### FILE: FUNCS.txt
#   Updated:   Mon Jun 13 20:36:03 EDT 2005
#   Previous updates:
#              Thu Jun  9 16:27:56 EDT 2005
#              Thu May 26 22:57:37 EDT 2005
#              Tue May 24 12:53:23 EDT 2005
#              Thu May 12 12:18:09 EDT 2005
#              Mon May  9 23:22:25 EDT 2005
#
# This file contains some maple programs for computing
# partitions, and q-series.
#

with(combinat):

 
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:

 
ptnOP:=proc(ptn)
local np,k;
   #Returns true if ptn is a partition into odd parts
   np:=nops(ptn); # np is number of parts of ptn
      for k from 1 to np do
	 if modp(ptn[k],2)=0 then
	    RETURN(false)
	 fi:
     od:
   RETURN(true):
end:

 
standptn:=proc(ptn)
local np,x,j,y;
   #prints a ptn in standard form
   np:=nops(ptn); # np is number of parts of ptn
   if np=1 then
     RETURN(ptn[1]):
   else
     x:=convert(ptn[np],string):
     for j from (np-1) by -1 to 1 do
        x:=cat(x," + ",convert(ptn[j],string));
     od:
     y:=convert(x,symbol):	
     RETURN(y);
   fi:
end:

 
printptns:=proc(ptns)
local nps,L,j;
   # Prints a list of ptns in standard form
   nps:=nops(ptns):
   L:=map(standptn,ptns):
   for j from nps by -1 to 1 do
     lprint(L[j]);
   od:
   RETURN():
end:
   
pmake:=proc(QS,q,T)
local b,X,m,g;

	#Computes formal product prod( (1-q^n)^b[n],n=1..T) 
        #of qseries QS O(q^(T+1)
	
	b[1]:=-coeff(QS,q,1):
        X:=(1-q)^b[1]:
	for m from 2 to T do
	    g := coeff(series(X,q,m+1),q,m):
            b[m]:= g - coeff(QS,q,m):
	    X:=X*(1-q^m)^b[m]:
	od:
	RETURN(X);
end:

 
P:=proc(n)
	#Computes p(n) the hard way using partition(n)
        RETURN(nops(combinat[partition](n)));
end:

 
PDP:=proc(n)
local ptns,ptns1;
	#Computes p(D,n) the hard way
        ptns:=combinat[partition](n):
	ptns1:=select(ptnDP,ptns):
	RETURN(nops(ptns1));
end:

 
 
ptnRR:=proc(ptn)
local np,k,d;
   #Returns true if ptn is a partition in which difference between parts
   #is at least 2.
   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<=1 then
	    RETURN(false)
	 fi:
     od:
   fi:
   RETURN(true):
end:

 
PRR:=proc(n)
local ptns,ptns1;
	#Computes the number of ptns of n in which difference between parts
        #is at least 2.

        ptns:=combinat[partition](n):
	ptns1:=select(ptnRR,ptns):
	RETURN(nops(ptns1));
end:

ptnCC:=proc(ptn)
local np,k,d;
   #Returns true if ptn is a partition in which difference between parts
   #is at least 3.
   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<=2 then
	    RETURN(false)
	 fi:
     od:
   fi:
   RETURN(true):
end:

ptnSCHURC:=proc(ptn)
local np,k,d;
   #Returns true if ptn is a partition in which difference between parts
   #is at least 3.
   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]:
	 pm:=modp(ptn[k+1],3):
	 if d<=2 then
	    RETURN(false)
	 else
	    if d=3 and pm<>0 then
	      RETURN(false):
	    fi:
	 fi:
     od:
   fi:
   RETURN(true):
end:

 
PSCHURC:=proc(n)
local ptns,ptns1;
	#Computes the number of ptns of n in which difference between parts
        #is at least 3.

        ptns:=combinat[partition](n):
	ptns1:=select(ptnSCHURC,ptns):
	RETURN(nops(ptns1));
end:
 
dissectq:=proc(f,q,m)
local S,r,dg,j;
	# For a polynomial f =sum a[n]*q^n
	# this proc returns the m-dissection of f
	# ie f = f[0] + f[1] + ... f[m-1]
        # f[j] = sum a[n]*q^n where n=j mod m
	S:=array(0..m-1):
	for r from 0 to m-1 do
	   S[r]:=0:
	od:
	dg:=degree(f,q):
	for j from 0 to dg do
	   r:=modp(j,m):
	   S[r]:=S[r]+coeff(f,q,j)*q^j:
	od:
	for r from 0 to m-1 do
	  S[r]:=factor(S[r]):
	od:
	RETURN(add(S[r],r=0..m-1));
end:


vecptns:=proc(n)
local j,PTNS,DPTNS,V,i,k,T,vp;
	#GENERATES A LIST OF VECTOR PARTITIONS OF n
	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:

 
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:

 
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:

vecptnsC:=proc(n,k,t)
       #GENERATES vector partitions of n with crank congruent to
       #k  mod 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;	

aqprod:=proc(a,q,n)
   #aqprod(a,q,n)=(a;q)[n]=(1-a)*(1-a*q)*...*(1-a*q^(n-1))
   local x,i,x1,x2:
   if type(n,integer) then
        x:=1:
      if n>0 then
         for i from 1 to n do
           x := x * (1-a*q^(i-1)):
         od:
      else
         m:=-n:
         x1:=mul((1-q/a*q^(i-1)),i=1..(m)):
	 x2:=(-q/a)^m*q^(m*(m-1)/2):
         x:=x2/x1:
      fi:	
   else
      x:=``(a,q)[n];
   fi:
   RETURN(x):
end:

qbin:=proc(q,m,n)
   #qbin(q,m,n) is the q-binomial coefficient or Gaussian polynomial
   # (q)[n]/((q)[m](q)[n-m]) for 0<=m<=n
        if whattype(m)=integer and whattype(n)=integer then
	  if m>=0 and m<=n then
	    RETURN(normal(aqprod(q,q,n)/aqprod(q,q,m)/aqprod(q,q,n-m)));
          else
            RETURN(0);
          fi:
          else
             ERROR(` m and n must be integers.`); 
        fi:
end:

qsplit:=proc(f::polynom)
    local tc,ld,td,f2;
    f2:=expand(f):##<--- this line added 10/19/99
                        ##If left unexpanded tcoeff may return incorrect value.
    tc := tcoeff(f2, q);
    ld := ldegree(f2, q);
    td := degree(f2, q);
    RETURN([f,tc,ld,td]);
end;

#######################################################################
# FUNCTION :  qfactor - factor a q polynom into factors (1-q^i)          
# CALLING SEQUENCE :  qfactor(f,T,s)       
#                      qfactor(f)           
#                      qfactor(f,T)         
#                      qfactor(f,s)         
# PARAMETERS :   f   - a polynomial in q
#                T - upper bound for i (optional)
#                s - string (optional)
# SYNOPSIS :   
# qfactor(f) attempts to factor the poly q into factors (1-q^i)
# if it fails it simply returns the input f.
# If the third (or second) argument s=test is given then qfactor
# will return FAIL if f can not be q-factored.
#######################################################################
##
 
 
 
qfactor:=proc()
	local m1, m2, LBITS, tc, ld, g, gl, gl1, gl2, gl3, gd, ggc, gprod, gser, gno,f;
	f:=args[1]:
	if nargs=1 or (nargs=2 and type(args[2],integer)<>true) then
	  m1:=4:
	  m2:=3:
	else
	  m1:=0:
	  m2:=args[2]
	fi:
	LBITS:=qsplit(f);
	tc:=LBITS[2]:
	ld:=LBITS[3]:
        g:=convert(f/tc/q^ld,polynom);
        gl:=convert(g,list);
	gl1:=subs(q=1,gl);
	gl2:=map(whattype,gl1);
	gl3:=convert(gl2,set);
	if gl3 = '{integer}' then
	   gd:=degree(g,q);
	   ggc:=coeff(g,q,0);
	   gprod:=pmake(g,q,m1*gd+m2);
	   gser:=convert(series(g-gprod,q,m1*gd+m2+5),polynom);
	   if gser<>0 then
	        if nargs=2 and args[2]=test then
		   RETURN(FAIL);
		fi:
		if nargs=3 and args[3]=test then
		   RETURN(FAIL);
		fi:
		if nargs=1 or (nargs=2 and args[2]<>test) then
		   RETURN(f);
	        fi:
	   else
	        gno:=normal(gprod-g);
	  	if gno = 0 then
		   RETURN(tc*gprod*q^ld);
		else
	            if nargs=2 and args[2]=test then
		       RETURN(FAIL);
		    fi:
		    if nargs=3 and args[3]=test then
		       RETURN(FAIL);
		    fi:
		    if nargs=1 or (nargs=2 and args[2]<>test) then
		       RETURN(f);
	            fi:
	         fi:
	    fi:
	else
	        if nargs=2 and args[2]=test then
		   RETURN(FAIL);
		fi:
		if nargs=3 and args[3]=test then
		   RETURN(FAIL);
		fi:
		if nargs=1 or (nargs=2 and args[2]<>test) then
		   RETURN(f);
	        fi:
	fi:
end;


beta2alpha:=proc(beta,a,q,n)
	# given sequence beta[n]
	# This proc returns alpha[n] such that (alpha[n],beta[n]) is
	# a Bailey pair
	option remember:
	local x:
	if n=0 then
	   beta(a,q,0):
	else
	   x:=add(beta2alpha(beta,a,q,r)/aqprod(q,q,n-r)/aqprod(a*q,q,n+r),r=0..n-1):
	   y:=aqprod(a*q,q,2*n)*(beta(a,q,n) - x):
	   RETURN(y):
	fi:
end: