#LEGO: A Maple package for computing generating functions for LEGO #configurations, accompanying the paper `Automatic Counting of LEGO #Towers', by Doron Zeilberger #Written by Doron Zeilberger, Temple University ,zeilberg@math.temple.edu. #Please report bugs to zeilberg@math.temple.edu print(`LEGO: a package for computing generating functions for LEGO`): print(`Configurations`): print(`Version of Jan. 31, 1997`): print(`written by Doron Zeilberger(zeilberg@math.temple.edu).`): print(`This program accompanies the paper "Automatic Counting of LEGO`): print(`Towers" by Doron Zeilberger`): print(`The most current version of the paper and program`): print(` are available from`): print(`http://www.math.temple.edu/~zeilberg`): print(`For a list of the procedures type ezra(), for help with`): print(`a specific procedure, type ezra(procedure_name)`): print(``): print(`For a list of procedures that are not mentioned explicity`): print(`in the paper, type ezra1();`): print(`and for specific help ezra1(procedure_name)`): ezra:=proc() if args=NULL then print(`Contains the following procedures:`): print(`LEGO,muLEGO, MIGDAL, LEGOmul, muLEGOmul `): fi: if nops([args])=1 and op(1,[args])=`MIGDAL` then print(`MIGDAL(R,t): g.f., in t, whose coeff. of t^n is the number of`): print(`LEGO towers in which the pieces are a by i, a>=1 1<=i<=R, and`): print(`such that the a sides are all parallel to each other`): fi: if nops([args])=1 and op(1,[args])=`LEGO` then print(`LEGO(L,p,a,b,t): `): print(`Given a linear polynomial L(a) in a, and a polynomial`): print(`p(a,b), computes the generating function sum(F(a),a=1..infinity)`): print(`where F(a)=t^L(a)+t^L(a)*(sum(p(a,b)*F(b),b=1..infinity))`): print(``): print(`For examaple, try LEGO(a,a+b-1,a,b,t)`): fi: if nops([args])=1 and op(1,[args])=`LEGOmul` then print(`LEGOmul(L,p,a,b,t): `): print(`Given lists of variables, a and b s.t. nops(a)=nops(b)=s, a var. t`): print(`and a polynomial p in the variables a_1, ..., a_s, b_1, ..., b_s`): print(`and a linear form L in a_1, ..., a_s`): print(`computes the generating function`): print(` sum(F(a_1, ..., a_s),a_1, ..., a_s=1..infinity)`): print(`where F(a_1, ..., a_s)=t^L(a_1, ..., a_s)+`): print(`t^L(a_1, ..., a_s)*(sum(p(a_1, ..., a_s,b_1, ..., b_s)`): print(`*F(b_1, ..., b_s),b_1, ..., b_s=1..infinity))`): print(``): print(`For examaple, try `): print(`LEGOmul(a1+a2,(a1+b1-1)*(a2+b2-1),[a1,a2],[b1,b2],t)`): fi: if nops([args])=1 and op(1,[args])=`muLEGO` then print(`muLEGO(Ls,pols,a,b,t): `): print(`Given a list of affine-linear polynomials`): print(` Ls=[L[1],...,L[s]] in the multi-variables a`): print(` and a list of lists of polynomials, pols: `): print(`[[p[1,1],p[1,2],...,p[1,s]],[p[2,1],p[2,2],...,p[2,s]], ... `): print(`[p[s,1],...,p[s,s]]], in the variables a and b`): print(` computes the generating function`): print(`sum(sum(F[i](a),a=1..infinity),i=1..s) where the F[i] satisfy`): print(`the system, for i=1..s:`): print(`where F[i](a)=t^L[i](a)+t^L[i](a)*`): print(`sum(sum(p[i,j](a,b)*F[j](b),b=1..infinity), j=1..s)`): print(``): print(`For examaple, try muLEGO([2*a,2*a-1],[[2*a-1,2*a],`): print(`[2*a-2,2*a-1]] ,a,b,t);`): lprint(): print(`When it fails, it returns 0`): fi: if nops([args])=1 and op(1,[args])=`muLEGOmul` then print(`muLEGOmul(Dims,Ls,pols,a,b,t): `): print(`Given a list of dimensions Dims [d_1, ..., d_s]`): print(`and given a list of affine-linear polynomials`): print(` Ls=[L[1],...,L[s]], where L[i] depends on a[1], ..., a[d_i]`): print(` and a list of lists of polynomials, pols: `): print(`[[p[1,1],p[1,2],...,p[1,s]],[p[2,1],p[2,2],...,p[2,s]], ... `): print(`[p[s,1],...,p[s,s]]], where p[i,j] depends on the d_i+d_j vars.`): print(`a[1], ..., a[d_i], b[1], ..., b[d_j]`): print(` computes the generating function`): print(`sum(sum(F[i](a_1, ..., a_d_i),a_1, ..., a_d_i=1..infinity),i=1..s)`): print(` where the F[i] satisfy the system, for i=1..s:`): print(`where F[i](a_1, ..., a_d_i)=t^L[i](a_1, ..., a_d_i)+t^L[i](a)*`): print(`sum(sum(p[i,j](a_1, ..., a_d_i,b_1, ..., b_d_j)*`): print(`F[j](b_1, ..., b_d_j),b_1, ..., b_d_j=1..infinity), j=1..s)`): print(``): print(` a and b are names for the indexed variables, and t is the`): print(`designated variable for the generating function`): print(`For examaple, try muLEGOmul([2,1],[2*a[1]+2*a[2],2*a[1]-1]`): print(`,[[(a[1]+a[2]-1)*(b[1]+b[2]-1), a[1]*a[2]*b[1]],`): print(`[b[1]*a[1]+a[2],a[1]*b[1]]] ,a,b,t)`): print(`or try`): print(`muLEGOmul([2,2],[a[1]+a[2],2*a[1]+2*a[2]],`): print(`[[(a[1]+b[1]-1)*(a[2]+b[2]-1), (a[1]+b[1]-1)*(a[2]+b[2]-1)],`): print(`[(a[1]+b[1]-1)*(a[2]+b[2]-1), (a[1]+b[1]-1)*(a[2]+b[2]-1)]]`): print(`,a,b,t);`): print(`or try`): print(`muLEGOmul([1,2],[a[1]+1,a[1]+a[2]],`): print(`[[a[1]+b[1]-1,(a[1]+b[1]-1)*b[2]],[a[2]*(a[1]+b[1]-1),`): print(`(a[2]+b[2]-1)*(a[1]+b[1]-1)]],a,b,t);`): print(`if it is unsuccesful, it returns 0`): fi: end: ezra1:=proc() if args=NULL then print(`The procedures not directly refered to in the paper are: `): print(` Compo, COMP, wt, empirLEGO,sempirLEGO, FLar, FLr `): print(` sempirLEGOpq, SeqsempirLEGOpq, SeqsempirLEGO, Findrec `): print(` gfmul, gfpmul, Empir `): fi: if nops([args])=1 and op(1,[args])=`FLr` then print(`FLr:=proc(L,r,p,a,b) :given a linear form L(a), a polynomial`): print(`p in the variables a and b, and integer r, computes the`): print(`sum of the wt of all compositions [a_1 ,..., a_s]`): print(`such that L(a_1)+ ... +L(a_s)=r`): fi: if nops([args])=1 and op(1,[args])=`FLar` then print(`FLar:=proc(L,a1,r,p,a,b) :given a linear form L(a), a polynomial`): print(`p in the variables a and b, and integers a1 and r, computes the`): print(`sum of the wt of all compositions [a_1 ,..., a_s]`): print(`such that L(a_1)+ ... +L(a_s)=r`): fi: if nops([args])=1 and op(1,[args])=`empirLEGO` then print(`empirLEGO(L,p,a,b,t,NUTERMS) finds empirically the first`): print(`NUTERMS terms of the generation function for weighted`): print(`compositions according to the linear form a and the polynomial`): print(`p(a,b)`): print(``): print(`For examaple, try empirLEGO(a,a+b-1,a,b,t,10)`): fi: if nops([args])=1 and op(1,[args])=`sempirLEGO` then print(`sempirLEGO(L,p,a,b,t,NUTERMS) finds semi-empirically the first`): print(`NUTERMS terms of the generation function for weighted`): print(`compositions according to the linear form a and the polynomial`): print(`p(a,b)`): print(``): print(`For examaple, try sempirLEGO(a,a+b-1,a,b,t,10)`): fi: if nops([args])=1 and op(1,[args])=`wt` then print(`wt(com1,p,a,b): finds the weight of the composition com1 according`): print(`to p(a,b) i.e. wt([a_1,a_2, .., a_s])=`): print(`p(a_1,a_2)*p(a_2,a_3)*...*p(a_(s-1),a_s)`): fi: if nops([args])=1 and op(1,[args])=`COMP` then print(`COMP(r,L,a): gives the set of all compositions`): print(` [a_1, ..., a_s] `): print(`s.t. L(a_1)+ ... +L(a_s)=r, where L has the form`): print(`L(a)=c1*a+c2 where c1 c2 are constants`): print(`a is the name of the variable`): fi: if nops([args])=1 and op(1,[args])=`Compo` then print(`Compo(a1,r,L,a): gives the set of all compositions [a_1, ..., a_s]`): print(` that start`): print(`with a_1=a1 and s.t. L(a_1)+ ... +L(a_s)=r, where L has the form`): print(`L(a)=c1*a+c2 where c1 c2 are constants`): print(`a is the name of the variable`): fi: if nops([args])=1 and op(1,[args])=`gfmul` then print(`gfmul(L,a,t,z); Given a linear polynomial L(a),`): print(`in the list of variables a, and a list of variables`): print(`z, of the same length as a, computes`): print(`the rational function `): print(`sum(t^L(a_1, ..., a_s)*z_1^a_1*z_2^a2* ... *z_s^a_s,`): print(`a_1, ..., a_s=1..infinity), here s is the length of a or z`): fi: if nops([args])=1 and op(1,[args])=`Empir` then print(`Empir(gu,x,P) empirically finds an algebraic equation`): print(` F(P(x),x)=0 for the `): print(`formal power series P(x):=sum_{i=0} gu[i]*x^i, where gu is a list`): print(`for example Empir([seq(1,i=1..20)],x,P) should yield`): print(` P-xP-1=0`): print(`If there is not enough data, it reurns 0. You should then`): print(`try it again with a longer sequence`): fi: if nops([args])=1 and op(1,[args])=`Findrec` then print(`Findrec(f,n,N) finds empirically an ordi. linear recurrence`): print(` with polynomial coeffs. The input is a sequence f given as a list`): print(`STARTING at f[1],i.e. f[0] is not considered`): print(` in n and N, where N is the forward unit shift: Nf(n):=f(n+1).`): print(`For example Findrec([1,1,2,3,5,8,13,21,34],n,N) should yield`): print(`N^2-N-1 , and findrec([1,2,5,14,42,132,429],m,M) should yield`): print(`(m+1)*M-(4*m-2). If there is not enough data, you will get 0`): print(`Of course, if a sequence is not P-recursive there would never`): print(`enough data`): fi: if nops([args])=1 and op(1,[args])=`sempirLEGOpq` then print(`sempirLEGOpq(L,p,q,a,b,t,NUTERMS) finds semi-empirically the first`): print(`NUTERMS terms of the generation function for weighted`): print(`compositions according to the linear form a and the polynomials`): print(`p(a,b) and q(a,b).`): print(` Output as a truncated Taylor series.`): fi: if nops([args])=1 and op(1,[args])=`SeqsempirLEGOpq` then print(`sempirSeqLEGOpq(L,p,q,a,b,NUTERMS) finds semi-empirically the first`): print(`NUTERMS terms of the generation function for weighted`): print(`compositions according to the linear form a and the polynomials`): print(`p(a,b) and q(a,b).`): print(`Output as a sequence`): fi: if nops([args])=1 and op(1,[args])=`SeqsempirLEGO` then print(`sempirSeqLEGO(L,p,a,b,NUTERMS) finds semi-empirically the first`): print(`NUTERMS terms of the generation function for weighted`): print(`compositions according to the linear form a and the polynomials`): print(`p(a,b).`): print(`Output as a sequence`): fi: if nops([args])=1 and op(1,[args])=`gfpmul` then print(`gfpmul(L,a,r,t,z); Given a linear polynomial L(a) in the list of`): print(`variables a, and a list of nonegative integers r, and a variable t`): print(`and a list of variables z s.t. nops(a)=nops(r)=nops(z),`): print(` computes the rational function `): print(`sum(a_1^r_1*a_2^r_2*...*a_s^r_s*t^L(a_1, .., a_s)*`): print(`z_1^a_1*z_2^a_2* ... *z_s^a_s,a_1, ..., a_s=1..infinity)`): print(`where s is the length of a (=nops(r)=nops(z))`): fi: end: #gfp(L,a,r,t,z); Given a linear polynomial L(a) computes #the rational function sum(a^r*t^L(a)*z^a,a=1..infinity) gfp:=proc(L,a,r,t,z) local gu: if not type(r,integer) or r<0 then ERROR(`Bad intput`): fi: if degree(L,a)>1 then ERROR(`the first input much be a poly of degree <=1 in`,a): fi: if r=0 then RETURN(gf(L,a,t,z)): fi: gu:=gfp(L,a,r-1,t,z): gu:=z*diff(gu,z): normal(gu): end: #gf(L,a,t,z); Given a linear polynomial L(a) computes #the rational function sum(t^L(a)*z^a,a=1..infinity) gf:=proc(L,a,t,z) local k0,k1: if degree(L,a)>1 then ERROR(`the first input much be a poly of degree <=1 in`,a): fi: k0:=coeff(L,a,0): k1:=coeff(L,a,1): t^(k0+k1)*z/(1-t^k1*z): end: #gfpol(L,a,POL,t,z); Given a linear polynomial L(a), and a polynomial #POL, in a, computes #the rational function sum(POL(a)*t^L(a)*z^a,a=1..infinity) gfpol:=proc(L,a,POL,t,z) local gu,r: gu:=0: for r from 0 to degree(POL,a) do gu:=gu+coeff(POL,a,r)*gfp(L,a,r,t,z): gu:=normal(gu): od: gu: end: LEGO:=proc(L,p,a,b,t) local gu,mu,r,G,eq,var,z: gu:=gf(L,a,t,z): var:={}: for r from 0 to degree(p,b) do gu:=gu+G[r]*gfpol(L,a,coeff(p,b,r),t,z): var:=var union {G[r]}: od: eq:={G[0]=subs(z=1,gu)}: mu:=gu: for r from 1 to degree(p,b) do mu:=z*diff(mu,z): eq:=eq union {G[r]=subs(z=1,mu)}: od: var:=solve(eq,var): subs(var,G[0]): factor(normal(%)): end: muLEGO:=proc(Ls,pols,a,b,t) local deg,i,j,gu,mu,r,G,eq,var,p,hakhi,z: for j from 1 to nops(Ls) do hakhi[j]:=0: for i from 1 to nops(Ls) do deg:=degree(op(j,op(i,pols)),b): if deg>hakhi[j] then hakhi[j]:=deg: fi: od: od: if nops(Ls)<>nops(pols) then ERROR(`Ls and pols must have the same length`): fi: for i from 1 to nops(Ls) do gu[i]:=gf(op(i,Ls),a,t,z): var:={}: for j from 1 to nops(Ls) do p:=op(j,op(i,pols)): for r from 0 to degree(p,b) do gu[i]:=gu[i]+G[j,r]*gfpol(op(i,Ls),a,coeff(p,b,r),t,z): var:=var union {G[j,r]}: od: od: od: eq:={}: for i from 1 to nops(Ls) do eq:=eq union {G[i,0]=subs(z=1,gu[i])}: mu:=gu[i]: for r from 1 to hakhi[i] do mu:=z*diff(mu,z): eq:=eq union {G[i,r]=subs(z=1,mu)}: od: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: mu:=0: for i from 1 to nops(Ls) do mu:=mu+subs(var,G[i,0]): od: normal(mu): end: #Compo(a1,r,L,a): gives the set of all compositions [a_1, ..., a_s] that start #with a_1=a1 and s.t. L(a_1)+ ... +L(a_s)=r, where L has the form #L(a)=c1*a+c2 where c1 c2 are constants #a is the name of the variable Compo:=proc(a1,r,L,a) local i,gu,a2,mu: option remember: if degree(L,a)>1 then ERROR(`L is a linear form in`, a): fi: if r1 then ERROR(`L is a linear form in`, a): fi: gu:={}: for a1 from 1 while subs(a=a1,L)<=r do gu:=gu union Compo(a1,r,L,a): od: gu: end: #wt(com1,p,a,b): finds the weight of the composition com1 according #to p(a,b) i.e. wt([a_1,a_2, .., a_s])= #p(a_1,a_2)*p(a_2,a_3)*...*p(a_(s-1),a_s) wt:=proc(com1,p,a,b) local gu,i: gu:=1: for i from 1 to nops(com1)-1 do gu:=gu*subs({a=op(i,com1), b=op(i+1,com1)},p): od: gu: end: #empirLEGO(L,p,a,b,t,NUTERMS) finds empirically the first #NUTERMS terms of the generation function for weighted #compositions according to the linear form a and the polynomial #p(a,b) empirLEGO:=proc(L,p,a,b,t,NUTERMS) local r,gu,mu,i: mu:=0: for r from 1 to NUTERMS do gu:=COMP(r,L,a): for i from 1 to nops(gu) do mu:=mu+wt(op(i,gu),p,a,b)*t^r: od: od: mu: end: #FLar:=proc(L,a1,r,p,a,b) :given a linear form L(a), a polynomial #p in the variables a and b, and integers a1 and r, computes the #sum of the wt of all compositions [a_1 ,..., a_s] #such that L(a_1)+ ... +L(a_s)=r FLar:=proc(L,a1,r,p,a,b) local gu,b1: option remember: if subs(a=a1,L)>r then RETURN(0): fi: if subs(a=a1,L)=r then RETURN(1): fi: gu:=0: for b1 from 1 while subs(a=b1,L)+subs(a=a1,L)<=r do gu:=gu+subs({a=a1,b=b1},p)*FLar(L,b1,r-subs(a=a1,L),p,a,b): od: gu: end: #FLr:=proc(L,r,p,a,b) :given a linear form L(a), a polynomial #p in the variables a and b, and integer r, computes the #sum of the wt of all compositions [a_1 ,..., a_s] #such that L(a_1)+ ... +L(a_s)=r FLr:=proc(L,r,p,a,b) local gu,a1: gu:=0: for a1 from 1 while subs(a=a1,L)<=r do gu:=gu+FLar(L,a1,r,p,a,b): od: gu: end: #sempirLEGO(L,p,a,b,t,NUTERMS) finds semi-empirically the first #NUTERMS terms of the generation function for weighted #compositions according to the linear form a and the polynomial #p(a,b) sempirLEGO:=proc(L,p,a,b,t,NUTERMS) local gu,r: gu:=0: for r from 1 to NUTERMS do gu:=gu+FLr(L,r,p,a,b)*t^r: od: gu: end: #MIGDAL(R,t): g.f., in t, whose coeff. of t^n is the number of #LEGO towers in which the pieces are a by i, a>=1 1<=i<=R, and #such that the a sides are all parallel to each other MIGDAL:=proc(R,t) local i,j,Ls,pols,a,b,pols1: Ls:=[]: pols:=[]: for i from 1 to R do Ls:=[op(Ls),i*a]: pols1:=[]: for j from 1 to R do pols1:=[op(pols1),(i+j-1)*(a+b-1)]: od: pols:=[op(pols),pols1]: od: muLEGO(Ls,pols,a,b,t): end: #FLarpq:=proc(L,a1,r,p,q,a,b) :given a linear form L(a), a #polynomial p in the variables a and b, and integers a1 and r, #computes the sum of the wt of all compositions [a_1 ,..., a_s] #such that L(a_1)+ ... +L(a_s)=r and the weight is the product #of R(a_i,a_{i+1}) where R(a,b)=p(a,b) if ba FLarpq:=proc(L,a1,r,p,q,a,b) local gu,b1: option remember: if subs(a=a1,L)>r then RETURN(0): fi: if subs(a=a1,L)=r then RETURN(1): fi: gu:=0: for b1 from 1 while subs(a=b1,L)+subs(a=a1,L)<=r do if b1nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if nops(kv)>1 then RETURN(0): fi: for i from 1 to nops(kv) do ope:=subs(op(i,kv)=1,ope): od: ope: end: Findrec:=proc(f,n,N) local i1,kal,DEGREE, ORDER,L,gu,gu1,i: for L from 1 to nops(f)-5 do for ORDER from 1 to min(L,trunc((nops(f)-3)/2)-1) do for DEGREE from 0 to min(L-ORDER, trunc((nops(f)-3)/2/(1+ORDER))-2) do gu:=findrec(f,DEGREE,ORDER,n,N): if gu<>0 then gu1:=0: for i from 0 to degree(gu,N) do gu1:=gu1+factor(coeff(gu,N,i))*N^i: od: kal:=checkrec(gu1,n,N,f): if kal=[seq(0,i1=1..nops(kal))] then RETURN(gu1): fi: fi: od: od: od: 0: end: #checkrec(ope,n,N,Seq1) applies the operator ope(N,n) to the sequence gu checkrec:=proc(ope,n,N,Seq1) local i,j,lu,ru: lu:=[]: for i from 1 to nops(Seq1)-degree(ope,N) do ru:=0: for j from 0 to degree(ope,N) do ru:=ru+subs(n=i,coeff(ope,N,j))*op(i+j,Seq1): od: lu:=[op(lu),ru]: od: lu: end: empir:=proc(gu,degx,degP,x,P) local i1,i2,F,a,cand,lu,eq,var,mu,flo,pip: if (1+degx)*(1+degP) > nops(gu)-3 then RETURN(`sequence too small`): fi: F:=0: var:={}: for i1 from 0 to degx do for i2 from 0 to degP do F:=F+a[i1,i2]*x^i1*P^i2: var:=var union {a[i1,i2]}: od: od: cand:=0: for i1 from 0 to nops(gu)-1 do cand:=cand+op(i1+1,gu)*x^i1: od: lu:=subs(P=cand,F): lu:=taylor(lu,x=0,nops(gu)-1): eq:={}: for i1 from 0 to nops(gu)-2 do eq:=eq union {coeff(lu,x,i1)=0} od: mu:=solve(eq,var): F:=subs(mu,F): if F=0 then RETURN(0): fi: flo:=degree(F,P): pip:=coeff(F,P,flo): flo:=degree(pip,x): pip:=coeff(pip,x,flo): F:=F/pip: normal(F): end: Empir:=proc(gu,x,P) local degx,degP,L,lu: for L from 1 to (nops(gu)-3)/3 do for degP from 1 to L do for degx from 0 to min(trunc((nops(gu)-3)/(1+degP))-1,L-degP) do lu:=empir(gu,degx,degP,x,P): if lu<>0 then RETURN(lu): fi: od: od: od: 0: end: #gfmul(L,a,t,z); Given a linear polynomial L(a), #in the list of variables a, and a list of variables #z, of the same length as a, computes #the rational function #sum(t^L(a_1, ..., a_s)*z_1^a_1*z_2^a2* ... *z_s^a_s, #a_1, ..., a_s=1..infinity), here s is the length of a or z gfmul:=proc(L,a,t,z) local gu,i,s,k0: if not type(a,list) or not type(z,list) or nops(a)<>nops(z) then ERROR(`The second and fourth arguments of gfmul must be lists of same len.`): fi: s:=nops(a): for i from 1 to s do if degree(L,op(i,a))>1 then ERROR(`the first input much be a poly of degree <=1 in`,op(i,a)): fi: od: k0:=L: for i from 1 to s do k0:=subs(op(i,a)=0,k0): od: gu:=t^k0: for i from 1 to s do gu:=gu*t^(coeff(L,op(i,a),1))*op(i,z)/(1-t^(coeff(L,op(i,a),1))*op(i,z)): od: gu: end: #gfpmul(L,a,r,t,z); Given a linear polynomial L(a) in the list of #variables a, and a list of nonegative integers r, and a variable t #and a list of variables z s.t. nops(a)=nops(r)=nops(z), # computes the rational function #sum(a_1^r_1*a_2^r_2*...*a_s^r_s*t^L(a_1, .., a_s)* #z_1^a_1*z_2^a_2* ... *z_s^a_s,a_1, ..., a_s=1..infinity) #where s is the length of a (=nops(r)=nops(z)) gfpmul:=proc(L,a,r,t,z) local i,ri,ri1,s,gu: if not type(a,list) or not type(z,list) or not type(r,list) or nops(a)<>nops(z) or nops(a)<>nops(r) or nops(r)<>nops(z) then ERROR(`The 2nd,3rd and 5th arguments of gfpmul must be lists of same len.`): fi: s:=nops(a): gu:=gfmul(L,a,t,z): for i from 1 to s do ri:=op(i,r): if not type(ri,integer) or ri<0 then ERROR(`All the components of the 3rd argument must be >=0 integers`): fi: if ri>0 then for ri1 from 1 to ri do gu:=op(i,z)*diff(gu,op(i,z)): od: fi: od: normal(gu): end: #`LEGOmul(L,p,a,b,t): `): #`Given lists of variables, a and b s.t. nops(a)=nops(b)=s, a var. t`): #(`and a polynomial p in the variables a_1, ..., a_s, b_1, ..., b_s`): #(`and a linear form L in a_1, ..., a_s`): #(`computes the generating function`): #(` sum(F(a_1, ..., a_s),a_1, ..., a_s=1..infinity)`): #(`where F(a_1, ..., a_s)=t^L(a_1, ..., a_s)+`): #(`t^L(a_1, ..., a_s)*(sum(p(a_1, ..., a_s,b_1, ..., b_s)`): #(`*F(b_1, ..., b_s),b_1, ..., b_s=1..infinity))`): LEGOmul:=proc(L,p,a,b,t) local rj,rj1,j1,i1,s,gu,j,r1,r,G,eq,var,z,i,Z,Pol,mono,gu1,Listr: if not type(a,list) or not type(b,list) or nops(a)<>nops(b) then ERROR(`the 3rd and 4th arguments are lists of the same length`): fi: s:=nops(a): Z:=[seq(z[i],i=1..s)]: gu:=gfmul(L,a,t,Z ): var:={G[seq(0,i1=1..s)]}: Listr:={[seq(0,i1=1..s)]}: Pol:=expand(p): if type(Pol, `+`) then for i from 1 to nops(Pol) do mono:=op(i,Pol): r:=[]: for j from 1 to s do r:=[op(r),degree(mono,op(j,b))]: mono:=normal(mono/op(j,b)^degree(mono,op(j,b))): od: r1:=[]: for j from 1 to s do r1:=[op(r1),degree(mono,op(j,a))]: mono:=normal(mono/op(j,a)^degree(mono,op(j,a))): od: gu:=gu+G[op(r)]*gfpmul(L,a,r1,t,Z)*mono: var:=var union {G[op(r)]}: Listr:=Listr union {r}: od: fi: if not type(Pol, `+`) then mono:=Pol: r:=[]: for j from 1 to s do r:=[op(r),degree(mono,op(j,b))]: mono:=normal(mono/op(j,b)^degree(mono,op(j,b))): od: r1:=[]: for j from 1 to s do r1:=[op(r1),degree(mono,op(j,a))]: mono:=normal(mono/op(j,a)^degree(mono,op(j,a))): od: gu:=gu+G[op(r)]*gfpmul(L,a,r1,t,Z)*mono: var:=var union {G[op(r)]}: Listr:=Listr union {r}: fi: eq:={}: for j from 1 to nops(Listr) do r:=op(j,Listr): gu1:=gu: for j1 from 1 to s do rj:=op(j1,r): for rj1 from 1 to rj do gu1:=z[j1]*diff(gu1,z[j1]): od: od: for i from 1 to s do gu1:=subs(z[i]=1,gu1): od: eq:=eq union {G[op(r)]=gu1}: od: var:=solve(eq,var): subs(var,G[seq(0,i1=1..s)]): normal(%): end: muLEGOmul:=proc(Dims,Ls,pols,a,b,t) local i,j,gu,mu,r,G,eq,var,p,Listr,z, gu1, i1, j1, k1, mono, r1, rj, rj1,i2,i3: var:={}: for i from 1 to nops(Ls) do Listr[i]:={[seq(0,k1=1..op(i,Dims))]}: od: for i from 1 to nops(Ls) do gu[i]:=gfmul(op(i,Ls),[seq(a[i1],i1=1..op(i,Dims))],t, [seq(z[i1],i1=1..op(i,Dims))] ): for j from 1 to nops(Ls) do p:=op(j,op(i,pols)): p:=expand(p): if type(p, `+`) then for i1 from 1 to nops(p) do mono:=op(i1,p): r:=[]: for j1 from 1 to op(j,Dims) do r:=[op(r),degree(mono,b[j1])]: mono:=normal(mono/b[j1]^degree(mono,b[j1])): od: r1:=[]: for j1 from 1 to op(i,Dims) do r1:=[op(r1),degree(mono,a[j1])]: mono:=normal(mono/a[j1]^degree(mono,a[j1])): od: gu[i]:=gu[i]+ G[j,r]*gfpmul(op(i,Ls),[seq(a[i2],i2=1..op(i,Dims))],r1,t, [seq(z[i2],i2=1..op(i,Dims))] )*mono: var:=var union {G[j,r]}: Listr[j]:=Listr[j] union {r}: od: fi: if not type(p, `+`) then mono:=p: r:=[]: for j1 from 1 to op(j,Dims) do r:=[op(r),degree(mono,b[j1])]: mono:=normal(mono/b[j1]^degree(mono,b[j1])): od: r1:=[]: for j1 from 1 to op(i,Dims) do r1:=[op(r1),degree(mono,a[j1])]: mono:=normal(mono/a[j1]^degree(mono,a[j1])): od: gu[i]:=gu[i]+G[j,r]*gfpmul(op(i,Ls),[seq(a[i1],i1=1..op(i,Dims))],r1,t, [seq(z[i1],i1=1..op(i,Dims))] )*mono: var:=var union {G[j,r]}: Listr[j]:=Listr[j] union {r}: fi: od: od: eq:={}: for i from 1 to nops(Ls) do for j from 1 to nops(Listr[i]) do r:=op(j,Listr[i]): gu1:=gu[i]: for j1 from 1 to nops(r) do rj:=op(j1,r): for rj1 from 1 to rj do gu1:=z[j1]*diff(gu1,z[j1]): od: od: for i3 from 1 to op(i,Dims) do gu1:=subs(z[i3]=1,gu1): od: eq:=eq union {G[i,r]=gu1}: od: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: mu:=0: for i from 1 to nops(Ls) do mu:=mu+subs(var,G[i,[seq(0,k1=1..op(i,Dims))]]): od: normal(mu): end: