####################################################################### ## GoodDyson: Save this file as GoodDyson. To use it, stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read GoodDyson: # ## Then follow the instructions given there # ## # ## Written by Drew Sills and Doron Zeilberger # ## (Rutgers University) # ## asills at math dot rutgers dot edu, # ## zeilberg at math dot rutgers dot edu. # ####################################################################### ## First draft: 19 Feb 2004 print(`Tested on Maple 8 and 9`): #print(`This is the version of 18 May 2005.`): print(`This is the version of 2 November 2005.`): print(` `): print(`This is GoodDyson, a Maple package`): print(`that conjectures formulas inspired by variations of`): print(`Freeman Dyson's 1962 constant term conjecture`); print(`and automatically supplies proofs inspired by I. J. Good's 1970 proof`): print(`of Dyson's conjecture.`); print(``): print(`It accompanies the article : Disturbing the Dyson Conjecture`): print(`(In a GOOD Way) `): print(`by Drew Sills and Doron Zeilberger`): print(` `): print(`The most current version is available on WWW at:`): print(` http://www.math.rutgers.edu/~asills/papers.html`): print(`Please report all bugs to: asills at math dot rutgers dot edu and/or`); print(`zeilberg at math dot rutgers dot edu .`): print(): print(`For general help, and a list of the available functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name)" `): ezra:=proc() if args=NULL then print(`This is GoodDyson`): print(`a Maple package`): print(): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`The main procedures are`); print(`DysonExpansion, GetDysonCoeff, GuessRat, GuessDysonCoeff, TurboDyson and WritePaper.`): print(`The supporting procedures are: DysonProd, GenPol, GenRat, Boundary.`): print(): elif nargs=1 and args[1]=GenPol then print(`GenPol(Resh,a,deg): Given a list of variables Resh, a symbol a,`): print(`and a non-negative integer deg, outputs a generic polynomial`): print(`of total degree deg,in the variables Resh, indexed by the letter`): print(` followed by the set of coeffs. For example, try `): print(` GenPol([x,y],a,1); `): elif nargs=1 and args[1]=GenRat then print(`GenRat(Resh,a,b,degT,degB): The generic form of a rational function`): print(`in the list of variables Resh indexed by a on top and b at bottom`): print(`of total degree degT on top and degB at bottom`): print(`for example, try: GenRat([x,y],a,b,1,2)`): elif nargs=1 and args[1]=GuessRat then print(`GuessRat(F,Resh,deg,m): Given a function or procedure F with`): print(`nops(Resh) arguments (that are integers), guesses a rational`): print(`function, using data that starts at m`): print(`with the sum of total degrees of top and bottom <=deg`): print(`that seems fits it, or returns FAIL`): print(`For example, try: GuessRat(a->a/(a+1),[a],2,1);`): elif nargs=1 and args[1]=DysonExpansion then print(`DysonExpansion(n,C) prints the partial expansion of the`); print(`n-dimensional Dyson product including all terms of complexity `): print(`at most C, where the complexity of a term is the sum of the`); print(`positive exponents in the monomial.`); print(`For example,`); print(`> DysonExpansion(2,3)`); DysonExpansion(2,3); elif nargs=1 and args[1]=GuessDysonCoeff then print(`GuessDysonCoeff(b1,b2,b3,...,bn):`); print(`Given that b1,b2,b3,...,bn are n specific integers,`); print(`this procedure attempts to find a rational function`); print(` R = R(a1,a2,a3,...,an) such that`): print(`the coefficient of x[1]^b1 * x[2]^b2 * .... * x[n]^bn in`); print(`the Dyson product`); #print( Product(Product( (1- x[i]/x[j])^a[i] * (1-x[j]/x[i])^a[j], # j=i+1..n),i=1..n)); print( Product(Product( (1- x[i]/x[j])^a[j] * (1-x[j]/x[i])^a[i], j=i+1..n),i=1..n)); print(` is R * (a1+a2+ ... + an)!/ (a1! a2! ... an!).`); print(`If no such R exists, the procedure returns 'FAIL'.`); elif nargs=1 and args[1]=DysonProd then print(`DysonProd(x,a,b): For x and a symbolic or a list of symbols,`); print(`and b a list such that n=nops(b),`); print(`DysonProd returns the n-dimensional Dyson product`); #print( # Product( x[i]^(-b[i]), i=1..n) * # Product(Product( (1-x[i]/x[j])^a[i] * (1-x[j]/x[i])^a[j],j=i+1..n),i=1..n)); print( Product( x[i]^(-b[i]), i=1..n) * Product(Product( (1-x[i]/x[j])^a[j] * (1-x[j]/x[i])^a[i],j=i+1..n),i=1..n)); elif nargs=1 and args[1]=Boundary then print(`Boundary(x,a,b): `); print(`For x and a (lists of symbols) and b (a list of n integers),`); print(`this procedure returns a list LP of n Laurent polynomials`); print(`where the k-th item on the list LP is the coefficient`); print(`which arises as a multiple of the (n-1)-dimensional Dyson product`); print(`in the series expansion of the a[k]=0 case of the`); print(`full n-dimensional Dyson product about the variable x[k]`); elif nargs=1 and args[1]=WritePaper then print(`WritePaper(b1,b2,...,bn) states and proves a theorem which gives`); print(`an explicit expression, as a funtion of a1,a2,...,an,`); print(`for the coefficient of x1^b1 * x2^b2 * ... * xn^bn`); print(`in the expansion of the Dyson product`) elif nargs=1 and args[1]=TurboDyson then print(`TurboDyson(n,C) calculates and stores for future use all`); print(`coefficients of x1^b1 * x2^b2 * ... * xn^bn`); print(`in the expansion of the Dyson product`); print(`which have complexity less than or equal to C`); print(`The complexity C of (b1,b2,b3,...,bn) is the sum of the`); print(`positive components of the vector (b1,b2,...,bn)`); elif nargs=1 and args[1]=GetDysonCoeff then print(`GetDysonCoeff(b1,b2,...,bn) returns the same output as`); print(`GuessDysonCoeff...the only difference is that GetDysonCoeff`); print(`first checks to see if the calculation has already been done`); print(`and retrieves this prior result when available. If not, then`); print(`the calculation is done and stored.`); else print(`There is no Ezra for`, args): fi: end: Ezra:=proc() ezra() end: Help:=proc() ezra() end: rf:=proc(a,n) local i; if type(n,integer) then product(a+i,i=0..n-1) else 'rf'(a,n) fi end: FallFact:=proc(a,n) local i; if type(n,integer) then product(a-i,i=0..n-1) else 'FallFact'(a,n) fi end: DysonProd:=proc(x,a,b) local i,j,n; n:=nops(b); product(x[i]^(-b[i]), i=1..n)* product(product( (1-x[i]/x[j])^a[j] *(1-x[j]/x[i])^a[i], j=i+1..n),i=1..n-1) end: WritePaper:=proc() local b,x,a,n,i,j,nextlist,xpwrs; b:=[args]; x:=[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20]; a:=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20]; n:=nops(b); xpwrs:=product(x[i]^b[i],i=1..n); if add(b[i],i=1..n) <> 0 then if xpwrs=1 then print(`The constant term`) else print(`The coefficient of`, xpwrs) fi; print(`in the expansion of`); print(DysonProd(x,a,[seq(0,i=1..n)])); print( `is trivially zero. `); return() fi; if n=0 then ERROR(`Empty argument`) elif n=1 then print(`The n=1 case is trivial.`); return() elif n=2 then print(`By the binomial theorem,`); if xpwrs=1 then print(`The constant term`) else print(`The coefficient of`, xpwrs) fi; print(`in the expansion of`); print(DysonProd(x,a,[seq(0,i=1..n)])); print(` is :`); print( (-1)^(b[1]) * binomial(a[1]+a[2],a[1]+b[1])); return() else GoodThm(b,THEOREM) fi; return() end: GoodThm:=proc(b,title) local B,i,j,R,n,x,a,xx,aa,bb,xsub,asub,bsub,LP, xxm,aam,bbm,one,bzero,IC,nextlist,xpwrs; x:=[x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20]; a:=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20]; B:=[b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,b16,b17,b18,b19,b20]; n:=nops(b); xx:=[seq(x[i],i=1..n)]; xxm:=[seq(x[i],i=1..n-1)]; aa:=[seq(a[i],i=1..n)]; aam:=[seq(a[i],i=1..n-1)]; bb:=[seq(B[i],i=1..n)]; bbm:=[seq(B[i],i=1..n-1)]; one:=[seq(1,i=1..n)]; if type(_LL[op(b)],indexed) then R:=GuessDysonCoeff(op(b)) else R:=_LL[op(b)] fi; LP:=Boundary(xx,aa,b); if title=THEOREM then print(`A THEOREM BASED ON A VARIATION OF THE`); print(N=n); print(`CASE OF THE 1962 DYSON CONJECTURE`); print(`by Shalosh B. Ekhad`); fi; print(`=======================================`); print(title); xpwrs:=product(x[i]^b[i],i=1..n); if xpwrs=1 then print(`The constant term`) else print(`The coefficient of`, xpwrs) fi; print(`in the expansion of`); print(DysonProd(x,a,[seq(0,i=1..n)])); print( `is: `); print(R*M(aa)); print(`=========================================`); print(`PROOF:`); print(`Define:`); print(CT(X),`= the constant term of X`); print(F[n](xx,aa,bb)= DysonProd(xx,aa,bb)); print(F[n-1](xxm,aam,bbm) = DysonProd(xxm,aam,bbm)); print(c[n](aa,bb)= CT(F[n](xx,aa,bb))); print(c[n-1](aam,bbm) = CT(F[n-1](xxm,aam,bbm))); print(d[n](aa,b)= R*M(aa)); print(`----------------------------------------`); print(`It is easily verified that we have the following recursion, for`); print(seq(a[i]>0,i=1..n),` `); print(F[n](xx,aa,bb) = add(subs(a[i]=a[i]-1,F[n](xx,aa,bb)),i=1..n)); print(`and thus, by applying the constant term operator to both sides,`); print(c[n](aa,bb) = add(subs(a[i]=a[i]-1,c[n](aa,bb)),i=1..n)); print(`---------------------------------------`); print(`By explicit series expansion, we obtain the boundary conditions:`); for i from 1 to n do xsub:=[seq(xx[j],j=1..i-1),seq(xx[j],j=i+1..n)]; asub:=[seq(aa[j],j=1..i-1),seq(aa[j],j=i+1..n)]; bzero:=[seq(0,j=1..n-1)]; print(CT(subs(a[i]=0,F[n](xx,aa,b))) = CT(F[n-1](xsub,asub,bzero)*LP[i])) od; print(`and thus`); nextlist:={}; for i from 1 to n do xsub:=[seq(xx[j],j=1..i-1),seq(xx[j],j=i+1..n)]; asub:=[seq(aa[j],j=1..i-1),seq(aa[j],j=i+1..n)]; if type(LP[i],`+`) then nextlist:= nextlist union {seq( [seq(-degree(op(j,LP[i]),x[l]),l=1..i-1), seq(-degree(op(j,LP[i]),x[l]),l=i+1..n)], j=1..nops(LP[i])) }; print(subs(a[i]=0,c[n](aa,b)) = add( c[n-1](asub, [seq(-degree(op(j,LP[i]),x[l]),l=1..i-1), seq(-degree(op(j,LP[i]),x[l]),l=i+1..n)]) *op(j,LP[i])/ convert([seq(x[l]^(degree(op(j,LP[i]),x[l])),l=1..n)],`*`) , j=1..nops(LP[i]) ) ) else nextlist:=nextlist union {[seq(-degree(LP[i],x[l]),l=1..i-1), seq(-degree(LP[i],x[l]),l=i+1..n)]}; print(subs(a[i]=0,c[n](aa,b)) = c[n-1](asub, [seq(-degree(LP[i],x[l]),l=1..i-1), seq(-degree(LP[i],x[l]),l=i+1..n)]) *LP[i]/ convert([seq(x[l]^(degree(LP[i],x[l])),l=1..n)],`*`) ) fi od; print(`----------------------------------------`); print(`Initial condition:`); if {op(b)}={0} then IC:=1 else IC:=0 fi; print(c[n]([seq(0,i=1..n)],b) = IC); print(`----------------------------------------`); print(`Since the`, d[n](aa,b), `satisfy the same recursions, boundary conditions,`); print(`and initial condition as the`, c[n](aa,b),`the result follows.`); print(`==========================================================`); nextlist:=nextlist minus {[seq(infinity,i=1..n-1)]}; #print(`NEXTLIST`=nextlist); if n>3 then print(`...once the following lemmas are established:`); for i from 1 to nops(nextlist) do GoodThm(nextlist[i],LEMMA) od fi; return() end: #GenPol(Resh,a,deg): Given a list of variables Resh, a symbol a, #and a non-negative integer deg, outputs a generic polynomial #of total degree deg,in the variables Resh, indexed by the letter a of #followed by the set of coeffs. For example, try: #GenPol([x,y],a,1); GenPol:=proc(Resh,a,deg) local var,POL1,i,x,Resh1,POL: if not type(Resh,list) or nops(Resh)=0 then ERROR(`Bad Input`): fi: x:=Resh[1]: if nops(Resh)=1 then RETURN(convert([seq(a[i]*x^i,i=0..deg)],`+`),{seq(a[i],i=0..deg)}): fi: Resh1:=[op(2..nops(Resh),Resh)]: var:={}: POL:=0: for i from 0 to deg do POL1:=GenPol(Resh1,a[i],deg-i): var:=var union POL1[2]: POL:=expand(POL+POL1[1]*x^i): od: POL,var: end: #GenRat(Resh,a,b,degT,degB): The generic form of a rational function #in the list of variables Resh indexed by a on top and b at bottom #of total degree degT on top and degB at bottom #for example, try: GenRat([x,y],a,b,1,2) GenRat:=proc(Resh,a,b,degT,degB) local TOP,BOT: TOP:=GenPol(Resh,a,degT): BOT:=GenPol(Resh,b,degB): TOP[1]/BOT[1],TOP[2] union BOT[2]: end: #Simp(k,L): all the weakly-decreasing k-vectors of non-negative integers with #entries<=L, for example Simp(2,4); Simp:=proc(k,L) local gu,mu,i,j,i1: if k=1 then RETURN({seq([i],i=0..L)}): fi: gu:={}: for i from 0 to L do mu:=Simp(k-1,i): gu:=gu union {seq(seq([op(mu[i1]),j],j=i..L),i1=1..nops(mu))}: od: gu: end: #Simp1(k,m,L): all the weakly-decreasing k-vectors of non-negative #integers with #m<=entries<=L, for example Simp(2,3,4); Simp1:=proc(k,m,L) local gu,mu,i,j,i1: if k=1 then RETURN({seq([i],i=m..L)}): fi: gu:={}: for i from m to L do mu:=Simp1(k-1,m,i): gu:=gu union {seq(seq([op(mu[i1]),j],j=i..L),i1=1..nops(mu))}: od: gu: end: #GuessRat1(F,Resh,degT,degB,m): Given a function or procedure F with #nops(Resh) arguments (that are integers), guesses a rational function #with total degree degT and degB for numerator and denominator respectively #that seems fits it, using data of integer k-tuples that start at m GuessRat1:=proc(F,Resh,degT,degB,m) local eq,var,Rat,a,b,i1,j1,mu,k,L: k:=nops(Resh): Rat:=GenRat(Resh,a,b,degT,degB): var:=Rat[2]: Rat:=Rat[1]: eq:={}: for L from m while nops(Simp1(k,m,L))a/(a+b+1),[a,b],2,0); GuessRat:=proc(F,Resh,Deg,m) local lu,L: for L from 0 to Deg do #lu:=GuessRat1(F,Resh,L,Deg-L,m): lu:=GuessRat1(F,Resh,Deg-L,L,m): if lu<>FAIL then RETURN(factor(lu)): fi: od: FAIL: end: #VerifyMN(Rat,Resh):Given a rational function Rat in the variables #in the list Resh, verifies that Rat*MultiNomial(Resh) satisfies the #same the famous Pascal-Like recurrence of the Mult-nomial coeff. VerifyMN:=proc(Rat,Resh) local i,su: su:=convert(Resh,`+`): normal(su*Rat- convert([seq(subs(Resh[i]=Resh[i]-1,Rat)*Resh[i],i=1..nops(Resh))],`+`)): end: GuessDysonCoeff1:=proc(b, TotDeg) local a,c,ff,i,n,r,ss,vars; n:=nops(b); if add(b[i],i=1..n)<>0 then #print(`Trivial case!`); return(0) fi; a:=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20]; vars:= [seq(a[i],i=1..n)]; c:=proc() local x,a,gu,i,n,lu,j,f,s; n:=nargs; gu:=1; for i to n do lu:=1; #for j to i-1 do lu:=lu*(1-x[i]/x[j])^args[i] od; for j to i-1 do lu:=lu*(1-x[i]/x[j])^args[j] od; #for j from i+1 to n do lu:=lu*(1-x[i]/x[j])^args[i] od; for j from i+1 to n do lu:=lu*(1-x[i]/x[j])^args[j] od; gu:= gu*lu od; gu:=expand(gu); for i to n do gu:=coeff(gu,x[i],b[i]) od; ### find extra factors for accelerated ansatz: f:=1; s:=add(args[i],i=1..n); for i to n do if b[i]<0 then f:=f*FallFact(args[i],ceil(-b[i]/2))/rf(1+s-args[i],-b[i]) fi od; # f:=f*args[i]/rf(1+s-args[i],-b[i]) fi od; # for i to n do if b[i]>0 then f:=f/rf(1+args[i],floor(b[i]/2)) fi od; gu:=gu/M([args])/f; gu end proc: r:=GuessRat(c,vars,TotDeg,max(op(b))); ff:=1; ss:=add(a[i],i=1..n); for i to n do if b[i]<0 then ff:=ff*FallFact(a[i],ceil(-b[i]/2))/rf(1+ss-a[i],-b[i]) fi od; # ff:=ff*a[i]/rf(1+ss-a[i],-b[i]) fi od; # for i to n do if b[i]>0 then ff:=ff/rf(1+a[i],floor(b[i]/2)) fi od; return(r*ff) end: GuessDysonCoeff:=proc() local b,deg,r; b:=[args]; if add(b[i],i=1..nops(b))<>0 then print(`Trivial Case!`); return(0) fi; #if type(_LL[args],indexed) then for deg from max(args) do r:=GuessDysonCoeff1(b,deg); if r<>FAIL then #FindAndRecordSymmetric(args); return(r) fi od #else return(_LL[args]) #fi end: Multinomial:=proc(a) local i,n; n:=nops(a); add(a[i],i=1..n)!/product( a[i]!,i=1..n) end: M2:=proc(a,n) local i; add(a[i],i=1..n)! /product(a[i]!,i=1..n) end: M:=proc(a) Multinomial(a) end: Boundary:=proc(x,a,B) local b, i,j,k; #b:=[seq(-B[i],i=1..nops(B))]; b:=B; [seq( combine( expand(simplify( coeff( series( convert([seq( (x[i]-x[k])^a[i]/x[i]^(b[i]+a[i]),i=1..k-1)],`*`) *convert([seq( (x[j]-x[k])^a[j]/x[j]^(b[j]+a[j]),j=k+1..nops(b))],`*`) , x[k]=0, abs(b[k])+1) ,x[k],b[k]) ) ) ,power,symbolic) ,k=1..nops(b))] end: CanonicalVectors:=proc(n,T) local CP,i,m,negs,pos,s,tmp,t; if n-1>T then m:=T else m:=n-1 fi; pos:=combinat[partition](T,m); pos:=[seq(combinat[conjpart](pos[i]),i=1..nops(pos))]; negs:=[seq(-pos[i],i=1..nops(pos))]; CP:=combinat[cartprod]([negs,pos]); s:=[]; while not CP[finished] do tmp:=CP[nextvalue](); s:=[op(s),([op(tmp[1]),op(tmp[2])])] od; t:=[]; for i from 1 to nops(s) do if nops(s[i])=n then t:=[op(t), sort(s[i],`<`)] elif nops(s[i])`,_LL[subscript]) fi od; return() end: ElemVector:=proc(k,n) local i; [seq(0,i=1..k-1),1,seq(0,i=k+1..n)] end: OffsetVector:=proc(L,n) local i,v; v:=[seq(0,i=1..n-1),nops(L)]; for i from 1 to nops(L) do v:=v-ElemVector(L[i],n) od end: TurboDyson:=proc(n,C) local i,a; global _LL; # Finds and stores all n-Dyson coeffs of complexity <= C a:=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20]; _LL[seq(0,i=1..n)]:=1; print([seq(0,i=1..n)],`-->`,_LL[seq(0,i=1..n)]); _LL[-1,1,seq(0,i=1..n-2)]:=-a[1]/(1+add(a[i],i=2..n)); print([-1,1,seq(0,i=1..n-2)],`-->`,_LL[-1,1,seq(0,i=1..n-2)]); FindAndRecordSymmetric(-1,1,seq(0,i=1..n-2)); for i from 2 to C do TurboDyson1(n,i) od end: TurboDyson1:=proc(n,C) local a,b,CV,i,j,k,tmp,MV,NeededVectors; global _LL; a:=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20]; CV:=CanonicalVectors(n,C); print(`Canonical vectors of complexity`,C,CV); for i from 1 to nops(CV) do if type(_LL[op(CV[i])],indexed) then MV:=[seq(1,i=1..n-1),-n+1]; b:=CV[i]+MV; # print('B'=b); # print('MainVec'=MV); NeededVectors:=[b]; for j from 1 to n-2 do tmp:=combinat[choose](n-1,j); for k from 1 to nops(tmp) do # print('temp'=tmp); NeededVectors:=[op(NeededVectors),b+OffsetVector(tmp[k],n)] od od; print(`In order to get`,CV[i],`we will need`,NeededVectors); for j from 1 to nops(NeededVectors) do if type(_LL[op(NeededVectors[j])],indexed) then _LL[op(NeededVectors[j])]:=GuessDysonCoeff(op(NeededVectors[j])); print([op(NeededVectors[j])],`-->`,_LL[op(NeededVectors[j])]); FindAndRecordSymmetric(op(NeededVectors[j])) fi od; ## Calculate _LL[op(CV[i])] _LL[op(CV[i])]:=factor(normal( (-1)^(n+1)*(1+add(a[i],i=1..n))/(1+a[n]) *subs(a[n]=a[n]+1,_LL[op(NeededVectors[1])]) +(-1)^n*_LL[op(NeededVectors[1])] +add((-1)^(n+ op(n,NeededVectors[j])-b[n] ) * _LL[op(NeededVectors[j])],j=2..nops(NeededVectors)))); print([op(CV[i])],`-->`,_LL[op(CV[i])]); FindAndRecordSymmetric(op(CV[i])) fi od; return() end: GetDysonCoeff:=proc() local R; global _LL; if type(_LL[args],indexed) then R:=GuessDysonCoeff(args); _LL[args]:=R; return(R) else return(_LL[args]) fi end: DysonExpansion:=proc(n,C) local A,AA,c,i,L; AA:=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20]; A:=[ seq(AA[i],i=1..n)]; L:=[seq(op(AllVectors(n,c)),c=0..C)]; print( DysonProd(x,a,[seq(0,i=1..n)]) = M(A)* add( GetDysonCoeff(op(L[i]))*product(x[j]^L[i][j],j=1..n) ,i=1..nops(L) ) + `. . .` ); end: