#######################################################################
## GoodDyson: Save this file as GoodDyson. To use it, stay in the     #
## same directory, get into Maple (by typing: maple <Enter> )         #
## and then type:  read GoodDyson: <Enter>                            #
## 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))<nops(var)+10 do
od:
mu:=Simp1(k,m,L):

eq:=
{seq(
numer(normal(F(op(mu[i1]))-subs({seq(Resh[j1]=mu[i1][j1],j1=1..k)},Rat))),
i1=1..nops(mu))}:

eq:=eq minus {0}:

if nops(eq)<nops(var)+5 then
#print(`Too many zeroes, try to raise L `):
RETURN(FAIL)
fi:

#var:=SolveTools[Linear](eq,var):
var:=solve(eq,var):
if subs(var,denom(Rat))=0 or var=NULL then
  RETURN(FAIL)
else
 RETURN(normal(subs(var,Rat)))
fi
end:


#GuessRat(F,Resh,deg): Given a function or procedure F with
#nops(Resh)  arguments (that are integers), guesses a rational function
#with total degrees of sums of degrees of top and bottom up to
#deg that seems fits it, or  RETURNS FAIL
#For example, try: GuessRat1((a,b)->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])<n then
    tmp:= sort([op(s[i]),seq(0,i=1..n-nops(s[i]))],`<`);
    t:=[op(t), tmp]
  fi
od;
t
end:

AllVectors:=proc(n,c) local L,i;
L:=CanonicalVectors(n,c);
[seq(op( combinat[permute](op(i,L))), i=1..nops(L) )];
end:


FindAndRecordSymmetric:=proc()
local a,b,i,j,n,L,R,subscript,substitutions,varlist;
global _LL;
n:=nargs;
a:=[a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20];
if type(_LL[args],indexed) then
   R:=GuessDysonCoeff(args);
   _LL[args]:=R
else R:=_LL[args]
fi;
L:=combinat[permute](n);
varlist:=[seq(a[i],i=1..n)];
for i from 2 to nops(L) do
#print(`permutation in use:`,L[i]);
subscript:= seq( args[L[i][j] ],j=1..n); 
if type(_LL[subscript],indexed) then 
  #print(subscript);
  substitutions:={seq(a[L[i][j]]=a[j],j=1..n) };
  #print(substitutions);
  _LL[subscript]:=subs(substitutions,R);
  #print(`In FindAndRecordSymmetric`,[subscript],`-->`,_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: