#######################################################
## make-win-package ETA
## Sat, Jun 29, 2019  9:28:11 PM
#######################################################
#######################################################

printf("BEGIN ETA package\n");
printf("THIS VERSION DATED Sat, Jun 29, 2019  9:28:11 PM \n");
ETA:=table();
ETA[cuspmake]:=proc()
##		 Computes a set of inequivalent cusps for GAMMA_0(N)
		local S,SoD,c,a,lasta, SSc,lastd,gcN,d,md,N:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("cuspmake(N)                                                  \n");
     printf("   N     =  postive integer                                  \n");
     printf("   Returns a set of inequivalent cusps for GAMMA[0](N).     \n"); 
     printf("                                                             \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 1 then
        N:=args[1]:
		SoD:=numtheory[divisors](N):
		SoD := SoD minus {1}:
	 	S:={0}:
		for c in SoD do
                   SSc:={}:
		   lastd:=c-1:  
		   gcN:=igcd(c,N/c):
		   for d from 1 to lastd do
			md:=modp(d,gcN):
			if igcd(d,c)=1 and member(md,SSc)=false then
			   S:= S union {d/c}:
			   SSc:= SSc union {md}:
			fi:
		   od:
		od:
		RETURN(S):
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 2`);
  fi:
end:
ETA[cuspord]:=proc()
		# Computes the invariant order at a cusp 
		# of the given etaproduct.
		# 06.27.19: 3 args to 2 args
		# 06.27.19: gcd -> igcd         
		local GP,ngp,S,s,ord,i,t,r,c,etaprod,N,cusp:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("cuspord(etaprod,cusp)                                      \n");
     printf("   etaprod = eta-quotient                                    \n");
     printf("   cusp  =  rational                                          \n");
     printf("   Computes the invariant order at a cusp                    \n");
     printf("   of the given etaproduct.                                   \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 2 then
        etaprod:=args[1]:
        cusp:=args[2]:
		GP:=ETA[GPmake](etaprod):
		ngp:=nops(GP):
		  ord:=0:
		  for i from 1 to (ngp/2) do
		     t:=GP[2*i-1]:
		     r:=GP[2*i]:
		     c:=denom(cusp):
		     ord:=ord+igcd(t,c)^2/t*r/24:
		  od:
		RETURN(ord):
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 2`);
  fi:
end:
ETA[cuspORD]:=proc()
  local etaprod, N, cusp:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("cuspORD(etaprod,N,cusp)                                      \n");
     printf("   etaprod = eta-quotient                                    \n");
     printf("   N     =  postive integer                                  \n");
     printf("   cusp  =  rational                                          \n");
     printf("   Computes the order at a cusp z with respect to the group \n");
     printf("   G= GAMMA_0(N) of the given etaproduct f; ie. ORD[G](f,z) \n"); 
     printf("-------------------------------------------------------------\n");
     elif nargs = 3 then
        etaprod:=args[1]:
        N:=args[2]:
        cusp:=args[3]:
	RETURN(cuspord(etaprod,cusp)*fanwidth(cusp,N)):
  fi:
end:
ETA[cuspORDS]:=proc()
		# Computes the order at each cusp (GAMMA_0(N))
		# of the given etaproduct and return result as a list.
		#
		local GP,ngp,S,s,ords,i,t,r,c, TOTORD,ORDS,cuspORDSLIST,ORD:
		local etaprod,CUSPS,N:
		global toterror:  #for error checking
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("cuspORDS(etaprod,CUSPS,N)                                    \n");
     printf("   etaprod = eta-quotient                                    \n");
     printf("   CUSPS =  list of cusps of GAMMA0(N)                        \n");
     printf("   N     =  postive integer                                  \n");
     printf("   Computes the order at each cusp z with respect to the group\n");
     printf("   G= GAMMA_0(N) of the given etaproduct f.                 \n"); 
     printf("   A list [[cusp1, ORD1], [cusp2, ORD2], ...] is returned.  \n"); 
     printf("   Global variable toterror (total ORD) should be zero.  \n"); 
     printf("-------------------------------------------------------------\n");
     elif nargs = 3 then
        etaprod:=args[1]:
        CUSPS:=args[2]:
        N:=args[3]:
                cuspORDSLIST:=[]:
		GP:=ETA[GPmake](etaprod):
		ngp:=nops(GP):
                TOTORD:=0:
		if type(etaprod,constant) then
                    cuspORDSLIST:=[seq([s,0],s in CUSPS)]:
                else
		for s in CUSPS do
		  ords:=0:
		  for i from 1 to (ngp/2) do
		     t:=GP[2*i-1]:
		     r:=GP[2*i]:
		     c:=denom(s):
		     ords:=ords+igcd(t,c)^2/t*r/24:
		  od:
                  ORD:=ords*fanwidth(s,N):
                  cuspORDSLIST:=[op(cuspORDSLIST),[s,ORD]]:
                  TOTORD:=TOTORD + ORD:
		od:	
		fi:
                toterror:=TOTORD:
		RETURN(cuspORDSLIST):
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 2`);
  fi:
end:
ETA[cuspORDSnotoo]:=proc()
		# Computes the order at each cusp<> oo  (GAMMA_0(N))
		# of the given etaproduct and return result as a list.
		#
		local GP,ngp,S,s,ords,i,t,r,c,TOTORD,ORDS, cset,cuspORDSLIST;
	        local ORD, etaprod, CUSPS, N:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("cuspORDSnotoo(etaprod,CUSPS,N)                               \n");
     printf("   etaprod = eta-quotient                                    \n");
     printf("   CUSPS =  list of cusps of GAMMA0(N) excluding oo ~ 1/N    \n");
     printf("   N     =  postive integer                                  \n");
     printf("   Computes the order at each cusp z with respect to the group\n");
     printf("   G= GAMMA_0(N) of the given etaproduct f.                 \n"); 
     printf("   A list [[cusp1, ORD1], [cusp2, ORD2], ...] is returned.  \n"); 
     printf("-------------------------------------------------------------\n");
     elif nargs = 3 then
        etaprod:=args[1]:
        CUSPS:=args[2]:
        N:=args[3]:
                cset:=convert(CUSPS,set):
                if member(1/N,CUSPS) then
                  ERROR("oo in CUSPS");
                fi:
                cuspORDSLIST:=[]:
		if type(etaprod,constant) then
                    cuspORDSLIST:=[seq([s,0],s in CUSPS)]:
                else
   		   GP:=ETA[GPmake](etaprod):
   		   ngp:=nops(GP):
   		   for s in CUSPS do
   		     ords:=0:
   		     for i from 1 to (ngp/2) do
   		        t:=GP[2*i-1]:
   		        r:=GP[2*i]:
   		        c:=denom(s):
   		        ords:=ords+igcd(t,c)^2/t*r/24:
   		     od:
                     ORD:=ords*fanwidth(s,N):
                     cuspORDSLIST:=[op(cuspORDSLIST),[s,ORD]]:
   		   od:	
     		fi:
		RETURN(cuspORDSLIST):
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 3`);
  fi:
end:
ETA[ETAchanges]:=proc()
printf("**************************************************************\n");
printf("*\n");
printf("*\n");
printf("*  ETA package version 0.3  - Sat, Jun 29, 2019  8:19:25 PM\n");
printf("*  ETA package version 0.2  - Fri, Jun 21, 2019  1:04:35 PM\n");
printf("*  ETA package version 0.1  - Thu Jul 11 13:20:48 EDT 2013\n");
printf("*\n");
printf("*  Current version 0.3 tested on MAPLE 2019\n");
printf("*\n");
printf("*  Changes since version 0.2\n");
printf("*\n");
printf("*     * New function qetacombo                           \n");
printf("*       computes q-expansion of lin combo of eta-products\n");
printf("*\n");
printf("*     * New function  provemodfuncGAMMA0idBATCH           \n");
printf("*       BATCH version of  provemodfuncGAMMA0id           \n");
printf("*\n");
printf("*     * Fixed some minor typos in mintotGAMMA0ORDS       \n");
printf("*     * Fixed minor bug in provemodfuncGAMMA0id       \n");
printf("*     * Fixed minor typo in printETAIDORDStable    \n");
printf("*\n");
printf("*     * Minor changes to cuspord and CUSPORD    \n");
printf("*       cuspord now has 2 args not 3.           \n");
printf("*       USEAGE: cuspord(etaprod,cusp)           \n");
printf("*\n");
printf("*     * Fixed MAJOR BUG in UpLB                 \n");
printf("*       Should use cuspord NOT cuspORD          \n");
printf("*\n");
printf("*     * Fixed MAJOR BUGS in provemodfuncGAMMA0UpETAid \n");
printf("*\n");
printf("*     * New function provemodfuncGAMMA0UpETAidBATCH \n");
printf("*\n");
printf("*  Changes since version 0.1\n");
printf("*\n");
printf("*     * UpLB now accepts nargs = 0           \n");
printf("*\n");
printf("*     * LIST of new functions:\n");
printf("*         ETAchanges \n");
printf("*         ETApversion \n");
printf("*         provemodfuncGAMMA0UpETAid \n");
printf("*\n");
printf("*     *   The function provemodfuncGAMMA0UpETAid \n");
printf("          PROVES the id U[p](EP) = etacombo           \n");
printf("          where EP is an eta-product on Gamma[0](N*p), p is prime,\n");
printf("          etacombo is a linear combination of eta-products on\n");
printf("          Gamma[0](N) where p divides N.         \n");
printf("          USAGE:                                     \n");
printf("          provemodfuncGAMMA0UpETAid()                   \n");
printf("          provemodfuncGAMMA0UpETAid(EP,p,etacombo,N)   \n");
printf("*\n");
printf("*\n");
printf("**************************************************************\n");
RETURN():
end:
ETA[etaCOF]:=proc(EP)
   RETURN(1):
end:
ETA[etaCONSTANT]:=proc()
  local epterm,cc:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("etaCONSTANT(epterm)                                          \n");
     printf("   epterm = c*eta(t[1]*tau)^r[1]* ... *eta(t[s]*tau)^r[s]    \n");
     printf("            ie. an eta-quotient times a constant. \n"); 
     printf("   Returns the constant c.        \n");
     printf("                                                             \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 1 then
        epterm:=args[1]: 
        if type(epterm,constant) then
           RETURN(epterm):
        else
           cc:=eval(subs(eta=etaCOF,epterm)):
           RETURN(cc):
        fi:
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 1`);
  fi:
end:
	ETA[etamult]:=proc(M)
		# computes the etamultiplier of eta(M tau)
		# for a given M in SL_2(Z).
		local a,b,c,d,cm,dm:
		a:=M[1,1]:
		b:=M[1,2]:
		c:=M[2,1]:
		d:=M[2,2]:
		cm:=modp(c,2):
		dm:=modp(d,2):	
		if cm=1 then
		   ETA[jactopstar](d,c)*exp(Pi*I*((a+d)*c-b*d*(c^2-1)-3*c)/12):
		   RETURN(%):
		else
		   if dm=0 then 
		     print(` d= `, d);
		     print(`ERROR: d is even `);
		   fi:
		   ETA[jacbotstar](c,d)*exp(Pi*I*((a+d)*c-b*d*(c^2-1)+3*d-3-3*c*d)/12):
		   RETURN(%):
		fi:
	end:
ETA[etanormalid]:=proc()
     local epcombo,K,VINF,N,j,VINF0,N0;
     global xprint:
        if not(type(xprint,boolean)) then
           xprint:=false:
        fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("etanormalid(epcombo)                                          \n");
     printf("   Renormalizes a sum of etaprods by dividing by the term \n");
     printf("   with the lowest v[oo]. \n");
     printf("-------------------------------------------------------------\n");
  elif nargs = 1 and type(args[1],`+`) then
      epcombo:=args[1]:
        K:=1:
        VINF:=vetainf(op(1,epcombo)):
        N:=nops(op(1,epcombo)):
        ## find term with smallest degree 
        ## and smallest nops if there is a tie
        for j from 2 to nops(epcombo) do
           VINF0:=vetainf(op(j,epcombo)):
           N0:=nops(op(j,epcombo)):
           if VINF0<VINF then
              K:=j:
              VINF:=VINF0:
              N:=N0:
           fi:
           if VINF0=VINF and N0<N then
              K:=j:
              N:=N0:
           fi:
        od:
        if xprint then
	   printf("We divide the eta-prod identity by %a.\n",op(K,epcombo));
        fi:
        RETURN(expand(epcombo/op(K,epcombo)));
     else
        ERROR("epcombo not of type + or nargs not 0 or 1");
     fi:
end:
	ETA[etaprodtoqseries]:=proc(etaprod,T)
		# Computes q-series expansion of
		# the given etaproduct to O(q^T)
		#
	#VERSION: 06.21.2019
        #         Corrected problem with negative ord at oo
		local GP,ngp,xm,qinf,i,t,r,vinf,T1,s:
           if etaprod=1 then
                RETURN(1):
           else
		GP:=ETA[GPmake](etaprod):
                ngp:=nops(GP):
                xm:=1:
		qinf:=0:
                vinf:=vetainf(etaprod):
                if vinf<0 then 
                   T1:=T-2*floor(vinf):
                else
                   T1:=T:
                fi:
                for i from 1 to (ngp/2) do
                  t:=GP[2*i-1]:
                  r:=GP[2*i]:
                  xm:=xm*ETA[POWERq](qseries[etaq](q,t,T),r,T1):       
                  qinf:=qinf+r*t/24:
                od:
                if type(qinf,integer) and qinf<0 then
                     s:=expand(q^(qinf)*convert(series(xm,q,T1),polynom)):
 	 	else
		     s:=series(q^(qinf)*xm,q,T1):
                fi:
		RETURN(s):
           fi:
	end:
	ETA[etaprodtoqseries2]:=proc(etaprod,T)
		# Computes q-series expansion of
		# the given etaproduct to O(q^T)
		# This version of etaprodtoqseries
	  	# omits the q^(qinf) term
		#
		local GP,ngp,xm,qinf,i,t,r:
		global xprint:
          if not(type(xprint,boolean)) then
               xprint:=false:
          fi:
           if etaprod=1 then
                RETURN(1):
           else
		GP:=ETA[GPmake](etaprod):
                ngp:=nops(GP):
                xm:=1:
		qinf:=0:
                for i from 1 to (ngp/2) do
                  t:=GP[2*i-1]:
                  r:=GP[2*i]:
		  if xprint then
	            lprint("computing ",[t,r],eta(t*tau)^r);
                  fi:			
                  xm:=xm*ETA[POWERq](qseries[etaq](q,t,T),r,T):       
                  xm:=series(xm,q,T):
                  if xprint then
                    lprint(i," of ", ngp/2," series computed"):
                  fi:
                  qinf:=qinf+r*t/24:
                od:
	        if xprint then 
                  lprint("NOTE term ",q^(qinf), "is missing"); 
		fi:
		series(xm,q,T):
		RETURN(%):
           fi:
	end:
	ETA[etaprodtoqseriesMODP]:=proc(etaprod,T,p)
		# Computes q-series expansion mod p of
		# the given etaproduct to O(q^T)
		#
		local GP,ngp,xm,qinf,i,t,r:
		GP:=ETA[GPmake](etaprod):
                ngp:=nops(GP):
                xm:=1:
		qinf:=0:
                for i from 1 to (ngp/2) do
                  t:=GP[2*i-1]:
                  r:=GP[2*i]:
                  xm:=xm*ETA[POWERqMODP](qseries[etaq](q,t,T),r,T,p):       
                  xm:=modp(series(xm,q,T),p):
                  qinf:=qinf+r*t/24:
                od:
		series(q^(qinf)*xm,q,T):
                modp(%,p):
		RETURN(%):
	end:
	ETA[etaprodWe]:=proc(etaprod,W,ee,N)
		# Computes that action of W_e on the given etaproduct
		# where  W_e is an
		# Atkin-Lehner involution of GAMMA_0(N)
		local GP,ngp,xm,i,t,r:
		GP:=ETA[GPmake](etaprod):
		ngp:=nops(GP):
		xm:=1:
		for i from 1 to (ngp/2) do
		  t:=GP[2*i-1]:
		  r:=GP[2*i]:
		  xm:=xm*ETA[etaWe](t,W,ee,N)^r:
		od:
		RETURN(simplify(xm)):
	end:
ETA[ETApversion]:=proc()
printf("****************************************************\n");
printf("*\n");
printf("*CURRENT VERSION: \n");
printf("*  ETA   package version 0.3  \n");
printf("*  Sat, Jun 29, 2019  8:19:25 PM\n");
printf("*  This version tested on MAPLE 2019\n"); 
printf("*\n");
printf("*PREVIOUS VERSIONS: \n");
printf("*  ETA   package version 0.2 Fri, Jun 21, 2019  1:04:35 PM\n");
printf("*  misc  package version 0.1 Thu Jul 11 13:20:48 EDT 2013\n");
printf("*\n");
printf("*  Please report any problems to fgarvan@ufl.edu\n");
printf("        \n");
printf("****************************************************\n");
RETURN():
end:
	ETA[etaWe]:=proc(t,W,ee,N)
		# Computes eta(t W_e tau)^r where W_e is an
		# Atkin-Lehner involution of GAMMA_0(N)
		local detW,a,b,c,d,M,delta,detM:
		detW:=linalg[det](W):
		if detW=ee then
		   print(`det W = `, ee);
		else
		   print(`det W is not `, ee);
		   print(`det W = `, detW);
		   print(` W is not an Atkin-Lehner involution`);
		fi:
		if modp(N,ee)=0 and ee<N then 
		   print(` e||N `);
		else
		   print(` e is not a proper divisor of N `);
		   print(` W is not an Atkin-Lehner involution`);
		fi:
		b:=W[1,2]:
		a:=W[1,1]/ee:
		c:=W[2,1]/N:
		d:=W[2,2]/ee:
		if modp(N,t)=0 then
		   print(` t|N `);
		else
		   print(` t= `,t,` N= `,N);
		   print(`ERROR:  t does not divide N`);
		fi:
		delta:=gcd(ee,t):
		M:=array(1..2,1..2,[[a*delta, b*t/delta],[c*N*delta/ee/t,d*ee/delta]]):
		print(` M= `, M );
                detM:=linalg[det](M):
		if detM<>1 then
		   print(` ERROR: M is not in SL_2(Z) `);
		fi:
		ETA[etamult](M)*( (c*N*tau+d*ee)/delta )^(1/2)*eta(ee*t/delta^2*tau):
		RETURN(%):
	end:
	ETA[fanwidth]:=proc(rat,N)
	  local s,fw:
	  s:=denom(rat):
          fw:=N/igcd(N,s^2):
	end:
	ETA[Ffind]:=proc(etaprod,N)
                # Returns 1 if etaprod is in F ie. only
                # negative valence at ioo = 1/N
		# (GAMMA_0(N))
		local S,FS,find,cusp,ord:
                S:=ETA[cuspmake](N):	
                FS:= S minus {1/N}:
                find:=1:
                for cusp in FS do
                    ord:=ETA[cuspord](etaprod,N,cusp):
                    if ord < 0 then
                       find:=0:
                    fi:
                od:
		RETURN(find):
        end:

	ETA[Fricke]:=proc(GP,N)
	  #Fricke involution tau->-1/(N*tau)
	  local NM,NGP,i,m,r:
          NM:=[seq(GP[2*i-1]/N,i=1..nops(GP)/2)]:
	  NGP:=[]:
          for i from 1 to nops(GP)/2 do
             m:=GP[2*i-1]/N:
             r:=GP[2*i]:
             NGP:=[op(NGP),m,r]:
          od:
	  RETURN(NGP):
	end:
	ETA[gamma0FORMCHECK]:=proc(GP,N)
	##gammacheckk:=proc(GP,N)
		# This proc checks whether the eta-function
                # with generalised permutation GP is invariant under
		# GAMMA_0(N) (via Newman's Theorem)
		# Here GP=[t1,r1,t2,r2,...]
		# 
		local ngp,w1,w2,w3,fail,p,c4,i,r,t,S,s,FL,k:
                FL:=[]:
		ngp:=nops(GP):
		k:=add(GP[2*i],i=1..(ngp/2))/2:
		fail:=0:
		##
		## Condition (1): k in Z.
		##
		if type(k,integer) then
            	  # OK
        	else
 		  fail:=fail+1:
                  FL:=[op(FL),(1)]:
		fi:
		w2:=sum(GP[2*i-1]*GP[2*i],i=1..(ngp/2)):
		##
		## Condition (2): w2 = 0 mod 24.
		##
		if modp(w2,24)=0 then
		   ##print(`Condition (2) holds`);
	        else
		   ##print(`Condition (2) does not hold`);
		   fail:=fail+1:
                   FL:=[op(FL),(2)]:
		fi:
	        s:=product(GP[2*i-1]^abs(GP[2*i]),i=1..(ngp/2)):
##		if type(eval(sqrt(p)),integer)=true then
##		   ##print(`Condition (3) holds`);
##	        else
##		   ##print(`Condition (3) does not hold`);
##		   fail:=1:
##		fi:
		##
		## Condition (3): Each t is a divisor of N.
		##
		c4:=1:
		for i from 1 to (ngp/2) do
			r:=GP[2*i]:
                        t:=GP[2*i-1]:
		  if r<>0 and modp(N,t)<>0 then
		    c4:=0:
		  fi:
		od:
		if c4=1 then
		   ##print(`Condition (4) holds`);
	        else
		   ##print(`Condition (4) does not hold`);
		   fail:=fail+1:
                   FL:=[op(FL),(3)]:
		fi:
		##
		## Condition (4): w3 = sum N*r/t = 0 mod 24.  
		##
		w3:=sum((N/GP[2*'i'-1])*GP[2*'i'],'i'=1..(ngp/2)):
		if modp(w3,24)=0 then
		   ##print(`Condition (5) holds`);
	        else
		   ##print(`Condition (5) does not hold`);
		   fail:=fail+1:
                   FL:=[op(FL),(4)]:
		fi:
		if fail=0 then
		   ##print(`function is invariant`);
                   ["N=",N,"weight=",k,"character=",[(-1)^k*ifactor(s),d]]:
		else
                   ["NOT a modular form","Conditions that failed=",FL]:
		fi:
		end:
	ETA[gammacheck]:=proc(GP,N)
		# This proc checks whether the eta-function
                # with generalised permutation GP is invariant under
		# GAMMA_0(N) (via Newman's Theorem)
		# Here GP=[t1,r1,t2,r2,...]
		# 
		local ngp,w1,w2,w3,fail,p,c4,i,r,t,S,s,pp:
		global xprint:
        if not(type(xprint,boolean)) then
           xprint:=false:
        fi:
		ngp:=nops(GP):
		w1:=sum(GP[2*i],i=1..(ngp/2)):
		fail:=0:
		if w1=0 then
                   if xprint then
		   print(`Condition (1) holds`);
		   fi:
	        else
                   if xprint then
		   print(`Condition (1) does not hold`);
		   fi:
		   fail:=1:
		fi:
		w2:=sum(GP[2*i-1]*GP[2*i],i=1..(ngp/2)):
		if modp(w2,24)=0 then
                   if xprint then
		   print(`Condition (2) holds`);
		   fi:
	        else
                   if xprint then
		   print(`Condition (2) does not hold`);
		   fi:
		   fail:=1:
		fi:
	        p:=product(GP[2*i-1]^abs(GP[2*i]),i=1..(ngp/2)):
                pp:=sqrt(p):
		if type(pp,integer)=true then
                   if xprint then
		   print(`Condition (3) holds`);
		   fi:
	        else
                   if xprint then
		   print(`Condition (3) does not hold`);
		   fi:
		   fail:=1:
		fi:
		c4:=1:
		for i from 1 to (ngp/2) do
			r:=GP[2*i]:
                        t:=GP[2*i-1]:
		  if r<>0 and modp(N,t)<>0 then
		    c4:=0:
		  fi:
		od:
		if c4=1 then
                   if xprint then
		   print(`Condition (4) holds`);
		   fi:
	        else
                   if xprint then
		   print(`Condition (4) does not hold`);
		   fi:
		   fail:=1:
		fi:
		w3:=sum((N/GP[2*'i'-1])*GP[2*'i'],'i'=1..(ngp/2)):
		if modp(w3,24)=0 then
                   if xprint then
		   print(`Condition (5) holds`);
		   fi:
	        else
                   if xprint then
		   print(`Condition (5) does not hold`);
		   fi:
		   fail:=1:
		fi:
		if fail=0 then
                   if xprint then
		   print(`function is invariant`);
		   fi:
                   RETURN(1):
		else
                   if xprint then
		   print(`function is NOT invariant`);
		   fi:
                   RETURN(0):
		fi:
		end:
	ETA[gammacheckM]:=proc(GP,N)
	##gammacheckk:=proc(GP,N)
		# This proc checks whether the eta-function
                # with generalised permutation GP is invariant under
		# GAMMA_0(N) (via Newman's Theorem)
		# Here GP=[t1,r1,t2,r2,...]
		# 
		local ngp,w1,w2,w3,fail,p,c4,i,r,t,S,s:
		ngp:=nops(GP):
		w1:=sum(GP[2*i],i=1..(ngp/2)):
		fail:=0:
##		if w1=0 then
##		   ##print(`Condition (1) holds`);
##	        else
##		   ##print(`Condition (1) does not hold`);
##		   fail:=1:
##		fi:
		w2:=sum(GP[2*i-1]*GP[2*i],i=1..(ngp/2)):
		if modp(w2,24)=0 then
		   ##print(`Condition (2) holds`);
	        else
		   ##print(`Condition (2) does not hold`);
		   fail:=1:
		fi:
	        p:=product(GP[2*i-1]^abs(GP[2*i]),i=1..(ngp/2)):
		if type(eval(sqrt(p)),integer)=true then
		   ##print(`Condition (3) holds`);
	        else
		   ##print(`Condition (3) does not hold`);
		   fail:=1:
		fi:
		c4:=1:
		for i from 1 to (ngp/2) do
			r:=GP[2*i]:
                        t:=GP[2*i-1]:
		  if r<>0 and modp(N,t)<>0 then
		    c4:=0:
		  fi:
		od:
		if c4=1 then
		   ##print(`Condition (4) holds`);
	        else
		   ##print(`Condition (4) does not hold`);
		   fail:=1:
		fi:
		w3:=sum((N/GP[2*'i'-1])*GP[2*'i'],'i'=1..(ngp/2)):
		if modp(w3,24)=0 then
		   ##print(`Condition (5) holds`);
	        else
		   ##print(`Condition (5) does not hold`);
		   fail:=1:
		fi:
		if fail=0 then
		   ##print(`function is invariant`);
		else
		   ##print(`function is NOT invariant`);
		fi:
		RETURN(modp(fail+1,2)):
		end:
	ETA[gp2etaprod]:=proc(gp)
		# Converts a gen perm into an eta-product
		local ngp,xm,i,t,r:
                ngp:=nops(gp):
                xm:=1:
                for i from 1 to (ngp/2) do
                  t:=gp[2*i-1]:
                  r:=gp[2*i]:
                  xm:=xm*eta(tau*t)^r:       
                od:
		RETURN(xm):
	end:
	ETA[GPmake]:=proc(etaprod)
		# This proc finds the GP corresponding to the given etaproduct
		# 06/21/00: bug fix. Now works for single eta function
		#           like eta(40*tau)^12
		local L1,L1n,GP,i,r,p,t:
	        if whattype(etaprod)=`^` then
		   L1:=[etaprod]:
	     	else
		   if whattype(etaprod)=function then
			RETURN([op(etaprod)/tau,1]);
		   else
		        L1:=[op(etaprod)]:
	           fi:
		fi:
		L1n:=nops(L1):
		GP:=NULL:
		for i from 1 to L1n do
                    r:=degree(L1[i]):
		    if r=1 then
                       p:=op(L1[i]):
	               t:=p/tau:
		    else
                       p:=op(L1[i]):
	               t:=op(p[1])/tau:
		    fi:
		    GP:=GP,t,r:
		od:
		RETURN([GP]);
	end:	
	ETA[jacbotstar]:=proc(c,d)
		local dd,m,mm:
		dd:=abs(d):
		if c=0 then
		   if dd=1 then
		      RETURN(1):
		   fi:
		else
		   m:=(csgn(c)-1)/2*(csgn(d)-1)/2:
                   mm:=(-1)^m:
          	   RETURN(mm*numtheory[jacobi](c,dd)):
		fi:
	end:
	ETA[jactopstar]:=proc(c,d)
		local dd:
		dd:=abs(d):
		if c=0 then
		   if dd=1 then
		      RETURN(1):
		   fi:
		else
          	   RETURN(numtheory[jacobi](c,dd)):
		fi:
	end:
ETA[mintotGAMMA0ORDS]:=proc()
# 09/06/09: version of mintotORD2 (mintotORD2 does not work if L in an array)
# 06/21/19: num should be N                                                      
local nL,numORDS,L,num,N,i,j;
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 = (c1 f1 + c2 f2 + ...)
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("mintotORDS(L,N)                                            \n");
     printf("   L is a list (or array) [[ORDS(f1), ORDS(f2), ...].        \n");
     printf("   N   = nops(L) (list) or dim(L) (array)                    \n"); 
     printf("   This proc calculates a lower bound for sum ORD g where    \n");
     printf("   g = (c1 f1 + c2 f2 + ...)                             \n");
     printf("   Here sum is over cusps not equivalent to oo.              \n");
     printf("   (Usually) each ORDS(f) was produced by cuspORDSnotoo   \n");
     printf("   (Usually) 1 is one term of g (but not necessarily)     \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 2 then
        L:=args[1]: num:=args[2]: 
	nL:=num:
	numORDS:=nops(L[1]):
        ## cusp oo = 1/N has been removed.
        ##for i from 1 to numORDS do
        ##print(i,seq(L[j][i][2],j=1..nL),"*",min(seq(L[j][i][2],j=1..nL)));
        ##od;
        RETURN(add(min(seq(L[j][i][2],j=1..nL)),i=1..numORDS));
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 2`);
  fi:
end:
ETA[POWERPq]:=proc(a,N,M)  #compute a^N to O(q^M) (where N>=0)
	local j,Q,Y,n,b:
        global xprint:
	j:=1:
	Q:=1:
	Y:=1:
	n:=N:
        if not(type(xprint,boolean)) then
           xprint:=false:
        fi:
	while n>0 do
if xprint then
print("n=",n);
fi:
	  if modp(n,2)=1 then
	    b:=1:
	    if n=N then
		Q:=series(a,q,M): Y:=Q:                  
	    else
	        Q:=series(Q^2,q,M): Y:=series(Q*Y,q,M):
	    fi:
	  else 
	    b:=0:
	    if n=N then
		Q:=series(a,q,M):
	    else
	        Q:=series(Q^2,q,M): 
	    fi:
	  fi:
	  n:=trunc(n/2):
	od:
if xprint then
print("POWERPq done.");
fi:
	RETURN(Y);
end:
ETA[POWERPqMODP]:=proc(a,N,M,p)  #compute a^n to O(q^M) mod p
                                 #here a>=0
	local j,Q,Y,n,b:
	j:=1:
	Q:=1:
	Y:=1:
	n:=N:
	while n>0 do
print("N=",N,"n=",n);
	  if modp(n,2)=1 then
	    b:=1:
	    if n=N then
		Q:=modp(series(a,q,M),p): Y:=Q:
	    else
	        Q:=modp(series(Q^2,q,M),p): Y:=modp(series(Q*Y,q,M),p):
	    fi:
	  else 
	    b:=0:
	    if n=N then
		Q:=modp(series(a,q,M),p): 
	    else
	        Q:=modp(series(Q^2,q,M),p): 
	    fi:
	  fi:
	  n:=trunc(n/2):
	od:
	RETURN(Y);
end:
ETA[POWERq]:=proc(a,N,M)  #compute a^n to O(q^M)
	local mN:
	if N>=0 then
           RETURN(ETA[POWERPq](a,N,M))
        else
           mN:=-N:
           RETURN(ETA[POWERPq](1/a,mN,M))
        fi:
end:
ETA[POWERqMODP]:=proc(a,N,M,p)  #compute a^n to O(q^M) mod p
	local mN:
	if N>=0 then
           RETURN(ETA[POWERPqMODP](a,N,M,p))
        else
           mN:=-N:
           RETURN(ETA[POWERPqMODP](1/a,mN,M,p))
        fi:
end:
	ETA[printcuspords]:=proc(etaprod,N)
		# Computes the order at each cusp (GAMMA_0(N))
		# of the given etaproduct and prints result.
		#
		local GP,ngp,S,s,ords,i,t,r,c:
		GP:=ETA[GPmake](etaprod):
		ngp:=nops(GP):
		S:=ETA[cuspmake](N):
		for s in S do
		  ords:=0:
		  for i from 1 to (ngp/2) do
		     t:=GP[2*i-1]:
		     r:=GP[2*i]:
		     c:=denom(s):
		     ords:=ords+gcd(t,c)^2/t*r/24:
		  od:
		  print(` Order at cusp `,s, ` is `,ords);
		od:	
		RETURN():
	end:
	ETA[printcuspORDS]:=proc(etaprod,N)
		# Computes the order at each cusp (GAMMA_0(N))
		# of the given etaproduct and prints result.
		#
		local GP,ngp,S,s,ords,i,t,r,c, TOTORD,ORDS :
		GP:=ETA[GPmake](etaprod):
		ngp:=nops(GP):
		S:=ETA[cuspmake](N):
                TOTORD:=0:
		for s in S do
		  ords:=0:
		  for i from 1 to (ngp/2) do
		     t:=GP[2*i-1]:
		     r:=GP[2*i]:
		     c:=denom(s):
		     ords:=ords+gcd(t,c)^2/t*r/24:
		  od:
                  ORDS:=ords*fanwidth(s,N):
                  TOTORD:=TOTORD+ORDS:
                  print(`cusp =`,s,` order =`,ords,` ORD=`,ORDS);
		od:	
                print(`TOT ORD = `,TOTORD);
		RETURN():
	end:
 
ETA[printETAIDORDStable]:=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("printETAIDORDStable()                                        \n");
     printf("   Print a table of ORDS for each term in a eta-prod-identity \n");
     printf("   using global data produced by the function \n");
     printf("   provemodfuncGAMMA0id.\n");
     printf("   Table is stored in the matrix bigmat which is returned.   \n");
     printf("-------------------------------------------------------------\n");
        numfuncs:= ArrayNumElems(_ORDS): ## added 7/26/13
        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-1]): 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][2]:
            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(consL[j]*_F[j],j=1..numfuncs)=0);
        printf("where\n");
        for k from 1 to numfuncs do
           printf("_F[%a] = \n",k-1);
           print(etaPRODL[k]);
        od:
        sumcolG:=add(bigmat[j,numfuncs+2],j=2..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 we see that \n");
        printf("the identity can be proved by showing that the coefficients\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 provemodfuncGAMMA0id.\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:
#++++++++++++++++++++++++++++++++++++++++++++++++++++
ETA[provemodfuncGAMMA0id]:=proc()
local CW;
local etaid,N,CUSPS,CUSPS1,nn,j,XX,cc,ordtmp,mftmp,num2;
local num3,proofq,ff1,TT,qsj,ss;
# 06.26.19 -mintottmp --> -mintottmp+1
    global qcheck,modfunccheck, totcheck, _ORDS,jptmp,jpqd,eptmp,gltmp:
    global etaPRODL,GPL,COFS,conpres,CONTERMS,mintottmp,consL,MFLB:
    ##
    ## proving jacid as a modfunc on Gamma1(N)
    ##
 global xprint,_CUSPS:
        if not(type(xprint,boolean)) then
          xprint:=false:
        fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("provemodfuncGAMMA0id(etaid,N)                                \n");
     printf("   etaid = sum of modular functions on Gamma[0](N)           \n");
     printf("           Each term in the sum is a eta-quotient to base N. \n"); 
     printf("   CUSPS = Set of inequivalent cusps for Gamma[0](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]:
        etaid:=args[1]: N:=args[2]:
        CUSPS:=convert(cuspmake(N),list):
        CUSPS1:=cuspmake(N) minus {1/N}:
        CUSPS1:=convert(CUSPS1,list):
        _CUSPS:=CUSPS1:
    conpres:=0:  #conpres=0 if no terms are constants
    qcheck:=[]:
    modfunccheck:=[]:
    totcheck:=[]:
    CONTERMS:=[]:
    nn:=nops(etaid):
    _ORDS:=Array(1..nn):
    etaPRODL:=array(1..nn):
    consL:=array(1..nn):
    GPL:=array(1..nn):
    COFS:=array(1..nn):
    for j from 1 to nn do
      lprint("TERM ",j,"of ",nn," *****************");
      XX:=op(j,etaid):
      if xprint then
        print("XX=",XX);
      fi:
      cc:=eval(subs(eta=etaCOF,XX)):
      COFS[j]:=cc:
      eptmp:=eval(XX/cc):
      if eptmp=1 then
          etaPRODL[j]:=1:
          consL[j]:=cc: 
          GPL[j]:=[]:
          ordtmp:=cuspORDSnotoo(eptmp,CUSPS1,N):
          _ORDS[j]:=ordtmp:
          conpres:=1:
          CONTERMS:=[op(CONTERMS),j]:
      else
          etaPRODL[j]:=eptmp:
          gltmp:=GPmake(eptmp):            
          consL[j]:=cc:
          GPL[j]:=gltmp:
          cuspORDS(eptmp,CUSPS,N):
          ordtmp:=cuspORDSnotoo(eptmp,CUSPS1,N):
          _ORDS[j]:=ordtmp:
          if toterror<>0 then
            lprint("WARNING: totcheck: ");
            totcheck:=[op(totcheck),j]:
          fi:
          mftmp:=gammacheck(gltmp,N):
          if mftmp<>1 then
            lprint("WARNING: modfunccheck: ");
            modfunccheck:=[op(modfunccheck),j]:
          fi:
      fi:
   od:
   mintottmp:=mintotGAMMA0ORDS(_ORDS,nn): 
   lprint("mintotord = ",mintottmp);
   lprint("TO PROVE the identity we need to show that v[oo](ID) > ",-mintottmp);
   MFLB:=-mintottmp:
   num2:=nops(modfunccheck):
   num3:=nops(totcheck):
   if 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("      Gamma0(%a). \n",N);
     printf("*** o We also checked that the total order of\n",N);
     printf("      each term was zero.\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+1);
   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:=0:
         TT:=-mintottmp+2*N:
         for j from 1 to nn do
             if etaPRODL[j]<>1 then
                qsj:=etaprodtoqseries(etaPRODL[j],TT):
                qsj:=convert(qsj,polynom):
                ff1 := ff1 + consL[j]*qsj:
              else
                ff1 := ff1 + consL[j]*1:                
             fi:
         od:
         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();
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 2`);
  fi:
end:
ETA[provemodfuncGAMMA0idBATCH]:=proc()
local CW;
local etaid,N,CUSPS,CUSPS1,nn,j,XX,cc,ordtmp,mftmp,num2;
local num3,proofq,ff1,TT,qsj,ss,symid,pp,ss0;
# 06.26.19 -mintottmp --> -mintottmp+1
    global qcheck,modfunccheck, totcheck, _ORDS,jptmp,jpqd,eptmp,gltmp:
    global etaPRODL,GPL,COFS,conpres,CONTERMS,mintottmp,consL,MFLB:
    ##
    ## proving etaid as a modfunc on Gamma0(N)
    ##
 global xprint, overrideproofq,_CUSPS,noprint:
 global qthreshold,proveit,printlocalsymid:
        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("provemodfuncGAMMA0idBATCH(symid,etaid,N)                     \n");
     printf("provemodfuncGAMMA0idBATCH(etaid,N)                           \n");
     printf("   This is a BATCH version of provemodfuncGAMMA0id \n");
     printf("   It uses the global var qthreshold .\n");
     printf("   symid = symbolic form of the identity. \n");
     printf("   etaid = sum of modular functions on Gamma[0](N)           \n");
     printf("           Each term in the sum is a eta-quotient to base N. \n"); 
     printf("   CUSPS = Set of inequivalent cusps for Gamma[0](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("-------------------------------------------------------------\n");
     elif nargs = 2 or nargs = 3 then
        if nargs = 2 then
           etaid:=args[1]: N:=args[2]:
           symid:=etaid=0:
        else
           symid:=args[1]: etaid:=args[2]: N:=args[3]:
        fi:
        CUSPS:=convert(cuspmake(N),list):
        CUSPS1:=cuspmake(N) minus {1/N}:
        CUSPS1:=convert(CUSPS1,list):
        _CUSPS:=CUSPS1:
    conpres:=0:  #conpres=0 if no terms are constants
    qcheck:=[]:
    modfunccheck:=[]:
    totcheck:=[]:
    CONTERMS:=[]:
    nn:=nops(etaid):
    _ORDS:=Array(1..nn):
    etaPRODL:=array(1..nn):
    consL:=array(1..nn):
    GPL:=array(1..nn):
    COFS:=array(1..nn):
    for j from 1 to nn do
      ##lprint("TERM ",j,"of ",nn," *****************");
      XX:=op(j,etaid):
      if xprint then
        print("XX=",XX);
      fi:
      cc:=eval(subs(eta=etaCOF,XX)):
      COFS[j]:=cc:
      eptmp:=eval(XX/cc):
      if eptmp=1 then
          etaPRODL[j]:=1:
          consL[j]:=cc: 
          GPL[j]:=[]:
          ordtmp:=cuspORDSnotoo(eptmp,CUSPS1,N):
          _ORDS[j]:=ordtmp:
          conpres:=1:
          CONTERMS:=[op(CONTERMS),j]:
      else
          etaPRODL[j]:=eptmp:
          gltmp:=GPmake(eptmp):            
          consL[j]:=cc:
          GPL[j]:=gltmp:
          cuspORDS(eptmp,CUSPS,N):
          ordtmp:=cuspORDSnotoo(eptmp,CUSPS1,N):
          _ORDS[j]:=ordtmp:
          if toterror<>0 then
            lprint("WARNING: totcheck: ");
            totcheck:=[op(totcheck),j]:
          fi:
          mftmp:=gammacheck(gltmp,N):
          if mftmp<>1 then
            lprint("WARNING: modfunccheck: ");
            modfunccheck:=[op(modfunccheck),j]:
          fi:
      fi:
   od:
   mintottmp:=mintotGAMMA0ORDS(_ORDS,nn): 
   ##lprint("mintotord = ",mintottmp);
   ##lprint("TO PROVE the identity we need to show that v[oo](ID) > ",-mintottmp);
   MFLB:=-mintottmp:
   num2:=nops(modfunccheck):
   num3:=nops(totcheck):
   if num2+num3>0 then
     lprint("***********************************");
     lprint("See qcheck, modfunccheck, totcheck.");
     error("There were errors.");
   else
     if noprint=false then
     printf("*** There were NO errors. \n");
     printf("*** o Each term was modular function on\n");
     printf("      Gamma0(%a). \n",N);
     printf("*** o We also checked that the total order of\n",N);
     printf("      each term was zero.\n");
     fi:
   fi:
   if conpres=1 and xprint then
     lprint("*** WARNING: some terms were constants. ***");
     lprint("See array CONTERMS.");
   fi:                    
#HERE
   if noprint=false then
     printf("To prove the identity we will need to verify if up to \n");
     printf("q^(%a).\n",-mintottmp+1);
   fi:
         ff1:=0:
         TT:=-mintottmp+2*N:
         for j from 1 to nn do
             if etaPRODL[j]<>1 then
                qsj:=etaprodtoqseries(etaPRODL[j],TT):
                qsj:=convert(qsj,polynom):
                ff1 := ff1 + consL[j]*qsj:
              else
                ff1 := ff1 + consL[j]*1:                
             fi:
         od:
         ss0:=series(ff1,q,-mintottmp+2*N);
         ss:=simplify(ss0);
         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,mintottmp,ss0]):
   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,mintottmp]):
   fi:
   ## 06.26.16
   if printlocalsymid then
     print(symid);
   fi:
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0,2 or 3`);
  fi:
end:
##OLD         if ss=0 then
##OLD           printf("RESULT: The identity holds to O(q^(\%a)).\n",-mintottmp+2*N);
##OLD           printf("CONCLUSION: This proves the identity since we had only\n");
##OLD           printf("            to show that v[oo](ID) > %a.\n",MFLB);  
##OLD         else 
##OLD           printf("WARNING: ss <> 0. Something is wrong.\n");
##OLD           printf("ss = %a. Something is wrong.\n",ss);
##OLD         fi:
##OLD      else
##OLD         printf("You did not enter yes.\n");
##OLD      fi:
##OLD    RETURN();
##OLD  else
##OLD     printf("nargs=%a\n",nargs);
##OLD     ERROR(`nargs must be 0 or 2`);
##OLD  fi
ETA[ provemodfuncGAMMA0UpETAid]:=proc()
local CW;
local etaid,N,CUSPS,CUSPS1,nn,j,XX,cc,ordtmp,mftmp,num2;
local num3,proofq,ff1,TT,qsj,ss;
local EP,p,etacombo,n,ff0,sff0,chkEP,vinf,r,ss0;
    global qcheck,modfunccheck, totcheck, _ORDS,jptmp,jpqd,eptmp,gltmp:
    global etaPRODL,GPL,COFS,conpres,CONTERMS,mintottmp,consL,MFLB:
    ##
    ## proving jacid as a modfunc on Gamma1(N)
    ##
 global xprint,_CUSPS,UpORDL,_ORDS2:
    ## CREATED 06.21.2019
        # _CUSPS list of cusps except oo
        # UpORDL list of cusp ORDS of U[p](EP)
        # _ORDS2 list of cusp ORDS of all funcs involved
        if not(type(xprint,boolean)) then
          xprint:=false:
        fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("provemodfuncGAMMA0UpETAid(EP,p,etacombo,N)                   \n");
     printf("      EP = one eta-product                                   \n");
     printf("       p = prime                                             \n");
     printf("etacombo = sum of modular functions on Gamma[0](N)           \n");
     printf("           Each term in the sum is a eta-quotient to base N. \n"); 
     printf("       N = Positive integer multiple of p                    \n");
     printf("                                                             \n");
     printf("   This function PROVES the id U[p](EP) = etacombo           \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("-------------------------------------------------------------\n");
     elif nargs = 4 then
     ###############################################################
     # 06.27.19
     EP:=args[1]: p:=args[2]: etacombo:=args[3]: N:=args[4]:
     ###############################################################
     # first we check that EP is a MF in Gamma[0](N*p)
     ###############################################################
     # 06.27.19
     if modp(N,p)<>0 then
       ERROR("N not divisible by p");
     fi:
     ###############################################################
     chkEP:=gammacheck(GPmake(EP),p*N):
     if chkEP = 1 then
        ##jacid:=args[1]: CUSPS:=args[2]: WIDS:=args[3]: N:=args[4]:
        CUSPS:=convert(cuspmake(N),list):
        CUSPS1:=cuspmake(N) minus {1/N}:
        CUSPS1:=convert(CUSPS1,list):
        _CUSPS:=CUSPS1:
    conpres:=0:  #conpres=0 if no terms are constants
    qcheck:=[]:
    modfunccheck:=[]:
    totcheck:=[]:
    CONTERMS:=[]:
     ###############################################################
     # 06.27.19 fixed nn=1 problem
    if type(etacombo,`*`) then
      nn:=1:
    else
      nn:=nops(etacombo):
    fi:
     ###############################################################
    _ORDS:=Array(1..nn):
    _ORDS2:=Array(1..nn+1):
    etaPRODL:=array(1..nn):
    consL:=array(1..nn):
    GPL:=array(1..nn):
    COFS:=array(1..nn):
    for j from 1 to nn do
      if nn=1 then
        XX:=etacombo:
      else
        XX:=op(j,etacombo):
      fi:
      if xprint then
        lprint("TERM ",j,"of ",nn," *****************");
        print("XX=",XX);
      fi:
      cc:=eval(subs(eta=etaCOF,XX)):
      COFS[j]:=cc:
      eptmp:=eval(XX/cc):
      if eptmp=1 then
          etaPRODL[j]:=1:
          consL[j]:=cc: 
          GPL[j]:=[]:
          ordtmp:=cuspORDSnotoo(eptmp,CUSPS1,N):
          _ORDS[j]:=ordtmp:
          conpres:=1:
          CONTERMS:=[op(CONTERMS),j]:
      else
          etaPRODL[j]:=eptmp:
          gltmp:=GPmake(eptmp):            
          consL[j]:=cc:
          GPL[j]:=gltmp:
          cuspORDS(eptmp,CUSPS,N):
          ordtmp:=cuspORDSnotoo(eptmp,CUSPS1,N):
          _ORDS[j]:=ordtmp:
          if toterror<>0 then
            if xprint then
              lprint("WARNING: totcheck: ");
            fi:
            totcheck:=[op(totcheck),j]:
          fi:
          mftmp:=gammacheck(gltmp,N):
          if mftmp<>1 then
            if xprint then
              lprint("WARNING: modfunccheck: ");
            fi:
            modfunccheck:=[op(modfunccheck),j]:
          fi:
      fi:
   od:
   ## NOTE: _ORDS is list of cuspORDS of funcs in etacombo
   ## Now compute lower bound for each ORD(U[p](EP,cusp,N)
   ######################################################################
   ## UpORDL:=[seq([r,ceil(UpLB(EP,r,N,p)*fanwidth(r,N))],r in CUSPS1)]:
   # line above corrected 06.27.2019
   UpORDL:=[seq([r,ceil(UpLB(EP,r,N*p,p))],r in CUSPS1)]:
   ######################################################################
   _ORDS2[1]:=UpORDL:            
   for n from 2 to nn+1 do _ORDS2[n]:=_ORDS[n-1]: od:
   ##print("_ORDS2=",_ORDS2);
   mintottmp:=mintotGAMMA0ORDS(_ORDS2,nn): 
   ##lprint("mintotord = ",mintottmp);
   ##lprint("TO PROVE the identity we need to show that v[oo](ID) > ",-mintottmp);
   MFLB:=-mintottmp:
   num2:=nops(modfunccheck):
   num3:=nops(totcheck):
   if num2+num3>0 then
     lprint("***********************************");
     lprint("See qcheck, modfunccheck, totcheck.");
     error("There were errors.");
   else
     printf("*** There were NO errors. \n");
     printf("*** o EP is an MF on Gamma[0](%a)\n",N*p);
     printf("*** o Each term in the etacombo is a  modular function on\n");
     printf("      Gamma0(%a). \n",N);
     printf("*** o We also checked that the total order of\n",N);
     printf("      each term etacombo was zero.\n");
   fi:
   if conpres=1 then
     lprint("*** WARNING: some terms were constants. ***");
     lprint("See array CONTERMS.");
   fi:                    
     printf("*** To prove the identity U[%a](EP)=etacombo we need to show\n",p);
     printf("    that v[oo](ID) > %a",-mintottmp);
     printf("    This means checking up to q^(%a).\n",-mintottmp+1);
   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:=0:
         TT:=-mintottmp+2*N:
         for j from 1 to nn do
             if etaPRODL[j]<>1 then
                qsj:=etaprodtoqseries(etaPRODL[j],TT):
                qsj:=convert(qsj,polynom):
                ff1 := ff1 + consL[j]*qsj:
              else
                ff1 := ff1 + consL[j]*1:                
             fi:
         od:
         ## now calculate q-expansion of left side
	 vinf:=vetainf(EP):
         if vinf<0 then
            ff0:=etaprodtoqseries(EP,p*TT-vinf*p):
	 else
            ff0:=etaprodtoqseries(EP,p*TT+p+1):
         fi:
         sff0:=sift(ff0,q,p,0,p*TT+p):
         ss0:=series(sff0-ff1,q,-mintottmp+2*N);
         printf("We find that LHS - RHS is \n");
         print(ss0);
         ss:=simplify(series(sff0-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();
  else
    ERROR("EP is not a MF on Gamma[0](p*N)");
  fi:
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 4`);
  fi:
end:
ETA[ provemodfuncGAMMA0UpETAidBATCH]:=proc()
local CW;
local etaid,N,CUSPS,CUSPS1,nn,j,XX,cc,ordtmp,mftmp,num2;
local num3,proofq,ff1,TT,qsj,ss;
local EP,p,etacombo,n,ff0,sff0,chkEP,vinf,r,ss0,pp;
    global qcheck,modfunccheck, totcheck, _ORDS,jptmp,jpqd,eptmp,gltmp:
    global etaPRODL,GPL,COFS,conpres,CONTERMS,mintottmp,consL,MFLB:
    ##
    ## proving jacid as a modfunc on Gamma1(N)
    ##
 global xprint,_CUSPS,UpORDL,_ORDS2,noprint:
    ## CREATED 06.21.2019
        # _CUSPS list of cusps except oo
        # UpORDL list of cusp ORDS of U[p](EP)
        # _ORDS2 list of cusp ORDS of all funcs involved
        if not(type(xprint,boolean)) then
          xprint:=false:
        fi:
        if not(type(noprint,boolean)) then
          noprint:=false:  
        fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("provemodfuncGAMMA0UpETAidBATCH(EP,p,etacombo,N)              \n");
     printf("This a BATCH version of provemodfuncGAMMA0UpETAid            \n");
     printf("      EP = one eta-product                                   \n");
     printf("       p = prime                                             \n");
     printf("etacombo = sum of modular functions on Gamma[0](N)           \n");
     printf("           Each term in the sum is a eta-quotient to base N. \n"); 
     printf("       N = Positive integer multiple of p                    \n");
     printf("                                                             \n");
     printf("   This function PROVES the id U[p](EP) = etacombo           \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("-------------------------------------------------------------\n");
     elif nargs = 4 then
        EP:=args[1]: p:=args[2]: etacombo:=args[3]: N:=args[4]:
     # first we check that EP is a MF in Gamma[0](N*p)
     ###############################################################
     # 06.27.19
     if modp(N,p)<>0 then
       ERROR("N not divisible by p");
     fi:
     ###############################################################
     chkEP:=gammacheck(GPmake(EP),p*N):
     if chkEP = 1 then
        ##jacid:=args[1]: CUSPS:=args[2]: WIDS:=args[3]: N:=args[4]:
        CUSPS:=convert(cuspmake(N),list):
        CUSPS1:=cuspmake(N) minus {1/N}:
        CUSPS1:=convert(CUSPS1,list):
        _CUSPS:=CUSPS1:
    conpres:=0:  #conpres=0 if no terms are constants
    qcheck:=[]:
    modfunccheck:=[]:
    totcheck:=[]:
    CONTERMS:=[]:
      ###############################################################
      # 06.27.19 fixed nn=1 problem
     if type(etacombo,`*`) then
       nn:=1:
     else
       nn:=nops(etacombo):
     fi:
      ###############################################################
    _ORDS:=Array(1..nn):
    _ORDS2:=Array(1..nn+1):
    etaPRODL:=array(1..nn):
    consL:=array(1..nn):
    GPL:=array(1..nn):
    COFS:=array(1..nn):
    for j from 1 to nn do
      ###############################################################
       if nn=1 then
         XX:=etacombo:
       else
         XX:=op(j,etacombo):
       fi:
      if xprint then
        lprint("TERM ",j,"of ",nn," *****************");
        print("XX=",XX);
      fi:
      cc:=eval(subs(eta=etaCOF,XX)):
      COFS[j]:=cc:
      eptmp:=eval(XX/cc):
      if eptmp=1 then
          etaPRODL[j]:=1:
          consL[j]:=cc: 
          GPL[j]:=[]:
          ordtmp:=cuspORDSnotoo(eptmp,CUSPS1,N):
          _ORDS[j]:=ordtmp:
          conpres:=1:
          CONTERMS:=[op(CONTERMS),j]:
      else
          etaPRODL[j]:=eptmp:
          gltmp:=GPmake(eptmp):            
          consL[j]:=cc:
          GPL[j]:=gltmp:
          cuspORDS(eptmp,CUSPS,N):
          ordtmp:=cuspORDSnotoo(eptmp,CUSPS1,N):
          _ORDS[j]:=ordtmp:
          if toterror<>0 then
            if xprint then
              lprint("WARNING: totcheck: ");
            fi:
            totcheck:=[op(totcheck),j]:
          fi:
          mftmp:=gammacheck(gltmp,N):
          if mftmp<>1 then
            if xprint then
              lprint("WARNING: modfunccheck: ");
            fi:
            modfunccheck:=[op(modfunccheck),j]:
          fi:
      fi:
   od:
   ## NOTE: _ORDS is list of cuspORDS of funcs in etacombo
   ## Now compute lower bound for each ORD(U[p](EP,cusp,N)
   ######################################################################
   ## UpORDL:=[seq([r,ceil(UpLB(EP,r,N,p)*fanwidth(r,N))],r in CUSPS1)]:
   # line above corrected 06.27.2019
   UpORDL:=[seq([r,ceil(UpLB(EP,r,N*p,p))],r in CUSPS1)]: 
   ######################################################################
   _ORDS2[1]:=UpORDL:            
   for n from 2 to nn+1 do _ORDS2[n]:=_ORDS[n-1]: od:
   ##print("_ORDS2=",_ORDS2);
   mintottmp:=mintotGAMMA0ORDS(_ORDS2,nn): 
   ##lprint("mintotord = ",mintottmp);
   ##lprint("TO PROVE the identity we need to show that v[oo](ID) > ",-mintottmp);
   MFLB:=-mintottmp:
   num2:=nops(modfunccheck):
   num3:=nops(totcheck):
   if num2+num3>0 then
     lprint("***********************************");
     lprint("See qcheck, modfunccheck, totcheck.");
     error("There were errors.");
   else
     if not noprint then
      printf("*** There were NO errors. \n");
      printf("*** o EP is an MF on Gamma[0](%a)\n",N*p);
      printf("*** o Each term in the etacombo is a  modular function on\n");
      printf("      Gamma0(%a). \n",N);
      printf("*** o We also checked that the total order of\n",N);
      printf("      each term etacombo was zero.\n");
    fi:
   fi:
   if conpres=1 and xprint then 
     lprint("*** WARNING: some terms were constants. ***");
     lprint("See array CONTERMS.");
   fi:                    
     if noprint=false then
      printf("*** To prove the identity U[%a](EP)=etacombo we need to show\n",p);
      printf("    that v[oo](ID) > %a",-mintottmp);
      printf("    This means checking up to q^(%a).\n",-mintottmp+1);
     fi:
   ##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:=0:
         TT:=-mintottmp+2*N:
         for j from 1 to nn do
             if etaPRODL[j]<>1 then
                qsj:=etaprodtoqseries(etaPRODL[j],TT):
                qsj:=convert(qsj,polynom):
                ff1 := ff1 + consL[j]*qsj:
              else
                ff1 := ff1 + consL[j]*1:                
             fi:
         od:
         ## now calculate q-expansion of left side
	 vinf:=vetainf(EP):
         if vinf<0 then
            ff0:=etaprodtoqseries(EP,p*TT-vinf*p):
	 else
            ff0:=etaprodtoqseries(EP,p*TT+p+1):
         fi:
         sff0:=sift(ff0,q,p,0,p*TT+p):
         ss0:=series(sff0-ff1,q,-mintottmp+2*N);
         if not noprint then
         printf("We find that LHS - RHS is \n");
         print(ss0);
         fi:
         ss:=simplify(series(sff0-ff1,q,-mintottmp+2*N));
         pp:=convert(ss,polynom):
         if pp=0 then
            RETURN([1,mintottmp,ss0]):
         else 
            RETURN([0,mintottmp]):
         fi:
      ##else
      ##   printf("You did not enter yes.\n");
      ##fi:
  else
    ERROR("EP is not a MF on Gamma[0](p*N)");
  fi:
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 4`);
  fi:
end:
ETA[qetacombo]:=proc()
local etaid,TT,BB,nn,etaPRODL,COFS,j,XX,cc,eptmp,GPL,conpres,gltmp:
local ff1,qsj,ss,etacombo;
global consL:
##X########################################################################
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("qetacombo()                                \n");
     printf("qetacombo(etacombo,TT)                                \n");
     printf("qetacombo(etacombo,TT,BB)                                \n");
     printf("   etacombo = lin combo of eta-products                      \n");
     printf("   TT = positive integer                                     \n");
     printf("   BB = nonnegative integer (default = 0)                    \n"); 
     printf("                                                             \n"); 
     printf("   Compute q-expansions of etacombo up to q^TT   \n"); 
     printf("   Increase BB if not enough terms               \n"); 
     printf("-------------------------------------------------------------\n");
##X########################################################################
  else 
#####X#####################################################################
     if nargs = 2 or nargs = 3 then
        etacombo:=args[1]: TT:=args[2]: 
########X##################################################################
        if nargs=3 then
          BB:=args[3]:
        else
          BB:=0:
        fi:
#####X#####################################################################
      # 06.27.19 fixed nn=1 problem
     if type(etacombo,`*`) then
       nn:=1:
     else
       nn:=nops(etacombo):
     fi:
      ###############################################################
     etaPRODL:=array(1..nn):
     consL:=array(1..nn):
     COFS:=array(1..nn):
     for j from 1 to nn do
      ###############################################################
       if nn=1 then
         XX:=etacombo:
       else
         XX:=op(j,etacombo):
       fi:
       cc:=eval(subs(eta=etaCOF,XX)):
       COFS[j]:=cc:
       eptmp:=eval(XX/cc):
#######X###################################################################
       if eptmp=1 then
           etaPRODL[j]:=1:
           consL[j]:=cc: 
           GPL[j]:=[]:
           conpres:=1:
#######X###################################################################
       else
           etaPRODL[j]:=eptmp:
           consL[j]:=cc:
       fi:
#######X###################################################################
     od:
     ff1:=0:
     for j from 1 to nn do
#########X#################################################################
         if etaPRODL[j]<>1 then
            qsj:=etaprodtoqseries(etaPRODL[j],TT):
            qsj:=convert(qsj,polynom):
            ff1 := ff1 + consL[j]*qsj:
#########X#################################################################
         else
            ff1 := ff1 + consL[j]*1:                
#########X#################################################################
         fi:
      od:
      ss:=simplify(series(ff1,q,TT+BB));
      ss:=convert(ss,polynom):
      RETURN(ss);
   else
      printf("nargs=%a\n",nargs);
      ERROR(`nargs must be 0,2,or 3`);
   fi:
  fi:
##X########################################################################
end:
ETA[UpLB]:=proc()
###
### See
### math/reports/referee/bkim/revised-paper/prog1a
###
### 06.27.19: fixed major bug (use cuspord NOT cuspORD)
    local N,d,LB,EP,r,N1,p,b:
    global xprint, minSET:
        if not(type(xprint,boolean)) then
           xprint:=false:
        fi:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("UpLB(EP,r,N1,p)                                              \n");
     printf("      EP = etaproduct                                        \n");
     printf("       r = rational                                          \n");
     printf("      N1 = positive integer divisible by prime p             \n"); 
     printf("       p = prime                                             \n");
     printf(" EP assumed to be a modular function on level divisible by p^2\n");
     printf(" Let N = N1/p (assumed p^2|N1)                   \n");
     printf(" So if EP is MF on GAMMA0(N1) then Up(EP) is MF on GAMMA0(N) \n");
     printf(" This function returns a lower bound for ORD[G](Up(EP),r)       \n");
     printf(" where r is a cusp of GAMMA0(N)=G.                           \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 4 then
        ##jacid:=args[1]: CUSPS:=args[2]: WIDS:=args[3]: N:=args[4]:
        EP:=args[1]: r:=args[2]: N1:=args[3]: p:=args[4]:
    if modp(N1,p)<>0 then
      ERROR("p does not divide N1");
    fi:
    # Suppose EP is an etaprod on Gamma0(N1)
    # p|N1 and r=c/d is a cusp Gamma0(N) where N1=p*N if p^2|N1 and N1=N otherwise
    #
    if modp(N1,p^2)=0 then 
       N:=N1/p:
    else
       ERROR("N1 not divisible by p^2");
    fi:
    d:=denom(r):
    if modp(N,d)<>0 then
      ERROR("d does not divide N");
    fi:
    if vp(d,p) >= vp(N,p)/2 then
       LB:=  cuspORD(EP,N1,r/p)/p:  
if xprint then
       print("case 1",LB);
fi:
    fi:
    if vp(d,p) > 0 and vp(d,p)< vp(N,p)/2 then
       LB:=  cuspORD(EP,N1,r/p):
if xprint then
       print("case 2",LB);
fi:
    fi:
    if vp(d,p)=0 then
       minSET:=[seq(cuspORD(EP,N1,(r+b)/p),b=0..p-1)]:
       LB:=min(op(minSET)):
if xprint then
       print("case 3",LB);
fi:
    fi:
    RETURN(LB);
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 4`);
  fi:
end:
ETA[vetainf]:=proc()
  local epterm,cc,eprod,GP,ngp,w2:
  if nargs=0 then
     printf("-------------------------------------------------------------\n");
     printf("vetainf(epterm)                                          \n");
     printf("   epterm = c*eta(t[1]*tau)^r[1]* ... *eta(t[s]*tau)^r[s]    \n");
     printf("            ie. an eta-quotient times a constant. \n"); 
     printf("   Returns the v[oo](epterm).     \n");
     printf("                                                             \n");
     printf("-------------------------------------------------------------\n");
     elif nargs = 1 then
        epterm:=args[1]: 
        if type(epterm,constant) then
           RETURN(0):
        else
           eprod:=epterm/ETA[etaCONSTANT](epterm):
           GP:=GPmake(eprod):
           ngp:=nops(GP):
	   w2:=sum(GP[2*i-1]*GP[2*i],i=1..(ngp/2)):
           RETURN(w2/24):
        fi:
  else
     printf("nargs=%a\n",nargs);
     ERROR(`nargs must be 0 or 1`);
  fi:
end:
ETA[vp]:=proc(n,p)
   local x,n1:
   x:=0:
   if n=0 then
       RETURN(oo):
   else
   n1:=n:
   while type(n1,integer) do
      n1:=n1/p:
      x:=x+1:
   od:
   RETURN(x-1);
   fi:
end:

printf("TABLE TYPE ETA = %a\n",type(ETA,table));
savelib( ETA 
,
"c:\\cygwin64/home/fgarv/maple/mylib/ETA.mla");
printf("END ETA package\n");
## mylib  above must be changed to the name of the directory
## in which you want stuff stored