#Code presented by Frank Garvan at his amazing guest-lecture
#April 25, 2013, Doron Zeilberger's Experimental Math
#Class

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

#freqtab(L): Inputs a list of numbers and outputs
#and rewrites it in frequency notation
#For example
#freqtab([2,2,2,4,4]);  outputs
#2,3
#4,2
freqtab:=proc(L)
local fset,i,M,j:
        fset:={}:
        for i from 1 to nops(L) do
           fset:=fset union {L[i]}:
           if type(M[L[i]],integer) then
              M[L[i]]:=M[L[i]]+1:
           else
              M[L[i]]:=1:
           fi:
        od:
        for j in fset do
        lprint(j,M[j]);
        od;
end:

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

#etaq(q,i,T): Inputs a symbol q, pos. integer i and
#a pos. T and outputs the Maclaurin polynomial of degree T of
#the infinite product
#(1-q^i)*(1-q^(2*i)*(1-q^(3*i)*.... .... 
etaq:=proc(q,i,T)
        local k,x,z1,z,w:
        z1:=(i + sqrt( i*i + 24*T*i ) )/(6*i):
        z:=1+trunc( evalf(z1) ):
        x:=0:
        for k from -z to z do
          w:=i*k*(3*k-1)/2:
          if w<=T then
           x:=x+ q^( w )*(-1)^k:
          fi:
        od:
        RETURN(x):
     end:

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

#prodmake(): inputs a polynomial f of q (arising, for out application
#from the begining of a q-series), a variable q, and a pos. integer T
#(and optionally, also, the name list)
#and outputs a product of the from 
#(1-q)^(a1)*(1-q^2)^(a2)*...*(1-q^T)^(aT)*...*
#that agrees with f up to the q^T power
#if list is given as the optional 4th argument, it returns the LIST
#[a1,a2, ..., aT]
	prodmake:=proc()
        local ft,f0,_a,_b,i,n,j,d,_A,_B,sum1,sum2,divj,divjb,m,prd,f,q,T,bseq,lst;
	####
	#USAGE: prodmake(f,q,T)       returns q-product
	#       prodmake(f,q,T,list)  returns q-product as a list of exponents
	#
	####
        #04/01/00: added nargs bit
        #03/30/07: add expand for symbolic prods
        if nargs>4 then
          ERROR(` number of arguments must be 3 or 4.`);
        fi:
        f:=args[1]:
        q:=args[2]:
        T:=args[3]:
	####
	ft:=series(f,q,T+5):
	if whattype(ft)=series  then
	   f0:=coeff(ft,q,0):
	   if f0=1 then
	      _b:=series(f,q,T):
	      _b:=convert(_b,polynom):
	      for i from 1 to T-1 do
	         _B[i]:=coeff(_b,q,i):
	      od:
	      _A[1]:=_B[1]:
 	      for n from 2 to T-1 do
                 sum2:=0:
	         for j from 1 to n-1 do
	            divj:=numtheory[divisors](j):
                    sum1:=0:
	            for d in divj do
	               sum1:=expand(sum1+d*(_A[d])):
                    od:
                    sum2:=expand(sum2 + (_B[n-j])*sum1):
                 od:
                 sum1:=0:
	            divjb:=numtheory[divisors](j) minus {n}:
                    for d in divjb do
	               sum1:=expand(sum1+d*(_A[d])):
	            od:
	         sum2:=expand(n*_B[n] - sum2 - sum1):
	         _A[n]:=sum2/n:
	      od:
	      ####
	      #04/01/00: added option of returning list
	      if nargs=3 then		
	         prd:=product((1-q^m)^(-_A[m]),m=1..(T-1)):
	         RETURN(prd):
	      else
	         bseq:=[seq(-_A[m],m=1..(T-1))]:
                 lst:=convert(bseq,list):
                 RETURN(lst):
	      fi:
	   else
	      ERROR(`coeff of q^0 must be 1`);
	   fi:
	else
	   ERROR(`f must be a series`);
	fi:
	end:
########################################################################


#etamake(): inputs a polynomial f of q (arising, for out application
#from the begining of a q-series), a variable q, and a pos. integer last
#(and optionally, also, the name list)
#and outputs a product of the from 
#q^b*eta(tau)^a1*eta(2*tau)^a2*...
#(tau is global variable such that q=exp(2*Pi*I*tau))
etamake:=proc(f,q,last) 
#
#VERSION: 01/30/10
# Use series(`leadterm`(expr),q) to get leadterm instead of tc etc.
# This is handle BBG c(q) function
#This procedure computes an eta-product expansion
#of the series f to order O(q^last).              
	   local fp,tc,exq,g,aa,ld,h,hh,i,cf,etaprod,alast,sf:
	   global ebase:
           sf:=series(f,q,last+10):
	   fp:=convert(sf,polynom):
	   tc:=tcoeff(fp,q):
	   exq:=ldegree(fp,q):
	   g:=normal(fp/tc/q^exq):
	   aa:=tc:
	   ld:=1:
	   ebase:=1:
	   alast:=last-exq:
	   while ld>0 do
	      h:=series(g-1,q=0,alast+1):
	      hh:=0:
	      for i from 1 to alast do 
	         hh:=hh+coeff(h,q,i)*q^i:
	      od:
	      h:=hh:
	      ##
	      ##08/20/99: bug fix ldegree(0,q) returns infinity
	      ##          in maple V release 5
	      if h=0 then
	        ld:=0:
	      else
	        ld:=ldegree(h,q):
	      fi:
	      cf:=coeff(h,q,ld):
	      if ld>0 then
	         exq:=exq+(ld*cf)/24:
	         aa:=eta(ld*tau)^(-cf)*aa:
	         g:=g*etaq(q,ld,alast)^cf:
                 ebase:=ilcm(ebase,ld):
	      fi:
	   od:
	   etaprod:=q^(exq)*aa:
	   RETURN(etaprod):
	end:
########################################################################

 #findcong(): inputs a polynomial QS, in q (for our applications
#coming from the beginning of a q-series), and a pos. integer T
#and outputs a list of conjenctured linear congruences
#If it has three arguments then the third (optional) argument, LM
#looks for congruences of the form
#f(a*n+b)==0 (mod c) with a<=LM (the default is sqrt(T)
#and an optional 4th argument, XSET, tells you a set of c's
#that you want to exclude (being chaff rather than wheat)
 findcong:=proc()
 local s,j,c,S,p,r,QS2,S1,QS,T,LM,LM1,M,xxT,xxT1,CFS1,IGCD1,CFS,IGCD,f,IFX,HH,h,P,R,m,L,ra,Ma,pa,XSET;
 global xprint:
        ## findcong(QS,T)
        ## findcong(QS,T,LM)
        ## findcong(QS,T,LM,XSET)
        ## find linear congruences for q-series QS up to q^T
        ## EXAMPLE:
        ## P:=series(1/etaq(q,1,1001),q,1001):
        ## > findcong(P,1000);
        ##                        [4, 5, 5]
        ##                        [5, 7, 7]
        ##                       [6, 11, 11]
        ##                       [24, 25, 25]
        ##    {[5, 7, 7], [4, 5, 5], [24, 25, 25], [6, 11, 11]}
        ## NOTE: XSET set of moduli to exclude
        if nargs>4 or nargs<2 then
           ERROR(`findcong takes 2, 3 or 4 args.`);
        fi:
        if not(type(xprint,boolean)) then
          xprint:=false:
        fi:
        QS:=args[1]:
        T:=args[2]:
        LM1:=2:
        XSET:={}:
        if nargs=2 then
           LM:=trunc(sqrt(T)):
        else
           LM:=args[3]: 
        fi:
        if nargs=4 then
           XSET:=args[4]:
        fi:
        S:={}:
      for M from LM1 to LM  do
      for r from 0 to M-1 do
         xxT:=trunc((T-r)/M):
         xxT1:=min(xxT,100):
         CFS1:=[seq(coeff(QS,q,M*n+r),n=0..xxT1)]:
         if xprint then
            lprint("xxT1=",xxT1,"CFS1=",CFS1);
         fi:
         IGCD1:=igcd(op(CFS1));
         if xprint then
            lprint("M=",M,"r=",r,"IGCD1=",IGCD1);
         fi:
         if IGCD1<>1 and not member(IGCD1,XSET) then
         CFS:=[seq(coeff(QS,q,M*n+r),n=0..trunc((T-r)/M))]:
         IGCD:=igcd(op(CFS));
         if IGCD<>1 then
           # new one?
           f:=1:
           IFX:=ifactors(IGCD):
           HH:=IFX[2]:
           for h from 1 to nops(HH) do
              P:=HH[h][1]:
              R:=HH[h][2]:
           for m from 1 to nops(S) do
              L:=S[m]:
              ra:=L[1]: Ma:=L[2]: pa:=L[3]:
              if modp(r,Ma)=ra and pa=P^R and modp(M,Ma)=0 then
                 f:=0:
              fi:
           od:
           if f=1 then
              S1:= S union {[r,M,P^R]}: 
              if nops(S1) > nops(S) then 
                print([r,M,P^R]); 
                S:=S1:
              fi:
           fi:
           od:
          fi:
          fi:
       od:od:
 RETURN(S);
 end:           


########################################################################
#sift(s,q,n,k,T), inputs  a polynomial s in q (of degree T),
#and pos. integers n,k,T
#and outputs sum a_[ni+k] q^i
#if s= sum a_i q^i
	sift:=proc(s,q,n,k,T)
	        #
		# sum a_i q^i --> sum a_[ni+k] q^i
	        #
	 	local y,i,st,lasti:
                st:=series(s,q,T+5):
        if whattype(st)=series  then
	 	y:=0:
		lasti:=floor((T-k)/n):
	 	for i from 0 to lasti do
	 	  ##y:=y+coeff(s,q,(n*i+k))*q^i:
                  ##fixed 11-16-06
	 	  y:=y+coeff(st,q,(n*i+k))*q^i:
	 	od:
		RETURN(y):
	else
	      ERROR(`s must be a series`);
	fi:
	end:

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

with(combinat):

#ptnDP(ptn): inputs a partition, ptn, and returns true (false)
#if all its parts are distinct (not so)
ptnDP:=proc(ptn)
local np,k,d;
   #Returns true if ptn is a partition into distinct parts
   np:=nops(ptn); # np is number of parts of ptn
   if np=1 then
        RETURN(true):
   else
      for k from 1 to np-1 do
         d:=ptn[k+1]-ptn[k]:
         if d=0 then
            RETURN(false)
         fi:
     od:
   fi:
   RETURN(true):
end:
vpw:=proc(vptn)
local P1,np;
    # weight of vector ptn [P1,P2,P3] = (-1)^(#(P1))
    P1:=vptn[1]:
    np:=nops(P1):
    RETURN( (-1)^(np) ):
end:
 
#vecptns(n): inputs a pos. integer n and outputs the set of triples
#[DPTNS[i],PTNS[j],PTNS[k]], 
#obeying some conditions (see the code below)
vecptns:=proc(n)
local j,PTNS,DPTNS,V,i,k,T,vp,vp1,vp2,vp3,s1,s2,s3;
        #GENERATES A LIST OF VECTOR PARTITIONS OF n (cor. to spt)
        for j from 0 to n do
           PTNS[j]:=partition(j):
           DPTNS[j]:=select(ptnDP,PTNS[j]):
        od:
        V:=[]:
        for i from 0 to n do
           for j from 0 to n-i do
                k:=n-i-j:
                T:=cartprod([DPTNS[i],PTNS[j],PTNS[k]]):
                while not T[finished] do 
                  vp:=T[nextvalue]():
                  V:=[op(V),vp]: 
                od:
           od:
        od:
        RETURN(V):
end:

#lp(ptn): largest part of a partion ptn
lp:=proc(ptn)
  if nops(ptn)=0 then
     0:
  else
    ptn[nops(ptn)]:
 fi:
end:

#vpcrank(vptn): The famous Garvan original crank 
#of a vector partition
vpcrank:=proc(vptn)
local P2,np2,P3,np3;
    # crank of vector ptn [P1,P2,P3] is #(P2)-#(P3)
    P2:=vptn[2]:
    np2:=nops(P2):	
    P3:=vptn[3]:
    np3:=nops(P3):	
    RETURN(np2 - np3):
end:


#inputs pos. integers n,k,t
#GENERATES new vector partitions of n with crank congruent to
#k  mod t
vecptnsC:=proc(n,k,t)
       local localCS,vptns:
       localCS:=proc(vptn)
          if modp(vpcrank(vptn),t)=k then
             true:
          else:
             false:
          fi:
       end:
       vptns:=vecptns(n):
       select(localCS,vptns);
end:

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

#nu2(F): the highest power of 2 divisible by all the coeff. of F
nu2:=proc(F)
local n,F1:
n:=0:
F1:=0:
while F1=0 do
  if modp(F,2^(n+1))=0 then
    n:=n+1: 
  else
    F1:=1:
   fi:
od:
RETURN(n):
end:

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