#######################################################
## make-win-package thetaids
## Sat Jul 21 12:00:05 EDT 2018
#######################################################
#######################################################

printf("BEGIN thetaids package\n");
printf("THIS VERSION DATED Sat Jul 21 12:00:05 EDT 2018 \n");
thetaids:=table();
thetaids[Acmake]:=proc()
local S1,S,sp,s,c,N;
    # c is a positive divisor of N
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf(" Acmake(c,N)                                             \n");
     printf("    Returns the set A[c] where c is a positive           \n");
     printf("    divisor of N.                                        \n");
     printf("---------------------------------------------------------\n");
  elif nargs = 2 then
     #check input:
      c:=args[1]: N:=args[2]:
      if map(x->type(x,posint),[c,N])<>[true,true] then
         ERROR(`c,N must be positive integers`);
      fi:
      if modp(N,c)<>0 then
         ERROR(`c must be a positive divisor of N`);
      fi:
    if (N/c)=1 or (N/c)=2 then
       S1:=sset(c):
    else
       S1:=phiset(c):
    fi:
    S:={}:
    for sp in S1 do
        s:=srep(sp,c,N):
        S:=S union {s}:
    od:
    RETURN(S):
  else
     ERROR(`nargs should be 0 or 2`);
  fi:
end:
thetaids[Bord]:=proc()
	local enc,nr,t,g,n,r,a,c;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("Bord(n,r,a,c)                                                    \n");
     printf("   Calculates ord(f[n,r],a/c)\n");                      
     printf("   where f[n,r] = q^( (n-2r)^2/(8n) )*(q^r,q^(n-r),q^n;q^n)oo\n");
     printf("   as defined on p.277 of Biagioli. See Lemma 3.2 of Biagioli\n");
     printf("   Assumptions: (a,c)=1 and n does not divide r.\n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 4 then
        n:=args[1]: r:=args[2]: a:=args[3]: c:=args[4]:
	if gcd(a,c)<>1 then
	   ERROR(` gcd(a,c) must be 1.`);
	fi:	
	enc:=gcd(n,c):
	nr:=modp(r,n):
	if nr=0 then
	   ERROR(` n can not divide r.`);
	fi:
	t:=floor(a*r/enc):
	RETURN(  enc^2/2/n*(a*r/enc - t - 1/2)^2 ):
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 4`);
  fi:
end:
thetaids[cuspequiv1]:=proc()
local checkgcd,n,mm1,mm2,a1,c1,a2,c2,N,check1,epc;
    ##global check1,epc;
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("cuspequiv1(a1,c1,a2,c2,N)                                \n");
     printf("   Determines whether the cusps a1/c1 and a2/c2 are      \n");
     printf("   Gamma[1](N)-equivalent.                               \n");
     printf("---------------------------------------------------------\n");
  elif nargs = 5 then
     #check input:
      a1:=args[1]: c1:=args[2]: a2:=args[3]: c2:=args[4]: N:=args[5]:
      if map(x->type(x,integer),[a1,c1,a2,c2,N])<>[true,true,true,true,true] then
         ERROR(`args must be integers`);
      fi:
      if not type(N,posint) then
         ERROR(`N must be a positive integer.`);
      fi:
      if igcd(a1,c1)<>1 or igcd(a2,c2)<>1 then
         ERROR(`Need gcd(a1,c1)=gcd(a2,c2)=1`);
      fi:
    #global check1,epc;
    # Determine whether the cusp a/c -s Gamma1(N)-equivalent
    # to a1/c1
   check1:=0:
## checkgcd:=igcd(a1,c1)*igcd(a2,c2):
## if checkgcd<>1 then
##    lprint("gcd not in cuspequiv1 with args",a1,c1,a2,c2,N);
##    error "gcd not in cuspequiv1" ;
## fi:
    if modp(c1-c2,N)<>0 and modp(c1+c2,N)<>0 then
       RETURN(false):
   else
      if modp(c1-c2,N)=0 then
         epc:=1:
      else
         epc:=-1:
      fi:
      #
      #NOTE: c1==c2 and c1==-c2 mod N implies 2c1==0 mod N.
      #
      if modp(2*c1,N)=0 then
          for n from 0 to N-1 do
             mm1:=modp(a2-a1-n*c1,N):
             mm2:=modp(-a2-a1-n*c1,N):
             if mm1=0 or mm2=0 then     
               check1:=1:
               RETURN(true): 
             fi:
          od:
       else
          for n from 0 to N-1 do
             if modp(epc*a2-a1-n*c1,N)=0 then
               check1:=1:
               RETURN(true): 
             fi:
          od:
      fi:
   fi:
   if check1=0 then
       RETURN(false):
   fi:
  else
     ERROR(`invalid input type`);
  fi:
end:
thetaids[cuspmake1]:=proc()
local divN,xx,cusps,c,SSET,ASET,sci,acj,x,y,NXY,nx1,ny1,N;
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("cuspmake1(N)                                             \n");
     printf("   Returns a set of inequivalent cusps                   \n");
     printf("   for Gamma[1](N).                                      \n");
     printf("   Each cusp a/c in the list is represented by [a,c]     \n");
     printf("   so that oo is represented by [1,0]                    \n");
     printf("---------------------------------------------------------\n");
  elif nargs = 1 and type(args[1],posint) then
    N:=args[1]:
    divN:=numtheory[divisors](N):
    xx:=0:
    cusps:={}:
    for c in divN do
       #make Sc set and Ac set
        SSET:=Scmake(c,N):
        ASET:=Acmake(c,N):
        for sci in SSET do
        for acj in ASET do
           x:=modp(c*sci,N):
           y:=modp(acj,N):
           NXY:=newxy(x,y,N):
           nx1:=NXY[1]:
           ny1:=NXY[2]:
           if nx1<>0 then
              ny1:=modp(ny1,nx1):
           fi:  
           ## cusp = y/x
           ##cusps:=cusps union {[x,y,igcd(x,y)]}:
           cusps:=cusps union {[ny1,nx1]}:
        od:
        od:
    od:
    RETURN(cusps):
  else
     ERROR(`invalid input type`);
  fi:
end:
thetaids[CUSPSANDWIDMAKE1]:=proc()
     local cuspsN,widsN,CUSPSN,N:
     global CUSPMAKEERROR:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("CUSPSANDWIDMAKE(N) \n");
     printf("   N is a positive integer\n");                      
     printf("   Returns a set of inequivalent cusps for Gamma[1](N),\n");  
     printf("   and corresponding widths.\n");                      
     printf("   Output has the form [CUSPLIST,WIDTHLIST].\n");     
     printf("-------------------------------------------------------------\n");
  elif nargs = 1 and type(args[1],posint) then
     N:=args[1]:
     cuspsN:=cuspmake1(N):
     if member([1,0],cuspsN) then
        cuspsN:=cuspsN minus {[1,0]}:
        cuspsN:=convert(cuspsN,list):
        widsN:=map(x->cuspwid1(x[1],x[2],N),cuspsN);
        widsN:=[1,op(widsN)];
        CUSPSN:=map(x->x[1]/x[2],cuspsN);
        CUSPSN:=[oo,op(CUSPSN)];
        RETURN([CUSPSN,widsN]):
     else
        CUSPMAKEERROR:=1:
        ERROR("[1,0] missing in cusp set");
     fi:
  else
     ERROR(`invalid input type`);
  fi:
end:
thetaids[cuspsetinequiv1]:=proc()
local j,c,cm,c1,checkinequiv,nn,choose2,kk,k,f1,f2,cusp1,cusp2,checkequiv,cset,N;
     global DSET,equivpairs:
     # Determine whether the cusps in the given set
     # are inequivalent mod Gamma1(N).
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("cuspsetinequiv1(cset,N)                                  \n");
     printf("   Determines whether the cusps in in cset are           \n");
     printf("   Gamma[1](N)-inequivalent.                             \n");
     printf("   This proc returns true or false. If equivalent pairs  \n");
     printf("   are found they are stored in the global variable      \n");
     printf("   equivpairs.                                           \n");
     printf("---------------------------------------------------------\n");
  elif nargs = 2 and type(args[2],posint) and type(args[1],set) then
     #
     # First we breakup the cset according the denom mod N
     cset:=args[1]:  N:=args[2]:
     equivpairs:={}:
     DSET:=array(0..trunc(N/2)):     
     # Initialize DSET
     for j from 0 to trunc(N/2) do
       DSET[j]:={}:
     od:
     for j from 1 to nops(cset) do
       c:=cset[j][2]:
       cm:=modp(c,N):
       if 2*cm>N then
         c1:=N-cm:
       else
         c1:=cm:
       fi:
       DSET[c1]:=DSET[c1] union {cset[j]}:
    od:
    checkinequiv:=1:
    # now check pairs in each DSET for equivalence
    for j from 0 to trunc(N/2) do
       nn:=nops(DSET[j]):
       lprint("DSET",j,"nops=",nn);
       if nn>1 then
         choose2:=combinat[choose](nn,2):
         kk:=nops(choose2):
         for k from 1 to kk do
             f1:=choose2[k][1]:
             f2:=choose2[k][2]:
             # f1,f2 is a pair from the set {1,2,...,nn}
             cusp1:=DSET[j][f1]:
             cusp2:=DSET[j][f2]:
             checkequiv:=cuspequiv1(cusp1[1],cusp1[2],cusp2[1],cusp2[2],N):
             if checkequiv then
                lprint("WARNING: ",cusp1,cusp2,"are equivalent");
                equivpairs:=equivpairs union {[cusp1,cusp2]}:
                checkinequiv:=0:
             fi:
         od:
       fi:
    od:
    if checkinequiv=1 then
      printf("All cusps in the set are inequivalent.");
      RETURN(true):
    else
      lprint("WARNING: Set contains equivalent cusps.");
      lprint("SEE global equivpairs.");
      RETURN(false):
    fi:
  else
     ERROR(`invalid input type`);
  fi:
end:
thetaids[cuspwid1]:=proc()
   local a,c,N:
   # cusp width of a/c in Gamma1(N)
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("cuspwid1(a,c,N)                                          \n");
     printf("   Returns the width of the cusp a/c                     \n");
     printf("   for the group Gamma[1](N).                            \n");
     printf("---------------------------------------------------------\n");
  elif nargs = 3 then
     #check input:
      a:=args[1]: c:=args[2]: N:=args[3]:
      if map(x->type(x,integer),[a,c,N])<>[true,true,true] then
         ERROR(`a,c,N must be integers`);
      fi:
      if N<=0 then
         ERROR(`N must be a positive integer`);         
      fi:
   if N=4 and c=2 then
      RETURN(1):
   else
      RETURN(N/igcd(c,N)):
   fi:
  else
     ERROR(`nargs should be 0 or 3`);
  fi:
end:
thetaids[DivCheck]:=proc()
local divs,chk,nL,j,v,n,L,N;
  #check that all [n,m,r] in L satisfy n divisor of N
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("DivCheck(L,N)                                               \n");
     printf("   L is a getalist; i.e. L=[[b1,a1,c1],....]\n");
     printf("   Checks that n is a divisor of N for each [n,m,r] in L.\n");
     printf("   Returns 1 if it checks.\n");                      
     printf("-------------------------------------------------------------\n");
  elif nargs = 2 and type(args[2],posint) then
  L:=args[1]: N:=args[2]:
  divs:=numtheory[divisors](N):
  chk:=1:
  nL:=nops(L):
  for j from 1 to nL do
      v:=L[j]:
      n:=v[1]:
      if member(n,divs) then
         ##OK
      else
	 chk:=0:
      fi:
  od:
  RETURN(chk):
  else
     ERROR(`invalid input type`);
  fi:
end:
thetaids[eeta]:=proc(b)
        eta(b*tau):
end:
thetaids[egeta]:=proc(d,g)
        eta[d,g](tau):
end:
thetaids[ejac]:=proc(a,b,c)
# 07/19/13: added power of q
local qpow1:
if a=0 then
	EETA(b)/q^(b/24):
else
	qpow1:=getacuspord(b,a,1,0):
	GETA(b,a)*EETA(b)/q^(qpow1+b/24):
fi:
end:
thetaids[eta2geta]:=proc()
  local etaexp,EE:                  
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("eta2geta(etaexp)                                       \n");
     printf("   This function is a notation conversion.                  \n");
     printf("   It does the following                                   \n");
     printf("        ETAnm(d,g) ->  eta[d,g](tau)                       \n");
     printf("        ETAn(b) ->  eta(b*tau)                       \n");
     printf("-------------------------------------------------------------\n");
#############################################################################################
## EXAMPLE:
## > EP:=ETAnm(5,1)/ETAnm(5,2):
## > EP1:=ETAn(1):
## > eta2jac(EP);
##                                                          1/5
##                                     JAC(1, 5, infinity) q
##                                     ------------------------
##                                       JAC(2, 5, infinity)
## 
## > eta2jac(EP1);
##                                    (1/24)
##                                   q       JAC(0, 1, infinity)
## 
## > eta2geta(EP);
##                                          eta[5, 1](tau)
##                                          --------------
##                                          eta[5, 2](tau)
## 
## > eta2geta(EP1);
##                                             eta(tau)
## 
#############################################################################################
  elif nargs = 1 then
     etaexp:=args[1]:
     EE:=eval(subs({ETAnm=_etanm,ETAn=_etan},etaexp)):
     RETURN(EE):
  else
     ERROR(`nargs should be 0 or 1`);
  fi:
end:
thetaids[eta2jac]:=proc()
  local GP:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("eta2jac(getaexp)                                            \n");
     printf("   getaexp is expression involving ETAnm(d,g) and ETAn(b).\n");
     printf("   ETAnm(d.g)=eta[d,g](tau) amd ETAn(b)=eta(b*tau). \n");
     printf("   This function converts getaexp into JAC-prod form;ie.  \n");
     printf("   in terms of JAC(a,b,infinity).                         \n");
     printf("-------------------------------------------------------------\n");
#############################################################################################
## EXAMPLE:
## > EP:=ETAnm(5,1)/ETAnm(5,2):
## > EP1:=ETAn(1):
## > eta2jac(EP);
##                                                          1/5
##                                     JAC(1, 5, infinity) q
##                                     ------------------------
##                                       JAC(2, 5, infinity)
## 
## > eta2jac(EP1);
##                                    (1/24)
##                                   q       JAC(0, 1, infinity)
## 
## > eta2geta(EP);
##                                          eta[5, 1](tau)
##                                          --------------
##                                          eta[5, 2](tau)
## 
## > eta2geta(EP1);
##                                             eta(tau)
## 
#############################################################################################
  elif nargs = 1 then
     GP:=args[1]:
     eval(subs({ETAnm = `_ETAnm`,ETAn = `_ETAn`}, GP)):
  else
     ERROR(`nargs should be 0 or 1`);
  fi:
end:
thetaids[_etan]:=n->eta(n*tau);                     
thetaids[_ETAn]:=n->q^(n/24)*JAC(0,n,infinity);
thetaids[_etanm]:=(n,m)->eta[n,m](tau);
thetaids[_ETAnm]:=(n,m)->JAC(m,n,infinity)*q^(QP2(m/n)*n/2)/JAC(0,n,infinity);
thetaids[fpart]:=t->t-trunc(t):
thetaids[Gamma1ModFunc]:=proc() ##checks whether geta prod is a modfunc on
## EXAMPLE
## > L2;
##                                  [[25, 5, -1], [25, 10, 1]]
## > Gamma1ModFunc(L2,25);
##                                 "All n are divisors of ", 25
##                                          "val0=", 0
##                                       "which is even."
##                                        "valinf=", -2
##                                       "which is even."
##                            "It IS a modfunc on Gamma1(", 25, ")"
##                                              1
local mfunc,val0,valinf,L,N;
global xprint:
        if not(type(xprint,boolean)) then
          xprint:=false:
        fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("Gamma1ModFunc(L,N)                                               \n");
     printf("   L is a getalist; i.e. L=[[b1,a1,c1],....]\n");
     printf("   and assumes each a1, a2, etc is not zero.\n");
     printf("   Let G be the corresponding generalized eta-function.\n");
     printf("   Checks whether G is a modular function of Gamma[1](N).\n");
     printf("   Returns 1 if it is a modular function otherwise 0.\n");
     printf("   It also prints diagnostics.\n");                               
     printf("-------------------------------------------------------------\n");
  elif nargs = 2 and type(args[2],posint) then
  L:=args[1]: N:=args[2]:
  if xprint then
     printf("* starting Gamma1ModFunc with L=%a and N=%a \n",L,N);
  fi:
  mfunc:=1:
  if DivCheck(L,N)=1 then
     if xprint then
       printf("All n are divisors of %a\n",N);
     fi:
  else
     if xprint then
       printf("All n are NOT divisors of %a\n",N);
       printf("Hence NOT a modfunc on Gamma1(%a)\n",N); 
     fi:
     mfunc:=0:
  fi:
  val0:=v0(L,N):
  if xprint then
    printf("val0=%a\n",val0);
  fi:
  if type(val0,integer) then
    #OK
    if modp(val0,2)=0 then
      if xprint then
        printf("which is even.\n");
      fi:
    else
      if xprint then
        printf("which is odd. Hence NOT a modfunc on Gamma1(%a)\n",N);
      fi:
      mfunc:=0
    fi:
  else
    if xprint then
      printf("which is NOT an integer. Hence NOT a modfunc on Gamma1(%a)\n",N);
    fi:
    mfunc:=0:
  fi:
  valinf:=vinf(L,N):
  if xprint then
     printf("valinf=%a\n",valinf);
  fi:
  if type(valinf,integer) then
    #OK
    if modp(valinf,2)=0 then
      if xprint then
        printf("which is even.\n");
      fi:
    else
      if xprint then
        printf("which is odd. Hence NOT a modfunc on Gamma1(%a)\n",N);
      fi:
      mfunc:=0:
    fi:
  else
    if xprint then
      printf("which is NOT an integer. Hence NOT a modfunc on Gamma1(%a)\n",N);
    fi:
    mfunc:=0:
  fi:
  if mfunc=1 then
     if xprint then
       printf("It IS a modfunc on Gamma1(%a)\n",N);
     fi:
  fi:
  RETURN(mfunc):
  else
     ERROR(`invalid input type`);
  fi:
end:
thetaids[getacuspord]:=proc()
local x1,x2,n,r,a,c,nr;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("getacuspord(n,r,a,c)                                         \n");
     printf("   Calculates order of generalized eta-product         \n");
     printf("   eta[n,r](tau) at the cusp a/c.\n");
     printf("   Assumptions: (a,c)=1 and n does not divide r.\n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 4 then
        n:=args[1]: r:=args[2]: a:=args[3]: c:=args[4]:
	if gcd(a,c)<>1 then
	   ERROR(` gcd(a,c) must be 1.`);
	fi:	
	nr:=modp(r,n):
	if nr=0 then
	   ERROR(` n can not divide r.`);
	fi:
   if r=0 then
        x1:=igcd(n,c)^2/n/24:
	x2:=Bord(3*n,n,a,c):
	if x1<>x2 then
	  ERROR(`x1 not x2`):
	fi:
	x1:
   else
	Bord(n,r,a,c)-Bord(3*n,n,a,c):
   fi:
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 4`);
  fi:
end:
thetaids[getalist2jacprod]:=proc()
global glist,JP,gterm: 
## EXAMPLE: Rogers-Ramanujan Continued Fraction
## M := 100: a[ M ] := 1+q; 
## for n from M-1 by -1 to 1 do 
## a[ n ] := 1 + q^n/a[ n+1 ]:
## od:
## > series(a[1],q,100):
## > jacprodmake(%,q,50);
##                                     JAC(2, 5, infinity)
##                                     -------------------
##                                     JAC(1, 5, infinity)
## > jac2eprod(%);
##                                         GETA(5, 2)
##                                         ----------
##                                         GETA(5, 1)
## > GETAP2getalist(%);
##                                   [[5, 1, -1], [5, 2, 1]]
		local L1,L1n,GP,i,p,t,chk,L1i,s,r:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("getalistjacprod(getalist)                                   \n");
     printf("   Each item [b,a,c] in getalist corresponds to \n");
     printf("   GETA(b,a)^c.  The list is converted into a product\n");
     printf("   of jacprods.             \n");
     printf("-------------------------------------------------------------\n");
  elif nargs = 1 then
     glist:=args[1]:
     JP:=mul(gterm2jacprod(gterm), gterm in glist):
     RETURN(JP):
  fi:
end:
thetaids[GETAP2getalist]:=proc()
global geprod:
## EXAMPLE: Rogers-Ramanujan Continued Fraction
## M := 100: a[ M ] := 1+q; 
## for n from M-1 by -1 to 1 do 
## a[ n ] := 1 + q^n/a[ n+1 ]:
## od:
## > series(a[1],q,100):
## > jacprodmake(%,q,50);
##                                     JAC(2, 5, infinity)
##                                     -------------------
##                                     JAC(1, 5, infinity)
## > jac2eprod(%);
##                                         GETA(5, 2)
##                                         ----------
##                                         GETA(5, 1)
## > GETAP2getalist(%);
##                                   [[5, 1, -1], [5, 2, 1]]
		local L1,L1n,GP,i,p,t,chk,L1i,s,r:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("GETAP2getalist(geprod)                                   \n");
     printf("   geprod is a quotient of generalized eta-functions encoded\n");
     printf("   in terms of GETA(b,a).                                   \n");
     printf("   The quotient is converted to a list.\n");
     printf("   Each term GETA(b,a)^c in the quotient is converted to an\n");
     printf("   item [b,a,c] in the list.\n");
     printf("-------------------------------------------------------------\n");
  elif nargs = 1 then
     geprod:=args[1]:
	        if whattype(geprod)=`^` then
		   L1:=[geprod]:
	     	else
		   if whattype(geprod)=function then
			RETURN([op(geprod),1]);
		   else
		        L1:=[op(geprod)]:
	           fi:
		fi:
		L1n:=nops(L1):
		GP:=[]:   
		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:
	            chk:=simplify(L1i-GETA(s)^r):
	            if chk<>0 then
		      error "ERROR in GETAP2getalist chk<>0";
	            fi:  
		    GP:=[op(GP),[s,r]]:
		od:
		RETURN(GP);
  else
     ERROR(`invalid input type`);
  fi:
	end:	
 
thetaids[getaprodcuspord]:=proc()
local tot,nc,nL,ORDS,j,a,c,kco,k,gg,n,r,mm,L,cusps,wids,xprint,cusp;
  #05/07/10: added xprint
  if not(type(xprint,boolean)) then
     xprint:=false:
  fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("getaprodcuspord(L,cusp)                                     \n");
     printf("   Let G be a generalized-etaproduct corresponding to the    \n");
     printf("   getalist L. This proc calculates invariant order ord(G,z)\n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 2 then
        L:=args[1]: cusp:=args[2]:  
        nL:=nops(L):
	if cusp=oo then
	   a:=1:
	   c:=0:
	else
	   a:=numer(cusp):
	   c:=denom(cusp):
	fi:
	kco:=0:
	##print("cusp=",a,c);
	for k from 1 to nL do
	  gg:=L[k]:
          n:=gg[1]:
	  r:=gg[2]:
          mm:=gg[3]:
	  kco:=	kco+getacuspord(n,r,a,c)*mm:
	od:
        ##print(a,c,kco);
        RETURN(kco);
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 3`);
  fi:
end:
#++++++++++++++++++++++++++++++++++++++++++++++++++++
thetaids[getaprodcuspORDS]:=proc()
local tot,nc,nL,ORDS,j,a,c,kco,k,gg,n,r,mm,L,cusps,wids;
global toterror,xprint:
  #05/07/10: added xprint
  if not(type(xprint,boolean)) then
     xprint:=false:
  fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("getaprodcuspORDS(L,cusps,wids)                               \n");
     printf("   Let G be a generalized-etaproduct corresponding to the    \n");
     printf("   getalist L. This proc calculates ORD(G,z) with respect to \n");
     printf("   Gamma[1](N) for each cusp z in cusps. Here cusps is a list\n");
     printf("   of inequivalent cusps of Gamma[1](N) and wids is the list \n");
     printf("   of their corresponding widths. The cusp at infinity is    \n");
     printf("   repesented by oo. The total ORD should be 0.              \n");
     printf("   Global var toterror = total ORD (for error-checking).     \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 3 then
        L:=args[1]: cusps:=args[2]: wids:=args[3]: 
  tot:=0:
  nc:=nops(cusps):
  nL:=nops(L):
  ORDS:=[]:	
  for j from 1 to nc do
	if cusps[j]=oo then
	   a:=1:
	   c:=0:
	else
	   a:=numer(cusps[j]):
	   c:=denom(cusps[j]):
	fi:
	kco:=0:
	##print("cusp=",a,c);
	for k from 1 to nL do
	  gg:=L[k]:
          n:=gg[1]:
	  r:=gg[2]:
          mm:=gg[3]:
	  kco:=	kco+getacuspord(n,r,a,c)*mm*wids[j]:
	od:
        ##print(a,c,kco);
	ORDS:=[op(ORDS),kco]:
	tot:=tot+kco:
  od:
  if xprint then
    printf("Cusp ords: \n");
    print([seq([cusps[j],ORDS[j]],j=1..nc)]);
    printf("TOTAL ORD = %a\n",tot);
  fi:
  toterror:=tot:
RETURN(ORDS):
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 3`);
  fi:
end:
thetaids[gterm2jacprod]:=proc()
local gterm,a,b,c,JP,qexp:   
## EXAMPLE: Rogers-Ramanujan Continued Fraction
## M := 100: a[ M ] := 1+q; 
## for n from M-1 by -1 to 1 do 
## a[ n ] := 1 + q^n/a[ n+1 ]:
## od:
## > series(a[1],q,100):
## > jacprodmake(%,q,50);
##                                     JAC(2, 5, infinity)
##                                     -------------------
##                                     JAC(1, 5, infinity)
## > jac2eprod(%);
##                                         GETA(5, 2)
##                                         ----------
##                                         GETA(5, 1)
## > GETAP2getalist(%);
##                                   [[5, 1, -1], [5, 2, 1]]
		local L1,L1n,GP,i,p,t,chk,L1i,s,r:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("gterm2jacprod(gterm)                                   \n");
     printf("   gterm = [b,a,c] corresponds to \n");
     printf("   GETA(b,a)^c.  This is converted to a jacprod.    \n");
     printf("-------------------------------------------------------------\n");
  elif nargs = 1 then
     gterm:=args[1]:
	        b:=gterm[1]: a:=gterm[2]: c:=gterm[3]:
                qexp:=b*QP2(a/b)/2:
                if a=0 then
		   JP:=JAC(0,b,infinity)^(2*c)*q^(b*c/12):
                else
                   JP:=JAC(a,b,infinity)^c*q^(qexp*c)/JAC(0,b,infinity)^c:
	     	fi:
     RETURN(JP):
  fi:
end:
thetaids[jac2eprod]:=proc()
  local JJ:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("jac2eprod(JJ)                                            \n");
     printf("   JJ is a quotient of theta functions encoded in terms of \n");
     printf("   JAC(a,b,infinity). This proc converts this quotient into\n");
     printf("   a quotient of eta-products and generalized eta-products\n");
     printf("   The output is terms of EETA(b) and GETA(b,a).          \n");
     printf("   EETA(b) corresponds to eta(b*tau), and                 \n");
     printf("   GETA(b,a) corresponds to the generalized eta-function \n");
     printf("   eta[b,a](tau).\n");
     printf("-------------------------------------------------------------\n");
# EXAMPLE:
# > with(thetaids):
# > with(qseries):
# > X:=add( (-1)^n*q^(n*(5*n-1)/2),n=-20..20):
# > J:=jacprodmake(X,q,100);
#                            J := JAC(2, 5, infinity)
# 
# > jac2eprod(J);
#                               GETA(5, 2) EETA(5)
#                               ------------------
#                                      1/40
#                                     q
# 
  elif nargs = 1 then
     JJ:=args[1]:
     simplify(eval(subs(JAC = `ejac`, JJ)),symbolic):
  else
     ERROR(`nargs should be 0 or 1`);
  fi:
end:
thetaids[jac2getacombo]:=proc()
  local JJ,xx,nn,n,JJ2;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("jac2getacombo(JJ)                                            \n");
     printf("   This function converts a combination of jacprods into a  \n");
     printf("   combination of getaprods. \n");
     printf("-------------------------------------------------------------\n");
  elif nargs = 1 then
     JJ:=args[1]:
     if type(JJ,`+`) then
         xx:=0:
         nn:=nops(JJ):
         for n from 1 to nn do
             JJ2:=jac2getaprod(op(n,JJ)):
             xx := xx + JJ2:
         od:
     else
         xx:=jac2getaprod(JJ):
     fi:
     RETURN(xx):
  else
     ERROR(`nargs should be 0 or 1`);
  fi:
end:
thetaids[jac2getaprod]:=proc()
  local JJ,JJ1,JJ2,JJ3:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("jac2getaprod(JJ)                                            \n");
     printf("   This function is a version of jac2eprod.                 \n");
     printf("   JJ is a quotient of theta functions encoded in terms of \n");
     printf("   JAC(a,b,infinity). This proc converts this quotient into\n");
     printf("   a quotient of eta-products and generalized eta-products\n");
     printf("   using the notation eta(b*tau) and eta[b,a](tau).       \n");
     printf("   Instead of eta(b*tau) and eta[b,a](tau) jac2eprod uses \n");
     printf("   the notation EETA(b) and GETA(b,a).\n");
     printf("-------------------------------------------------------------\n");
#EXAMPLE:
#> X:=add((-1)^n*q^(n*(5*n-1)/2),n=-10..10);
#          255    207    164    126    93    65    42    24    11    3    2    9
#X := 1 + q    - q    + q    - q    + q   - q   + q   - q   + q   - q  - q  + q
#
#       21    38    60    87    119    156    198    245
#    - q   + q   - q   + q   - q    + q    - q    + q
#> JJ:=jacprodmake(X,q,100);
#                           JJ := JAC(2, 5, infinity)
#> jac2getaprod(JJ);
#                               GETA(5, 2) EETA(5)
  elif nargs = 1 then
     JJ:=args[1]:
     JJ1:=jac2eprod(JJ):
     JJ2:=eval(subs({GETA=egeta,EETA=eeta},JJ1)):
     JJ3:= simplify(JJ2,symbolic):
     RETURN(JJ3):
  else
     ERROR(`nargs should be 0 or 1`);
  fi:
end:
thetaids[jacbase]:=proc()
local jacterm,cc,jbases,jbase,jprod; 
global jacterm1, jacterm2,geprod,jaclist:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("jacbase(jacterm)                                        \n");
    printf("   jacterm is  weight zero quotient of theta functions written\n");
     printf("   in terms of JAC(a,b,infinity). It may also involve a      \n");
     printf("   power of q. It returns the base if jacprodmake is used    \n");
     printf("   correctly.                                  \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 1 then
        jacterm:=args[1]:  
        cc:=eval(subs({JAC=JACCOF,q=1},jacterm)):
        jacterm1:=jacterm/cc:
        jacterm2:=jacterm1/q^lqdegree0(jacterm1):
        if jacterm2=1 then
           RETURN(1):
        else
           jprod:=jacterm2:                 
           jaclist:=JACP2jaclist(jprod):
           jbases:=seq(jaclist[j][2],j=1..nops(jaclist)):
           jbase:=ilcm(jbases):
           RETURN(jbase):
        fi:
  fi:
end:
thetaids[JACCOF]:=(a,b,c)->1:
thetaids[jacnormalid]:=proc()
     local K,DG,N,j,DG0,N0,jcombo;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("jacnormalid(jcombo)                                          \n");
     printf("   Renormalizes a sum of jacprods by dividing by the term \n");
     printf("   with the lowest power of q. \n");
     printf("-------------------------------------------------------------\n");
  elif nargs = 1 and type(args[1],`+`) then
      jcombo:=args[1]:
        K:=1:
        DG:=qjacdegree(op(1,jcombo)):
        N:=nops(op(1,jcombo)):
        ## find term with smallest degree 
        ## and smallest nops if there is a tie
        for j from 2 to nops(jcombo) do
           DG0:=qjacdegree(op(j,jcombo)):
           N0:=nops(op(j,jcombo)):
           if DG0<DG then
              K:=j:
              DG:=DG0:
              N:=N0:
           fi:
           if DG0=DG and N0<N then
              K:=j:
              N:=N0:
           fi:
        od:
        RETURN(expand(jcombo/op(K,jcombo)));
     else
        ERROR("jcombo not of type + or nargs not 0 or 1");
     fi:
end:
thetaids[JACP2jaclist]:=proc()
## EXAMPLE: Rojacrs-Ramanujan Continued Fraction
## M := 100: a[ M ] := 1+q; 
## for n from M-1 by -1 to 1 do 
## a[ n ] := 1 + q^n/a[ n+1 ]:
## od:
## > series(a[1],q,100):
## > jacprodmake(%,q,50);
##                                     JAC(2, 5, infinity)
##                                     -------------------
##                                     JAC(1, 5, infinity)
## > jac2eprod(%);
##                                         GETA(5, 2)
##                                         ----------
##                                         GETA(5, 1)
## > GETAP2jaclist(%);
##                                   [[5, 1, -1], [5, 2, 1]]
		local L1,L1n,JP,i,p,t,chk,L1i,s,r:
global jprod:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("JACP2jaclist(jprod)                                   \n");
     printf("   jprod is a quotient of theta-functions encoded\n");
     printf("   in terms of JAC(a,b,infinity).                            \n");
     printf("   The quotient is converted to a list.\n");
     printf("   Each term JAC(a,b,infinity)^c in the quotient is converted to an\n");
     printf("   item [a,b,infinity,c] in the list.\n");
     printf("-------------------------------------------------------------\n");
  elif nargs = 1 then
     jprod:=args[1]:
	        if whattype(jprod)=`^` then
		   L1:=[jprod]:
	     	else
		   if whattype(jprod)=function then
			RETURN([[op(jprod),1]]);
		   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:
	            chk:=simplify(L1i-JAC(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:	
thetaids[jcombobase]:=proc()
  local jacprod,JP,BASESET,jcombo,S,j,jterm,jb:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("jcombobase(jcombo)                                   \n");
     printf("   jcombo is a sum of jacprods. \n");
     printf("   jacbase of each term is computed and the common lcm     \n");
     printf("   is returned.             \n");
     printf("-------------------------------------------------------------\n");
## EXAMPLE:
## > with(qseries):
## > jcomb:= q^(1/5)*JAC(1,5,infinity)/JAC(2,5,infinity) + 
##           q^(2/5)*JAC(2,10,infinity)/JAC(4,10,infinity);
##                    1/5                        2/5
##                   q    JAC(1, 5, infinity)   q    JAC(2, 10, infinity)
##          jcomb := ------------------------ + -------------------------
##                     JAC(2, 5, infinity)        JAC(4, 10, infinity)
## 
## > jac2getacombo(jcomb);                                                        
##                        eta[1, 5](tau)   eta[2, 10](tau)
##                        -------------- + ---------------
##                        eta[2, 5](tau)   eta[4, 10](tau)
## 
## > jcombobase(jcomb); 
##                                       10
## 
  elif nargs = 1 then
     jcombo:=args[1]:
     if type(jcombo,`+`) then
        S:={}:
        for j from 1 to nops(jcombo) do
           jterm:=rmcofjac(op(j,jcombo)):
           jb:=jacbase(jterm):
           S := S union {jb}:
        od:
        RETURN(ilcm(op(S))):
     else
        RETURN(jacbase(jcombo)):
     fi:
  else
     ERROR(`invalid input type`);
  fi:
end:	
thetaids[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:
thetaids[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:
thetaids[mintotORDS]:=proc()
# 09/06/09: version of mintotORD2 (mintotORD2 does not work if L in an array)
local nL,numORDS,L,num;
global xprint:
  #05/07/10: added xprint
  if not(type(xprint,boolean)) then
     xprint:=false:
  fi:
	#L is a list [ORDS(f1), ORDS(f2), ...]
	#It calculates
	#a lower bounds for sum ORD g where g = (1 + c1 f1 + c2 f2 + ...)
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("mintotORDS(L,num)                                            \n");
     printf("   L is a list (or array) [[ORDS(f1), ORDS(f2), ...].        \n");
     printf("   num = nops(L) (list) or dim(L) (array)                    \n"); 
     printf("   This proc calculates a lower bound for sum ORD g where    \n");
     printf("   g = (1 + c1 f1 + c2 f2 + ...)                             \n");
     printf("   Here sum is over cusps not equivalent to oo.              \n");
     printf("   (Usually) each ORDS(f) was produced by getaprodcuspORDS   \n");
     printf("   and the first cusp in the list is oo.                     \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 2 then
        L:=args[1]: num:=args[2]: 
	nL:=num:    
	numORDS:=nops(L[1]):
        if xprint then
          printf("min inf ord=%a\n",min(seq(L[j][1],j=1..nL)));
        fi:
        RETURN(add(min(seq(L[j][i],j=1..nL),0),i=2..numORDS));
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 2`);
  fi:
end:
thetaids[mixedjac2jac]:=proc()
local xx,n,j,jacterm,cc,jacprodtmp,jacprodtmp2,jacexp,T,curperiod;
 global xprint,jmxperiod:
        if not(type(xprint,boolean)) then
          xprint:=false:
        fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("mixedjac2jac(jacexp)                                       \n");
     printf("   jacexp is a sum of quotients of theta functions written   \n");
     printf("   in terms of JAC(a,b,infinity). jac2series and             \n"); 
     printf("   jacprodmake is used to convert each quotient into a       \n");
     printf("   a quotient of the same b.                                 \n");
     printf("   T is chosen large enough for conversion to work.          \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 1 then
        jacexp:=args[1]:  
        jmxperiod:=1:
       if type(jacexp,`+`) then
         xx:=0:
          n:=nops(jacexp):
          for j from 1 to n do
             if xprint then
               lprint("term ",j,"of ",n, "in mixedjac2jac");
             fi:
             jacterm:=op(j,jacexp):
             T:=jacbase(jacterm):
             cc:=eval(subs(JAC=JACCOF,jacterm)):
             jacprodtmp:=eval(jacterm/cc):
             jac2series(jacprodtmp,3*T+10):
             jacprodtmp2:=jacprodmake(%,q,2*T+5): 
             if jacprodtmp2=1 then
                curperiod:=1:      
             else
                curperiod:=jacp[2]:
             fi:
             jmxperiod:=ilcm(jmxperiod,curperiod):
             xx:=xx+cc*jacprodtmp2:
          od:
          RETURN(xx):
      else
         lprint("WARNING: no change");
         RETURN(jacexp):
      fi:
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 1`);
  fi:
end:
thetaids[newxy]:=proc()
local found,a,b,A,x1,y1,x,y,N;
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("newxy(x,y,N)                                             \n");
     printf("   Suppose (x,y,N)=1. This proc returns [x1,y1] such     \n");
     printf("   (x1,y1), x1 == x mod N and y1 == y mod N.             \n");
     printf("---------------------------------------------------------\n");
  elif nargs = 3 then
     #check input:
      x:=args[1]: y:=args[2]: N:=args[3]:
      if map(x->type(x,integer),[x,y,N])<>[true,true,true] then
         ERROR(`x,y,N must be integers`);
      fi:
      if not type(N,posint) then
         ERROR(`N must be a positive integer.`);
      fi:
      if igcd(x,y,N)<>1 then
         ERROR(`gcd(x,y,N) must be 1`);
      fi:
    # given (x,y,N)=1
    # find x1,y1 such that 
    # x1 == x mod N, y1 == y mod N
    # and (x1,y1)=1.
    if igcd(x,y)=1 then
        RETURN([x,y]):
    else
        found:=0:
        a:=0:
        b:=0:
        A:=0:
        while found=0 do
           x1:=x+a*N:
           y1:=y+b*N:
           if igcd(x1,y1)=1 then
              found:=1:
              RETURN([x1,y1]):
           else
              if b>0 then
                 a:=a+1:
                 b:=b-1:
              else
                 A:=A+1:
                 b:=A:
                 a:=0:
              fi:
           fi:
        od:
     fi:
  else
     ERROR(`nargs should be 0 or 3`);
  fi:
end:
thetaids[numcuspequiv1]:=proc()
local N,dd,xx,d;
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("numcuspequiv1(N)                                         \n");
     printf("   Returns the number of inequivalent cusps              \n");
     printf("   of Gamma[1](N).                                       \n");   
     printf("---------------------------------------------------------\n");
  elif nargs = 1 and type(args[1],posint) then
    N:=args[1]:
    dd:=numtheory[divisors](N):
    xx:=0:
    for d in dd do
      xx:=numtheory[phi](d)*numtheory[phi](N/d)+xx:
    od:
    RETURN(xx/2):
  else
     ERROR(`invalid input type`);
  fi:
end:
thetaids[phiset]:=proc()
local SP,j,m;
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("phiset(m)                                                \n");
     printf("   Thus proc returns the integers between 1 and m that are\n");
     printf("   relatively prime to m.                                \n");
     printf("---------------------------------------------------------\n");
  elif nargs = 1 and type(args[1],posint) then
     m:=args[1]:
     SP:={}:
     for j from 1 to m do
         if igcd(j,m)=1 then
           SP:=SP union {j}:
         fi:
     od:
     RETURN(SP):
  else
     ERROR(`invalid input type`);
  fi:
end:
 
thetaids[printJACIDORDStable]:=proc()
local jacterm,cc,jbases,jbase,jprod,bigmat,sumcolG,numfuncs,numcusps,k,j,MFLB2; 
global jacterm1, jacterm2,geprod,jaclist:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("printJACIDORDStable()                                        \n");
     printf("   Print a table of ORDS for each term in a jacprod-identity \n");
     printf("   using global data produced by the function provemodfuncid.\n");
     printf("   Table is stored in the matrix bigmat which is returned.   \n");
     printf("-------------------------------------------------------------\n");
        ## numfuncs:=nops(op(_ORDS)):
        ## 07.14.18: fixed numfuncs
        numfuncs:=nops([indices(_ORDS)]):
        numcusps:=nops(_ORDS[1]):
        bigmat:=matrix(numcusps+1,numfuncs+2):
        bigmat[1,1]:=cusp:
        for k from 1 to numfuncs do bigmat[1,k+1]:=ORD(_F[k]): od:
        bigmat[1,numfuncs+2]:=ORD(_G):
        for j from 1 to numcusps do
            bigmat[j+1,1]:=_CUSPS[j]:
            for k from 1 to numfuncs do
              bigmat[j+1,k+1]:=_ORDS[k][j]:
            od:
            bigmat[j+1,numfuncs+2]:=min(seq(_ORDS[k][j],k=1..numfuncs)):
        od:
        printf("ORDS Table for the jacprod identity\n");
        printf("_G = ");
        print(add(COFS[j]*_F[j],j=1..numfuncs)=0);
        printf("where\n");
        for k from 1 to numfuncs do
           printf("_F[%a] = \n",k);
           print(EPRODL[k]);
        od:
        ##sumcolG:=add(bigmat[j,numfuncs+2],j=2..numcusps+1):
        ##07.14.18: Fixed sumcolG
        sumcolG:=add(bigmat[j,numfuncs+2],j=3..numcusps+1):
        MFLB2:=-sumcolG:
        print(op(bigmat));
        printf("The last column of the table gives a lower bound for\n");
        printf("ORDS of _G. By summing this last column (except for oo)\n");
        printf("we see that the identity can be proved by showing that\n");
        printf("the coefficients of\n");
        printf("q^0, q^1, ... q^%a are all zero.\n",MFLB2+1);
        if MFLB2=MFLB then
           printf("This confirms the calculation done by provemodfuncid.\n");
        else
           printf("WARNING: Previous statement suspect since \n");
           printf("         MFLB2 = -sumcolG = %a and MFLB = %a.\n",MFLB2,MFLB);
           printf("         These two quantities should be equal.\n");    
        fi:
        RETURN(bigmat):
 fi:
end:
#++++++++++++++++++++++++++++++++++++++++++++++++++++
thetaids[processjacid]:=proc()
local xx,n,j,jacterm,cc,jacprodtmp,jacprodtmp2,jacexp,T,ratjac,ratjac2;
local newjacid,nj,tdq,tdq1,tj,newjacid0;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("processjacid(ratjac)                                       \n");
     printf("   ratjac is rational function of theta functions.           \n");
     printf("   This function calls mixedjac2jac.\n");
     printf("   This proc expands the numerator a jacprodcombo.   \n");
     printf("   Then uses mixedjac2jac and                        \n");
     printf("   renormalizes new jacprodcombo by dividing by the term \n");
     printf("   with the lowest power of q. \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 1 then
        ratjac:=args[1]:  
        ratjac2:= normal(ratjac):
        jacexp:=expand(numer(ratjac2)):
        ##OLD: newjacid:=mixedjac2jac(jacexp,T):
        newjacid:=mixedjac2jac(jacexp):
        nj:=nops(newjacid):
        tdq:=1000:
        for j from 1 to nj do
            ##tdq1:=qdegree(op(j,newjacid)):
            ##Changed tdq1 12/19/12 to handle fractional exponents of q
            tdq1:=lqdegree0(op(j,newjacid)):
            if tdq1<tdq then
               tdq:=tdq1:
               tj:=j:
            fi:
        od:
        newjacid0:=expand(newjacid/op(tj,newjacid)):
        newjacid:=mixedjac2jac(newjacid0):
        RETURN(newjacid):
  fi:
end:
thetaids[provemodfuncid]:=proc()
local nn,j,XX,cc,EPROD,ordtmp,inford,mftmp,num1,num2,num3,proofq,ff1,ss,desq,jacid,N,CUSPS,WIDS,CW;
    local cc2:
    global qcheck,modfunccheck, totcheck, _ORDS,jptmp,jpqd,eptmp,gltmp,EPRODL,GETAL,COFS,conpres,CONTERMS,mintottmp:
    global MFLB:
    ##
    ## proving jacid as a modfunc on Gamma1(N)
    ##
 global xprint, overrideproofq,_CUSPS:
        if not(type(xprint,boolean)) then
          xprint:=false:
        fi:
        if not(type(overrideproofq,boolean)) then
          overrideproofq:=false:
        fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("provemodfuncid(jacid,N)                                      \n");
     printf("   jacid = sum of modular functions on Gamma[1](N)           \n");
     printf("           Each term in the sum is a JAC-quotient to base N. \n"); 
     printf("   CUSPS = Set of inequivalent cusps for Gamma[1](N).        \n");
     printf("   WIDS  = List of corresponding widths.                     \n");
     printf("   global vars (can be used for error-checking):             \n");
     printf("   qcheck, modfunccheck, totcheck, _ORDS, jptmp, jpqd, eptmp,\n");
     printf("   gltmp, EPRODL, GETAL, COFS, conpres, CONTERMS, mintottmp  \n");
     printf("                                                             \n");
     printf("** Do you want a DETAILED DESCRIPTION of this proc? (yes/no) ");
     desq:=readline(terminal);
     printf("\n");
     if desq = "yes" then
     printf("DESCRIPTION:                                                 \n");
     printf("  (I)  We cycle through the terms of jacid.                  \n");
     printf("       Let j be term number.                                 \n");
     printf("       (1) Test if term is a constant.                       \n");
     printf("           If constant then conpres=1, j added to CONTERMS   \n");
     printf("           list, EPROD[j]=1, GETAL[j]=[], and                \n");
     printf("           _ORDS[j]=[0,0..].                                 \n");
     printf("       (2) Assuming term is not a constant.                  \n");
     printf("           Let jpqd = power of q in jacterm.                 \n");
     printf("           Use jac2eprod to convert jacterm to GETA-prod.    \n");
     printf("           Use GETAP2getalist to convert eprod to getalist.  \n");
     printf("           Use getaprodcuspORDS to calculate ORDS of jacterm.\n");
     printf("           Result is stored in the array _ORDS as _ORDS[j].  \n");
     printf("       (3) Check that the power of q matches ORD at oo.      \n");
     printf("           If not, j is added to qcheck list.                \n");
     printf("       (4) Use Gamma1ModFunc to check whether GETA-prod is a \n");
     printf("           modular function on Gamma[1](N).                  \n");
     printf("           If not, j is added to modfunccheck list.          \n");
     printf("  (II) We now should have a complete array _ORDS.            \n");
     printf("       (5) Final error check made. Each of the arrays qcheck,\n");
     printf("           modfunccheck, and totcheck should be empty. If not\n");
     printf("           an error message is returned which terminates the \n");
     printf("           proc.                                             \n");
     printf("       (6) A WARNING is printed if any terms were constants. \n");
     printf("       (7) We use mintotORDS to min power of q to check      \n");
     printf("           identity.                                         \n");
     printf("           A query is made whether to check now.             \n");
     printf("           If not (suggested), this min power is returned.   \n");
     fi:
     printf("-------------------------------------------------------------\n");
     elif nargs = 2 then
        ##jacid:=args[1]: CUSPS:=args[2]: WIDS:=args[3]: N:=args[4]:
        jacid:=args[1]: N:=args[2]:
        CW:=CUSPSANDWIDMAKE1(N):
        CUSPS:=CW[1]: WIDS:=CW[2]:
        _CUSPS:=CUSPS:
    conpres:=0:  #conpres=0 if no terms are constants
    qcheck:=[]:
    modfunccheck:=[]:
    totcheck:=[]:
    CONTERMS:=[]:
    nn:=nops(jacid):
    _ORDS:=array(1..nn):
    EPRODL:=array(1..nn):
    GETAL:=array(1..nn):
    COFS:=array(1..nn):
    for j from 1 to nn do
      lprint("TERM ",j,"of ",nn," *****************");
      XX:=op(j,jacid):
      if xprint then
        print("XX=",XX);
      fi:
      ##cc:=eval(subs(JAC=JACCOF,XX)):
      cc:=eval(subs(JAC=JACCOF,XX)):
      cc2:=eval(subs({JAC=JACCOF,q=1},XX)):
      COFS[j]:=cc2:
      jptmp:=eval(XX/cc2):
      #11/03/09: added condition to check term is a constant
      if jptmp=1 then
          EPRODL[j]:=1:
          GETAL[j]:=[]:
          _ORDS[j]:=[seq(0,k=1..nops(CUSPS))]:
          conpres:=1:
          CONTERMS:=[op(CONTERMS),j]:
      else
          jpqd:=degree(cc,q);
          eptmp:=jac2eprod(jptmp):
          EPRODL[j]:=eptmp:
          gltmp:=GETAP2getalist(eptmp):
          GETAL[j]:=gltmp:
          ordtmp:=getaprodcuspORDS(gltmp,CUSPS,WIDS):
          _ORDS[j]:=ordtmp:
          if toterror<>0 then
            lprint("WARNING: totcheck: ");
            totcheck:=[op(totcheck),j]:
          fi:
          inford:=ordtmp[1]:
          if inford<>jpqd then
            lprint("WARNING: qcheck: ");
            qcheck:=[op(qcheck),j]:
          else
            if xprint then
             lprint("POWER of q CORRECT");
            fi:
          fi:
          mftmp:=Gamma1ModFunc(gltmp,N):
          if mftmp<>1 then
            lprint("WARNING: modfunccheck: ");
            modfunccheck:=[op(modfunccheck),j]:
          fi:
      fi:
   od:
   mintottmp:=mintotORDS(_ORDS,nn): 
   lprint("mintotord = ",mintottmp);
   lprint("TO PROVE the identity we need to show that v[oo](ID) > ",-mintottmp);
   num1:=nops(qcheck):
   num2:=nops(modfunccheck):
   num3:=nops(totcheck):
   if num1+num2+num3>0 then
     lprint("***********************************");
     lprint("See qcheck, modfunccheck, totcheck.");
     error("There were errors.");
   else
     printf("*** There were NO errors. \n");
     printf("*** o Each term was modular function on\n");
     printf("      Gamma1(%a). \n",N);
     printf("*** o We also checked that the total order of\n",N);
     printf("      each term was zero.\n");
     printf("*** o We also checked that the power of q was correct in\n");
     printf("      each term.\n");
   fi:
   if conpres=1 then
     lprint("*** WARNING: some terms were constants. ***");
     lprint("See array CONTERMS.");
   fi:                    
   printf("To prove the identity we will need to verify if up to \n");
   printf("q^(%a).\n",-mintottmp);
   printf("Do you want to prove the identity? (yes/no)\n");
   if overrideproofq then
      proofq := "yes":
   else
      proofq:=readline(terminal);
   fi:
   MFLB:=-mintottmp:
     if proofq = "yes" then
        printf("You entered yes.\n");
        printf("We verify the identity to O(q^(\%a)).\n",-mintottmp+2*N);
        ff1:=jac2series(jacid,-mintottmp+2*N):
         ss:=simplify(series(ff1,q,-mintottmp+2*N));
         ss:=convert(ss,polynom):
         if ss=0 then
           printf("RESULT: The identity holds to O(q^(\%a)).\n",-mintottmp+2*N);
           printf("CONCLUSION: This proves the identity since we had only\n");
           printf("            to show that v[oo](ID) > %a.\n",MFLB);  
         else 
           printf("WARNING: ss <> 0. Something is wrong.\n");
           printf("ss = %a. Something is wrong.\n",ss);
         fi:
     else
        printf("You did not enter yes.\n");
     fi:
   RETURN(-mintottmp);
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 2`);
  fi:
end:
thetaids[provemodfuncidBATCH]:=proc()
local nn,j,XX,cc,EPROD,ordtmp,inford,mftmp,num1,num2,num3,proofq,ff1,ss,desq,jacid,N,CUSPS,WIDS,CW;
    local cc2,pp,symid:
    global qcheck,modfunccheck, totcheck, _ORDS,jptmp,jpqd,eptmp,gltmp,EPRODL,GETAL,COFS,conpres,CONTERMS,mintottmp:
    global MFLB:
    ##
    ## proving jacid as a modfunc on Gamma1(N)
    ##
 global xprint, overrideproofq,_CUSPS,noprint:
 global qthreshold,proveit,printlocalsymid:
 ## 06.26.16 added global var printlocalsymid
 ## 07.15.18 added global var noprint
        if not(type(xprint,boolean)) then
          xprint:=false:
        fi:
        if not(type(noprint,boolean)) then
          noprint:=false:
        fi:
        if not(type(qthreshold,integer)) then
          qthreshold:=1000: 
        fi:
        if not(type(proveit,boolean)) then
          proveit:=false:
        fi:
        if not(type(overrideproofq,boolean)) then
          overrideproofq:=false:
        fi:
        if not(type(printlocalsymid,boolean)) then
          printlocalsymid:=false:
        fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("provemodfuncidBATCH(symid,jacid,N)                           \n");
     printf("provemodfuncidBATCH(jacid,N)                           \n");
     printf("   This is a BATCH version od provemodfuncid, \n");
     printf("   It uses the global var qthreshold .\n");
     printf("   symid = symbolic form of the identity. \n");
     printf("   jacid = sum of modular functions on Gamma[1](N)           \n");
     printf("           Each term in the sum is a JAC-quotient to base N. \n"); 
     printf("   CUSPS = Set of inequivalent cusps for Gamma[1](N).        \n");
     printf("   WIDS  = List of corresponding widths.                     \n");
     printf("   global vars (can be used for error-checking):             \n");
     printf("   qcheck, modfunccheck, totcheck, _ORDS, jptmp, jpqd, eptmp,\n");
     printf("   gltmp, EPRODL, GETAL, COFS, conpres, CONTERMS, mintottmp,  \n");
     printf("   printlocalsymid                                           \n");
     printf("                                                             \n");
     printf("** Do you want a DETAILED DESCRIPTION of this proc? (yes/no) ");
     desq:=readline(terminal);
     printf("\n");
     if desq = "yes" then
     printf("DESCRIPTION:                                                 \n");
     printf("  (I)  We cycle through the terms of jacid.                  \n");
     printf("       Let j be term number.                                 \n");
     printf("       (1) Test if term is a constant.                       \n");
     printf("           If constant then conpres=1, j added to CONTERMS   \n");
     printf("           list, EPROD[j]=1, GETAL[j]=[], and                \n");
     printf("           _ORDS[j]=[0,0..].                                 \n");
     printf("       (2) Assuming term is not a constant.                  \n");
     printf("           Let jpqd = power of q in jacterm.                 \n");
     printf("           Use jac2eprod to convert jacterm to GETA-prod.    \n");
     printf("           Use GETAP2getalist to convert eprod to getalist.  \n");
     printf("           Use getaprodcuspORDS to calculate ORDS of jacterm.\n");
     printf("           Result is stored in the array _ORDS as _ORDS[j].  \n");
     printf("       (3) Check that the power of q matches ORD at oo.      \n");
     printf("           If not, j is added to qcheck list.                \n");
     printf("       (4) Use Gamma1ModFunc to check whether GETA-prod is a \n");
     printf("           modular function on Gamma[1](N).                  \n");
     printf("           If not, j is added to modfunccheck list.          \n");
     printf("  (II) We now should have a complete array _ORDS.            \n");
     printf("       (5) Final error check made. Each of the arrays qcheck,\n");
     printf("           modfunccheck, and totcheck should be empty. If not\n");
     printf("           an error message is returned which terminates the \n");
     printf("           proc.                                             \n");
     printf("       (6) A WARNING is printed if any terms were constants. \n");
     printf("       (7) We use mintotORDS to min power of q to check      \n");
     printf("           identity.                                         \n");
     printf("           A query is made whether to check now.             \n");
     printf("           If not (suggested), this min power is returned.   \n");
     printf("  (III) RETURNS 1 if identity proved.                        \n");
     fi:
     printf("-------------------------------------------------------------\n");
     elif nargs = 2 or nargs = 3 then
        ##jacid:=args[1]: CUSPS:=args[2]: WIDS:=args[3]: N:=args[4]:
        if nargs = 2 then
           jacid:=args[1]: N:=args[2]:
           symid:=jacid=0:
        else
           symid:=args[1]: jacid:=args[2]: N:=args[3]:
        fi:
        CW:=CUSPSANDWIDMAKE1(N):
        CUSPS:=CW[1]: WIDS:=CW[2]:
        _CUSPS:=CUSPS:
    conpres:=0:  #conpres=0 if no terms are constants
    qcheck:=[]:
    modfunccheck:=[]:
    totcheck:=[]:
    CONTERMS:=[]:
    nn:=nops(jacid):
    _ORDS:=array(1..nn):
    EPRODL:=array(1..nn):
    GETAL:=array(1..nn):
    COFS:=array(1..nn):
    for j from 1 to nn do
      ##lprint("TERM ",j,"of ",nn," *****************");
      XX:=op(j,jacid):
      if xprint then
        print("XX=",XX);
      fi:
      ##cc:=eval(subs(JAC=JACCOF,XX)):
      cc:=eval(subs(JAC=JACCOF,XX)):
      cc2:=eval(subs({JAC=JACCOF,q=1},XX)):
      COFS[j]:=cc2:
      jptmp:=eval(XX/cc2):
      #11/03/09: added condition to check term is a constant
      if jptmp=1 then
          EPRODL[j]:=1:
          GETAL[j]:=[]:
          _ORDS[j]:=[seq(0,k=1..nops(CUSPS))]:
          conpres:=1:
          CONTERMS:=[op(CONTERMS),j]:
      else
          jpqd:=degree(cc,q);
          eptmp:=jac2eprod(jptmp):
          EPRODL[j]:=eptmp:
          gltmp:=GETAP2getalist(eptmp):
          GETAL[j]:=gltmp:
          ordtmp:=getaprodcuspORDS(gltmp,CUSPS,WIDS):
          _ORDS[j]:=ordtmp:
          if toterror<>0 then
            lprint("WARNING: totcheck: ");
            totcheck:=[op(totcheck),j]:
          fi:
          inford:=ordtmp[1]:
          if inford<>jpqd then
            lprint("WARNING: qcheck: ");
            qcheck:=[op(qcheck),j]:
          else
            if xprint then
             lprint("POWER of q CORRECT");
            fi:
          fi:
          mftmp:=Gamma1ModFunc(gltmp,N):
          if mftmp<>1 then
            lprint("WARNING: modfunccheck: ");
            modfunccheck:=[op(modfunccheck),j]:
          fi:
      fi:
   od:
   mintottmp:=mintotORDS(_ORDS,nn): 
   ##lprint("mintotord = ",mintottmp);
   ##lprint("TO PROVE the identity we need to show that v[oo](ID) > ",-mintottmp);
   if -mintottmp > qthreshold then
     printf(" abs(mintotord)=%a, which is too large \n",mintottmp);
     printf(" Try increasing the global var qthreshold. \n");
     RETURN():
   fi:
   num1:=nops(qcheck):
   num2:=nops(modfunccheck):
   num3:=nops(totcheck):
   if num1+num2+num3>0 then
     lprint("***********************************");
     lprint("See qcheck, modfunccheck, totcheck.");
     error("There were errors.");
   else
     if noprint=false then
       printf("*** There were NO errors.  Each term was modular function on\n");
       printf("    Gamma1(%a). Also -mintotord=%a. To prove the identity\n",N,-mintottmp);
       printf("    we need to  check up to O(q^(\%a)).\n",-mintottmp+2);
       printf("    To be on the safe side we check up to O(q^(\%a)).\n",-mintottmp+2*N);
     fi:
   fi:
   if conpres=1 and xprint then
     lprint("*** WARNING: some terms were constants. ***");
     lprint("See array CONTERMS.");
   fi:                    
   MFLB:=-mintottmp:
   ff1:=jac2series(jacid,-mintottmp+2*N):
   ss:=simplify(series(ff1,q,-mintottmp+2*N));
   pp:=convert(ss,polynom);
   if pp=0 then
     if noprint=false then
       printf("*** The identity below is PROVED!\n");
     fi:
   if printlocalsymid then
     print(symid);
   fi:
   ##added 07.12.16
     RETURN(1):
   else
     printf("*** The identity below is FALSE!\n");
     printf("    It does NOT hold up to O(q^(\%a)).\n",-mintottmp+2*N);
     if printlocalsymid then
       print(symid);
     fi:
     RETURN(0):
   fi:
   ## 06.26.16
   if printlocalsymid then
     print(symid);
   fi:
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 2`);
  fi:
end:
##         if ss=0 then
##           printf("RESULT: The identity holds to O(q^(\%a)).\n",-mintottmp+2*N);
##           printf("CONCLUSION: This proves the identity since we had only\n");
##           printf("            to show that v[oo](ID) > %a.\n",MFLB);  
##         else 
##           printf("WARNING: ss <> 0. Something is wrong.\n");
##           printf("ss = %a. Something is wrong.\n",ss);
##         fi:
##     else
##        printf("You did not enter yes.\n");
##     fi:
##   RETURN(-mintottmp);
##  else
##     printf("nargs=%a\n",nargs);
##     ERROR(`nargs must be 0 or 2`);
##  fi:
##end:
thetaids[qjacdegree]:=proc()
local cc,jacterm,X,R:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("qjacdegree(jacterm)                                            \n");
     printf("   Computes degree of jacterm with respect to q.             \n");
     printf("   jacterm is a quotient of JAC(a,b,infinity)s               \n");
     printf("   In other words,                                           \n");
     printf("   c*q^r*JAC(a1,b1)*JAC(a2,b2)*.. --> r    \n");
     printf("   Here r is a rational number.            \n");
     printf("-------------------------------------------------------------\n");
## EXAMPLE:
## > jacterm:=q^(749/24)*JAC(0, 10, infinity)*JAC(0, 120, infinity)^3
##            *JAC(0, 360, infinity)*JAC(1, 3, infinity)*JAC(2, 6, infinity)
##            *JAC(3, 10, infinity)*JAC(4, 12, infinity)*JAC(8, 20, infinity)^3;
##             /749\
##             |---|
##             \24 /                                           3
## jacterm := q      JAC(0, 10, infinity) JAC(0, 120, infinity)
## 
##     JAC(0, 360, infinity) JAC(1, 3, infinity) JAC(2, 6, infinity)
## 
##                                                                   3
##     JAC(3, 10, infinity) JAC(4, 12, infinity) JAC(8, 20, infinity)
## 
## > qjacdegree(jacterm);
##                                       749
##                                       ---
##                                       24
## 
     elif nargs = 1 then
       jacterm:=args[1]:
       if degree(jacterm,q)=FAIL then
          X:=eval(subs(JAC=JACCOF,jacterm));
          if X=0 then
            RETURN(0):
          else
            R:=simplify(eval(q*diff(X,q)/X)):
            if type(R,rational) then
               RETURN(R);
            else
               ERROR(`qjacdegree is not rational`);
            fi:
          fi:
       else
          RETURN(degree(X,q));
       fi:
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 1`);
  fi:
end:
thetaids[QP2]:=x->fpart(x)^2 - fpart(x) + 1/6:
thetaids[rmcofjac]:=proc()
local cc,jacterm;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("rmcofjac(jacterm)                                            \n");
     printf("   Removes constant coefficient of a jacterm.                \n");
     printf("   jacterm is a quotient of JAC(a,b,infinity)s               \n");
     printf("   In other words,                                           \n");
     printf("   c*JAC(a1,b1)*JAC(a2,b2)*.. --> JAC(a1,b1)*JAC(a2,b2)*..   \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 1 then
        jacterm:=args[1]: 
             cc:=eval(subs(JAC=JACCOF,jacterm)):
             RETURN(eval(jacterm/cc));
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 1`);
  fi:
end:
thetaids[rmcofnotqjac]:=proc()
local cc,jacterm;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("rmcofnotqjac(jacterm)                                         \n");
     printf("   Removes only constant coefficient of a jacterm and leave   \n");
     printf("   power of q intact.                                          \n");
     printf("   jacterm is a quotient of JAC(a,b,infinity)s               \n");
     printf("   In other words,                                           \n");
     printf("   c*q^d*JAC(a1,b1)*JAC(a2,b2)*..-->q^d*JAC(a1,b1)*JAC(a2,b2)*..\n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 1 then
        jacterm:=args[1]: 
             cc:=eval(subs({JAC=JACCOF,q=1},jacterm)):
             RETURN(eval(jacterm/cc));
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 1`);
  fi:
end:
thetaids[Scmake]:=proc()
local S1,S,sp,s,c,N;
    # c is a positive divisor of N
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("Scmake(c,N)                                              \n");
     printf("   Returns the set S[c] where c is a positive            \n");
     printf("   divisor of N.                                         \n");
     printf("---------------------------------------------------------\n");
  elif nargs = 2 then
     #check input:
      c:=args[1]: N:=args[2]:
      if map(x->type(x,posint),[c,N])<>[true,true] then
         ERROR(`c,N must be positive integers`);
      fi:
      if modp(N,c)<>0 then
         ERROR(`c must be a positive divisor of N`);
      fi:
    S1:=sset(N/c):
    S:={}:
    for sp in S1 do
        s:=srep(sp,N/c,N):
        S:=S union {s}:
    od:
    RETURN(S):
  else
     ERROR(`nargs should be 0 or 2`);
  fi:
end:
thetaids[srep]:=proc()
local find,j,news,s,c,N;
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("srep(s,c,N)                                              \n");
     printf("   Let c be a positive divisor of N. Given s in (Z/cZ)^x \n");
     printf("   this proc returns an s1 in (Z/NZ)^x such that         \n");
     printf("   s1 == s mod c.                                        \n"); 
     printf("---------------------------------------------------------\n");
  elif nargs = 3 and map(x->type(x,posint),[args])=[true,true,true] then
     # c is a positive divisor of N
     # given s in Z[c]x find s1 in Z[n]x such that s1=s mod c
     s:=args[1]: c:=args[2]: N:=args[3]:
     if modp(N,c)<>0 then
        ERROR(` srep: c must be a divisor of N`);
     fi:
     if igcd(s,c)<>1 then
        ERROR(` srep: s, c must be relatively prime`);
     fi:
     find:=0:
     for j from 0 to trunc(N/c) do
        news:=j*c+s:
        if igcd(news,N)=1 then
           find:=1:
           RETURN(news):
        fi:
     od:
     if find=0 then
     lprint("ERROR in srep with args",s,c,N):
     error "found=0";
    fi:
  else
     ERROR(`invalid input type`);
  fi:
end;
thetaids[sset]:=proc()
local SP,j,m;
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("sset(m)                                                  \n");
     printf("   Returns a set of representives for the following      \n");
     printf("   equivalence relation on the set of integers relatively\n");
     printf("   prime to m. We say to two such integers n1, n2 are    \n");
     printf("   equivalent if n1 == +/- n2 mod m.                     \n");
     printf("---------------------------------------------------------\n");
  elif nargs = 1 and type(args[1],posint) then
     m:=args[1]:
     if m=1 or m=2 then
        SP:={1}:
     else
        SP:={}:
        for j from 1 to trunc(m/2) do
            if igcd(j,m)=1 then
              SP:=SP union {j}:
            fi:
       od:
     fi:
     RETURN(SP):
  else
     ERROR(`invalid input type`);
  fi:
end:
thetaids[thetaidschanges]:=proc()
printf("**************************************************************\n");
printf("*\n");
printf("*\n");
printf("*   thetaids package version 0.8 -  Sun Jun 26 14:08:15 EDT 2016\n");
printf("*   thetaids package version 0.9 -  Tue Jul 12 13:39:16 EDT 2016\n");
printf("*   thetaids package version 1.0 -  Sat Jul 21 12:00:05 EDT 2018\n");
printf("*  This version tested on MAPLE 2017\n");
printf("*\n");
printf("*\n");
printf("*  Changes since previous version 0.9\n");
printf("*\n");
printf("*     * Fixed two bugs in printJACIDORDStable\n");
printf("*       (was giving incorrect numfuncs)\n");
printf("*       (need to omit oo entry for sumcolG)\n");
printf("*\n");
printf("*     * Added noprint global var to\n");
printf("*       provemodfuncidBATCH so when noprint is true\n");
printf("*       it only RETURNS 1 if identity is true.     \n");
printf("*\n");
printf("*     * Updated thetaids/ramarobinsids tutorial\n");
printf("*\n");
printf("*\n");
printf("*  Changes since previous version 0.8\n");
printf("*\n");
printf("*     * Changed pversion to thetaidspversion\n");
printf("*\n");
printf("*     * provemodmfundidBATCH: \n");            
printf("*       RETURNS 1 if identity proved\n");
printf("*\n");
printf("*\n");
printf("*\n");
printf("*  Changes since previous version 0.7\n");
printf("*\n");
printf("*      UNDER CONSTRUCTION\n");
printf("*\n");
printf("**************************************************************\n");
RETURN():
end:
thetaids[thetaidspversion]:=proc()
printf("****************************************************\n");
printf("*\n");
printf("*  thetaids package version 1.0  \n");
printf("*  Sat Jul 21 12:00:05 EDT 2018\n");
printf("*  This version tested on MAPLE 2017\n"); 
printf("*\n");
printf("*  Please report any problems to fgarvan@ufl.edu\n");
printf("*  Previous versions:\n");
printf("        0.9 - Jul 2016 (MAPLE 2015)\n");
printf("        0.8 - Jun 2016 (MAPLE 2015)\n");
printf("        .........\n");
printf("****************************************************\n");
RETURN():
end:
thetaids[v0]:=proc()
local nL,x,j,v,n,m,r,L,N;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("v0(L,N)                                               \n");
     printf("   L is a getalist; i.e. L=[[b1,a1,c1],....]\n");
     printf("   Let G be the corresponding generalized eta-function.\n");
     printf("   v0(L,N) returns ord(G,0)\n");                      
     printf("-------------------------------------------------------------\n");
  elif nargs = 2 and type(args[2],posint) then
  L:=args[1]: N:=args[2]:
  nL:=nops(L):
  x:=0:
  for j from 1 to nL do
      v:=L[j]:
      n:=v[1]:m:=v[2]:r:=v[3]:
      x := x + N*QP2(0)*r/n:
  od:
  RETURN(x):
  else
     ERROR(`invalid input type`);
  fi:
end:
thetaids[vinf]:=proc()
local nL,x,j,v,n,m,r,L,N;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("vinf(L,N)                                               \n");
     printf("   L is a getalist; i.e. L=[[b1,a1,c1],....]\n");
     printf("   Let G be the corresponding generalized eta-function.\n");
     printf("   vinf(L,N) returns ord(G,oo)\n");                      
     printf("-------------------------------------------------------------\n");
  elif nargs = 2 and type(args[2],posint) then
  L:=args[1]: N:=args[2]:
  nL:=nops(L):
  x:=0:
  for j from 1 to nL do
      v:=L[j]:
      n:=v[1]:m:=v[2]:r:=v[3]:
      x := x + n*QP2(m/n)*r:
  od:
  RETURN(x):
  else
     ERROR(`invalid input type`);
  fi:
end:
printf("TABLE TYPE thetaids = %a\n",type(thetaids,table));
savelib( thetaids 
,
"c:\\cygwin64/home/Owner/maple/mylib/thetaids.mla");
##"H:\\maple\\newmylib\\thetaids.mla");
printf("END thetaids package\n");
## mylib  above must be changed to the name of the directory
## in which you want stuff stored