##
## 12/28/2011
##
with(qseries):
 
 
CRANKB:=proc(z,q,K,T)
	local y,n:
	y:=0:
	for n from  -T to T do
	   y:=y + (-1)^n*q^(K*n*(n+3)/2)/(1-z*q^n);
	od:
	RETURN(y):
end:

 
 
SIGk:=proc(z,q,K,T)
	local y,n:
	y:=0:
	for n from  -T to T do
	   y:=y + (-1)^n*q^(K*n*(n+1)/2)/(1-z*q^n);
	od:
	RETURN(y):
end:
 
 
Sk:=proc(zeta,z,q,K,T)
	local y,n,y1:
	y:=0:
	for n from  -T to T do
	   y1:=(-1)^n*q^(K*n*(n+1)/2)*(
                zeta^(-K*n)/(1-z/zeta*q^n) +
                zeta^(K*n+K)/(1-z*zeta*q^n)): 
           y:=y+normal(y1):
	od:
	RETURN(y):
end:
with(numtheory):

 
 
EISEN:=proc(k,q,T)
	local y,n:
	y:=0:
	for n from 1 to T do
	y:=y+sigma[k](n)*q^n:
	od:
RETURN(y):
end:

 
 
POWERq:=proc(a,N,M)  #compute a^n to O(q^M)
        local j,Q,Y,n,b:
        j:=1:
        Q:=1:
        Y:=1:
        n:=N:
        while n>0 do
print(n);
          if modp(n,2)=1 then
            b:=1:
            if n=N then
                Q:=a: Y:=a:
            else
                Q:=normal(series(Q^2,q,M)): Y:=normal(series(Q*Y,q,M)):
            fi:
          else
            b:=0:
            if n=N then
                Q:=a:
            else
                Q:=normal(series(Q*Q,q,M)):
            fi:
          fi:
          n:=trunc(n/2):
        od:
        RETURN(Y);
end:


ZOP:=X->z*diff(X,z):
XOP:=X->x*diff(X,x):
QOP:=X->q*diff(X,q):
ZETAOP:=X->diff(X,zeta):
NEWZETAOP:=X->zeta*diff(X,zeta):
ZETAOP2:=X->zeta*diff(X,zeta):
#######################################################################
 
 
HkOP:=proc(X,k)
	local ans:
	ans:=k*z*diff(X,z) + 2*k*q*diff(X,q) + z*diff(z*diff(X,z),z):
	RETURN(ans):
end:
 
 
HkOPS:=proc(X,k,n)
	local a,m:
	a:=X:
	for m from 1 to n do
	  a:=HkOP(a,k):
	od:
	RETURN(a):
end:
#######################################################################
 
 
ZOPS:=proc(X,n)
	local a,k:
	a:=X:
	for k from 1 to n do
	a:=ZOP(a):
	od:
	RETURN(a):
end:
 
 
XOPS:=proc(X,n)
	local a,k:
	a:=X:
	for k from 1 to n do
	a:=XOP(a):
	od:
	RETURN(subs(x=1,a)):
end:
 
 
ZETAOPS:=proc(X,n)
	local a,k:
	a:=X:
	for k from 1 to n do
	a:=ZETAOP(a):
	od:
	RETURN(subs(zeta=1,a)):
end:
 
NEWZETAOPS:=proc(X,n)
	local a,k:
	a:=X:
	for k from 1 to n do
	a:=normal(NEWZETAOP(a)):
	od:
	RETURN(subs(zeta=1,a)):
end:
 
ZETAOPS2:=proc(X,n)
	local a,k:
	a:=X:
	for k from 1 to n do
	a:=ZETAOP2(a):
	od:
	RETURN(subs(zeta=1,a)):
end:
 
 
QOPS:=proc(Q,n)
	local a,k:
	a:=Q:
	for k from 1 to n do
	a:=QOP(a):
	od:
	RETURN(a):
end:
###############################################################
with(combinat):
 
 
PellPOLY:=proc(ell,x,k)
	local ans,j,m;
        ans:=0:
        for j from 1 to ell do
          for m from 0 to trunc(j/2) do
              ans:= ans + 
              stirling1(ell,j)*j*(j-m-1)!/(j-2*m)!/m!*x^m*k^(j-2*m):
 	  od:
        od:
	RETURN(ans):
end:      
###############################################################
 
 
checkDSPHid:=proc(ell,k,T,TT)
local L1,R1,poly,dg,r,j,cj;
   global L2,R2:
   L1:=normal(series(Sk(zeta,z,q,k,T),q,TT)):
   L1:=convert(L1,polynom):
   L2:=ZETAOPS(L1,ell):
   R1:=normal(series(SIGk(z,q,k,T),q,TT)):
   poly:=PellPOLY(ell,x,k):
   dg:=degree(poly,x):
   r:=coeff(poly,x,0)*R1:
   lprint("POLY has degree ",dg);
   lprint("POLY:",poly);
   for j from 1 to dg do
       cj:=coeff(poly,x,j):
       r := r + cj*HkOPS(R1,k,j):
   od:
   R2:=r:
   RETURN(L2-R2):
end:

 
findPDE:=proc(k,T,TT)
local C1,j,FJ,kj,m,wt,func,sfunc;
   global R1,C2,FUNCS,SFUNCS:
   C1:=normal(series(SIGk(z,q,1,TT),q,TT)):
   C2:=POWERq(C1,k,TT):
   lprint("C2 computed");
   R1:=normal(series(SIGk(z,q,k,T),q,TT)):
   makequasilist(k-1,TT): ## quasimodular forms of wt <= N
   FUNCS:=[]:
   SFUNCS:=[]:
   for j from 1 to nops(MLIST) do
      FUNCS:=[op(FUNCS),MLIST[j]*R1]:
      SFUNCS:=[op(SFUNCS),MSLIST[j]]:
   od:
   for j from 1 to (k-1)/2 do
      FJ:=HkOPS(R1,k,j):
      kj:=(k-1-j):
      if modp(kj,2)=1 then
        kj:=kj-1:
      fi:
      for m from 1 to nops(MLIST) do
        wt:=WLIST[m]:
        if wt<=kj then
           func:=MLIST[m]*FJ:
           sfunc:=MSLIST[m]*HK^j:
           FUNCS:=[op(FUNCS),func]:
           SFUNCS:=[op(SFUNCS),sfunc]:
        fi:
      od:
   od:
   lprint("number of funcs =",nops(FUNCS));
end:
 
findPDE2:=proc(bump,p,k,T,TT)
local C1,j,FJ,kj,m,wt,func,sfunc;
   global R1,C2,FUNCS,SFUNCS:
   C1:=normal(series(SIGk(z,q,1,TT),q,TT)):
   C2:=POWERq(C1,k,TT):
   C2:=normal(series(C1/etaq(q,1,TT)^p,q,TT)):
   lprint("C2 computed");
   R1:=normal(series(SIGk(z,q,k,T)/etaq(q,1,TT)^p,q,TT)):
   makequasilist(k-1+bump,TT): ## quasimodular forms of wt <= N
   FUNCS:=[]:
   SFUNCS:=[]:
   for j from 1 to nops(MLIST) do
      FUNCS:=[op(FUNCS),MLIST[j]*R1]:
      SFUNCS:=[op(SFUNCS),MSLIST[j]]:
   od:
   for j from 1 to (k-1)/2 do
      FJ:=HkOPS(R1,k,j):
      kj:=(k-1-j):
      if modp(kj,2)=1 then
        kj:=kj-1:
      fi:
      for m from 1 to nops(MLIST) do
        wt:=WLIST[m]:
        if wt<=kj then
           func:=MLIST[m]*FJ:
           sfunc:=MSLIST[m]*HK^j:
           FUNCS:=[op(FUNCS),func]:
           SFUNCS:=[op(SFUNCS),sfunc]:
        fi:
      od:
   od:
   lprint("number of funcs =",nops(FUNCS));
end:


###############################################################

 
 
makequasilist:=proc(N,TT) ## quasimodular forms of wt <= N
local a,b,c,func,sfunc,phi1,phi3,phi5;
    global MLIST,MSLIST,WLIST:
   phi1:=add(sigma[1](n)*q^n,n=1..TT):
   phi3:=add(sigma[3](n)*q^n,n=1..TT):
   phi5:=add(sigma[5](n)*q^n,n=1..TT):
    MLIST:=[]:
    MSLIST:=[]:
    WLIST:=[]:
    for a from 0 to N/2 do
    for b from 0 to (N/2 - a)/2 do
    for c from 0 to (N/2 - a - 2*b)/3 do
      func:=phi1^a*phi3^b*phi5^c:
      sfunc:=PHI[1]^a*PHI[3]^b*PHI[5]^c:
      MLIST:=[op(MLIST),func]:
      MSLIST:=[op(MSLIST),sfunc]:
      WLIST:=[op(WLIST),2*a+3*b+6*c]:
    od:
    od:
    od:
RETURN():
end:




## PHILIST:=[]:
## PHISLIST:=[]:
## for a from 0 to 12 do
## for b from 0 to 6 do
## for c from 0 to 4 do
##   wt:=2*a+4*b+6*c:
##   if wt>0 and wt<=12 then
##     print("a ,b ,c =",a,b,c);
##     phifunc:=series(PHIQ(1)^a*PHIQ(3)^b*PHIQ(5)^c,q,TT):
##     PHILIST:=[op(PHILIST),phifunc]:
##     PHISLIST:=[op(PHISLIST), PHI[1]^a*PHI[3]^b*PHI[5]^c]:
##   fi:
## od:
## od:
## od:



 
 
makediffRlist:=proc(K,N,TT)
local a,b,c,xx,XX;
global FLIST,FSLIST,R1:
FLIST:=[]:
FSLIST:=[]:
R1:=normal(series(CRANK(z,q,K,20),q,TT)):
R1:=convert(R1,polynom):
for a from 0 to N do
for b from 0 to N/2 do
    c:=a+b:
if c<=N then
   print("a = ",a," b = ",b);
     xx:=ZOPS(R1,a):
     xx:=QOPS(xx,b):
     XX:=normal(series(xx,q,TT)):
     FLIST:=[op(FLIST),XX]:
     FSLIST:=[op(FSLIST),DELZ^a*DELQ^b]:
fi:
od:
od:
RETURN():
end:



##D_{k-1}\, S_k(\zeta,z,q)
##= 2 \left( \tfrac{1}{2}(k-1)k + \mathcal{H}_k\right)
##  \prod_{a=2}^{\frac{1}{2}(k-1)} \left(a(k-a) +\mathcal{H}_k\right) \,
##\Sigma^{(k)}(z,q).
## k\delta_z + 2k \delta_q + \delta_z^2
HH:=k->k*DELZ + 2*k*DELQ+DELZ^2:
V:=(k,j)->2*(k*(k-1)/2 + HH(k))*mul(a*(k-a)+HH(k),a=2..j):

## NEWPOLY:=(x,k,ell)->add(ell*(ell-m-1)!/(ell-2*m)!/m!*k^(ell-2*m)*x^m,
##          m=0..floor(ell/2)):
NEWPOLY:=proc(x,k,ell)
 if ell=0 then
    2:
 else
    add(ell*(ell-m-1)!/(ell-2*m)!/m!*k^(ell-2*m)*x^m,
    m=0..floor(ell/2)):
 fi:
end:
NEWPOLY0:=(x,k,ell)->expand( (k/2-1/2*sqrt(k^2+4*x))^ell
                            + (k/2+1/2*sqrt(k^2+4*x))^ell):

checkNEWDSPHid:=proc(ell,k,T,TT)
local L1,R1,poly,dg,r,j,cj;
   global L2,R2:
   L1:=normal(series(Sk(zeta,z,q,k,T),q,TT)):
   L1:=convert(L1,polynom):
   L2:=NEWZETAOPS(L1,ell):
   R1:=normal(series(SIGk(z,q,k,T),q,TT)):
   R2:=k^ell*R1 + 
       add(ell*(ell-m-1)!/(ell-2*m)!/m!*k^(ell-2*m)*HkOPS(R1,k,m),
       m=1..floor(ell/2)):
   RETURN(L2-R2):
end:

F1:=(zeta,q,m,T)->tripleprod(zeta^(m+1),q,T)/tripleprod(zeta^m,q,T)*zeta^m:
G1:=(zeta,m)->normal(m + m*zeta^m/(1-zeta^m) - (m+1)*zeta^(m+1)/(1-zeta^(m+1))):

DjL:=proc(j,m)
     if modp(j,2)=0 then
        NEWZETAOPS(G1(zeta,m),j):
     else
        NEWZETAOPS(G1(zeta,m),j) + 2*( m^(j+1) - (m+1)^(j+1))*PHI[j]:
     fi:
end:

## DaF1:=proc(a,m)
##     option remember:
##     if a=0 then
##       (m+1)/m:
##     else
##       xx:=add(binomial(a-1,2*j-1)*DjL(2*j-1,m)*_DaF1(a-2*j,m),j=1..trunc(a/2)):
##       + _DaF1(a-1,m)*(m+1/2):
##     fi:
## end:
##DaF1:=proc(a,m)
##    option remember:
##    if a=0 then
##      (m+1)/m:
##    else
##      xx:=add(binomial(a-1,j)*DjL(j,m)*DaF1(a-1-j,m),j=0..(a-1)):          
##    fi:
##end:
DaF1:=proc(a,m)
    local xx:
    option remember:
    if a=0 then
      (m+1)/m:
    else
      ##xx:= DjL(0,m)*DaF1(a-1,m)+
      xx:= (m + 1/2)*DaF1(a-1,m)+
           add(binomial(a-1,2*j-1)*DjL(2*j-1,m)*DaF1(a-2*j,m),j=1..floor(a/2)):  
    fi:
end:
phi1:=add(sigma[1](n)*q^n,n=1..100):
phi3:=add(sigma[3](n)*q^n,n=1..100):
phi5:=add(sigma[5](n)*q^n,n=1..100):
phi7:=add(sigma[7](n)*q^n,n=1..100):
phi9:=add(sigma[9](n)*q^n,n=1..100):
phi11:=add(sigma[11](n)*q^n,n=1..100):
phi13:=add(sigma[13](n)*q^n,n=1..100):
   
#++++++++++++++++++++++++++++++++++++++++++++++++++++

G:=(m,j,zeta)->normal(add( (m-k)*zeta^(k-m)/(1-zeta^(k-m)) 
+ (m+k+1)*zeta^(m+k+1)/(1-zeta^(m+k+1)),k=1..j)):


#++++++++++++++++++++++++++++++++++++++++++++++++++++


DaLj:=proc(a,m,j)
     if modp(a,2)=0 then
        NEWZETAOPS(G(m,j,zeta),a):
     else
        NEWZETAOPS(G(m,j,zeta),a) + 
        2*add( (m+k+1)^(a+1)-(m-k)^(a+1),k=1..j)*PHI[a]:
     fi:
end:
DaFj:=proc(a,m,j)
    local xx:
    option remember:
    if a=0 then
      (-1)^j*mul( (m-k)/(m+k+1), k=1..j):
    else
      xx:= DaLj(0,m,j)*DaFj(a-1,m,j)+
   add(binomial(a-1,2*k-1)*DaLj(2*k-1,m,j)*DaFj(a-2*k,m,j),k=1..floor(a/2)):  
    fi:
end:

Fj:=(zeta,q,m,j,T)->
mul(tripleprod(zeta^(-(m-k)),q,T)/tripleprod(zeta^(m+k+1),q,T),
   k=1..j):

#############################################################
## 12/28/2011: New versions of DaF0m etc
symbDaF0m:=proc(a,m)
     local xx:
     option remember:
     if a=0 then
        (m+1)/m:
     else
        xx:=(m+ 1/2)*symbDaF0m(a-1,m) +
        add(2*binomial(a-1,2*i-1)*(m^(2*i)-(m+1)^(2*i))*G[2*i]*symbDaF0m(a-2*i,m),
        i=1..floor(a/2)):
     fi:
end:
G2Phi:=a->-bernoulli(a)/a/2+Phi[a-1]:
G2Phiq:=a->-bernoulli(a)/a/2+EISEN(a-1,q,100):
symbDaF0mv2:=proc(a,m)
     local xx:
     option remember:
     if a=0 then
        (m+1)/m:
     else
        xx:=(m+ 1/2)*symbDaF0mv2(a-1,m) +
    add(2*binomial(a-1,2*i-1)*(m^(2*i)-(m+1)^(2*i))*G2Phi(2*i)*symbDaF0mv2(a-2*i,m),
        i=1..floor(a/2)):
     fi:
end:
qDaF0m:=proc(a,m)
     local xx:
     option remember:
     if a=0 then
        (m+1)/m:
     else
        xx:=(m+ 1/2)*qDaF0m(a-1,m) +
        add(2*binomial(a-1,2*i-1)*(m^(2*i)-(m+1)^(2*i))*G2Phiq(2*i)*qDaF0m(a-2*i,m),
        i=1..floor(a/2)):
     fi:
end:
symbDaFjm:=proc(a,j,m)
     local xx,i,k:
     option remember:
     if a=0 then
        (-1)^j*mul((m-i)/(m+i+1),i=1..j):
     else
        xx:=-j*(m+ 1/2)*symbDaFjm(a-1,j,m) +
        add(
        add(2*binomial(a-1,2*i-1)*((m+k+1)^(2*i)-(m-k)^(2*i))*G[2*i]
            *symbDaFjm(a-2*i,j,m), k=1..j),
        i=1..floor(a/2)):
     fi:
end:
symbDaFjmv2:=proc(a,j,m)
     local xx,i,k:
     option remember:
     if a=0 then
        (-1)^j*mul((m-i)/(m+i+1),i=1..j):
     else
        xx:=-j*(m+ 1/2)*symbDaFjmv2(a-1,j,m) +
        add(
        add(2*binomial(a-1,2*i-1)*((m+k+1)^(2*i)-(m-k)^(2*i))*G2Phi(2*i)
            *symbDaFjmv2(a-2*i,j,m), k=1..j),
        i=1..floor(a/2)):
     fi:
end:
qDaFjm:=proc(a,j,m)
     local xx,i,k:
     option remember:
     if a=0 then
        (-1)^j*mul((m-i)/(m+i+1),i=1..j):
     else
        xx:=-j*(m+ 1/2)*qDaFjm(a-1,j,m) +
        add(
        add(2*binomial(a-1,2*i-1)*((m+k+1)^(2*i)-(m-k)^(2*i))*G2Phiq(2*i)
            *qDaFjm(a-2*i,j,m), k=1..j),
        i=1..floor(a/2)):
     fi:
end:
#############################################################
CKCOF:=m->(-1)^(m+1)*(m+1)!*(m-1)!*(2*m)!;
#############################################################
SYMBOLPDEGENv0:=proc(k)
     local m,LHS,RHS,RHS1,RHS2,RHS3:
     m:=1/2*(k-1):
     LHS:=CKCOF(m)*CS^k*QINF^k:
     RHS1:=NEWPOLY(_H[k],k,2*m):
     RHS2:=add(add(binomial(2*m,a)*DaFj(a,m,j)*(j+1)^(2*m-a)*NEWPOLY(_H[k],k,2*m-a),
           a=0..2*m),j=1..m-1):
     RHS3:=-DaF1(2*m,m):
     LHS=collect(expand(RHS1+RHS2+RHS3),_H[k]):
end:
##symbDaF0mv2:=proc(a,m)
##symbDaFjmv2:=proc(a,j,m)
SYMBOLPDEGEN:=proc(k)
     local m,LHS,RHS,RHS1,RHS2,RHS3:
     m:=1/2*(k-1):
     LHS:=CKCOF(m)*CS^k*QINF^k:
     RHS1:=NEWPOLY(_H[k],k,2*m):
     RHS2:=add(add(binomial(2*m,a)*symbDaFjmv2(a,j,m)*(j+1)^(2*m-a)*NEWPOLY(_H[k],k,2*m-a),
           a=0..2*m),j=1..m-1):
     RHS3:=-symbDaF0mv2(2*m,m):
     LHS=collect(expand(RHS1+RHS2+RHS3),_H[k]):
end:
SYMBOLPDEGEN2v0:=proc(k)
     local m,LHS,RHS,RHS1,RHS2,RHS3:
     m:=1/2*(k-1):
     LHS:=CKCOF(m)*CS^k*QINF^k:
     RHS1:=NEWPOLY(_H[k],k,2*m):
     RHS2:=add(add(binomial(2*m,a)*_DaFj(a,m,j)*(j+1)^(2*m-a)*NEWPOLY(_H[k],k,2*m-a),
           a=0..2*m),j=1..m-1):
     RHS3:=-_DaF1(2*m,m):
     LHS=RHS1 + RHS2 + RHS3:                       
end:
SYMBOLPDEGEN2:=proc(k)
     local m,LHS,RHS,RHS1,RHS2,RHS3:
     m:=1/2*(k-1):
     LHS:=CKCOF(m)*CS^k*QINF^k:
     RHS1:=NEWPOLY(_H[k],k,2*m):
     RHS2:=add(add(binomial(2*m,a)*_D[a](_F[j,m])*(j+1)^(2*m-a)*NEWPOLY(_H[k],k,2*m-a),
           a=0..2*m),j=1..m-1):
     RHS3:=-_D[2*m](_F[0,m]):
     LHS=RHS1 + RHS2 + RHS3:                       
end:
SYMBOLPDEGEN3v0:=proc(k)
     local m,LHS,RHS,RHS1,RHS2,RHS3:
     m:=1/2*(k-1):
     LHS:=CKCOF(m)*CS^k*QINF^k:
     RHS1:=NEWPOLY(_H[k],k,2*m):
     RHS2:=add(add(binomial(2*m,a)*{expand(DaFj(a,m,j))}*(j+1)^(2*m-a)*NEWPOLY(_H[k],k,2*m-a),
           a=0..2*m),j=1..m-1):
     RHS3:=-[expand(DaF1(2*m,m))]:
     LHS=RHS1 + RHS2 + RHS3:                       
end:
SYMBOLPDEGEN3v1:=proc(k)
     local m,LHS,RHS,RHS1,RHS2,RHS3:
     m:=1/2*(k-1):
     LHS:=CKCOF(m)*CS^k*QINF^k:
     RHS1:=NEWPOLY(_H[k],k,2*m):
     RHS2:=add(add({expand(binomial(2*m,a)*expand(DaFj(a,m,j))*(j+1)^(2*m-a))}*NEWPOLY(_H[k],k,2*m-a),
           a=0..2*m),j=1..m-1):
     RHS3:=-{expand(DaF1(2*m,m))}:
     LHS=RHS1 + RHS2 + RHS3:                       
end: 
##symbDaF0mv2:=proc(a,m)
##symbDaFjmv2:=proc(a,j,m)
SYMBOLPDEGEN3:=proc(k)
     local m,LHS,RHS,RHS1,RHS2,RHS3:
     m:=1/2*(k-1):
     LHS:=CKCOF(m)*CS^k*QINF^k:
     RHS1:=NEWPOLY(_H[k],k,2*m):
     RHS2:=add(add({expand(binomial(2*m,a)*expand(symbDaFjmv2(a,j,m))*(j+1)^(2*m-a))}*NEWPOLY(_H[k],k,2*m-a),
           a=0..2*m),j=1..m-1):
     RHS3:=-{expand(symbDaF0mv2(2*m,m))}:
     LHS=RHS1 + RHS2 + RHS3:                       
end: 
## SYMBOLPDEGEN(3);
##                         3     3                          
##                     4 CS  QINF  = 4 + 2 _H[3] + 12 PHI[1]
## 
checkPDE:=proc(K,T,TT)
local C1,j,FJ,kj,m,wt,func,sfunc,PDE,LS,RS,LRS,cofj,XX;
   global R1,C2,FUNCS,SFUNCS:
   PDE:=SYMBOLPDEGEN(K):
   print(PDE);
   C1:=normal(series(SIGk(z,q,1,TT),q,TT)):
   C2:=POWERq(C1,K,TT):
   lprint("C2 computed");
   R1:=normal(series(SIGk(z,q,K,T),q,TT)):
   m:=(K-1)/2:
   LS:=series(CKCOF(m)*C2,q,TT):
   RS:=0:
   LRS:=degree(rhs(PDE),_H[K]):
   for j from 0 to LRS do
      cofj:=coeff(rhs(PDE),_H[K],j):
      RS := RS + cofj*HkOPS(R1,K,j):
   od:
##   XX:=normal(series( subs({PHI[1]=phi1,PHI[3]=phi3,PHI[5]=phi5,PHI[7]=phi7,PHI[9]=phi9,PHI[11]=phi11,PHI[13]=phi13},LS-RS),q,TT));
   XX:=normal(series( subs({Phi[1]=phi1,Phi[3]=phi3,Phi[5]=phi5,Phi[7]=phi7,Phi[9]=phi9,Phi[11]=phi11,Phi[13]=phi13},LS-RS),q,TT));
   RETURN(XX):
end:

CKCOF:=m->(-1)^(m+1)*(m+1)!*(m-1)!*(2*m)!;
CKCOF2:=m->(-1)^(m+1)*(m+1)!*(m-1)!;
flistgen:=proc(k) ## calculate f[j] in Corollary 4.5 [maincor]
     local m,LHS,RHS1,RHS2,RHS3,j,LHSa:
     global RHS:
     m:=1/2*(k-1):
     LHSa:=CKCOF(m)*CS^k*QINF^k:
     LHS:=(2*m)!*CS^k*QINF^k:
     RHS1:=NEWPOLY(_H[k],k,2*m):
     RHS2:=add(add(expand(binomial(2*m,a)*expand(symbDaFjmv2(a,j,m))*(j+1)^(2*m-a))*NEWPOLY(_H[k],k,2*m-a),
           a=0..2*m),j=1..m-1):
     RHS3:=-expand(symbDaF0mv2(2*m,m)):
     RHS:=collect(expand((RHS1+RHS2+RHS3)/CKCOF2(m)),[_H[k]]); 
     print(LHS=RHS);                      
     for j from 0 to m do
     print(_f[m-j],"=",coeff(RHS,_H[k],j));
     od:
     RETURN():
end: