#######################################################
## make-win-package qseries
## Sun, Feb  9, 2020  6:37:42 PM
#######################################################
#######################################################

printf("BEGIN qseries package\n");
printf("THIS VERSION DATED Sun, Feb  9, 2020  6:37:42 PM \n");
qseries:=table();
qseries[aqprod]:=proc(a,q,n)
   #11/25/99: added else bit when n not an integer
   #          (existed in an earlier version but got lost somehow)
   #09.07.14: can handle negative n
   local x,i,n1,m1,x1:
   if type(n,integer) then
      if n >=0 then
         x:=1:
         for i from 1 to n do
            x := x * (1-a*q^(i-1)):
         od:
      else
         n1:=-n:
         m1:=mul(1-q^(j+1)/a,j=0..n1-1):
         x1:=(-q/a)^n1*q^(n1*(n1-1)/2):
         RETURN(x1/m1):
      fi:
   else
      x:=``(a,q)[n];
   fi:
   RETURN(x):
end:
qseries[changes]:=proc()
printf("**************************************************************\n");
printf("*\n");
printf("*\n");
printf("*  qseries package version 1.1 -  Mon Jul 16 15:24:26 EDT 2012\n");
printf("*  qseries package version 1.2  - Thu Dec 20 11:52:01 EST 2012\n");
printf("*  qseries package version 1.3  - Fri Aug 12 15:07:08 EDT 2016\n");
printf("*  qseries package version 1.3a - Tue Nov 28 13:48:19 EST 2017\n");
printf("*  qseries package version 1.3b - Wed, Mar  6, 2019  8:15:20 PM\n");
printf("*  qseries package version 1.3c - Wed, Jul 17, 2019  9:41:24 AM\n");
printf("*  qseries package version 1.3d - Sat, Oct 19, 2019  9:58:53 AM\n");
printf("*  qseries package version 1.3e - Fri, Jan 24, 2020 11:06:28 PM\n");
printf("*  qseries package version 1.3f - Sat, Feb  8, 2020 10:33:00 PM\n");
printf("*\n");
printf("*  This version tested on MAPLE 2019\n");
printf("*\n");
printf("*  Changes since version 1.3\n");
printf("*\n");
printf("*    *   New function J2jaclist \n");
printf("*         Takes 0 or 1 args. Converts theta-quotient in J-notation \n");
printf("*         to a jaclist. J[b,a]^c and J[b]^c  correspond to \n");
printf("*         [b,a,c] or [b,c] in the list.  \n");
printf("*\n");
printf("*    *   New function jaclist2JACPROD \n");
printf("*         Takes 0 or 2 args.Converts theta-quotient as a jaclist to\n");
printf("*         JAC-notation. \n");
printf("*\n");
printf("*    *   New function dilatejaclist \n");
printf("*         Takes 0 or 2 args.\n");
printf("          Does q -> q^r in jacprod corresponding to\n");
printf("*         a jaclist. It is used by Jterm2JACPRPOD. \n");
printf("*\n");
printf("*    *   New function Jterm2JACPROD \n");
printf("*           Takes 0 or 2 args.Converts jterm to JACPROD and also  \n");
printf("*           does q -> q^r using dilatejaclist.                \n");
printf("*    *   New function jac2J\n");                         
printf("*         This functions takes 0 or 3 args. \n");
printf("*         This is just a variant of etamake (or qetamake) \n");
printf("*         which returns J[m] instead of eta(m*tau) etc.\n");
printf("*\n");
printf("*      * New function Jetamake\n");                         
printf("*         This functions takes 0 or 3 args. \n");
printf("*         This is just a variant of etamake (or qetamake) \n");
printf("*         which returns J[m] instead of eta(m*tau) etc.\n");
printf("*\n");
printf("*      *  Testing new function version findcongbeta\n");
printf("*         Now takes up to 5 args try findcongbeta() for help\n");
printf("*\n");
printf("*      *  Removed print from findprod \n");
printf("*\n");
printf("*      *  Fixed small bug in findmaxind \n");
printf("*\n");
printf("*      *  Fixed bug in findhom (coeffs now computed from constant\n");
printf("*         term instead of from q^1) \n");
printf("*\n");
printf("*      *  Added global var _xprint to findhomcombo\n");
printf("*\n");
printf("*\n");
printf("*  Changes since version 1.2\n");
printf("*\n");
printf("*      *  Added global vars _etalisttop, _etalistbot, \n");
printf("*         _etaconst to etamake (used by latex progs in \n");
printf("*         new ramarobinsids package)\n");
printf("*\n");
printf("*      *  New function findmaxind \n");
printf("*\n");
printf("*      *  Changed taylor to series in both findhom and findhommodp\n");
printf("*         so they can handle negative powers of q.                \n");
printf("*\n");
printf("*      *  New version of sift can now handle negative exponents\n");
printf("*         See oldsift for old version. The new sift needs the\n");
printf("*         functions lqdegree and lqdegree0 from thetaids package.\n");
printf("*         So these have also been added to qseries package.\n");
printf("*\n");
printf("*      *  Added xprint option to findhom to avoid extra printing in large loops\n");
printf("*         Default is xprint=false (no extra printing)                          \n");
printf("*\n");
printf("*      *  jac2series can now handle JAC(a,b,infinity) when \n");
printf("*         a<0 or a>b.                                      \n");
printf("*\n");
printf("*      *  fixed major bug in findcong\n");
printf("*         The old version missed some congruences.\n");
printf("*\n");
printf("*      *  findhom works correctly when nops(L)=1\n");
printf("*\n");
printf("*      *  added new function checkmult\n");
printf("*         This function checks with a give qseries has multiplicative coefficients.\n");
printf("          USAGE:                                     \n");
printf("             checkmult()                                         \n");
printf("             checkmult(QS,T)                                         \n");
printf("             checkmult(QS,T,yes)                                     \n");
printf("*\n");
printf("*      *  aqprod(a,q,n) can now handle negative n\n");
printf("*\n");
printf("*      *  added new function checkprod\n");
printf("*         This function checks a qseries has a nice product.\n");
printf("          USAGE:                                     \n");
printf("             checkprod(f,M,Q)                                    \n");
printf("*\n");
printf("*      *  added new function findprod\n");
printf("*         This function searches for Z-linear combinations of \n");
printf("*         the functions in FL which are probable products. \n");
printf("*         This function also uses checkprod.               \n");
printf("          USAGE:                                     \n");
printf("             findprod(FL,T,M,Q)                                 \n");
printf("*\n");
printf("*      *  aqprod(a,q,n) can now handle negative n\n");
printf("*\n");
printf("*  Changes since version 1.1\n");
printf("*\n");
printf("*      *  Added global var jacp to jacmake.\n");
printf("*\n");
printf("*      *  In jacprodmake change failure from ERROR to WARNING:\n");
printf("          In jacprodmake changed handling of half-integer exponents.\n");
printf("*\n");
printf("*      *  In etamake added  var ebase (to be used by jacprodmake\n");
printf("*          when converting an etaprod to a jacprod.\n");
printf("*\n");
printf("*\n");
printf("*  Changes since version 1.0\n");
printf("*\n");
printf("*      *  Added new function qdegree\n");
printf("*         qdegree(QS) returns degree in q.\n");
printf("*\n");
printf("*      *  Added new function findcong\n");
printf("          This function finds congruences c(A*n+B) == 0 mod C\n");
printf("          for a give qseries f=add(c(n)*q^n).\n");
printf("*\n");
printf("*      *  Added new function qs2jaccombo\n");
printf("          This function converts a sum of q-series to a sum\n");
printf("          of jacprods if possible. \n");
printf("*\n");
printf("*      *  Local vars added to procs where necessary\n");
printf("*\n");
printf("*      *  jacprodmake now takes 3 or 4 args\n");
printf("*         jacprodmake(f,q,T) works as before.\n");
printf("*         jacprodmake(f,q,T,P) returns a jacprod (if it exists) on\n");
printf("          base q^b where b is a divisor of P. ERROR is returned \n");
printf("          if f is not a jacprod.\n");
printf("          This new version also handles the case where the coeff \n");
printf("          of q^0 is not 1.\n");
printf("          PLEASE report any bugs.\n");
printf("*\n");
printf("*      *  New hidden function periodfind2 added. \n");
printf("          This is used by jacprodmake to find period assuming it\n");
printf("          is a divisor of a given integer.\n");
printf("*\n");
printf("*      *  New hidden function jacmake now takes 2 or 3 args.\n");
printf("*         It is used by jacprodmake.\n");
printf("*\n");
printf("*      *  Added instructions for windows installation\n");
printf("*         Tested with Maple 11\n");
printf("*\n");
printf("*      *  Corrected email and homepage listed below\n");
printf("*         New email: fgarvan@ufl.edu\n");
printf("*         New homepage:  http://www.math.ufl.edu/~fgarvan\n");
printf("*\n");
printf("*  Changes since version 0.9\n");
printf("*\n");
printf("*      *  Updated packageversion with new email address and homepage.\n");
printf("*\n");
printf("*      *  qbin accepts q=root of unity etc (Suggested by Michael Somos)\n");
printf("*\n");
printf("*      *  Added function findlincombomodp\n");
printf("*\n");
printf("*      *  Fixed bug in findlincombo\n");
printf("*\n");
printf("*      *  Local username changed from frank to fgarvan\n");
printf("*         New email: fgarvan@ufl.edu\n");
printf("*         New homepage:  http://www.math.ufl.edu/~fgarvan\n");
printf("*\n");
printf("*      *  Corrected typo in online help for findlincombo\n");
printf("*\n");
printf("*  Changes since version 0.8\n");
printf("*\n");
printf("*      *  Added online FUNCTIONS page:\n");
printf("*         http://www.math.ufl.edu/~fgarvan/qmaple/functions/\n");
printf("*\n");
printf("*      *  Added gooo to global list in findhomcombo\n");
printf("*      *  Added new function findlincombo\n");
printf("*\n");
printf("*      *  Fixed bugs in findhomcombo and modp version\n");
printf("*         There was unnecessary etamake comp when etaopt=no\n");
printf("*\n");
printf("*      *  updated prodmake so it can have symbolic exponents\n");
printf("*         Used series instead of taylor\n");
printf("*\n");
printf("*      *  new function findhomcombomodp\n");
printf("*         This is a modp version of  findhomcombo\n");
printf("*         Useage: findhomcombomodp(f,L,p,q,n,topshift,etaoption))\n");
printf("*\n");
printf("*      *  Fixed a bug in findhomcombo when nops(L)=1\n");
printf("*\n");
printf("*      *  Fixed bug in findpoly \n");
printf("*\n");
printf("*      *  In qetamake change E(q) to _E(q) \n");
printf("*\n");
printf("*      *  Fixed minor bug in sift \n");
printf("*\n");
printf("*  Changes since version 0.7\n");
printf("*\n");
printf("*      *  changed name of function version to packageversion\n");
printf("*         since version is maple builtin function\n");
printf("*         Useage: packageversion()\n");
printf("*\n");
printf("*      *  new function findhommodp\n");
printf("*         This is a modp version of  findhom\n");
printf("*         Useage: findhommodp(L,p,q,n,topshift)\n");
printf("*\n");
printf("*      *  new function mprodmake\n");
printf("*         Finds product expansion (1+q^n1)*(1+q^n2)*...\n");
printf("*         Useage: mprodmake(f,q,last)\n");
printf("*\n");
printf("*      *  fixed bug in qetamake         \n");
printf("*\n");
printf("*  Changes since version 0.6 (Nov 2002):\n");
printf("*\n");
printf("*      *  new function qetamake          \n");
printf("*         This is just a variant of etamake which returns E(q^m)\n");
printf("*         instead of eta(m*tau).\n");
printf("*\n");
printf("*      *  fixed bug in jac2series        \n");
printf("*\n");
printf("*  Changes since version 0.5:\n");
printf("*\n");
printf("*      *  new function zqfactor          \n");
printf("*         Takes 4 or 5 args: \n");
printf("*         zqfactor(F,z,q,N,buglim) \n");
printf("*         zqfactor(F,z,q,N) \n");
printf("*            factors F into a product of terms \n");
printf("*            of the form (1 - a[i,j] z^i q^j).\n");
printf("*            Factor F into a product of terms (1 - a[i,j] z^i q^j)\n");
printf("*            N = largest j \n");
printf("*            buglim = max number of terms (default is 1000).\n");
printf("*\n");
printf("*  Changes since version 0.4:\n");
printf("*\n");
printf("*      *  prodmake now takes 3 or 4 args.\n");
printf("*         prodmake(f,q,t) returns a q-product.\n");
printf("*         prodmake(f,q,t,list) returns a q-product as a list of exponents.\n");
printf("*      *  new function changes() - returns changes since previous release.\n");
printf("*\n");
printf("*  Please report any problems to fgarvan@ufl.edu\n");
printf("*  See \n");
printf("*  http://www.math.ufl.edu/~fgarvan/qmaple/qmaple.html \n");
printf("*  for documentation and help.\n");
printf("*\n");
printf("**************************************************************\n");
RETURN():
end:
qseries[checkmult]:=proc()
 local NICE,m,n,c1,c2,c3,QS,T:
 NICE:=1:
 if not member(nargs,{0,2,3}) then
     ERROR("checkmult takes 0, 2, or 3 args");
 fi:
 if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf(" checkmult(QS,T)                                         \n");
     printf(" checkmult(QS,T,yes)                                     \n");
     printf("    Returns 1 if QS up to Q^T is multiplicative           \n");
     printf("    If nargs=3 and arg3=yes then m,n where QS fails to    \n");
     printf("    be multiplicative are printed.                        \n");
     printf("    Otherwise only the first instance of failure is       \n");
     printf("    printed.                                              \n");
     printf("---------------------------------------------------------\n");
else
 QS:=args[1]:
 T:=args[2]:
 if nargs=2 then
    for n from 2 to trunc(T/2) do
    for m from 2 to trunc(T/2) do
    if igcd(m,n)=1 and m*n<=T then
       c1:=coeff(QS,q,m):
       c2:=coeff(QS,q,n):
       c3:=coeff(QS,q,m*n):
       if c3<>c1*c2 then
          lprint("NOT MULTIPLICATIVE");
          lprint("failure at m=",m,"and n=",n);
          RETURN(0):
       fi:
    fi:
    od:
    od:
    lprint(" MULTIPLICATIVE");
    RETURN(1):
 else
    for n from 2 to trunc(T/2) do
    for m from 2 to trunc(T/2) do
    if igcd(m,n)=1 and m*n<=T then
       c1:=coeff(QS,q,m):
       c2:=coeff(QS,q,n):
       c3:=coeff(QS,q,m*n):
       if c3<>c1*c2 then
          lprint("NOT MULTIPLICATIVE");
          lprint("failure at m=",m,"and n=",n);
          NICE:=0:
       fi:
    fi:
    od:
    od:
    if NICE=1 then
    lprint(" MULTIPLICATIVE");
    fi:
    RETURN(NICE):
 fi:
fi:
end:
qseries[checkprod]:=proc(f,M,Q)
local ss,ssp,LL,nss,c0,cb0,nss1,xx,mx:
    ss:=series(f,q,Q+1):
    ssp:=convert(ss,polynom):
    LL:=ldegree(ssp,q):
    nss:=expand(ssp/q^LL):
    c0:=coeff(nss,q,0):
    cb0:=abs(c0):
    ##if coeff(nss,q,0)=1 then
if ssp<>0 then
    nss1:=expand(nss/c0):
    if modp(ssp,cb0)=0 then
       xx:=qseries[prodmake](nss1,q,Q-LL,list):
       mx:=max(op(map(abs,xx))):
       if mx<M then
        RETURN([LL,1]):
       else
        RETURN([LL,mx]):
       fi:
    else
        RETURN([[LL,coeff(nss,q,0)],-1]):
    fi:
else
   RETURN([ZEROTERM,oo]):
fi:
end:
qseries[dilatejaclist]:=proc()
   local njaclist,L,jaclist,r:
   # 02.08.2020
#############################################################################
##  dilatejaclist(jaclist,r) does q->q^r in the jacprod corresponding to
##  jaclist.
## > jlist:=[[25, 5, -1], [25, 10, 1], [25, 1], [50, 5, 1]];
## > jaclist2JACPROD(jlist);
##         JAC(10, 25, infinity) JAC(0, 25, infinity) JAC(5, 50, infinity)
##         ---------------------------------------------------------------
##                              JAC(5, 25, infinity)
## > jlist2:=dilatejaclist(jlist,1/5);
##              jlist2 := [[5, 1, -1], [5, 2, 1], [5, 1], [10, 1, 1]]
## > jaclist2JACPROD(jlist2);
##           JAC(2, 5, infinity) JAC(0, 5, infinity) JAC(1, 10, infinity)
##           ------------------------------------------------------------
##                               JAC(1, 5, infinity)
#############################################################################
 if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf(" dilatejaclist(jaclist,r)                                \n");
     printf("    Does q -> q^r in the jacprod corresponding to jaclist.\n");
     printf("                                                          \n");
     printf(" EXAMPLE:                                                 \n");
     printf(" jlist:=[[25, 5, -1], [25, 10, 1], [25, 1], [50, 5, 1]];\n");
     printf(" jlist2:=dilatejaclist(jlist,1/5); \n");
     printf("--------------------------------------------------------\n");
else
   if nargs=2 then
      njaclist:=[]:
      jaclist:=args[1]:
      r:=args[2]:
      for L in jaclist do
         if nops(L)=2 then
            njaclist:=[op(njaclist),[L[1]*r,L[2]]]:
         else
    
           njaclist:=[op(njaclist),[L[1]*r,r*L[2],L[3]]]:
         fi:
     od:
   RETURN(njaclist):
  else
    ERROR("nargs = 0 or 2");
  fi:
fi:
end:

qseries[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).              
#12/20/12: added global var ebase (to be used by jacprodmake when
#          converting an etaprod to a jacprod
	   local fp,tc,exq,g,aa,ld,h,hh,i,cf,etaprod,alast,sf:
           global ebase,_etalisttop,_etalistbot,_etaconst:
#07/11/16: added global vars _etalisttop, _etalistbot
           ebase:=1:
           _etalisttop:=[]:
           _etalistbot:=[]:
           sf:=series(f,q,last+10):
	   fp:=convert(sf,polynom):
	   tc:=tcoeff(fp,q):
           _etaconst:=tc:
	   exq:=ldegree(fp,q):
	   g:=normal(fp/tc/q^exq):
	   aa:=tc:
	   ld:=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(ld,ebase):
                 if cf<0 then
                   _etalisttop:=[op(_etalisttop),[ld,-cf]]
                 else
                   _etalistbot:=[op(_etalistbot),[ld,cf]]
                 fi:
	      fi:
	   od:
	   etaprod:=q^(exq)*aa:
	   RETURN(etaprod):
	end:
     qseries[etaq]:=proc(q,i,trunk)
#       This proc returns (q^i:q^i)_inf up to q^trunk
        local k,x,z1,z,w:
        z1:=(i + sqrt( i*i + 24*trunk*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<=trunk then
           x:=x+ q^( w )*(-1)^k:
          fi:
        od:
        RETURN(x):
     end:

 
 qseries[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,fH;
 global xprint:
 ## VERSION: 09/09/2013
 ## PREVIOUS VERSION: 07/25/2012
        ## 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:={}:
####
##"xxT1=", 11, "CFS1=", [660, 91740, 3593480, 76836100, 1128188020, 
## 12733555780, 117951349040, 935191373080, 6532841375940, 41057165666720, 
## 235820387905900, 1253017582965800]
##"M=", 40, "r=", 27, "IGCD1=", 20
####
      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:
         #next line fixed 09/09/2013
         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 and not member(IGCD,XSET) then
           # new one?
           f:=0:
           IFX:=ifactors(IGCD):
           HH:=IFX[2]:
           if xprint then
              lprint("IGCD=",IGCD,"HH=",HH);
           fi:
           for h from 1 to nops(HH) do
              P:=HH[h][1]:
              R:=HH[h][2]:
              #BUG FIX 10.20.214 (In the old version the variable fH was 
              #        called f and did not appear here but was outside
              #        this loop and hence some congruences were missed)
              fH:=1:
              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
                     fH:=0: 
                  fi:
              od:
              if fH=1 then
              # next line added 7/27/13
                 if not member(P^R,XSET) then
                    S1:= S union {[r,M,P^R]}: 
                    if nops(S1) > nops(S) then 
                       print([r,M,P^R]); 
                       S:=S1:
                    fi:
                  fi:
               fi:
           od:
          fi:
          fi:
       od:od:
 RETURN(S);
 end:           
#++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 qseries[findcongbeta]:=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,fH,xxTLIM;
 global xprint:
 ## VERSION: 09/09/2013
 ## PREVIOUS VERSION: 07/25/2012
        ## findcongbeta(QS,T)
        ## findcongbeta(QS,T,LM)
        ## findcongbeta(QS,T,LM,XSET)
        ## findcongbeta(QS,T,LM,XSET,xxTLIM)
        ## find linear congruences for q-series QS up to q^T
        ## If QS = sum a(n) q^n then findcong searches for congruences
        ##    a(M*n + r) == 0 (mod K)
        ## and returns possible set of congruences : {[r1,M1,K1], ..... }
        ## EXAMPLE:
        ## P:=series(1/etaq(q,1,1001),q,1001):
        ## > findcongbeta(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]}
        ## NOTE1: LM is maximum M (default = [sqrt(T)]
        ## NOTE2: XSET set of moduli to exclude
        ## NOTE3: xxTLIM max number of coeffs to test gcd (default = 100)
        if nargs>5 or nargs=1 then
           ERROR(`findcong takes 0, 2, 3 or 4  or 5 args.`);
        fi:
 if nargs=0 then
   printf("  findcongbeta(QS,T)\n");
   printf("  findcongbeta(QS,T,LM)\n");
   printf("  findcongbeta(QS,T,LM,XSET)\n");
   printf("  findcongbeta(QS,T,LM,XSET,xxTLIM)\n");
   printf("  find linear congruences for q-series QS up to q^T\n");
   printf("  If QS = sum a(n) q^n then findcong searches for congruences\n");
   printf("     a(M*n + r) == 0 (mod K)\n");
   printf("  and returns possible set of congruences : {[r1,M1,K1], ..... }\n");
   printf("  EXAMPLE:\n");
   printf("  P:=series(1/etaq(q,1,1001),q,1001):\n");
   printf("  > findcongbeta(P,1000);\n");
   printf("                         [4, 5, 5]\n");
   printf("                         [5, 7, 7]\n");
   printf("                        [6, 11, 11]\n");
   printf("                        [24, 25, 25]\n");
   printf("     {[5, 7, 7], [4, 5, 5], [24, 25, 25], [6, 11, 11]}\n");
   printf("  NOTE1: LM is maximum M (default = [sqrt(T)]\n");
   printf("  NOTE2: XSET set of moduli to exclude\n");
   printf("  NOTE3: xxTLIM max number of coeffs to test gcd (default = 100)\n");
else
        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:
        if nargs=5 then
          xxTLIM:=args[5]:
        else
          xxTLIM:=100:
        fi:
        S:={}:
####
##"xxT1=", 11, "CFS1=", [660, 91740, 3593480, 76836100, 1128188020, 
## 12733555780, 117951349040, 935191373080, 6532841375940, 41057165666720, 
## 235820387905900, 1253017582965800]
##"M=", 40, "r=", 27, "IGCD1=", 20
####
      for M from LM1 to LM  do
      for r from 0 to M-1 do
         xxT:=trunc((T-r)/M):
         xxT1:=min(xxT,xxTLIM):
         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:
         #next line fixed 09/09/2013
         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 and not member(IGCD,XSET) then
           # new one?
           f:=0:
           IFX:=ifactors(IGCD):
           HH:=IFX[2]:
           if xprint then
              lprint("IGCD=",IGCD,"HH=",HH);
           fi:
           for h from 1 to nops(HH) do
              P:=HH[h][1]:
              R:=HH[h][2]:
              #BUG FIX 10.20.214 (In the old version the variable fH was 
              #        called f and did not appear here but was outside
              #        this loop and hence some congruences were missed)
              fH:=1:
              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
                     fH:=0: 
                  fi:
              od:
              if fH=1 then
              # next line added 7/27/13
                 if not member(P^R,XSET) then
                    S1:= S union {[r,M,P^R]}: 
                    if nops(S1) > nops(S) then 
                       print([r,M,P^R]); 
                       S:=S1:
                    fi:
                  fi:
               fi:
           od:
          fi:
          fi:
       od:od:
 RETURN(S);
fi:
 end:           
#++++++++++++++++++++++++++++++++++++++++++++++++++++

qseries[findhom]:=proc(L,q,n,topshift)
#This program generates a list of potential homogeneous relations
#of degree n among the list of q-series L.
#HISORY:
#  11.28.17: Coefficients now computed from constant term
#            findhom([1,q],1,0) gave  {X[1]} in error
#
#  05/07/13: now handles case nops(L)=1 correctly.
#  04.18.13: added global xprint to annoying extra prints
#  10.28.15: changed taylor to series
   local xx,x,y,z,nx,n2,n3,gg,goo,gooo,c,dh,jj,g,k,w,d,i,ct,D,ii,h,G,ANS,yy;
   global X,xprint:
   if not(type(xprint,boolean)) then
          xprint:=false:
   fi:
if whattype(n)=integer and whattype(topshift)=integer then
if n>0 and topshift>-20 then

   D:=nops(L);
   gg[1]:={};
   xx:=expand(sum(x[i],i=1..D)^(n) ):
   nx:=nops(xx):
   n2:=nx+20+topshift:#step this up if some solutions are spurious
   n3:=n2+1:
   ################################################################
   ##05/07/13: handle case when nops(L)=1
   if nops(L)=1 then
      yy:=convert(series(L[1],q,n3),polynom):
      if yy=0 then
         RETURN({1}):
      else
         RETURN({{}}):
      fi:
   fi:
   ################################################################
   if xprint then
      print(`# of terms `,n3);
   fi:
   c:=array(1..nx,1..n2);
   dh:=array(1..nx,1..D);
   for jj from 1 to nx do
      G:=op(jj,xx):
      for i from 1 to D do
      h[i]:=degree(G,x[i]):
      dh[jj,i]:=h[i]:
      od:
      w:=product(L[ii]^h[ii],ii=1..D):
      d:=series(w,q,n3+5);
      for i from 1 to n2 do
         c[jj,i]:=coeff(d,q,i-1);
         ##11.28.2017: old c[jj,i]:=coeff(d,q,i);
      od;
   od:
   ct:=  linalg[transpose](c):
   gooo:=(linalg[kernel](ct));
   for k from 1 to nops(gooo) do
      goo:=gooo[k];
      d:=op(dh):
      g[k]:=0:
      for i from 1 to nx do
         g[k]:=g[k]+goo[i]*product(X[ii]^d[i,ii],ii=1..D):
      od;
      gg[k]:=g[k];
  od;
  if xprint then
     print(`-----RELATIONS----of order---`,n);  
  fi:
   ANS:=convert(gg,set);
   RETURN(ANS);
else
   ERROR(`n must be a positive integer and topshift must apositive integer greater then -20`);
fi:
else
   ERROR(`n and topshift must be integers.`);
fi:
end:

 
qseries[findhomcombo]:=proc(f,L,q,n,topshift,etaoption)
#This program tries to express f as a homogeneous polynomial     
#of degree n among the list of q-series in L.
#
#01/26/06: Fixed bug when nops(L)=1
#
   local xx,x,y,z,nx,n2,n3,gg,goo,c,dh,jj,g,k,w,d,i,D,ii,h,G,ANS,func,fq,ANSETA,geta;
   global X,ct,gooo,xprint:
   #06/20/07: added gooo to global
   #03/03/17: added xprint to global
   if not(type(xprint,boolean)) then
          xprint:=false:
   fi:
   X:='X':
if whattype(n)=integer and whattype(topshift)=integer then
if n>0 and topshift>-20 then

   D:=nops(L);
   gg[1]:={};
   xx:=expand(sum(x[i],i=1..D)^(n) ):
   nx:=nops(xx):
   if D=1 then
      nx:=1:
   fi:
   n2:=nx+20+topshift:#step this up if some solutions are spurious
   n3:=n2+1:
   if xprint then
    print(`# of terms `,n3);
   fi:
   c:=array(1..nx+1,1..n2);
   dh:=array(1..nx,1..D);
   if D>1 then
      for jj from 1 to nx do
         G:=op(jj,xx):
         for i from 1 to D do
         h[i]:=degree(G,x[i]):
         dh[jj,i]:=h[i]:
         od:
         w:=product(L[ii]^h[ii],ii=1..D):
         d:=series(w,q,n3);
         func[jj]:=d:
         for i from 1 to n2 do
            c[jj,i]:=coeff(d,q,i);
         od;
      od:
   else
      d:=series(L[1]^n,q,n3):
      dh[1,1]:=n:
      func[1]:=d:
         for i from 1 to n2 do
            c[1,i]:=coeff(d,q,i);
         od;
   fi:
   fq:=series(f,q,n3):
   for i from 1 to n2 do c[nx+1,i]:=coeff(fq,q,i):od:
   ct:=  linalg[transpose](c):
   gooo:=(linalg[kernel](ct)):
   ##warning added 09/13/05
   if nops(gooo)>1 then
      print(`WARNING: dim ker =`,nops(gooo));
   fi:
   ANS:={}: 
   ANSETA:={}: 
   for k from 1 to nops(gooo) do
       goo:=gooo[k]:
       d:=op(dh):
       g[k]:=0:
       geta[k]:=0:
       if goo[nx+1]<>0 then
           for i from 1 to nx do
              g[k]:=g[k]+goo[i]*product(X[ii]^d[i,ii],ii=1..D):
              if etaoption = `yes` then
                geta[k]:=geta[k]+goo[i]*qseries[etamake](func[i],q,n3-2): 
              fi:
           od:
           ANS:=ANS union {-g[k]/goo[nx+1]}:
           if etaoption = `yes` then
             ANSETA:=ANSETA union {-geta[k]/goo[nx+1]}:
           fi:
        fi:
   od:
 if xprint then
 print(`-----possible linear combinations of degree------`,n):  
 fi:
if etaoption = `yes` then
  print(ANSETA);
fi:
RETURN(ANS);
else
   ERROR(`n must be a positive integer and topshift must apositive integer greater then -20`);
fi:
else
   ERROR(`n and topshift must be integers.`);
fi:
end:
#++++++++++++++++++++++++++++++++++++++++++++++++++++
 
qseries[findhomcombomodp]:=proc(f,L,p,q,n,topshift,etaoption)
#This program tries to express f (mod p) as a homogeneous polynomial     
#of degree n among the list of q-series in L.
#
#
   local xx,x,y,z,nx,n2,n3,gg,goo,gooo,c,ct,dh,jj,g,k,w,d,i,D,ii,h,G,ANS,func,fq,ANSETA,geta;
   global X:
   X:='X':
if whattype(n)=integer and whattype(topshift)=integer then
if n>0 and topshift>-20 then

   D:=nops(L);
   gg[1]:={};
   xx:=expand(sum(x[i],i=1..D)^(n) ):
   nx:=nops(xx):
   if D=1 then
      nx:=1:
   fi:
   n2:=nx+20+topshift:#step this up if some solutions are spurious
   n3:=n2+1:
   print(`# of terms `,n3);
   c:=array(1..nx+1,1..n2);
   dh:=array(1..nx,1..D);
   if D>1 then
      for jj from 1 to nx do
         G:=op(jj,xx):
         for i from 1 to D do
         h[i]:=degree(G,x[i]):
         dh[jj,i]:=h[i]:
         od:
         w:=product(L[ii]^h[ii],ii=1..D):
         d:=series(w,q,n3);
         func[jj]:=d:
         for i from 1 to n2 do
            c[jj,i]:=coeff(d,q,i);
         od;
      od:
   else
      d:=series(L[1]^n,q,n3):
      dh[1,1]:=n:
      func[1]:=d:
         for i from 1 to n2 do
            c[1,i]:=coeff(d,q,i);
         od;
   fi:
   fq:=series(f,q,n3):
   for i from 1 to n2 do c[nx+1,i]:=coeff(fq,q,i):od:
   ct:=  linalg[transpose](c):
   gooo:=Nullspace(ct) mod p:  
   print("Nullspace computed");
   if nops(gooo)>1 then
      print(`WARNING: dim ker =`,nops(gooo));
   fi:
   ANS:={}: 
   ANSETA:={}: 
   for k from 1 to nops(gooo) do
       goo:=gooo[k]:
       d:=op(dh):
       g[k]:=0:
       geta[k]:=0:
       if goo[nx+1]<>0 then
           for i from 1 to nx do
              g[k]:=g[k]+goo[i]*product(X[ii]^d[i,ii],ii=1..D):
              #added if to stop unnecessary slowdown using etamake
              if etaoption = `yes` then
                geta[k]:=geta[k]+goo[i]*qseries[etamake](func[i],q,n3-2): 
              fi:
           od:
           ANS:=ANS union {-g[k]/goo[nx+1]}:
           if etaoption = `yes` then
             ANSETA:=ANSETA union {-geta[k]/goo[nx+1]}:
           fi:
        fi:
   od:
   ANS:=modp(ANS,p):
 print(`-----possible linear combinations of degree------`,n,`---- mod `,p):  
if etaoption = `yes` then
  print(ANSETA);
fi:
RETURN(ANS);
else
   ERROR(`n must be a positive integer and topshift must apositive integer greater then -20`);
fi:
else
   ERROR(`n and topshift must be integers.`);
fi:
end:
#++++++++++++++++++++++++++++++++++++++++++++++++++++

qseries[findhommodp]:=proc(L,p,q,n,topshift)
# HISTORY:
##05/18/05 This is a modp version of findhom.
##10.30.15: changed taylor to series
#This program generates a list of potential homogeneous relations mod p
#of degree n among the list of q-series L.
#
   local xx,x,y,z,nx,n2,n3,gg,goo,gooo,c,dh,jj,g,k,w,d,i,D,ii,h,G,ANS;
   global X,ct:
if whattype(n)=integer and whattype(topshift)=integer then
if n>0 and topshift>-20 then
   D:=nops(L);
   gg[1]:={};
   xx:=expand(sum(x[i],i=1..D)^(n) ):
   nx:=nops(xx):
   n2:=nx+20+topshift:#step this up if some solutions are spurious
   n3:=n2+1:
   print(`# of terms `,n3);
   c:=array(1..nx,1..n2);
   dh:=array(1..nx,1..D);
   for jj from 1 to nx do
      G:=op(jj,xx):
      for i from 1 to D do
      h[i]:=degree(G,x[i]):
      dh[jj,i]:=h[i]:
      od:
      w:=product(L[ii]^h[ii],ii=1..D):
      d:=series(w,q,n3+5);
      for i from 1 to n2 do
         c[jj,i]:=coeff(d,q,i);
      od;
   od:
   ct:=  linalg[transpose](c):
   ##gooo:=(linalg[kernel](ct));
   gooo:=Nullspace(ct) mod p;  
   for k from 1 to nops(gooo) do
      goo:=gooo[k];
      d:=op(dh):
      g[k]:=0:
      for i from 1 to nx do
         g[k]:=g[k]+goo[i]*product(X[ii]^d[i,ii],ii=1..D):
      od;
      gg[k]:=g[k];
  od;
 print(`-----RELATIONS----of order---`,n,`---- mod `,p);  
   ANS:=convert(gg,set);
   RETURN(ANS);
else
   ERROR(`n must be a positive integer and topshift must apositive integer greater then -20`);
fi:
else
   ERROR(`n and topshift must be integers.`);
fi:
end:

qseries[findlincombo]:=proc(f,L,SL,q,topshift)
#04/27/08
#This program tries to express f as a linear combo of funcs in L 
#SL is list of names of funcs in L.          
#
   local fq, nx, n2, n3, c, jj, d, i, goo, ANS;
   global ct,gooo:
   #06/20/07: added gooo to global
if whattype(topshift)=integer then
if topshift>-20 then
   nx:=nops(L):
   print(`nx = `,nx);
   n2:=nx+20+topshift:#step this up if some solutions are spurious
   n3:=n2+1:
   print(`# of terms `,n3);
   c:=array(1..nx+1,1..n2);
      for jj from 1 to nx do
         d:=series(L[jj],q,n3);
   ## 02/14/09: changed loop  [i from 1 to n2] TO [for i from 0 to n2-1] 
         for i from 0 to n2-1 do
            c[jj,i+1]:=coeff(d,q,i);
         od;
      od:
   fq:=series(f,q,n3):
   ## 02/14/09: changed loop  [i from 1 to n2] TO [for i from 0 to n2-1] 
   for i from 0 to n2-1 do c[nx+1,i+1]:=coeff(fq,q,i):od:
   ct:=  linalg[transpose](c):
   gooo:=(linalg[kernel](ct)):
   ##warning added 09/13/05
   if nops(gooo)>1 or nops(gooo)=0 then
      print(`WARNING: dim ker =`,nops(gooo));
      if nops(gooo)=0 then
           print("NOT A LINEAR COMBO.");
      fi:
      RETURN():
   fi:
       goo:=gooo[1]:
       if goo[nx+1]<>0 then
           ANS:=add(-goo[k]/goo[nx+1]*SL[k],k=1..nx);
           RETURN(ANS);
       else
           print("NOT A LINEAR COMBO.");
           RETURN():
       fi:
else
   ERROR(`topshift must apositive integer greater then -20`);
fi:
else
   ERROR(`topshift must be integers.`);
fi:
end:
#####################################################################
##EXAMPLE:
##>    f1:=1+q+q^3;
##>    f2:=1+q-q^3-q^5-q^7;
##>    f3:=12*f1+13*f2;
##> qseries[findlinhomcombo](f3,[f1,f2],[F1,F2],q,0);
##                                    nx = , 2
##
##                                # of terms , 23
##
##                                 12 F1 + 13 F2
##
##> qseries[findlinhomcombo](f3,[f1,f2+q],[F1,F2],q,0);
##                                    nx = , 2
##
##                                # of terms , 23
##
##                             WARNING: dim ker =, 0
##
##                             "NOT A LINEAR COMBO."
##
##
#####################################################################
qseries[findlincombomodp]:=proc(f,L,SL,p,q,topshift)
#02/15/09 (This is a mod p version of findlincombo)
#This program tries to express f as a linear combo modp of funcs in L 
#SL is list of names of funcs in L.          
#
   local fq, nx, n2, n3, c, jj, d, i, goo, ANS;
   global ct,gooo:
   #06/20/07: added gooo to global
if whattype(topshift)=integer then
if topshift>-20 then
   nx:=nops(L):
   print(`nx = `,nx);
   n2:=nx+20+topshift:#step this up if some solutions are spurious
   n3:=n2+1:
   print(`# of terms `,n3);
   c:=array(1..nx+1,1..n2);
      for jj from 1 to nx do
         d:=series(L[jj],q,n3);
   ## 02/14/09: changed loop  [i from 1 to n2] TO [for i from 0 to n2-1] 
         for i from 0 to n2-1 do
            c[jj,i+1]:=coeff(d,q,i);
         od;
      od:
   fq:=series(f,q,n3):
   ## 02/14/09: changed loop  [i from 1 to n2] TO [for i from 0 to n2-1] 
   for i from 0 to n2-1 do c[nx+1,i+1]:=coeff(fq,q,i):od:
   ct:=  linalg[transpose](c):
   gooo:=Nullspace(ct) mod p: 
   ##warning added 09/13/05
   if nops(gooo)>1 or nops(gooo)=0 then
      print(`WARNING: dim ker =`,nops(gooo));
      if nops(gooo)=0 then
           print("NOT A LINEAR COMBO MOD ",p,".");
      fi:
      RETURN():
   fi:
       goo:=gooo[1]:
       if goo[nx+1]<>0 then
           ANS:=add(-goo[k]/goo[nx+1]*SL[k],k=1..nx);
           RETURN(ANS);
       else
           print("NOT A LINEAR COMBO MOD ",p,".");
           RETURN():
       fi:
else
   ERROR(`topshift must apositive integer greater then -20`);
fi:
else
   ERROR(`topshift must be integers.`);
fi:
end:
#####################################################################
#####################################################################
##EXAMPLE:
##
## > E4:=1+240*add(n^3*q^n/(1-q^n),n=1..300):
## > E6:=1-504*add(n^5*q^n/(1-q^n),n=1..300):
## > P:=series(1/etaq(q,1,300),q,300):       
## > ETA1:=etaq(q,1,300):
## > F1:=sift(E4^16*E6*P,q,13,6,299):
## > BB13:=[seq(ETA1^11*E6*E4^(3*j+1),j=0..5)]:
## > SBB13:=[seq(_ETA1^11*_E6*_E4^(3*j+1),j=0..5)]:
## > findlincombomodp(F1,BB13,SBB13,13,q,-8);
##                                             nx = , 6
## 
##                                         # of terms , 19
## 
##                                                11
##                                        -6 _ETA1   _E6 _E4
## 
##   Let E4(q), E6(q) be the standard Eisenstein series of
##   weights 4 and 6, and let E(q) = (q;q)[infinity].
##   Let F = sum a(n) q^n = E4^16*E6/E(q).
##   We suspect that
##   sum a(13*n+6) q^n = linear combination of
##   E(q)^11*E6*E4^(3*j+1) (j=0..5) mod 13.
##   Using findlincombomodp it seems that
##   sum a(13*n+6) q^n = -6*E(q)^11*E4(q)*E6(q) mod 13.
## 
##
#####################################################################
qseries[findmaxind]:=proc(XFL,T)
  local P,NXFL,j,NP,NNXFL:
  ## 11.28.17: findhom like fixed was findhom(NP,q,1,0,T)
  ## Created NOV 2015
  if member(nargs,{0,2}) then
     if nargs=0 then
       printf("-------------------------------------------------------------\n");
       printf("findmaxind(XFL,T)                                             \n");
       printf("   XFL is list of q-polynomials.           \n");
       printf("   T is an integer (ussually 0). Increase value of T if you \n");
       printf("   a problem with not computing enough coefficients.        \n");
       printf("   Returns [P,NXFL],   \n");
       printf("   where P is maximal independent subset of XFL, and \n");
       printf("   and NXFL is a list of indices.\n");
       printf("-------------------------------------------------------------\n");
     else
       P:=[XFL[1]]: NXFL:=[1]:
       for j from 2 to nops(XFL) do
          NP:=[op(P),XFL[j]]: NNXFL:=[op(NXFL),j]:
          if findhom(NP,q,1,T)={{}} then
             P:=NP: NXFL:=NNXFL: 
          fi:
       od:
       RETURN([P,NXFL]):
     fi:
  else
     ERROR(` number of arguments must be 0 or 2.`);
  fi:
end:

qseries[findnonhom]:=proc(L,q,n,topshift)
#This program generates a list of potential nonhomogeneous relations
#of degree n among the list of q-series L.
#
   local xx,x,y,z,nx,n2,n3,gg,goo,gooo,c,dh,jj,g,k,w,d,i,ct,D,ii,h,G,ANS;
   global X:
if whattype(n)=integer and whattype(topshift)=integer then
if n>0 and topshift>-20 then
   X:='X';
   D:=nops(L);
   gg[1]:={};
   xx:=expand((sum(x[i],i=1..D)+1)^(n) ):
   nx:=nops(xx):
   n2:=nx+20+topshift:#step this up if some solutions are spurious
   n3:=n2+1:
   print(`# of terms `,n3);
   c:=array(1..nx,1..n2+1);
   dh:=array(1..nx,1..D);
   for jj from 1 to nx do
      
      G:=op(jj,xx):
      
      for i from 1 to D do
         h[i]:=degree(G,x[i]):
         dh[jj,i]:=h[i]:
      od:
      w:=product(L[ii]^h[ii],ii=1..D):
      d:=series(w,q=0,n3+4);
      for i from 0 to n2 do
         c[jj,i+1]:=coeff(d,q,i);
      od;
      
   od:
   ct:=  linalg[transpose](c):
   gooo:=(linalg[kernel](ct));
  
   for k from 1 to nops(gooo) do
      goo:=gooo[k];
      d:=op(dh):
      g[k]:=0:
   for i from 1 to nx do
      g[k]:=g[k]+goo[i]*product(X[ii]^d[i,ii],ii=1..D):
   od;
      gg[k]:=g[k];
  od;
 print(`-----RELATIONS----of order---`,n);  
   ANS:=convert(gg,set);
   RETURN(ANS);
else
   ERROR(`n must be a positive integer and topshift must apositive integer greater then -20`);
fi:
else
   ERROR(`n and topshift must be integers.`);
fi:
end:

qseries[findnonhomcombo]:=proc()
#This program tries to express f as a polynomial     
#of degree n among the list of q-series in L.
#
   local xx,x,y,z,nx,n2,n3,gg,goo,gooo,c,dh,jj,g,k,w,d,i,ct,D,ii,h,G,ANSf,L,q,n,etaoption,
       topshift,f,func,fq,ANS,ANSETA,geta;
   global X:
   X:='X':
   f:=args[1]: L:=args[2]: q:=args[3]: n:=args[4]:
   if nargs>4 and whattype(args[5])=integer then
      topshift:=args[5]:
      if nargs=6 then
         etaoption:=args[6]:
      else
         etaoption:=no:      
      fi:
   fi:
   if nargs=4 then
      topshift:=0:
      etaoption:=no:
   fi:
   if nargs=5 and whattype(args[5])<>integer then
      topshift:=0:
      etaoption:=args[5]:
   fi:
  ##print(` n= `,n, `topshift =`,topshift); 
if whattype(n)=integer and whattype(topshift)=integer then
if n>0 and topshift>-20 then

   D:=nops(L);
   gg[1]:={};
   xx:=expand((sum(x[i],i=1..D)+1)^(n) ):
   nx:=nops(xx):
   n2:=nx+20+topshift:#step this up if some solutions are spurious
   n3:=n2+1:
   print(`# of terms `,n3);
   c:=array(1..nx+1,1..n2+1);
   print(`matrix is `,nx+1,` x `,n2+1);
   dh:=array(1..nx,1..D);
   for jj from 1 to nx do
      G:=op(jj,xx):
      for i from 1 to D do
      h[i]:=degree(G,x[i]):
      dh[jj,i]:=h[i]:
      od:
      w:=product(L[ii]^h[ii],ii=1..D):
      ##d:=taylor(w,q=0,n3+4);
      ##09-17-10: changed taylor to series
      d:=series(w,q=0,n3+4);
      func[jj]:=d:
      for i from 0 to n2 do
         c[jj,i+1]:=coeff(d,q,i);
      od;
   od:
   fq:=series(f,q,n3+4):
   for i from 0 to n2 do 
      c[nx+1,i+1]:=coeff(fq,q,i):
   od:
   ct:=  linalg[transpose](c):
   gooo:=(linalg[kernel](ct)):
   ANS:={}: 
   ANSETA:={}: 
   for k from 1 to nops(gooo) do
       goo:=gooo[k]:
       d:=op(dh):
       g[k]:=0:
       geta[k]:=0:
       if goo[nx+1]<>0 then
           for i from 1 to nx do
              g[k]:=g[k]+goo[i]*product(X[ii]^d[i,ii],ii=1..D):
              if etaoption=`yes` then
              geta[k]:=geta[k]+goo[i]*qseries[etamake](func[i],q,n3-2): 
              fi:
           od:
           ANS:=ANS union {-g[k]/goo[nx+1]}:
              if etaoption=`yes` then
           ANSETA:=ANSETA union {-geta[k]/goo[nx+1]}:
              fi:
        fi:
   od:
 print(`-----possible linear combinations of degree------`,n):  
if etaoption = `yes` then
  print(ANSETA);
fi:
RETURN(ANS);
else
   ERROR(`n must be a positive integer and topshift must apositive integer greater then -20`);
fi:
else
   ERROR(`n and topshift must be integers.`);
fi:
end:
qseries[findpoly]:=proc() #This proc looks for a polynomial
       local x,y,q,deg1,deg2,ARGLIST,TYPELIST,CORRECTL,check,dim1,dim2,A,num,k,j,B,qq,l,kk,POLY,kkk,i,POLYg,POLYfunc,ss,numtyp,tA,opkk:
       global X,Y:
       lprint(`WARNING: X,Y are global.`);
       #relation between x,y of degree deg1 in x, and degree deg2 in y.
       #The relation is checked to order O(q^check).
       if nargs<5 then
          ERROR(` number of arguments must be 5 or 6.`);
       fi:
       x:=args[1]:
       y:=args[2]:
       q:=args[3]:
       deg1:=args[4]:
       deg2:=args[5]:
       ARGLIST:=[q,deg1,deg2];
       TYPELIST:=map(whattype,ARGLIST);
       ##05/15/99: change string->symbol here.
       ##OLD: CORRECTL:=[string,integer,integer]:
       CORRECTL:=[symbol,integer,integer]:
       if TYPELIST=CORRECTL then
         if nargs=6 then
            check:=args[6]:
         fi:
         if nargs>6 then
            ERROR(` findpoly can at most 6 arguments`);
         fi:
	 dim1:=(deg1+1)*(deg2+1):
	 dim2:=dim1+10:
	 print(` dims `, dim1, dim2);
	 A:=array(1..dim1,1..dim2):
	 num:=0:
	 for k from 0 to deg1 do
            for j from 0 to deg2 do
	       num:=num+1:
               B[num]:=X^k*Y^j:
	       qq:=series(x^k*y^j,q,dim2+3):
               for l from 0 to (dim2-1) do
	          A[num,l+1]:=coeff(qq,q,l):
	       od:
	     od:
	 od:
	#removed 10/19/06: with(linalg):
	tA:=linalg[transpose](A):
	kk:=linalg[kernel](tA):
	if kk={} then 
	  lprint(` NO polynomial relation found. `);
	else	
	 POLY:=0;
	 opkk:=op(kk);
	 kkk:=opkk;
	 for i from 1 to num do
	   POLY:=POLY+B[i]*kkk[i]:
	 od:
	 POLYg:=0;
	for j from 0 to deg2 do
	  POLYg:=POLYg+factor(coeff(POLY,Y,deg2-j))*Y^(deg2-j):
	od:
        if nargs=6 then
          print(`The polynomial is`);
	  print(POLYg);
	  POLYfunc:=unapply(POLYg,X,Y):
	  ss:=series(POLYfunc(x,y),q,check+1):
          print(`Checking to order`,check);
	  print(ss);
	  RETURN(POLYg):
        else 
          print(`The polynomial is`);
	fi:
        RETURN(POLY):
        fi:
      else
        numtyp:=nops(TYPELIST):
        ERROR(` Wrong type of argument. ARGLIST has type `,seq(TYPELIST[i],i=1..numtyp),
        `It should have type`,seq(CORRECTL[i],i=1..numtyp));
      fi:
	end:
qseries[findprod]:=proc(FL,T,M,Q)      ## find z-linear combination of the
## 09.25.15: This function uses checkprod.
## 03.06.19: removed print statement
        local COMBS,nfl,S,BS,CP,vec,mcomb,cc,LL:
        COMBS:=[]:            ## functions in FL which is a prob.product
	nfl:=nops(FL):        ## T=max(abs) of coeffs in z-lin.combo
        S:=[seq(j,j=-T..T)]:  ## M=max(abs(multiplicity))) in prod
        BS:=[seq(S,k=1..nfl)]:## Q is highest power of q
        CP:=combinat[cartprod](BS):
        while not CP[finished] do
           vec:=CP[nextvalue]();
           if igcd(op(vec))=1 then
           mcomb:=add(vec[k]*FL[k],k=1..nfl):
           cc:=qseries[checkprod](mcomb,M,Q):
           ##print(vec,cc);
             if cc[2]=1 then
                ##print(cc);
                LL:=cc[1]:
                COMBS:=[op(COMBS),[LL,vec]]:
             fi:
           fi:
        od;
RETURN(COMBS):
end:
qseries[J2jaclist]:=proc()
## EXAMPLE: 
## > with(qseries):
## >    z1:=J[25, 5]*J[25]*J[50, 15]/J[25, 10];
##                                J[25, 5] J[25] J[50, 15]
##                          z1 := ------------------------
##                                       J[25, 10]
## 
## >    J2jaclist(z1);
##                 [[25, 5, 1], [25, 1], [50, 15, 1], [25, 10, -1]]
## 
## >    x1:=J[25,1];
##                                  x1 := J[25, 1]
## 
## >    J2jaclist(x1);
##                                   [[25, 1, 1]]
## 
## >    x2:=J[25]^2;
##                                              2
##                                   x2 := J[25]
## 
## >    J2jaclist(x2);
##                                    [[25, 2]]
## 
local L1,L1n,JP,i,r,L1i,s,chk;
global Jprod:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("J2jaclist(Jprod)                                   \n");
     printf("   Jprod is a quotient of theta-functions encoded\n");
     printf("   in terms of J[b,a], J[b];                                 \n");
     printf("   The quotient is converted to a list.\n");
     printf("   Each term J[b,a]^c or J[b]^c in the quotient is converted to an\n");
     printf("   item [b,a,c] or [b,c] in the list.\n");
     printf("-------------------------------------------------------------\n");
  elif nargs = 1 then
     Jprod:=args[1]:
	        if whattype(Jprod)=`^` then
		   L1:=[Jprod]:
	     	else
		   if whattype(Jprod)=indexed then
                       if op(0,Jprod)='J' then
			RETURN([[op(Jprod),1]]);
                       else
                        ERROR(`not a J-prod`);
                       fi:
		   else
		        L1:=[op(Jprod)]:
	           fi:
		fi:
		L1n:=nops(L1):
		JP:=[]:   
		for i from 1 to L1n do
                    r:=degree(L1[i]):
	            L1i:=L1[i]:
	            if whattype(L1i)=`^` then
 	               r:=op(2,L1i):
		       s:=op(op(1,L1i)):
	            else
		       r:=1:
		       s:=op(L1i):
	            fi:
		    #check:
		    ##print("L1i=",L1i,"J[s]^r=",J[s]^r);
	            chk:=simplify(L1i-J[s]^r):
	            if chk<>0 then
		      error "chk<>0";
	            fi:  
		    JP:=[op(JP),[s,r]]:
		od:
		RETURN(JP);
  else
     ERROR(`invalid input type`);
  fi:
	end:	
#######################################################################
#  qseries[jac2J]
#######################################################################

macro(jacJ = `jac2J/jacJ`);

jacJ:=proc(a::integer, b::integer, c::infinity)
    local Q;
    if a<b then
          if a>0 then
	     Q:=J[b,a];
             RETURN(Q):
          else
             RETURN(J[b]);
          fi:
    else
       ERROR(` a,b must satisy a<b `);
       RETURN():
    fi:
end:
#######################################################################
#  qseries[jac2J]
#######################################################################
#######################################################################
# FUNCTION :  jac2J    - converts a product of theta functions into
#                          J-product form.
# CALLING SEQUENCE :  jac2prod(jacexpr)
# PARAMETERS :   jacexpr   - product of JAC(i,j,oo)
#                            where i,j are integers
#                            and 0 < i < j.
#                                  =
# SYNOPSIS :   
# jac2prod(jacexpr) returns a J-product;ie a product of the
# of functions of the form J[b,a] or J[b].
#######################################################################
qseries[jac2J]:=proc()
   local JJ;
 if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf(" jac2J(JEXP)                                         \n");
     printf(" jac2J is a variant of jac2prod.  \n");
     printf(" It converts JAC notation to J[b,a] and J[b]\n");
     printf("---------------------------------------------------------\n");
else
   if nargs=1 then
      JJ:=args[1]:
      RETURN(eval(subs(JAC=jacJ,JJ))):
   else
     ERROR(`jac2J takes 0 or 1 args`);      
  fi:
fi:
end:

macro(jacJ = jacJ):     
#######################################################################
#  qseries[jac2prod]
#######################################################################

macro(jac = `jac2prod/jac`);

jac:=proc(a::integer, b::integer, c::infinity)
    local Q;
    if a<b then
          if a>0 then
Q:=` `(q^a,q^b)[infinity]*` `(q^(b-a),q^b)[infinity]*` `(q^b,q^b)[infinity];
             RETURN(Q):
          else
             RETURN(` `(q^b,q^b)[infinity]);
          fi:
    else
       ERROR(` a,b must satisy a<b `);
       RETURN():
    fi:
end:
#######################################################################
#  qseries[jac2prod]
#######################################################################
#######################################################################
# FUNCTION :  jac2prod   - converts a product of theta functions into
#                          q-product form.
# CALLING SEQUENCE :  jac2prod(jacexpr)
# PARAMETERS :   jacexpr   - product of JAC(i,j,oo)
#                            where i,j are integers
#                            and 0 < i < j.
#                                  =
# SYNOPSIS :   
# jac2prod(jacexpr) returns a q-product;ie a product of the
#                            a   b
# of functions of the form (q , q ).
#######################################################################
qseries[jac2prod]:=proc(JEXPR)
   local JJ;
   JJ:=JEXPR:
   RETURN(eval(subs(JAC=jac,JJ))):
end:

macro(jac = jac):     
#######################################################################
#  qseries[jac2series]
#######################################################################

macro(tripleprod2 = `jac2series/tripleprod2`, qjac = `jac2series/qjac`);

tripleprod2:=proc(z,q,T)
    local x,lasti,i;
        x := 0;
        lasti := T;
        for i from -lasti to lasti do
            x := x+(-1)^i*z^i*q^(1/2*i*(i-1))
        od;
        RETURN(x)
end:
#########
##
##
##OCT 24, 1997 - fixed a bug in qjac. Previously it called tripleprod
##               which did not give enough terms. This been a replaced
##               by a variant called tripleprod2 (see above).
##
#########
qjac:=proc(a::integer, b::integer, T::integer)
    local Q,t,c;
    option remember:
    ##08.27.14: now can handle a<0 and a>b
    if T>0 then
          if a<>0 then
             if a<0 then
                RETURN(qjac(-a,b,T+b)*(-q^(a))):
             else
                if a>b then
                   RETURN(qjac(a-b,b,T+b)*(-q^(b-a))):
                else
                   t:=max(T,abs(b)):
                   ##c:=trunc(1/2*(b-2*a+(4*a^2-4*a*b+b^2+8*b*t)^(1/2))/b)+1: 
	           ##Nov 15, 2002 added evalf to trunc
                   c:=trunc(evalf(1/2*(b-2*a+(4*a^2-4*a*b+b^2+8*b*t)^(1/2))/b))+1: 
                   RETURN(tripleprod2(q^a,q^b,c));   
                 fi:
             fi:
          else
             RETURN(qseries[etaq](q,b,T));
          fi:
    else
       ERROR(`  T must be positive. `);
       RETURN():
    fi:
end:
#######################################################################
#  qseries[jac2series]
#######################################################################
#######################################################################
# FUNCTION :  jac2series   - converts a product of theta functions into
#                            q-product form.
# CALLING SEQUENCE :  jacseries(jacexpr)
# PARAMETERS :   jacexpr   - product of JAC(i,j,oo)
#                            where i,j are integers
#                            and 0 < i < j.
#                                  =
# SYNOPSIS :   
# jacprod(jacexpr) returns a q-product;ie a product of the
#                            a   b
# of functions of the form (q , q ).
#######################################################################
qseries[jac2series]:=proc(JEXPR, T::integer)
   local JJ;
   JJ:=JEXPR:
   RETURN(eval(subs({JAC=qjac,infinity=T},JJ))):
end:

macro(tripleprod2=tripleprod2, qjac = qjac):     
qseries[jaclist2JACPROD]:=proc()
## EXAMPLE: 
## > with(qseries):
## >    z1:=J[25, 5]*J[25]*J[50, 15]/J[25, 10]:
## >    L1:=J2jaclist(z1):
## >    jaclist2JACPROD(L1);
##         JAC(5, 25, infinity) JAC(0, 25, infinity) JAC(15, 50, infinity)
##         ---------------------------------------------------------------
##                              JAC(10, 25, infinity)
## 
## >    x1:=J[25,1]:
## >    L2:=J2jaclist(x1):
## >    jaclist2JACPROD(L2);
##                               JAC(1, 25, infinity)
## 
## >    x2:=J[25]^2:
## >    L3:=J2jaclist(x2):
## >    jaclist2JACPROD(L3);
##                                                  2
##                              JAC(0, 25, infinity)
## 
local jaclist,xj,L;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("jaclist2JACPROD(jaclist,T)                                   \n");
     printf("  jaclist is a list of terms [b,a,c] or [b,c]   \n");
     printf("  corresponding to J[b,a]^c or J[b]^c                      \n");
     printf("  These are converted to JAC(a,b,infinity)^c or JAC(0,b,infinity)^c");
     printf("-------------------------------------------------------------\n");
  elif nargs = 1 then
     jaclist:=args[1]:
     xj:=1:
     for L in jaclist do
        if nops(L)=2 then
             xj:=xj*JAC(0,L[1],infinity)^L[2]:
        else
             xj:=xj*JAC(L[2],L[1],infinity)^L[3]:
        fi:
     od:
     RETURN(xj):
 fi:
end:
qseries[jacprod]:=proc(a,b,q,T)
	local x,i,lasti:
	tripleprod(q^a,q^b,T)/tripleprod(q^(b),q^(3*b),T):
end:
#######################################################################
#  qseries[jacprodmake]
#######################################################################

macro(jaccheck = `jacprodmake/jaccheck`,
      periodfind = `jacprodmake/periodfind`,
      periodfind2 = `jacprodmake/periodfind2`,
      jacmake = `jacprodmake/jacmake`);


#######################################################################
# FUNCTION :  periodfind - finds the smallest period of a finite list
# CALLING SEQUENCE :  periodfind(A,T)
# PARAMETERS :   A   - list of integers
#                P   - positive integer
# SYNOPSIS :   
# periodfind(A,T) returns P if 
#      A[i] = A[i+P]   for each  i. 
#######################################################################
periodfind:=proc(A,T)
   local num,n,perfi,pp,dud,i,j,lastj,diffc;
   num:=T:       
   n:=trunc(num/2):
   perfi:=0:
   pp:=2:
   while perfi=0 do
     dud:=0:
     while dud=0 do
      for i from 1 to pp do
        lastj:=trunc(num/pp)-2:
        for j from 0 to lastj do
           diffc:=A[i+pp*j]-A[i+pp*j+pp]:
           if diffc<>0 then
              dud:=1:
           fi: 
        od:
      od:
      if dud=0 then
         perfi:=1:
         RETURN(pp):
         dud:=1:
      fi:
     od:
      pp:=pp+1:
  od:
end;
periodfind2:=proc(A,T,PP)
   global xprint:
   local num,n,perfi,pp,dud,i,j,lastj,diffc,pset,k,np;
   num:=T:       
   n:=trunc(num/2):
   perfi:=0:
   pset:=numtheory[divisors](PP) minus {1}:
   pp:=pset[1]:
   if xprint then
     lprint("TRYING period pp=",pp);
   fi:
   k:=1:
   np:=nops(pset):
   while perfi=0 do
     dud:=0:
     while dud=0 do
      for i from 1 to pp do
        lastj:=trunc(num/pp)-2:
        for j from 0 to lastj do
           diffc:=A[i+pp*j]-A[i+pp*j+pp]:
           if diffc<>0 then
              dud:=1:
           fi: 
        od:
      od:
      if dud=0 then
         perfi:=1:
         if xprint then
           lprint("PERIOD FOUND: ",pp);
         fi:
         RETURN(pp):
         dud:=1:
      fi:
     od:
      if k < np then
        k:=k+1:
       pp:=pset[k]:
       if xprint then
         lprint("TRYING period pp=",pp);
       fi:
     else
       RETURN( trunc(T/2)+2 ):
       dud:=1:
     fi:
  od:
end;
#######################################################################
# FUNCTION :  jaccheck - checks whether an array corresponds a 
#                        theta-product.
# CALLING SEQUENCE :  jaccheck(A,P):
# PARAMETERS :   A   - list of integers
#               P   - positive integer
# SYNOPSIS :   
# jaccheck(A,P) returns 1 if
#     A[i] = A[P-i]   for each  i, 1 < i < (P/2);
# otherwise 0 is returned. If jaccheck(A,P)=1 and A is periodic
# with period P, then A corresponds to a certain theta product.
#######################################################################
jaccheck:=proc(A,P)
     local n,jac,i:
     n:=trunc(P/2):
     jac:=1:
     for i from 1 to n do
         if A[i]<>A[P-i] then
            jac:=0:
            RETURN(0):
         fi:
     od:
     if jac=1 then
        RETURN(1):
     fi:  
     RETURN(-1):
end:
#######################################################################
# FUNCTION :  jacmake  - connverts a list of integers into a 
#                        theta-product if possible.
# CALLING SEQUENCE :  jacmake(A,T):
#                     jacmake(A,T,P):
# PARAMETERS :   A   - list of integers
#                T   - positive integer
#                P   - positive integer
# SYNOPSIS :                                    T
# If f = product((1-q^i)^(-A[i]),i=1..T-1) + O(q )
# then
# jacmake(A,T) will return a product of the form
#
#    product( JAC(j,B,oo), j in S)
#                                     T
# if this product agrees with f to O(q ).
#
# Here JAC(j,B,oo) corresponds to the theta-function
# tripleprod(q^j,q^B,oo) when  0 < j < B
# and
# JAC(0,B,oo) corresponds to the eta-product etaq(q,B,oo).
# jacmake(A,T,P) does the same except it assumes the base (period) B is
# a divisor of P.
#######################################################################
jacmake:=proc()
   # 04/23/10 NEW VERSION
   global xprint,jacp:
   local A,T,PP,pp,y,mp,i,y1,y2:
   # 12/17/12: jacp:=[jaccheck(A,pp),pp]
   #           if jaccheck=1 then pp=theta period
   if xprint then
     lprint("jacmake started with args",args);
   fi:
   if nargs=2 or nargs=3 then
     # OK
   else
      ERROR("jacmake takes 2 or 3 args");
   fi:
   A:=args[1]:
   T:=args[2]:
   if nargs=3 then
     PP:=args[3]:
   fi:
   if nargs=2 then
      pp:=periodfind(A,T):
   else
      pp:=periodfind2(A,T,PP):
   fi:
   if nargs=3 then
     if  modp(PP,pp)<>0 then
     print(`ERROR: PP should be multiple of pp`);
     RETURN():
     fi:
   fi:
   if pp>(T/2) then
     print(` no period found `);
     print(` This means either that there is no period,`);
     print(` or a larger value of T should be chosen.`);
     jacp:=[0,pp]:
     RETURN():
   else
     ## 12/17/12: added global jacp
     jacp:=[jaccheck(A,pp),pp]:
     if jaccheck(A,pp)<>1 then
        print(` not a theta function `);
        RETURN():          
     else
     ##if nargs=3 then
     ##  pp:=PP:
     ##fi:
     y:=1:
     mp:=modp(pp,2):
     if mp=0 then
        for i from 1 to ((pp)/2-1) do
           y:=y/(JAC(i,pp,infinity)/JAC(0,pp,infinity))^A[i]:
        od:
        ##OLD:y:=y/(JAC(pp/2,pp,infinity)/JAC(0,pp,infinity))^(A[pp/2]/2):
        y1:=JAC(pp/2,pp,infinity):
        y2:=JAC(0,pp,infinity):
        y:=y/y1^(A[pp/2]/2)*y2^(A[pp/2]/2):
     else
        for i from 1 to (pp-1)/2 do
           y:=y/(JAC(i,pp,infinity)/JAC(0,pp,infinity))^A[i]:
        od:
      fi:
      y:=y/JAC(0,pp,infinity)^A[pp]:
      RETURN(y):
      fi:
      fi:
end:
#######################################################################
#  THE MAIN FUNCTION qseries[jacprodmake]
#######################################################################
        # 04/23/10 NEW VERSION
        # 05/03/10 added global var xprint
	qseries[jacprodmake]:=proc()
          global xprint,ft0,fixjacp,jacp:
          local f,q,T,ft,f0,_a,_b,i,n,j,d,_B,sum1,sum2,divj,divjb,m,_A,LT,PP,newT,je;
          fixjacp:=args:
          if not(type(xprint,boolean)) then
               xprint:=false:
          fi:
          if nargs=3 or nargs=4 then
            # OK
          else
            ERROR("nargs in jacprodmake should be 3 or 4");
          fi:
          f:=args[1]:
          q:=args[2]:
          T:=args[3]:
          if nargs=4 then
            PP:=args[4]:
          fi:
          ##05/03/10 added leadterm stuff
	  ft0:=series(f,q,T+5):
          if nops(ft0)=2 and op(1,ft0)=O(1) then
             newT:=trunc(op(2,ft0))+5:
	     ft:=series(f,q,newT):
          else
             ft:=ft0:
             newT:=T+5:
          fi:
          LT:=series(`leadterm`(ft),q,newT):
          LT:=convert(LT,polynom):
	  if whattype(ft)=series  then
	   f0:=coeff(ft,q,0):
	   ##if f0=1 then
	      _b:=series(ft/LT,q,newT):
	      _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:=sum1+d*(_A[d]):
                    od:
                    sum2:=sum2 + (_B[n-j])*sum1:
                 od:
                 sum1:=0:
	            divjb:=numtheory[divisors](j) minus {n}:
                    for d in divjb do
	               sum1:=sum1+d*(_A[d]):
	            od:
	         sum2:=n*_B[n] - sum2 - sum1:
	         _A[n]:=sum2/n:
	      od:
              if nargs=3 then
                 je:=jacmake(_A,T-1):
              else
                 if xprint then
                   lprint("_A done");
                 fi:
                 je:=jacmake(_A,T-1,PP):
              fi:
              #05/04/10: added error and LT stuff
              if nops([je])=0 then
                 fixjacp:=args:
                 #12/12/12: changed following error to warning
                 printf("WARNING: jacprodmake failed (see gobal var fixjacp)");
                 RETURN():
              else
                 RETURN(LT*je):
              fi:
	   ##else
	   ##   ERROR(`coeff of q^0 must be 1`);
	   ##fi:
	else
	   ERROR(`f must be a series`);
	fi:
	end:

macro(jaccheck = jacprod,
      periodfind = periodfind,
      periodfind2 = periodfind2,
      jacmake = jacmake);

qseries[Jetamake]:=proc() 
#
#10.19.19: This is version of qetamake which returns eta-products using
#          J-notation; ie. eta(p*tau) <-> J[p]
#This is a version of etamake which returns _E(q^m) instead of eta(m*tau)
#This procedure computes an eta-product expansion
#of the series f to order O(q^last).              
	   local f,last,fp,tc,exq,g,aa,ld,h,hh,i,cf,Jetaprod,alast,sf:
           global q:
if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf(" Jetamake(f,q,T)                                    \n");
     printf(" qetamake is a variant of etamake.   \n");
     printf(" etamake(f,q,T) converts the q-series f into an eta-product\n");
     printf(" expansion that agrees with f to O(q^T).   \n");
     printf(" Jetamake returns a product of the functions of the form J[m]\n");
     printf(" where J[m] = (q^m;q^m)oo = (1-q^m)*(1-q^(2m))*(1-q^(3m))...   \n");
     printf(" In other words, J[m] = eta(m*tau)/q^(m/24).   \n");
     printf("---------------------------------------------------------\n");
else
	f:=args[1]:
	q:=args[2]:
	last:=args[3]:
           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:
	   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
	         ##aa:=_E(q^ld)^(-cf)*aa:
	         aa:=J[ld]^(-cf)*aa:
	         g:=g*etaq(q,ld,alast)^cf:
	      fi:
	   od:
	   ##etaprod:=aa
           ##Fixed 05/18/05
	   Jetaprod:=aa*q^exq:
	   RETURN(Jetaprod):
fi:
	end:
qseries[Jterm2JACPROD]:=proc()
  local LJT,selcon,jterm2,JL,JP,jterm,dilfactor:
#############################################################################
## > y1:=2*J[25, 5]*J[25]*J[50, 15]/J[25, 10];
##                               2 J[25, 5] J[25] J[50, 15]
##                         y1 := --------------------------
##                                       J[25, 10]
## 
## > Jterm2JACPROD(y1,1);
##        2 JAC(5, 25, infinity) JAC(0, 25, infinity) JAC(15, 50, infinity)
##        -----------------------------------------------------------------
##                              JAC(10, 25, infinity)
## 
## > Jterm2JACPROD(y1,1/5);
##          2 JAC(1, 5, infinity) JAC(0, 5, infinity) JAC(3, 10, infinity)
##          --------------------------------------------------------------
##                               JAC(2, 5, infinity)
#############################################################################
 if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf(" Jterm2JACPRPOD(jterm,dilfactor)                         \n");
     printf("   Converts jterm to JACPROD and also does q -> q^dilfactor\n");
     printf("                                                         \n");
     printf(" EXAMPLE: \n");
     printf(" y1:=2*J[25, 5]*J[25]*J[50, 15]/J[25, 10];\n");
     printf(" Jterm2JACPROD(y1,1);\n");
     printf(" Jterm2JACPROD(y1,1/5);\n");
     printf("---------------------------------------------------------\n");
else
  if nargs=2 then
     jterm:=args[1]:
     dilfactor:=args[2]:
     if type(jterm,`*`) then
        LJT:=[op(jterm)]:
        selcon:=select(x->if type(x,constant) then true else false fi, LJT):
        if nops(selcon)=0 then 
           jterm2:=jterm:
        else
           jterm2:=expand(jterm/selcon[1]):
        fi:
        JL:=dilatejaclist(J2jaclist(jterm2),dilfactor):
        JP:=2*jaclist2JACPROD(JL):
        RETURN(JP):
     else
        RETURN(jterm):
     fi:
  else
     ERROR("nargs = 0 or 2 ");
  fi:
fi:
end:

qseries[lqdegree]:=proc()
 local FT,FTC,FT1,FT2,M1,MS1,qexp:
##> Y:=1/2*JAC(1,2,infinity) + q*JAC(3,4,infinity) + JAC(1,5,infinity)/q^(11/2);
##                                                         JAC(1, 5, infinity)
##  Y := 1/2 JAC(1, 2, infinity) + q JAC(3, 4, infinity) + -------------------
##                                                                 11/2
##                                                                q
##> lqdegree(Y);
##                                     -11/2
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("lqdegree(qexp)                                               \n");
     printf("   qexp is a sum of q-monomials.           \n");
     printf("   Degrees in q are allowed to be negative and fractional. \n");
     printf("   Returns the lowest degree in q that occurs.  \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 1 then
       qexp:=args[1]:  
       if convert(qexp,polynom) <>FAIL and ldegree(convert(qexp,polynom),q)<> FAIL 
       then
         ldegree(convert(qexp,polynom),q);
       else
         FT:=qexp:
         M1:=map(lqdegree0,[op(FT)]):
         MS1:=convert(M1,set) minus{FAIL}:
         min(MS1);
       fi:
  fi:
end:
qseries[lqdegree0]:=proc()
 local FT,FTC,FT1,FT2,qexp:
#########################################################################
##> B:=sqrt(q)*JAC(0,1,infinity);
##                               1/2
##                         B := q    JAC(0, 1, infinity)
##> lqdegree0(B);
##                                      1/2
##> B1:=(q)*JAC(0,1,infinity);   
##                          B1 := q JAC(0, 1, infinity)
##> lqdegree0(B1);            
##                                       1
##> B2:=JAC(0,1,infinity);    
##                           B2 := JAC(0, 1, infinity)
##> lqdegree0(B2);        
##                                       0
##> B3:=JAC(0,1,infinity)/q^3;
##                                 JAC(0, 1, infinity)
##                           B3 := -------------------
##                                          3
##                                         q
##> lqdegree0(B3);            
##                                      -3
##> B4:=JAC(0,1,infinity)/q^(1/3);
##                                 JAC(0, 1, infinity)
##                           B4 := -------------------
##                                         1/3
##                                        q
##> lqdegree0(B4);                
##                                     -1/3
#########################################################################
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("lqdegree0(qexp)                                       \n");
     printf("   qexp is q-monomial.           \n");
     printf("   This function returns the degree in q.\n");
     printf("   It also handles fractional degrees.                       \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 1 then
       qexp:=args[1]:  
     if ldegree(qexp,q)<>FAIL then
       ldegree(convert(qexp,polynom),q);
     else
       FT:=qexp:
       FTC:=subs(q=1,FT):
       FT1:=FT/FTC:
       FT2:=[op(FT1)]:
       if nops(FT2)=2 and FT2[1]=q then
         FT2[2];
       else
         RETURN(FAIL);
       fi:
     fi:
  fi:
end:
qseries[mprodmake]:=proc(f,q,last) 
#05/18/05
#This procedure computes a product expansion of the form
# (1+q^n1)*(1+q^n2)....
#of the series f to order O(q^last).              
	   local fp,tc,exq,g,aa,ld,h,hh,i,cf,etaprod,alast,sf,qqprod:
           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:
	   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
	         aa:=(1+q^ld)^(cf)*aa:
	         g:=g*(1+q^ld)^(-cf):
	      fi:
	   od:
	   qqprod:=q^(exq)*aa:
	   RETURN(qqprod):
	end:
	qseries[oldsift]:=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:
qseries[packageversion]:=proc()
printf("****************************************************\n");
printf("*\n");
printf("*  qseries package version 1.3f  \n");
printf("*  Sat, Feb  8, 2020 10:33:00 PM\n");
printf("*\n");
printf("*  This version tested on MAPLE 2019\n"); 
printf("*\n");
printf("*  Please report any problems to fgarvan@ufl.edu\n");
printf("*  See \n");
printf("*  http://qseries.org/fgarvan/qmaple/qmaple.html \n");
printf("*  for documentation and help.\n");
printf("*\n");
printf("*  Previous versions:\n");
printf("        1.3f -Feb 2020 (MAPLE 2019)\n");
printf("        1.3e -Jan 2020 (MAPLE 2019)\n");
printf("        1.3d -Oct 2019 (MAPLE 2019)\n");
printf("        1.3c -Jul 2019 (MAPLE 2019)\n");
printf("        1.3b -Mar 2019 (MAPLE 2018)\n");
printf("        1.3a -Nov 2017 (MAPLE 2017)\n");
printf("        1.3 - Aug 2016 (MAPLE 2015)\n");
printf("        1.2 - Dec 2012 (MAPLE 16)\n");
printf("        1.1 - Ju1 2012 (MAPLE 13)\n");
printf("        1.0 - Jun 2009 (MAPLE 10)\n");
printf("        0.9 - Apr 2008 (MAPLE 10)\n");
printf("        0.8 - May 2005 (MAPLE 9)\n");
printf("        0.7 - Mar 2004\n");
printf("        0.6 - Nov 2002\n");
printf("        0.5 - May 2000\n");
printf("        0.4 - Jan 2000\n");
printf("        0.3 - Nov 1999\n");
printf("        0.2 - Dec 1998\n");
printf("        0.1 - Dec 1997\n");
printf("****************************************************\n");
RETURN():
end:

	qseries[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:
qseries[qbin]:=proc(q,m,n)
        #06/21/09: allow q to be root of unity etc
        if whattype(m)=integer and whattype(n)=integer then
	  if m>=0 and m<=n then
            if type(q,symbol) then
	      RETURN(normal(aqprod(q,q,n)/aqprod(q,q,m)/aqprod(q,q,n-m)));
            else
	      RETURN(subs(_x=q,normal(aqprod(_x,_x,n)/aqprod(_x,_x,m)/aqprod(_x,_x,n-m))));
            fi: 
          else
            RETURN(0);
          fi:
          else
             ERROR(` m and n must be integers.`); 
        fi:
end:
 
 qseries[qdegree]:=proc(QS)
       # 10/14/10
       RETURN(degree(convert(QS,polynom),q));
 end:
qseries[qetamake]:=proc(f,q,last) 
#
#03/16/04
#Bug fixed 05/18/05
#This is a version of etamake which returns _E(q^m) instead of eta(m*tau)
#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:
           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:
	   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
	         aa:=_E(q^ld)^(-cf)*aa:
	         g:=g*etaq(q,ld,alast)^cf:
	      fi:
	   od:
	   ##etaprod:=aa
           ##Fixed 05/18/05
	   etaprod:=aa*q^exq:
	   RETURN(etaprod):
	end:
#######################################################################
#  qseries[qfactor]
#######################################################################
##01/08/00: 
##BUGS FIXED: 
## (1) Fix problem when numerator=1
## Example OLD version: 
## >    qfactor(1/(1-q));
## Error, (in qfactor/bigqpfac) unable to convert
##
## Example NEW version: 
## >    qfactor(1/(1-q));
##                                         1
##                                     - ------
##                                       -1 + q
## (2) Fix problem when input is a polynom
## Example OLD version: 
## >    qfactor(1+q+q^2);
##                                    2
##                                   q  + q + 1
## Example NEW version: 
## >    qfactor(1+q+q^2);
##                                          3
##                                     1 - q
##                                     ------
##                                     1 - q
##NOTE: Problems fixed by fixing qfopexp.
##
#######################################################################
##10/24/99: did maple-update and maple-local-make
#######################################################################
##10/19/99:  BUGS FIXED:  (1) added line in qsplit to make sure tcoeff
##			      returns correct value.
##                        (2) added extra condition statement to qfactor
##			      so that qfactor(f)=f if f has type symbol
##			      or rational. 
#######################################################################

macro(qsplit = `qfactor/qsplit`, qpfac = `qfactor/qpfac`, qfac=`qfactor/qfac`,qfopexp=`qfactor/qfopexp`, bigqpfac=`qfactor/bigqpfac`);


#######################################################################
# FUNCTION :  qsplit - returns a list [f,tc,ld,td]                      
# CALLING SEQUENCE :  qsplit(f)       
# PARAMETERS :   f   - a polynomial in q
# SYNOPSIS :   
# qsplit(f) returns the list [f,tc,ld,td]
#      where f = tc*(q^ld + ........+()q^td.
#######################################################################
##
 
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);
### WARNING: ldegree(0,x) now returns infinity
    ld := ldegree(f2, q);
### WARNING: degree(0,x) now returns -infinity
    td := degree(f2, q);
    RETURN([f,tc,ld,td]);
end;

#######################################################################
# FUNCTION :  qpfac - factor a q polynom into factors (1-q^i)          
# CALLING SEQUENCE :  qpfac(f,T,s)       
#                      qpfac(f)           
#                      qpfac(f,T)         
#                      qpfac(f,s)         
# PARAMETERS :   f   - a polynomial in q
#                T - upper bound for i (optional)
#                s - string (optional)
# SYNOPSIS :   
# qpfac(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 qpfac
# will return FAIL if f can not be q-factored.
#######################################################################
##
 
#++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
qpfac:=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
### WARNING: degree(0,x) now returns -infinity
	   gd:=degree(g,q);
	   ggc:=coeff(g,q,0);
	   gprod:=prodmake(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;
#######################################################################
# FUNCTION :  qfopexp - factors a SINGLE q-expression and factors 
#                       as much as possible into terms of the form (1-q^i)
# CALLING SEQUENCE :  qfopexp(f,T,s)       
#                      qfopexp(f)           
#                      qfopexp(f,T)         
#                      qfopexp(f,s)         
# PARAMETERS :   f   - a single q-expression (ie. a factor of a factored
#                      expression.
#                T - upper bound for i (optional)
#                s - string (optional)
# SYNOPSIS :   
#   An expression that cant be q-factored is left alone.
#######################################################################
##
 
#++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
 qfopexp:=proc()
local expr,wh,n,whL,cond,s;
    expr:=args[1];
    wh:=whattype(expr);
    n:=nops([op(expr)]);
    whL:=map(whattype,[op(expr)]);
    if n=2 and whL[1]=`+` and whL[2]=integer then
       cond:=1:
    else
       cond:=0;
    fi;
    if wh=`+` or cond=1 then
       s:=subs(q=2,expr);
       if type(s,rational)=true then
          RETURN(qfac(args));
        else
           RETURN(expr);
        fi:
    else
        RETURN(expr);
    fi;
 end;
#######################################################################
# FUNCTION :  bigqpfac - factor a q polynom into factors (1-q^i)          
# CALLING SEQUENCE :  bigqpfac(f,T,s)       
#                      bigqpfac(f)           
#                      bigqpfac(f,T)         
#                      bigqpfac(f,s)         
# PARAMETERS :   f   - a polynomial in q
#                T - upper bound for i (optional)
#                s - string (optional)
# SYNOPSIS :   
# bigqpfac(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 qpfac
# will return FAIL if f can not be q-factored.
# This is a version of qpfac. This version does not fail if it comes
# across a term with other variables beside q - these are simply 
# left alone.
#######################################################################
##
 
#++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
 bigqpfac:=proc()
        local k,expr,na,qtmp,rst,L1,L2,new,conc;
    expr:=args[1]:
###
###01/08/00: add following if .. fi
###
if type(expr,symbol)=true or type(expr,rational)=true then
 RETURN(expr);
else
    na:=nargs;
    if na=1 then
       qtmp:=x->qfopexp(x);
    else
        rst:=seq(args[k],k=2..nargs);
        qtmp:=x->qfopexp(x,rst);
    fi:
###
###01/08/00: add following if .. fi
###
    conc:=0:
    if type(expr,polynom)=true and type(expr,`+`)=true then
	RETURN(qtmp(expr)):
    else	
       L1:=[op(expr)];
       L2:=map(qtmp,L1);
       if type(expr,`^`)=true then
          new:=L2[1]^L2[2];    
          if nops(L2)>2 then
            ERROR(`nops(L2)>2 in bigqfac`);
          fi:
       else
          ##new:=convert(L2,`*`);
          new:=convert(L2,whattype(expr));
       fi:
       RETURN(new);
    fi:
fi:    
### WARNING: `rst` is a lexically scoped local
 end;
#######################################################################
#  2ND MAIN FUNCTION qseries[qfac]
#######################################################################
 
#++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
	        qfac:=proc()
		local f,T;
		f:=args[1]:
	        if nargs=1 then
		   if type(f,ratpoly)=true then
		      if type(f,polynom) then
		         RETURN(qpfac(f));
		      else
	                 RETURN(qpfac(numer(f))/qpfac(denom(f)));
		      fi:
		   else
	              ERROR(`f must be ratpoly`);		
		   fi:
		else
		   T:=args[2]:
		   if type(f,ratpoly)=true then
		      if type(f,polynom) then
		         RETURN(qpfac(f,T));
		      else
	                 RETURN(qpfac(numer(f),T)/qpfac(denom(f),T));
		      fi:
		   else
	              ERROR(`f must be ratpoly`);		
		   fi:
		fi:
	end:
#######################################################################
#  THE MAIN FUNCTION qseries[qfactor]
#######################################################################
		   		
##
 
#++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
	qseries[qfactor]:=proc()
		local f,T,k,n,newargs1,newargs2;
		f:=factor(args[1]):
                ## This part added 10/19/99
### WARNING: the definition of the type `symbol` has changed'; see help page for
### details
                if type(f,symbol)=true or type(f,rational)=true then
                  RETURN(f);
                else

		   if type(f,ratpoly)=true then
		      if type(f,polynom) then
		         RETURN(bigqpfac(f));
		      else
		         n:=nargs;
                         if n=1 then
	                    RETURN(bigqpfac(numer(f))/bigqpfac(denom(f)));
                         else
                            newargs1:=numer(f),seq(args[k],k=2..n);
                            newargs2:=denom(f),seq(args[k],k=2..n);
	                    RETURN(bigqpfac(newargs)/bigqpfac(newargs));
		         fi:
                      fi:
		   else
	              ERROR(`f must be ratpoly`);		
		   fi:
		fi:
	end:
macro(qsplit = qsplit, qpfac = qpfac, qfac=qfac,qfopexp=qfopexp, bigqpfac=bigqpfac);

#++++++++++++++++++++++++++++++++++++++++++++++++++++
     qseries[qs2jaccombo]:=proc()
#       This proc tries to write qseries as a linear combination
#       of jacprods.
        local xx,nn,j,yy,qs,q,T,PP;
        if nargs<>3 and nargs<>4 then
           ERROR("nargs = 3 or 4");
        fi:
        qs:=args[1]:
        q:=args[2]:
        T:=args[3]:
        if nargs=4 then
         PP:=args[4]:
        fi:
        if not(type(qs,`+`)) then
           ERROR("qs not of + type");
        fi:
        xx:=0:
        nn:=nops(qs):
        if nargs=3 then
           for j from 1 to nn do
             yy:=jacprodmake(op(j,qs),q,T):
             xx:=xx+yy:
           od:
        else
           for j from 1 to nn do
             yy:=jacprodmake(op(j,qs),q,T,PP):
             xx:=xx+yy:
           od:
        fi:
       RETURN(xx):
     end:
qseries[quinprod]:=proc(z,q,T)
        local x,i,j,lasti,leftid,x1,x2,rightid;
	# If T = prodid
        # then a product form of the identity is returned.
	# If T = seriesid
        # then a series form of the identity is returned.
	# If T is a positive integer then a series to O(q^T)
	# is returned.
	# Determine output form
	if type(T,posint) then	
	   x:=0:
	   lasti:=trunc( 7/6+1/6*(73+24*T)^(1/2))+1:
	   for i from -lasti to lasti do
	     x:=x + ( (-z)^(-3*i)-(-z)^(3*i+1) )
                    *q^(i*(3*i+1)/2):
	   od:
	   RETURN(x):
	else
           if T=prodid or T=seriesid then
	      leftid:=` `(-z,q)[infinity]
	              *` `(-q/z,q)[infinity]
	              *` `(z^2*q,q^2)[infinity]
	              *` `(1/z^2*q,q^2)[infinity]
	              *` `(q,q)[infinity]:
	      if T=prodid then
	      x1:=` `(q^2/z^3,q^3)[infinity]
	          *` `(q*z^3,q^3)[infinity]
	          *` `(q^3,q^3)[infinity]:
	      x2:=` `(q/z^3,q^3)[infinity]
	          *` `(q^2*z^3,q^3)[infinity]
	          *` `(q^3,q^3)[infinity]:
	      x:= x1+z*x2:
	      RETURN(leftid=x):
              else
	     rightid:=sum(( (-z)^(-3*m)-(-z)^(3*m-1) )
                    *q^(m*(3*m+1)/2),m=-infinity..infinity):
	      RETURN(leftid=rightid):
              fi:
           else
	   fi:
	fi:
end:
	qseries[sift]:=proc(s,q,n,k,T)
	        #
		# sum a_i q^i --> sum a_[ni+k] q^i
	        # 09.04.15: new version of sift can now handle negative exponents
                # See oldsift for oldversion
	 	local y,i,st,lasti,sp,firsti:
                st:=series(s,q,T+5):
        if whattype(st)=series  then
                sp:=convert(st,polynom):
	 	y:=0:
		lasti:=floor((T-k)/n):
                ## NOTE: 09.23.15 lqdegree and lqdegree0 borrowed
                ##       from thetaids package.
                firsti:=ceil((lqdegree(sp)-k)/n):
                ##print("firsti=",firsti,"lasti=",lasti);
	 	for i from firsti to lasti do
	 	  ##y:=y+coeff(s,q,(n*i+k))*q^i:
                  ##fixed 11-16-06
                  ##print("**** i=",i,"y=",y);
	 	  y:=y+coeff(st,q,(n*i+k))*q^i:
	 	od:
		RETURN(y):
	else
	      ERROR(`s must be a series`);
	fi:
	end:
qseries[theta]:=proc(z,q,T)
	local x,i:
	x:=0:
	for i from -T to T do
	  x := x + z^i*q^(i*i):
	od:
	RETURN(x):
end:
qseries[theta2]:=proc(q,T)
        local N:
        N:=trunc(sqrt(T))+2:
        RETURN(qseries[theta](q,q,N)*q^(1/4)):
end:
qseries[theta3]:=proc(q,T)
        local N:
        N:=trunc(sqrt(T))+1:
        RETURN(qseries[theta](1,q,N)):
end:
qseries[theta4]:=proc(q,T)
        local N:
        N:=trunc(sqrt(T))+1:
        RETURN(qseries[theta](-1,q,N)):
end:
qseries[tripleprod]:=proc(z,q,T)
        local x,lasti,i,leftid,rightid;
	# If T = seriesid
        # then a series form of the identity is returned.
	# If T is a positive integer then a series to O(q^T)
	# is returned.
	# Determine output form
	if type(T,posint) then	
   	   x:=0:
   	   lasti:=trunc( sqrt(2*T+1/4)+1/2)+1:
   	   for i from -lasti to lasti do
   	     x:=x + (-1)^i*z^i*q^(i*(i-1)/2):
   	   od:
   	   RETURN(x):
	else
           if T=seriesid then
	      leftid:=` `(z,q)[infinity]
	              *` `(q/z,q)[infinity]
	              *` `(q,q)[infinity]:
	      rightid:=sum((-1)^m*z^m*q^(m*(m-1)/2),   
                   m=-infinity..infinity):
	      RETURN(leftid=rightid):
           else
              ERROR(`Invalid T. T must either be a positive integer or the string seriesid`);
	   fi:
	fi:
end:
qseries[winquist]:=proc(a,b,q,T)
	local x,i,j,lasti:
	x:=0:
	lasti:=trunc( 7/6+1/6*(25+24*T)^(1/2))+1:
	for i from 0 to lasti do
	for j from -lasti to lasti do
	  x:=x + (-1)^(i+j)*( (a^(-3*i)-a^(3*i+3))*(b^(-3*j)-b^(3*j+1))
	                     + (a^(-3*j+1)-a^(3*j+2))*(b^(3*i+2)-b^(-3*i-1)))
                           *q^( 3*i*(i+1)/2 + j*(3*j+1)/2 ):
	od:
	od:
	RETURN(x):
end:
 
#++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
 
 qseries[zqfactor]:=proc()
        local NF,G,numt,numlim,c,FF,TT,LL,ld,tt,zqf,buglim,F,z,q,N;
	if nargs>5 or nargs<4 then
	   ERROR(`zqfactor takes 4 or 5 args.`);
	fi:
	if nargs=5 then 
	   buglim:=args[5]:
	else
	   buglim:=1000:
	fi:
	F:=args[1]:
	z:=args[2]:
	q:=args[3]:
	N:=args[4]:
	#
	#Factor F into a product of terms (1 - a[i,j] z^i q^j)
	#N = largest j 
	#buglim = max number of terms (default is 1000 but can be assigned)
	#
	##
	##
	##EXAMPLE: L1:=tripleprod(-z^2*q^3,q^4,40);
	##         L2:=tripleprod(-z^2*q,q^4,40);
	##         zqfactor(L1-L2/z,q,20);     
	##
        G:=1:
        NF:=F:
        numt:=0:
        numlim:=buglim:
        while NF<>1 and numt<=numlim do
             c:=coeff(NF,q,0):
             if c<>1 then
                 NF:=NF/c:
                 G:=G/c:
 	          numt:=numt+1:
             fi:
             FF:=normal(series(NF,q,N+1)):
 	      FF:=convert(FF,polynom):
             TT:=expand(tcoeff(FF-1,q)):
             LL:=[op(TT)]:
             ld:=ldegree(FF-1,q):
             for tt in LL do
                if type(tt,integer) and tt<>(-1) then
                   zqf:=(1-q^ld)^(-tt):
                else
                   zqf:=1+tt*q^ld:
                fi:
                numt:=numt+1:
                G:=G/zqf: NF:=NF/zqf:
                NF:=normal(series(NF,q,N+1));
                 NF:=convert(NF,polynom):
             od;
          od;
          RETURN(1/G);
 end;
#++++++++++++++++++++++++++++++++++++++++++++++++++++
printf("TABLE TYPE qseries = %a\n",type(qseries,table));
savelib( qseries, 
`jacprodmake/jaccheck`
, 
`jacprodmake/periodfind` 
, 
`jacprodmake/periodfind2` 
, 
`jacprodmake/jacmake`
, 
`jac2prod/jac`
, 
`jac2J/jacJ`
, 
`jac2series/tripleprod2`
, 
`jac2series/qjac`
, 
`qfactor/qsplit`
, 
`qfactor/qpfac`
, 
`qfactor/qfac`
, 
`qfactor/qfopexp`
, 
`qfactor/bigqpfac`
,
"c:\\cygwin64/home/fgarv/maple/mylib/qseries.mla");
printf("END qseries package\n");
## mylib  above must be changed to the name of the directory
## in which you want stuff stored