You are not logged in.
BeginPackage["DedekindEta`"];
Unprotect[DedekindEta, LogDedekindEta, DedekindSum, WeberF];
ClearAll[LogDedekindEta, DedekindSum, WeberF]
Begin["`Private`"];
DedekindEta /: HoldPattern[Derivative[1][DedekindEta]] :=
DedekindEta[#]*LogDedekindEta'[#] &;
LogDedekindEta[z_] /; InexactNumberQ[z] && Im[z] > 0 :=
Block[{t = z, n, pr, f},
n = Round[Re[t]]; t -= n; pr = Precision[t];
f = If[7*Im[z] < 6,
n*Pi*I/12 + LogDedekindEta[SetPrecision[-1/t, pr]] - Log[-I*t]/2,
t = z*Pi*I;
t/12 + Log[QPochhammer[Exp[2*t]]],
Null
];
f /; NumberQ[f]
];
LogDedekindEta /: HoldPattern[Derivative[1][LogDedekindEta]] :=
With[{m = ModularLambda[#]},
With[{t = EllipticK[m]},
I*t*(EllipticE[m] + t*(m - 2)/3)/Pi
]
] &;
LogDedekindEta /: HoldPattern[E^LogDedekindEta[z_]] := DedekindEta[z];
LogDedekindEta /: MakeBoxes[LogDedekindEta[z_], TraditionalForm] :=
RowBox[{InterpretationBox["log\[Eta]", LogDedekindEta, Editable -> False, Selectable -> False, Tooltip -> "LogDedekindEta"], "(", ToBoxes[z], ")"}];
dedekindSum[x_Rational] :=
Module[{h, k, kk, s = 1, s1 = 0, s2, p = 1, pp = 0, r, nexth, a},
{h, k} = NumeratorDenominator[x];
h = Mod[h, k];
kk = k;
s2 = h;
While[h != 0,
{a, nexth} = QuotientRemainder[kk, h];
If[h == 1, s2 += p*s];
s1 += a*s;
s = -s;
kk = h; h = nexth;
r = a*p + pp; pp = p; p = r;
];
If[s < 0, s1 -= 3];
(k*s1 + s2)/(12*k)
];
dedekindSum[x_Integer] := 0;
dedekindSum[___] := $Failed;
DedekindSum[h_Integer, k_Integer] :=
With[{s = dedekindSum[h/k]},
s /; NumberQ[s]
];
DedekindSum[a_Integer, b_Integer, c_Integer] :=
Module[{g, rb, rc, d, s},
{g, {rb, rc}} = ExtendedGCD[b, c];
d = GCD[a, g];
s = dedekindSum[a*g*rb/(c*d)];
d*s /; NumberQ[s]
];
WeberF[z_, p_, n_] /; InexactNumberQ[z] && Im[z] > 0 && PrimeQ[p] && IntegerQ[n] && n >= 0 && n < p :=
With[{f = Exp[-Pi*I*n/(12*p)]*DedekindEta[(z + n)/p]/DedekindEta[z]},
f /; NumberQ[f]
];
WeberF[z_, p_] /; InexactNumberQ[z] && Im[z] > 0 && p >= 2 && PrimeQ[p] :=
With[{f = Sqrt[p]*DedekindEta[p*z]/DedekindEta[z]},
f /; NumberQ[f]
];
WeberF[z_, 1] /; InexactNumberQ[z] && Im[z] > 0 :=
With[{f = WeberF[z, 2, 0]},
f /; NumberQ[f]
];
WeberF[z_] /; InexactNumberQ[z] && Im[z] > 0 :=
With[{f = WeberF[z, 2, 1]},
f /; NumberQ[f]
];
End[];
SetAttributes[{LogDedekindEta, DedekindSum, WeberF}, {Listable, NumericFunction, ReadProtected}];
SetAttributes[DedekindSum, {NHoldAll}];
SetAttributes[WeberF, {NHoldRest}];
Protect[DedekindEta, LogDedekindEta, DedekindSum, WeberF];
EndPackage[];Test cases for Dedekind sum:
Module[{p, q, u, v, s, t, a, b, p1, q1, u1, v1, r, g},
r = 10^10;
p = RandomInteger[{-r, r}];
q = 1 + RandomInteger[r];
{g, {p1, q1}} = ExtendedGCD[p, q];
{p, q} /= g;
u = RandomInteger[{-r, r}];
v = 1 + RandomInteger[r];
{g, {u1, v1}} = ExtendedGCD[u, v];
{u, v} /= g;
(* Addition theorem of Dedekind sum *)
s = DedekindSum[p, q] + DedekindSum[u, v];
t = p*v + q*u;
If[t == 0,
Assert[s == 0],
a = p*u1 - q*v1;
b = u*p1 - v*q1;
r = q*v*t;
r = (1/12)*(q^2 + v^2 + t^2)/r - Sign[r]/4;
Assert[Mod[a*b - 1, Abs[t]] == 0];
Assert[s == DedekindSum[a, t] + r == DedekindSum[b, t] + r]
];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Identities:
DedekindSum[h, k] == DedekindSum[1, h, k]
DedekindSum[h, k] == 0 /; Mod[h^2 + 1, k] == 0
DedekindSum[a, b, c] == DedekindSum[b, a, c]
DedekindSum[-a, b, c] == -DedekindSum[a, b, c]
DedekindSum[a, b, -c] == -DedekindSum[a, b, c]
DedekindSum[a*d, b*d, c*d] == d*DedekindSum[a, b, c]
DedekindSum[a, b*d, c*d] == DedekindSum[a, b, c] /; GCD[a, d] == 1
DedekindSum[a*d, b*d, c] == DedekindSum[a, b, c] /; GCD[c, d] == 1
DedekindSum[a, b, c] + DedekindSum[b, c, a] + DedekindSum[c, a, b] == (1/12)*(a^2 + b^2 + c^2)/d - Sign[d]/4 /; d == a*b*c != 0 && GCD[a, b] == 1 && GCD[b, c] == 1 && GCD[c, a] == 1
WeberF[z + k, p] == Exp[(p - 1)*k*I*Pi/12]*WeberF[z, p] /; p >= 2
WeberF[z + k, p, n] == Exp[(1/p - 1)*k*I*Pi/12]*WeberF[z, p, Mod[n + k, p]] /; p >= 2
WeberF[-1/z, p] == WeberF[z, p, 0] /; p >= 2
WeberF[-1/z, p, 0] == WeberF[z, p] /; p >= 2
WeberF[-1/z, p, n] == Exp[-I*Pi*DedekindSum[n, p]]*WeberF[z, p, PowerMod[-n, -1, p]] /; p >= 2 && 1 <= n <= p - 1BeginPackage["KleinInvariantJ`"];
Unprotect[InverseKleinInvariantJ];
ClearAll[InverseKleinInvariantJ];
Begin["`Private`"];
InverseKleinInvariantJ[j_] :=
Module[{a},
Which[
PossibleZeroQ[j],
Return[(I*Sqrt[3] - 1)/2 + 4*I*(-j)^(1/3)*(Pi/(Gamma[1/3]^2))^3, Module],
PossibleZeroQ[a = j - 1],
Return[I + 8*I*Pi^2*Sqrt[a/3]/Gamma[1/4]^4, Module],
InexactNumberQ[j],
a = Tan[ArcCsc[Sqrt[j]]/3];
a *= 2/(a + Sqrt[3]);
a = I*ArithmeticGeometricMean[1, Sqrt[1 - a]]/ArithmeticGeometricMean[1, Sqrt[a]];
If[InexactNumberQ[a], Return[a, Module]],
Head[j] === DirectedInfinity,
Return[I*Infinity, Module]
];
Null /; False
];
InverseKleinInvariantJ /:
HoldPattern[KleinInvariantJ[InverseKleinInvariantJ[j_]]] := j;
InverseKleinInvariantJ /:
HoldPattern[Derivative[1][InverseKleinInvariantJ]] :=
I/(4*Sqrt[3]*Pi*DedekindEta[InverseKleinInvariantJ[#]]^4*#^(2/3)*Sqrt[# - 1]) &
InverseKleinInvariantJ /:
MakeBoxes[InverseKleinInvariantJ[j_], TraditionalForm] :=
RowBox[{InterpretationBox[SuperscriptBox["J", RowBox[{"-", "1"}]], InverseKleinInvariantJ, Editable -> False, Selectable -> False, Tooltip -> "InverseKleinInvariantJ"], "(", ToBoxes[j], ")"}];
End[];
SetAttributes[InverseKleinInvariantJ, {Listable, NumericFunction, ReadProtected}];
Protect[InverseKleinInvariantJ];
EndPackage[];BeginPackage["EisensteinE`"];
Unprotect[EisensteinE];
ClearAll[EisensteinE];
Begin["`Private`"];
EisensteinE[2, q_] /; InexactNumberQ[q] && Abs[q] < 1 :=
Block[{q2 = q^2, qn = -1, qn2 = 1, a = 1, b = 1, k = 1, t},
While[qn *= q2; qn2 *= qn; k += 2; t = a + k^3*qn2; b += k*qn2; t != a, a = t];
t / b
];
EisensteinE[n_, q_] /; EvenQ[n] && n >= 4 && InexactNumberQ[q] && Abs[q] < 1 :=
Block[{q2 = q^2, qn = 1, k = 1, r = n - 1, s = 0, t},
While[qn *= q2; t = s + k^r*qn/(1 - qn); k++; t != s, s = t];
1 - 2*n*t/BernoulliB[n]
];
End[];
SetAttributes[EisensteinE, {Listable, NHoldFirst, NumericFunction, ReadProtected}];
Protect[EisensteinE];
EndPackage[];Test cases:
Module[{q},
q = RandomReal[WorkingPrecision -> 50];
Assert[ResourceFunction["EisensteinE"][#, q] == EisensteinE[#, q]] & /@ Range[2,8,2];
] & // Block[{$AssertFunction = Automatic}, #[]] &;BeginPackage["Lemniscate`"];
Unprotect[
LemniscateSin, LemniscateCos, LemniscateTan, LemniscateCot, LemniscateSinh, LemniscateCosh, LemniscateSinhPrime, LemniscateCoshPrime, LemniscateTanh, LemniscateCoth,
InverseLemniscateSin, InverseLemniscateCos, InverseLemniscateTan, InverseLemniscateCot, InverseLemniscateSinh, InverseLemniscateCosh, InverseLemniscateTanh, InverseLemniscateCoth
];
ClearAll[
LemniscateSin, LemniscateCos, LemniscateTan, LemniscateCot, LemniscateSinh, LemniscateCosh, LemniscateSinhPrime, LemniscateCoshPrime, LemniscateTanh, LemniscateCoth,
InverseLemniscateSin, InverseLemniscateCos, InverseLemniscateTan, InverseLemniscateCot, InverseLemniscateSinh, InverseLemniscateCosh, InverseLemniscateTanh, InverseLemniscateCoth
];
Begin["`Private`"];
LemniscateSin[z_] /; (z == 0) := z;
LemniscateSin[z_?InexactNumberQ] := JacobiSN[z, -1];
LemniscateSin /: HoldPattern[LemniscateSin'] :=
LemniscateCos[#]*(1 + LemniscateSin[#]^2) &;
LemniscateSin /: MakeBoxes[LemniscateSin[z_], TraditionalForm] :=
RowBox[{InterpretationBox["sl", LemniscateSin, Editable -> False, Selectable -> False, Tooltip -> "LemniscateSin"], "(", ToBoxes[z], ")"}]
LemniscateCos[z_] /; (z == 0) := 1 - z^2;
LemniscateCos[z_?InexactNumberQ] := JacobiCD[z, -1]
LemniscateCos /: HoldPattern[LemniscateCos'] :=
-LemniscateSin[#]*(1 + LemniscateCos[#]^2) &;
LemniscateCos /: MakeBoxes[LemniscateCos[z_], TraditionalForm] :=
RowBox[{InterpretationBox["cl", LemniscateCos, Editable -> False, Selectable -> False, Tooltip -> "LemniscateCos"], "(", ToBoxes[z], ")"}]
LemniscateTan[z_] /; (z == 0) := z;
LemniscateTan[z_?InexactNumberQ] :=
JacobiSC[z, -1]*Sqrt[JacobiCD[z, -1]];
LemniscateTan /: HoldPattern[LemniscateTan'] :=
LemniscateCot[#]^3 &;
LemniscateTan /: MakeBoxes[LemniscateTan[z_], TraditionalForm] :=
RowBox[{InterpretationBox["tl", LemniscateTan, Editable -> False, Selectable -> False, Tooltip -> "LemniscateTan"], "(", ToBoxes[z], ")"}]
LemniscateCot[z_] /; (z == 0) := 1 + z^4/4;
LemniscateCot[z_?InexactNumberQ] :=
JacobiNC[z, -1]*Sqrt[JacobiCD[z, -1]];
LemniscateCot /: HoldPattern[LemniscateCot'] :=
LemniscateTan[#]^3 &;
LemniscateCot /: MakeBoxes[LemniscateCot[z_], TraditionalForm] :=
RowBox[{InterpretationBox["ctl", LemniscateCot, Editable -> False, Selectable -> False, Tooltip -> "LemniscateCot"], "(", ToBoxes[z], ")"}]
LemniscateSinh[z_] /; (z == 0) := z;
LemniscateSinh[z_?InexactNumberQ] :=
JacobiSN[z, 1/2]*JacobiDC[z, 1/2];
LemniscateSinh /: HoldPattern[LemniscateSinh'] :=
LemniscateSinhPrime;
LemniscateSinh /: MakeBoxes[LemniscateSinh[z_], TraditionalForm] :=
RowBox[{InterpretationBox["slh", LemniscateSinh, Editable -> False, Selectable -> False, Tooltip -> "LemniscateSinh"], "(", ToBoxes[z], ")"}]
LemniscateSinhPrime[z_] /; (z == 0) := 1 + z^4/2;
LemniscateSinhPrime[z_?InexactNumberQ] :=
With[{t = JacobiCN[z, 1/2]^2},
(t + t^(-1))/2
];
LemniscateSinhPrime /: HoldPattern[LemniscateSinhPrime'] :=
2*LemniscateSinh[#]^3 &;
LemniscateSinhPrime /: MakeBoxes[LemniscateSinhPrime[z_], TraditionalForm] :=
RowBox[{InterpretationBox[SuperscriptBox["slh", "\[Prime]"], LemniscateSinhPrime, Editable -> False, Selectable -> False, Tooltip -> "LemniscateSinhPrime"], "(", ToBoxes[z], ")"}]
LemniscateCosh[z_] /; (z == 0) := z^(-1);
LemniscateCosh[z_?InexactNumberQ] :=
JacobiNS[z, 1/2]*JacobiCD[z, 1/2];
LemniscateCosh /: HoldPattern[LemniscateCosh'] :=
LemniscateCoshPrime;
LemniscateCosh /: MakeBoxes[LemniscateCosh[z_], TraditionalForm] :=
RowBox[{InterpretationBox["clh", LemniscateCosh, Editable -> False, Selectable -> False, Tooltip -> "LemniscateCosh"], "(", ToBoxes[z], ")"}]
LemniscateCoshPrime[z_] /; (z == 0) := -z^(-2);
LemniscateCoshPrime[z_?InexactNumberQ] :=
With[{t = JacobiSD[z, 1/2]^2},
-(t^(-1) + t/4)
];
LemniscateCoshPrime /: HoldPattern[LemniscateCoshPrime'] :=
2*LemniscateCosh[#]^3 &;
LemniscateCoshPrime /: MakeBoxes[LemniscateCoshPrime[z_], TraditionalForm] :=
RowBox[{InterpretationBox[SuperscriptBox["clh", "\[Prime]"], LemniscateCoshPrime, Editable -> False, Selectable -> False, Tooltip -> "LemniscateCoshPrime"], "(", ToBoxes[z], ")"}]
LemniscateTanh[z_] /; (z == 0) := z;
LemniscateTanh[z_?InexactNumberQ] :=
With[{t = JacobiSD[z, 1/2]},
t/Sqrt[1 + t^4/4]
];
LemniscateTanh /: HoldPattern[LemniscateTanh'] :=
LemniscateCoth[#]^3 &;
LemniscateTanh /: MakeBoxes[LemniscateTanh[z_], TraditionalForm] :=
RowBox[{InterpretationBox["tlh", LemniscateTanh, Editable -> False, Selectable -> False, Tooltip -> "LemniscateTanh"], "(", ToBoxes[z], ")"}]
LemniscateCoth[z_] /; (z == 0) := 1 - z^4/4;
LemniscateCoth[z_?InexactNumberQ] :=
JacobiCD[z, 1/2]*JacobiND[z, 1/2]/Sqrt[1 + JacobiSD[z, 1/2]^4/4];
LemniscateCoth /: HoldPattern[LemniscateCoth'] :=
-LemniscateTanh[#]^3 &;
LemniscateCoth /: MakeBoxes[LemniscateCoth[z_], TraditionalForm] :=
RowBox[{InterpretationBox["ctlh", LemniscateCoth, Editable -> False, Selectable -> False, Tooltip -> "LemniscateCoth"], "(", ToBoxes[z], ")"}]
InverseLemniscateSin[z_] := With[{t = z*Hypergeometric2F1[1/4, 1/2, 5/4, z^4]},
t /; FreeQ[t, Hypergeometric2F1]
];
InverseLemniscateSin /: HoldPattern[InverseLemniscateSin'] :=
(1 - #^4)^(-1/2) &;
InverseLemniscateSin /: HoldPattern[LemniscateSin[InverseLemniscateSin[z_]]] :=
z;
InverseLemniscateSin /: HoldPattern[LemniscateCos[InverseLemniscateSin[z_]]] :=
Sqrt[1 - z^2]/Sqrt[1 + z^2];
InverseLemniscateSin /: HoldPattern[LemniscateTan[InverseLemniscateSin[z_]]] :=
z/(1 - z^4)^(1/4);
InverseLemniscateSin /: HoldPattern[LemniscateCot[InverseLemniscateSin[z_]]] :=
(1 - z^4)^(-1/4);
InverseLemniscateCos[z_] := With[{t = InverseLemniscateSin[1] - InverseLemniscateSin[z]},
t /; FreeQ[t, InverseLemniscateSin]
];
InverseLemniscateCos /: HoldPattern[InverseLemniscateCos'] :=
-(1 - #^4)^(-1/2) &;
InverseLemniscateCos /: HoldPattern[LemniscateSin[InverseLemniscateCos[z_]]] :=
Sqrt[1 - z^2]/Sqrt[1 + z^2];
InverseLemniscateCos /: HoldPattern[LemniscateCos[InverseLemniscateCos[z_]]] :=
z;
InverseLemniscateCos /: HoldPattern[LemniscateTan[InverseLemniscateCos[z_]]] :=
Sqrt[1 - z^2]/(Sqrt[2]*Sqrt[z]);
InverseLemniscateCos /: HoldPattern[LemniscateCot[InverseLemniscateCos[z_]]] :=
Sqrt[1 + z^2]/(Sqrt[2]*Sqrt[z]);
InverseLemniscateTan[z_] := With[{t = z*Hypergeometric2F1[1/4, 3/4, 5/4, -z^4]},
t /; FreeQ[t, Hypergeometric2F1]
];
InverseLemniscateTan /: HoldPattern[InverseLemniscateTan'] :=
(1 + #^4)^(-3/4) &;
InverseLemniscateCot[z_] := With[{t = z*(1 - z^(-4))^(1/4)*Hypergeometric2F1[1/4, 3/4, 5/4, 1 - z^4]},
t /; FreeQ[t, Hypergeometric2F1]
];
InverseLemniscateCot /: HoldPattern[InverseLemniscateCot'] :=
(#^4 - 1)^(-3/4) &;
InverseLemniscateSinh[z_] := With[{t = z*Hypergeometric2F1[1/4, 1/2, 5/4, -z^4]},
t /; FreeQ[t, Hypergeometric2F1]
]
InverseLemniscateSinh /: HoldPattern[InverseLemniscateSinh'] :=
(1 + #^4)^(-1/2) &
InverseLemniscateSinh /: HoldPattern[LemniscateSinh[InverseLemniscateSinh[z_]]] :=
z;
InverseLemniscateSinh /: HoldPattern[LemniscateCosh[InverseLemniscateSinh[z_]]] :=
z^(-1);
InverseLemniscateSinh /: HoldPattern[LemniscateSinhPrime[InverseLemniscateSinh[z_]]] :=
Sqrt[1 + z^4];
InverseLemniscateSinh /: HoldPattern[LemniscateCoshPrime[InverseLemniscateSinh[z_]]] :=
Sqrt[1 + z^4]/z^2;
InverseLemniscateSinh /: HoldPattern[LemniscateTanh[InverseLemniscateSinh[z_]]] :=
z/(1 + z^4)^(1/4);
InverseLemniscateSinh /: HoldPattern[LemniscateCoth[InverseLemniscateSinh[z_]]] :=
(1 + z^4)^(-1/4);
InverseLemniscateCosh[z_] := With[{t = 2*InverseLemniscateSinh[1] - InverseLemniscateSinh[z]},
t /; FreeQ[t, InverseLemniscateSinh]
];
InverseLemniscateCosh /: HoldPattern[InverseLemniscateCosh'] :=
(1 + #^4)^(-1/2) &
InverseLemniscateCosh /: HoldPattern[LemniscateSinh[InverseLemniscateCosh[z_]]] :=
z^(-1);
InverseLemniscateCosh /: HoldPattern[LemniscateCosh[InverseLemniscateCosh[z_]]] :=
z;
InverseLemniscateCosh /: HoldPattern[LemniscateSinhPrime[InverseLemniscateCosh[z_]]] :=
-Sqrt[1 + z^4]/z^2;
InverseLemniscateCosh /: HoldPattern[LemniscateCoshPrime[InverseLemniscateCosh[z_]]] :=
-Sqrt[1 + z^4];
InverseLemniscateCosh /: HoldPattern[LemniscateTanh[InverseLemniscateCosh[z_]]] :=
(1 + z^4)^(-1/4);
InverseLemniscateCosh /: HoldPattern[LemniscateCoth[InverseLemniscateCosh[z_]]] :=
z/(1 + z^4)^(1/4);
InverseLemniscateTanh[z_] := With[{t = z*Hypergeometric2F1[1/4, 3/4, 5/4, z^4]},
t /; FreeQ[t, Hypergeometric2F1]
];
InverseLemniscateTanh /: HoldPattern[InverseLemniscateTanh'] :=
(1 + #^4)^(-3/4) &;
InverseLemniscateCoth[z_] := With[{t = InverseLemniscateTanh[1] - InverseLemniscateTanh[z]},
t /; FreeQ[t, InverseLemniscateTanh]
];
InverseLemniscateCoth /: HoldPattern[InverseLemniscateCoth'] :=
-(1 + #^4)^(-3/4) &;
End[];
SetAttributes[{
LemniscateSin, LemniscateCos, LemniscateTan, LemniscateCot, LemniscateSinh, LemniscateCosh, LemniscateSinhPrime, LemniscateCoshPrime, LemniscateTanh, LemniscateCoth,
InverseLemniscateSin, InverseLemniscateCos, InverseLemniscateTan, InverseLemniscateCot, InverseLemniscateSinh, InverseLemniscateCosh, InverseLemniscateTanh, InverseLemniscateCoth
}, {Listable, NumericFunction, ReadProtected}];
Protect[
LemniscateSin, LemniscateCos, LemniscateTan, LemniscateCot, LemniscateSinh, LemniscateCosh, LemniscateSinhPrime, LemniscateCoshPrime, LemniscateTanh, LemniscateCoth,
InverseLemniscateSin, InverseLemniscateCos, InverseLemniscateTan, InverseLemniscateCot, InverseLemniscateSinh, InverseLemniscateCosh, InverseLemniscateTanh, InverseLemniscateCoth
];
EndPackage[];Special values:
Module[{w1, w2, w3, g2, g3, Δ, p, ζ, case},
w1 = RandomReal[WorkingPrecision -> 50];
z = w1*RandomReal[WorkingPrecision -> 50];
p[z_] := WeierstrassP[z, {g2, g3}];
ζ[z_] := WeierstrassZeta[z, {g2, g3}];
case[t_] := (w3 = t*w1; w2 = -w1 - w3; {g2, g3} = WeierstrassInvariants[{w1, w3}]; Δ = g2^3 - 27*g3^2);
case[I]; (* Lemniscatic cases *)
Assert[ζ[w1] == Pi/(4*w1)];
Assert[ζ[w3] == -I*Pi/(4*w1)];
Assert[p[w1] == Gamma[1/4]^4/(32*Pi*w1^2)];
Assert[p[w2] == 0];
Assert[p[w3] == -Gamma[1/4]^4/(32*Pi*w1^2)];
Assert[g2 == Gamma[1/4]^8/(16*Pi*w1^2)^2];
Assert[g3 == 0];
Assert[Δ == Gamma[1/4]^24/(16*Pi*w1^2)^6];
Assert[g2^3/Δ == 1];
Assert[p[(2/5)*w1] - p[(4/5)*w1] == Sqrt[GoldenRatio]*Gamma[1/4]^4/(16*Pi*w1^2)];
With[{x = p[z], y = p'[z]},
Assert[p[(1 + I)*z] == I*(g2/x - 4*x)/8];
Assert[p'[(1 + I)*z] == -(1 + I)*y*(g2/x^2 + 4)/16];
Assert[p[(1 + I)*z/2] == -I*x*(1 + Sqrt[4 - g2/x^2]/2)];
Assert[z == x^(-1/2)*Hypergeometric2F1[1/4, 1/2, 5/4, g2/(4*x^2)]];
];
case[Exp[2*Pi*I/3]]; (* Equianharmonic cases *)
Assert[ζ[w1] == Pi/(2*Sqrt[3]*w1)];
Assert[ζ[w2] == Exp[2*Pi*I/3]*Pi/(2*Sqrt[3]*w1)];
Assert[ζ[w3] == Exp[-2*Pi*I/3]*Pi/(2*Sqrt[3]*w1)];
Assert[p[w1] == Gamma[1/3]^6/(2^(14/3)*(Pi*w1)^2)];
Assert[p[w2] == Exp[-2*Pi*I/3]*Gamma[1/3]^6/(2^(14/3)*(Pi*w1)^2)];
Assert[p[w3] == Exp[2*Pi*I/3]*Gamma[1/3]^6/(2^(14/3)*(Pi*w1)^2)];
Assert[g2 == 0];
Assert[g3 == Gamma[1/3]^18/(4*Pi*w1)^6];
Assert[Δ == -27*Gamma[1/3]^36/(4*Pi*w1)^12];
Assert[g2^3/Δ == 0];
Assert[p[(2/5)*w1] - p[(4/5)*w1] == Sqrt[3]*5^(-1/12)*(2 + Sqrt[5])^(1/6)*Gamma[1/3]^6/(4*Pi*w1)^2];
With[{x = p[z], y = p'[z]},
Assert[p[I*Sqrt[3]*z] == 3^(-1)*(g3/x^2 - x)];
Assert[p'[I*Sqrt[3]*z] == I*3^(-3/2)*y*(1 + 8*g3/(y^2 + g3))];
Assert[p[I*z/Sqrt[3]] == -2^(-1)*x^(-1/2)*g3^(1/2)*Csc[ArcCsc[2*x^(3/2)*g3^(-1/2)]/3]];
Assert[z == x^(-1/2)*Hypergeometric2F1[1/6, 1/2, 7/6, g3/(4*x^3)]];
];
case[I*Infinity]; (* Degenerate cases *)
Assert[ζ[w1] == Pi^2/(12*w1)];
Assert[p[w1] == (Pi/w1)^2/6];
Assert[g2 == (Pi/w1)^4/12];
Assert[g3 == (Pi/w1)^6/216];
Assert[Δ == 0];
Assert[p[(2/5)*w1] - p[(4/5)*w1] == 5^(-1/2)*(Pi/w1)^2]
case[I*Sqrt[2]];
Assert[ζ[w1] == (Gamma[1/8]*Gamma[3/8])^2/(384*Pi*w1) + Pi/(2^(5/2)*w1)];
Assert[ζ[w3] == I*(Sqrt[2]*(Gamma[1/8]*Gamma[3/8])^2/(384*Pi*w1) - Pi/(4*w1))];
Assert[p[w1] == (3 + Sqrt[2])*(Gamma[1/8]*Gamma[3/8])^2/(2^(13/2)*3*Pi*w1^2)];
Assert[p[w2] == -(3 - Sqrt[2])*(Gamma[1/8]*Gamma[3/8])^2/(2^(13/2)*3*Pi*w1^2)];
Assert[p[w3] == -(Gamma[1/8]*Gamma[3/8])^2/(2^5*3*Pi*w1^2)];
Assert[g2 == 5*(Gamma[1/8]*Gamma[3/8])^4/(2^11*3*(Pi*w1^2)^2)];
Assert[g3 == 7*(Gamma[1/8]*Gamma[3/8])^6/(2^16*(3*Pi*w1^2)^3)];
Assert[Δ == (Gamma[1/8]*Gamma[3/8])^12/(2^33*(Pi*w1^2)^6)];
Assert[g2^3/Δ == 125/27];
case[I*Sqrt[3]];
Assert[ζ[w1] == Gamma[1/3]^6/(2^(20/3)*Pi^2*w1) + Pi/(4*Sqrt[3]*w1)];
Assert[ζ[w3] == I*(Sqrt[3]*Gamma[1/3]^6/(2^(20/3)*Pi^2*w1) - Pi/(4*w1))];
Assert[p[w1] == (2*Sqrt[3] + 1)*Gamma[1/3]^6/(2^(20/3)*(Pi*w1)^2)];
Assert[p[w2] == -Gamma[1/3]^6/(2^(17/3)*(Pi*w1)^2)];
Assert[p[w3] == -(2*Sqrt[3] - 1)*Gamma[1/3]^6/(2^(20/3)*(Pi*w1)^2)];
Assert[g2 == 15*Gamma[1/3]^12/(2^(34/3)*(Pi*w1)^4)];
Assert[g3 == 11*Gamma[1/3]^18/(2^17*(Pi*w1)^6)];
Assert[Δ == 27*Gamma[1/3]^36/(2^32*(Pi*w1)^12)];
Assert[g2^3/Δ == 125/4];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Half period and Quarter period values:
Module[{w1, w2, w3, w4, g2, g3, Δ, e1, e2, e3, η1, η2, a, b, c, ϕ, p, ζ, σ},
On[Assert];
w1 = 1;
w2 = I*Sqrt[Tan[RandomReal[WorkingPrecision -> 50]*Pi/2]];
w3 = w1 + w2;
w4 = w1 - w2;
{g2, g3} = WeierstrassInvariants[{w1, w2}];
Δ = g2^3 - 27*g3^2;
p[z_] := WeierstrassP[z, {g2, g3}];
ζ[z_] := WeierstrassZeta[z, {g2, g3}];
σ[z__] := Times @@ WeierstrassSigma[{z}, {g2, g3}];
(* Half period values *)
{e1, e2, e3} = p[{w1, w2, w3}];
{η1, η2} = ζ[{w1, w2}];
Assert[4*#^3 - g2*# - g3 == 0] & /@ {e1, e2, e3};
Assert[e1 + e2 + e3 == 0];
Assert[g2 == 2*(e1^2 + e2^2 + e3^2) == -4*(e1*e2 + e1*e3 + e2*e3)];
Assert[g3 == 4*e1*e2*e3 == (4/3)*(e1^3 + e2^3 + e3^3)];
Assert[Δ == 16*(e1 - e2)^2*(e1 - e3)^2*(e2 - e3)^2];
Assert[Δ == #^2*(4*#^2 - 3*g2)^2] & /@ (Subtract @@@ Permutations[{e1, e2, e3}, {2}]);
Assert[4*(e1 - e2)*(e1 - e3) == 12*e1^2 - g2];
Assert[4*(e2 - e3)*(e2 - e1) == 12*e2^2 - g2];
Assert[4*(e3 - e1)*(e3 - e2) == 12*e3^2 - g2];
Assert[p'[#] == 0] & /@ {w1, w2, w3};
a = Sqrt[g2/3];
ϕ = ArcSin[g3/a^3]/3;
Assert[e1 == a*Cos[Pi/6 - ϕ]];
Assert[e2 == -a*Cos[Pi/6 + ϕ]];
Assert[e3 == -a*Sin[ϕ]];
Assert[ζ[w3] == η1 + η2];
Assert[ζ[w4] == η1 - η2];
Assert[σ[w1] == Sqrt[2]*E^((1/2)*w1*η1)/(12*e1^2 - g2)^(1/4)];
Assert[σ[w2] == I*Sqrt[2]*E^((1/2)*w2*η2)/(12*e2^2 - g2)^(1/4)];
(* Quarter period values *)
a = Sqrt[e1 - e2];
b = Sqrt[e1 - e3];
c = Sqrt[e3 - e2];
Assert[p[(1/2)*w1] == e1 + a*b];
Assert[p[(1/2)*w2] == e2 - a*c];
Assert[p[(1/2)*w3] == e3 - I*b*c];
Assert[p[(1/2)*w4] == e3 + I*b*c];
Assert[p[(1/2)*w1 + w2] == e1 - a*b];
Assert[p[(1/2)*w2 + w1] == e2 + a*c];
Assert[p'[(1/2)*w1] == -2*a*b*(a + b)];
Assert[p'[(1/2)*w2] == -2*I*a*c*(a + c)];
Assert[p'[(1/2)*w3] == 2*b*c*(b + I*c)];
Assert[p'[(1/2)*w4] == 2*b*c*(b - I*c)];
Assert[p'[(1/2)*w1 + w2] == 2*a*b*(a - b)];
Assert[p'[(1/2)*w2 + w1] == 2*I*a*c*(a - c)];
Assert[ζ[(1/2)*w1] == (1/2)*(ζ[w1] + (a + b))];
Assert[ζ[(1/2)*w2] == (1/2)*(ζ[w2] - I*(a + c))];
Assert[σ[(1/2)*w1] == E^((1/8)*w1*η1)/(2^(1/4)*(a + b)^(1/4)*(a*b)^(3/8))];
Assert[σ[(1/2)*w2] == I*E^((1/8)*w2*η2)/(2^(1/4)*(a + c)^(1/4)*(a*c)^(3/8))];
(* Dedekind eta function representations *)
With[{τ = w2/w1, η = DedekindEta},
Assert[e1 == (Pi/w1)^2*((η[τ/2]^8 + 8*η[2*τ]^8)/η[τ]^4)/6];
Assert[e2 == -(Pi/w1)^2*((η[τ/2]^8 + 32*η[2*τ]^8)/η[τ]^4)/12];
Assert[e3 == -(Pi/w1)^2*((η[τ/2]^8 - 16*η[2*τ]^8)/η[τ]^4)/12];
];
] & // Block[{$AssertFunction = Automatic}, #[]] &;One-third period values:
Module[{w1, w2, w3, w4, g2, g3, Δ, x, y, η1, η2, cbrtΔ, a0, a1, a2, a3, a4, p, ζ, σ},
w1 = 1;
w2 = I*Sqrt[Tan[RandomReal[WorkingPrecision -> 50]*Pi/2]];
w3 = w1 + w2;
w4 = w1 - w2;
{g2, g3} = WeierstrassInvariants[{w1, w2}];
Δ = g2^3 - 27*g3^2;
p[z_] := WeierstrassP[z, {g2, g3}];
ζ[z_] := WeierstrassZeta[z, {g2, g3}];
σ[z__] := Times @@ WeierstrassSigma[{z}, {g2, g3}];
{x, y} = Through[{p, p'}[(2/3)*{w1, w2, w3, w4}]];
Assert[#^4 - (1/2)*g2*#^2 - g3*# - (1/48)*g2^2 == 0] & /@ x;
Assert[#^8 - 8*g3*#^6 - (2/3)*Δ*#^4 - (1/27)*Δ^2 == 0] & /@ Join[y, -y];
Assert[g2 == Total[x^2]];
Assert[g3 == Total[x^3]/3 == Total[y^2]/8];
Assert[Δ == -I*3^(3/2)*(Times @@ y)];
cbrtΔ = CubeRoot[Δ];
a0 = (g2 - cbrtΔ)/3;
a1 = Sign[g3]*Sqrt[a0];
a2 = 2*If[a1 == 0, g2/Sqrt[3], g3/a1];
a3 = g2 - a0;
a4 = Sqrt[a2 - a0 - cbrtΔ];
Assert[x[[1]] == (a1 + Sqrt[a2 + a3])/2];
Assert[x[[2]] == (a1 - Sqrt[a2 + a3])/2];
Assert[x[[3]] == (-a1 - I*Sqrt[a2 - a3])/2];
Assert[x[[4]] == (-a1 + I*Sqrt[a2 - a3])/2];
Assert[y[[1]] == -cbrtΔ/(Sqrt[3]*Sqrt[a4 - a1])];
Assert[y[[2]] == -I*cbrtΔ/(Sqrt[3]*Sqrt[a4 + a1])];
{η1, η2} = ζ[{w1, w2}];
Assert[ζ[(2/3)*w1] == (2/3)*η1 + Sqrt[x[[1]]]/Sqrt[3]];
Assert[ζ[(2/3)*w2] == (2/3)*η2 - I*Sqrt[-x[[2]]]/Sqrt[3]];
Assert[σ[(2/3)*w1] == E^((2/9)*w1*η1)/(-y[[1]])^(1/3)];
Assert[σ[(2/3)*w2] == I*E^((2/9)*w2*η2)/(I*y[[2]])^(1/3)];
(* Dedekind eta function representations *)
With[{τ = w2/w1, η = DedekindEta},
Assert[x[[1]] == (Pi*(η[τ/3]^3 + 3*η[3*τ]^3)/(w1*η[τ]))^2/4];
Assert[x[[2]] == -(Pi*(η[τ/3]^3 + 9*η[3*τ]^3)/(w1*η[τ]))^2/12];
Assert[x[[3]] == -(Pi*(η[τ/3]^3 + I*3^(3/2)*η[3*τ]^3)/(w1*η[τ]))^2/12];
Assert[x[[4]] == -(Pi*(η[τ/3]^3 - I*3^(3/2)*η[3*τ]^3)/(w1*η[τ]))^2/12];
Assert[y[[1]] == -(Pi*η[τ]^3/(Sqrt[3]*w1*η[3*τ]))^3];
Assert[y[[2]] == -I*(Pi*η[τ]^3/(w1*η[τ/3]))^3];
Assert[y[[3]] == Sqrt[ I]*(Pi*η[τ]^3/(w1*η[(τ+1)/3]))^3];
Assert[y[[4]] == Sqrt[-I]*(Pi*η[τ]^3/(w1*η[(τ-1)/3]))^3];
];
] & // Block[{$AssertFunction = Automatic}, #[]] &;One-fifth period values:
Module[{w1, w2, g2, g3, Δ, p, ζ, σ},
w1 = 1;
w2 = I*Sqrt[Tan[RandomReal[WorkingPrecision -> 50]*Pi/2]];
{g2, g3} = WeierstrassInvariants[{w1, w2}];
Δ = g2^3 - 27*g3^2;
p[z_] := WeierstrassP[z, {g2, g3}];
ζ[z_] := WeierstrassZeta[z, {g2, g3}];
σ[z__] := Times @@ WeierstrassSigma[{z}, {g2, g3}];
Module[{x, y, sqrtx, t, a},
{x, y} = p[(2/5)*#] + {-1, 1}*p[(4/5)*#];
t = Sign[#];
sqrtx = Sqrt[t^2*x]/t;
t = Δ/x^6;
a = 3*y/x + (11 + t)/2;
Assert[x^12 - (12/5)*g2*x^10 + 2*Δ*x^6 + (1/5)*Δ^2 == 0];
Assert[y == x*Sqrt[5^3 + 22*t + t^2]/6];
Assert[p[(2/5)*#] == (y + x)/2 == x*(2/a + 17 + t)/12];
Assert[p[(4/5)*#] == (y - x)/2 == x*(2/a + 5 + t)/12];
Assert[p'[(2/5)*#] == -sqrtx^3*a^( 1/2)];
Assert[p'[(4/5)*#] == -sqrtx^3*a^(-1/2)];
Assert[ζ[(2/5)*#] == (2/5)*ζ[#] + sqrtx*(1 + 3*a)/(10*Sqrt[a])];
Assert[ζ[(4/5)*#] == (4/5)*ζ[#] - sqrtx*(3 - 1*a)/(10*Sqrt[a])];
Assert[σ[(2/5)*#] == E^((2/25)*#*ζ[#])*a^(-1/10)/sqrtx];
Assert[σ[(4/5)*#] == E^((8/25)*#*ζ[#])*a^( 1/10)/sqrtx];
(* https://mathworld.wolfram.com/AbelsImpossibilityTheorem.html *)
] & /@ {w1, w2};
(* Dedekind eta function representations *)
With[{τ = w2/w1, η = DedekindEta},
Assert[p[(2/5)*w1] - p[(4/5)*w1] == (Pi/w1)^2*η[τ]^5/(Sqrt[5]*η[5*τ])];
Assert[p[(2/5)*w2] - p[(4/5)*w2] == -(Pi/w1)^2*η[τ]^5/η[τ/5]];
];
(* Hypergeometric series representations *)
With[{t = g3*(3/g2)^(3/2), f = (Times @@ Hypergeometric2F1[1/6, 5/6, {4/5, 6/5}, #]) &},
Assert[p[(2/5)*w1] - p[(4/5)*w1] == 2*Sqrt[(3/5)*g2]/f[(1 - t)/2]];
Assert[p[(2/5)*w2] - p[(4/5)*w2] == -2*Sqrt[(3/5)*g2]/f[(1 + t)/2]];
];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Sigma function notation:
Module[{z, a, b, c0, c1, c2, w1, w2, g2, g3, d, p, ζ, σ},
w1 = 1;
w2 = I*Sqrt[Tan[RandomReal[WorkingPrecision -> 50]*Pi/2]];
{g2, g3} = WeierstrassInvariants[{w1, w2}];
z = w1*(2*RandomReal[WorkingPrecision -> 50] - 1);
a = w1*(2*RandomReal[WorkingPrecision -> 50] - 1);
b = w1*(2*RandomReal[WorkingPrecision -> 50] - 1);
c0 = a - b;
p[z_] := WeierstrassP[z, {g2, g3}];
ζ[z_] := WeierstrassZeta[z, {g2, g3}];
σ[z__] := Times @@ WeierstrassSigma[{z}, {g2, g3}];
With[{wp = p[c0], pd = p'[c0]},
c1 = InverseWeierstrassP[{(-wp + Sqrt[g2 - 3*wp^2])/2, pd}, {g2, g3}];
c2 = -c0 - c1;
];
Assert[ζ[z + a] - ζ[z - a] == 2*ζ[a] - σ[2*a, z, z]/σ[z + a, z - a, a, a]];
d = ζ[a - b];
Assert[ζ[z + a] - ζ[z + b] == d - σ[a - b, z + a + c1, z + b - c1]/σ[z + a, z + b, c1, c2]];
Assert[ζ[z - a] - ζ[z - b] == -d + σ[a - b, z - a - c1, z - b + c1]/σ[z - a, z - b, c1, c2]];
Assert[ζ[z + a] - ζ[z + b] == d - (σ[2*(z + b), z + a, z + a]/σ[z + b, z + b] + σ[2*(z + a), z + b, z + b]/σ[z + a, z + a])/(2*σ[a - b, 2*z + a + b])];
Assert[ζ[z - a] - ζ[z - b] == -d + (σ[2*(z - b), z - a, z - a]/σ[z - b, z - b] + σ[2*(z - a), z - b, z - b]/σ[z - a, z - a])/(2*σ[a - b, 2*z - a - b])];
d = ζ[a - b] + (3*p[a - b]^2 - g2/4)/p'[a - b];
Assert[ζ[z + a] - ζ[z + b] == d + σ[z - a + 2*b, z - b + 2*a]/σ[z + a, z + b, 2*(a - b)]];
Assert[ζ[z - a] - ζ[z - b] == -d - σ[z + a - 2*b, z + b - 2*a]/σ[z - a, z - b, 2*(a - b)]];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Relations to Jacobi elliptic functions:
Module[{z, p, q, s, c, d, n, ϕ1, ϕ2, ϕ3, e1, e2, e3, g2, g3, m, w1, w2, w3, eq, mod, nearestSqrt},
eq = Module[{l},
l = Cases[{##}, _?NumericQ];
Equal @@ SetAccuracy[l, Accuracy[l]]
] &;
mod[a_, b_, m_] := With[{z = Mod[a - b, m, -m/2]},
!NumericQ[z] || PossibleZeroQ[z]
];
nearestSqrt[z_, zt_] := If[ReIm[z].ReIm[zt] < 0, -z, z, z];
m = RandomReal[WorkingPrecision -> 50];
w1 = EllipticK[m];
w3 = I*EllipticK[1 - m];
w2 = -w1 - w3;
z = w1*RandomReal[WorkingPrecision -> 50];
e1 = (2 - m)/3;
e2 = -(1 - 2*m)/3;
e3 = -(1 + m)/3;
With[{t = m*(1 - m)},
g2 = (4/3)*(1 - t);
g3 = (4/27)*(1 - 2*m)*(2 + t);
];
Assert[WeierstrassInvariants[{w1, w3}] == {g2, g3}];
p = WeierstrassP[z/2, {g2, g3}];
q = WeierstrassPPrime[z/2, {g2, g3}];
s = q;
c = 1 - m - (p - e1)^2;
d = m*(m - 1) - (p - e2)^2;
n = m - (p - e3)^2;
ϕ1 = -2*ArcTan[2*Sqrt[m]*(p - e1)/q];
ϕ2 = -2*ArcTan[2*(p - e2)/q];
ϕ3 = -2*ArcTan[2*Sqrt[1 - m]*(p - e3)/q];
Assert[eq[JacobiSN[z, m], s/n]];
Assert[eq[JacobiCN[z, m], c/n]];
Assert[eq[JacobiDN[z, m], d/n]];
Assert[mod[JacobiAmplitude[z, m], ϕ2, 2*Pi]];
Assert[mod[JacobiAmplitude[w1 - z, m], Pi/2 - ϕ3, 2*Pi]];
Assert[eq[ϕ1, ArcSin[Sqrt[m]*Sin[ϕ2]], ArcTan[Sqrt[m]*Sin[ϕ3]/Sqrt[1 - m]]]];
Assert[mod[EllipticF[ϕ2, m], z, 2*w1]];
With[{ζ = WeierstrassZeta[#, {g2, g3}] &, σ = Times @@ WeierstrassSigma[{##}, {g2, g3}] &},
Assert[eq[JacobiZeta[ϕ2, m], ζ[z + w3] - ζ[w3] - z*ζ[w1]/w1]];
Assert[eq[EllipticE[m], ζ[w1] + w1*e1]];
Assert[eq[I*EllipticE[1 - m], -ζ[w3] - w3*e3]];
Assert[eq[NevilleThetaS[z, m], σ[z]*E^(-z^2*ζ[w1]/(2*w1))]];
Assert[eq[NevilleThetaC[z, m], σ[z + w1]/σ[w1]*E^(-z^2*ζ[w1]/(2*w1) - z*ζ[w1])]];
Assert[eq[NevilleThetaD[z, m], σ[z + w2]/σ[w2]*E^(-z^2*ζ[w1]/(2*w1) - z*ζ[w2])]];
Assert[eq[NevilleThetaN[z, m], σ[z + w3]/σ[w3]*E^(-z^2*ζ[w1]/(2*w1) - z*ζ[w3])]];
];
p = WeierstrassP[z, {g2, g3}];
q = WeierstrassPPrime[z, {g2, g3}];
s = 1;
c = p - e1;
d = p - e2;
n = p - e3;
Assert[eq[p, Csc[ϕ2]^2 + e3, Cot[ϕ2]^2 + e1]];
Assert[eq[q, -2*Csc[ϕ2]^2*Cot[ϕ2]*Sqrt[1 - m*Sin[ϕ2]^2]]];
Assert[eq[JacobiSN[z, m]^2, s/n]];
Assert[eq[JacobiCN[z, m]^2, c/n]];
Assert[eq[JacobiDN[z, m]^2, d/n]];
Assert[mod[ϕ2, -ArcCot[nearestSqrt[Sqrt[c], q/(c + Sqrt[1 - m])]], Pi]];
] & // Block[{$AssertFunction = Automatic},
Quiet[#[], {Power::infy, Infinity::indet}]
] &;(* https://resources.wolframcloud.com/PacletRepository/resources/JanM/Dixon/ *)
<< "JanM`Dixon`";
Hold[{f, g, u, v, α, β, s, c, m, x, y, z, b, γ, w1, w2, λ, λ0, λ1, λ2, λd, makeλ, κ, κ0, κ1, κ2, makeκ, μ, ϕ, τ, η, p, q, pr, sqrt},
pr = 50;
α = -RandomReal[WorkingPrecision -> pr];
b = (9*(1 - α + α^2))^(1/3);
γ = CubeRoot[1 + α^3];
sqrt[z_] := RandomChoice[{1, -1}]*Sqrt[z];
w1 = Exp[2*sqrt[-1]*Pi/3];
w2 = Conjugate[w1];
cis[z_] := Exp[I*Sign@Im[w1]*z];
λ0 = Gamma[1/3]^3/(2*Sqrt[3]*Pi); (* λ(0) *)
κ0 = 4*Pi^2/(3*Gamma[1/3]^3); (* κ(0) *)
λd = 2*Pi*3^(-3/2); (* λ(-1) *)
Assert[N[DixonF[λ0, 0] == κ0, pr]];
Assert[N[λ0 * κ0 == λd, pr]];
makeλ[α_] := Block[{x, y, z},
x = λ0*Hypergeometric2F1[1/3, 1/3, 2/3, -α^3];
y = κ0*α*Hypergeometric2F1[2/3, 2/3, 4/3, -α^3];
z = x + y;
If[PossibleZeroQ[Im[α]], z = Re[z]];
{z, w1*x + w2*y, w2*x + w1*y}
];
makeκ[α_] := Block[{x, y, z},
x = κ0*Hypergeometric2F1[-1/3, 2/3, 1/3, -α^3];
y = λ0*α^2*Hypergeometric2F1[1/3, 4/3, 5/3, -α^3]/2;
z = x + y;
If[PossibleZeroQ[Im[α]], z = Re[z]];
α + {x + y, w2*x + w1*y, w1*x + w2*y}
];
{λ, λ1, λ2} = Quiet@makeλ[α];
{κ, κ1, κ2} = Quiet@makeκ[α];
Which[
VectorQ[{λ, κ}, NumericQ] (* γ != 0 *),
Assert[λ + λ1 + λ2 == 0];
With[{m = 3*ReIm[{λ, λ1, λ2}].Inverse@ReIm@DixonPeriods[α]},
Assert[MatrixQ[m, # == Round[#] &]];
Assert[VectorQ[Cross @@ Transpose[m], #^2 == 1 &]];
],
α < 0,
With[{b = 1 + α^(-3)},
{λ, κ} = λd*{
-Hypergeometric2F1[1/3, 2/3, 1, b]/α,
Hypergeometric2F1[2/3, 4/3, 1, b]/α^2
};
κ += α;
]
];
{u, v} = RandomReal[{0, λ}, 2, WorkingPrecision -> pr];
x = DixonSM[0, α];
y = DixonCM[0, α];
Assert[x == 0];
Assert[y == 1];
x = DixonSM[λ, α];
y = DixonCM[λ, α];
Assert[x == 1];
Assert[y == 0];
x = DixonSM[λ/2, α];
y = DixonCM[λ/2, α];
Assert[x == y];
Assert[2*x^3 == 1 + 3*α*x^2];
x = DixonSM[λ/3, α];
y = DixonCM[λ/3, α];
Assert[x^2*y + x - y^2 == 0];
Assert[y^2*x + y - x^2 == b*x*y];
Assert[x^3 + 1 == b*x];
Assert[y^3 + 1 == b*y^2];
x = DixonSM[λ/4, α];
y = DixonCM[λ/4, α];
Assert[x + x^3 == y^3*(1 - x)];
Assert[DixonF[λ, α] == κ];
If[α != -1,
Assert[DixonF[λ1, α] == κ1];
Assert[DixonF[λ2, α] == κ2];
];
s = DixonSM[u, α];
c = DixonCM[u, α];
Assert[s^3 + c^3 == 1 + 3*α*s*c];
Assert[DixonSM[λ - u, α] == c];
Assert[DixonCM[λ - u, α] == s];
Assert[DixonSM[u + λ, α] == DixonCM[-u, α] == 1/c];
Assert[DixonCM[u + λ, α] == DixonSM[-u, α] == -s/c];
Assert[DixonSM[u - λ, α] == DixonCM[-u - λ, α] == -c/s];
Assert[DixonCM[u - λ, α] == DixonSM[-u - λ, α] == 1/s];
Assert[DixonSM[2*u, α] == s*(1 + c^3)/(c*(1 + s^3))];
Assert[DixonCM[2*u, α] == (c^3 - s^3)/(c*(1 + s^3))];
Assert[DixonSM[3*u, α] == s*c*(1 + s^3 + c^3 + s^6 + c^6 - s^3*c^3)/(c^3 - s^6 + 3*s^3*c^3 + c^6*s^3)];
Assert[DixonCM[3*u, α] == -(s^3 - c^6 + 3*s^3*c^3 + s^6*c^3)/(c^3 - s^6 + 3*s^3*c^3 + c^6*s^3)];
x = Null;
With[{v = λ1 - λ2}, If[NumericQ[v],
Assert[DixonSM[u + v, α] == w2*s];
Assert[DixonCM[u + v, α] == w1*c];
x = -3^(-3/2)*Pi*v/(w1 - w2) + (1/2)*α^2*HypergeometricPFQ[{1, 1, 1}, {4/3, 5/3}, -α^3];
]];
(* Power series at α = -1 *)
f[α_] :=
Module[{x, y, z, n, b, c, sum},
b = c = 1 + SetAccuracy[α, Accuracy[α] + 5];
x = 0;
y = 1/3;
sum = x;
n = 2;
While[
z = sum;
sum += c*y;
z != sum,
c *= b;
z = ((1 - 3*n + 3*n^2)*y - (n - 1)^2*x)/(3*n^2);
x = y;
y = z;
n++;
];
SetAccuracy[sum, Accuracy[sum] - 5]
];
(* Differential equations:
w(α) + 7*α*w'(α) + 6*α^2*w''(α) + (1 + α^3)*w'''(α) = 0
Wronskian(λ(α), λ'(α), w(α)) = constant * (1 + α^3)^(-2)
*)
If[Abs[α + 1] < 1,
y = f[α] + (Pi/Sqrt[3] - Sqrt[3]*PolyGamma[1, 1/3]/(2*Pi))*λ;
If[NumericQ[x], Assert[x == y], x = y];
];
With[{b = SetPrecision[α, Infinity], μ = SetPrecision[λ, Infinity]},
y = NIntegrate[Log[DixonSM[u, b]/u], {u, 0, μ}, WorkingPrecision -> pr] + λ*(Log[λ] - 1)
If[NumericQ[x], Assert[x == y]];
];
Module[{s1, s2, s3, s4, c1, c2, c3, c4, f, g, z},
(* Addition theorem *)
{s1, s2, s3, s4} = DixonSM[{u, v, u + v, u - v}, α];
{c1, c2, c3, c4} = DixonCM[{u, v, u + v, u - v}, α];
f = DixonF[#, α] &;
g = Times @@ DixonG[{##}, α] &;
z = c2 + s1*c1*s2^2;
Assert[s3*z == s1 + s2*c2*c1^2];
Assert[s4*z == s1*c2^2 - s2*c1^2];
Assert[c3*z == c1*c2^2 - s1^2*s2];
Assert[c4*z == c1 + s2*c2*s1^2];
Assert[s3*s4*z == s1^2*c2 - s2^2*c1];
Assert[s3*c4*z == s2 + s1*c1*c2^2];
Assert[c3*s4*z == s1*c1 - s2*c2];
Assert[c3*c4*z == c1^2*c2 - s1*s2^2];
Assert[f[u + v] == f[u] + f[v] + s1*s2*c3];
Assert[f[u - v] == f[u] - f[v] - c1*s2*s4];
Assert[f[-u - v] == -f[u] - f[v] + (s1^3 - s2^3)/(s1*c1 - s2*c2)];
Assert[g[u + v, u - v, 0, 0]/g[u, v]^2 == E^(-α*v)*z];
];
With[{v = sqrt[-1]*Log[RandomReal[WorkingPrecision -> pr]]},
Assert[Abs[DixonCM[v, α]] == 1];
];
Assert[Derivative[1, 0][DixonSM][u, α] == c^2 - α*s];
Assert[Derivative[1, 0][DixonCM][u, α] == -s^2 + α*c];
f[u_, α_] := DixonSM[u, α];
x = f[u, α];
y = Derivative[1, 0][f][u, α];
Assert[y^3 - 3*α*x*y^2 - x^6 + (2 + 4*α^3)*x^3 - 1 == 0];
If[α < 0 && x > 0, Assert[y*Hypergeometric2F1[1/3, 2/3, 1/2, (1 - x^3)^2/(4*(α*x)^3)] == -α*x]];
f[u_, α_] := DixonCM[u, α];
x = f[u, α];
y = Derivative[1, 0][f][u, α];
Assert[y^3 + 3*α*x*y^2 + x^6 - (2 + 4*α^3)*x^3 + 1 == 0];
If[α < 0 && x > 0, Assert[y*Hypergeometric2F1[1/3, 2/3, 1/2, (1 - x^3)^2/(4*(α*x)^3)] == α*x]];
f[u_, α_] := DixonSM[u, α] + DixonCM[u, α];
x = f[u, α];
y = Derivative[1, 0][f][u, α];
Assert[3*y^2 + (x + α)*(x^3 - 3*α*x^2 - 4) == 0];
f[u_, α_] := DixonSM[u, α]*DixonCM[u, α];
x = f[u, α];
y = Derivative[1, 0][f][u, α];
Assert[y^2 + 4*x^3 == (1 + 3*α*x)^2];
f[u_, α_] := If[PossibleZeroQ[1 + α], -1/3, 1/(DixonSM[u, α] + DixonCM[u, α] + α)];
x = f[u, α];
y = Derivative[1, 0][f][u, α];
Assert[3*y^2 == 4*(1 + α^3)*x^3 - (1 - 3*α*x)^2];
If[α < 0,
f[u_, α_] := With[{s = DixonSM[u, α], c = DixonCM[u, α]}, ArcTan[(1 - c - 2*α*s)/(Sqrt[3]*(1 + c))]];
x = f[u, α];
y = Derivative[1, 0][f][u, α];
Assert[y*Hypergeometric2F1[1/3, 2/3, 1/2, (1 + α^(-3))*Sin[x]^2] == -Sqrt[3]*α/2];
];
If[α < 2,
f[u_, α_] := With[{s = DixonSM[u, α], c = DixonCM[u, α]}, ArcTan[(1 - c + 2*s)/(Sqrt[3]*(1 + c))]];
x = f[u, α];
y = Derivative[1, 0][f][u, α];
Assert[y*Hypergeometric2F1[1/3, 2/3, 1/2, -((1 + α)/(2 - α))^3*Sin[3*x]^2] == (2 - α)/(2*Sqrt[3])];
];
Clear[f];
If[Context[CarlsonRF] === "System`",
f[u_, v_] :=
Block[{x, y, z, b = Abs[α]^(3/2), ϕ},
ϕ = 2*I*Pi*{1, -1, 0};
ϕ = Which[
α > 0, (ArcCsch[b] + ϕ)/3,
α < 0, (ArcSech[b] + ϕ)/3
];
b = v*α;
{x, y, z} = Which[
b > 0, 1 + Csch[ϕ]^2/4,
b < 0, 1 - Sech[ϕ]^2/4,
True, Table[1, 3]
];
Assert[y == Conjugate[x] && z >= 3/4];
Assert[u == Re[CarlsonRF[x, y, z]]];
Assert[u == Sqrt[2]*CarlsonRF[Re[x] + Abs[x], z + Abs[x] + Abs[z - x], z + Abs[x] - Abs[z - x]]];
]
];
If[α != -1, With[{u = HypergeometricPFQ[{1/2, 1/2, 1}, {5/6, 7/6}, -α^3], s = sqrt[3*α]},
(* Differential equations:
u(α) + 3*α*u'(α) + (1 + α^3)/α*u''(α) + sqrt(3)/(4*α^(5/2)) = 0
5*α*u(α) + 23*α^2*u'(α) + 3*(1 + 5*α^3)*u''(a) + 2*α*(1 + α^3)*u'''(α) = 0
Wronskian(λ(α), λ'(α), u(α)) = π*i/(2*sqrt(3)*α^(3/2)*(1 + α^3)^2)
*)
Assert[DixonSM[s*u, α] == s];
Assert[DixonCM[s*u, α] == 1];
Assert[DixonF[s*u, α] == (3/2)*α + s*(3/5)*α^2*HypergeometricPFQ[{1/2, 1, 3/2}, {7/6, 11/6}, -α^3]];
Assert[DixonF[s*u + λ, α] == κ - (3/2)*α + s*(3/5)*α^2*HypergeometricPFQ[{1/2, 1, 3/2}, {7/6, 11/6}, -α^3]];
Assert[PossibleZeroQ[α] ||
DixonF[s*u - λ, α] == -κ + (3/2)*α + s^(-1)*HypergeometricPFQ[{-1/2, 1/2, 1}, {1/6, 5/6}, -α^3]];
If[α > 0, Assert[s*(u + (1/9)*α^(-3)*HypergeometricPFQ[{5/6, 7/6, 1}, {3/2, 3/2}, -α^(-3)]) == Sign[s]*(3/4)*λ]];
If[α > -1, With[{v = u*Sqrt[1 + α^3]},
With[{b = SetPrecision[(3*α)^3/4, Infinity]},
Assert[u == (1/2)*NIntegrate[(x^3 + b*(x - 1)^2)^(-1/2), {x, 1, Infinity}, WorkingPrecision -> pr]]; (* hard case at x = 3 *)
Assert[u == NIntegrate[(1 + b*x^2*(1 - x^2)^2)^(-1/2), {x, 0, 1}, WorkingPrecision -> pr]]; (* hard case at x = 3^(-1/2) *)
];
If[!PossibleZeroQ[α], With[{b = SetPrecision[1 + α^(-3), Infinity], ϕ = SetPrecision[ArcTan[α^(3/2)], Infinity]},
Assert[v == NIntegrate[(4*x^3 - (1 - 3*x)^2/b)^(-1/2), {x, 1, Infinity}, WorkingPrecision -> pr]];
Assert[v == NIntegrate[(1 - x^2*(3 - x^2)^2/(4*b))^(-1/2), {x, 0, 1}, WorkingPrecision -> pr]];
Assert[u == (2/3)*α^(-3/2)*N[NIntegrate[Hypergeometric2F1[1/3, 2/3, 1/2, b*Sin[x]^2], {x, 0, ϕ},
WorkingPrecision -> pr, PrecisionGoal -> pr/2
], pr/2]];
]];
f[v, -1];
]];
]];
If[α > -1, With[{u = HypergeometricPFQ[{1/2, 1/2, 1}, {5/6, 7/6}, α^3/(1 + α^3)]/Sqrt[1 + α^3], s = sqrt[α]},
(* Differential equations:
4*α^(3/2)*sqrt(1 + α^3)*(α*u(α) + 3*α^2*u'(α) + (1 + α^3)*u''(α)) + 1 = 0
α*(5 + 8*α^3)*u(α) + α^2*(23 + 32*α^3)*u'(α) + 3*(1 + α^3)*(1 + 6*α^3)*u''(α) + 2*α*(1 + α^3)^2*u'''(α) = 0
Wronskian(λ(α), λ'(α), u(α)) = π*i/(6*α^(3/2)*(1 + α^3)^(5/2))
*)
Assert[DixonSM[s*u, α] == s*E^((1/3)*ArcSinh[s^3])];
Assert[DixonCM[s*u, α] == E^((2/3)*ArcSinh[s^3])];
Assert[DixonSM[2*s*u, α] == 2*s*E^(-(1/3)*ArcSinh[s^3])];
Assert[DixonCM[2*s*u, α] == E^(-(2/3)*ArcSinh[s^3])];
Assert[Derivative[1, 0][DixonSM][s*u, α] == Sqrt[1 + α^3]*E^((1/3)*ArcSinh[s^3])];
Assert[Derivative[1, 0][DixonCM][s*u, α] == 0];
Assert[DixonF[s*u, α] == (1/2)*α + (s*(2/5)*α^2/Sqrt[1 + α^3])*HypergeometricPFQ[{1/2, 1/2, 1}, {7/6, 11/6}, α^3/(1 + α^3)]];
Assert[DixonF[s*u + λ, α] == κ - (1/2)*α + (s*(2/5)*α^2/Sqrt[1 + α^3])*HypergeometricPFQ[{1/2, 1/2, 1}, {7/6, 11/6}, α^3/(1 + α^3)]];
Assert[PossibleZeroQ[α] ||
DixonF[s*u - λ, α] == -κ + (3/2)*α + (Sqrt[1 + α^3]/(9*s))*(7 + 2*HypergeometricPFQ[{-1/2, -1/2, 1}, {1/6, 5/6}, α^3/(1 + α^3)])];
If[α < 0, Assert[s*(u - (1/9)*(Sqrt[1 + α^3]/α^3)*HypergeometricPFQ[{5/6, 7/6, 1}, {3/2, 3/2}, 1 + α^(-3)]) == Sign[s]*(1/4)*Abs[λ1 - λ2]]];
If[!PossibleZeroQ[α], With[{b = SetPrecision[α^(-3), Infinity], ϕ = SetPrecision[ArcSinh[α^(3/2)], Infinity]},
Assert[u == NIntegrate[(4*x^3 + (1 - 3*x)^2/b)^(-1/2), {x, 1, Infinity}, WorkingPrecision -> pr]];
Assert[u == NIntegrate[(1 + x^2*(3 - x^2)^2/(4*b))^(-1/2), {x, 0, 1}, WorkingPrecision -> pr]];
Assert[u == (2/3)*α^(-3/2)*N[NIntegrate[Hypergeometric2F1[1/3, 2/3, 1/2, b*Sinh[x]^2], {x, 0, ϕ},
WorkingPrecision -> pr, PrecisionGoal -> pr/2
], pr/2]];
]];
f[u, 1];
]];
If[α <= 0, With[{b = SetPrecision[α^3, Infinity]},
Assert[λ == NIntegrate[(1 - b*x^3)^(-1/3)*(1 + x^3)^(-2/3), {x, 0, Infinity}, WorkingPrecision -> pr]];
Assert[λ == NIntegrate[x*(x^3 - b)^(-1/3)*(1 + x^3)^(-2/3), {x, 0, Infinity}, WorkingPrecision -> pr]];
]];
f[u_, α_] := Through[{DixonSM, DixonCM}[u, α]]; (* ComapApply *)
If[!PossibleZeroQ[γ],
Assert[f[w2*u, w1*α] == {w2*s, c}];
Assert[makeλ[w1*α] == w2*{λ1, λ2, λ}];
Assert[makeκ[w1*α] == w1*{κ1, κ2, κ}];
Assert[f[w1*u, w2*α] == {w1*s, c}];
Assert[makeλ[w2*α] == w1*{λ2, λ, λ1}];
Assert[makeκ[w2*α] == w2*{κ2, κ, κ1}];
β = (2 - α)/(1 + α);
Assert[f[u*(1 + α)/(w2 - w1), β] == -{s + c - 1, s + w2*c - w1}/(s + w1*c - w2)];
Assert[DixonF[u*(1 + α)/(w2 - w1), β] ==
(-(1 - α + α^2)*u/3 + (w1 - w2)*(1 + α*w1)/3 + (w2*c - w1*s)*(s + w1*c + w2*α)/(s + w1*c - w2) + DixonF[u, α])*(w2 - w1)/(1 + α)];
Assert[DixonG[u*(1 + α)/(w2 - w1), β]/DixonG[0, β] ==
E^(-(1 - α + α^2)*u^2/6 + u*(w1 - w2)*(1 + α*w1)/3)*(s + w1*c - w2)*DixonG[u, α]/((w1 - w2)*DixonG[0, α])];
β = -α/γ;
Assert[f[3*u*γ/(w2 - w1), β] == -{3*γ*s*c, s^3 + w2*c^3 - w1}/(s^3 + w1*c^3 - w2)];
Assert[DixonG[3*u*γ/(w2 - w1), β]/DixonG[0, β] ==
E^(-3*α^2*u^2/2 + (w2 - 1)*α*u)*(s^3 + w1*c^3 - w2)*DixonG[u, α]^3/((w1 - w2)*DixonG[0, α]^3)];
Assert[f[u*γ/(w2 - w1), β]^3 == -{α*s - c + 1, α*s - w2*c + w1}/(α*s - w1*c + w2)];
Assert[DixonG[u*γ/(w2 - w1), β]^3/DixonG[0, β]^3 ==
E^(-α^2*u^2/2 + w2*α*u)*(α*s - w1*c + w2)*DixonG[u, α]/((w2 - w1)*DixonG[0, α])];
β = (α - 2)/b;
Assert[f[b*u, β] == {b*s*c, -(s - c^2 + c*s^2)}/(c - s^2 + s*c^2)];
Assert[DixonG[b*u, β]/DixonG[0, β] ==
E^((3*(1 - α)*u^2/2 - (1 + α)*u))*(c - s^2 + s*c^2)*DixonG[u, α]^3/DixonG[0, α]^3];
Assert[f[b*u/3, β]^3 == {α*s - c + 1, α*c - s + 1}/(s + c + α)];
Assert[DixonF[b*u/3, β] ==
((1 - α)*u - s + c - 1 + DixonF[u, α])/b];
Assert[DixonG[b*u/3, β]^3/DixonG[0, β]^3 ==
E^((1 - α)*u^2/2 - u)*(s + c + α)*DixonG[u, α]/((1 + α)*DixonG[0, α])];
β = (2*γ + α)/(γ - α);
Assert[f[(γ - α)*u, β] == {s^3 + (γ - α)*s*c, c^3 + (γ - α)*s*c}/(1 - (γ - α)*s*c)];
Assert[DixonF[(γ - α)*u, β] ==
(γ + 2*α)*u + 1 + (s^3 - c^3)/(1 - (γ - α)*s*c) + 3*DixonF[u, α]/(γ - α)];
Assert[DixonG[(γ - α)*u, β]/DixonG[0, β] ==
E^((γ - α)*((γ + 2*α)*u^2/2 + u))*(1 - (γ - α)*s*c)*DixonG[u, α]^3/DixonG[0, α]^3];
];
If[γ > 0,
β = (2 - α)/(1 + α);
Assert[makeλ[β] == (1 + α)/(w2 - w1)*{λ2 - λ1, -λ2, λ1}];
Assert[makeκ[β] == β + (3*{κ1 - κ2, κ2 - α, α - κ1} + (1 - α + α^2)*{λ2 - λ1, -λ2, λ1})/((w2 - w1)*(1 + α))];
β = -α/γ;
Assert[makeλ[β] == γ/(w2 - w1)*{λ2 - λ1, λ - λ2, λ1 - λ}];
Assert[makeκ[β] == β + ({κ1 - κ2, κ2 - κ, κ - κ1} + α^2*{λ2 - λ1, λ - λ2, λ1 - λ})/((w2 - w1)*γ)];
β = (α - 2)/b;
Assert[makeλ[β] == b*(λ/3 + {0, λ1, λ2})];
Assert[makeκ[β] == β + (κ - α + 3*{0, κ1 - α, κ2 - α} + (1 - α)*(λ + 3*{0, λ1, λ2}))/b];
β = (2*γ + α)/(γ - α);
Assert[makeλ[β] == (γ - α)*{λ, (λ1 - λ)/3, (λ2 - λ)/3}];
Assert[makeκ[β] == β + {3*(κ - α), κ1 - κ, κ2 - κ}/(γ - α) + (γ + 2*α)*{λ, (λ1 - λ)/3, (λ2 - λ)/3}];
];
Clear[f];
f = QPochhammer;
g[q_] := 1 - 24*QPolyGamma[1, 1, q]/Log[q]^2;
g[z_, q_] := (QPolyGamma[Log[q, z], q] + Log[1 - q])/Log[q];
If[NumericQ[λ1],
τ = I*Abs@Im[λ1/λ]/3;
q = -Sign[γ]*Exp[2*I*Pi*τ];
p = CubeRoot[q];
TimeConstrained[
Assert[α == -1 - 9*q*(f[q^9]/f[q])^3];
Assert[γ == -3*p*(f[q^3]/f[q])^4];
Assert[α/γ == 1 + (3*p)^(-1)*(f[p]/f[q^3])^3];
Assert[b == 3*f[q^3]^4/(f[q]^3*f[q^9])];
Assert[γ - α == f[p]^3*f[q^3]/f[q]^4];
Assert[λ/λd == f[q]^3/f[q^3]];
Assert[(κ - α)/λ - (3/4)*α^2 == (1/4)*f[q^3]^2*g[q]/f[q]^6];
Assert[(κ - α)/λ - (1/4)*α^2 == (3/4)*f[q^3]^2*g[q^3]/f[q]^6];
Assert[DixonSM[(3/2)*λ, α] == 2*f[q^2]^3*f[q^3]/(f[q]^3*f[q^6])];
If[q > 0,
Assert[DixonSM[(3/2)*I*Im[λ1], α] == -f[q^(1/2)]^3*f[q^3]/(f[q]^3*f[q^(3/2)])];
Assert[DixonSM[(3/2)*(λ + I*Im[λ1]), α] == -f[-q^(1/2)]^3*f[q^3]/(f[q]^3*f[-q^(3/2)])],
Assert[DixonSM[(3/2)*λ1, α] == -f[q^2]^3*f[q^3]/(f[q]^3*f[q^6]) + Sqrt[-3*q]*(w1 - w2)*f[q^6]^3/(f[q]*f[q^2]*f[q^3])]
];
Assert[DixonSM[(2/3)*I*Im[λ1], α] == -w2*f[p]/f[w2*p]];
Assert[DixonCM[(2/3)*I*Im[λ1], α] == -w1*f[w1*p]/f[w2*p]];
With[{z = w1^(u/(2*λ))},
Assert[DixonSM[u, α] ==
-(Im[z]/Im[w2*z])*Abs[f[q*z^2, q]/f[q*w1*z^2, q]]^2];
Assert[DixonCM[u, α] ==
-(Im[w1*z]/Im[w2*z])*Abs[f[q*w2*z^2, q]/f[q*w1*z^2, q]]^2];
Assert[(α - κ)*u/λ - α/2 + DixonF[u, α] ==
λd*Im[w1]*(4*Im[g[q*w1*z^2, q]] + Re[w2*z]/Im[w2*z])/λ];
Assert[Exp[(α - κ)*u^2/(2*λ) - α*u/2]*DixonG[u, α]/DixonG[0, α] ==
(Im[w2*z]/Im[w2])*Abs[f[q*w1*z^2, q]/f[q*w1, q]]^2];
],
0.1
];
With[{n = 1 + Boole[γ > 0]},
q = Exp[-2*I*Pi/(3*n*τ)];
p = (-1)^(n - 1)*q^(3/n);
TimeConstrained[
With[{q = Surd[p, 9]},
Assert[α == -1 - (3*q)^(-1)*(f[q]/f[p])^3];
Assert[γ == -(3*q)^(-1)*(f[q^3]/f[p])^4];
Assert[α/γ == 1 + 9*q^3*(f[p^3]/f[q^3])^3];
Assert[b == Abs[q]^(-2/3)*f[q^3]^4/(f[p]^3*f[q])];
Assert[γ - α == 3*q^2*f[p^3]^3*f[q^3]/f[p]^4];
Assert[Abs@Im[λ1] == (2*Pi/n)*Abs[q]*f[p]^3/f[q^3]];
Assert[Im[κ1]/Im[λ1] - (3/4)*α^2 == -(1/12)*q^(-2)*f[q^3]^2*g[p]/f[p]^6];
Assert[Im[κ1]/Im[λ1] - (1/4)*α^2 == -(1/36)*q^(-2)*f[q^3]^2*g[q^3]/f[p]^6];
(* ... *)
];
With[{z = q^(u/λ)},
Assert[DixonSM[u, α] ==
q^(1/(3*n))*z^(-n/3)*f[z, p]*f[p/z, p]/(f[q*z, p]*f[p/(q*z), p])];
Assert[DixonCM[u, α] ==
z^(1 - 2*n/3)*f[q/z, p]*f[p*z/q, p]/(f[q*z, p]*f[p/(q*z), p])];
Assert[Exp[-Im[κ1]*u^2/(2*Im[λ1]) - α*u/2]*DixonG[u, α]/DixonG[0, α] ==
z^(n/3 - 1/2)*f[q*z, p]*f[p/(q*z), p]/(f[q, p]*f[p/q, p])];
If[n == 1, Assert[-Im[κ1]*u/Im[λ1] - α/2 + DixonF[u, α] ==
-((1/6)*Log[p] + QPolyGamma[1/3 + u/(3*λ), p] - QPolyGamma[2/3 - u/(3*λ), p])/(3*λ)]];
],
0.1
];
];
];
Block[{u, v},
{u, v} = RandomReal[{0, λ0}, 2, WorkingPrecision -> pr];
With[{s1 = DixonSM[u], c1 = DixonCM[u], s2 = DixonSM[v], c2 = DixonCM[v]},
Assert[DixonSM[u + w1*v] == (s1*c2*(c2 + s1*c1*s2^2) + w1*s2*c1*(c1 + s2*c2*s1^2))/(1 - s1^3*s2^3)];
Assert[DixonCM[u + w1*v] == (c1*(c2 + s1*c1*s2^2) + w1*s1*s2*(s2*c1^2 - s1*c2^2))/(1 - s1^3*s2^3)];
Assert[s1*AppellF1[1/3, 2/3, 1, 4/3, s1^3, s1^3*s2^3] ==
u + s2*(u*DixonF[v] - (Log[DixonG[u + v]] + w1*Log[DixonG[u + w1*v]] + w2*Log[DixonG[u + w2*v]])/3)/c2^2];
Assert[s1^2*AppellF1[2/3, 1/3, 1, 5/3, s1^3, s1^3*s2^3] ==
(2/3)*(Log[DixonG[u + v]] + w2*Log[DixonG[u + w1*v]] + w1*Log[DixonG[u + w2*v]])/(s2*c2)];
Assert[c2*s1*AppellF1[1/3, 2/3, 1, 4/3, s1^3, s1^3*s2^3] - s2^2*s1^2*AppellF1[2/3, 1/3, 1, 5/3, s1^3, s1^3*s2^3]/2 ==
c2*u - (s2/c2)*(Log[DixonG[u + v]] - Log[DixonG[u]] - Log[DixonG[v]] + Log[DixonG[0]] - Log[1 - s1^3*s2^3]/3 - u*DixonF[v])];
];
];
Clear[f];
f[α_, x_, y_] := With[{λ = N[makeλ[α][[1]], pr]},
If[α < 2, Assert[λ == (2*Pi/(Sqrt[3]*(2 - α)))*Hypergeometric2F1[1/3, 2/3, 1, -((1 + α)/(2 - α))^3]]];
With[{b = α^3},
(Sqrt[α]*{Sqrt[3]*x, y/Sqrt[1 + b]} . HypergeometricPFQ[{1/2, 1/2, 1}, {5/6, 7/6}, {-b, b/(1 + b)}]) / λ
]
];
With[{ϕ = GoldenRatio},
Assert[f[1/3, 1, 0] == 1/2];
Assert[f[1/(3*ϕ^2), 1, 0] == 1/3];
Assert[f[ϕ^2/3, 1, 0] == 2/3];
Assert[f[Root[#^3 + 2*# - 1 &, 1]^2/3, 1, 0] == 1/4];
Assert[f[1/2, 0, 1] == 1/3];
Assert[f[80^(-1/3), 0, 1] == 1/4];
Assert[f[1/(5 + Sqrt[33]), 0, 1] == 1/6];
Assert[f[3^(-1/3), 1, 1] == 1];
Assert[f[5^(-2/3)/3, 1, 1] == 1/2];
Assert[f[48^(-1/3), 1, 2] == 1];
Assert[f[(2/25)^(1/3), 1, 3] == 3/2];
Assert[f[(7 - Sqrt[33])/8, 1, 3] == 1];
Assert[f[(7 + Sqrt[33])/8, 1, 3] == 2];
Assert[f[Root[4*#^3 - 10*#^2 + 12*# - 3 &, 1]^2/3, 1, 3] == 1/2];
Assert[f[Root[2*#^3 - 4*#^2 + 3*# - 3 &, 1]^2/3, -1, 3] == 1/2];
];
] /.
HoldPattern[_[vars:{___Symbol}, expr_]] :>
Block[{$AssertFunction = Automatic},
With[{temp = Unique[Unevaluated[vars], Temporary]},
WithCleanup[
Unevaluated[expr] /. Dispatch[Thread@Hold[vars, temp] /. Hold[lhs_, rhs_] :> (HoldPattern[lhs] -> rhs)],
Remove @@ Unevaluated[temp]
] /; ListQ[temp]
]
];Theta functions:
Module[{u, q, a, b},
u = RandomReal[WorkingPrecision -> 50];
q = RandomReal[WorkingPrecision -> 50];
a = ((EllipticTheta[4, q]/EllipticTheta[4, Pi/3, q])^2 + 2*EllipticTheta[4, Pi/3, q]/EllipticTheta[4, q])/3;
b = EllipticTheta[2, q]*EllipticTheta[3, q]*EllipticTheta[4, q]/(2*Sqrt[3]*EllipticTheta[1, Pi/3, q]);
Assert[EllipticTheta[1, u, q^3]^3 == b*(a*EllipticTheta[1, u, q] - EllipticTheta[1, u + Pi/3, q] + EllipticTheta[1, u + 2*Pi/3, q])];
Assert[EllipticTheta[2, u, q^3]^3 == b*(a*EllipticTheta[2, u, q] - EllipticTheta[2, u + Pi/3, q] + EllipticTheta[2, u + 2*Pi/3, q])];
Assert[EllipticTheta[3, u, q^3]^3 == b*(a*EllipticTheta[3, u, q] + EllipticTheta[3, u + Pi/3, q] + EllipticTheta[3, u + 2*Pi/3, q])];
Assert[EllipticTheta[4, u, q^3]^3 == b*(a*EllipticTheta[4, u, q] + EllipticTheta[4, u + Pi/3, q] + EllipticTheta[4, u + 2*Pi/3, q])];
Assert[EllipticTheta[1, u, q^3] == EllipticTheta[1, u/3, q]*EllipticTheta[1, (u + Pi)/3, q]*EllipticTheta[1, (u + 2*Pi)/3, q]/(3*b)];
Assert[EllipticTheta[2, u, q^3] == -EllipticTheta[2, u/3, q]*EllipticTheta[2, (u + Pi)/3, q]*EllipticTheta[2, (u + 2*Pi)/3, q]/(3*b)];
Assert[EllipticTheta[3, u, q^3] == EllipticTheta[3, u/3, q]*EllipticTheta[3, (u + Pi)/3, q]*EllipticTheta[3, (u + 2*Pi)/3, q]/(3*b)];
Assert[EllipticTheta[4, u, q^3] == EllipticTheta[4, u/3, q]*EllipticTheta[4, (u + Pi)/3, q]*EllipticTheta[4, (u + 2*Pi)/3, q]/(3*b)];
Assert[EllipticTheta[1, u, q^9] == -(EllipticTheta[1, u/3, q] - EllipticTheta[1, (u + Pi)/3, q] + EllipticTheta[1, (u + 2*Pi)/3, q])/3];
Assert[EllipticTheta[2, u, q^9] == (EllipticTheta[2, u/3, q] - EllipticTheta[2, (u + Pi)/3, q] + EllipticTheta[2, (u + 2*Pi)/3, q])/3];
Assert[EllipticTheta[3, u, q^9] == (EllipticTheta[3, u/3, q] + EllipticTheta[3, (u + Pi)/3, q] + EllipticTheta[3, (u + 2*Pi)/3, q])/3];
Assert[EllipticTheta[4, u, q^9] == (EllipticTheta[4, u/3, q] + EllipticTheta[4, (u + Pi)/3, q] + EllipticTheta[4, (u + 2*Pi)/3, q])/3];
With[{s = EllipticTheta[1, u, q], m = EllipticTheta[1, u + Pi/3, q], c = EllipticTheta[1, u + 2*Pi/3, q]},
Assert[EllipticTheta[1, u+Pi/9, q]^3 == ( EllipticTheta[1, Pi/9, q]^3*c*m^2 + EllipticTheta[1, 2*Pi/9, q]^3*s*c^2 + EllipticTheta[1, 4*Pi/9, q]^3*m*s^2)/EllipticTheta[1, Pi/3, q]^3]
Assert[EllipticTheta[1, u-Pi/9, q]^3 == (-EllipticTheta[1, Pi/9, q]^3*m*c^2 + EllipticTheta[1, 2*Pi/9, q]^3*s*m^2 - EllipticTheta[1, 4*Pi/9, q]^3*c*s^2)/EllipticTheta[1, Pi/3, q]^3];
];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Inverse of Dixon's sm function:
(* Numerical computation *)
InverseDixonSM[z_, alf_] /;
VectorQ[{z, alf}, NumericQ] && Precision[{z, alf}] < Infinity :=
Module[{g2, g3, c, t},
{g2, g3} = {alf/12 (alf^3 - 8), (8 - alf^3 (20 + alf^3))/216};
t = 1 + alf*z;
c = Nearest[c /. {ToRules[NRoots[c^3 + z^3 == 1 + 3*alf*c*z, c]]}, t][[1]];
t = (1 + alf^3)/(3*(c - t));
t = InverseWeierstrassP[{-z*t - alf^2/4, (1 + c)*t}, {g2, g3}];
t /; InexactNumberQ[t]
];
(* Series expansion *)
InverseDixonSM[z, α] == z*(
+ Sum[Pochhammer[n + 2/3, n]/((3*n+1)*n!)*(α*z)^(3*n)*Hypergeometric2F1[n + 1/3, 2*(n + 1/3), n + 4/3, z^3], {n, 0, Infinity}]
- Sum[Pochhammer[n + 4/3, n]/((3*n+2)*n!)*(α*z)^(3*n+1)*Hypergeometric2F1[n + 2/3, 2*(n + 2/3), n + 5/3, z^3], {n, 0, Infinity}]
)
(* Branch points *)
BranchPoints[InverseDixonSM[z, α], z] == {
E^((2/3)*ArcSinh[α^(3/2)]),
E^(-(2/3)*ArcSinh[α^(3/2)]),
E^((2/3)*(ArcSinh[α^(3/2)]+I*Pi)),
E^(-(2/3)*(ArcSinh[α^(3/2)]+I*Pi)),
E^((2/3)*(ArcSinh[α^(3/2)]-I*Pi)),
E^(-(2/3)*(ArcSinh[α^(3/2)]-I*Pi)),
ComplexInfinity
}(* Addition *)
D[f[x] + g[x], x] == f'[x] + g'[x]
D[f[x] + g[x], {x, 2}] == f''[x] + g''[x]
(* Subtraction *)
D[f[x] - g[x], x] == f'[x] - g'[x]
D[f[x] - g[x], {x, 2}] == f''[x] - g''[x]
(* Multiplication by constant *)
D[a * f[x], x] == a * f'[x]
D[a * f[x], {x, 2}] == a * f''[x]
(* Dot product *)
D[f[x] . g[x], x] == f[x] . g'[x] + f'[x] . g[x]
D[f[x] . g[x], {x, 2}] == f[x] . g''[x] + 2 * f'[x] . g'[x] + f''[x] . g[x]
(* Hadamard product *)
D[f[x] * g[x], x] == f[x] * g'[x] + f'[x] * g[x]
D[f[x] * g[x], {x, 2}] == f[x] * g''[x] + 2 * f'[x] * g'[x] + f''[x] * g[x]
(* Inverse *)
D[Inverse[f[x]], x] == -Inverse[f[x]] . f'[x] . Inverse[f[x]]
D[Inverse[f[x]], {x, 2}] == Inverse[f[x]] . (2 * f'[x] . Inverse[f[x]] . f'[x] - f''[x]) . Inverse[f[x]]
(* Transpose *)
D[Transpose[f[x]], x] == Transpose[f'[x]]
D[Transpose[f[x]], {x, 2}] == Transpose[f''[x]]
(* Trace *)
D[Tr[f[x]], x] == Tr[f'[x]]
D[Tr[f[x]], {x, 2}] == Tr[f''[x]]
D[Tr[MatrixPower[f[x], n]], x] == n * Tr[MatrixPower[f[x], n - 1] . f'[x]]
D[Tr[MatrixExp[f[x]]], x] == Tr[MatrixExp[f[x]] . f'[x]]
D[Tr[MatrixLog[f[x]]], x] == Tr[Inverse[f[x]] . f'[x]]
D[Tr[MatrixFunction[p, f[x]]], x] == Tr[MatrixFunction[p', f[x]] . f'[x]]
(* Determinant *)
D[Det[f[x]], x] == Det[f[x]] * Tr[Inverse[f[x]] . f'[x]]
D[Det[f[x]], {x, 2}] == Det[f[x]] * (Tr[Inverse[f[x]] . f'[x]]^2 + Tr[Inverse[f[x]] . (f''[x] - f'[x] . Inverse[f[x]] . f'[x])])
D[Det[MatrixPower[f[x], n]], x] == n * Det[MatrixPower[f[x], n]] * Tr[Inverse[f[x]] . f'[x]]
D[Det[MatrixExp[f[x]]], x] == Det[MatrixExp[f[x]]] * Tr[f'[x]]
D[Det[MatrixLog[f[x]]], x] == Det[MatrixLog[f[x]]] * Tr[f'[x] . Inverse[f[x]] . Inverse[MatrixLog[f[x]]]]
D[Det[MatrixFunction[p, f[x]]], x] == Det[MatrixFunction[p, f[x]]] * Tr[f'[x] . MatrixFunction[p', f[x]] . Inverse[MatrixFunction[p, f[x]]]]
(* Logarithm of determinant *)
D[Log[Det[f[x]]], x] == Tr[Inverse[f[x]] . f'[x]]
D[Log[Det[f[x]]], {x, 2}] == Tr[Inverse[f[x]] . (f''[x] - f'[x] . Inverse[f[x]] . f'[x])]
D[Log[Det[MatrixPower[f[x], n]]], x] == n * Tr[Inverse[f[x]] . f'[x]]
D[Log[Det[MatrixExp[f[x]]]], x] == Tr[f'[x]]
D[Log[Det[MatrixLog[f[x]]]], x] == Tr[f'[x] . Inverse[f[x]] . Inverse[MatrixLog[f[x]]]]
D[Log[Det[MatrixFunction[p, f[x]]]], x] == Tr[f'[x] . MatrixFunction[p', f[x]] . Inverse[MatrixFunction[p, f[x]]]]{
my(
f=default(format),
p=default(realprecision),
bp=default(realbitprecision),
w,w1,w2,w3,g2,g3,z1,z2,z3,z4,p1,p2,p3,p4,P
);
default(realbitprecision,bp+64);
default(format,f);
setrand(getwalltime());
f=strprintf("%%.%dg\n",p);
w1=tan(random(1.)*Pi/2);
w2=I*tan(random(1.)*Pi/2);
w3=-w1-w2;
w=2*[w1,w2];
g2=elleisnum(w,4)/12;
g3=-elleisnum(w,6)/216;
P=((z,n)->my(a,b,x,y,t);[x,y]=ellwp(w,z,1);t=vector(n);a='x;t[1]=x;for(i=2,n,t[i]=if(i%2,a=(4*'x^3-g2*'x-g3)*deriv(b,'x)+(6*'x^2-g2/2)*b;subst(a,'x,x),b=deriv(a,'x);y*subst(b,'x,x)));return(t));
z1=random(1.)*w1;p1=P(z1,6);
z2=random(1.)*w1;p2=P(z2,4);
z3=random(1.)*w1;p3=P(z3,4);
z4=random(1.)*w1;p4=P(z4,4);
printf(f, matdet([
1,p1[1];
1,p2[1]
]));
printf(f, ellsigma(w,z1-z2)*ellsigma(w,z1+z2)/(ellsigma(w,z1)*ellsigma(w,z2))^2);
printf(f, matdet([
1,p1[1],p1[2];
1,p2[1],p2[2];
1,p3[1],p3[2]
]));
printf(f, -2!*ellsigma(w,z1-z2)*ellsigma(w,z1-z3)*ellsigma(w,z2-z3)*ellsigma(w,z1+z2+z3)/(ellsigma(w,z1)*ellsigma(w,z2)*ellsigma(w,z3))^3);
printf(f, matdet([
1,p1[1],p1[2],p1[3];
1,p2[1],p2[2],p2[3];
1,p3[1],p3[2],p3[3];
1,p4[1],p4[2],p4[3]
]));
printf(f, -2!*3!*ellsigma(w,z1-z2)*ellsigma(w,z1-z3)*ellsigma(w,z1-z4)*ellsigma(w,z2-z3)*ellsigma(w,z2-z4)*ellsigma(w,z3-z4)*ellsigma(w,z1+z2+z3+z4)/(ellsigma(w,z1)*ellsigma(w,z2)*ellsigma(w,z3)*ellsigma(w,z4))^4);
printf(f, matdet([
1,p1[1];
0,p1[2]
]));
printf(f, -ellsigma(w,2*z1)/ellsigma(w,z1)^4);
printf(f, matdet([
1,p1[1],p1[2];
0,p1[2],p1[3];
0,p1[3],p1[4]
]));
printf(f, 2!^2*ellsigma(w,3*z1)/ellsigma(w,z1)^9);
printf(f, matdet([
1,p1[1],p1[2],p1[3];
0,p1[2],p1[3],p1[4];
0,p1[3],p1[4],p1[5];
0,p1[4],p1[5],p1[6]
]));
printf(f, -(2!*3!)^2*ellsigma(w,4*z1)/ellsigma(w,z1)^16);
printf(f, matdet([
1,p1[1],p1[2];
1,p2[1],p2[2];
1,p3[1],p3[2]
])/matdet([
1,p1[1];
1,p2[1]
]));
printf(f, -2!*ellsigma(w,z1-z3)*ellsigma(w,z2-z3)*ellsigma(w,z1+z2+z3)/(ellsigma(w,z1)*ellsigma(w,z2)*ellsigma(w,z1+z2)*ellsigma(w,z3)^3));
printf(f, matdet([
1,p1[1],p1[2],p1[3];
1,p2[1],p2[2],p2[3];
1,p3[1],p3[2],p3[3];
1,p4[1],p4[2],p4[3]
])/matdet([
1,p1[1],p1[2];
1,p2[1],p2[2];
1,p3[1],p3[2]
]));
printf(f, 3!*ellsigma(w,z1-z4)*ellsigma(w,z2-z4)*ellsigma(w,z3-z4)*ellsigma(w,z1+z2+z3+z4)/(ellsigma(w,z1)*ellsigma(w,z2)*ellsigma(w,z3)*ellsigma(w,z1+z2+z3)*ellsigma(w,z4)^4));
printf(f, matdet([
1,p1[1],p1[2];
0,p1[2],p1[3];
1,p2[1],p2[2]
])/matdet([
1,p1[1];
0,p1[2]
]));
printf(f, -2!*ellsigma(w,z1-z2)^2*ellsigma(w,2*z1+z2)/(ellsigma(w,z1)^2*ellsigma(w,2*z1)*ellsigma(w,z2)^3));
printf(f, matdet([
1,p1[1],p1[2],p1[3];
0,p1[2],p1[3],p1[4];
0,p1[3],p1[4],p1[5];
1,p2[1],p2[2],p2[3]
])/matdet([
1,p1[1],p1[2];
0,p1[2],p1[3];
0,p1[3],p1[4]
]));
printf(f, 3!*ellsigma(w,z1-z2)^3*ellsigma(w,3*z1+z2)/(ellsigma(w,z1)^3*ellsigma(w,3*z1)*ellsigma(w,z2)^4));
printf(f, (1/2)*matdet([
1,p1[2];
1,p2[2]
])/matdet([
1,p1[1];
1,p2[1]
]));
printf(f, ellzeta(w,z1+z2)-ellzeta(w,z1)-ellzeta(w,z2));
printf(f, (1/3)*matdet([
1,p1[1],p1[3];
1,p2[1],p2[3];
1,p3[1],p3[3]
])/matdet([
1,p1[1],p1[2];
1,p2[1],p2[2];
1,p3[1],p3[2]
]));
printf(f, ellzeta(w,z1+z2+z3)-ellzeta(w,z1)-ellzeta(w,z2)-ellzeta(w,z3));
printf(f, (1/4)*matdet([
1,p1[1],p1[2],p1[4];
1,p2[1],p2[2],p2[4];
1,p3[1],p3[2],p3[4];
1,p4[1],p4[2],p4[4]
])/matdet([
1,p1[1],p1[2],p1[3];
1,p2[1],p2[2],p2[3];
1,p3[1],p3[2],p3[3];
1,p4[1],p4[2],p4[3]
]));
printf(f, ellzeta(w,z1+z2+z3+z4)-ellzeta(w,z1)-ellzeta(w,z2)-ellzeta(w,z3)-ellzeta(w,z4));
default(realbitprecision,bp);
}BeginPackage["CubicTheta`"];
Unprotect[CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime];
ClearAll[CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime];
Begin["`Private`"];
CubicTheta[{n:(-1 | 0 | 1), m:(-1 | 0 | 1)}, {u_, v_}, q_] /;
(VectorQ[{u, v, q}, NumericQ] && Abs[q] < 1 && Precision[{u, v, q}] < Infinity) :=
Module[{p = q^3, t, a = (u + v)*Pi, b = (u - v + m/3)*Pi},
If[n != 0, a += n*I*Log[q]];
t = EllipticTheta[3, a, p]*EllipticTheta[3, b, q];
If[!PossibleZeroQ[p], t += EllipticTheta[2, a, p]*EllipticTheta[2, b, q]*q/(p^(1/4)*q^(1/4))];
If[n != 0, t *= q^(-n^2/3)*E^(-2*n*I*a/3)];
Return[t];
];
CubicThetaA[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
(QPochhammer[q]^3 + 9*q*QPochhammer[q^9]^3) / QPochhammer[q^3];
CubicThetaA /: HoldPattern[Derivative[n_][CubicThetaA]] /; n > 0 := Derivative[n - 1][CubicThetaAPrime];
CubicThetaB[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
QPochhammer[q]^3 / QPochhammer[q^3];
CubicThetaB /: HoldPattern[Derivative[n_][CubicThetaB]] /; n > 0 := Derivative[n - 1][CubicThetaBPrime];
CubicThetaC[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
3 * q^(1/3) * QPochhammer[q^3]^3 / QPochhammer[q];
CubicThetaC /: HoldPattern[Derivative[n_][CubicThetaC]] /; n > 0 := Derivative[n - 1][CubicThetaCPrime];
CubicThetaAPrime[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
With[{a = CubicThetaA[q], c = CubicThetaC[q]},
(c^3/3 - a^3/12 + a*WeierstrassZeta[1, WeierstrassInvariants[{1, -I*Log[q]/(2*Pi)}]]/Pi^2)/q
];
CubicThetaAPrime /: HoldPattern[Derivative[1][CubicThetaAPrime]] :=
(With[{a = CubicThetaA[#], b = CubicThetaB[#], c = CubicThetaC[#], ap = CubicThetaAPrime[#]},
2*b^3*c^3/(a*(3*#)^2) - ap/# + 2*ap^2/a
] &);
CubicThetaBPrime[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
With[{a = CubicThetaA[q], b = CubicThetaB[q]},
b*(-a^2/12 + WeierstrassZeta[1, WeierstrassInvariants[{1, -I*Log[q]/(2*Pi)}]]/Pi^2)/q
];
CubicThetaBPrime /: HoldPattern[Derivative[1][CubicThetaBPrime]] :=
(With[{a = CubicThetaA[#], b = CubicThetaB[#], c = CubicThetaC[#], bp = CubicThetaBPrime[#]},
-a*b*c^3/(3*#)^2 - bp/# + 2*bp^2/b
] &);
CubicThetaCPrime[q_] /;
(NumericQ[q] && Abs[q] < 1 && Precision[q] < Infinity) :=
With[{a = CubicThetaA[q], c = CubicThetaC[q]},
c*(a^2/4 + WeierstrassZeta[1, WeierstrassInvariants[{1, -I*Log[q]/(2*Pi)}]]/Pi^2)/q
];
CubicThetaCPrime /: HoldPattern[Derivative[1][CubicThetaCPrime]] :=
(With[{a = CubicThetaA[#], b = CubicThetaB[#], c = CubicThetaC[#], cp = CubicThetaCPrime[#]},
-a*b^3*c/(3*#)^2 - cp/# + 2*cp^2/c
] &);
End[];
SetAttributes[{CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime}, {Listable, NumericFunction, ReadProtected}];
SetAttributes[{CubicTheta}, {NHoldFirst, NumericFunction, ReadProtected}];
Protect[CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime];
EndPackage[];Test cases for 0 < q < 1:
Module[{a, b, c, q},
q = RandomReal[WorkingPrecision -> 20];
a = CubicThetaA[q];
b = CubicThetaB[q];
c = CubicThetaC[q];
Assert[TrueQ[And[
CubicThetaA[q^2] == a*Cos[(2/3)*ArcTan[c^(3/2)/b^(3/2)]],
CubicThetaB[q^2] == a^(1/2)*b^(1/2)*Cos[(1/3)*ArcCos[b^(3/2)/a^(3/2)]],
CubicThetaC[q^2] == a^(1/2)*c^(1/2)*Sin[(1/3)*ArcSin[c^(3/2)/a^(3/2)]],
CubicThetaA[q^3] == (1/3)*(a+2*b),
CubicThetaB[q^3] == ((1/3)*b*(a^2+a*b+b^2))^(1/3),
CubicThetaC[q^3] == (1/3)*(a-b),
CubicThetaA[q^4] == (1/4)*(a+6^(1/2)*b^(3/4)*c^(3/4)*a^(-1/2)*Sin[(2/3)*ArcTan[c^(3/2)/b^(3/2)]]^(-1/2)),
CubicThetaB[q^4] == (1/4)*(b+3^(1/2)*b^(1/4)*c^(3/4)*Tan[(1/3)*ArcTan[c^(3/2)/b^(3/2)]]^(-1/2)),
CubicThetaC[q^4] == (1/4)*(c-3^(1/2)*b^(3/4)*c^(1/4)*Tan[(1/3)*ArcTan[c^(3/2)/b^(3/2)]]^(1/2)),
CubicThetaA[q^(1/2)] == 2*a*Cos[(2/3)*ArcTan[b^(3/2)/c^(3/2)]],
CubicThetaB[q^(1/2)] == 2*a^(1/2)*b^(1/2)*Sin[(1/3)*ArcSin[b^(3/2)/a^(3/2)]],
CubicThetaC[q^(1/2)] == 2*a^(1/2)*c^(1/2)*Cos[(1/3)*ArcCos[c^(3/2)/a^(3/2)]],
CubicThetaA[q^(1/3)] == a+2*c,
CubicThetaB[q^(1/3)] == a-c,
CubicThetaC[q^(1/3)] == (9*c*(a^2+a*c+c^2))^(1/3),
CubicThetaA[q^(1/4)] == (a+6^(1/2)*b^(3/4)*c^(3/4)*a^(-1/2)*Sin[(2/3)*ArcTan[b^(3/2)/c^(3/2)]]^(-1/2)),
CubicThetaB[q^(1/4)] == (b-3^(1/2)*b^(1/4)*c^(3/4)*Tan[(1/3)*ArcTan[b^(3/2)/c^(3/2)]]^(1/2)),
CubicThetaC[q^(1/4)] == (c+3^(1/2)*b^(3/4)*c^(1/4)*Tan[(1/3)*ArcTan[b^(3/2)/c^(3/2)]]^(-1/2))
]]];
(* CubicTheta[{0, 0}, {0, 0}, q]^3 == CubicThetaA[q] *)
(* CubicTheta[{0, 1}, {0, 0}, q]^3 == CubicThetaB[q] *)
(* CubicTheta[{1, 0}, {0, 0}, q]^3 == CubicThetaC[q] *)
(* CubicTheta[{1, 1}, {0, 0}, q]^3 == 0 *)
(* CubicTheta[{0, 0}, {x, y}, q]^3 + CubicTheta[{1, 1}, {x, y}, q]^3 == CubicTheta[{0, 1}, {x, y}, q]^3 + CubicTheta[{1, 0}, {x, y}, q]^3 *)
] & // Block[{$AssertFunction = Automatic}, #[]] &;Level 5 and 13 analogues:
Module[{q, f, α, β},
f[a_List, q_] := Times @@ (QPochhammer[#, q] & /@ a);
q = RandomReal[WorkingPrecision -> 50];
(* Level 5 *)
α = GoldenRatio;
β = 1/α;
Assert[(f[q^{2,3},q^5]^5 - α^5*q*f[q^{1,4},q^5]^5)*f[{q^5},q^5]^8 == f[(-1)^((2/5)*{1,4})*q,q]^5*f[{q},q]^8];
Assert[(f[q^{2,3},q^5]^5 + β^5*q*f[q^{1,4},q^5]^5)*f[{q^5},q^5]^8 == f[(-1)^((2/5)*{2,3})*q,q]^5*f[{q},q]^8];
Assert[f[q^{10,15,25},q^25] - α*q*f[q^{5,20,25},q^25] == f[(-1)^((2/5)*{1,4,0})*q,q]];
Assert[f[q^{10,15,25},q^25] + β*q*f[q^{5,20,25},q^25] == f[(-1)^((2/5)*{2,3,0})*q,q]];
(* Level 13 *)
α = (3 + Sqrt[13])/2;
β = 1/α;
Assert[f[q^{2,5,6,7,8,11,13,13},q^13] - α*q*f[q^{1,3,4,9,10,12,13,13},q^13] == f[(-1)^((2/13)*{1,3,4,9,10,12,0,0})*q,q]];
Assert[f[q^{2,5,6,7,8,11,13,13},q^13] + β*q*f[q^{1,3,4,9,10,12,13,13},q^13] == f[(-1)^((2/13)*{2,5,6,7, 8,11,0,0})*q,q]];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Fourier series (with q = exp(-π K(k')/K(k))):
Specific values for equianharmonic (α=0) case: