BeginPackage["GoodDyson`"]

<<Combinatorica`;

Print["GoodDyson Mathematica Package, Version 1.102, 1 March 2009"];
Print["written by Andrew Sills"];
Print["inspired by Sills & Zeilberger's Maple package `GoodDyson.`"];
Print["and Doron Zeilberger's Maple package `GuessRat`."];
Print["The package accompanies the article"];
Print["`Disturbing the Dyson Conjecture (in a GOOD way)` by Sills & Zeilberger"];
Print["Experimental Math. 15(2006), 187-191."]
Print["Latest version available at http://math.georgiasouthern.edu/~asills."];
Print["Please report bugs and/or send comments to asills at georgiasouthern dot edu"];
Print["For a listing of available functions, enter 'GoodDysonFunctions[]'."];

GoodDysonFunctions::usage "GoodDysonFunctions[] lists available functions in
the GoodDyson Mathematica package."

GuessRat::usage = "GuessRat[F,Resh,deg,M]:  Given a function F with argument 
list Resh, guesses a rational function representation P/Q 
(with degree(P)+degree(Q) at most deg) of F using internally
generated data that uses the integer M as a seed.  Returns 'FAIL' if
unsuccessful."

Freeman3::usage = "Freeman3[b1,b2,T]: guesses a closed form formula for
the coefficient of x^b1 y^b2 z^(-b1-b2) in the expansion of
((1-x/y)(1-x/z))^a ((1-y/x)(1-y/z))^b ((1-z/x)(1-z/y))^c in the form
P(a,b,c)/Q(a,b,c) Multinomial[a,b,c], where P and Q are polynomials
such that deg(P)+deg(Q)<= T.  If no such closed form expression is
found, the function returns 'Fail'."

GuessRatFunctions::usage = "GuessRatFunctions[] prints the list of
functions in this package which are available to the user."

GuessDysonCoeff::usage="GuessDysonCoeff[b]:  Given that b is a list of n 
integers {b1,b2,...,bn}, this procedure attempts to find a closed form 
expression for the coefficient of x1^b1 x2^b2 ... xn^bn in the Dyson product.
The result returned is not the closed form expression itself, but rather 
the rational function multiple of the multinomial
coefficient. "

GetDysonCoeff::usage="GetDysonCoeff[b]: Similar to `GetDysonCoeff`, except
that this procedure first checks to see if the result has already been
calculated, and if so, returns it immediately.  If not, it automatically
calls `GuessDysonCoeff`. "

WritePaper::usage="WritePaper[b], where b=={b1,b2,...bn} is a list of
integers states and proves a theorem which gives an explicit closed form
expression for the coefficient of x1^b1 x2^b2 ... xn^bn in the expansion
of the Dyson product. "

IntegerPartition::usage="IntegerPartition[n,m] returns a list of partitions
of the integer n into parts less than or equal to m. "

IntegerPartRNP::usage="IntegerPartitionRNP[n,m,k] returns a list
of partitions of the integer in into at most k parts, each less than or
equal to m."

TurboDyson::usage="TurboDyson(n,C) calculates and stores for future use all
coefficients of x1^b1 * x2^b2 * ... * xn^bn
in the expansion of the Dyson product
which have complexity less than or equal to C.
The complexity C of (b1,b2,b3,...,bn) is the sum of the
positive components of the vector (b1,b2,...,bn).";

Begin["`Private`"]

GoodDysonFunctions[]:=
  Module[{},
    Print["The GoodDyson package includes the following functions:
GetDysonCoeff, GuessDysonCoeff, TurboDyson, WritePaper.  
Enter ?<function name> for help on a specific function."]
  ];

GuessRatFunctions[]:=
  Module[{},
   Print["Available functions: GuessRat, Freeman3"];
   Print["Enter '?<function name>' for a description of the function <function name>"] 
  ];

MC[a_list]:=
  Module[{i,n},
    n=Length[a];
    Return[Factorial[Sum[a[[i]],{i,n}] ]/ Product[ Factorial[ a[[i]] ],{i,n}]]
  ];

IntegerPartition[n_, m_]:=
  Module[{i,j,p,t},
   If[n==0, Return[ {{}} ]];
   If[m==1, Return[ {Table[ 1,{i,n}]}]];
   p:={};
   For[ i=1, i<=m, i++,
     t = IntegerPartition[ n-i, Min[n-i,i]];
     p = Union[ p, Table[ Append[ t[[j]],i], {j,Length[t]}]]
   ];
   Return[p]
  ];

IntegerPartRNP[n_, m_, k_]:=
  Module[{i,p,o},
    p = IntegerPartition[n,m];
    o = {};
    For[ i=1, i<=Length[p], i++,
      If[ Length[ p[[i]] ]<=k, o= Append[o, p[[i]] ]
      ]
    ];
    o
  ]; 

GenPol[Resh_, a_, deg_]:=
  Module[{var,POL1,i,x,Resh1,POL},
   If[Length[Resh]==0, 
    Return[ {Sum  [ Subscript[a,i] Resh^i, {i,0,deg}],
             Table[ Subscript[a,i]       , {i,0,deg}]}
          ]
   ];
   x = Resh[[1]];
   If[Length[Resh]==1, 
    Return[ {Sum  [ Subscript[a,i] x^i, {i,0,deg}],
             Table[ Subscript[a,i],     {i,0,deg}]}
          ]
   ];
   Resh1 = Delete[Resh,1];
   var = {};
   POL = 0;
   For[ i=0, i<=deg, i++,
        POL1 = GenPol[ Resh1, Subscript[a,i], deg-i];
        var = Union[ var, POL1[[2]] ];
        POL = Expand[ POL + POL1[[1]] x^i ];
   ];
   {POL,var}
  ];


GenRat[ Resh_, a_, b_, degT_, degB_]:=
 Module[{TOP,BOT},
   TOP = GenPol[Resh,a,degT];
   BOT = GenPol[Resh,b,degB];
   {TOP[[1]]/BOT[[1]], Union[ TOP[[2]], BOT[[2]] ] }
 ];


OffsetList[ L_, a_ ]:=
  Module[{i},
   Table[ L[[i]]+a, {i,Length[L]}]
  ];

Simp[k_, M_, L_]:=
  Module[{n,i},
    p = {};
    For[n=0, n<= k(L-M), n++, 
         p = Union[p, IntegerPartRNP[n,L-M,k]]
    ];
    Table[ OffsetList[PadLeft[p[[i]],k],M],{i,Length[p]}]
  ];
  
GuessRat1[F_, Resh_, degT_, degB_, M_ ]:=
  Module[{eq,var,sol,Rat,a,b,i,j,mu,k,l,L,r},
    k = Length[Resh];
    Rat = GenRat[ Resh, a, b, degT, degB ];
    var = Rat[[2]];
    Rat = Rat[[1]];
    For[ L=M, Length[mu]< Length[var]+10, L++, mu = Simp[k,M,L]];
    eq = Table[
           Numerator[
             Together[
               (* Apply[ F,  mu[[i]] ] *)
               F[ mu[[i]] ]
               - (Rat/.Table[ Resh[[j]] -> mu[[i,j]], {j,k} ])
             ]
           ]==0,
         {i,Length[mu]}
         ];
  
    eq = Complement[eq,{0==0,False,True}];  
    If[ Length[eq] < Length[var] + 5, 
         Print["Too Many Zeros.  Try to raise L"];
         Return[FAIL]
      ];
    sol = Solve[ eq, var ];
    r = Denominator[Rat]/.sol;
    If[ r == {0} || sol =={} , Return[FAIL]]; 
    r = Simplify[Rat/.sol];
    If[ Length[r] == 1, Return[ r[[1]] ], Return[r] ]
  ]; 


GuessRat[F_, Resh_, deg_, M_]:=
  Module[{lu,m},
    For[m=0, m<=deg, m++,
        lu = GuessRat1[F,Resh,deg-m,m,M]; 
        If [ lu =!= FAIL , Return[lu] ] 
      ];
    Return[FAIL]
   ];

Dyson3[b1_, b2_]:=
  Module[{k,L,F},
    L= Max[ Abs[b1], Abs[b2], Abs[-b1-b2] ];
    F[a1_, a2_, a3_]:=
      Block[{x,y,z,xpwrlist},
        xpwrlist=Complement[ {x^b1, y^b2, z^(-b1-b2)},{1}];
        Fold[ Coefficient, 
              ((1-x/y)(1-x/z))^a1 ((1-y/x)(1-y/z))^a2 ((1-z/x)(1-z/y))^a3/
              Multinomial[a1,a2,a3],
             xpwrlist 
            ]
      ];
    Return[F]
];

Freeman3[b1_, b2_, T_, x_:"x", y_:"y", z_:"z", a_:"a", b_:"b", c_:"c"]:=
   Block[{lu, deg},
     lu = FAIL;
     For[deg=1, deg<=T && lu==FAIL, deg++,
       lu = GuessRat[Dyson3[b1,b2],{a,b,c},deg,Max[b1,b2]];
       Print["Dyson3[b1,b2]=",Dyson3[b1,b2] ];
       Print["deg=",deg,"  lu=",lu]
     ];
     If[ lu =!= FAIL,
      Print["The coefficient of ",x^b1 y^b2 z^(-b1-b2)," in the expansion of"];
      Print[((1-x/y)(1-x/z))^a ((1-y/x)(1-y/z))^b ((1-z/x)(1-z/y))^c];
      Print["is"];
      Return[Multinomial[a,b,c] Factor[lu]],
      Return[Fail]
     ]
   ];


FallingFact[a_, n_]:= Module[{i},Product[a-i,{i,0,n-1}]];
       
DysonProduct[x_, a_, b_List]:=
  Module[{i,j,n},
    n = Length[b];
    Product[ Subscript[x,i]^(-b[[i]]), {i,1,n}] *
    Product[ (1-Subscript[x,i]/Subscript[x,j])^Subscript[a,j] *
             (1-Subscript[x,j]/Subscript[x,i])^Subscript[a,i],
        {i,n-1},{j,i+1,n}
           ]   
  ];

WritePaper[b_List]:=
   Module[{n,i,j,xpwrs},
     n = Length[b];
     xpwrs = Product[ Subscript[x,i]^b[[i]], {i,n}];
     If[ Sum[ b[[i]],{i,n}] != 0,
        If[xpwrs === 1, 
           Print["The constant term"], 
           Print["The coefficient of ",xpwrs] 
        ];
        Print["in the expansion of"];
        Print[DysonProduct[x,a, Table[0,{n}] ] ];
        Print["is trivially zero."];
        Return[]
     ];
     Which[ n==0, Print["Empty argument"]; Return[],
            n==1, Print["The case n=1 is trivial."]; Return[],
            n==2, Print["By the binomial theorem,"];
                  Print["The coefficient of ", xpwrs," in the expansion of"];
                  Print[DysonProduct[x,a,Table[0,{n}] ] ];
                  Print["is"];
                  Print[ (-1)^b[[1]] *
                    Binomial[Subscript[a,1]+Subscript[a,2],
                      Subscript[a,1]+b[[1]]
                    ]
                  ];
                  Return[],
            n>2, GoodThm[b,"THEOREM"]; Return[]
      ]
   ];

DysonN[b_List]:=
  Module[{k,L,F,n},
    L= Apply[ Max, b ];
    n = Length[b];
    F[a_List]:=
      Block[{f,s,i,j,x,xpwrlist},
       xpwrlist=Complement[Table[Subscript[x,i]^b[[i]],{i,n}],{1}];
       s = Sum[ a[[i]], {i,n}];
       f = 1;
       Do[ 
           If[ b[[i]]<0, 
              f = f*FallingFact[a[[i]], Ceiling[ -b[[i]]/2] ] /
                    Pochhammer[1 + s - a[[i]], -b[[i]] ]
           ]
       ,{i,n}
       ];
       Fold[ Coefficient,
         Product[(1-Subscript[x,i]/Subscript[x,j])^a[[j]] *
             (1-Subscript[x,j]/Subscript[x,i])^a[[i]],
             {i,n-1},{j,i+1,n}
         ]/ Apply[Multinomial,a]/ f,
            xpwrlist
        ]
      ];
    Return[F]
];

GuessDysonCoeff1[b_List, TotDeg_Integer, a_List]:=
   Block[{lu, deg, n, i,ff,ss },
      n = Length[b];
      If[ Sum[b[[i]], {i,n}] != 0, Return[0] ];
      If[ Union[b] == {0}, Return[ 1 ] ];
      lu = FAIL;
      For[deg=1, deg<=TotDeg && lu==FAIL , deg++,
          lu = GuessRat[DysonN[b], a, deg, Apply[Max,b] ];
      ];
      If[lu =!= FAIL,
         ss = Sum[ a[[i]], {i,n}];
         ff = 1;
         Do[
            If[ b[[i]]<0 ,
                ff = ff*FallingFact[ a[[i]], Ceiling[-b[[i]]/2] ] /
                        Pochhammer[1 + ss - a[[i]], -b[[i]] ]
            ]
         ,{i,n}
         ]; 
         Return[ff Factor[lu]], 
         Return[Fail]
      ]
   ];

GuessDysonCoeff[b_List,a_List]:=
  Module[{deg,i,r},
    If[ Sum[b[[i]],{i,Length[b]}] != 0, Return[0] ];
    r = Fail;
    For[ deg=Apply[Max,b],r==Fail, deg++,
      r= GuessDysonCoeff1[b,deg,a];
    ];
    Return[r]
  ];

x="x";
a="a";

GoodThm[b_List, title_]:=
  Module[{an,anm,asub,avars,bn,bnm,bzero,i,IC,j,l,LP,n,nextlist,R,varsubs,xn,xnm,xpwrs,xsub,xvars},
    n = Length[b];
    an = Table[Subscript["a",i],{i,n}];
    anm= Table[Subscript["a",i],{i,n-1}];
    bn = Table[Subscript["b",i],{i,n}];
    bnm= Table[Subscript["b",i],{i,n-1}];
    xn = Table[Subscript["x",i],{i,n}];
    xnm= Table[Subscript["x",i],{i,n-1}];
    xvars = {x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15};
    xvars = Table[xvars[[i]],{i,n}];
    avars = {a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15};
    avars = Table[avars[[i]],{i,n}];
    varsubs= Join[Table[xvars[[i]]->xn[[i]],{i,n}],
             Table[avars[[i]]->an[[i]],{i,n}] ];
    If[ Head[ LL[b] ]===LL,
        R = GuessDysonCoeff[b, an],
        R = LL[b]
    ];
    (*  R = GuessDysonCoeff[b, Table[ Subscript[a,i], {i,n}] ]; *)
    LP = Boundary[xvars, avars, b]/.varsubs;
    If[ title=="THEOREM",
        CellPrint[
          Cell[
             TextData[{ "A theorem based on a variation of the n = ",
                Cell[ BoxData[
                        FormBox[ n, TraditionalForm] 
                      ]
                ],
             " case of the 1962 Dyson conjecture"            
             }],
             "Title"
          ]
        ];
    ];
    CellPrint[
          Cell[
             TextData[{ " ",
                Cell[ BoxData[
                        FormBox[ title , TraditionalForm]
                      ]
                ],
             " "
             }],
             "Section"
          ]
        ];

    xpwrs = Product[ Subscript[x,i]^b[[i]], {i,n}];
    If[ xpwrs===1, 
        Print["The constant term"], 
        Print["The coefficient of ",xpwrs]
    ];
    Print["in the expansion of"];
    Print[ DysonProduct[x,a, Table[0,{n}] ] ]; 
    Print["is"];
    Print[R Apply[Multinomial,an]//TraditionalForm ];
    Print["where ",Apply[Multinomial,an]//TraditionalForm," = ", Apply[Multinomial,an] ];
    Print["===================================================="];
    Print[StyleForm["Proof:", FontSlant->"Italic", FontWeight->"Bold", FontSize->18]];
    Print["----------------------------------------"];
    Print["Define the following notation:"];
    Print["CT[X] = the constant term of X"];
    Print[Subscript["F",n][xn,an,bn],"=", DysonProduct[x,a,bn] ];
    Print[Subscript["c",n][an,bn],"=", "CT"[Subscript["F",n][xn,an,bn]] ];
    Print[Subscript["d",n][an,b],"=",R Apply[Multinomial,an]//TraditionalForm  ];
    Print["---------------------------------------"];
    Print["It is easily verified that we have the following recursion:"];
    Print["For ",an,">", Table[0,{n}]];
    Print[Subscript["F",n][xn,an,bn],"==", 
          Sum[(Subscript["F",n][xn,an,bn])/.Subscript[a,i]->Subscript[a,i]-1
           ,{i,n}
         ]
    ];
    Print["and thus, by applying the constant term operator to both sides,"];
    Print[Subscript["c",n][an,bn],"==",
           Sum[(Subscript["c",n][an,bn])/.Subscript[a,i]->Subscript[a,i]-1
           ,{i,n}
           ]
    ];
    Print["----------------------------------"];
    Print["By explicit series expansion, we obtain the boundary conditions"];
    Do[
      xsub=Join[Table[xn[[j]],{j,i-1}],Table[xn[[j]],{j,i+1,n}] ];
      asub=Join[Table[an[[j]],{j,i-1}],Table[an[[j]],{j,i+1,n}] ];
      bzero=Table[0,{n-1}];
      Print["CT"[ (Subscript["F",n][xn,an,b]/. an[[i]]->0) ],"==",
            "CT"[ Subscript["F",n-1][xsub,asub,bzero] * LP[[i]] ]],
      {i,n}
    ];
    Print["and thus,"];
    nextlist={};
    Do[  xsub =Join[Table[xn[[j]],{j,i-1}],Table[xn[[j]],{j,i+1,n}] ];
         asub =Join[Table[an[[j]],{j,i-1}],Table[an[[j]],{j,i+1,n}] ];
         If[ Head[ LP[[i]] ]=== Plus,
            nextlist = Union[nextlist,
                             Table[ 
                                Join[
                                  Table[-Exponent[ LP[[i,j]], 
                                         Subscript["x",l]
                                         ],
                                  {l,i-1}
                                  ], 
                                  Table[-Exponent[ LP[[i,j]], 
                                         Subscript["x",l]
                                         ],
                                  {l,i+1,n}
                                  ]
                                ],
                             {j,Length[LP[[i]] ]}
                             ]
                       ];
            Print[ Subscript["c",n][an,b]/.Subscript["a",i]->0 , "==",
                   Sum[ Subscript["c",n-1][asub,
                                           Join[ 
                                              Table[-Exponent[ LP[[i,j]],
                                                               Subscript["x",l]
                                                     ], 
                                              {l,i-1}
                                              ],
                                              Table[-Exponent[ LP[[i,j]],
                                                               Subscript["x",l]
                                                     ],
                                              {l,i+1,n} 
                                              ]
                                           ] 
                                          ] * LP[[i,j]] /
                                              Product[ Subscript["x",l]^
                                                       (Exponent[ LP[[i,j]],
                                                        Subscript["x",l] 
                                                        ]),
                                              {l,n}
                                              ], 
                   {j,Length[LP[[i]] ]}
                   ]  
            ]  
        , 
          nextlist = Union[ nextlist,
                       {Join[
                         Table[Exponent[LP[[i]],Subscript["x",l] ],
                         {l,i-1}
                         ],
                         Table[Exponent[LP[[i]],Subscript["x",l] ],
                         {l,i+1,n}
                         ]
                       ]} 
                   ]; 
          nextlist = Complement[nextlist, {Table[Infinity,{n-1}]} ];
          If[ MemberQ[ 
               Join[
                 Table[-Exponent[LP[[i]],Subscript["x",l]],
                 {l,i-1}
                 ],
                 Table[-Exponent[LP[[i]],Subscript["x",l]],
                 {l,i+1,n}
                 ]
               ], 
               Infinity 
              ],
            Print[ Subscript["c",n][an,b]/.Subscript["a",i]->0, " == ", 0]
            , 
            Print[ Subscript["c",n][an,b]/.Subscript["a",i]->0, " == ",
                 Subscript["c",n-1][asub,
                                   {Join[  
                                     Table[ 
                                       -Exponent[LP[[i]], Subscript["x",l]],
                                     {l,i-1}
                                     ],
                                     Table[  
                                       -Exponent[LP[[i]], Subscript["x",l]],
                                     {l,i+1,n}
                                     ]
                                  ]} 
                                 ] * LP[[i]] /
                                     Product[
                                           Subscript["x",l]^
                                          (Exponent[LP[[i]],Subscript["x",l]]),
                                     {l,n}
                                     ] 
            ]
          ]   
        ] 
    ,{i,n}
    ]; 
    Print["----------------------------------"];
    Print["Initial Condition:"];
    If[ Union[b]=={0}, IC= 1, IC=0];
    Print[Subscript["c",n][ Table[0,{n}], b]," == ",IC];
    Print["----------------------------------"];
    Print["Since the ", Subscript["d",n][an,b]," satisfy the same recursions, boundary conditions, and initial condition as the ", Subscript["c",n][an,b], ", the result follows."];
    Print["==========================================="];
    If[ n>3,
        Print["...once the following lemmas are established."];
        Do[ GoodThm[ nextlist[[i]], "LEMMA"],
        {i, Length[nextlist]}
        ]
    ];
Return[] 
];

Boundary[X_List,A_List,b_List]:=
  Module[{i,j,k,n},
     n = Length[b];
     Table[
       Expand[
       Simplify[
        Coefficient[
         Series[ 
           Product[ (X[[i]]-X[[k]])^A[[i]]/
                     X[[i]]^(b[[i]]+A[[i]]),{i,k-1}]*
           Product[ (X[[j]]-X[[k]])^A[[j]]/
                     X[[j]]^(b[[j]]+A[[j]]),{j,k+1,n}],
         {X[[k]],0, Abs[b[[k]]+1]} 
         ],
        X[[k]], b[[k]] 
        ]
       ]
       ]
     ,
     {k,n}
     ]
  ];

CanonicalVectors[n_,T_]:=
Module[{CP,i,m,negs,pos,s,t},
  If[ n-1 > T, m = T, m = n-1 ];
  pos = IntegerPartition[T,m];
  pos = Table[ Reverse[TransposePartition[Reverse[pos[[i]]]]],{i,Length[pos]}];
  (*negs = Table[ -pos[[i]], {i,Length[pos]} ];*)
  negs=-pos;
  CP = CartesianProduct[ negs, pos ];
  t=   Table[ Join[CP[[i,1]], CP[[i,2]] ],{i,Length[CP]} ];
  s = {};
  Do[
    If[ Length[ t[[i]] ] <= n,
        s = Union[s, {Sort[ PadLeft[ t[[i]], n] ]} ]
    ],
  {i,Length[t]} 
  ];
  s 
];

FindAndRecordSymmetric[b_List]:=
Module[{a,i,j,n,L,R,args,subst,varlist},
n = Length[b];
a = Table[Subscript["a",i],{i,n}];
If[ Head[ LL[b ]] === LL, 
    R=GuessDysonCoeff[b,a]; LL[b]=R, 
    R=LL[b] ];
L = Permutations[n];
Do[
   args = Table[ b[[ L[[i,j]] ]], {j,n}];
   If[ Head[ LL[args] ] === LL, 
       subst = Table[ a[[L[[i,j]] ]] -> a[[j]], {j,n} ];
       (*Print["Recording ", args];*)
       LL[args] = R/. subst;
   ], 
{i,Length[L]}
];
Return[]
];

GetDysonCoeff[b_List]:=
Module[{i,R},
  If[ Head[ LL[b] ] === LL, 
      R=GuessDysonCoeff[b,Table[ Subscript["a",i],{i,Length[b]}] ]; 
      LL[b]=R; Return[R],
      Return[ LL[b] ]
  ]
];

ElemVector[k_,n_]:=Module[{i},  PadRight[PadLeft[{1},k],n] ];

OffsetVector[L_,n_]:=
Module[{i,v},
 v = PadLeft[{Length[L]},n];
 Do[ v=v-ElemVector[ L[[i]], n], {i,Length[L]}];
 v
];

TurboDyson[n_,C_]:=
Module[{i,a},
 a = Table[ Subscript["a",i], {i,n}];
 LL[ Table[0,{n}] ]=1;
 Print[ Table[0,{n}],"-->", LL[ Table[0,{n}] ] ];
 LL[ PadRight[{-1,1},n] ] = -a[[1]] /(1+ Sum[ a[[i]],{i,2,n}]); 
 FindAndRecordSymmetric[ PadRight[ {-1,1}, n] ];
 Do[ TurboDyson1[n,i], {i,2,C} ];
 Return[]
];



TurboDyson1[n_, C_]:=
Module[{a,b,CV,i,j,k,l,tmp,MV,NeededVectors},
  a = Table[Subscript["a",i],{i,n}];
  CV = CanonicalVectors[n,C];
  Print["Canonical vectors of complexity ",C," : ", CV];
  Do[
     If[ Head[ LL[ CV[[i]] ] ]===LL,
         MV = Join[Table[ 1,{n-1}], {1-n}];
         b = CV[[i]] + MV;
      (*   Print["B=",b];  *)
      (*   Print["MainVec=",MV]; *)
         NeededVectors = {b};
         Do[ tmp = KSubsets[ Table[l,{l,n-1}], j];
             (*Print["tmp=",tmp];*)
             Do[ NeededVectors = Append[ NeededVectors, b + OffsetVector[ tmp[[k]],n]];
              ,
             {k,Length[tmp]}
             ],
         {j,n-2}
         ];
         Print["In order to get ", CV[[i]], " we will need ",NeededVectors]; 
         Do[ 
             If[ Head[ LL[ NeededVectors[[j]] ] ]=== LL,
                 LL[ NeededVectors[[j]] ] = 
                    GuessDysonCoeff[NeededVectors[[j]], a ];
                 Print[ NeededVectors[[j]],"--->", LL[ NeededVectors[[j]] ] ];
                 FindAndRecordSymmetric[ NeededVectors[[j]] ]
             ],
         {j,Length[NeededVectors]}  
         ];
         (*Print["Last part"]; *)
         LL[ CV[[i]] ] = Factor[Together[  
                           (-1)^(n+1) * (1+ Sum[ a[[i]], {i,n} ])/ (1+a[[n]] )*
                           ( LL[ NeededVectors[[1]] ]/. a[[n]] -> 1+a[[n]] )
                        +(-1)^n * LL[ NeededVectors[[1]] ]
                        + Sum[ (-1)^(n+ NeededVectors[[j,n]] - b[[n]] ) 
                         *LL[ NeededVectors[[j]]  ], {j,2,Length[NeededVectors]}]]];
        Print[ CV[[i]],"-->",LL[ CV[[i]] ] ];
        FindAndRecordSymmetric[ CV[[i]] ]
     ],
  {i,Length[CV]}
  ]
];


End[]

EndPackage[]