This site is supported by donations to The OEIS Foundation.

User:R. J. Mathar/transforms3

From OeisWiki

Jump to: navigation, search
# This is a file with Maple programs. Usage from within another Maple program:
# read("transforms3") ;

# Read the b-file, named in the argument, and return the list of the
# numbers of its 2nd column. The information on the offset of the sequence
# is not transmitted: the first number in the list that is returned is just
# the topmost value in the b-file, which may be a(0), a(1) or another
# entry depending on the offset.
#
# Example to get a list of prime numbers after having downloaded b000040.txt
# from the OEIS web page:
#     a := BFILETOLIST("b000040.txt") ;
#
# @param bfilename The string with the input file which is read/scanned.
# @since 2010-09-06 Delete non-ASCII characters in the input
# Richard J. Mathar, 2010-02-20
BFILETOLIST := proc(bfilename)
        local a,f,bline,b,nold,n ;
        # List intially empty
        nold := -999999 ;
        a := [] ;
        f := fopen(bfilename,READ) :
        bline := readline(f) :
        # Loop over the entire b-file
        while bline <> 0 do
                # There have been files with non-ASCII characters in the past
                # (BOM markers of UTF codes in particular) which are difficult to detect.
                # The following is an attempt to treat this error case gently,
                # deleting non-ASCII bytes completely.
                bline := StringTools[Select](StringTools[IsASCII],bline) ;
                bline := StringTools[TrimLeft](bline) ;
                # Skip lines with comments
                if StringTools[FirstFromLeft]("#",bline ) = 1 then
                        ;
                else
                        # interpret the line contents  as two integers (index and value)
                        b := sscanf( bline,"%d %d") ;
                        # grab the second component of each line, the value; append
                        if nops(b) >1 then
                                a := [op(a),op(2,b)] ;
                                n := op(1,b) ;
                                if n <> nold+1 and nold <> -999999 then
                                        printf("# discontinuous first column %d %d\n",nold,n) ;
                                end if;
                                nold := n ;
                        end if;
                end if ;
                bline := readline(f) ;
        end do:
        fclose(f) ;
        RETURN(a) ;
end proc:

# Write the list "L" to the file "bfilename" with the first number
# given the offset "offs".
#
# Warning: existing files are overwritten.
#
# Example:
#     L := [seq(ithprime(i),i=1..20)] ;
#     LISTTOBFILE("B000040.txt",L,1) ; # note capital B in the example, not in real life.
#
# Richard J. Mathar, 2008-03-12
LISTTOBFILE := proc(bfilename,L,offs)
        local n,f ;
        f := fopen(bfilename,WRITE) :
        # Loop over the list members
        for n from 1 to nops(L) do
                fprintf(f,"%d %d\n",n+offs-1, op(n,L) ) ;
        end do:
        fclose(f) ;
        printf("%%I A%6s\n", substring(bfilename,2..7) ) ;
        printf("%%H A%6s F. Lastname, <a href=\"/A%6s/%s\">Table of n, a(n) for n=%d..%d</a>\n",substring(bfilename,2..7),substring(bfilename,2..7), bfilename,offs, nops(L)+offs-1) ;
        RETURN() ;
end proc:

# Read a b-file and replace tabs by spaces and <cr><lf> combinations
# by <lf>. Remove empty lines.
# The argument bfilename should be a readable b-file name b??????.txt
# (no directory component). The result will be placed in a file named
# B??????.txt (so that's just the first letter replaced by a capital-B), which is
# overwritten if existing.
# Richard J. Mathar, 2009-03-10
BFILE2BFILE := proc(bfilename)
        local f,fo,bline,b,n,boutname ;
        f := fopen(bfilename,READ) :

        boutname := sprintf("B%6s.txt", substring(bfilename,2..7) ) ;
        fo := fopen(boutname,WRITE) :

        bline := readline(f) :
        # Loop over the entire b-file; eliminate trailing blank lines
        while bline <> 0 do
                # copy lines with comments
                if StringTools[FirstFromLeft]("#",bline ) <> 0 then
                        fprintf(fo,"%s\n",bline ) ;
                else
                        # interpret the line contents  as two integers (index and value)
                        b := sscanf( bline,"%d %d") ;
                        if nops(b) >=2 then
                                # grab the second component of each line, the value; append
                                fprintf(fo,"%d %d\n",op(1,b), op(2,b) ) ;
                        fi;
                end if ;
                bline := readline(f) ;
        end do:
        fclose(f) ;
        fclose(fo) ;
        RETURN() ;
end:

# Convert a floating point constant to a sequence of decimal digits.
# @param c the floating point constant
# @return the list with zeros retained according to OEIS guidelines.
#           The length of the list depends on the value of Digits
#           at the time of the call. Two digits are chopped off as we assume
#           that the value of the constant would typically be rounded, not
#           truncated.
#
# Examples:
#       CONSTTOLIST(3.14159) ; # returns [3,1,4,1,5,9,0,0,...]
#       CONSTTOLIST(59) ; # returns [5,9,0,0,0,..]
#
# Richard J. Mathar, 2008-09-19
CONSTTOLIST := proc(c)
        local absc,a,dgs;
        # ensure it is the positive value, if negative
        absc := abs(c) ;
        a := [] ;
        # for numbers < 0.1, start with the zeros to maintain offset 0
        while absc < 0.1 do
                a := [op(a),0] ;
                absc := 10*absc ;
        end do:

        # dgs is the number of digits of absc before the decimal point, the number that remains to be printed
        dgs := 1+ilog10(absc) ;
        # dgs is the number of places we want to shift absc to the right, truncating two more digits than DIGITS
        dgs := Digits-2-dgs ;
        absc := floor(absc*10^dgs) ;
        a := [op(a),op(ListTools[Reverse](convert(absc,base,10)))] ;
        RETURN(a) ;
end proc:
# Convert a floating point constant to a sequence of digits in base b.
# @param c the floating point constant
# @param b the base (2, 3,... )
# @return the list with zeros retained according to OEIS guidelines.
#           The length of the list depends on the value of Digits
#           at the time of the call. Two digits are chopped off as we assume
#           that the value of the constant would typically be rounded, not
#           truncated.
#
# Examples:
#       CONSTTOLIST(evalf(Pi/3),2) ; # returns 1,0,0,0,0,1,... (See A019670)
#       CONSTTOLIST(59,2) ;
#
# Richard J. Mathar, 2009-02-06
CONSTTOLISTB := proc(c,b)
        local absc,a,dgs;
        # ensure it is the positive value, if negative
        absc := abs(c) ;
        a := [] ;
        # for numbers < 0.1, start with the zeros to maintain offset 0
        while absc < 1/b do
                a := [op(a),0] ;
                absc := b*absc ;
        end do:

        # dgs is the number of digits of absc before the decimal point, the number that remains to be printed
        dgs := 1+floor(log(absc)/log(b)) ;
        # dgs is the number of places we want to shift absc to the right, truncating two more digits than DIGITS
        dgs := Digits-2-dgs ;
        absc := floor(absc*b^dgs) ;
        a := [op(a),op(ListTools[Reverse](convert(absc,base,b)))] ;
        RETURN(a) ;
end proc:

# Convert a sequence to a floating point constant.
# The arguments are the list and the (correct) offset.
# Richard J. Mathar, 2009-02-02
LISTTOCONST := proc(L,offs)
        local c,n,dgs,realdigs ;
        c := 0 ;
        dgs := 0 ;
        for n from 1 to nops(L) do
                if dgs = 0 then
                        if op(n,L) <> 0 then
                                dgs := 1 ;
                        end if;
                else
                        dgs := dgs+1  ;
                end if;
                c := 10*c+op(n,L) ;
        end do:

        # Set effective precision for the floating point conversion to
        # what the list's length (after removing leading zeros) offers.
        realdigs := Digits ;
        Digits := dgs ;
        c := evalf( c*10^(offs-nops(L)) ) ;
        Digits := realdigs ;
        RETURN(c) ;
end proc:
# Motzkin transform of the sequence provided by the list L.
# This is equivalent to replacing the independent variable x in the
# generating function for L by x*M(x), where M(x) is the generating
# function for the Motzkin numbers.
#
# Examples:
#    L := [0,1,0,0,0,0,0,0,0] ; Motzkin(L) ; # returns [0,1,1,2,4,9,21,51,27]
#    L := [0,1,0,1,0,1,0,1,0,1] ; Motzkin(L) ; # returns [0,1,1,3,7,19,51,141,393,1107,...]
#    L := [1,0,1,0,1,0,1,0,1] ; Motzkin(L) ; # returns [1,0,1,2,6,16,45,126,357,1016,2907,...]
#
# Richard J. Mathar, 2008-11-07
MOTZKIN := proc(L)
        local d, bgf, mL, maxd;
        # output list intially empty

        # maximum degree of output is limited by initial length of L
        maxd := nops(L)-1 ;
        # generate an auxiliary list of Motzkin numbers
        mL := taylor( (1-x-sqrt(1-2*x-3*x^2))/2/x, x=0,maxd+2) ;

        bgf := op(1,L) ;
        for d from 1 to maxd do
                bgf := bgf+taylor(op(d+1,L)*mL^d,x=0,maxd+1) ;
        end do:
        [seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:

# Inverse Motzkin transform of the sequence provided by the list L.
# The inverse of the operation given by MOTZKIN above.
#
# Richard J. Mathar, 2008-11-08
MOTZKINi := proc(L)
        local d, bgf, mL, maxd;
        # output list intially empty

        # maximum degree of output is limited by initial length of L
        maxd := nops(L)-1 ;
        # generate an auxiliary list of Motzkin numbers
        mL := taylor( x/(1+x+x^2), x=0,maxd+2) ;

        bgf := op(1,L) ;
        for d from 1 to maxd do
                bgf := bgf+taylor(op(d+1,L)*mL^d,x=0,maxd+1) ;
        end do:
        [seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:

# Catalan transform of the sequence provided by the list L.
# This is equivalent to replacing the independent variable x in the
# generating function for L by x*C(x), where C(x) is the generating
# function for the Catalan numbers.
#
# Richard J. Mathar, 2008-12-17
CATALAN := proc(L)
        local d, bgf, cL, maxd;
        # output list intially empty

        # maximum degree of output is limited by initial length of L
        maxd := nops(L)-1 ;
        # generate an auxiliary list of Catalan numbers
        cL := taylor( (1-sqrt(1-4*x))/2, x=0,maxd+2) ;

        bgf := op(1,L) ;
        for d from 1 to maxd do
                bgf := bgf+taylor(op(d+1,L)*cL^d,x=0,maxd+1) ;
        end do:
        [seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:
# Inverse Catalan transform of the sequence provided by the list L.
# The inverse of the operation given by CATALAN above.
# Richard J. Mathar, 2008-12-11
CATALANi := proc(L)
        local d, bgf, cL, maxd;
        # output list intially empty

        # maximum degree of output is limited by initial length of L
        maxd := nops(L)-1 ;
        cL := taylor( x*(1-x), x=0,maxd+2) ;

        bgf := op(1,L) ;
        for d from 1 to maxd do
                bgf := bgf+taylor(op(d+1,L)*cL^d,x=0,maxd+1) ;
        end do:
        # bug in Maple 9: the following line is necessary to get correct results
        bgf := taylor(bgf,x=0,maxd+1) ;
        [seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:

# Pierce sequence generated from a constant.
# The input is a number between 0 and 1, the output the list of Pierce
# expansion coefficients.
# Richard J. Mathar, 2009-01-21
PIERCE := proc(c)
        local resid,L,i,an ;
        resid := c ;
        L := [] ;
        if whattype(c) = `float` then
                if c >0 and c <= 1 then
                        for i from 1 do
                                an := floor(1./resid) ;
                                L := [op(L),an] ;
                                resid := evalf(1.-an*resid) ;
                                if ilog10( mul(i,i=L)) > 0.7*Digits then
                                        break ;
                                end if ;
                        end do:
                end if;
        end if;
        L ;
end proc:

# Engel sequence generated from a constant.
# The input is a number between 0 and 1, the output the list of Engel
# expansion coefficients.
# Richard J. Mathar, 2009-07-30
ENGEL := proc(c)
        local resid,L,i,an ;
        resid := c ;
        L := [] ;
        if whattype(c) = `float` then
                if c > 1 then
                        RETURN( [1,op(procname(c-1.))] ) ;
                elif c >0 and c <= 1 then
                        for i from 1 do
                                an := ceil(1./resid) ;
                                L := [op(L),an] ;
                                resid := evalf(an*resid-1.) ;
                                if ilog10( mul(i,i=L)) > 0.7*Digits then
                                        break ;
                                end if ;
                        end do:
                end if;
        end if;
        L ;
end proc:

# Beatty sequence generated from a constant.
# Richard J. Mathar, 2009-05-20
BEATTY := proc(c)
        local L,n,x ;
        L := [] ;
        x := evalf(c) ;
        for n from 1 to 0.8*Digits do
                L := [op(L),floor(n*x)] ;
        end do:
        L ;
end proc:


# L-Schroeder transform of the sequence provided by the list L.
# This is equivalent to replacing the independent variable x in the
# generating function for L by x*S(x), where S(x) is the generating
# function for the large Schroeder numbers.
#
# Richard J. Mathar, 2008-11-13
SCHROEDER := proc(L)
        local d, bgf, mL, maxd;
        # output list intially empty

        # maximum degree of output is limited by initial length of L
        maxd := nops(L)-1 ;
        # generate an auxiliary list of Schroeder numbers
        mL := taylor( (1-x-sqrt(1-6*x+x^2))/2, x=0,maxd+2) ;

        bgf := op(1,L) ;
        for d from 1 to maxd do
                bgf := bgf+taylor(op(d+1,L)*mL^d,x=0,maxd+1) ;
        end do:
        [seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:

# Inverse L-Schroeder transform of the sequence provided by the list L.
# This is equivalent to replacing the independent variable x in the
# generating function for L by x*(1-x)/(1+x).
#
# Example: if L is A001850, the funtion returns a list of 1 followed
#   by A001333.
# Richard J. Mathar, 2008-11-13
SCHROEDERi := proc(L)
        local d, bgf, mL, maxd;
        # output list intially empty

        # maximum degree of output is limited by initial length of L
        maxd := nops(L)-1 ;
        # generate an auxiliary list of Schroeder numbers
        mL := taylor( (1-x)*x/(1+x), x=0,maxd+2) ;

        bgf := op(1,L) ;
        for d from 1 to maxd do
                bgf := bgf+taylor(op(d+1,L)*mL^d,x=0,maxd+1) ;
        end do:
        [seq(coeftayl(bgf,x=0,d),d=0..maxd)] ;
end proc:

# Reduce the list given by the argument L to a list in which immediately repeated
# numbers (ie, clusters of the same numbers) are reduced to a single occurrence.
#
# Example: if L contains [3,5,6,7,7,8,3,3,3,10],
# then the list returned is [3,5,6,7,8,3,10], where "7,7" has shrunk to "7"
# and "3,3,3" has shrunk to "3".
#
# Richard J. Mathar, 2008-05-03
LISTNODUPL := proc(L)
        local b, a, aprev;
        # output list intially empty
        b := [] ;
        # Note that aprev is not initialized to avoid a coincidental equality with L[1].
        for a in L do
                if a <> aprev then
                        b := [op(b),a] ;
                end if ;
                aprev := a ;
        end do:
        RETURN(b) ;
end proc:

# Reduce the list given by the argument "L" to a list which only contains clusters
# of 2 or more repeated numbers. The argument "multpl" is either a number
# larger than 0 which specifies the multiplicity (cluster length) which the
# repetitions need to have (exactly) to enter the new list, or is 0 (zero)
# to indicate that any cluster length of 2 or more will have one representative
# be put on the list returned.
#
# Example: If L contains [3,5,7,7,8,3,3,3,10,7,7,7,11]
# then the list returned by LISTDUPL(L,2) is [7], which represents the
# first of the two clusters with 7's.
# Example: If L contains [3,5,7,7,8,3,3,3,10,7,7,7,11]
# then the list returned by LISTDUPL(L,3) is [3,7], which represents the
# cluster with the 3's and the 2nd cluster with the 7's, each showing
# numbers repeated 3 times.
# Example: If L contains [3,5,7,7,8,3,3,3,10,7,7,7,11,11]
# then the list returned by LISTDUPL(L,0) is [7,3,7,11], which represents the
# cluster with the two 7's, the 2nd cluster with the 3's, the 2nd cluster
# with 7's, and the final cluster with 11's.
# Example: If L contains [3,5,7,7,8,3,3,3,10,7,7,7,11,11]
# then the list returned by LISTDUPL(L,1) is [3,5,8,10], which represents the
# isolated numbers (numbers in clusters of length 1).
#
# To obtain results which register repeated numbers without regard to
# their placement in the original list, these may be sorted before calling
# the routine: LISTDUPL(sort(L),multpl) ;
#
# Richard J. Mathar, 2008-05-03
LISTDUPL := proc(L,multpl)
        local b, a, aprev, cnt;
        if multpl < 0 then
                ERROR("negative multiplicity ",args[2]) ;
        end if ;
        # output list intially empty
        b := [] ;
        cnt :=0 ;
        aprev := op(1,L) ;
        for a in L do
                if a <> aprev then
                        if ( cnt = multpl and multpl >0 ) or (cnt > 1 and multpl = 0 ) then
                                b := [op(b),aprev] ;
                        end if ;
                        cnt := 1 ;
                else
                        cnt := cnt+1 ;
                end if ;
                aprev := a ;
        end do:
        # At that time, aprev represents a number in the last cluster. Since
        # the lists in the OEIS are generally incomplete at the end, the following
        # check on the length of this last cluster is not executed, because
        # the length of this cluster might actually be larger than the list
        # representatives tell; this may lead to misinterpretations if the lists
        # are "full"!
        # if ( cnt = multpl and multpl >0 ) or (cnt > 1 and multpl = 0 ) then
        #               b := [op(b),aprev] ;
        #       end if ;
        if cnt > 1 and multpl = 0 then
                b := [op(b),aprev] ;
        end if ;
        RETURN(b) ;
end proc:

# Lehmer's cotangent representation of an irrational number.
# Example
# CONSTTOCOT(evalf(Pi)); # see A002667
# CONSTTOCOT(evalf(log(2))); # see A081785
# Fritz Schweiger, Metrische Ergebnisse ueber den Kotangensalgorithmus, Acta Arith. 26 (1975) 217
CONSTTOCOT := proc(x)
        local L,xcut ;
        xcut := x ;
        L := [floor(xcut)] ;
        while type(L[-1],'infinity') = false do
                arccot(L[-1])-arccot(xcut) ;
                xcut := cot(evalf(%) ) ;
                L := [op(L),floor(xcut)] ;
        end do;
        return subsop(-1=NULL,L) ;
end proc:
# LOG transformation operation on o.g.f.
# The input list L represents [a0=1,a1,a2,a3,a4....] of an ordinary
# generating function, g=1+a1*x+a2*x^2+.... This is to say that the first
# element of L is tacitly assumed to be equal to 1.
#
# The output contains the coefficients of log(g) = sum_{i=1..inf} b(i)*x^i/i .
# Richard J. Mathar, 2009-03-13
OLOG := proc(L)
        local x,gf,gfl,i,a;
        gf := add( op(i,L)*x^(i-1),i=2..nops(L)) ;
        gfl := 0 ;
        for i from 1 to nops(L) do
                gfl := taylor(gfl+gf^i/i,x=0,nops(L)) ;
        end do:
        a := [] ;
        for i from 1 to nops(L)-1 do
                a := [op(a), i*coeftayl(gfl,x=0,i)] ;
        end do:
        a;
end proc:

# Akiyama-Tanigawa algorithm, one step
# The input list L with L[0], L[1] is transformed to the output list
# A[0], A[1] as described in
# D. Merlini, R. Sprugnoli, M. C. Verri, The Akiyama-Tanigawa Transformation, Integers, 5 (1) (2005) #A05
# Richard J. Mathar 2010-12-02
AKIYATANI := proc(L)
        local a,k ;
        a := [] ;
        for k from 1 to nops(L)-1 do
                a := [op(a),k*(op(k,L)-op(k+1,L))] ;
        end do;
        a ;
end proc:

# Akiyama-Tanigawa algorithm, generate first column.
# The input list L is iteratively transformed to construct the first columnof the table.
# D. Merlini, R. Sprugnoli, M. C. Verri, The Akiyama-Tanigawa Transformation, Integers, 5 (1) (2005) #A05
# Richard J. Mathar 2010-12-02
AKIYAMATANIGAWA := proc(L)
        local a,r,img ;
        a := [op(1,L)] ;
        img := L ;
        for r from 1 to nops(L)-1 do
                img := AKIYATANI(img) ;
                a := [op(a),op(1,img)] ;
        end do;
        a ;
end proc:

# Inverse Akiyama-Tanigawa tabulation: generate top row given the column k=0.
# Richard J. Mathar 2010-12-02
AKIYAMATANIGAWAi := proc(L)
        local a,r,img,k,row ;
        img := [op(-1,L)] ;
        for r from nops(L)-1 to 1 by -1 do
                row := [op(r,L)] ;
                for k from 2 to 1+nops(L)-r do
                        op(k-1,row)-op(k-1,img)/(k-1) ;
                        row := [op(row),%] ;
                end do;
                img := row ;
        end do;
        img ;
end proc:

# Given the left entry d(x) and the right entry h(x) of the Riordan array,
# compute the (n,k) element of the array.
# @param d a function of x
# @param h a function of x
# @param n the row index >=0
# @param k the column index >=0
# Richard J. Mathar 2011-01-09
RIORDAN := proc(d,h,n,k)
        d*h^k ;
        expand(%) ;
        coeftayl(%,x=0,n) ;
end proc:

RIORDANEXP := proc(d,h,n,k)
        RIORDAN(d,h,n,k)*n!/k! ;
end proc:

# Generate a list of multinomial partitions.
# @param n The sum of the first moment of the parts, sum_j j*p_j = n
# @param k The sum of the parts, sum_j p_j = k
# @param nparts The number of parts, sum_j 1 = nparts
# @param zIsPart If true, the constrained is 0<=p_j, other wise 1<=p_j.
# @return A list of the form [[p_1,p_2,p_3,...],[p_1,p_2,p_3,...]]
#   where each of the sublists is a partition in nparts obeying the two constraints.
# There no other restrictions but p_j<=k. Some of the p_j may be equal.
# Richard J. Mathar, 2011-04-06
multinomial2q := proc(n::integer,k::integer,nparts::integer,zIsPart::boolean)
        local lpar ,res, constrp;
        res := [] ;
        if n< 0 or nparts <= 0 then
                ;
        elif nparts = 1 then
                if n = k then
                        if n > 0 or zIsPart then
                                return [[n]] ;
                        end if;
                end if;
        else
                for lpar from `if`(zIsPart,0,1) do
                        if lpar*nparts > n or lpar > k then
                                break;
                        end if;
                        # recursive call, assuming that lpar is the last part
                        for constrp in procname(n-nparts*lpar,k-lpar,nparts-1,zIsPart) do
                                if nops(constrp) > 0 then
                                        res := [op(res),[op(constrp),lpar]] ;
                                end if;
                        end do:
                end do:
        end if ;
        return res ;
end proc:

# Generate a list of multinomial partitions.
# @param n The sum of the first moment of the parts, sum_j j*p_j = n
# @param k The sum of the parts, sum_j p_j = k
# @param zIsPart If true, the constrained is 0<=p_j, other wise 1<=p_j.
# @return A list of the form [[p_1,p_2,p_3,...],[p_1,p_2,p_3,...]]
#   where each of the sublists is a partition obeying the two constraints.
# There are no other restrictions but 1<=p_j<=k. Some of the p_j may be equal.
# The Bell polynomials of the second kind are an application of this.
# Note that ordering in the sublists matters.
# Empty partitions (in case n=k=0) are not returned as empty sublists but not returned at all.
# Example:
#   multinomial2(10,5,false) returns [[1,3,1],[2,1,2]] where 1+3+1=5, 1+2*3+3*1= 10
#   for the first sublist and 2+1+2=5, 2+2*1+3*2=10 for the second sublist.
# Richard J. Mathar, 2011-04-06
multinomial2 := proc(n::integer,k::integer,zIsPart::boolean)
        local npart,res,constrp ;
        res := [] ;
        for npart from 1 do
                # The smallest part is 1. So sum_j j*p_j >= npart*(npart+1)/2 =n means
                # that the number of parts is restricted by n, only to a lesser degree
                # by k.
                if npart*(npart+1)/2 > n or npart > k then
                        return res;
                end if;
                # loop over the partitions with fixed number of parts and
                # attach those solultions to the result vector
                for constrp in multinomial2q(n,k,npart,zIsPart) do
                        if nops(constrp) > 0 then
                                res := [op(res),constrp] ;
                        end if ;
                end do:
        end do:
end proc:

# Bell polynomial of the second kind.
# @param n
# @param k
# @param x The vector of x1,x2,x3...
# @return sum_{j1+j2+..=k, j1+2*j2+.. =n} n!/(j1!j2!..) (x1/1!)^j1 (x2/2!)^j2 ...
# Richard J. Mathar, 2011-04-06
bell2 := proc(n::integer,k::integer,x::list)
        local p,npart,b,m,j;
        npart := nops(x) ;
        b := 0 ;
        for p in multinomial2q(n,k,npart,true) do
                m := 1 ;
                for j from 1 to nops(p) do
                        m := m*(op(j,x)/j!)^op(j,p)/op(j,p)! ;
                end do:
                b := b+ n!*m ;
        end do:
        return b;
end proc:
Personal tools