####################################################################
##GEPENR.txt: Save this file as  GEPNER.txt                          #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read `GEPNER.txt`<Enter>                           #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: 
 
print(`Created:  March 31, 2016`):
print(` This is GEPNER.txt `):
print(`It is the Maple package that accompanies the article `):
print(`   Doron Gepner's Statistics on Words in {1,2,3}^* is (Most Probably) Asymptotically Logistic`):
print(`by Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AveAndMoms, Chop, DGn, DoronSeq12, F1abc, F1abcD, Hafokh, G1abc, G1abcT, gep, GepMomsD, GP, Milim , Mishkal `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: DGseq, Fabc, Gabc,  gepPmoms, GepMoms, GepSeq `):
 print(` `):



elif nops([args])=1 and op(1,[args])=AveAndMoms then
print(`AveAndMoms(f,x,N): Given a probability generating function`):
print(`f(x) (a polynomial, rational function and possibly beyond) `):
print(` returns a list whose first entry is the average  `):
print(` (alias expectation, alias mean) `):
print(` followed by the variance, the third moment (about the mean) and `):
print(` so on, until the N-th moment (about the mean). `):
print(` If f(1) is not 1, than it first normalizes it by dividing `):
print(` by f(1) (if it is not 0) . `):
print(` For example, try: `):
print(` AveAndMoms(((1+x)/2)^100,x,4); `):


elif nops([args])=1 and op(1,[args])=Chop then
print(`Chop(P,var,K): chops the polynomial P of the list of variables var to its terms of (total) degree<=K`):
print(`Try: Chop((x+y)^5,[x,y],2);`):

elif nops([args])=1 and op(1,[args])=DGseq then
print(`DGseq(q,N): The list, of length N, of weight-enumerator of the set of n-permutations according to the variable q and the gep-statistics. `):
print(`from n=1 to n=N. Try:`):
print(`DGseq(q,6);`):

elif nops([args])=1 and op(1,[args])=DGn then
print(`DGn(q,n): The weight-enumerator of the set of n-permutations according to the variable q and the gep-statistics. Try:`):
print(`DGn(q,5);`):

elif nops([args])=1 and op(1,[args])=DoronSeq12 then
print(`DoronSeq12(q): The pre-computed sequence of the first 12 Gepner polynomials. In other words the output of GepSeq(12,q) `):

elif nops([args])=1 and op(1,[args])=Fabc then
print(`Fabc(a,b,c,q): the weight-enumerator of the set of words with a 1's b 2's c 3's `):
print(`according to the gep statistics q^(#321+#132+#213) done via a recurrence`):
print(`Try:`):
print(`Fabc(2,2,2,q);`):

elif nops([args])=1 and op(1,[args])=F1abc then
print(`F1abc(a,b,c,t32,t13,t21,q): the weight-enumerator of the set of words with a 1's b 2's c 3's done via the functional-recurrence equation)`):
print(`according to the Gepner statistics q^(#321+#132+#213)*t32^(#t32)*t13^(#13)*t21^(#21) done via a recurrence`):
print(`Try:`):
print(`F1abc(2,2,2,t32,t13,t21,q); `):

elif nops([args])=1 and op(1,[args])=F1abcD then
print(`F1abcD(a,b,c,t32,t13,t21,q): the weight-enumerator of the set of words with a 1's b 2's c 3's `):
print(`according to the Gepner statistics q^(#321+#132+#213)*t32^(#t32)*t13^(#13)*t21^(#21) done via a recurrence`):
print(`done directly (for checking purposes only!) `):
print(`Try:`):
print(`F1abcD(2,2,2,t32,t13,t21,q); `):



elif nops([args])=1 and op(1,[args])=gepPmoms then
print(`gepPmoms(n,r): the polynomial expressions for (i) the average (ii) the variance and the third through the r-th moments about the mean`):
print(`as polynomial expessions in n for the Gepner statistics on n-permutations. Warning: entire brute force, very slow`):
print(`gepPmoms(n,2);`):

elif nops([args])=1 and op(1,[args])=Hafokh then
print(`Hafokh(P,x,N): inputs a polynomial P whose coefficient of x^i is the factorial moment generating function outputs`):
print(`the sequence of moments about the mean up to the N-th moment. Try:`):
print(`Hafokh(((1+z)/2)^20,z,4);`):

elif nops([args])=1 and op(1,[args])=Gabc then
print(`Gabc(a,b,c,p,K): the weight-enumerator of the set of words with a 1's b 2's c 3's `):
print(`according to the gep statistics (1+p)^(#321+#132+#213) truncated to degree K, done via a recurrence`):
print(`Try: `):
print(` Gabc(2,2,2,p,3); `):

elif nops([args])=1 and op(1,[args])=G1abc then
print(`G1abc(a,b,c,x,y,z,p): the weight-enumerator of the set of words with a 1's b 2's c 3's done via the functional-recurrence equation)`):
print(`according to the Gepner statistics (1+p)^(#321+#132+#213)*(1+x)^(#t32)*(1+y)^(#13)*(1+z)^(#21) done via a recurrence`):
print(`Try:`):
print(`G1abc(2,2,2,x,y,z,p); `):

elif nops([args])=1 and op(1,[args])=G1abcT then
print(`G1abcT(a,b,c,x,y,z,p,K): the truncated, to degree K,`):
print(`weight-enumerator of the set of words with a 1's b 2's c 3's done via the functional-recurrence equation)`):
print(`according to the Gepner statistics (1+p)^(#321+#132+#213)*(1+x)^(#t32)*(1+y)^(#13)*(1+z)^(#21) done via a recurrence`):
print(`Try:`):
print(`G1abcT(2,2,2,x,y,z,p,2); `):


elif nops([args])=1 and op(1,[args])=gep then
print(`gep(w): The Gepner statistics on a permutation or word w. Try:`):
print(`gep([2,1,3,4]); `):

elif nops([args])=1 and op(1,[args])=GepMoms then
print(`GepMoms(n,r): the polynomial expressions for (i) the average (ii) the variance and the third through the r-th moments about the mean`):
print(` as polynomial expessions in n for the Gepner statistics on words with n occurences of each of 1,2,3. `):
print(` done via factorial moments and Gabc (i.e. by truncation). Try: `):
print(` GepMoms(n,2); `):

elif nops([args])=1 and op(1,[args])=GepMomsD then
print(`GepMomsD(n,r): the polynomial expressions for (i) the average (ii) the variance and the third through the r-th moments about the mean`):
print(` as polynomial expessions in n for the Gepner statistics on words with n occurences of each of 1,2,3. `):
print(` done directly (and hence slower than GepMoms(n,r). For checking purposes only. Try: `):
print(` GepMomsD(n,2); `):

elif nops([args])=1 and op(1,[args])=GepSeq then
print(`GepSeq(N,q): the first N terms of the Gepner polynomials Fabc(i,i,i,q), from i=1 to N. Try:`):
print(`GepSeq(6,q);`):

elif nops([args])=1 and op(1,[args])=GP then
print(`GP(L,n): tries to guess a polynomial of degree<=nops(L)-3 in n fitting L`):

elif nops([args])=1 and op(1,[args])=Milim then
print(`Milim(a,b,c): the set of words in {1,2,3} with a 1's, b 2's, and c 3's. Try`):
print(`Milim(2,2,2);`):

elif nops([args])=1 and op(1,[args])=Mishkal then
print(`Mishkal(w,t32,t13,t21,q): The weight of a word w in {1,2,3} defined by q^(#321+#132+#213)*t32^#32*t13^#13*t21^#21. Try:`):
print(` Mishkal([2,1,1,3],t32,t13,t21,q); `):

else
print(`There is no ezra for`,args):
fi:
 
end:

#GP1(L,n,d): inputs a sequence of rational numbers, a symbol n, and a non-neg.
#integer d, outputs a polynomial p(n) such that p(i)=L[i] for all i from 1 to nops(L)
GP1:=proc(L,n,d) local i,a,p,var, var1, eq:

if nops(L)<d+2 then
 print(`Make the list longer`):
 RETURN(FAIL):
fi:

p:=add(a[i]*n^i, i=0..d):

var:={ seq(a[i],i=0..d)}:

eq:={seq( subs(n=i, p)=L[i],i=1..d+1)}:

var1:=solve(eq,var):

if var1=NULL then

 RETURN(FAIL):
fi:


 p:=subs(var1,p):

if {seq(subs(n=i,p)-L[i],i=1..nops(L))}={0} then
 p:
else
 FAIL:
fi:

end:


#GP1t(L,n,d): inputs a sequence of rational numbers, a symbol n, and a non-neg.
#integer d, outputs a polynomial p(n) such that p(i)=L[i] for all i from 1 to nops(L)
GP1t:=proc(L,n,d) local i,a,p,var, var1, eq:

if nops(L)<d+1 then
 print(`Make the list longer`):
 RETURN(FAIL):
fi:

p:=add(a[i]*n^i, i=0..d):

var:={ seq(a[i],i=0..d)}:

eq:={seq( subs(n=i, p)=L[i],i=1..d+1)}:

var1:=solve(eq,var):

if var1=NULL then

 RETURN(FAIL):
fi:


 p:=subs(var1,p):

if {seq(subs(n=i,p)-L[i],i=1..nops(L))}={0} then
 p:
else
 FAIL:
fi:

end:


#GP(L,n): tries to guess a polynomial of degree<=nops(L)-3 in n fitting L
GP:=proc(L,n) local d,p:

for d from 0 to nops(L)-2 do
 p:=GP1(L,n,d):
 if p<>FAIL then
   RETURN(p):
 fi
od:

FAIL:

end:

#DoronSeq12(q): The pre-computed sequence of the first 12 Gepner polynomials
DoronSeq12:=proc(q):
[3*q+3, 6*q^8+21*q^6+36*q^4+21*q^2+6, 9*q^27+27*q^24+108*q^21+264*q^18+432*q^15
+432*q^12+264*q^9+108*q^6+27*q^3+9, 12*q^64+36*q^60+108*q^56+336*q^52+870*q^48+
2016*q^44+4041*q^40+6252*q^36+7308*q^32+6252*q^28+4041*q^24+2016*q^20+870*q^16+
336*q^12+108*q^8+36*q^4+12, 15*q^125+45*q^120+135*q^115+330*q^110+900*q^105+
2115*q^100+4710*q^95+10230*q^90+21195*q^85+40290*q^80+69423*q^75+102780*q^70+
126210*q^65+126210*q^60+102780*q^55+69423*q^50+40290*q^45+21195*q^40+10230*q^35
+4710*q^30+2115*q^25+900*q^20+330*q^15+135*q^10+45*q^5+15, 18*q^216+54*q^210+
162*q^204+396*q^198+918*q^192+2160*q^186+4788*q^180+9990*q^174+20250*q^168+
40563*q^162+79272*q^156+149769*q^150+276054*q^144+491886*q^138+830196*q^132+
1299720*q^126+1840536*q^120+2293326*q^114+2473020*q^108+2293326*q^102+1840536*q
^96+1299720*q^90+830196*q^84+491886*q^78+276054*q^72+149769*q^66+79272*q^60+
40563*q^54+20250*q^48+9990*q^42+4788*q^36+2160*q^30+918*q^24+396*q^18+162*q^12+
54*q^6+18, 21*q^343+63*q^336+189*q^329+462*q^322+1071*q^315+2268*q^308+4956*q^
301+10206*q^294+20412*q^287+39501*q^280+74781*q^273+139608*q^266+257397*q^259+
463827*q^252+822906*q^245+1442028*q^238+2493267*q^231+4228497*q^224+7018851*q^
217+11330508*q^210+17574606*q^203+25776663*q^196+35139888*q^189+43727985*q^182+
48966519*q^175+48966519*q^168+43727985*q^161+35139888*q^154+25776663*q^147+
17574606*q^140+11330508*q^133+7018851*q^126+4228497*q^119+2493267*q^112+1442028
*q^105+822906*q^98+463827*q^91+257397*q^84+139608*q^77+74781*q^70+39501*q^63+
20412*q^56+10206*q^49+4956*q^42+2268*q^35+1071*q^28+462*q^21+189*q^14+63*q^7+21
, 24*q^512+72*q^504+216*q^496+528*q^488+1224*q^480+2592*q^472+5304*q^464+10728*
q^456+21096*q^448+40248*q^440+75240*q^432+137520*q^424+246576*q^416+437400*q^
408+766872*q^400+1324752*q^392+2259804*q^384+3810888*q^376+6367428*q^368+
10550712*q^360+17305884*q^352+28091832*q^344+45167472*q^336+71848872*q^328+
112716510*q^320+173693976*q^312+261480672*q^304+381306840*q^296+532682913*q^288
+703991232*q^280+869087388*q^272+991121088*q^264+1036403964*q^256+991121088*q^
248+869087388*q^240+703991232*q^232+532682913*q^224+381306840*q^216+261480672*q
^208+173693976*q^200+112716510*q^192+71848872*q^184+45167472*q^176+28091832*q^
168+17305884*q^160+10550712*q^152+6367428*q^144+3810888*q^136+2259804*q^128+
1324752*q^120+766872*q^112+437400*q^104+246576*q^96+137520*q^88+75240*q^80+
40248*q^72+21096*q^64+10728*q^56+5304*q^48+2592*q^40+1224*q^32+528*q^24+216*q^
16+72*q^8+24, 27*q^729+81*q^720+243*q^711+594*q^702+1377*q^693+2916*q^684+5967*
q^675+11583*q^666+22437*q^657+42120*q^648+77598*q^639+140049*q^630+248886*q^621
+434160*q^612+747063*q^603+1271241*q^594+2142531*q^585+3569670*q^576+5889888*q^
567+9623448*q^558+15584724*q^549+25045524*q^540+39981114*q^531+63353016*q^522+
99669528*q^513+155742345*q^504+241849395*q^495+373299822*q^486+572455431*q^477+
871460208*q^468+1316250243*q^459+1969993791*q^450+2914546779*q^441+4247322777*q
^432+6068178666*q^423+8445890493*q^414+11360462187*q^405+14632926291*q^396+
17880194097*q^387+20549444877*q^378+22068832563*q^369+22068832563*q^360+
20549444877*q^351+17880194097*q^342+14632926291*q^333+11360462187*q^324+
8445890493*q^315+6068178666*q^306+4247322777*q^297+2914546779*q^288+1969993791*
q^279+1316250243*q^270+871460208*q^261+572455431*q^252+373299822*q^243+
241849395*q^234+155742345*q^225+99669528*q^216+63353016*q^207+39981114*q^198+
25045524*q^189+15584724*q^180+9623448*q^171+5889888*q^162+3569670*q^153+2142531
*q^144+1271241*q^135+747063*q^126+434160*q^117+248886*q^108+140049*q^99+77598*q
^90+42120*q^81+22437*q^72+11583*q^63+5967*q^54+2916*q^45+1377*q^36+594*q^27+243
*q^18+81*q^9+27, 30*q^1000+90*q^990+270*q^980+660*q^970+1530*q^960+3240*q^950+
6630*q^940+12870*q^930+24300*q^920+45090*q^910+81990*q^900+146070*q^890+256470*
q^880+443340*q^870+755730*q^860+1271370*q^850+2112930*q^840+3478050*q^830+
5675190*q^820+9174780*q^810+14709330*q^800+23397750*q^790+36932040*q^780+
57872520*q^770+90108750*q^760+139475265*q^750+214619130*q^740+328353225*q^730+
499600230*q^720+756154665*q^710+1138855200*q^700+1707397830*q^690+2548079460*q^
680+3785051340*q^670+5596720770*q^660+8237946645*q^650+12069383340*q^640+
17593605870*q^630+25498496970*q^620+36707592060*q^610+52419707076*q^600+
74101488885*q^590+103388691210*q^580+141823576980*q^570+190331660160*q^560+
248398892601*q^550+313112464830*q^540+378509540220*q^530+435889298910*q^520+
475564087230*q^510+489782289096*q^500+475564087230*q^490+435889298910*q^480+
378509540220*q^470+313112464830*q^460+248398892601*q^450+190331660160*q^440+
141823576980*q^430+103388691210*q^420+74101488885*q^410+52419707076*q^400+
36707592060*q^390+25498496970*q^380+17593605870*q^370+12069383340*q^360+
8237946645*q^350+5596720770*q^340+3785051340*q^330+2548079460*q^320+1707397830*
q^310+1138855200*q^300+756154665*q^290+499600230*q^280+328353225*q^270+
214619130*q^260+139475265*q^250+90108750*q^240+57872520*q^230+36932040*q^220+
23397750*q^210+14709330*q^200+9174780*q^190+5675190*q^180+3478050*q^170+2112930
*q^160+1271370*q^150+755730*q^140+443340*q^130+256470*q^120+146070*q^110+81990*
q^100+45090*q^90+24300*q^80+12870*q^70+6630*q^60+3240*q^50+1530*q^40+660*q^30+
270*q^20+90*q^10+30, 33*q^1331+99*q^1320+297*q^1309+726*q^1298+1683*q^1287+3564
*q^1276+7293*q^1265+14157*q^1254+26730*q^1243+48807*q^1232+88011*q^1221+155232*
q^1210+269742*q^1199+461439*q^1188+779823*q^1177+1300398*q^1166+2145231*q^1155+
3498759*q^1144+5650293*q^1133+9045729*q^1122+14369454*q^1111+22643181*q^1100+
35428734*q^1089+55052316*q^1078+84988662*q^1067+130374090*q^1056+198808434*q^
1045+301479486*q^1034+454855302*q^1023+682811415*q^1012+1020064155*q^1001+
1516842789*q^990+2245480182*q^979+3309707445*q^968+4858468032*q^957+7104448098*
q^946+10349809041*q^935+15022056555*q^924+21725209242*q^913+31309592028*q^902+
44970070464*q^891+64379070855*q^880+91864454352*q^869+130646851665*q^858+
185163978087*q^847+261492449988*q^836+367873957440*q^825+515326255983*q^814+
718301536227*q^803+995292706287*q^792+1369119307050*q^781+1866317716824*q^770+
2514870581487*q^759+3339309642075*q^748+4352211409722*q^737+5541912014499*q^726
+6858669928770*q^715+8204908310511*q^704+9437745326625*q^693+10390499432775*q^
682+10912156743147*q^671+10912156743147*q^660+10390499432775*q^649+
9437745326625*q^638+8204908310511*q^627+6858669928770*q^616+5541912014499*q^605
+4352211409722*q^594+3339309642075*q^583+2514870581487*q^572+1866317716824*q^
561+1369119307050*q^550+995292706287*q^539+718301536227*q^528+515326255983*q^
517+367873957440*q^506+261492449988*q^495+185163978087*q^484+130646851665*q^473
+91864454352*q^462+64379070855*q^451+44970070464*q^440+31309592028*q^429+
21725209242*q^418+15022056555*q^407+10349809041*q^396+7104448098*q^385+
4858468032*q^374+3309707445*q^363+2245480182*q^352+1516842789*q^341+1020064155*
q^330+682811415*q^319+454855302*q^308+301479486*q^297+198808434*q^286+130374090
*q^275+84988662*q^264+55052316*q^253+35428734*q^242+22643181*q^231+14369454*q^
220+9045729*q^209+5650293*q^198+3498759*q^187+2145231*q^176+1300398*q^165+
779823*q^154+461439*q^143+269742*q^132+155232*q^121+88011*q^110+48807*q^99+
26730*q^88+14157*q^77+7293*q^66+3564*q^55+1683*q^44+726*q^33+297*q^22+99*q^11+
33, 36*q^1728+108*q^1716+324*q^1704+792*q^1692+1836*q^1680+3888*q^1668+7956*q^
1656+15444*q^1644+29160*q^1632+53244*q^1620+95040*q^1608+166644*q^1596+287460*q
^1584+487836*q^1572+817560*q^1560+1353168*q^1548+2214864*q^1536+3587544*q^1524+
5754960*q^1512+9146088*q^1500+14411412*q^1488+22534740*q^1476+34987032*q^1464+
53945568*q^1452+82646280*q^1440+125850888*q^1428+190536300*q^1416+286876512*q^
1404+429643332*q^1392+640181988*q^1380+949354200*q^1368+1401517836*q^1356+
2060029800*q^1344+3015259956*q^1332+4395743532*q^1320+6383377368*q^1308+
9234911358*q^1296+13311747684*q^1284+19121540454*q^1272+27375736968*q^1260+
39067069230*q^1248+55576489176*q^1236+78822046008*q^1224+111459662280*q^1212+
157157269146*q^1200+220973148504*q^1188+309867399936*q^1176+433384307172*q^1164
+604573080996*q^1152+841234458816*q^1140+1167593474790*q^1128+1616531094540*q^
1116+2232529588566*q^1104+3075490504764*q^1092+4225641894207*q^1080+
5789806926876*q^1068+7909195714848*q^1056+10768563081660*q^1044+14606120253465*
q^1032+19722874911120*q^1020+26488879398834*q^1008+35341438519800*q^996+
46766568139830*q^984+61251628264668*q^972+79195888985256*q^960+100768591023084*
q^948+125717122046880*q^936+153157712093964*q^924+181423677004200*q^912+
208079161805040*q^900+230194620348606*q^888+244895878206288*q^876+
250058196381780*q^864+244895878206288*q^852+230194620348606*q^840+
208079161805040*q^828+181423677004200*q^816+153157712093964*q^804+
125717122046880*q^792+100768591023084*q^780+79195888985256*q^768+61251628264668
*q^756+46766568139830*q^744+35341438519800*q^732+26488879398834*q^720+
19722874911120*q^708+14606120253465*q^696+10768563081660*q^684+7909195714848*q^
672+5789806926876*q^660+4225641894207*q^648+3075490504764*q^636+2232529588566*q
^624+1616531094540*q^612+1167593474790*q^600+841234458816*q^588+604573080996*q^
576+433384307172*q^564+309867399936*q^552+220973148504*q^540+157157269146*q^528
+111459662280*q^516+78822046008*q^504+55576489176*q^492+39067069230*q^480+
27375736968*q^468+19121540454*q^456+13311747684*q^444+9234911358*q^432+
6383377368*q^420+4395743532*q^408+3015259956*q^396+2060029800*q^384+1401517836*
q^372+949354200*q^360+640181988*q^348+429643332*q^336+286876512*q^324+190536300
*q^312+125850888*q^300+82646280*q^288+53945568*q^276+34987032*q^264+22534740*q^
252+14411412*q^240+9146088*q^228+5754960*q^216+3587544*q^204+2214864*q^192+
1353168*q^180+817560*q^168+487836*q^156+287460*q^144+166644*q^132+95040*q^120+
53244*q^108+29160*q^96+15444*q^84+7956*q^72+3888*q^60+1836*q^48+792*q^36+324*q^
24+108*q^12+36]:



end:

###Start from HISTABRUT
#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:


#AlphaSeq(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AlphaSeq(((1+x)/2)^100,x,4);
AlphaSeq:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:
###end from HISTABRUT


#Milim(a,b,c): the set of words in {1,2,3} with a 1's, b 2's, and c 3's. Try
#Milim(2,2,2);
Milim:=proc(a,b,c) local w,gu1,gu2,gu3:
option remember:

if a=0 and b=0 and c=0 then
 RETURN({[]}):
fi:

if a<0 or b<0 or c<0 then
  RETURN({}):
fi:


gu1:=Milim(a-1,b,c):
gu2:=Milim(a,b-1,c):
gu3:=Milim(a,b,c-1):


{seq([op(w),1], w in gu1), seq([op(w),2], w in gu2), seq([op(w),3], w in gu3)}:
end:

#Mishkal(w,t32,t13,t21,q): The weight of a word w in {1,2,3} defined by q^(#321+#132+#213)*t32^#32*t13^#13*t21^#21. Try:
#Mishkal([2,1,1,3],t32,t13,t21,q);
Mishkal:=proc(w,t32,t13,t21,q) local gu,i,j,k:
gu:=1:

for i from 1 to nops(w)-1 do
  for j from i+1 to nops(w) do
    if w[i]=3 and w[j]=2 then
      gu:=gu*t32:
     elif w[i]=1 and w[j]=3 then
       gu:=gu*t13:
     elif w[i]=2 and w[j]=1 then
       gu:=gu*t21:
    fi:
 od:
od:

for i from 1 to nops(w)-2 do
 for j from i+1 to nops(w)-1 do
   for k from j+1 to nops(w) do
      if ( (w[i]=1 and w[j]=3 and w[k]=2) or (w[i]=2 and w[j]=1 and w[k]=3) or (w[i]=3 and w[j]=2 and w[k]=1)) then
       gu:=gu*q:
      fi:
   od:
 od:
od:
    

gu:

end:



#F1abcD(a,b,c,t32,t13,t21,q): the weight-enumerator of the set of words with a 1's b 2's c 3's done directly (for checking purposes only!)
#Try:
#F1abcD(2,2,2,t32,t13,t21,q);
F1abcD:=proc(a,b,c,t32,t13,t21,q) local gu, w:
gu:=Milim(a,b,c):
add(Mishkal(w,t32,t13,t21,q), w in gu):
end:


#F1abc(a,b,c,t32,t13,t21,q): the weight-enumerator of the set of words with a 1's b 2's c 3's done via a recurrence
#Try:
#F1abc(2,2,2,t32,t13,t21,q);
F1abc:=proc(a,b,c,t32,t13,t21,q)
option remember:

if a=0 and b=0 and c=0 then
  RETURN(1):
fi:

if a<0 or b<0 or c<0 then
  RETURN(0):
fi:

expand(
t21^b*subs(t32=q*t32,F1abc(a-1,b,c,t32,t13,t21,q))+
t32^c*subs(t13=q*t13,F1abc(a,b-1,c,t32,t13,t21,q))+
t13^a*subs(t21=q*t21,F1abc(a,b,c-1,t32,t13,t21,q))
):

end:


#Fabc(a,b,c,q): the weight-enumerator of the set of words with a 1's b 2's c 3's 
#according to the gep statistics q^(#321+#132+#213) done via a recurrence
#Try:
#Fabc(2,2,2,q);
Fabc:=proc(a,b,c,q) local gu,t32,t13,t21:
gu:=F1abc(a,b,c,t32,t13,t21,q):
subs({t32=1,t13=1,t21=1},gu):
end:


#GepSeq(N,q): the first N terms of the Gepner polynomials Fabc(i,i,i,q), from i=1 to N. Try:
#GepSeq(6,q);
GepSeq:=proc(N,q) local gu,i:
if N<=12 then
 gu:=DoronSeq12(q):
 RETURN([op(1..N,gu)]):
else
 RETURN([seq(Fabc(i,i,i,q),i=1..N)]):
fi:
end:





#Chop(P,var,K): chops the polynomial P of the list of variables var to its terms of (total) degree<=K
#Try: Chop((x+y)^5,[x,y],2);
Chop:=proc(P,var,K) local i,t,lu:
lu:=subs({seq(var[i]=t*var[i],i=1..nops(var))},P):
lu:=expand(add(coeff(lu,t,i)*t^i,i=0..K)):
subs(t=1,lu):
end:



#G1abc(a,b,c,x,y,z,p): the weight-enumerator of the set of words with a 1's b 2's c 3's done via a recurrence
#according to the gep statistics (1+x)*q^(#321+#132+#213) done via a recurrence`):
#Try:
#G1abc(2,2,2,x,y,z,p);
G1abc:=proc(a,b,c,x,y,z,p)
option remember:

if a=0 and b=0 and c=0 then
  RETURN(1):
fi:

if a<0 or b<0 or c<0 then
  RETURN(0):
fi:

expand(
(1+z)^b*subs(x=p+x+p*x,G1abc(a-1,b,c,x,y,z,p))+
(1+x)^c*subs(y=p+y+p*y,G1abc(a,b-1,c,x,y,z,p))+
(1+y)^a*subs(z=p+z+p*z,G1abc(a,b,c-1,x,y,z,p))
):

end:


#G1abcT(a,b,c,x,y,z,p,K): the truncated, to degree K, weight-enumerator of the set of words with a 1's b 2's c 3's done via a recurrence
#according to the gep statistics (1+x)*q^(#321+#132+#213) done via a recurrence`):
#Try:
#G1abcT(2,2,2,x,y,z,p,2);
G1abcT:=proc(a,b,c,x,y,z,p,K)
option remember:

if a=0 and b=0 and c=0 then
  RETURN(1):
fi:

if a<0 or b<0 or c<0 then
  RETURN(0):
fi:

Chop(
(1+z)^b*subs(x=p+x+p*x,G1abcT(a-1,b,c,x,y,z,p,K))+
(1+x)^c*subs(y=p+y+p*y,G1abcT(a,b-1,c,x,y,z,p,K))+
(1+y)^a*subs(z=p+z+p*z,G1abcT(a,b,c-1,x,y,z,p,K)),[x,y,z,p],K):

end:

#Gabc(a,b,c,p,K): the weight-enumerator of the set of words with a 1's b 2's c 3's 
#according to the gep statistics (1+p)^(#321+#132+#213) truncated to degree K, done via a recurrence
#Try:
#Gabc(2,2,2,p,3);
Gabc:=proc(a,b,c,p,K) local gu,x,y,z:
gu:=G1abcT(a,b,c,x,y,z,p,K):
gu:=subs({x=0,y=0,z=0},gu):
end:


#GepMomsD(n,r): the polynomial expressions for (i) the average (ii) the variance and the third through the r-th moments about the mean
#as polynomial expessions in n for the Gepner statistics on words with n occurences of each of 1,2,3.
#done directly (and hence slower than GepMoms(n,r). For checking purposes only. Try:
#GepMomsD(n,2);
GepMomsD:=proc(n,r) local i,q,gu,mu,j,S,lu,ho:
 gu:=GepSeq(2*r+2,q):

S:=[]:

mu:=[seq(AveAndMoms(gu[i],q,r),i=1..nops(gu))]:

S:=[]:

for j from 1 to r do
 lu:=[seq(mu[i][j],i=1..nops(mu))]:
 ho:=GP(lu,n):
  if ho=FAIL then
    RETURN(S):
  else
   S:=[op(S),ho]:
  fi:
od:

factor(S):
 
end:


#Hafokh(P,x,N): inputs a polynomial P whose coefficient of x^i is the factorial moment generating function outputs
#the sequence of moments about the mean. Try:
#Hafokh(((1+z)/2)^20,z);
Hafokh:=proc(P,x,N) local gu,P1,i,j,av:
P1:=P/subs(x=0,P):
gu:=[seq(i!*coeff(P1,x,i),i=1..N)]:

gu:=[seq(add(stirling2(i,j)*gu[j],j=1..i),i=1..N)]:

av:=gu[1]:

[av,seq(av^i+add((-1)^j*binomial(i,j)*gu[j]*av^(i-j),j=1..i),i=2..N)]:

#[ave1, seq(simplify(add(stirling2(i,j)*gu[j],j=1..i)),i=2..N)]:

end:


#GepMoms(n,r): the polynomial expressions for (i) the average (ii) the variance and the third through the r-th moments about the mean
#as polynomial expessions in n for the Gepner statistics on words with n occurences of each of 1,2,3.
#done via factorial moments and the truncated version. Try:
#GepMoms(n,4);
GepMoms:=proc(n,r) local p,i,gu,mu,j,S,lu,ho:
gu:=[seq(Gabc(i,i,i,p,r),i=1..2*r+2)]:
gu:=[seq(gu[i]/coeff(gu[i],p,0),i=1..nops(gu))]:
mu:=[seq(Hafokh(gu[i],p,r),i=1..nops(gu))]:

S:=[]:

for j from 1 to r do
 lu:=[seq(mu[i][j],i=1..nops(mu))]:
 ho:=GP(lu,n):
  if ho=FAIL then
    RETURN(S):
  else
   S:=[op(S),ho]:
  fi:
od:

factor(S):


end:






#gep(w): The Gepner statistics on a permutation or word p
#try:
#gep([2,1,3,4]);
gep:=proc(w) local n,i,j,k,co:
n:=nops(w):
co:=0:

for i from 1 to n do
 for j from i+1 to n do
  for k from j+1 to n do
   if ( (w[i]<w[k] and w[k]<w[j]) or (w[k]<w[j] and w[j]<w[i]) or (w[j]<w[i] and w[i]<w[k]) ) then
    co:=co+1:
   fi:
  od:
 od:
od:
co:
end:


#DGn(q,n): The weight-enumerator of the set of n-permutations according to the variable q and the gep-statistics. Try:
#DGn(q,5);
DGn:=proc(q,n) local gu,w:
gu:=permute(n):
add(q^gep(w), w in gu):
end:

#DGseq(q,N): The sequence of the the weight-enumerator of the set of n-permutations according to the variable q and the gep-statistics
#from n=1 to n=N
#DGseq(q,5);

DGseq:=proc(q,N) local n:
[seq(DGn(q,n),n=1..N)]:
end:

#gepPmoms(n,r): the polynomial expressions for (i) the average (ii) the variance and the third through the r-th moments about the mean
#as polynomial expessions in n for the Gepner statistics on n-permutations. Warning: entire brute force, very slow
#gepPmoms(n,2);
gepPmoms:=proc(n,r) local i,q,gu,mu,j,S,lu,ho:
 gu:=DGseq(q,2*r+2):

S:=[]:

mu:=[seq(AveAndMoms(gu[i],q,r),i=1..nops(gu))]:

S:=[]:

for j from 1 to r do
 lu:=[seq(mu[i][j],i=1..nops(mu))]:
 ho:=GP(lu,n):
  if ho=FAIL then
    RETURN(S):
  else
   S:=[op(S),ho]:
  fi:
od:

factor(S):
 
end: