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:
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:
CUSPSANDWIDMAKE1:=proc()
     local cuspsN,widsN,CUSPSN:
     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:
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:
	GETAP2getalist:=proc()
## 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,geprod:
  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 "chk<>0";
	            fi:  
		    GP:=[op(GP),[s,r]]:
		od:
		RETURN(GP);
  else
     ERROR(`invalid input type`);
  fi:
	end:	
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;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("Gamma1ModFunc(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("   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]:
  mfunc:=1:
  if DivCheck(L,N)=1 then
     print("All n are divisors of ",N);
  else
     print("All n are NOT divisors of ",N);
     print("Hence NOT a modfunc on Gamma1(",N,")");
     mfunc:=0:
  fi:
  val0:=v0(L,N):
  print("val0=",val0);
  if type(val0,integer) then
    #OK
    if modp(val0,2)=0 then
      print("which is even.");
    else
      print("which is odd. Hence NOT a modfunc on Gamma1(",N,")");
      mfunc:=0:
    fi:
  else
    print("which is NOT an integer. Hence NOT a modfunc on Gamma1(",N,")");
    mfunc:=0:
  fi:
  valinf:=vinf(L,N):
  print("valinf=",valinf);
  if type(valinf,integer) then
    #OK
    if modp(valinf,2)=0 then
      print("which is even.");
    else
      print("which is odd. Hence NOT a modfunc on Gamma1(",N,")");
      mfunc:=0:
    fi:
  else
    print("which is NOT an integer. Hence NOT a modfunc on Gamma1(",N,")");
    mfunc:=0:
  fi:
  if mfunc=1 then
     print("It IS a modfunc on Gamma1(",N,")");
  fi:
  RETURN(mfunc):
  else
     ERROR(`invalid input type`);
  fi:
end:
JACCOF:=(a,b,c)->1:
QP2:=x->fpart(x)^2 - fpart(x) + 1/6:
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:
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:
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:
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 unequivalent 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:
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:
ejac:=proc(a,b,c)
if a=0 then
	EETA(b):
else
	GETA(b,a)*EETA(b):
fi:
end:
fpart:=t->t-trunc(t):
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:
getaprodcuspORDS:=proc()
local tot,nc,nL,ORDS,j,a,c,kco,k,gg,n,r,mm,L,cusps,wids;
global toterror:
  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:
  print("Cusp ords: ");
  print([seq([cusps[j],ORDS[j]],j=1..nc)]);
  print("TOTAL ORD = ",tot);
  toterror:=tot:
RETURN(ORDS):
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 3`);
  fi:
end:
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:
#> 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)
#> jac2eprod(JJ);
#                               GETA(5, 2) EETA(5)
  elif nargs = 1 then
     JJ:=args[1]:
     eval(subs(JAC = `ejac`, JJ)):
  else
     ERROR(`nargs should be 0 or 1`);
  fi:
end:
lw:=proc(N)
local found,a,b,A;
    # given (x,y,N)=1
    # find x1,y1 such that 
    # x1 == x mod N, y1 == y mod N
    # and (x1,y1)=1.
        found:=0:
        a:=0:
        b:=0:
        A:=0:
        while found=0 do
        lprint(a,b);
              if b>0 then
                 a:=a+1:
                 b:=b-1:
              else
                 A:=A+1:
                 b:=A:
                 a:=0:
              fi:
        if (a+b)>N then
           found:=1:
        fi:
        od:
end:
mintotORDS:=proc()
# 09/06/09: version of mintotORD2 (mintotORD2 does not work if L in an array)
local nL,numORDS,L,num;
	#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]):
        print("min inf ord=",min(seq(L[j][1],j=1..nL)));
        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:
mixedjac2jac:=proc()
local xx,n,j,jacterm,cc,jacprodtmp,jacprodtmp2,jacexp,T;
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("mixedjac2jac(jacexp,T)                                       \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 = 2 then
        jacexp:=args[1]: T:=args[2]: 
       if type(jacexp,`+`) then
         xx:=0:
          n:=nops(jacexp):
          for j from 1 to n do
             lprint("term ",j,"of ",n);
             jacterm:=op(j,jacexp):
             cc:=eval(subs(JAC=JACCOF,jacterm)):
             jacprodtmp:=eval(jacterm/cc):
             jac2series(jacprodtmp,3*T):
             jacprodtmp2:=jacprodmake(%,q,T): 
             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 2`);
  fi:
end:
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:
numcos0:=proc(N)
local d,pd,xx,p;
   d:=numtheory[divisors](N):
   pd:=select(isprime,d):
   xx:=1:
   for p in pd do
     xx:=xx*(1+1/p):
   od:
   RETURN(xx*N):
end:      
numcos01:=proc(N)
local d,pd,xx,p;
   d:=numtheory[divisors](N):
   pd:=select(isprime,d):
   xx:=1:
   for p in pd do
     xx:=xx*(1-1/p):
   od:
   RETURN(xx*N):
end:      
numcuspequiv0:=proc()
local N,dd,xx,d;
  if nargs=0 then
     printf("---------------------------------------------------------\n");
     printf("numcuspquiv0(N) returns the number of inequivalent cusps\n");
     printf("of Gamma[0](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:=xx+numtheory[phi](igcd(d,N/d)):                    
     od:
     RETURN(xx):
  else
     ERROR(`invalid input type`);
  fi:
end:
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:
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:
provemodfuncid:=proc()
local nn,j,XX,cc,EPROD,ordtmp,inford,mftmp,num1,num2,num3,proofq,ff1,ss,desq,jacid,CUSPS,WIDS,N;
    global qcheck,modfunccheck, totcheck, _ORDS,jptmp,jpqd,eptmp,gltmp,EPRODL,GETAL,COFS,conpres,CONTERMS,mintottmp:
    ##
    ## proving jacid as a modfunc on Gamma1(N)
    ##
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("provemodfuncid(jacid,CUSPS,WIDS,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 = 4 then
        jacid:=args[1]: CUSPS:=args[2]: WIDS:=args[3]: N:=args[4]:
    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):
      print("XX=",XX);
      cc:=eval(subs(JAC=JACCOF,XX)):
      COFS[j]:=cc:
      jptmp:=eval(XX/cc):
      #11/03/09: added condition to check term is a constant
      if jptmp=1 then
          EPROD[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
            lprint("POWER of q CORRECT");
          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
     lprint("*** There were NO errors. ***");
   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");
   proofq:=readline(terminal);
     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));
        print(ss);
     else
        printf("You did not enter yes.\n");
     fi:
   RETURN(-mintottmp);
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 4`);
  fi:
end:
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:
setcusps:=proc(N)
local SC,a,b;
     SC:=[]:
     for a from 0 to N-1 do
     for b from 0 to N-1 do
       if igcd(a,b,N)=1 then
          SC:=[op(SC),[a,b]]:
       fi:
     od:
     od:
RETURN(SC):
end:
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;
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:
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:
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: