qseries:=table(); qseries[aqprod]:=proc(a,q,n) local x,i: if type(n,nonnegint) then x:=1: for i from 1 to n do x := x * (1-a*q^(i-1)): od: else x:=``(a,q)[n]; fi: RETURN(x): end: qseries[etamake]:=proc(f,q,last) #This procedure computes an eta-product expansion #of the series f to order O(q^last). local fp,tc,exq,g,aa,ld,h,hh,i,cf,etaprod,alast: series(f,q,last+10): fp:=convert(",polynom): tc:=tcoeff(fp,q): exq:=ldegree(fp,q): g:=normal(fp/tc/q^exq): aa:=tc: ld:=1: alast:=last-exq: while ld>0 do h:=series(g-1,q=0,alast+1): hh:=0: for i from 1 to alast do hh:=hh+coeff(h,q,i)*q^i: od: h:=hh: ld:=ldegree(h,q): cf:=coeff(h,q,ld): if ld>0 then exq:=exq+(ld*cf)/24: aa:=eta(ld*tau)^(-cf)*aa: g:=g*etaq(q,ld,alast)^cf: fi: od: etaprod:=q^(exq)*aa: RETURN(etaprod): end: qseries[etaq]:=proc(q,i,trunk) # This proc returns (q^i:q^i)_inf up to q^trunk local k,x,z1,z,w: z1:=(i + sqrt( i*i + 24*trunk*i ) )/(6*i): z:=1+trunc( evalf(z1) ): x:=0: for k from -z to z do w:=i*k*(3*k-1)/2: if w<=trunk then x:=x+ q^( w )*(-1)^k: fi: od: RETURN(x): end: qseries[findhom]:=proc(L,q,n,topshift) #This program generates a list of potential homogeneous relations #of degree n among the list of q-series L. # local xx,x,y,z,nx,n2,n3,gg,goo,gooo,c,dh,jj,g,k,w,d,i,ct,D,ii,h,G,ANS; global X: if whattype(n)=integer and whattype(topshift)=integer then if n>0 and topshift>-20 then D:=nops(L); gg[1]:={}; xx:=expand(sum(x[i],i=1..D)^(n) ): nx:=nops(xx): n2:=nx+20+topshift:#step this up if some solutions are spurious n3:=n2+1: print(`# of terms `,n3); c:=array(1..nx,1..n2); dh:=array(1..nx,1..D); for jj from 1 to nx do G:=op(jj,xx): for i from 1 to D do h[i]:=degree(G,x[i]): dh[jj,i]:=h[i]: od: w:=product(L[ii]^h[ii],ii=1..D): d:=taylor(w,q=0,n3); for i from 1 to n2 do c[jj,i]:=coeff(d,q,i); od; od: ct:= linalg[transpose](c): gooo:=(linalg[kernel](ct)); for k from 1 to nops(gooo) do goo:=gooo[k]; d:=op(dh): g[k]:=0: for i from 1 to nx do g[k]:=g[k]+goo[i]*product(X[ii]^d[i,ii],ii=1..D): od; gg[k]:=g[k]; od; print(`-----RELATIONS----of order---`,n); ANS:=convert(gg,set); RETURN(ANS); else ERROR(`n must be a positive integer and topshift must apositive integer greater then -20`); fi: else ERROR(`n and topshift must be integers.`); fi: end: qseries[findhomcombo]:=proc(f,L,q,n,topshift,etaoption) #This program tries to express f as a homogeneous polynomial #of degree n among the list of q-series in L. # local xx,x,y,z,nx,n2,n3,gg,goo,gooo,c,dh,jj,g,k,w,d,i,ct,D,ii,h,G,ANS; global X: X:='X': if whattype(n)=integer and whattype(topshift)=integer then if n>0 and topshift>-20 then D:=nops(L); gg[1]:={}; xx:=expand(sum(x[i],i=1..D)^(n) ): nx:=nops(xx): n2:=nx+20+topshift:#step this up if some solutions are spurious n3:=n2+1: print(`# of terms `,n3); c:=array(1..nx+1,1..n2); dh:=array(1..nx,1..D); for jj from 1 to nx do G:=op(jj,xx): for i from 1 to D do h[i]:=degree(G,x[i]): dh[jj,i]:=h[i]: od: w:=product(L[ii]^h[ii],ii=1..D): d:=taylor(w,q=0,n3); func[jj]:=d: for i from 1 to n2 do c[jj,i]:=coeff(d,q,i); od; od: fq:=series(f,q,n3): for i from 1 to n2 do c[nx+1,i]:=coeff(fq,q,i):od: ct:= linalg[transpose](c): gooo:=(linalg[kernel](ct)): ANS:={}: ANSETA:={}: for k from 1 to nops(gooo) do goo:=gooo[k]: d:=op(dh): g[k]:=0: geta[k]:=0: if goo[nx+1]<>0 then for i from 1 to nx do g[k]:=g[k]+goo[i]*product(X[ii]^d[i,ii],ii=1..D): geta[k]:=geta[k]+goo[i]*qseries[etamake](func[i],q,n3-2): od: ANS:=ANS union {-g[k]/goo[nx+1]}: ANSETA:=ANSETA union {-geta[k]/goo[nx+1]}: fi: od: print(`-----possible linear combinations of degree------`,n): if etaoption = `yes` then print(ANSETA); fi: RETURN(ANS); else ERROR(`n must be a positive integer and topshift must apositive integer greater then -20`); fi: else ERROR(`n and topshift must be integers.`); fi: end: qseries[findnonhom]:=proc(L,q,n,topshift) #This program generates a list of potential nonhomogeneous relations #of degree n among the list of q-series L. # local xx,x,y,z,nx,n2,n3,gg,goo,gooo,c,dh,jj,g,k,w,d,i,ct,D,ii,h,G,ANS; global X: if whattype(n)=integer and whattype(topshift)=integer then if n>0 and topshift>-20 then X:='X'; D:=nops(L); gg[1]:={}; xx:=expand((sum(x[i],i=1..D)+1)^(n) ): nx:=nops(xx): n2:=nx+20+topshift:#step this up if some solutions are spurious n3:=n2+1: print(`# of terms `,n3); c:=array(1..nx,1..n2+1); dh:=array(1..nx,1..D); for jj from 1 to nx do G:=op(jj,xx): for i from 1 to D do h[i]:=degree(G,x[i]): dh[jj,i]:=h[i]: od: w:=product(L[ii]^h[ii],ii=1..D): d:=series(w,q=0,n3+4); for i from 0 to n2 do c[jj,i+1]:=coeff(d,q,i); od; od: ct:= linalg[transpose](c): gooo:=(linalg[kernel](ct)); for k from 1 to nops(gooo) do goo:=gooo[k]; d:=op(dh): g[k]:=0: for i from 1 to nx do g[k]:=g[k]+goo[i]*product(X[ii]^d[i,ii],ii=1..D): od; gg[k]:=g[k]; od; print(`-----RELATIONS----of order---`,n); ANS:=convert(gg,set); RETURN(ANS); else ERROR(`n must be a positive integer and topshift must apositive integer greater then -20`); fi: else ERROR(`n and topshift must be integers.`); fi: end: qseries[findnonhomcombo]:=proc() #This program tries to express f as a polynomial #of degree n among the list of q-series in L. # local xx,x,y,z,nx,n2,n3,gg,goo,gooo,c,dh,jj,g,k,w,d,i,ct,D,ii,h,G,ANSf,L,q,n,etaoption,topshift; global X: X:='X': f:=args[1]: L:=args[2]: q:=args[3]: n:=args[4]: if nargs>4 and whattype(args[5])=integer then topshift:=args[5]: if nargs=6 then etaoption:=args[6]: else etaoption:=no: fi: fi: if nargs=4 then topshift:=0: etaoption:=no: fi: if nargs=5 and whattype(args[5])<>integer then topshift:=0: etaoption:=args[5]: fi: ##print(` n= `,n, `topshift =`,topshift); if whattype(n)=integer and whattype(topshift)=integer then if n>0 and topshift>-20 then D:=nops(L); gg[1]:={}; xx:=expand((sum(x[i],i=1..D)+1)^(n) ): nx:=nops(xx): n2:=nx+20+topshift:#step this up if some solutions are spurious n3:=n2+1: print(`# of terms `,n3); c:=array(1..nx+1,1..n2+1); print(`matrix is `,nx+1,` x `,n2+1); dh:=array(1..nx,1..D); for jj from 1 to nx do G:=op(jj,xx): for i from 1 to D do h[i]:=degree(G,x[i]): dh[jj,i]:=h[i]: od: w:=product(L[ii]^h[ii],ii=1..D): d:=taylor(w,q=0,n3+4); func[jj]:=d: for i from 0 to n2 do c[jj,i+1]:=coeff(d,q,i); od; od: fq:=series(f,q,n3+4): for i from 0 to n2 do c[nx+1,i+1]:=coeff(fq,q,i): od: ct:= linalg[transpose](c): gooo:=(linalg[kernel](ct)): ANS:={}: ANSETA:={}: for k from 1 to nops(gooo) do goo:=gooo[k]: d:=op(dh): g[k]:=0: geta[k]:=0: if goo[nx+1]<>0 then for i from 1 to nx do g[k]:=g[k]+goo[i]*product(X[ii]^d[i,ii],ii=1..D): if etaoption=`yes` then geta[k]:=geta[k]+goo[i]*qseries[etamake](func[i],q,n3-2): fi: od: ANS:=ANS union {-g[k]/goo[nx+1]}: if etaoption=`yes` then ANSETA:=ANSETA union {-geta[k]/goo[nx+1]}: fi: fi: od: print(`-----possible linear combinations of degree------`,n): if etaoption = `yes` then print(ANSETA); fi: RETURN(ANS); else ERROR(`n must be a positive integer and topshift must apositive integer greater then -20`); fi: else ERROR(`n and topshift must be integers.`); fi: end: qseries[findpoly]:=proc() #This proc looks for a polynomial local x,y,q,deg1,deg2,ARGLIST,TYPELIST,CORRECTL,check,dim1,dim2,A,num,k,j,B,qq,l,kk,POLY,kkk,i,POLYg,POLYfunc,ss,numtyp: global X,Y: lprint(`WARNING: X,Y are global.`); #relation between x,y of degree deg1 in x, and degree deg2 in y. #The relation is checked to order O(q^check). if nargs<5 then ERROR(` number of arguments must be 5 or 6.`); fi: x:=args[1]: y:=args[2]: q:=args[3]: deg1:=args[4]: deg2:=args[5]: ARGLIST:=[q,deg1,deg2]; TYPELIST:=map(whattype,ARGLIST); CORRECTL:=[string,integer,integer]: if TYPELIST=CORRECTL then if nargs=6 then check:=args[6]: fi: if nargs>6 then ERROR(` findpoly can at most 6 arguments`); fi: dim1:=(deg1+1)*(deg2+1): dim2:=dim1+10: print(` dims `, dim1, dim2); A:=array(1..dim1,1..dim2): num:=0: for k from 0 to deg1 do for j from 0 to deg2 do num:=num+1: B[num]:=X^k*Y^j: qq:=series(x^k*y^j,q,dim2+3): for l from 0 to (dim2-1) do A[num,l+1]:=coeff(qq,q,l): od: od: od: with(linalg): transpose(A): kk:=kernel("): if kk={} then lprint(` NO polynomial relation found. `); else POLY:=0; op(kk); kkk:="; for i from 1 to num do POLY:=POLY+B[i]*kkk[i]: od: POLYg:=0; for j from 0 to deg2 do POLYg:=POLYg+factor(coeff(POLY,Y,deg2-j))*Y^(deg2-j): od: if nargs=6 then print(`The polynomial is`); print(POLYg); POLYfunc:=unapply(POLYg,X,Y): ss:=series(POLYfunc(x,y),q,check+1): print(`Checking to order`,check); print(ss); RETURN(POLYg): else print(`The polynomial is`); fi: RETURN(POLY): fi: else numtyp:=nops(TYPELIST): ERROR(` Wrong type of argument. ARGLIST has type `,seq(TYPELIST[i],i=1..numtyp), `It should have type`,seq(CORRECTL[i],i=1..numtyp)); fi: end: qseries[jacprod]:=proc(a,b,q,T) local x,i,lasti: tripleprod(q^a,q^b,T)/tripleprod(q^(b),q^(3*b),T): end: ####################################################################### # qseries[jacprodmake] ####################################################################### macro(jaccheck = `jacprodmake/jaccheck`, periodfind = `jacprodmake/periodfind`, jacmake = `jacprodmake/jacmake`); ####################################################################### # FUNCTION : periodfind - finds the smallest period of a finite list # CALLING SEQUENCE : periodfind(A,T) # PARAMETERS : A - list of integers # P - positive integer # SYNOPSIS : # periodfind(A,T) returns P if # A[i] = A[i+P] for each i. ####################################################################### periodfind:=proc(A,T) local num,n,perfi,pp,dud,i,j,lastj,diffc; num:=T: n:=trunc(num/2): perfi:=0: pp:=2: while perfi=0 do dud:=0: while dud=0 do for i from 1 to pp do lastj:=trunc(num/pp)-2: for j from 0 to lastj do diffc:=A[i+pp*j]-A[i+pp*j+pp]: if diffc<>0 then dud:=1: fi: od: od: if dud=0 then perfi:=1: RETURN(pp): dud:=1: fi: od: pp:=pp+1: od: end; ####################################################################### # FUNCTION : jaccheck - checks whether an array corresponds a # theta-product. # CALLING SEQUENCE : jaccheck(A,P): # PARAMETERS : A - list of integers # P - positive integer # SYNOPSIS : # jaccheck(A,P) returns 1 if # A[i] = A[P-i] for each i, 1 < i < (P/2); # otherwise 0 is returned. If jaccheck(A,P)=1 and A is periodic # with period P, then A corresponds to a certain theta product. ####################################################################### jaccheck:=proc(A,P) local n,jac,i: n:=trunc(P/2): jac:=1: for i from 1 to n do if A[i]<>A[P-i] then jac:=0: RETURN(0): fi: od: if jac=1 then RETURN(1): fi: RETURN(-1): end: ####################################################################### # FUNCTION : jacmake - connverts a list of integers into a # theta-product if possible. # CALLING SEQUENCE : jacmake(A,T): # PARAMETERS : A - list of integers # T - positive integer # SYNOPSIS : T # If f = product((1-q^i)^(-A[i]),i=1..T-1) + O(q ) # then # jacmake(A,T) will return a product of the form # # product( JAC(j,P,oo), j in S) # T # if this product agrees with f to O(q ). # # Here JAC(j,P,oo) corresponds to the theta-function # tripleprod(q^j,q^P,oo) when 0 < j < P # and # JAC(0,P,oo) corresponds to the eta-product etaq(q,P,oo). ####################################################################### jacmake:=proc(A,T) local pp,y,mp,i: pp:=periodfind(A,T): if pp>(T/2) then print(` no period found `); print(` This means either that there is no period,`); print(` or a larger value of T should be chosen.`); RETURN(): else if jaccheck(A,pp)<>1 then print(` not a theta function `); RETURN(): else y:=1: mp:=modp(pp,2): if mp=0 then for i from 1 to ((pp)/2-1) do y:=y/(JAC(i,pp,infinity)/JAC(0,pp,infinity))^A[i]: od: y:=y/(JAC(pp/2,pp,infinity)/JAC(0,pp,infinity))^(A[pp/2]/2): else for i from 1 to (pp-1)/2 do y:=y/(JAC(i,pp,infinity)/JAC(0,pp,infinity))^A[i]: od: fi: y:=y/JAC(0,pp,infinity)^A[pp]: RETURN(y): fi: fi: end: ####################################################################### # THE MAIN FUNCTION qseries[jacprodmake] ####################################################################### qseries[jacprodmake]:=proc(f,q,T) local ft,f0,_a,_b,i,n,j,d,_B,sum1,sum2,divj,divjb,m,_A: ft:=series(f,q,T+5): if whattype(ft)=series then f0:=coeff(ft,q,0): if f0=1 then _b:=series(f,q,T): _b:=convert(_b,polynom): for i from 1 to T-1 do _B[i]:=coeff(_b,q,i): od: _A[1]:=_B[1]: for n from 2 to T-1 do sum2:=0: for j from 1 to n-1 do divj:=numtheory[divisors](j): sum1:=0: for d in divj do sum1:=sum1+d*(_A[d]): od: sum2:=sum2 + (_B[n-j])*sum1: od: sum1:=0: divjb:=numtheory[divisors](j) minus {n}: for d in divjb do sum1:=sum1+d*(_A[d]): od: sum2:=n*_B[n] - sum2 - sum1: _A[n]:=sum2/n: od: RETURN(jacmake(_A,T-1)): else ERROR(`coeff of q^0 must be 1`); fi: else ERROR(`f must be a series`); fi: end: macro(jaccheck = jacprod, periodfind = periodfind, jacmake = jacmake); ####################################################################### # qseries[jac2prod] ####################################################################### macro(jac = `jac2prod/jac`); jac:=proc(a::integer, b::integer, c::infinity) local Q; if a0 then Q:=` `(q^a,q^b)[infinity]*` `(q^(b-a),q^b)[infinity]*` `(q^b,q^b)[infinity]; RETURN(Q): else RETURN(` `(q^b,q^b)[infinity]); fi: else ERROR(` a,b must satisy a0 then if a<>0 then t:=max(T,abs(b)): c:=trunc(1/2*(b-2*a+(4*a^2-4*a*b+b^2+8*b*t)^(1/2))/b)+1: RETURN(tripleprod2(q^a,q^b,c)); else RETURN(qseries[etaq](q,b,T)); fi: else ERROR(` T must be positive. `); RETURN(): fi: end: ####################################################################### # qseries[jac2series] ####################################################################### ####################################################################### # FUNCTION : jac2series - converts a product of theta functions into # q-product form. # CALLING SEQUENCE : jacseries(jacexpr) # PARAMETERS : jacexpr - product of JAC(i,j,oo) # where i,j are integers # and 0 < i < j. # = # SYNOPSIS : # jacprod(jacexpr) returns a q-product;ie a product of the # a b # of functions of the form (q , q ). ####################################################################### qseries[jac2series]:=proc(JEXPR, T::integer) local JJ; JJ:=JEXPR: RETURN(eval(subs({JAC=qjac,infinity=T},JJ))): end: macro(tripleprod2=tripleprod2, qjac = qjac): qseries[prodmake]:=proc(f,q,T) local ft,f0,_a,_b,i,n,j,d,_A,_B,sum1,sum2,divj,divjb,m: ft:=series(f,q,T+5): if whattype(ft)=series then f0:=coeff(ft,q,0): if f0=1 then _b:=series(f,q,T): _b:=convert(_b,polynom): for i from 1 to T-1 do _B[i]:=coeff(_b,q,i): od: _A[1]:=_B[1]: for n from 2 to T-1 do sum2:=0: for j from 1 to n-1 do divj:=numtheory[divisors](j): sum1:=0: for d in divj do sum1:=sum1+d*(_A[d]): od: sum2:=sum2 + (_B[n-j])*sum1: od: sum1:=0: divjb:=numtheory[divisors](j) minus {n}: for d in divjb do sum1:=sum1+d*(_A[d]): od: sum2:=n*_B[n] - sum2 - sum1: _A[n]:=sum2/n: od: product((1-q^m)^(-_A[m]),m=1..(T-1)): RETURN("): else ERROR(`coeff of q^0 must be 1`); fi: else ERROR(`f must be a series`); fi: end: qseries[qbin]:=proc(q,m,n) if whattype(m)=integer and whattype(n)=integer then if m>=0 and m<=n then RETURN(normal(aqprod(q,q,n)/aqprod(q,q,m)/aqprod(q,q,n-m))); else RETURN(0); fi: else ERROR(` m and n must be integers.`); fi: end: ####################################################################### # qseries[qfactor] ####################################################################### macro(qsplit = `qfactor/qsplit`, qpfac = `qfactor/qpfac`, qfac=`qfactor/qfac`,qfopexp=`qfactor/qfopexp`, bigqpfac=`qfactor/bigqpfac`); ####################################################################### # FUNCTION : qsplit - returns a list [f,tc,ld,td] # CALLING SEQUENCE : qsplit(f) # PARAMETERS : f - a polynomial in q # SYNOPSIS : # qsplit(f) returns the list [f,tc,ld,td] # where f = tc*(q^ld + ........+()q^td. ####################################################################### print(`qsplit:=proc(f::polynom)`); qsplit:=proc(f::polynom) local tc,ld,td; tc := tcoeff(f, q); ld := ldegree(f, q); td := degree(f, q); RETURN([f,tc,ld,td]); end; ####################################################################### # FUNCTION : qpfac - factor a q polynom into factors (1-q^i) # CALLING SEQUENCE : qpfac(f,T,s) # qpfac(f) # qpfac(f,T) # qpfac(f,s) # PARAMETERS : f - a polynomial in q # T - upper bound for i (optional) # s - string (optional) # SYNOPSIS : # qpfac(f) attempts to factor the poly q into factors (1-q^i) # if it fails it simply returns the input f. # If the third (or second) argument s=test is given then qpfac # will return FAIL if f can not be q-factored. ####################################################################### print(`qpfac:=proc()`); qpfac:=proc() local m1, m2, LBITS, tc, ld, g, gl, gl1, gl2, gl3, gd, ggc, gprod, gser, gno; f:=args[1]: if nargs=1 or (nargs=2 and type(args[2],integer)<>true) then m1:=4: m2:=3: else m1:=0: m2:=args[2] fi: LBITS:=qsplit(f); tc:=LBITS[2]: ld:=LBITS[3]: g:=convert(f/tc/q^ld,polynom); gl:=convert(g,list); gl1:=subs(q=1,gl); gl2:=map(whattype,gl1); gl3:=convert(gl2,set); if gl3 = '{integer}' then gd:=degree(g,q); ggc:=coeff(g,q,0); gprod:=prodmake(g,q,m1*gd+m2); gser:=convert(series(g-gprod,q,m1*gd+m2+5),polynom); if gser<>0 then if nargs=2 and args[2]=test then RETURN(FAIL); fi: if nargs=3 and args[3]=test then RETURN(FAIL); fi: if nargs=1 or (nargs=2 and args[2]<>test) then RETURN(f); fi: else gno:=normal(gprod-g); if gno = 0 then RETURN(tc*gprod*q^ld); else if nargs=2 and args[2]=test then RETURN(FAIL); fi: if nargs=3 and args[3]=test then RETURN(FAIL); fi: if nargs=1 or (nargs=2 and args[2]<>test) then RETURN(f); fi: fi: fi: else if nargs=2 and args[2]=test then RETURN(FAIL); fi: if nargs=3 and args[3]=test then RETURN(FAIL); fi: if nargs=1 or (nargs=2 and args[2]<>test) then RETURN(f); fi: fi: end; ####################################################################### # FUNCTION : qfopexp - factors a SINGLE q-expression and factors # as much as possible into terms of the form (1-q^i) # CALLING SEQUENCE : qfopexp(f,T,s) # qfopexp(f) # qfopexp(f,T) # qfopexp(f,s) # PARAMETERS : f - a single q-expression (ie. a factor of a factored # expression. # T - upper bound for i (optional) # s - string (optional) # SYNOPSIS : # An expression that cant be q-factored is left alone. ####################################################################### print(` qfopexp:=proc()`); qfopexp:=proc() expr:=args[1]; wh:=whattype(expr); n:=nops([op(expr)]); whL:=map(whattype,[op(expr)]); if n=2 and whL[1]=`+` and whL[2]=integer then cond:=1: else cond:=0; fi; if wh=`+` or cond=1 then s:=subs(q=2,expr); if type(s,rational)=true then RETURN(qfac(args)); else RETURN(expr); fi: else RETURN(expr); fi; end; ####################################################################### # FUNCTION : bigqpfac - factor a q polynom into factors (1-q^i) # CALLING SEQUENCE : bigqpfac(f,T,s) # bigqpfac(f) # bigqpfac(f,T) # bigqpfac(f,s) # PARAMETERS : f - a polynomial in q # T - upper bound for i (optional) # s - string (optional) # SYNOPSIS : # bigqpfac(f) attempts to factor the poly q into factors (1-q^i) # if it fails it simply returns the input f. # If the third (or second) argument s=test is given then qpfac # will return FAIL if f can not be q-factored. # This is a version of qpfac. This version does not fail if it comes # across a term with other variables beside q - these are simply # left alone. ####################################################################### print(` bigqpfac:=proc()`); bigqpfac:=proc() local k; expr:=args[1]: na:=nargs; if na=1 then qtmp:=x->qfopexp(x); else rst:=seq(args[k],k=2..nargs); qtmp:=x->qfopexp(x,rst); fi: L1:=[op(expr)]; L2:=map(qtmp,L1); if type(expr,`^`)=true then new:=L2[1]^L2[2]; if nops(L2)>2 then ERROR(`nops(L2)>2 in bigqfac`); fi: else ##new:=convert(L2,`*`); new:=convert(L2,whattype(expr)); fi: RETURN(new); end; ####################################################################### # 2ND MAIN FUNCTION qseries[qfac] ####################################################################### qfac:=proc() local f,T; f:=args[1]: if nargs=1 then if type(f,ratpoly)=true then if type(f,polynom) then RETURN(qpfac(f)); else RETURN(qpfac(numer(f))/qpfac(denom(f))); fi: else ERROR(`f must be ratpoly`); fi: else T:=args[2]: if type(f,ratpoly)=true then if type(f,polynom) then RETURN(qpfac(f,T)); else RETURN(qpfac(numer(f),T)/qpfac(denom(f),T)); fi: else ERROR(`f must be ratpoly`); fi: fi: end: ####################################################################### # THE MAIN FUNCTION qseries[qfactor] ####################################################################### print(` qseries[qfactor]:=proc()`); qseries[qfactor]:=proc() local f,T,k; f:=factor(args[1]): if type(f,ratpoly)=true then if type(f,polynom) then RETURN(bigqpfac(f)); else n:=nargs; if n=1 then RETURN(bigqpfac(numer(f))/bigqpfac(denom(f))); else newargs1:=numer(f),seq(args[k],k=2..n); newargs2:=denom(f),seq(args[k],k=2..n); RETURN(bigqpfac(newargs)/bigqpfac(newargs)); fi: fi: else ERROR(`f must be ratpoly`); fi: end: macro(qsplit = qsplit, qpfac = qpfac, qfac=qfac,qfopexp=qfopexp, bigqpfac=bigqpfac); qseries[quinprod]:=proc(z,q,T) local x,i,j,lasti: # If T = prodid # then a product form of the identity is returned. # If T = seriesid # then a series form of the identity is returned. # If T is a positive integer then a series to O(q^T) # is returned. # Determine output form if type(T,posint) then x:=0: lasti:=trunc( 7/6+1/6*(73+24*T)^(1/2))+1: for i from -lasti to lasti do x:=x + ( (-z)^(-3*i)-(-z)^(3*i+1) ) *q^(i*(3*i+1)/2): od: RETURN(x): else if T=prodid or T=seriesid then leftid:=` `(-z,q)[infinity] *` `(-q/z,q)[infinity] *` `(z^2*q,q^2)[infinity] *` `(1/z^2*q,q^2)[infinity] *` `(q,q)[infinity]: if T=prodid then x1:=` `(q^2/z^3,q^3)[infinity] *` `(q*z^3,q^3)[infinity] *` `(q^3,q^3)[infinity]: x2:=` `(q/z^3,q^3)[infinity] *` `(q^2*z^3,q^3)[infinity] *` `(q^3,q^3)[infinity]: x:= x1+z*x2: RETURN(leftid=x): else rightid:=sum(( (-z)^(-3*m)-(-z)^(3*m-1) ) *q^(m*(3*m+1)/2),m=-infinity..infinity): RETURN(leftid=rightid): fi: else fi: fi: end: qseries[sift]:=proc(s,q,n,k,T) # # sum a_i q^i --> sum a_[ni+k] q^i # local y,i,st,lasti: st:=series(s,q,T+5): if whattype(st)=series then y:=0: lasti:=floor((T-k)/n): for i from 0 to lasti do y:=y+coeff(s,q,(n*i+k))*q^i: od: RETURN(y): else ERROR(`s must be a series`); fi: end: qseries[theta]:=proc(z,q,T) local x,i: x:=0: for i from -T to T do x := x + z^i*q^(i*i): od: RETURN(x): end: qseries[theta2]:=proc(q,T) local N: N:=trunc(sqrt(T))+2: RETURN(qseries[theta](q,q,N)*q^(1/4)): end: qseries[theta3]:=proc(q,T) local N: N:=trunc(sqrt(T))+1: RETURN(qseries[theta](1,q,N)): end: qseries[theta4]:=proc(q,T) local N: N:=trunc(sqrt(T))+1: RETURN(qseries[theta](-1,q,N)): end: qseries[tripleprod]:=proc(z,q,T) # If T = seriesid # then a series form of the identity is returned. # If T is a positive integer then a series to O(q^T) # is returned. # Determine output form if type(T,posint) then x:=0: lasti:=trunc( sqrt(2*T+1/4)+1/2)+1: for i from -lasti to lasti do x:=x + (-1)^i*z^i*q^(i*(i-1)/2): od: RETURN(x): else if T=seriesid then leftid:=` `(z,q)[infinity] *` `(q/z,q)[infinity] *` `(q,q)[infinity]: rightid:=sum((-1)^m*z^m*q^(m*(m-1)/2), m=-infinity..infinity): RETURN(leftid=rightid): else ERROR(`Invalid T. T must either be a positive integer or the string seriesid`); fi: fi: end: qseries[winquist]:=proc(a,b,q,T) local x,i,j,lasti: x:=0: lasti:=trunc( 7/6+1/6*(25+24*T)^(1/2))+1: for i from 0 to lasti do for j from -lasti to lasti do x:=x + (-1)^(i+j)*( (a^(-3*i)-a^(3*i+3))*(b^(-3*j)-b^(3*j+1)) + (a^(-3*j+1)-a^(3*j+2))*(b^(3*i+2)-b^(-3*i-1))) *q^( 3*i*(i+1)/2 + j*(3*j+1)/2 ): od: od: RETURN(x): end: save( qseries, `jacprodmake/jaccheck` , `jacprodmake/periodfind` , `jacprodmake/jacmake` , `jac2prod/jac` , `jac2series/tripleprod2` , `jac2series/qjac` , `qfactor/qsplit` , `qfactor/qpfac` , `qfactor/qfac` , `qfactor/qfopexp` , `qfactor/bigqpfac` , ###NOTE: `/home/fac29/frank/maple/mylib/qseries.m`); ### <--- this line should be edited to reflect your setup quit; ###