Required packages:with(LinearAlgebra):We use a precision of 50 digits for floating point calculations.Digits:=50;Rogers dilogarithm function:L:=proc(x):
if x=0 then: return 0; fi:
if x=1 then: return Pi^2/6; fi:
return dilog(1-x) + log(x)*log(1-x)/2;
end:ser(f,t,var)
Given a function `f` in variable `var`, find its Taylor series at var=0 up to `var^t` terms.
Try:
ser(1/(1-q),10,q);
ser:=proc(f,t,var):
return(convert(series(f,var=0,t+1),polynom)):
end:serexp(f,t:=6,var:=epsilon)
Given a function `f` in variable `var`, find the Taylor series of
1 + f + (f^2/2!) + ... + (f^t/t!) at var=0 up to `var^t` terms.
Default value of `t` is 6 and that of `var` is epsilon.
Try:
serexp(q+q^2,3,q);serexp:=proc(f,t:=6,var:=epsilon)
local ans:=1,i,temp:=1:
for i from 1 to t do:
temp:=ser(temp*f,t,var):
ans:=ans+ temp/factorial(i):
end:
end:findC(L)
Given a list L = (L[1],L[2],...L[n]), consider the periodic product
P = (1-q)^L[1] * (1-q^2)^L[2] * ... * (1-q^n)^L[n] * (1-q^(n+1))^L[1] ...
Finds a constant C such that q^C * P is modular.
Try (to find C corresponding to first Rogers-Ramanujan identity):
findC([-1,0,0,-1,0]);
findC:=proc(L)
local i, m, ans:=0:
m:=nops(L):
for i from 1 to m-1 do:
ans:=ans+L[i]*m/24 - i*(m-i)*L[i]/4/m;
end:
end:ppoly(p,var:=X)
Define a sequence of polynomials in variable `var` (default=X)
for p=1,2,....
such that:
ppoly(1,X)= X,
ppoly(p+1, X) = (X^2 + X) * d(ppoly(p,X))/dX
We have that:
Li_{2-r}(Q) = ppoly(r-1,Q/(1-Q)) for r \134geq 1.
Try:
ppoly(1,X);
subs(X=Q/(1-Q), ppoly(1,X));
ppoly(2,X);
subs(X=Q/(1-Q), ppoly(2,X));
ppoly:=proc(p,var:=X):
if p=1 then: return var: fi:
return diff(ppoly(p-1,var),var)*(var^2+var):
end:exparg(B,xi,T,J:=1,cutoff:=4,var:=epsilon)
This is the argument of exponential in equation (12) of the paper.
B = B_i,
xi = xi_i,
T = t_i,
J = J_i,
cutoff = we take the summation over `p` up to p=cutoff. Default is 4.
var = this is actually the sqrt(epsilon) from the paper. so, the epsilon here is the sqrt(epsilon) from the paper.
exparg:=proc(B,xi,T,J:=1,cutoff:=4,var:=epsilon)
local p:
return (B + xi/2)*T*var
- add(
J^(p-1)
*expand(bernoulli(p,T/var)*var^(2*p-2))
*subs(X=xi/J,ppoly(p-1,X))/factorial(p)
, p=3..cutoff):
end:int2all(f,Atinv,cutoff:=30,var1:=t1,var2:=t2)consider a function `f` in varialbes var1 and var2. (Default t1, t2, respectively.)
consider a matrix At (which is A-tilde from equations (9) and (14)) whose inverse is Atinverse
calculates the value from equation (14) with C_2p replaced by f.
we cutoff f by expanding it up to terms var1^cutoff * var2 ^ cutoff.int2all:=proc(f,Atinv,cutoff:=30,var1:=t1,var2:=t2)
option remember:
local xv,ip,ex,g,i,j,ans:=0,x1,x2:
xv:=Matrix(2,1,[x1,x2]);
ip:= ((Transpose(xv).Atinv.xv)[1][1])/2;
ex:=expand(ser(ser(exp(ip),cutoff+1,x1),cutoff+1,x2));
g:=ser(ser(f, cutoff+1,var1),cutoff+1,var2):
for i from 0 to cutoff do:
for j from 0 to cutoff do:
ans:=ans+coeff(coeff(g,var1,i),var2,j)
*factorial(i)*factorial(j)*coeff(coeff(ex,x1,i),x2,j):
od:
od:
return ans:
end:
We now look at the Capparelli identities
CapA:=Matrix(2,2,[4,6,6,12]);
Capxi1:=3;
Capxi2:=24;
CapAt:=CapA+DiagonalMatrix([Capxi1,Capxi2]);
CapAtinv:=MatrixInverse(CapAt);expargcutoff is the cutoff on the values of p in equation (12)
expcutoff is cutoff for the exponential function in equation (12)
expargcutoff:=10;
expcutoff:=expargcutoff-2;
Now we multiply the exponentials in equation (12) for i=1 and i=2 to to get the C_p polynomials in equation (13)Note that the epsilon here is sqrt(epsilon) from the paper.
Capwanted calculates the integrals (14) for all polynomials C_p. The terms corresponding to odd values of p are 0
and the ones corresponding to even values are the ones that lead to constants c_p in (14).
Capexptwovar:=serexp(ser(exparg(B1,Capxi1,T1,1,expargcutoff)+exparg(B2,Capxi2,T2,3,expargcutoff),expcutoff,epsilon), expcutoff):
Capwanted:=ser(int2all(Capexptwovar,CapAtinv,30,T1,T2),expargcutoff,epsilon):
Capc1 is the constant c1 from equation (14).
Similarly, Capc2 is c2, Capc3 is c3, Capc4 is c4.
Note again that epsilon here is the sqrt(epsilon) from the paper.
Capc1:=coeff(Capwanted,epsilon,2):
Capc2:=coeff(Capwanted,epsilon,4):
Capc3:=coeff(Capwanted,epsilon,6):
Capc4:=coeff(Capwanted,epsilon,8):
The gamma value is stored in Capg
Capg:=simplify(C + (2*Capxi1+1)/24 + (2*Capxi2+3)/24);
The expressions on the right-hand sides of (17) - (19) are stored below.
Capexpr1:=simplify(Capc1 - Capg):
Capexpr2:=simplify(Capc2 - Capg*Capc1 + Capg^2/2):
Capexpr3:=simplify(Capc3 - Capg*Capc2 + Capg^2/2*Capc1 - Capg^3/6):
Capexpr4:=simplify(Capc4 - Capg*Capc3 + Capg^2/2*Capc2 - Capg^3/6*Capc1 + Capg^4/24):
Let us check that the constraints are all equal to 0 for the first Capparelli identity
CapC:=findC([0,-1,-1,0,0,0,0,0,-1,-1,0,0]);
subs(C=CapC,B1=0,B2=0,Capexpr1);
subs(C=CapC,B1=0,B2=0,Capexpr2);
subs(C=CapC,B1=0,B2=0,Capexpr3);
subs(C=CapC,B1=0,B2=0,Capexpr4);
Now we consider general B1 and B2. We solve equation (16) for C and substitute this into (17) - (19).
CapCval:=simplify(solve(Capexpr1=0,C));
Capeq1:=expand(simplify(subs(C=CapCval,Capexpr2)));
Capeq2:=expand(simplify(subs(C=CapCval,Capexpr3)));
Capeq3:=expand(simplify(subs(C=CapCval,Capexpr4)));
We now solve these equations simultaneously:
solve({Capeq1=0,Capeq2=0,Capeq3=0},{B1,B2});Now for the second Capparelli identityHere, we have to take a sum of two expressions, as in equation (42).
The `beta` constants have some factors common.
Below, mp1 and mp2 denote those parts of beta1 and beta2 that are different
CapQ1, CapQ2 are the Q1 and Q2 constants corresponding to the matrix CapA.
CapQ1:=3/4;
CapQ2:=2/(3^(2/3));
mp1:=(CapQ1^B1*CapQ2^B2);
mp2:=(CapQ1^B1p*CapQ2^B2p);Now we calculate analogues of expressions appearing in (16) - (19) for this case.
We do not print out these expressions; they are quite big.
Cap2expr1:=simplify(mp1*Capexpr1 + mp2*subs(B1=B1p,B2=B2p, C=C+Cp, Capexpr1)):
Cap2expr2:=simplify(mp1*Capexpr2 + mp2*subs(B1=B1p,B2=B2p, C=C+Cp, Capexpr2)):
Cap2expr3:=simplify(mp1*Capexpr3 + mp2*subs(B1=B1p,B2=B2p, C=C+Cp, Capexpr3)):
Cap2expr4:=simplify(mp1*Capexpr4 + mp2*subs(B1=B1p,B2=B2p, C=C+Cp, Capexpr4)):Now we check if the sum-side for the second Capparelli identity gives all these expressions to be 0:Cap2C:=findC([-1, 1, -1, 0, -1, 0, -1, 0, -1, 1, -1, 0]);
subs(C=Cap2C,B1=1,B2=3,B1p=3,B2p=6,Cp=1,Cap2expr1);
subs(C=Cap2C,B1=1,B2=3,B1p=3,B2p=6,Cp=1,Cap2expr2);
subs(C=Cap2C,B1=1,B2=3,B1p=3,B2p=6,Cp=1,Cap2expr3);
subs(C=Cap2C,B1=1,B2=3,B1p=3,B2p=6,Cp=1,Cap2expr4);Now we check numerically if there are other values of B_1, B_2, B_1', B_2', C' that satisfy the constraints.lb:=0;
ub:=6;
for b1 from lb to ub do:
for b2 from lb to ub do:
for b1p from b1 to ub do:
for b2p from b2 to ub do:
for cp from lb to ub do:
cval:=solve(subs(B1=b1,B2=b2,B1p=b1p,B2p=b2p,Cp=cp,Cap2expr1)=0);
tempc2:=abs(subs(C=cval,B1=b1,B2=b2,B1p=b1p,B2p=b2p,Cp=cp, Cap2expr2)):
tempc3:=abs(subs(C=cval,B1=b1,B2=b2,B1p=b1p,B2p=b2p,Cp=cp, Cap2expr3)):
tempc4:=abs(subs(C=cval,B1=b1,B2=b2,B1p=b1p,B2p=b2p,Cp=cp, Cap2expr4)):
if (evalf(tempc2) < 0.001 and evalf(tempc3) < 0.001) and evalf(tempc4) < 0.001 then:
print("C=",cval," B1=",b1," B2=",b2," B1'=",b1p," B2'=",b2p," C'=",cp):
print("---"):
fi:
od:
od:
od:
od:
od:Now we look at the Mod-9 identities. Since the analysis follows similarly, we skip the explanations.Mod9A:=Matrix(2,2,[2,3,3,6]);
Mod9At:=Mod9A+DiagonalMatrix([xi1,xi2]);
Mod9Atinv:=MatrixInverse(Mod9At);Mod9Q1:=1-2*sin(Pi/18);
Mod9Q2:=(4*sin((1/18)*Pi)^2+4*sin((1/18)*Pi))^(1/3);
simplify(Mod9Q1^3 - 3*Mod9Q1^2 + 1);
simplify(Mod9Q2^9 - 6*Mod9Q2^6 + 3*Mod9Q2^3 + 1);
We shall substitute actual values for xi's later:
alias(zeta1=RootOf(x^3-3*x-1,x));
alias(zeta2=RootOf(x^3-9*x^2-54*x-27,x));
expargcutoff:=10;
expcutoff:=expargcutoff-2;exptwovar:=serexp(ser(exparg(B1,xi1,T1,1,expargcutoff)+exparg(B2,xi2,T2,3,expargcutoff),expcutoff,epsilon),expcutoff):
Mod9wanted:=ser(int2all(exptwovar,Mod9Atinv,20,T1,T2),expargcutoff,epsilon):Mod9c1:=coeff(Mod9wanted,epsilon,2):
Mod9c2:=coeff(Mod9wanted,epsilon,4):
Mod9c3:=coeff(Mod9wanted,epsilon,6):
Mod9c4:=coeff(Mod9wanted,epsilon,8):Mod9g:=simplify(C+ (2*zeta1+1)/24 + (2*zeta2+3)/24);
Mod9expr1:=simplify(subs(xi1=zeta1,xi2=zeta2, Mod9c1 - Mod9g)):
Mod9expr2:=simplify(subs(xi1=zeta1,xi2=zeta2, Mod9c2 - Mod9g*Mod9c1 + Mod9g^2/2)):
Mod9expr3:=simplify(subs(xi1=zeta1,xi2=zeta2, Mod9c3 - Mod9g*Mod9c2 + Mod9g^2/2*Mod9c1 - Mod9g^3/6)):Now we check that the three Mod-9 identities satisfy the required constraints.
Identity 1:id1C:=findC([-1,0,-1,0,0,-1,0,-1,0]);
evalf(subs(C=id1C,B1=0,B2=0,Mod9expr1));
evalf(subs(C=id1C,B1=0,B2=0,Mod9expr2));
evalf(subs(C=id1C,B1=0,B2=0,Mod9expr3));Identity 2:id2C:=findC([0,-1,-1,0,0,-1,-1,0,0]);
evalf(subs(C=id2C,B1=1,B2=3,Mod9expr1));
evalf(subs(C=id2C,B1=1,B2=3,Mod9expr2));
evalf(subs(C=id2C,B1=1,B2=3,Mod9expr3));Identity 3:id3C:=findC([0,0,-1,-1,-1,-1,0,0,0]);
evalf(subs(C=id3C,B1=2,B2=3,Mod9expr1));
evalf(subs(C=id3C,B1=2,B2=3,Mod9expr2));
evalf(subs(C=id3C,B1=2,B2=3,Mod9expr3));Now we proceed to the general search.
We do not print the output of the following commands, since they are quite big.Mod9cval:=solve(Mod9expr1=0,C):
Mod9eq1:=simplify(subs(C=Mod9cval,Mod9expr2)):
Mod9eq2:=simplify(subs(C=Mod9cval,Mod9expr3)):for b1 from -40 to 40 do:
for b2 from -40 to 40 do:
tempCval:=solve(subs(B1=b1,B2=b2,Mod9expr1)=0,C):
Mod9eq2eval:=evalf(subs(B1=b1,B2=b2,Mod9eq1)):
Mod9eq3eval:=evalf(subs(B1=b1,B2=b2,Mod9eq2)):
if abs(Mod9eq2eval)<10^(-3) and abs(Mod9eq3eval)<10^(-3) then:
print("C", identify(evalf(tempCval)), " B1", b1," B2", b2):
fi:
end:
end: Dilogarithm verification:evalf((L(1) - L(Mod9Q1)) + (L(1)-L(Mod9Q2^3))/3 - 2*Pi^2/27);