#DOMINO: A PACKAGE FOR COMPUTING and CONJECTURING Generating Functions #for the Dimer problem and generalizations, Written by Doron #Zeilberger, Temple University <zeilberg@math.temple.edu> #LAST UPDATE: Oct. 10, 1997 #The most current version is available from #http://www.math.temple.edu/~zeilberg #To use it get it first, call it `DOMINO` then go into Maple #by typing `maple`. Then, in Maple , type `read DOMINO;` #and see the help files by doing ezra() and ezra(function_name); #findrec(a,b,f) #program that empirically finds an ordinary linear recurrence with #polynomial coefficients. The input is a sequence f given as a list # and a:=the maximal degree of the coefficients #and b:=the order of the recurrence. The output is the operator in n #and E, where E is the forward unit shift: Ef(n):=f(n+1). #newts(pol,x): gives a list of the power sums of the roots pol(x) #from p_1 to p_(degree of pol(x)), using newtson's equations #newt(pol,x,R): gives a list of the power sums of the roots pol(x) #from p_1 to p_R, using newtson's equations #pol_in_terms_of_p(resh,x) :given the list of power sums [p_1,..,p_deg] #finds the monic polynomial, in x, whose roots are the power sums of #decomp(pol,x): given a polynomial in x, pol, decides whether it is #the algebraic product of a quadratic polynomial pol1, and another pol #pol2. If yes, it returns pol1 and pol2, if not, it returns 0 #decompx(pol,N,x): given a polynomial in the variables N and x #investigates whether, when viewed as a polynomial in N # it can be written as the `alg. product' (i.e. its roots are #the direct product of the sets of the roots of its `factors') #of a quadratic in N #of the form N^2+2*a*x*N-1 where a is an (algebraic) pure number #and another polynomial in N. with coeffs. that are polynomials in x print(`DOMINO: A PACKAGE FOR HANDLING THE DIMER AND MORE GENERAL`): print(`PROBLEMS`): print(`Version of Oct. 10, 1997`): print(`written by Doron Zeilberger(zeilberg@math.temple.edu).`): print(`The most current version is always 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)`): ezra:=proc() if args=NULL then print(`Contains the following procedures:`): print(`algprod, ,decomp,decompx,dimer,DimerBox,dimergf, `): print(`DimerSequence,Finch, findrec`): print(`Ge3D,Gm,Glam, Kast, MonoDimerrect, MonoDimerGe,MonoDimerSequence`): print(`monodimergf, newt,newts,pol_in_terms_of_p`): print(`ratios1,recun,trimergf`): print(``): print(`a specific procedure, type ezra(procedure_name)`): fi: if nops([args])=1 and op(1,[args])=`Finch` then print(`Finch(N): Finds the sequence describing the number of ways`): print(`of tiling a 2n by 2n by 2n box for n=0,1,2, ...,N.`): print(`It also prints out the intermediate values. Warning: It is `): print(`VERY SLOW!. It answers a question of Steve Finch`): print(` to the Domino mailgroup`): fi: if nops([args])=1 and op(1,[args])=`DimerBox` then print(`DimerBox(a,b,c): Number of ways to cover an`): print(`a by b by c box with domino tiles`): fi: if nops([args])=1 and op(1,[args])=`Ge3D` then print(`Ge3D(lam,x,y,z) Given a 3D board where the length in the`): print(`y-axis direction are described, floor by floor,`): print(`in terms of a list of lists lam`): print(`(with the property that lam_i[j]-min(lam)={0,1,2})`): print(`finds the sum of all the weights of domino-tilings`): print(`with the weight x^(tiles in the x direction)*`): print(`y^(tiles in the y direction)*z^(tiles in the z direction)`): fi: if nops([args])=1 and op(1,[args])=`MonoDimerSequence` then print(`MonoDimerSequence(L): finds the sequence whose (n)^term is`): print(` is the number of monomer-dimer tilings of an n by n square`): fi: if nops([args])=1 and op(1,[args])=`DimerSequence` then print(`DimerSequence(L): finds the sequence whose (n)^term is`): print(` is the number of dimer tilings of a 2n by 2n square`): fi: if nops([args])=1 and op(1,[args])=`MonoDimerrect(m,n)` then print(`MonoDimerrect(L): finds the number`): print(` of monomer-dimer tilings of an m by n square`): fi: if nops([args])=1 and op(1,[args])=`MonoDimerGe` then print(`MonoDimerGe(lam,x,y) Given a vector of integers lam, with`): print(`the property that lam_i-min(lam)={0,1,2}, finds the`): print(`number of monomer-dimer tilings of that board`): print(`weighted vy x^(#vertical dimers)*y^(#monomers)`): fi: if nops([args])=1 and op(1,[args])=`Kast` then print(`Kast(n,x): Given an integer n, empirically finds`): print(`Kasteleyn's formula for the number of n by n dimers`): print(`factored over the algebraic field Q(alpha), where`): print(`alpha is a primitive root of unity`): fi: if nops([args])=1 and op(1,[args])=`ratios1` then print(`ratios1(pol,t,L): `): print(`Given a polynomial in t, it determines the multiplicity of`): print(`of the ratios of its roots (if they are all real)`): print(`if all the multiplicities are 1, then the dimer miracle does`): print(`not happen, but if you get 2 and above, then there is hope`): print(`that the dimer miracle would happen`): print(`The precision is 1/10^(L+3)`): fi: if nops([args])=1 and op(1,[args])=`trimergf` then print(`trimergf(m,t,x): `): print(`computes the generating function, a rational function in t,and x, `): print(`sum(trimer(m,n)*t^n,n=0..infinity), where`): print(` trimer(m,n) is the weight enumerator for the tilings by `): print(`trimers of the m by n rectangular board, by the weight`): print(`x^(#of vertical trimers)`): fi: if nops([args])=1 and op(1,[args])=`monodimergf` then print(`monodimergf(m,t,x,y): `): print(`computes the generating function, a rational function in t,x, and y `): print(`sum(monodimer(m,n)*t^n,n=0..infinity), where`): print(` monodimer(m,n) is the weight enumerator for the tilings by monomers`): print(`and dimers of `): print(`the m by n rectangular board, by the weight`): print(`y^(#of monmers)*x^(#of vertical dimers)`): fi: if nops([args])=1 and op(1,[args])=`dimergf` then print(`dimergf(m,t,x): `): print(`computes the generating function, a rational function in t and x `): print(`sum(dimer(m,n)*t^n,n=0..infinity), where`): print(` dimer(m,n) is the weight enumerator for the tilings by dimers of `): print(`the m by n rectangular board, by the weights x^(#of vertical dimers)`): fi: if nops([args])=1 and op(1,[args])=`recun` then print(`recun(m,N,n,x): `): print(`computes the recurrence satisfied by the sequence a_m(n) `): print(` :=dimer(m,n) i.e. the g.f. for the tilings by dimers of `): print(`the m by n rectangular board`): fi: if nops([args])=1 and op(1,[args])=`dimer` then print(`dimer(m,n,x): `): print(`computes the weight enumerator, according to the weight `): print(` x^(nu. of vertical dominos) of the set of all tilings of`): print(`the m by n rectangular board`): fi: if nops([args])=1 and op(1,[args])=`Glam` then print(`Glam(lam,x): given a list lam=[lam_1,...,lam_n]`): print(`computes the weight enumerator, according to the weight `): print(` x^(nu. of vertical dominos) of the set of all tilings of`): print(`the histogram with vertical heights lam_1, ..., lam_n`): fi: if nops([args])=1 and op(1,[args])=`Gm` then print(`Gm(lam,x): given a decreasing list (partition) lam=[lam_1,...,lam_n]`): print(`computes the weight enumerator, according to the weight `): print(` x^(nu. of vertical dominos) of the set of all tilings of`): print(`the histogram with horizontal lengths lam_1, ..., lam_n`): fi: if nops([args])=1 and op(1,[args])=`decompx` then print(`decompx(pol,N,x): given a polynomial in the variables N and x`): print(`investigates whether, when viewed as a polynomial in N`): print(` it can be written as the alg. product (i.e. its roots are`): print(`the direct product of the sets of the roots of its " factors")`): print(`of a quadratic in N`): print(`of the form N^2+2*a*x*N-1 where a is an (algebraic) pure number`): print(`and another polynomial in N. with coeffs. that are polynomials in x`): fi: if nops([args])=1 and op(1,[args])=`decomp` then print(`decomp(pol,x): given a polynomial in x, pol, decides whether it is`): print(`the algebraic product of a quadratic polynomial pol1, and another pol`): print(`pol2. If yes, it returns pol1 and pol2, if not, it returns 0`): fi: if nops([args])=1 and op(1,[args])=`pol_in_terms_of_p` then print(`pol_in_terms_of_p(resh,x) :given the list of power sums [p_1,..,p_deg]`): print(`finds the monic polynomial, in x, whose roots are the power sums of`): fi: if nops([args])=1 and op(1,[args])=`newts` then print(`newts(pol,x): gives a list of the power sums of the roots pol(x)`): print(`from p_1 to p_(degree of pol(x)), using newtson's equations`): fi: if nops([args])=1 and op(1,[args])=`newt` then print(`newt(pol,x,R): gives a list of the power sums of the roots pol(x)`): print(`from p_1 to p_R, using newtson's equations`): fi: if nops([args])=1 and op(1,[args])=`algprod` then print(`algprod(P,Q,x): inputs polynomials P and Q in x , and outputs`): print(`a poly R, in x, of degree deg(P)deg(Q) whose roots are all the`): print(`possible products of the roots of P and Q`): fi: if nops([args])=1 and op(1,[args])=`findrec` then print(`findrec(DEG,ORDER,f,N,n) empir. finds an ordi. linear recurrence with`): print(`polynomial coefficients. The input is a sequence f given as a list`): print(`STARTING at f[1],i.e. f[0] is not considered`): print(` where DEG:=the maximal degree of the coefficients`): print(`and ORDER:=the order of the recurrence. The output is the operator`): print(` in n and N, where N is the forward unit shift: Nf(n):=f(n+1).`): print(`For example findrec(0,2,[1,1,2,3,5,8,13,21,34],N,n) should yield`): print(`N^2-N-1 , and findrec(1,1,[1,2,5,14,42,132,429],N,n) should yield`): print(`(n+2)*N-(4*n+2). If there is not enough data, you will get an`): print(`an error message. If there is no operator, you would get 0`): fi: end: pashet:=proc(p,N) local i,gu1,gu,p1,ra: p1:=normal(p): gu1:=denom(p1): ra:=degree(gu1,N): p1:=subs(n=n+ra,numer(p1)): p1:=expand(p1): gu:=0: for i from 0 to degree(p1,N) do gu:=gu+factor(coeff(p1,N,i))*N^i: od: RETURN(gu): end: #findrec(a,b,f) #program that empirically finds an ordinary linear recurrence with #polynomial coefficients. The input is a sequence f given as a list # and a:=the maximal degree of the coefficients #and b:=the order of the recurrence. The output is the linear #recurrence operator in n #and N, where N is the forward unit shift: Ef(n):=f(n+1). findrec:=proc(DEGREE,ORDER,f,N,n) local ope,var,eq,a,i,j,n0,kv,var1,eq1,mu: if (1+DEGREE)*(1+ORDER)+2+ORDER>nops(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)+2 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 print(`Overdetermined, either DEGREE or ORDER are too high`): print(`The output is not the minimal possible operator`): fi: for i from 1 to nops(kv) do ope:=subs(op(i,kv)=1,ope): od: pashet(ope,N): end: #algprod(P,Q,x): inputs polynomials P and Q in x , and outputs #a poly R, in x, of degree deg(P)deg(Q) whose roots are all the #possible products of the roots of P and Q algprod:=proc(P,Q,x) local degp,degq,deg,pol,var,eq,a,i,mu,lu,reshp,reshq,xp,xq,pol1,i1,i2: degp:=degree(P,x): degq:=degree(Q,x): deg:=degp*degq: reshp:=[]: for i from 0 to degp-1 do reshp:=[op(reshp),x^i]: od: mu:=P/coeff(P,x,degp): mu:=x^degp-mu: reshp:=[op(reshp),mu]: lu:=mu: for i from 1 to deg-degp do lu:=expand(x*lu): lu:=subs(x^(degp)=mu,lu): lu:=expand(lu): reshp:=[op(reshp),lu]: od: reshq:=[]: for i from 0 to degq-1 do reshq:=[op(reshq),x^i]: od: mu:=Q/coeff(Q,x,degq): mu:=x^degq-mu: reshq:=[op(reshq),mu]: lu:=mu: for i from 1 to deg-degq do lu:=expand(x*lu): lu:=subs(x^(degq)=mu,lu): lu:=expand(lu): reshq:=[op(reshq),lu]: od: reshp:=subs(x=xp,reshp): reshq:=subs(x=xq,reshq): pol:=0: pol1:=0: eq:={}: var:={}: for i from 0 to deg do pol:=pol+a[i]*x^i: var:=var union {a[i]}: pol1:=expand(pol1+a[i]*op(i+1,reshp)*op(i+1,reshq)): od: for i1 from 0 to degp-1 do for i2 from 0 to degq-1 do eq:=eq union {coeff(coeff(pol1,xp,i1),xq,i2)}: od: od: var:=solve(eq,var): pol:=subs(var,pol): normal(pol/coeff(pol,x,degree(pol,x))): end: newts:=proc(pol,x) local mu,deg,resh,pol1,r,i: deg:=degree(pol,x): pol1:=pol/coeff(pol,x,deg): pol1:=expand(pol1): resh:=[-coeff(pol1,x,deg-1)]: for r from 2 to deg do mu:=0: for i from 1 to r-1 do mu:=mu-coeff(pol1,x,deg-i)*op(r-i,resh): od: mu:=mu-r*coeff(pol1,x,deg-r): mu:=expand(mu): resh:=[op(resh),mu]: od: resh: end: newt:=proc(pol,x,R) local mu,deg,resh,pol1,r,i: if R<=degree(pol,x) then RETURN([op(1..R,newts(pol,x))]): fi: deg:=degree(pol,x): pol1:=pol/coeff(pol,x,deg): pol1:=expand(pol1): resh:=newts(pol,x): for r from deg+1 to R do mu:=0: for i from 1 to deg do mu:=mu-coeff(pol1,x,deg-i)*op(r-i,resh): od: mu:=expand(mu): resh:=[op(resh),mu]: od: resh: end: #pol_in_terms_of_p(resh) :given the list of power sums [p_1,..,p_deg] #finds the monic polynomial whose roots are the power sums of pol_in_terms_of_p:=proc(resh,x) local pol,mu,deg,r,i: deg:=nops(resh): pol:=x^deg: for r from 1 to deg do mu:=-op(r,resh): for i from 1 to r-1 do mu:=mu-coeff(pol,x,deg-i)*op(r-i,resh): od: mu:=mu/r: mu:=expand(mu): pol:=pol+mu*x^(deg-r): od: pol: end: decomp:=proc(pol,x) local i,deg,resh1,resh2,resh,a,pol1,pol2,eq,var,resh2a,eq1: deg:=degree(pol,x): deg:=deg/2: if not type(deg, integer) then RETURN(0): fi: pol1:=x^2+2*a*x-1: var:={a}: resh1:=newt(pol1,x,2*deg): resh:=newts(pol,x): resh2:=[]: for i from 1 to 2*deg do resh2:=[op(resh2),op(i,resh)/op(i,resh1)]: od: pol2:=pol_in_terms_of_p([op(1..deg,resh2)],x): resh2a:=newt(pol2,x,2*deg): eq:={}: for i from 1 to deg do eq1:=normal(op(i,resh2)-op(i,resh2a)): if eq1<>0 then ERROR(`Something is wrong`): fi: od: for i from deg+1 to 2*deg do eq1:=normal(op(i,resh2)-op(i,resh2a)): eq:=eq union {eq1}: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: pol1:=subs(var,pol1): pol2:=subs(var,pol2): pol1,pol2: end: decompx:=proc(pol,N,x) local i1,i,mu,deg,resh1,resh2,resh,a,pol1,pol2,eq,var,resh2a,eq1: deg:=degree(pol,N): deg:=deg/2: if not type(deg, integer) then RETURN(0): fi: pol1:=N^2+2*a*x*N-1: var:={a}: resh1:=newt(pol1,N,2*deg): resh:=newts(pol,N): resh2:=[]: for i from 1 to 2*deg do resh2:=[op(resh2), rem(op(i,resh),op(i,resh1),x)]: od: eq:={}: for i from 1 to 3 do mu:=op(i,resh2): mu:=expand(numer(normal(mu))): for i1 from 0 to degree(mu,x) do eq:=eq union {numer(normal(coeff(mu,x,i1)))=0}: od: od: var:=solve(eq,var): if var=NULL then RETURN(0): fi: var: end: Glam:=proc(lam,x) local i,i1,mu: option remember: for i from 1 to nops(lam) do if op(i,lam)<0 then RETURN(0): fi: od: if nops(lam)=0 then RETURN(1): fi: if op(1,lam)=0 then RETURN(Glam( [op(2..nops(lam),lam)],x)): fi: for i from 1 to nops(lam)-1 while op(i,lam)<op(i+1,lam) do od: i1:=i: mu:=0: if i1<nops(lam) and op(i1,lam)=op(i1+1,lam) then mu:=x* Glam([op(1..i1-1,lam),op(i1,lam)-1,op(i1+1,lam)-1,op(i1+2..nops(lam),lam)],x): fi: mu:=mu+ Glam([op(1..i1-1,lam),op(i1,lam)-2,op(i1+1..nops(lam),lam)],x): RETURN(expand(mu)): end: dimer:=proc(m,n,x) local i,lam: lam:=[]: for i from 1 to m do lam:=[op(lam),n]: od: Ge(lam,x): end: #recun finds a liner recurrence satisfied by the sequnce a_m(n):=dimer(m,n,x) #recurrence w.r.t. n recun:=proc(m,N,n,x) local i,ope,resh,ORDER: ope:=0: resh:=[]: for i from 1 to 6 do resh:=[op(resh),dimer(m,i,x)]: od: for ORDER from 0 while ope=0 do ope:=findrec(0,ORDER,resh,N,n): resh:=[op(resh),dimer(m,2*ORDER+7,x), dimer(m,2*ORDER+8,x)]: od: ope: end: checkdec:=proc(lam) local i: for i from 1 to nops(lam)-1 do if op(i+1,lam)>op(i,lam) then RETURN(0): fi: od: if op(nops(lam),lam)<=0 then RETURN(0): fi: 1: end: first_big_jump:=proc(lam) local i: for i from 1 to nops(lam)-1 do if op(i,lam)-op(i+1,lam)>=2 then RETURN(i): fi: od: if op(nops(lam),lam)>=2 then RETURN(nops(lam)): fi: nops(lam)+1: end: first_big_steady_jump:=proc(lam) local i: if nops(lam)=1 then if op(1,lam)>=2 then RETURN(1): else RETURN(2): fi: fi: if op(1,lam)-op(2,lam)>=2 then RETURN(1): fi: for i from 2 to nops(lam)-1 do if op(i,lam)-op(i+1,lam)>=2 and op(i-1,lam)=op(i,lam) then RETURN(i): fi: od: if op(nops(lam),lam)>=2 and op(nops(lam)-1,lam)=op(nops(lam),lam) then RETURN(nops(lam)): fi: nops(lam)+1: end: beg_staircase:=proc(lam,i) local j: for j from i-1 by -1 to 1 do if op(j+1,lam)=op(j,lam) then RETURN(j+1): fi: od: 1: end: chop:=proc(lam) local i,lamd: lamd:=[]: for i from 1 to nops(lam) do if op(i,lam)>0 then lamd:=[op(lamd),op(i,lam)]: fi: od: lamd: end: Gm:=proc(lam,x) local r,i,j,lam1,n,lam2: option remember: n:=nops(lam): if nops(lam)=0 then RETURN(1): fi: i:=first_big_steady_jump(lam): if i<=n then if i=n then lam1:=[op(1..n-1,lam),op(n,lam)-2]: fi: if i>1 then lam1:=[op(1..i-1,lam),op(i,lam)-2,op(i+1..n,lam)]: lam2:=[op(1..i-2,lam),op(i-1,lam)-1, op(i,lam)-1, op(i+1..n,lam)]: lam1:=chop(lam1): lam2:=chop(lam2): RETURN(expand(Gm(lam1,x)+x*Gm(lam2,x) )): else lam1:=[op(i,lam)-2,op(i+1..n,lam)]: lam1:=chop(lam1): RETURN(expand(Gm(lam1,x) )): fi: fi: i:=first_big_jump(lam): if i<=n then lam1:=[op(1..i-1,lam),op(i,lam)-2,op(i+1..n,lam)]: j:=beg_staircase(lam,i): if j=1 then lam1:=chop(lam1): RETURN(Gm(lam1,x)): fi: lam2:=[op(1..j-2,lam),op(j-1,lam)-1]: for r from j to i-1 do lam2:=[op(lam2),op(r,lam)-2]: od: lam2:=[op(lam2),op(i,lam)-1,op(i+1..n,lam)]: lam1:=chop(lam1): lam2:=chop(lam2): RETURN(expand(Gm(lam1,x)+x^(i-j+1)*Gm(lam2,x) )): fi: j:=beg_staircase(lam,n): if j=1 then RETURN(0): fi: lam2:=[op(1..j-2,lam),op(j-1,lam)-1]: for r from j to n-1 do lam2:=[op(lam2),op(r,lam)-2]: od: lam2:=[op(lam2),op(n,lam)-1]: lam2:=chop(lam2): RETURN(x^(n-j+1)*Gm(lam2,x) ): end: #Ge(lam) like Glam, but now the lambda are only allowed to have #the property that lam_i-min(lam)={0,1,2} Ge:=proc(lam,x) local i,n,mini,bet,maxi,lam1,lam2: option remember: n:=nops(lam): if nops(lam)=0 then RETURN(1): fi: for i from 1 to n while op(i,lam)=0 do od: if i=n+1 then RETURN(1): fi: mini:=min(op(lam)): if mini<0 then RETURN(0): fi: bet:=[]: for i from 1 to n do bet:=[op(bet),op(i,lam)-mini]: od: maxi:=max(op(bet)): if maxi>2 then ERROR(`Wrong input`): fi: if maxi=2 then for i from 1 to n while op(i,bet)<2 do od: lam1:=[op(1..i-1,lam),op(i,lam)-2,op(i+1..n,lam)]: if i=n or op(i+1,bet)<2 then RETURN(Ge(lam1,x)): else lam2:=[op(1..i-1,lam),op(i,lam)-1, op(i+1,lam)-1,op(i+2..n,lam)]: RETURN(expand(Ge(lam1,x)+x*Ge(lam2,x))): fi: fi: if maxi=1 then for i from 1 to n while op(i,bet)<1 do od: lam1:=[op(1..i-1,lam),op(i,lam)-2,op(i+1..n,lam)]: if i=n or op(i+1,bet)<1 then RETURN(Ge(lam1,x)): else lam2:=[op(1..i-1,lam),op(i,lam)-1, op(i+1,lam)-1,op(i+2..n,lam)]: RETURN(expand(Ge(lam1,x)+x*Ge(lam2,x))): fi: fi: if maxi=0 then lam1:=[op(1..n-1,lam),op(n,lam)-2]: if nops(lam)=1 then RETURN(Ge(lam1,x)): fi: lam2:=[op(1..n-2,lam),op(n-1,lam)-1,op(n,lam)-1]: RETURN(expand(Ge(lam1,x)+x*Ge(lam2,x))): fi: end: redu:=proc(resh) local i,katan,resh1: katan:=max(op(resh)): if katan>0 or katan<-2 then ERROR(`Something is wrong`): fi: resh1:=[]: for i from 1 to nops(resh) do resh1:=[op(resh1),op(i,resh)-katan]: od: resh1,-katan: end: equ:=proc(resh,t,x,A) local i,eq,var,resh1,resh2,n,gu: n:=nops(resh): for i from 1 to nops(resh) while op(i,resh)<>0 do od: resh1:=[op(1..i-1,resh),op(i,resh)-2,op(i+1..n,resh)]: gu:=redu(resh1): resh1:=gu[1]: eq:=A[op(resh)]-A[op(resh1)]*t^gu[2]: var:={A[op(resh1)]}: if i<n and op(i+1,resh)=0 then resh2:=[op(1..i-1,resh),op(i,resh)-1,op(i+1,resh)-1,op(i+2..n,resh)]: gu:=redu(resh2): resh2:=gu[1]: eq:=eq-A[op(resh2)]*t^gu[2]*x: var:=var union {A[op(resh2)]}: fi: if min(op(resh))=0 then eq:=eq-1: fi: eq,var: end: dimergf:=proc(m,t,x) local eq,var,i,RESHeq,resh,A,gu,mu: eq:={}: resh:=[seq(0,i=1..m)]: RESHeq:={A[op(resh)]}: var:={A[op(resh)]}: gu:=equ(resh,t,x,A): eq:={gu[1]}: var:=var union gu[2]: while var minus RESHeq<>{} do mu:=op(1,var minus RESHeq): gu:=equ([op(mu)],t,x,A): eq:=eq union {gu[1]}: var:=var union gu[2]: RESHeq:=RESHeq union {mu}: od: var:=solve(eq,var): gu:=subs(var,A[seq(0,i=1..m)]): normal(gu): end: equdm:=proc(resh,t,x,y,A) local i,eq,var,resh0,resh1,resh2,n,gu: n:=nops(resh): for i from 1 to nops(resh) while op(i,resh)<>0 do od: resh0:=[op(1..i-1,resh),op(i,resh)-1,op(i+1..n,resh)]: resh1:=[op(1..i-1,resh),op(i,resh)-2,op(i+1..n,resh)]: gu:=redu(resh1): resh1:=gu[1]: eq:=A[op(resh)]-A[op(resh1)]*t^gu[2]: var:={A[op(resh1)]}: gu:=redu(resh0): resh0:=gu[1]: eq:=eq-A[op(resh0)]*t^gu[2]*y: var:=var union {A[op(resh0)]}: if i<n and op(i+1,resh)=0 then resh2:=[op(1..i-1,resh),op(i,resh)-1,op(i+1,resh)-1,op(i+2..n,resh)]: gu:=redu(resh2): resh2:=gu[1]: eq:=eq-A[op(resh2)]*t^gu[2]*x: var:=var union {A[op(resh2)]}: fi: if min(op(resh))=0 then eq:=eq-1: fi: eq,var: end: #monodimergf(m,t,x,y): computes the gf sum(monodim(m,n)*t^n,n=0..infty) #where monodim(m,n) is the wt. enumerator according to the wt #x^(no. of vertical dimers)*y^(no. of monomers) of the tilings by #dimers and monomers of an m by n rect. board monodimergf:=proc(m,t,x,y) local eq,var,i,RESHeq,resh,A,gu,mu: eq:={}: resh:=[seq(0,i=1..m)]: RESHeq:={A[op(resh)]}: var:={A[op(resh)]}: gu:=equdm(resh,t,x,y,A): eq:={gu[1]}: var:=var union gu[2]: while var minus RESHeq<>{} do mu:=op(1,var minus RESHeq): gu:=equdm([op(mu)],t,x,y,A): eq:=eq union {gu[1]}: var:=var union gu[2]: RESHeq:=RESHeq union {mu}: od: var:=solve(eq,var): gu:=subs(var,A[seq(0,i=1..m)]): normal(gu): end: redu3:=proc(resh) local i,katan,resh1: katan:=max(op(resh)): if katan>0 or katan<-3 then ERROR(`Something is wrong`): fi: resh1:=[]: for i from 1 to nops(resh) do resh1:=[op(resh1),op(i,resh)-katan]: od: resh1,-katan: end: equ3:=proc(resh,t,x,A) local i,eq,var,resh1,resh2,n,gu: n:=nops(resh): for i from 1 to nops(resh) while op(i,resh)<>0 do od: resh1:=[op(1..i-1,resh),op(i,resh)-3,op(i+1..n,resh)]: gu:=redu3(resh1): resh1:=gu[1]: eq:=A[op(resh)]-A[op(resh1)]*t^gu[2]: var:={A[op(resh1)]}: if i<n-1 and op(i+1,resh)=0 and op(i+2,resh)=0 then resh2:=[op(1..i-1,resh),op(i,resh)-1,op(i+1,resh)-1, op(i+2,resh)-1,op(i+3..n,resh)]: gu:=redu3(resh2): resh2:=gu[1]: eq:=eq-A[op(resh2)]*t^gu[2]*x: var:=var union {A[op(resh2)]}: fi: if min(op(resh))=0 then eq:=eq-1: fi: eq,var: end: trimergf:=proc(m,t,x) local eq,var,i,RESHeq,resh,A,gu,mu: eq:={}: resh:=[seq(0,i=1..m)]: RESHeq:={A[op(resh)]}: var:={A[op(resh)]}: gu:=equ3(resh,t,x,A): eq:={gu[1]}: var:=var union gu[2]: while var minus RESHeq<>{} do mu:=op(1,var minus RESHeq): gu:=equ3([op(mu)],t,x,A): eq:=eq union {gu[1]}: var:=var union gu[2]: RESHeq:=RESHeq union {mu}: od: var:=solve(eq,var): gu:=subs(var,A[seq(0,i=1..m)]): normal(gu): end: kats:=proc(mis,L): trunc(mis*10^L)/10^L: end: ratios1:=proc(pol,t,L) local ku,gu,mu,i,j,pol1,lu,ru: readlib(realroot): ku:=realroot(pol,1/10^(L+3)): if nops(ku)<degree(pol,t) then ERROR(pol, `has complex roots`): fi: mu:=[]: for i from 1 to nops(ku) do lu:=op(i,ku): mu:=[op(mu),evalf(op(1,lu)/2+op(2,lu)/2)]: od: gu:=[]: for i from 1 to nops(mu) do for j from 1 to nops(mu) do if i<>j then gu:=[op(gu),kats(op(i,mu)/op(j,mu),L-2)]: fi: od: od: pol1:=0: for i from 1 to nops(gu) do pol1:=pol1+t^op(i,gu): od: mu:=[]: for i from 1 to nops(pol1) do ru:=op(1,op(i,pol1)): if ru=t then ru:=1: fi: mu:=[op(mu), ru]: od: sort(mu): end: #Kast(n,x): Given an integer n, empirically finds #Kasteleyn's formula for the number of n by n dimers #factored over the algebraic field Q(alpha), where #alpha is a primitive root of unity Kast:=proc(n,x) local p,Z,t: with(numtheory): p:=cyclotomic(n+1,t): Z:=RootOf(p): factor(dimer(n,n,x),Z); end: #MonoDimerGe(lam,x,y) Given a vector of integers lam, with #the property that lam_i-min(lam)={0,1,2}, finds the #number of monomer-dimer tilings of that board #weighted vy x^(#vertical dimers)*y^(#monomers) MonoDimerGe:=proc(lam,x,y) local i,n,mini,bet,maxi,lam1,lam2,lam0: option remember: n:=nops(lam): if nops(lam)=0 then RETURN(1): fi: for i from 1 to n while op(i,lam)=0 do od: if i=n+1 then RETURN(1): fi: mini:=min(op(lam)): if mini<0 then RETURN(0): fi: bet:=[]: for i from 1 to n do bet:=[op(bet),op(i,lam)-mini]: od: maxi:=max(op(bet)): if maxi>2 then ERROR(`Wrong input`): fi: if maxi=2 then for i from 1 to n while op(i,bet)<2 do od: lam1:=[op(1..i-1,lam),op(i,lam)-2,op(i+1..n,lam)]: lam0:=[op(1..i-1,lam),op(i,lam)-1,op(i+1..n,lam)]: if i=n or op(i+1,bet)<2 then RETURN(expand(MonoDimerGe(lam1,x,y)+y*MonoDimerGe(lam0,x,y))): else lam2:=[op(1..i-1,lam),op(i,lam)-1, op(i+1,lam)-1,op(i+2..n,lam)]: RETURN(expand(MonoDimerGe(lam1,x,y)+x*MonoDimerGe(lam2,x,y)+y*MonoDimerGe(lam0,x,y) )): fi: fi: if maxi=1 then for i from 1 to n while op(i,bet)<1 do od: lam0:=[op(1..i-1,lam),op(i,lam)-1,op(i+1..n,lam)]: lam1:=[op(1..i-1,lam),op(i,lam)-2,op(i+1..n,lam)]: if i=n or op(i+1,bet)<1 then RETURN(expand(MonoDimerGe(lam1,x,y)+y*MonoDimerGe(lam0,x,y))): else lam2:=[op(1..i-1,lam),op(i,lam)-1, op(i+1,lam)-1,op(i+2..n,lam)]: RETURN(expand(MonoDimerGe(lam1,x,y)+x*MonoDimerGe(lam2,x,y)+y*MonoDimerGe(lam0,x,y) )): fi: fi: if maxi=0 then lam0:=[op(1..n-1,lam),op(n,lam)-1]: lam1:=[op(1..n-1,lam),op(n,lam)-2]: if nops(lam)=1 then RETURN(expand(MonoDimerGe(lam1,x,y)+y*MonoDimerGe(lam0,x,y))): fi: lam2:=[op(1..n-2,lam),op(n-1,lam)-1,op(n,lam)-1]: RETURN(expand(MonoDimerGe(lam1,x,y)+x*MonoDimerGe(lam2,x,y)+y*MonoDimerGe(lam0,x,y) )): fi: end: MonoDimerrect:=proc(m,n) local i: MonoDimerGe([seq(m,i=1..n)],1,1): end: MonoDimerSequence:=proc(L) local gu,i: gu:=[]: for i from 1 to L do gu:=[op(gu),MonoDimerrect(i,i)]: #print(gu): od: gu: end: DimerSequence:=proc(L) local gu,i: gu:=[]: for i from 0 to L do gu:=[op(gu),dimer(2*i,2*i,1)]: #print(gu): od: gu: end: ez:=proc():print(`Uptsiot(lam,x,y,z,j,i), DimerBox(a,b,c)`):end: #Uptsiot(lam,x,y,z,j,i) Given a list of lists lam #and a row (j,i) (the ith item in the jth floor) #signalling the first place where it deviates by #2 from the minimum, lists the resulting #set of smaller lamdas, prefixed by the direction Uptsiot:=proc(lam,x,y,z,j,i) local lamx,lamy,lamz,k,n,gu,lamj,lamj1,lamjn,lamj1n: k:=nops(lam): if k=0 then RETURN({}): fi: n:=nops(op(1,lam)): gu:={}: lamj:=op(j,lam): if op(i,lamj)>=2 then lamjn:=[op(1..i-1,lamj),op(i,lamj)-2,op(i+1..nops(lamj),lamj)]: lamy:=[op(1..j-1,lam),lamjn,op(j+1..nops(lam),lam)]: gu:=gu union {[y,lamy]}: fi: if i<n then if op(i+1,lamj)>=op(i,lamj) then lamjn:=[op(1..i-1,lamj),op(i,lamj)-1,op(i+1,lamj)-1,op(i+2..nops(lamj),lamj)]: lamx:=[op(1..j-1,lam),lamjn, op(j+1..k,lam)]: gu:=gu union {[x,lamx]}: fi: fi: if j<k then lamj1:=op(j+1,lam): if op(i,lamj1)>=op(i,lamj) then lamjn:=[op(1..i-1,lamj),op(i,lamj)-1,op(i+1..nops(lamj),lamj)]: lamj1n:=[op(1..i-1,lamj1),op(i,lamj1)-1,op(i+1..nops(lamj1),lamj1)]: lamz:=[op(1..j-1,lam),lamjn, lamj1n,op(j+2..k,lam)]: gu:=gu union {[z,lamz]}: fi: fi: gu: end: #Gek(lam,x,y) like Glam, but now the lambda are only allowed to have #the property that lam_i-min(lam)={0,1,2} #with the weight x^(tiles in the x direction)* #y^(tiles in the y direction) Gek:=proc(lam,x,y) local i,n,mini,bet,maxi,lam1,lam2: option remember: n:=nops(lam): if nops(lam)=0 then RETURN(1): fi: for i from 1 to n while op(i,lam)=0 do od: if i=n+1 then RETURN(1): fi: mini:=min(op(lam)): if mini<0 then RETURN(0): fi: bet:=[]: for i from 1 to n do bet:=[op(bet),op(i,lam)-mini]: od: maxi:=max(op(bet)): if maxi>2 then ERROR(`Wrong input`): fi: if maxi=2 then for i from 1 to n while op(i,bet)<2 do od: lam1:=[op(1..i-1,lam),op(i,lam)-2,op(i+1..n,lam)]: if i=n or op(i+1,bet)<2 then RETURN(y*Gek(lam1,x,y)): else lam2:=[op(1..i-1,lam),op(i,lam)-1, op(i+1,lam)-1,op(i+2..n,lam)]: RETURN(expand(y*Gek(lam1,x,y)+x*Gek(lam2,x,y))): fi: fi: if maxi=1 then for i from 1 to n while op(i,bet)<1 do od: lam1:=[op(1..i-1,lam),op(i,lam)-2,op(i+1..n,lam)]: if i=n or op(i+1,bet)<1 then RETURN(y*Gek(lam1,x,y)): else lam2:=[op(1..i-1,lam),op(i,lam)-1, op(i+1,lam)-1,op(i+2..n,lam)]: RETURN(expand(y*Gek(lam1,x,y)+x*Gek(lam2,x,y))): fi: fi: if maxi=0 then lam1:=[op(1..n-1,lam),op(n,lam)-2]: if nops(lam)=1 then RETURN(y*Gek(lam1,x,y)): fi: lam2:=[op(1..n-2,lam),op(n-1,lam)-1,op(n,lam)-1]: RETURN(expand(y*Gek(lam1,x,y)+x*Gek(lam2,x,y))): fi: end: Metsa:=proc(bet,maxi) local n,k,i,j: k:=nops(bet): n:=nops(op(1,bet)): for j from 1 to k do for i from 1 to n do if op(i,op(j,bet))=maxi then RETURN(j,i): fi: od: od: end: #Ge3D(lam,x,y,z) Given a 3D board where the length in the #y-axis direction are described, floor by floor, #in terms of a list of lists lam #(with the property that lam_i[j]-min(lam)={0,1,2}) #finds the sum of all the weights of domino-tilings #with the weight x^(tiles in the x direction)* #y^(tiles in the y direction)*z^(tiles in the z direction) Ge3D:=proc(lam,x,y,z) local i,n,mini,bet,maxi,kak,alph,r1,mu,j,k,gu,ku,i1,j1,lu,k1: option remember: if nops(lam)=0 then RETURN(1): fi: k:=nops(lam): n:=nops(op(1,lam)): for i from 2 to k do if nops(op(i,lam))<>n then ERROR(`All members of the first input must be lists of the same length`): fi: od: ku:=0: for j from 1 to k do for i from 1 to n do ku:=ku+op(i,op(j,lam)): od: od: if ku=0 then RETURN(1): fi: mini:=min(op(op(1,lam))): for k1 from 2 to nops(lam) do kak:=min(op(op(k1,lam))): if kak<mini then mini:=kak: fi: od: if mini<0 then RETURN(0): fi: bet:=[]: for j from 1 to k do alph:=[]: for i from 1 to n do alph:=[op(alph),op(i,op(j,lam))-mini]: od: bet:=[op(bet),alph]: od: maxi:=max(op(op(1,bet))): for k1 from 2 to nops(bet) do kak:=max(op(op(k1,bet))): if kak>maxi then maxi:=kak: fi: od: if maxi>2 then ERROR(`Wrong input`): fi: if maxi=2 then lu:=Metsa(bet,maxi): j1:=lu[1]: i1:=lu[2]: gu:=Uptsiot(lam,x,y,z,j1,i1): mu:=0: for r1 from 1 to nops(gu) do mu:=mu+op(r1,gu)[1]*Ge3D(op(r1,gu)[2] ,x,y,z): od: RETURN(expand(mu)): fi: if maxi=1 then lu:=Metsa(bet,maxi): j1:=lu[1]: i1:=lu[2]: gu:=Uptsiot(lam,x,y,z,j1,i1): mu:=0: for r1 from 1 to nops(gu) do mu:=mu+op(r1,gu)[1]*Ge3D(op(r1,gu)[2] ,x,y,z): od: RETURN(expand(mu)): fi: if maxi=0 then gu:=Uptsiot(lam,x,y,z,1,1): mu:=0: for r1 from 1 to nops(gu) do mu:=mu+op(r1,gu)[1]*Ge3D(op(r1,gu)[2] ,x,y,z): od: RETURN(expand(mu)): fi: end: #DimerBox(a,b,c): Number of ways to cover an #a by b by c box with domino tiles DimerBox:=proc(a,b,c) local lam,i,j: lam:=[seq([seq(a,i=1..b)],j=1..c)]: Ge3D(lam,1,1,1): end: #Finch(N): Finds the sequence describing the number of ways #of tiling a 2n by 2n by 2n box for n=0,1,2, ...,N #It also prints out the intermediate answers Finch:=proc(N) local i,lu: lu:=[]: for i from 1 to N do lu:=[op(lu), DimerBox(2*i,2*i,2*i)]: print(lu): od: lu: end: