#!/usr/local/bin/maple # -*- maplev -*- # Nathaniel Shar # HW 25 # Experimental Mathematics # It is okay to link to this assignment on the course webpage. ############# # Problem 1 # ############# # The suggested formula is # exp(q) = mul((1-q^i)^(-mu(i)/i), i=2..infinity) # where mu is the Mobius function. # Proof: Take logs of both sides. The left side is q. # The right side is # sum(-mu(i)/i*(q^i + q^(2i)/2 + q^(3i)/3 + ...). # The coefficient of q^i is (1/i)sum(mu(d)), where the sum is # over d dividing i. The sum of the Mobius function over all # divisors of i is 1 if and only if i = 1, so the right hand side is # also q. Therefore, the formula is correct. ############# # Problem 2 # ############# # The ogf of the vector partitions is P(x)/P(x^2) * (P(x))^2 = # P(x)^3/P(x^2), where P(x) is the ogf of the partitions. This is # listed in OEIS as an alternative formula for the ogf of A182818. # OK, OK, I'll prove it myself. Take the log. We get # 3log(P(x)) - log(P(x^2)) # # = 3*[(q + q^2/2 + q^3/3 + ...) + (q^2 + q^4/2 + q^6/3 + ...) + ...] # - [(q^2 + q^4/2 + q^6/3 + ...) + (q^4 + q^8/2 + q^12/3 + ...) + ...] # # The coefficient of x^n is sum(3d/n, d|n, d odd) + sum(2d/n, d|n, d # even). # If n is odd, then we get 3/n*sum(d, d|n) = 3sigma(n)/n = # sigma(2)sigma(n)/n. Since sigma is multiplicative, this is # sigma(2n)/n. # If n is even, the sum of the even divisors of 2n is 2*sigma(n). # The sum of the odd divisors of 2n is the same as the sum of the odd # divisors of n. Adding these, the sum of the divisors of 2n is # 2*sigma(n) + sum(d, d|n, d odd) # = 3*sum(d, d|n, d odd) + 2*sum(d, d|n, d even), # which is n times the coefficient of x^n. # Hence the coefficient of x^n is equal to sigma(2n)/n, which proves # the result. # From PEP.txt: ######################################################################## freqtab:=proc(L) local fset,i,M,j: fset:={}: for i from 1 to nops(L) do fset:=fset union {L[i]}: if type(M[L[i]],integer) then M[L[i]]:=M[L[i]]+1: else M[L[i]]:=1: fi: od: for j in fset do lprint(j,M[j]); od; end: ######################################################################## etaq:=proc(q,i,T) # This proc returns (q^i:q^i)_inf up to q^T local k,x,z1,z,w: z1:=(i + sqrt( i*i + 24*T*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<=T then x:=x+ q^( w )*(-1)^k: fi: od: RETURN(x): end: ######################################################################## prodmake:=proc() local ft,f0,_a,_b,i,n,j,d,_A,_B,sum1,sum2,divj,divjb,m,prd,f,q,T,bseq,lst; #### #USAGE: prodmake(f,q,T) returns q-product # prodmake(f,q,T,list) returns q-product as a list of exponents # #### #04/01/00: added nargs bit #03/30/07: add expand for symbolic prods if nargs>4 then ERROR(` number of arguments must be 3 or 4.`); fi: f:=args[1]: q:=args[2]: T:=args[3]: #### ft:=series(f,q,T+5): if whattype(ft)=series then f0:=coeff(ft,q,0): if f0=1 then _b:=series(f,q,T): _b:=convert(_b,polynom): for i from 1 to T-1 do _B[i]:=coeff(_b,q,i): od: _A[1]:=_B[1]: for n from 2 to T-1 do sum2:=0: for j from 1 to n-1 do divj:=numtheory[divisors](j): sum1:=0: for d in divj do sum1:=expand(sum1+d*(_A[d])): od: sum2:=expand(sum2 + (_B[n-j])*sum1): od: sum1:=0: divjb:=numtheory[divisors](j) minus {n}: for d in divjb do sum1:=expand(sum1+d*(_A[d])): od: sum2:=expand(n*_B[n] - sum2 - sum1): _A[n]:=sum2/n: od: #### #04/01/00: added option of returning list if nargs=3 then prd:=product((1-q^m)^(-_A[m]),m=1..(T-1)): RETURN(prd): else bseq:=[seq(-_A[m],m=1..(T-1))]: lst:=convert(bseq,list): RETURN(lst): fi: else ERROR(`coeff of q^0 must be 1`); fi: else ERROR(`f must be a series`); fi: end: ######################################################################## etamake:=proc(f,q,last) # #VERSION: 01/30/10 # Use series(`leadterm`(expr),q) to get leadterm instead of tc etc. # This is handle BBG c(q) function #This procedure computes an eta-product expansion #of the series f to order O(q^last). local fp,tc,exq,g,aa,ld,h,hh,i,cf,etaprod,alast,sf: global ebase: sf:=series(f,q,last+10): fp:=convert(sf,polynom): tc:=tcoeff(fp,q): exq:=ldegree(fp,q): g:=normal(fp/tc/q^exq): aa:=tc: ld:=1: ebase:=1: alast:=last-exq: while ld>0 do h:=series(g-1,q=0,alast+1): hh:=0: for i from 1 to alast do hh:=hh+coeff(h,q,i)*q^i: od: h:=hh: ## ##08/20/99: bug fix ldegree(0,q) returns infinity ## in maple V release 5 if h=0 then ld:=0: else ld:=ldegree(h,q): fi: cf:=coeff(h,q,ld): if ld>0 then exq:=exq+(ld*cf)/24: aa:=eta(ld*tau)^(-cf)*aa: g:=g*etaq(q,ld,alast)^cf: ebase:=ilcm(ebase,ld): fi: od: etaprod:=q^(exq)*aa: RETURN(etaprod): end: ######################################################################## findcong:=proc() local s,j,c,S,p,r,QS2,S1,QS,T,LM,LM1,M,xxT,xxT1,CFS1,IGCD1,CFS,IGCD,f,IFX,HH,h,P,R,m,L,ra,Ma,pa,XSET; global xprint: ## findcong(QS,T) ## findcong(QS,T,LM) ## findcong(QS,T,LM,XSET) ## find linear congruences for q-series QS up to q^T ## EXAMPLE: ## P:=series(1/etaq(q,1,1001),q,1001): ## > findcong(P,1000); ## [4, 5, 5] ## [5, 7, 7] ## [6, 11, 11] ## [24, 25, 25] ## {[5, 7, 7], [4, 5, 5], [24, 25, 25], [6, 11, 11]} ## NOTE: XSET set of moduli to exclude if nargs>4 or nargs<2 then ERROR(`findcong takes 2, 3 or 4 args.`); fi: if not(type(xprint,boolean)) then xprint:=false: fi: QS:=args[1]: T:=args[2]: LM1:=2: XSET:={}: if nargs=2 then LM:=trunc(sqrt(T)): else LM:=args[3]: fi: if nargs=4 then XSET:=args[4]: fi: S:={}: for M from LM1 to LM do for r from 0 to M-1 do xxT:=trunc((T-r)/M): xxT1:=min(xxT,100): CFS1:=[seq(coeff(QS,q,M*n+r),n=0..xxT1)]: if xprint then lprint("xxT1=",xxT1,"CFS1=",CFS1); fi: IGCD1:=igcd(op(CFS1)); if xprint then lprint("M=",M,"r=",r,"IGCD1=",IGCD1); fi: if IGCD1<>1 and not member(IGCD1,XSET) then CFS:=[seq(coeff(QS,q,M*n+r),n=0..trunc((T-r)/M))]: IGCD:=igcd(op(CFS)); if IGCD<>1 then # new one? f:=1: IFX:=ifactors(IGCD): HH:=IFX[2]: for h from 1 to nops(HH) do P:=HH[h][1]: R:=HH[h][2]: for m from 1 to nops(S) do L:=S[m]: ra:=L[1]: Ma:=L[2]: pa:=L[3]: if modp(r,Ma)=ra and pa=P^R and modp(M,Ma)=0 then f:=0: fi: od: if f=1 then S1:= S union {[r,M,P^R]}: if nops(S1) > nops(S) then print([r,M,P^R]); S:=S1: fi: fi: od: fi: fi: od:od: RETURN(S); end: ######################################################################## 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: ##fixed 11-16-06 y:=y+coeff(st,q,(n*i+k))*q^i: od: RETURN(y): else ERROR(`s must be a series`); fi: end: ######################################################################## with(combinat): ptnDP:=proc(ptn) local np,k,d; #Returns true if ptn is a partition into distinct parts np:=nops(ptn); # np is number of parts of ptn if np=1 then RETURN(true): else for k from 1 to np-1 do d:=ptn[k+1]-ptn[k]: if d=0 then RETURN(false) fi: od: fi: RETURN(true): end: vpw:=proc(vptn) local P1,np; # weight of vector ptn [P1,P2,P3] = (-1)^(#(P1)) P1:=vptn[1]: np:=nops(P1): RETURN( (-1)^(np) ): end: vecptns:=proc(n) local j,PTNS,DPTNS,V,i,k,T,vp,vp1,vp2,vp3,s1,s2,s3; #GENERATES A LIST OF VECTOR PARTITIONS OF n (cor. to spt) for j from 0 to n do PTNS[j]:=partition(j): DPTNS[j]:=select(ptnDP,PTNS[j]): od: V:=[]: for i from 0 to n do for j from 0 to n-i do k:=n-i-j: T:=cartprod([DPTNS[i],PTNS[j],PTNS[k]]): while not T[finished] do vp:=T[nextvalue](): V:=[op(V),vp]: od: od: od: RETURN(V): end: lp:=proc(ptn) if nops(ptn)=0 then 0: else ptn[nops(ptn)]: fi: end: vpcrank:=proc(vptn) local P2,np2,P3,np3; # crank of vector ptn [P1,P2,P3] is #(P2)-#(P3) P2:=vptn[2]: np2:=nops(P2): P3:=vptn[3]: np3:=nops(P3): RETURN(np2 - np3): end: vecptnsC:=proc(n,k,t) #GENERATES new vector partitions of n with crank congruent to #k mod t local localCS,vptns: localCS:=proc(vptn) if modp(vpcrank(vptn),t)=k then true: else: false: fi: end: vptns:=vecptns(n): select(localCS,vptns); end: ######################################################################## nu2:=proc(F) local n,F1: n:=0: F1:=0: while F1=0 do if modp(F,2^(n+1))=0 then n:=n+1: else F1:=1: fi: od: RETURN(n): end: ########################################################################