You are not logged in.
Duplication theorems:
c[x_, y_] := CarlsonRC[x, y];
d[x_, y_, z_] := CarlsonRD[x, y, z];
f[x_, y_, z_] := CarlsonRF[x, y, z];
g[x_, y_, z_] := CarlsonRG[x, y, z] - (1/6)*(x + y + z)*CarlsonRF[x, y, z];
j[x_, y_, z_, p_] := CarlsonRJ[x, y, z, p];
i[x_, y_, z_] := With[{a = (x + y + z)/3},
WeierstrassSigma[CarlsonRF[x, y, z], {(2/3)*((x - y)^2 + (x - z)^2 + (y - z)^2), 4*(a - x)*(a - y)*(a - z)}]
];
{x, y, z, p} = RandomReal[{1, 2}, 4, WorkingPrecision -> 50];
λ = Sqrt[x]*Sqrt[y] + Sqrt[x]*Sqrt[z] + Sqrt[y]*Sqrt[z];
On[Assert];
Assert[c[x, y] == c[(Sqrt[x] + Sqrt[y])^2/4, (y + Sqrt[x]*Sqrt[y])/2]];
Assert[d[x, y, z] == (1/4)*d[(x + λ)/4, (y + λ)/4, (z + λ)/4] + 3/(Sqrt[z]*(z + λ))];
Assert[f[x, y, z] == f[(x + λ)/4, (y + λ)/4, (z + λ)/4]];
Assert[g[x, y, z] == 4*g[(x + λ)/4, (y + λ)/4, (z + λ)/4] - (1/2)*(Sqrt[x] + Sqrt[y] + Sqrt[z])];
Assert[i[x, y, z] == (1/8)*i[(x + λ)/4, (y + λ)/4, (z + λ)/4]^4*(Sqrt[x] + Sqrt[y])*(Sqrt[x] + Sqrt[z])*(Sqrt[y] + Sqrt[z])];
α = p*(Sqrt[x] + Sqrt[y] + Sqrt[z]) + Sqrt[x]*Sqrt[y]*Sqrt[z];
β = Sqrt[p]*(p + λ);
Assert[j[x, y, z, p] == (1/4)*j[(x + λ)/4, (y + λ)/4, (z + λ)/4, (p + λ)/4] + 3*CarlsonRC[α^2, β^2]];
δ = (p - x)*(p - y)*(p - z); (* β^2 - α^2 *)
a = Sqrt@Abs[δ]; (* sqrt(|β^2 - α^2|) *)
b = (Sqrt[p] + Sqrt[x])*(Sqrt[p] + Sqrt[y])*(Sqrt[p] + Sqrt[z]); (* β + α *)
(* Circular case *)
If[δ > 0, Assert[CarlsonRC[α^2, β^2] == 2*ArcTan[b, a]/a]];
(* Hyperbolic case *)
If[δ < 0, Assert[CarlsonRC[α^2, β^2] == Log[1 + a*(1 + a/b)/β]/a]];
(* Degenerate case *)
If[δ == 0, Assert[CarlsonRC[α^2, β^2] == 1/α == 1/β]];Addition theorems:
x = RandomReal[{0, 1}, WorkingPrecision -> 50];
{y, z, p} = x + RandomReal[{0, 1}, 3, WorkingPrecision -> 50];
a = y - x;
b = z - x;
c = p - x;
u = a*b/x;
v = a*z/x;
w = b*y/x;
q = u + c;
On[Assert];
Assert[CarlsonRC[x, y] + CarlsonRC[a^2/x, a*y/x] == Pi/(2*Sqrt[a])];
Assert[CarlsonRF[x, y, z] + CarlsonRF[u, v, w] == CarlsonRF[0, a, b]];
Assert[CarlsonRG[x, y, z] + CarlsonRG[u, v, w] == CarlsonRG[0, a, b]
+ (1/2)*((x - u)*CarlsonRF[x, y, z] + u*CarlsonRF[0, a, b] + Sqrt[y*z/x])];
Assert[CarlsonRJ[x, y, z, p] + CarlsonRJ[u, v, w, q] == CarlsonRJ[0, a, b, c]
- 3*CarlsonRC[c^2*y*z/x, c*p*q]];
Assert[CarlsonRD[x, y, z] + CarlsonRD[u, v, w] == CarlsonRD[0, a, b]
- 3/(b*Sqrt[y*z/x])];
Assert[CarlsonRD[x, z, y] + CarlsonRD[u, w, v] == CarlsonRD[0, b, a]
- 3/(a*Sqrt[y*z/x])];
Assert[a*b*(CarlsonRD[y, z, x] + CarlsonRD[v, w, u]) == -6*CarlsonRG[0, a, b]
+ 3*Sqrt[y*z/x]];Addition theorems (General case):
Module[{x, y, z, p, λ, μ, ν, γ, δ, rfλ, rfμ, rfν, rdλ, rdμ, rdν, rgλ, rgμ, rgν, rjλ, rjμ, rjν},
{x, y, z} = RandomReal[{0, 1}, 3, WorkingPrecision -> 50];
{λ, μ} = x + y + z + RandomReal[{0, 1}, 2, WorkingPrecision -> 50];
ν = ((Sqrt[(x + λ)*(y + λ)*(z + λ)] - Sqrt[(x + μ)*(y + μ)*(z + μ)])/(λ - μ))^2 - λ - μ - x - y - z;
Assert[Sqrt[z + ν]
== (Sqrt[(x + λ)*(y + λ)*(z + μ)] - Sqrt[(x + μ)*(y + μ)*(z + λ)])/(λ - μ)
== ((z + λ)*(z + μ) - (z - x)*(z - y))/(Sqrt[(x + λ)*(y + λ)*(z + μ)] + Sqrt[(x + μ)*(y + μ)*(z + λ)])
];
p = Min[x, y, z] + RandomReal[{0, 1}, WorkingPrecision -> 50];
rfλ = CarlsonRF[x + λ, y + λ, z + λ];
rfμ = CarlsonRF[x + μ, y + μ, z + μ];
rfν = CarlsonRF[x + ν, y + ν, z + ν];
rdλ = CarlsonRD[x + λ, y + λ, z + λ];
rdμ = CarlsonRD[x + μ, y + μ, z + μ];
rdν = CarlsonRD[x + ν, y + ν, z + ν];
rgλ = CarlsonRG[x + λ, y + λ, z + λ];
rgμ = CarlsonRG[x + μ, y + μ, z + μ];
rgν = CarlsonRG[x + ν, y + ν, z + ν];
rjλ = CarlsonRJ[x + λ, y + λ, z + λ, p + λ];
rjμ = CarlsonRJ[x + μ, y + μ, z + μ, p + μ];
rjν = CarlsonRJ[x + ν, y + ν, z + ν, p + ν];
Assert[rfλ + rfμ == rfν];
Assert[rdλ + rdμ == rdν - 3/Sqrt[(z + λ)*(z + μ)*(z + ν)]];
Assert[rgλ + rgμ == rgν + (1/2)*(λ*rfλ + μ*rfμ - ν*rfν + Sqrt[λ + μ + ν + x + y + z])];
γ = (p + λ)*(p + μ)*(p + ν);
δ = (p - x)*(p - y)*(p - z);
Assert[rjλ + rjμ == rjν - 3*CarlsonRC[γ - δ, γ]];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Carlson form ↔ Legendre form:
Module[{x, y, z, p, Δ, ϕ, m, m1, n, pr},
Δ = Sqrt[1 - m*Sin[#]^2] &;
pr = 50;
ϕ = RandomReal[{0, Pi/2}, WorkingPrecision -> pr];
m = RandomReal[{0, 1}, WorkingPrecision -> pr];
m1 = 1 - m;
n = 1 - 1/RandomReal[{0, 1}, WorkingPrecision -> pr];
Assert[EllipticF[ϕ, m] ==
+ Sin[ϕ]*CarlsonRF[Cos[ϕ]^2, Δ[ϕ]^2, 1]
];
Assert[EllipticE[ϕ, m] ==
+ Sin[ϕ]*CarlsonRF[Cos[ϕ]^2, Δ[ϕ]^2, 1]
- (1/3)*m*Sin[ϕ]^3*CarlsonRD[Cos[ϕ]^2, Δ[ϕ]^2, 1]
];
Assert[EllipticPi[n, ϕ, m] ==
+ Sin[ϕ]*CarlsonRF[Cos[ϕ]^2, Δ[ϕ]^2, 1]
+ (1/3)*n*Sin[ϕ]^3*CarlsonRJ[Cos[ϕ]^2, Δ[ϕ]^2, 1, 1 - n*Sin[ϕ]^2]
];
Assert[JacobiZeta[ϕ, m] ==
+ (1/3)*m*Sin[ϕ]*Cos[ϕ]*Δ[ϕ]*CarlsonRJ[0, m1, 1, Δ[ϕ]^2]/EllipticK[m]
];
{x, y, z, p} = {Cos[ϕ]^2, Δ[ϕ]^2, 1, 1 - n*Sin[ϕ]^2}*
RandomReal[{0, 1}, WorkingPrecision -> pr];
Assert[ϕ == ArcCos@Sqrt[x/z]];
Assert[m == (z - y)/(z - x)];
Assert[m1 == (y - x)/(z - x)];
Assert[n == (z - p)/(z - x)];
Assert[m - n == (p - y)/(z - x)];
Assert[1 - n == (p - x)/(z - x)];
Assert[CarlsonRC[x, z] ==
(z - x)^(-1/2)*ϕ];
Assert[CarlsonRF[x, y, z] ==
(z - x)^(-1/2)*EllipticF[ϕ, m]];
Assert[CarlsonRD[x, y, z] ==
(3/m)*(z - x)^(-3/2)*(EllipticF[ϕ, m] - EllipticE[ϕ, m])];
Assert[CarlsonRJ[x, y, z, p] ==
(3/n)*(z - x)^(-3/2)*(EllipticPi[n, ϕ, m] - EllipticF[ϕ, m])];
Assert[2*CarlsonRG[x, y, z] ==
(z - x)^(1/2)*(EllipticE[ϕ, m] + Cot[ϕ]^2*EllipticF[ϕ, m] + Cot[ϕ]*Δ[ϕ])];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Rhombic lattice → Rectangular lattice:
Module[{a0, a1, a2, ap, aq, b0, b1, b2, bp, bq, Δ, ϕ, ψ, m, m1, n, pr, rc, rd, rf, rg, rj, implies},
{rc, rd, rf, rg, rj} = Function[{f}, Module[{args, res},
args = {##};
res = f @@ SetAccuracy[args, Accuracy@args + 5];
SetAccuracy[res, Accuracy@res - 5]
] &] /@ {CarlsonRC, CarlsonRD, CarlsonRF, CarlsonRG, CarlsonRJ};
implies = Function[Null, !TrueQ[#1] || #2, HoldAll];
Δ = Sqrt[1 - m*Sin[#]^2] &;
pr = 50;
a0 = RandomReal[{0, 1}, WorkingPrecision -> pr];
a1 = RandomComplex[{0, 1 + I}, WorkingPrecision -> pr];
a2 = Conjugate[a1];
ap = RandomReal[{0, 1}, WorkingPrecision -> pr];
aq = a0 + Abs[a0 - a1]^2/(ap - a0);
b0 = (1/2)*(Abs[a1] + Re[a1]);
b1 = (1/2)*(Abs[a1] + Abs[a0 - a1] + a0);
b2 = (1/2)*(Abs[a1] - Abs[a0 - a1] + a0);
bp = (1/2)*(Abs[a1] + Abs[ap - a1] + ap);
bq = (1/2)*(Abs[a1] - Abs[ap - a1] + ap);
ϕ = 2*ArcTan[Abs[Sqrt[a0] + Sqrt[a1]], Sqrt@Abs[a0 - a1]];
ψ = 2*ArcTan[Sqrt[a0], Sqrt@Abs[a0 - a1]];
{m, m1} = ReIm@Sqrt[a0 - a1]^2/Abs[a0 - a1];
n = (b1 - bp)/Abs[a0 - a1];
Assert[0 <= b2 <= b0 <= b1];
Assert[0 <= bq <= b0 <= bp];
Assert[Sqrt[b0] == Re@Sqrt[a1]];
Assert[Sqrt[b1] == (1/2)*(Abs[Sqrt[a0] + Sqrt[a1]] + Abs[Sqrt[a0] - Sqrt[a1]])];
Assert[Sqrt[b2] == (1/2)*(Abs[Sqrt[a0] + Sqrt[a1]] - Abs[Sqrt[a0] - Sqrt[a1]])];
Assert[Sqrt[bp] == (1/2)*(Abs[Sqrt[ap] + Sqrt[a1]] + Abs[Sqrt[ap] - Sqrt[a1]])];
Assert[Sqrt[bq] == (1/2)*(Abs[Sqrt[ap] + Sqrt[a1]] - Abs[Sqrt[ap] - Sqrt[a1]])];
Assert[Sqrt[b1] + Sqrt[b2] == Abs[Sqrt[a0] + Sqrt[a1]]];
Assert[Sqrt[b1] - Sqrt[b2] == Abs[Sqrt[a0] - Sqrt[a1]]];
Assert[(bp - b0)*(b0 - bq) == (b1 - b0)*(b0 - b2)];
Assert[a0 == b1*b2/b0];
Assert[ap == bp*bq/b0];
Assert[Re[a1] == 2*b0 - Abs[a1]];
Assert[Abs@Im[a1] == 2*Sqrt[b1 - b0]*Sqrt[b0 - b2]];
Assert[Abs[a1] == b2 + b1 - a0];
Assert[Abs[a0 - a1] == b1 - b2];
Assert[Re@Sqrt[a0 - a1]^2 == b1 - b0];
Assert[Im@Sqrt[a0 - a1]^2 == b0 - b2];
Assert[Abs[ap - a1] == bp - bq];
Assert[Re@Sqrt[ap - a1]^2 == bp - b0];
Assert[Im@Sqrt[ap - a1]^2 == b0 - bq];
Assert[Sin[ϕ] == Sqrt[Abs[a0 - a1]/b1]];
Assert[Cos[ϕ] == Sqrt[b2/b1]];
Assert[Δ[ϕ] == Sqrt[b0/b1]];
Assert[Cos[ϕ]/Δ[ϕ] == Sqrt[a0/b1]];
Assert[Tan[ϕ]*Δ[ϕ] == Tan[ψ/2] == Sqrt[Abs[a0 - a1]/a0]];
Assert[Sin[ψ] == 2*Sqrt[a0*Abs[a0 - a1]]/(a0 + Abs[a0 - a1])];
Assert[Cos[ψ] == (a0 - Abs[a0 - a1])/(a0 + Abs[a0 - a1])];
Assert[Δ[ψ] == Abs[a1]/(a0 + Abs[a0 - a1])];
Assert[1 - m*(Sin[ϕ]*Cos[ϕ]/Δ[ϕ])^2 == Abs[a1]/b1];
Assert[Csc[ϕ]^2 - m*(Cos[ϕ]/Δ[ϕ])^2 == Abs[a1]/Abs[a0 - a1]];
Assert[Max[0, n] <= m <= 1];
Assert[m - n == (bp - b0)/Abs[a0 - a1]];
Assert[1 - n == (bp - b2)/Abs[a0 - a1]];
Assert[m*m1/(m - n) == (b0 - bq)/Abs[a0 - a1]];
Assert[m*(1 - n)/(m - n) == (b1 - bq)/Abs[a0 - a1]];
Assert[m1*n/(m - n) == (b2 - bq)/Abs[a0 - a1]];
(* CarlsonRF *)
Assert[rf[a1, a2, a0]
== rf[b1, b2, b0]
== Abs[a0 - a1]^(-1/2)*EllipticF[ϕ, m]
== (1/2)*Abs[a0 - a1]^(-1/2)*EllipticF[ψ, m]
];
(* CarlsonRD *)
Assert[Abs[a0 - a1]*rd[a1, a2, a0]
== 2*(b1 - b0)*rd[b0, b2, b1] - 3*rf[b1, b2, b0] + 3/Sqrt[a0]
== 2*(b0 - b2)*rd[b0, b1, b2] + 3*rf[b1, b2, b0] - 3/Sqrt[a0]
== 3*Abs[a0 - a1]^(-1/2)*(
- 2*EllipticE[ϕ, m]
+ EllipticF[ϕ, m]
+ Tan[ϕ]*Δ[ϕ]
)
== 3*Abs[a0 - a1]^(-1/2)*(
- EllipticE[ψ, m]
+ (1/2)*EllipticF[ψ, m]
+ Δ[ψ]*Tan[ψ/2]
)
];
(* CarlsonRG *)
Assert[rg[a1, a2, a0]
== 2*rg[b1, b2, b0] - (1/2)*(Abs[a1]*rf[b1, b2, b0] + Sqrt[a0])
== (1/12)*Im[a1]^2*rd[b1, b2, b0] + (1/2)*(Re[a1]*rf[b1, b2, b0] + Sqrt[a0])
== (1/2)*Abs[a0 - a1]^(1/2)*(
+ 2*EllipticE[ϕ, m]
+ (Cot[ϕ]^2 - m1/Δ[ϕ]^2)*EllipticF[ϕ, m]
+ Cot[ϕ]*(1 - 2*m*Sin[ϕ]^2)/Δ[ϕ]
)
== (1/2)*Abs[a0 - a1]^(1/2)*(
+ EllipticE[ψ, m]
+ Cot[ψ/2]*(Cot[ψ]*EllipticF[ψ, m] + Δ[ψ])
)
];
(* CarlsonRJ *)
Assert[Abs[ap - a1]*rj[a1, a2, a0, ap]
== 2*(bp - b0)*rj[b1, b2, b0, bp] - 3*rf[b1, b2, b0] + 3*rc[a0, ap]
== 2*(b0 - bq)*rj[b1, b2, b0, bq] + 3*rf[b1, b2, b0] - 3*rc[a0, ap]
== 3*Abs[a0 - a1]^(-1/2)*(
+ 2*(m/n - 1)*EllipticPi[n, ϕ, m]
+ (1 - 2*m/n)*EllipticF[ϕ, m]
) + 3*rc[a0, ap]
];
Assert@Implies[0 < aq, With[{γ = n*(1 - n), δ = (m - n^2)},
(ap - a0)*rj[a0, a1, a2, ap]
== (3/2)*Abs[a0 - a1]^(-1/2)*(
+ (1/2 - γ/δ)*EllipticPi[δ^2/(4*γ*(m - n)), ψ, m]
+ (γ/δ)*EllipticF[ψ, m]
) - (3/2)*rc[Abs[a1]^2/a0, ap*aq/a0]
]];
(* Branch cut of CarlsonRF *)
Assert@implies[a0 <= 0,
Re@rf[a1, a2, a0] == rf[0, Re@Sqrt[a1 - a0]^2, Abs[a1 - a0]] &&
Abs@Im@rf[a1, a2, a0] == Sqrt[-a0]*rf[Abs[a1]^2 - a0*a1, Abs[a1]^2 - a0*a2, Abs[a1 - a0]^2]
];
(* Branch cut of CarlsonRD *)
Assert@implies[a0 < 0,
Re@rd[a1, a2, a0] == -Abs[a1 - a0]^(-1)*(
+ 3*rf[0, Re@Sqrt[a1 - a0]^2, Abs[a1 - a0]]
- 2*Im@Sqrt[a1 - a0]^2*rd[0, Re@Sqrt[a1 - a0]^2, Abs[a1 - a0]]
) &&
Abs@Im@rd[a1, a2, a0] == Abs[
+ 3*Abs[a1]/(Sqrt[-a0]*Abs[a1 - a0]^2)
- (-a0)^(3/2)*rd[Abs[a1]^2 - a0*a1, Abs[a1]^2 - a0*a2, Abs[a1 - a0]^2]
]
];
(* Branch cut of CarlsonRG *)
Assert@implies[a0 <= 0,
Re@rg[a1, a2, a0] == (
+ (1/2)*(a0 + Abs[a1 - a0])*rf[0, Re@Sqrt[a1 - a0]^2, Abs[a1 - a0]]
- (1/3)*Abs[a1 - a0]*Im@Sqrt[a1 - a0]^2*rd[0, Re@Sqrt[a1 - a0]^2, Abs[a1 - a0]]
) &&
Abs@Im@rg[a1, a2, a0] == (1/2)*(-a0)^(3/2)*Abs[
+ rf[Abs[a1]^2 - a0*a1, Abs[a1]^2 - a0*a2, Abs[a1 - a0]^2]
- (1/3)*Abs[a1 - a0]^2*rd[Abs[a1]^2 - a0*a1, Abs[a1]^2 - a0*a2, Abs[a1 - a0]^2]
]
];
(* Inequalities *)
Assert@implies[a0 >= 0 && Re[a1] > 0,
rc[a0, Abs[a1]] <= Re@rf[a1, a2, a0] <= rc[a0, Re[a1]] &&
implies[0 < a0 < Re[a1],
rc[a0, Re[a1]] == Abs[a0 - a1]^(-1/2)*(m1 - m)^(-1/2)*ArcTan[Sqrt[m1 - m]*Tan[ϕ]*Δ[ϕ]]] &&
implies[a0 > Re[a1],
rc[a0, Re[a1]] == Abs[a0 - a1]^(-1/2)*(m - m1)^(-1/2)*ArcTanh[Sqrt[m - m1]*Tan[ϕ]*Δ[ϕ]]]
];
Assert@implies[0 < a0 < Re[a1],
Re@rf[a1, a2, a0] >= (2/Pi)*Abs[a0 - a1]^(-1/2)*EllipticK[m]*ArcCos@Sqrt[a0/Re[a1]]
];
Assert@implies[a0 >= Re[a1],
Re@rf[a1, a2, a0] >= a0^(-1/2)*(1 - (Im[a1]/a0)^2/10)
];
Assert@implies[0 < a0 <= Re[a1],
Re@rf[a1, a2, a0] <= a0^(-1/2)*(1 - Arg[a1]^2/10)
];
Assert[Re@rf[a1, a2, a0] < (Log[Pi/Abs@Arg[-a1]] + (1/2)*Log[E^Pi + 4*a0/Abs[a1]])/Sqrt@Abs[a0 - a1]];
] & // Block[{$AssertFunction = Automatic}, #[]] &;AGM method:
Module[{pr, ϕ, m, f},
pr = 50;
ϕ = RandomReal[{0, Pi/2}, WorkingPrecision -> pr];
m = RandomReal[{0, 1}, WorkingPrecision -> pr];
f = Module[{ϕ, m, m1, n, a, b, a1, b1, r, t},
ϕ = #1;
m = #2;
m1 = 1 - m;
n = Round[Re[ϕ]/Pi];
ϕ -= n*Pi;
t = Tan[ϕ/2];
r = t^2;
m1 = 1 - m;
a = 1 + (m1 - m)*r;
b = Sqrt[a^2 + 4*m*m1*r^2];
If[ReIm[a].ReIm[b] < 0, b = -b];
r = Sqrt[1 + 2*m*r/(a + b)];
a = Sqrt[m1];
b = 1;
While[
a != b || r != 1,
t *= r;
a1 = a;
b1 = b;
a = (a1 + b1)/2;
b = Sqrt[a1*b1];
r = Sqrt[(b1 + a1*r)/(a*(1 + r))];
];
Return[(2*ArcTan[a*t] + n*Pi)/a, Module];
] &;
Assert[EllipticF[ϕ, m] == f[ϕ, m]];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Module[{Δ, Δ1, m, k, kp, u, v, d, c, ϕu, ϕv, ϕP, ϕM, m1, am, sn, cn, dn, tn, Es, Ec, Ed, En, mod, pr},
mod[a_, b_, m_] := With[{t = (a - b)/m},
!NumericQ[t] || With[{tr = Re[t], ti = Im[t]},
PossibleZeroQ[ti] && PossibleZeroQ[tr - Round[tr]]
]
];
pr = 50;
Δ = Sqrt[1 - m*Sin[#]^2] &;
Δ1 = Sqrt[1 - m1*Sin[#]^2] &;
Es = EllipticE[#, m] + Cot[#]*Δ[#] &;
Ec = EllipticE[#, m] - Tan[#]*Δ[#] &;
Ed = EllipticE[#, m] - m*Sin[#]*Cos[#]/Δ[#] &;
En = EllipticE[#, m] &;
sn = JacobiSN[#, m] &;
cn = JacobiCN[#, m] &;
dn = JacobiDN[#, m] &;
tn = JacobiSC[#, m] &;
am = If[Quiet@TrueQ[Cos@N[JacobiAmplitude[2*I*EllipticK[1/2], 1/2], pr] == -1],
(* built-in *)
JacobiAmplitude[#, m] &,
(* fixup *)
Function[Null, Module[{z = #, w, n},
w = 2*EllipticK[m];
n = Round@Re[z/w];
z -= n*w;
Pi*n + ArcTan[cn[z], sn[z]]
], Listable]
];
m = RandomReal[WorkingPrecision -> pr];
m1 = 1 - m;
k = Sqrt[m];
kp = Sqrt[m1];
{u, v, c} = RandomReal[{0, EllipticK[m]}, 3, WorkingPrecision -> pr];
{ϕu, ϕv} = am@{u, v};
(* Addition formulas *)
{ϕP, ϕM} = am@{u + v, u - v};
Assert@mod[ϕP, ArcTan[Cos[ϕu], Sin[ϕu]*Δ[ϕv]] + ArcTan[Cos[ϕv], Sin[ϕv]*Δ[ϕu]], 2*Pi];
Assert@mod[ϕM, ArcTan[Cos[ϕu], Sin[ϕu]*Δ[ϕv]] - ArcTan[Cos[ϕv], Sin[ϕv]*Δ[ϕu]], 2*Pi];
Assert[ϕP == ϕu + ϕv - 2*ArcTan[1 + Δ[ϕu] + Δ[ϕv] + Δ[ϕP], m*Sin[ϕu]*Sin[ϕv]*Sin[ϕP]]];
Assert[ϕM == ϕu - ϕv + 2*ArcTan[1 + Δ[ϕu] + Δ[ϕv] + Δ[ϕM], m*Sin[ϕu]*Sin[ϕv]*Sin[ϕM]]];
d = 1 - m*Sin[ϕu]^2*Sin[ϕv]^2;
Assert[Sin[ϕP] == (Sin[ϕu]*Cos[ϕv]*Δ[ϕv] + Sin[ϕv]*Cos[ϕu]*Δ[ϕu])/d];
Assert[Sin[ϕM] == (Sin[ϕu]*Cos[ϕv]*Δ[ϕv] - Sin[ϕv]*Cos[ϕu]*Δ[ϕu])/d];
Assert[Cos[ϕP] == (Cos[ϕu]*Cos[ϕv] - Sin[ϕu]*Sin[ϕv]*Δ[ϕu]*Δ[ϕv])/d];
Assert[Cos[ϕM] == (Cos[ϕu]*Cos[ϕv] + Sin[ϕu]*Sin[ϕv]*Δ[ϕu]*Δ[ϕv])/d];
Assert[Δ[ϕP] == (Δ[ϕu]*Δ[ϕv] - m*Sin[ϕu]*Sin[ϕv]*Cos[ϕu]*Cos[ϕv])/d];
Assert[Δ[ϕM] == (Δ[ϕu]*Δ[ϕv] + m*Sin[ϕu]*Sin[ϕv]*Cos[ϕu]*Cos[ϕv])/d];
Assert[Sin[ϕP]*Sin[ϕM] == (Sin[ϕu]^2 - Sin[ϕv]^2)/d];
Assert[Sin[ϕP]*Cos[ϕM] == (Sin[ϕu]*Cos[ϕu]*Δ[ϕv] + Sin[ϕv]*Cos[ϕv]*Δ[ϕu])/d];
Assert[Sin[ϕP]*Δ[ϕM] == (Sin[ϕu]*Cos[ϕv]*Δ[ϕu] + Sin[ϕv]*Cos[ϕu]*Δ[ϕv])/d];
Assert[Cos[ϕP]*Sin[ϕM] == (Sin[ϕu]*Cos[ϕu]*Δ[ϕv] - Sin[ϕv]*Cos[ϕv]*Δ[ϕu])/d];
Assert[Cos[ϕP]*Cos[ϕM] == (1 - Sin[ϕu]^2 - Sin[ϕv]^2 + m*Sin[ϕu]^2*Sin[ϕv]^2)/d];
Assert[Cos[ϕP]*Δ[ϕM] == (Cos[ϕu]*Cos[ϕv]*Δ[ϕu]*Δ[ϕv] - m1*Sin[ϕu]*Sin[ϕv])/d];
Assert[Δ[ϕP]*Sin[ϕM] == (Sin[ϕu]*Cos[ϕv]*Δ[ϕu] - Sin[ϕv]*Cos[ϕu]*Δ[ϕv])/d];
Assert[Δ[ϕP]*Cos[ϕM] == (Cos[ϕu]*Cos[ϕv]*Δ[ϕu]*Δ[ϕv] + m1*Sin[ϕu]*Sin[ϕv])/d];
Assert[Δ[ϕP]*Δ[ϕM] == (1 - m*Sin[ϕu]^2 - m*Sin[ϕv]^2 + m*Sin[ϕu]^2*Sin[ϕv]^2)/d];
Assert[Sin[ϕP] + Sin[ϕM] == 2*Sin[ϕu]*Cos[ϕv]*Δ[ϕv]/d];
Assert[Sin[ϕP] - Sin[ϕM] == 2*Sin[ϕv]*Cos[ϕu]*Δ[ϕu]/d];
Assert[Cos[ϕP] + Cos[ϕM] == 2*Cos[ϕu]*Cos[ϕv]/d];
Assert[Cos[ϕP] - Cos[ϕM] == -2*Sin[ϕu]*Sin[ϕv]*Δ[ϕu]*Δ[ϕv]/d];
Assert[Δ[ϕP] + Δ[ϕM] == 2*Δ[ϕu]*Δ[ϕv]/d];
Assert[Δ[ϕP] - Δ[ϕM] == -2*m*Sin[ϕu]*Sin[ϕv]*Cos[ϕu]*Cos[ϕv]/d];
Assert[Sin[ϕP] + Cos[ϕM] == (Cos[ϕu] + Sin[ϕu]*Δ[ϕv])*(Cos[ϕv] + Sin[ϕv]*Δ[ϕu])/d];
Assert[Sin[ϕP] - Cos[ϕM] == -(Cos[ϕu] - Sin[ϕu]*Δ[ϕv])*(Cos[ϕv] - Sin[ϕv]*Δ[ϕu])/d];
Assert[k*Sin[ϕP] + Δ[ϕM] == (Δ[ϕu] + k*Sin[ϕu]*Cos[ϕv])*(Δ[ϕv] + k*Sin[ϕv]*Cos[ϕu])/d];
Assert[k*Sin[ϕP] - Δ[ϕM] == -(Δ[ϕu] - k*Sin[ϕu]*Cos[ϕv])*(Δ[ϕv] - k*Sin[ϕv]*Cos[ϕu])/d];
Assert[Sin[ϕP + ϕM] == Sin[2*ϕu]*Δ[ϕv]/d];
Assert[Sin[ϕP - ϕM] == Sin[2*ϕv]*Δ[ϕu]/d];
Assert[Cos[ϕP + ϕM] == (Cos[2*ϕu] + m*Sin[ϕu]^2*Sin[ϕv]^2)/d];
Assert[Cos[ϕP - ϕM] == (Cos[2*ϕv] + m*Sin[ϕu]^2*Sin[ϕv]^2)/d];
Assert[(1 + Sin[ϕP])*(1 + Sin[ϕM]) == (Cos[ϕv] + Sin[ϕu]*Δ[ϕv])^2/d];
Assert[(1 - Sin[ϕP])*(1 - Sin[ϕM]) == (Cos[ϕv] - Sin[ϕu]*Δ[ϕv])^2/d];
Assert[(1 + Sin[ϕP])*(1 - Sin[ϕM]) == (Cos[ϕu] + Sin[ϕv]*Δ[ϕu])^2/d];
Assert[(1 - Sin[ϕP])*(1 + Sin[ϕM]) == (Cos[ϕu] - Sin[ϕv]*Δ[ϕu])^2/d];
Assert[(1 + k*Sin[ϕP])*(1 + k*Sin[ϕM]) == (Δ[ϕv] + k*Sin[ϕu]*Cos[ϕv])^2/d];
Assert[(1 - k*Sin[ϕP])*(1 - k*Sin[ϕM]) == (Δ[ϕv] - k*Sin[ϕu]*Cos[ϕv])^2/d];
Assert[(1 + k*Sin[ϕP])*(1 - k*Sin[ϕM]) == (Δ[ϕu] + k*Sin[ϕv]*Cos[ϕu])^2/d];
Assert[(1 - k*Sin[ϕP])*(1 + k*Sin[ϕM]) == (Δ[ϕu] - k*Sin[ϕv]*Cos[ϕu])^2/d];
Assert[(1 + Cos[ϕP])*(1 + Cos[ϕM]) == (Cos[ϕu] + Cos[ϕv])^2/d];
Assert[(1 - Cos[ϕP])*(1 - Cos[ϕM]) == (Cos[ϕu] - Cos[ϕv])^2/d];
Assert[(1 + Cos[ϕP])*(1 - Cos[ϕM]) == (Sin[ϕu]*Δ[ϕv] - Sin[ϕv]*Δ[ϕu])^2/d];
Assert[(1 - Cos[ϕP])*(1 + Cos[ϕM]) == (Sin[ϕu]*Δ[ϕv] + Sin[ϕv]*Δ[ϕu])^2/d];
Assert[(1 + Δ[ϕP])*(1 + Δ[ϕM]) == (Δ[ϕu] + Δ[ϕv])^2/d];
Assert[(1 - Δ[ϕP])*(1 - Δ[ϕM]) == (Δ[ϕu] - Δ[ϕv])^2/d];
Assert[(1 + Δ[ϕP])*(1 - Δ[ϕM]) == m*Sin[ϕu - ϕv]^2/d];
Assert[(1 - Δ[ϕP])*(1 + Δ[ϕM]) == m*Sin[ϕu + ϕv]^2/d];
Assert[(Δ[ϕP] + Cos[ϕP])*(Δ[ϕM] + Cos[ϕM]) == (Cos[ϕu]*Δ[ϕv] + Cos[ϕv]*Δ[ϕu])^2/d];
Assert[(Δ[ϕP] - Cos[ϕP])*(Δ[ϕM] - Cos[ϕM]) == (Cos[ϕu]*Δ[ϕv] - Cos[ϕv]*Δ[ϕu])^2/d];
Assert[(Δ[ϕP] + Cos[ϕP])*(Δ[ϕM] - Cos[ϕM]) == m1*(Sin[ϕu] - Sin[ϕv])^2/d];
Assert[(Δ[ϕP] - Cos[ϕP])*(Δ[ϕM] + Cos[ϕM]) == m1*(Sin[ϕu] + Sin[ϕv])^2/d];
Assert[Δ[ϕP] + k*Cos[ϕP] == (Δ[ϕu]*Δ[ϕv] + k*Cos[ϕu]*Cos[ϕv])/(1 + k*Sin[ϕu]*Sin[ϕv])];
Assert[Δ[ϕP] - k*Cos[ϕP] == (Δ[ϕu]*Δ[ϕv] - k*Cos[ϕu]*Cos[ϕv])/(1 - k*Sin[ϕu]*Sin[ϕv])];
Assert[Δ[ϕM] + k*Cos[ϕM] == (Δ[ϕu]*Δ[ϕv] + k*Cos[ϕu]*Cos[ϕv])/(1 - k*Sin[ϕu]*Sin[ϕv])];
Assert[Δ[ϕM] - k*Cos[ϕM] == (Δ[ϕu]*Δ[ϕv] - k*Cos[ϕu]*Cos[ϕv])/(1 + k*Sin[ϕu]*Sin[ϕv])];
Assert[Sin[(ϕP + ϕM)/2] == Sin[ϕu]*Δ[ϕv]/Sqrt[d]];
Assert[Sin[(ϕP - ϕM)/2] == Sin[ϕv]*Δ[ϕu]/Sqrt[d]];
Assert[Cos[(ϕP + ϕM)/2] == Cos[ϕu]/Sqrt[d]];
Assert[Cos[(ϕP - ϕM)/2] == Cos[ϕv]/Sqrt[d]];
Assert[Tan[(ϕP + ϕM)/2] == Tan[ϕu]*Δ[ϕv]];
Assert[Tan[(ϕP - ϕM)/2] == Tan[ϕv]*Δ[ϕu]];
Assert[Tan[ϕP/2] == (Sin[ϕu]*Δ[ϕv] + Sin[ϕv]*Δ[ϕu])/(Cos[ϕu] + Cos[ϕv])];
Assert[Tan[ϕM/2] == (Sin[ϕu]*Δ[ϕv] - Sin[ϕv]*Δ[ϕu])/(Cos[ϕu] + Cos[ϕv])];
Assert[Sin[ϕP]/(1 + Δ[ϕP]) == Sin[ϕu + ϕv]/(Δ[ϕu] + Δ[ϕv])];
Assert[Sin[ϕM]/(1 + Δ[ϕM]) == Sin[ϕu - ϕv]/(Δ[ϕu] + Δ[ϕv])];
Assert[Tan[ϕP] == (Tan[ϕu]*Δ[ϕv] + Tan[ϕv]*Δ[ϕu])/(1 - Tan[ϕu]*Tan[ϕv]*Δ[ϕu]*Δ[ϕv])];
Assert[Tan[ϕM] == (Tan[ϕu]*Δ[ϕv] - Tan[ϕv]*Δ[ϕu])/(1 + Tan[ϕu]*Tan[ϕv]*Δ[ϕu]*Δ[ϕv])];
Assert[Es[ϕP] == Es[ϕu] + Es[ϕv] - Csc[ϕu]*Csc[ϕv]*Csc[ϕP] + Cot[ϕu]*Cot[ϕv]*Cot[ϕP]];
Assert[Es[ϕM] == Es[ϕu] - Es[ϕv] + Csc[ϕu]*Csc[ϕv]*Csc[ϕM] - Cot[ϕu]*Cot[ϕv]*Cot[ϕM]];
Assert[Ec[ϕP] == Ec[ϕu] + Ec[ϕv] - m1*Tan[ϕu]*Tan[ϕv]*Tan[ϕP]];
Assert[Ec[ϕM] == Ec[ϕu] - Ec[ϕv] + m1*Tan[ϕu]*Tan[ϕv]*Tan[ϕM]];
Assert[Ed[ϕP] == Ed[ϕu] + Ed[ϕv] + m*m1*Sin[ϕu]*Sin[ϕv]*Sin[ϕP]/(m1 + m*Cos[ϕu]*Cos[ϕv]*Cos[ϕP])];
Assert[Ed[ϕM] == Ed[ϕu] - Ed[ϕv] - m*m1*Sin[ϕu]*Sin[ϕv]*Sin[ϕM]/(m1 + m*Cos[ϕu]*Cos[ϕv]*Cos[ϕM])];
Assert[En[ϕP] == En[ϕu] + En[ϕv] - m*Sin[ϕu]*Sin[ϕv]*Sin[ϕP]];
Assert[En[ϕM] == En[ϕu] - En[ϕv] + m*Sin[ϕu]*Sin[ϕv]*Sin[ϕM]];
d = (1/2)*Log[(1 + m*Sin[ϕu]*Sin[ϕv]*sn[c]*sn[u + v + c])/(1 - m*Sin[ϕu]*Sin[ϕv]*sn[c]*sn[u + v - c])];
Assert[EllipticPi[m*sn[c]^2, ϕP, m] == EllipticPi[m*sn[c]^2, ϕu, m] + EllipticPi[m*sn[c]^2, ϕv, m] + (tn[c]/dn[c])*d];
With[{
x = m*sn[c]*cn[c]*dn[c]*Sin[ϕu]*Sin[ϕv]*Sin[ϕP],
y = dn[c]^2 + m*sn[c]^2*Cos[ϕu]*Cos[ϕv]*Cos[ϕP],
r = Sqrt[Times @@ ((1 - m*sn[c]^2*Sin[#]^2) & /@ {ϕu, ϕv, ϕP})]
},
Assert[d == +ArcTanh[x/y] == +ArcSinh[x/r]];
];
d = (1/2)*Log[(1 - m*Sin[ϕu]*Sin[ϕv]*sn[c]*sn[u - v + c])/(1 + m*Sin[ϕu]*Sin[ϕv]*sn[c]*sn[u - v - c])];
Assert[EllipticPi[m*sn[c]^2, ϕM, m] == EllipticPi[m*sn[c]^2, ϕu, m] - EllipticPi[m*sn[c]^2, ϕv, m] + (tn[c]/dn[c])*d];
With[{
x = m*sn[c]*cn[c]*dn[c]*Sin[ϕu]*Sin[ϕv]*Sin[ϕM],
y = dn[c]^2 + m*sn[c]^2*Cos[ϕu]*Cos[ϕv]*Cos[ϕM],
r = Sqrt[Times @@ ((1 - m*sn[c]^2*Sin[#]^2) & /@ {ϕu, ϕv, ϕM})]
},
Assert[d == -ArcTanh[x/y] == -ArcSinh[x/r]];
];
Clear[d, c];
(* Change of parameter, from https://dlmf.nist.gov/19.7.iii *)
Assert[Abs[ϕu] < Pi/2 \[Implies]
EllipticPi[m*Sin[ϕv]^2, ϕu, m] + EllipticPi[Csc[ϕv]^2, ϕu, m] ==
u + (1/2)*(Tan[ϕv]/Δ[ϕv])*Log@Abs[Sin[ϕP]/Sin[ϕM]]];
Assert[Abs[ϕu] < Pi/2 \[Implies]
Δ[ϕv]^2*EllipticPi[m*Sin[ϕv]^2, ϕu, m] - m1*Tan[ϕv]^2*EllipticPi[(Δ[ϕv]/Cos[ϕv])^2, ϕu, m] ==
u + (1/2)*(Tan[ϕv]*Δ[ϕv])*Log@Abs[Cos[ϕP]/Cos[ϕM]]];
Assert[Cos[ϕv]^2*EllipticPi[m*Sin[ϕv]^2, ϕu, m] + m1*(Sin[ϕv]/Δ[ϕv])^2*EllipticPi[m*(Cos[ϕv]/Δ[ϕv])^2, ϕu, m] ==
u + (1/2)*(Sin[ϕv]*Cos[ϕv]/Δ[ϕv])*Log[Δ[ϕP]/Δ[ϕM]]];
Assert[Cot[ϕv]*Δ[ϕv]*(EllipticPi[m*Sin[ϕv]^2, ϕu, m] - u) - Cot[ϕu]*Δ[ϕu]*(EllipticPi[m*Sin[ϕu]^2, ϕv, m] - v) ==
u*En[ϕv] - v*En[ϕu]];
Clear[ϕP, ϕM];
(* Double angle formulas *)
d = 1 - m*Sin[ϕu]^4;
Assert@mod[am[2*u], 2*ArcTan[Cos[ϕu], Sin[ϕu]*Δ[ϕu]], 4*Pi];
Assert[sn[2*u] == Sin[2*ϕu]*Δ[ϕu]/d];
Assert[cn[2*u] == (Cos[2*ϕu] + m*Sin[ϕu]^4)/d];
Assert[dn[2*u] == (m1 + m*Cos[ϕu]^4)/d];
Assert[1 + cn[2*u] == 2*Cos[ϕu]^2/d];
Assert[1 - cn[2*u] == 2*(Sin[ϕu]*Δ[ϕu])^2/d];
Assert[1 + dn[2*u] == 2*Δ[ϕu]^2/d];
Assert[1 - dn[2*u] == m*Sin[2*ϕu]^2/(2*d)];
Assert[En@am[2*u] == 2*(En[ϕu] - m*Sin[ϕu]^3*Cos[ϕu]*Δ[ϕu]/d)];
Assert[Cot[ϕu]*Δ[ϕu]*(EllipticPi[m*Sin[ϕu]^2, ϕu, m] - u) ==
u*JacobiZeta[am[u], m] - 2*Log@NevilleThetaN[u, m] - (1/2)*Log[d]];
Clear[d];
(* Triple angle formulas *)
d = m*(m*m1^3*Sin[ϕu]^8 - m*Cos[ϕu]^8 + Δ[ϕu]^8);
Assert[sn[3*u]/sn[u] == (-m^3*Cos[ϕu]^8 + Δ[ϕu]^8 - m1^3)/d];
Assert[cn[3*u]/cn[u] == m1*(m^3*m1*Sin[ϕu]^8 + Δ[ϕu]^8 - m1)/d];
Assert[dn[3*u]/dn[u] == m*m1*(-m*m1*Sin[ϕu]^8 + m*Cos[ϕu]^8 + m1)/d];
Assert[En@am[3*u] == 3*En[ϕu] - m^2*m1*(Sin[2*ϕu]*Δ[ϕu])^3/d];
Assert[Cot[ϕu]*Δ[ϕu]*(EllipticPi[m*Sin[ϕu]^2, am[2*u], m] - 2*u) ==
2*u*JacobiZeta[am[u], m] - 4*Log@NevilleThetaN[u, m] + (1/2)*Log[m*m1/d]];
Clear[d];
(* Half angle formulas *)
d = 1 + Δ[ϕu];
c = m1 - m*Cos[ϕu] + Δ[ϕu];
Assert[sn[u/2] == Sin[ϕu/2]*Sqrt[2/d]];
Assert[cn[u/2] == kp*Cos[ϕu/2]*Sqrt[2/c]];
Assert[dn[u/2] == kp*Sqrt[d/c]];
Assert[tn[u/2] == Tan[ϕu/2]*Sqrt[c/d]/kp];
Assert[En@am[u/2] == (1/2)*(En[ϕu] + m*Sin[ϕu]*(1 - Cos[ϕu])/(1 + Δ[ϕu]))];
Clear[d, c];
Assert[sn[(u + v)/2] == Sin[(ϕu + ϕv)/2]*Sqrt[2/(1 + m*Sin[ϕu]*Sin[ϕv] + Δ[ϕu]*Δ[ϕv])]];
Assert[cn[(u + v)/2] == Cos[(ϕu + ϕv)/2]*Sqrt[2*m1/(m1 - m*Cos[ϕu]*Cos[ϕv] + Δ[ϕu]*Δ[ϕv])]];
(* Complex arguments *)
v = RandomReal[EllipticK[m1]*{-1, 1}, WorkingPrecision -> pr];
ϕv = JacobiAmplitude[v, m1];
{ϕP, ϕM} = am@{u + I*v, u - I*v};
d = 1 - Sin[ϕv]^2 + m*Sin[ϕu]^2*Sin[ϕv]^2;
Assert@mod[ϕP, ArcTan[Cos[ϕu]*Cos[ϕv], Sin[ϕu]*Δ1[ϕv]] + I*ArcTanh[Δ[ϕu]*Sin[ϕv]], 2*Pi];
Assert[ArcSin[k*Sin[ϕP]] == ArcTan[Δ[ϕu]*Cos[ϕv], k*Sin[ϕu]] + I*ArcTanh[k*Cos[ϕu]*Sin[ϕv]/Δ1[ϕv]]];
Assert[Sin[ϕP] == (Sin[ϕu]*Δ1[ϕv] + I*Cos[ϕu]*Δ[ϕu]*Sin[ϕv]*Cos[ϕv])/d];
Assert[Cos[ϕP] == (Cos[ϕu]*Cos[ϕv] - I*Sin[ϕu]*Δ[ϕu]*Sin[ϕv]*Δ1[ϕv])/d];
Assert[Δ[ϕP] == (Δ[ϕu]*Cos[ϕv]*Δ1[ϕv] - I*m*Sin[ϕu]*Cos[ϕu]*Sin[ϕv])/d];
Assert[Abs@Sin[ϕP]^2 == (Sin[ϕu]^2 + Sin[ϕv]^2 - Sin[ϕu]^2*Sin[ϕv]^2)/d];
Assert[Abs@Cos[ϕP]^2 == (1 - Sin[ϕu]^2 + m1*Sin[ϕu]^2*Sin[ϕv]^2)/d];
Assert[Abs@Δ[ϕP]^2 == (1 - m*Sin[ϕu]^2 - m1*Sin[ϕv]^2)/d];
Assert[Tan[ϕP/2] == (Sin[ϕu]*Δ1[ϕv] + I*Δ[ϕu]*Sin[ϕv])/(1 + Cos[ϕu]*Cos[ϕv])];
Assert[Cot[ϕP/2] == (Sin[ϕu]*Δ1[ϕv] - I*Δ[ϕu]*Sin[ϕv])/(1 - Cos[ϕu]*Cos[ϕv])];
Assert[Tan[Pi/4 + ϕP/2] == (Cos[ϕu]*Cos[ϕv] + I*Δ[ϕu]*Sin[ϕv])/(1 - Sin[ϕu]*Δ1[ϕv])];
Assert[Tan[Pi/4 - ϕP/2] == (Cos[ϕu]*Cos[ϕv] - I*Δ[ϕu]*Sin[ϕv])/(1 + Sin[ϕu]*Δ1[ϕv])];
Assert[E^(+I*ϕP) == (Cos[ϕu]*Cos[ϕv] + I*Sin[ϕu]*Δ1[ϕv])/(1 + Δ[ϕu]*Sin[ϕv])];
Assert[E^(-I*ϕP) == (Cos[ϕu]*Cos[ϕv] - I*Sin[ϕu]*Δ1[ϕv])/(1 - Δ[ϕu]*Sin[ϕv])];
Assert[EllipticE[ϕP, m] == EllipticE[ϕu, m] + I*(v - EllipticE[ϕv, m1]) +
(m*Sin[ϕu]*Cos[ϕu]*Δ[ϕu]*Sin[ϕv]^2 + I*Δ[ϕu]^2*Sin[ϕv]*Cos[ϕv]*Δ1[ϕv])/d];
d = u/EllipticK[m];
c = -Cot[ϕv]^2;
Assert@mod[Arg@NevilleThetaS[u + I*v, m],
-2*Csc[2*ϕv]*Δ1[ϕv]*(EllipticPi[c, ϕu, m] - d*EllipticPi[c, m]) - Sign[v]*(Pi/2)*(d - 1), 2*Pi];
c = Δ1[ϕv]^2;
Assert@mod[Arg@NevilleThetaC[u + I*v, m],
-(m1*Sin[ϕv]*Cos[ϕv]/Δ1[ϕv])*(EllipticPi[c, ϕu, m] - d*EllipticPi[c, m]) - Sign[v]*(Pi/2)*d, 2*Pi];
c = m/Δ1[ϕv]^2;
Assert@mod[Arg@NevilleThetaD[u + I*v, m],
+(m1*Sin[ϕv]*Cos[ϕv]/Δ1[ϕv])*(EllipticPi[c, ϕu, m] - d*EllipticPi[c, m]), 2*Pi];
c = -m*Tan[ϕv]^2;
Assert@mod[Arg@NevilleThetaN[u + I*v, m],
+2*Csc[2*ϕv]*Δ1[ϕv]*(EllipticPi[c, ϕu, m] - d*EllipticPi[c, m]), 2*Pi];
Clear[d, c, ϕP, ϕM];
{ϕP, ϕM} = am[{u + I*v, u - I*v}/2];
d = Δ[ϕu]*Δ1[ϕv] - m*Cos[ϕu] + m1*Cos[ϕv];
Assert[Sin[ϕP]^2 == (Δ[ϕu]*Δ1[ϕv] - Cos[ϕu] + I*m1*Sin[ϕu]*Sin[ϕv])/d];
Assert[Cos[ϕP]^2 == m1*(Cos[ϕu] + Cos[ϕv] - I*Sin[ϕu]*Sin[ϕv])/d];
Assert[Δ[ϕP]^2 == m1*(Cos[ϕv] + Δ[ϕu]*Δ1[ϕv] - I*m*Sin[ϕu]*Sin[ϕv])/d];
Assert[Abs@Sin[ϕP]^4 == (Δ[ϕu]*Δ1[ϕv] - m*Cos[ϕu] - m1*Cos[ϕv])/(m*d)];
Assert[Abs@Cos[ϕP]^4 == m1*(Δ[ϕu]*Δ1[ϕv] + m*Cos[ϕu] - m1*Cos[ϕv])/(m*d)];
Assert[Abs@Δ[ϕP]^4 == m1*(Δ[ϕu]*Δ1[ϕv] + m*Cos[ϕu] + m1*Cos[ϕv])/d];
Assert[EllipticE[ϕP, m] == (1/2)*(EllipticE[ϕu, m] + I*(v - EllipticE[ϕv, m1]) +
(m*Sin[ϕu]*(Δ[ϕu] - Cos[ϕu]*Δ1[ϕv]) + I*m1*Sin[ϕv]*(Δ[ϕu]*Cos[ϕv] + Δ1[ϕv]))/d)];
Clear[d, ϕv, ϕP, ϕM];
(* Quadratic transformations *)
m1 = 4*k/(1 + k)^2;
v = (1 + k)*u;
d = 1 + k*Sin[ϕu]^2;
Assert[JacobiSN[v, m1] == (1 + k)*Sin[ϕu]/d];
Assert[JacobiCN[v, m1] == Cos[ϕu]*Δ[ϕu]/d];
Assert[JacobiDN[v, m1] == (1 - k*Sin[ϕu]^2)/d];
Clear[d];
Assert[JacobiSN[v/2, m1] == Sqrt[1 + k]*Sin[ϕu/2]*Sqrt[(1 - m - m*Cos[ϕu] + Δ[ϕu])/((1 + Δ[ϕu])*(Δ[ϕu] - k*Cos[ϕu]))]];
Assert[JacobiCN[v/2, m1] == Sqrt[1 - k]*Cos[ϕu/2]*Sqrt[(1 + Δ[ϕu])*(Δ[ϕu] + k*Cos[ϕu])/(1 - m - m*Cos[ϕu] + Δ[ϕu])]];
Assert[JacobiDN[v/2, m1] == (Δ[ϕu] + k*Cos[ϕu])/(1 + k)];
m1 = ((1 - kp)/(1 + kp))^2;
v = (1 + kp)*u;
Assert[JacobiSN[v, m1] == (1 + kp)*Sin[ϕu]*Cos[ϕu]/Δ[ϕu]];
Assert[JacobiCN[v, m1] == (1 - (1 + kp)*Sin[ϕu]^2)/Δ[ϕu]];
Assert[JacobiDN[v, m1] == (1 - (1 - kp)*Sin[ϕu]^2)/Δ[ϕu]];
Assert[JacobiSN[v/2, m1] == (1 + kp)*Sin[ϕu]/(1 + Δ[ϕu])];
Assert[JacobiCN[v/2, m1] == Cos[ϕu]*Sqrt[2]*Sqrt[1 + kp]/(Sqrt[1 + Δ[ϕu]]*Sqrt[kp + Δ[ϕu]])];
Assert[JacobiDN[v/2, m1] == Sqrt[2]*Sqrt[kp + Δ[ϕu]]/(Sqrt[1 + kp]*Sqrt[1 + Δ[ϕu]])];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Special values of Jacobi amplitude:
Module[{m, f, s, t, a, b, pr},
pr = 50;
m = RandomReal[{0, 1}, WorkingPrecision -> pr];
f[a_, b_, m_] := With[{z = a*EllipticK[m] + I*b*EllipticK[1 - m]},
JacobiAmplitude[z, m]
];
s = m^(1/4); (* sqrt(k) *)
t = (1 - m)^(1/4); (* sqrt(k') *)
(* Real arguments *)
a = ArcSin[(1 - t)/Sqrt[1 + t^2]];
b = ArcCsc[(1 + t)/Sqrt[1 + t^2]];
Assert[f[0 , 0, m] == 0];
Assert[f[1/4, 0, m] == (1/2)*(a + b)];
Assert[f[1/2, 0, m] == ArcCot[t]];
Assert[f[3/4, 0, m] == Pi/2 + (1/2)*(a - b)];
Assert[f[1 , 0, m] == Pi/2];
(* Complex arguments *)
a = ArcCosh[((1 + s)*Sqrt[1 + s^2] + 1)/s^2];
b = ArcCosh[((1 + s)*Sqrt[1 + s^2] - 1)/s^2];
Assert[f[0, 1/4, m] == I*(1/2)*(a - b)];
Assert[f[0, 1/2, m] == I*ArcCsch[s]];
Assert[f[0, 3/4, m] == I*(1/2)*(a + b)];
a = ArcCosh[(1 + (1 - s)*Sqrt[1 + s^2])/s^2];
b = ArcCosh[(1 - (1 - s)*Sqrt[1 + s^2])/s^2];
Assert[f[1, 1/4, m] == Pi/2 + I*(1/2)*(a - b)];
Assert[f[1, 1/2, m] == Pi/2 + I*ArcSech[s]];
Assert[f[1, 3/4, m] == Pi/2 + I*(1/2)*(a + b)];
a = ArcSinh[Sqrt[2*t]/(1 - t)];
b = ArcTanh[Sqrt[2*t]/(1 + t)];
Assert[f[1/4, 1, m] == Pi/2 + I*(1/2)*(a + b)];
Assert[f[1/2, 1, m] == Pi/2 + I*ArcTanh[t]];
Assert[f[3/4, 1, m] == Pi/2 + I*(1/2)*(a - b)];
Assert[f[1 , 1, m] == Pi/2 + I*ArcTanh[t^2]];
Assert[f[1/2, 1/2, m] == ArcTan[Sqrt[1 + s^2]/t] + I*ArcTanh[t/Sqrt[1 + s^2]]];
(* Symmetries *)
{a, b} = RandomReal[{0, 1}, 2, WorkingPrecision -> pr];
Assert[f[ a, -b, m] == Conjugate[f[a, b, m]]];
Assert[f[-a, b, m] == -Conjugate[f[a, b, m]]];
Assert[f[-a, -b, m] == -f[a, b, m]];
Assert[f[a + 2, b, m] == f[a, b, m] + Pi];
Assert[f[a - 2, b, m] == f[a, b, m] - Pi];
Assert[f[a, b + 2, m] == f[a, b - 2, m] == Pi*Mod[1, 2, a - 1] - f[a, b, m]];
Assert[f[a, b + 4, m] == f[a, b - 4, m] == f[a, b, m]];
Clear[f];
(* Degenerate cases *)
a = RandomReal[{0, Pi/2}, WorkingPrecision -> pr];
b = -Log@RandomReal[{0, 1}, WorkingPrecision -> pr];
Assert[JacobiAmplitude[a, 0] == a];
Assert[JacobiAmplitude[b, 1] == Gudermannian[b]];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Integration of elliptic function:
Module[{m, m1, u, a, n, f},
m = RandomReal[WorkingPrecision -> 50];
m1 = 1 - m;
{u, a} = EllipticK[m]*RandomReal[{0, 1}, 2, WorkingPrecision -> 50];
n = RandomInteger[{3, 10}];
If[Context[JacobiEpsilon] =!= "System`",
JacobiEpsilon = (EllipticE[JacobiAmplitude[#1, #2], #2]*JacobiDN[#1, #2]/Sqrt[1 - #2*JacobiSN[#1, #2]^2] &);
];
Assert[JacobiSN[u, m] == (-Log[JacobiDN[#, m] + Sqrt[m]*JacobiCN[#, m]]/Sqrt[m] &)'[u]];
Assert[JacobiCN[u, m] == (ArcTan[Sqrt[m]*JacobiSN[#, m]/JacobiDN[#, m]]/Sqrt[m] &)'[u]];
Assert[JacobiDN[u, m] == (ArcTan[JacobiSN[#, m]/JacobiCN[#, m]] &)'[u]];
Assert[JacobiCD[u, m] == (Log[JacobiND[#, m] + Sqrt[m]*JacobiSD[#, m]]/Sqrt[m] &)'[u]];
Assert[JacobiSD[u, m] == (-ArcTan[Sqrt[m]*JacobiCD[#, m]/(Sqrt[m1]*JacobiND[#, m])]/(Sqrt[m]*Sqrt[m1]) &)'[u]];
Assert[JacobiND[u, m] == (ArcTan[Sqrt[m1]*JacobiSD[#, m]/JacobiCD[#, m]]/(Sqrt[m1]) &)'[u]];
Assert[JacobiDC[u, m] == (Log[JacobiNC[#, m] + JacobiSC[#, m]] &)'[u]];
Assert[JacobiNC[u, m] == (Log[JacobiDC[#, m] + Sqrt[m1]*JacobiSC[#, m]]/Sqrt[m1] &)'[u]];
Assert[JacobiSC[u, m] == (Log[JacobiDC[#, m] + Sqrt[m1]*JacobiNC[#, m]]/Sqrt[m1] &)'[u]];
Assert[JacobiNS[u, m] == (-Log[JacobiDS[#, m] + JacobiCS[#, m]] &)'[u]];
Assert[JacobiDS[u, m] == (-Log[JacobiNS[#, m] + JacobiCS[#, m]] &)'[u]];
Assert[JacobiCS[u, m] == (-Log[JacobiNS[#, m] + JacobiDS[#, m]] &)'[u]];
Assert[JacobiSN[u, m]^2 == ((-JacobiEpsilon[#, m] + #)/m &)'[u]];
Assert[JacobiCN[u, m]^2 == ((JacobiEpsilon[#, m] - m1*#)/m &)'[u]];
Assert[JacobiDN[u, m]^2 == (JacobiEpsilon[#, m] &)'[u]];
Assert[JacobiCD[u, m]^2 == ((-JacobiEpsilon[#, m] + # + m*JacobiSN[#, m]*JacobiCD[#, m])/m &)'[u]];
Assert[JacobiSD[u, m]^2 == ((JacobiEpsilon[#, m] - m1*# - m*JacobiSN[#, m]*JacobiCD[#, m])/(m*m1) &)'[u]];
Assert[JacobiND[u, m]^2 == ((JacobiEpsilon[#, m] - m*JacobiSN[#, m]*JacobiCD[#, m])/m1 &)'[u]];
Assert[JacobiDC[u, m]^2 == (-JacobiEpsilon[#, m] + # + JacobiSN[#, m]*JacobiDC[#, m] &)'[u]];
Assert[JacobiNC[u, m]^2 == ((-JacobiEpsilon[#, m] + m1*# + JacobiSN[#, m]*JacobiDC[#, m])/m1 &)'[u]];
Assert[JacobiSC[u, m]^2 == ((-JacobiEpsilon[#, m] + JacobiSN[#, m]*JacobiDC[#, m])/m1 &)'[u]];
Assert[JacobiNS[u, m]^2 == (-JacobiEpsilon[#, m] + # - JacobiCN[#, m]*JacobiDS[#, m] &)'[u]];
Assert[JacobiDS[u, m]^2 == (-JacobiEpsilon[#, m] + m1*# - JacobiCN[#, m]*JacobiDS[#, m] &)'[u]];
Assert[JacobiCS[u, m]^2 == (-JacobiEpsilon[#, m] - JacobiCN[#, m]*JacobiDS[#, m] &)'[u]];
Assert[JacobiSN[2*u, m] == (ArcTanh[Sqrt[m]*JacobiSN[#, m]^2]/Sqrt[m] &)'[u]];
Assert[JacobiCN[2*u, m] == (ArcTan[Sqrt[m]*JacobiCN[#, m]*JacobiSD[#, m]]/Sqrt[m] &)'[u]];
Assert[JacobiDN[2*u, m] == (ArcTan[JacobiDN[#, m]*JacobiSC[#, m]] &)'[u]];
Assert[JacobiCD[2*u, m] == (ArcTanh[Sqrt[m]*JacobiCD[#, m]*JacobiSN[#, m]]/Sqrt[m] &)'[u]];
Assert[JacobiSD[2*u, m] == (ArcTan[Sqrt[m]*Sqrt[m1]*JacobiSD[#, m]^2]/(Sqrt[m]*Sqrt[m1]) &)'[u]];
Assert[JacobiND[2*u, m] == (ArcTan[Sqrt[m1]*JacobiND[#, m]*JacobiSC[#, m]]/Sqrt[m1] &)'[u]];
Assert[JacobiDC[2*u, m] == (ArcTanh[JacobiDC[#, m]*JacobiSN[#, m]] &)'[u]];
Assert[JacobiNC[2*u, m] == (ArcTanh[Sqrt[m1]*JacobiNC[#, m]*JacobiSD[#, m]]/Sqrt[m1] &)'[u]];
Assert[JacobiSC[2*u, m] == (ArcTanh[Sqrt[m1]*JacobiSC[#, m]^2]/Sqrt[m1] &)'[u]];
Assert[JacobiNS[2*u, m] == (Log[JacobiNC[#, m]*JacobiSD[#, m]]/2 &)'[u]];
Assert[JacobiDS[2*u, m] == (Log[JacobiSN[#, m]*JacobiDC[#, m]]/2 &)'[u]];
Assert[JacobiCS[2*u, m] == (Log[JacobiSN[#, m]*JacobiCD[#, m]]/2 &)'[u]];
f[u_] := JacobiSN[u, m]^(n - 3)*JacobiCN[u, m]*JacobiDN[u, m];
f[u_, n_] := JacobiSN[u, m]^n;
Assert[f[u, n] == 1/(m*(n - 1))*f'[u] + (1 + m)*(n - 2)/(m*(n - 1))*f[u, n - 2] - (n - 3)/(m*(n - 1))*f[u, n - 4]];
Clear[f];
f[u_] := JacobiCN[u, m]^(n - 3)*JacobiSN[u, m]*JacobiDN[u, m];
f[u_, n_] := JacobiCN[u, m]^n;
Assert[f[u, n] == 1/(m*(n - 1))*f'[u] + (m - m1)*(n - 2)/(m*(n - 1))*f[u, n - 2] + m1*(n - 3)/(m*(n - 1))*f[u, n - 4]];
Clear[f];
f[u_] := JacobiDN[u, m]^(n - 3)*JacobiSN[u, m]*JacobiCN[u, m];
f[u_, n_] := JacobiDN[u, m]^n;
Assert[f[u, n] == m/(n - 1)*f'[u] + (1 + m1)*(n - 2)/(n - 1)*f[u, n - 2] - m1*(n - 3)/(n - 1)*f[u, n - 4]];
Clear[f];
f[u_] := JacobiCD[u, m]^(n - 3)*JacobiSD[u, m]*JacobiND[u, m];
f[u_, n_] := JacobiCD[u, m]^n;
Assert[f[u, n] == -m1/(m*(n - 1))*f'[u] + (1 + m)*(n - 2)/(m*(n - 1))*f[u, n - 2] - (n - 3)/(m*(n - 1))*f[u, n - 4]];
Clear[f];
f[u_] := JacobiSD[u, m]^(n - 3)*JacobiCD[u, m]*JacobiND[u, m];
f[u_, n_] := JacobiSD[u, m]^n;
Assert[f[u, n] == -1/(m*m1*(n - 1))*f'[u] + (m - m1)*(n-2)/(m*m1*(n-1))*f[u, n - 2] + (n - 3)/(m*m1*(n - 1))*f[u, n - 4]];
Clear[f];
f[u_] := JacobiND[u, m]^(n - 3)*JacobiCD[u, m]*JacobiSD[u, m];
f[u_, n_] := JacobiND[u, m]^n;
Assert[f[u, n] == -m/(m1*(n - 1))*f'[u] + (1 + m1)*(n - 2)/(m1*(n - 1))*f[u, n - 2] - (n - 3)/(m1*(n - 1))*f[u, n - 4]];
Clear[f];
f[u_] := JacobiDC[u, m]^(n - 3)*JacobiNC[u, m]*JacobiSC[u, m];
f[u_, n_] := JacobiDC[u, m]^n;
Assert[f[u, n] == m1/(n - 1)*f'[u] + (1 + m)*(n - 2)/(n - 1)*f[u, n - 2] - m*(n - 3)/(n - 1)*f[u, n - 4]];
Clear[f];
f[u_] := JacobiNC[u, m]^(n - 3)*JacobiDC[u, m]*JacobiSC[u, m];
f[u_, n_] := JacobiNC[u, m]^n;
Assert[f[u, n] == 1/(m1*(n - 1))*f'[u] + (m1 - m)*(n - 2)/(m1*(n - 1))*f[u, n - 2] + m*(n - 3)/(m1*(n - 1))*f[u, n - 4]];
Clear[f];
f[u_] := JacobiSC[u, m]^(n - 3)*JacobiDC[u, m]*JacobiNC[u, m];
f[u_, n_] := JacobiSC[u, m]^n;
Assert[f[u, n] == 1/(m1*(n - 1))*f'[u] - (1 + m1)*(n - 2)/(m1*(n - 1))*f[u, n - 2] - (n - 3)/(m1*(n - 1))*f[u, n - 4]];
Clear[f];
f[u_] := JacobiNS[u, m]^(n - 3)*JacobiDS[u, m]*JacobiCS[u, m];
f[u_, n_] := JacobiNS[u, m]^n;
Assert[f[u, n] == -1/(n - 1)*f'[u] + (1 + m)*(n - 2)/(n - 1)*f[u, n - 2] - m*(n - 3)/(n - 1)*f[u, n - 4]];
Clear[f];
f[u_] := JacobiDS[u, m]^(n - 3)*JacobiNS[u, m]*JacobiCS[u, m];
f[u_, n_] := JacobiDS[u, m]^n;
Assert[f[u, n] == -1/(n - 1)*f'[u] + (m1 - m)*(n - 2)/(n - 1)*f[u, n - 2] + m*m1*(n - 3)/(n - 1)*f[u, n - 4]];
Clear[f];
f[u_] := JacobiCS[u, m]^(n - 3)*JacobiNS[u, m]*JacobiDS[u, m];
f[u_, n_] := JacobiCS[u, m]^n;
Assert[f[u, n] == -1/(n - 1)*f'[u] - (1 + m1)*(n - 2)/(n - 1)*f[u, n - 2] - m1*(n - 3)/(n - 1)*f[u, n - 4]];
Clear[f];
f[u_] := Cos[JacobiAmplitude[u, m]]*JacobiDN[u, m];
f[u_, n_] := Sin[JacobiAmplitude[u, m]]
Assert[f[u, n] == -4/(m*(n - 1))*f'[u] - 2*(1 + m1)*(n - 2)/(m*(n - 1))*f[u, n - 2] - (n - 3)/(n - 1)*f[u, n - 4]];
Clear[f];
f[u_] := Sin[JacobiAmplitude[u, m]]*JacobiDN[u, m];
f[u_, n_] := Cos[JacobiAmplitude[u, m]]
Assert[f[u, n] == 4/(m*(n - 1))*f'[u] - 2*(1 + m1)*(n - 2)/(m*(n - 1))*f[u, n - 2] - (n - 3)/(n - 1)*f[u, n - 4]];
Clear[f];
f[u_?InexactNumberQ] := -With[{k = Sqrt[m]}, Sqrt[2/k]*InverseJacobiSD[Sqrt[2*k/m1]*JacobiCN[u/2, m], (1 + k)/2]];
Assert[f'[u] == Sin[(1/2)*JacobiAmplitude[u, m]]];
Clear[f];
f[u_?InexactNumberQ] := With[{k = Sqrt[m]}, Sqrt[2/k]*InverseJacobiSD[Sqrt[2*k]*JacobiSD[u/2, m], (1 + k)/2]];
Assert[f'[u] == Cos[(1/2)*JacobiAmplitude[u, m]]];
Clear[f];
(* EllipticPi *)
f[u_] := EllipticPi[m*JacobiSN[a, m]^2, JacobiAmplitude[u, m], m];
Assert[f'[u] == 1/(1 - m*JacobiSN[u, m]^2*JacobiSN[a, m]^2)];
Clear[f];
f[u_] := With[{k = Sqrt[m]},
+ EllipticPi[m*JacobiSN[a, m]^2, JacobiAmplitude[u, m], m]
+ (1/2)*JacobiNC[a, m]*JacobiSD[a, m]*Log[(JacobiDN[a + u, m] + k*JacobiCN[a + u, m])/(JacobiDN[a - u, m] - k*JacobiCN[a - u, m])]
];
Assert[f'[u] == 1/(1 + Sqrt[m]*JacobiSN[u, m]*JacobiSN[a, m])];
Clear[f];
f[u_] := With[{k = Sqrt[m]},
+ EllipticPi[m*JacobiSN[a, m]^2, JacobiAmplitude[u, m], m]
- (1/2)*JacobiNC[a, m]*JacobiSD[a, m]*Log[(JacobiDN[a + u, m] + k*JacobiCN[a + u, m])/(JacobiDN[a - u, m] - k*JacobiCN[a - u, m])]
];
Assert[f'[u] == 1/(1 - Sqrt[m]*JacobiSN[u, m]*JacobiSN[a, m])];
Clear[f];
f[u_] := (
+ (1/2)*(JacobiSN[a, m]*JacobiCN[a, m]*JacobiDN[a, m])^(-1)*Log[JacobiSN[a - u, m]/JacobiSN[a + u, m]]
+ JacobiSN[a, m]^(-2)*(EllipticPi[m*JacobiSN[a, m]^2, JacobiAmplitude[u, m], m] - u)
);
Assert[f'[u] == 1/(JacobiSN[u, m]^2 - JacobiSN[a, m]^2)];
Clear[f];
f[u_] := (
- (1/2)*(JacobiCN[a, m]*JacobiDN[a, m])^(-1)*Log[(JacobiDN[a - u, m] + JacobiCN[a - u, m])/(JacobiDN[a + u, m] - JacobiCN[a + u, m])]
- JacobiSN[a, m]^(-1)*(EllipticPi[m*JacobiSN[a, m]^2, JacobiAmplitude[u, m], m] - u)
);
Assert[f'[u] == 1/(JacobiSN[u, m] + JacobiSN[a, m])];
Clear[f];
f[u_] := (
+ (1/2)*(JacobiCN[a, m]*JacobiDN[a, m])^(-1)*Log[(JacobiDN[a - u, m] - JacobiCN[a - u, m])/(JacobiDN[a + u, m] + JacobiCN[a + u, m])]
+ JacobiSN[a, m]^(-1)*(EllipticPi[m*JacobiSN[a, m]^2, JacobiAmplitude[u, m], m] - u)
);
Assert[f'[u] == 1/(JacobiSN[u, m] - JacobiSN[a, m])];
Clear[f];
f[u_] := (
+ (2*JacobiSN[a, m]*JacobiDN[a, m])^(-1)*Log[(1 + JacobiCN[a - u, m])/(1 + JacobiCN[a + u, m])]
+ (JacobiCN[a, m]/JacobiSN[a, m]^2)*(EllipticPi[m*JacobiSN[a, m]^2, JacobiAmplitude[u, m], m] - u)
);
Assert[f'[u] == 1/(JacobiCN[u, m] + JacobiCN[a, m])];
Clear[f];
f[u_] := (
- (2*JacobiSN[a, m]*JacobiDN[a, m])^(-1)*Log[(1 - JacobiCN[a - u, m])/(1 - JacobiCN[a + u, m])]
- (JacobiCN[a, m]/JacobiSN[a, m]^2)*(EllipticPi[m*JacobiSN[a, m]^2, JacobiAmplitude[u, m], m] - u)
);
Assert[f'[u] == 1/(JacobiCN[u, m] - JacobiCN[a, m])];
Clear[f];
f[u_] := (
+ (2*m*JacobiSN[a, m]*JacobiCN[a, m])^(-1)*Log[(1 + JacobiDN[a - u, m])/(1 + JacobiDN[a + u, m])]
+ (JacobiDN[a, m]/(m*JacobiSN[a, m]^2))*(EllipticPi[m*JacobiSN[a, m]^2, JacobiAmplitude[u, m], m] - u)
);
Assert[f'[u] == 1/(JacobiDN[u, m] + JacobiDN[a, m])];
Clear[f];
f[u_] := (
- (2*m*JacobiSN[a, m]*JacobiCN[a, m])^(-1)*Log[(1 - JacobiDN[a - u, m])/(1 - JacobiDN[a + u, m])]
- (JacobiDN[a, m]/(m*JacobiSN[a, m]^2))*(EllipticPi[m*JacobiSN[a, m]^2, JacobiAmplitude[u, m], m] - u)
);
Assert[f'[u] == 1/(JacobiDN[u, m] - JacobiDN[a, m])];
Clear[f];
(* Special cases for EllipticPi *)
f[u_] := With[{k = Sqrt[m]}, (1/2)*(u + (1 + k)^(-1)*ArcTan[(1 + k)*JacobiNC[u, m]*JacobiSD[u, m]])];
Assert[f'[u] == (1 + Sqrt[m]*JacobiSN[u, m]^2)^(-1)];
Clear[f];
f[u_] := With[{k = Sqrt[m]}, (1/2)*(u + (1 - k)^(-1)*ArcTan[(1 - k)*JacobiNC[u, m]*JacobiSD[u, m]])];
Assert[f'[u] == (1 - Sqrt[m]*JacobiSN[u, m]^2)^(-1)];
Clear[f];
f[u_] := With[{kp = Sqrt[m1]}, (2*kp)^(-1)*((kp + 1)*u + Log[JacobiDN[u + (1/2)*EllipticK[m], m]])];
Assert[f'[u] == (1 - (1 - Sqrt[m1])*JacobiSN[u, m]^2)^(-1)];
Clear[f];
f[u_] := With[{kp = Sqrt[m1]}, (2*kp)^(-1)*((kp - 1)*u + Log[JacobiSC[u + (1/2)*EllipticK[m], m]])];
Assert[f'[u] == (1 - (1 + Sqrt[m1])*JacobiSN[u, m]^2)^(-1)];
Clear[f];
Assert[EllipticPi[-Sqrt[m], m] == (1/4)*Pi*(1 + Sqrt[m])^(-1) + (1/2)*EllipticK[m]];
Assert[EllipticPi[+Sqrt[m], m] == (1/4)*Pi*(1 - Sqrt[m])^(-1) + (1/2)*EllipticK[m]];
Assert[EllipticPi[1 - Sqrt[1 - m], m] == (1/2)*(1 + (1 - m)^(-1/2))*EllipticK[m]];
Assert[EllipticPi[1 + Sqrt[1 - m], m] == (1/2)*(1 - (1 - m)^(-1/2))*EllipticK[m]]; (* Cauchy principal value *)
] & // Block[{$AssertFunction = Automatic}, #[]] &;Integration of zeta squared:
If[Context[JacobiZN] =!= "System`",
JacobiZN[u_, m_] := JacobiZeta[JacobiAmplitude[u, m], m];
];
f[u_, m_, n0_] :=
Module[{pr, q, z, m1, t, n, res = $Failed},
If[
VectorQ[{u, m}, NumericQ] && (pr = Accuracy[{u, m}]) < Infinity,
pr += 5;
m1 = SetAccuracy[m, pr];
q = EllipticNomeQ[m1];
z = Exp[I*Pi*SetAccuracy[u, pr]/EllipticK[m1]];
res = 1;
n = n0;
While[
t = res;
res *= ((1 - q^n*z)/(1 - q^n/z))^n;
t != res,
n += 2;
];
res = I*Log[res];
res = SetAccuracy[res, Accuracy[res] - 5];
];
res /; NumberQ[res]
];
On[Assert];
m = RandomReal[{1/3, 2/3}, WorkingPrecision -> 50];
u = RandomReal[{0, EllipticK[m]}, WorkingPrecision -> 50];
Assert[
Derivative[1, 0, 0][f][u, m, 1] ==
(EllipticK[m]/Pi*(JacobiZN[u, m]^2 - m*JacobiSN[u, m]^2 + (1/3)*(1 + m)) - Pi/(12*EllipticK[m]))
];
Assert[
Derivative[1, 0, 0][f][u, m, 2] ==
(EllipticK[m]/Pi*(JacobiZN[u, m]^2 + 2*JacobiZN[u, m]*JacobiCN[u, m]*JacobiDS[u, m] + m*JacobiSN[u, m]^2 - (2/3)*(1 + m)) + Pi/(6*EllipticK[m]))
];Repeated integration:
On[Assert];
g[z_] := HypergeometricPFQ[{1, 1, 1}, {3/2, 3/2}, z];
f[z_] :=
Module[{n = 1, a = 1, b = 1, c, t, s = 0, r},
c = 1 - SetPrecision[z, Precision[z] + 5];
t = c/4;
While[
r = s;
s += a*t;
r != s,
b *= n*(n + 1);
n += 2;
a = a*n^2 + b;
t *= c/(n + 1)^2;
];
SetPrecision[s, Precision[s] - 5]
];
m = RandomReal[{1/3, 2/3}, WorkingPrecision -> 50];
Assert[f[m] == Sqrt[m]*HypergeometricPFQ[{1, 1, 1}, {3/2, 3/2}, m] - (Pi/2)*EllipticK[m] + (4*Catalan/Pi)*EllipticK[1 - m]];
With[{z = SetPrecision[EllipticK[m], Infinity], m1 = SetPrecision[m, Infinity], kp = Sqrt[1 - m]},
(* ∫_0^K d(u') ∫_0^u' d(u) sn(u) *)
Assert[NIntegrate[(z - t)*JacobiSN[t, m1], {t, 0, z}, WorkingPrecision -> 50] == m^(-1/2)*ArcTanh[m^(1/2)]*z - (1 - m)^(-1)*g[m/(m - 1)]];
(* ∫_0^K d(u') ∫_0^u' d(u) cn(u) *)
Assert[NIntegrate[(z - t)*JacobiCN[t, m1], {t, 0, z}, WorkingPrecision -> 50] == g[m]];
(* ∫_0^K d(u') ∫_0^u' d(u) dn(u) *)
Assert[NIntegrate[(z - t)*JacobiDN[t, m1], {t, 0, z}, WorkingPrecision -> 50] == (1/4)*Pi*z + m/(1 + kp)^3*g[((1 - kp)/(1 + kp))^2]];
];
(* Functional equations *)
m += I*RandomReal[{-1/3, 1/3}, WorkingPrecision -> 50];
Assert[
+ Sqrt[m]*g[m] ==
+ (Sqrt[m] + I*Sqrt[1 - m])^3*g[(Sqrt[m] + I*Sqrt[1 - m])^4]
+ I*Sqrt[1 - m]*g[1 - m]
+ (Pi/4)*(EllipticK[m] - I*EllipticK[1 - m])
];
Assert[
+ Sqrt[m]*g[m] ==
+ (Sqrt[m] - I*Sqrt[1 - m])^3*g[(Sqrt[m] - I*Sqrt[1 - m])^4]
- I*Sqrt[1 - m]*g[1 - m]
+ (Pi/4)*(EllipticK[m] + I*EllipticK[1 - m])
];
Assert[
+ 2*Sqrt[m]*g[m] ==
+ (Sqrt[m] + I*Sqrt[1 - m])^3*g[(Sqrt[m] + I*Sqrt[1 - m])^4]
+ (Sqrt[m] - I*Sqrt[1 - m])^3*g[(Sqrt[m] - I*Sqrt[1 - m])^4]
+ (Pi/2)*EllipticK[m]
];
Assert[
+ 2*m*g[m] ==
+ (1 + Sqrt[m])*g[ (1 + Sqrt[m])^2/(4*Sqrt[m])]
+ (1 - Sqrt[m])*g[-(1 - Sqrt[m])^2/(4*Sqrt[m])]
+ Pi*Sqrt[m]*(-EllipticK[m] + I*Sign@Im[m]*EllipticK[1 - m])
];
Assert[
+ m*g[m] ==
+ (1 + Sqrt[m])*g[(1 + Sqrt[m])^2/(4*Sqrt[m])]
- m^(-1/2)*g[1/m]
+ Pi*Sqrt[m]*(-EllipticK[m] + I*(1/2)*Sign@Im[m]*EllipticK[1 - m])
];
Assert[
+ m*g[m] ==
+ (1 - Sqrt[m])*g[-(1 - Sqrt[m])^2/(4*Sqrt[m])]
+ m^(-1/2)*g[1/m]
+ I*(Pi/2)*Sign@Im[m]*Sqrt[m]*EllipticK[1 - m]
];References:
Transformations of the Jacobian amplitude function
Jacobi elliptic functions identities
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}, #[]] &;Duplication theorem:
c[x_, y_] := CarlsonRC[x, y];
d[x_, y_, z_] := CarlsonRD[x, y, z];
f[x_, y_, z_] := CarlsonRF[x, y, z];
g[x_, y_, z_] := CarlsonRG[x, y, z] - (1/6)*(x + y + z)*CarlsonRF[x, y, z];
j[x_, y_, z_, p_] := CarlsonRJ[x, y, z, p];
i[x_, y_, z_] := With[{a = (x + y + z)/3},
WeierstrassSigma[CarlsonRF[x, y, z], {(2/3)*((x - y)^2 + (x - z)^2 + (y - z)^2), 4*(a - x)*(a - y)*(a - z)}]
];
{x, y, z, p} = RandomReal[{1, 2}, 4, WorkingPrecision -> 50];
λ = Sqrt[x]*Sqrt[y] + Sqrt[x]*Sqrt[z] + Sqrt[y]*Sqrt[z];
On[Assert];
Assert[c[x, y] == c[(Sqrt[x] + Sqrt[y])^2/4, 2*(y + Sqrt[x]*Sqrt[y])/4]];
Assert[d[x, y, z] == (1/4)*d[(x + λ)/4, (y + λ)/4, (z + λ)/4] + 3/(Sqrt[z]*(z + λ))];
Assert[f[x, y, z] == f[(x + λ)/4, (y + λ)/4, (z + λ)/4]];
Assert[g[x, y, z] == 4*g[(x + λ)/4, (y + λ)/4, (z + λ)/4] - (1/2)*(Sqrt[x] + Sqrt[y] + Sqrt[z])];
Assert[i[x, y, z] == (1/8)*i[(x + λ)/4, (y + λ)/4, (z + λ)/4]^4*(Sqrt[x] + Sqrt[y])*(Sqrt[x] + Sqrt[z])*(Sqrt[y] + Sqrt[z])];
α = p*(Sqrt[x] + Sqrt[y] + Sqrt[z]) + Sqrt[x]*Sqrt[y]*Sqrt[z];
β = Sqrt[p]*(p + λ);
Assert[j[x, y, z, p] == (1/4)*j[(x + λ)/4, (y + λ)/4, (z + λ)/4, (p + λ)/4] + 3*CarlsonRC[α^2, β^2]];
δ = (p - x)*(p - y)*(p - z); (* β^2 - α^2 *)
a = Sqrt@Abs[δ]; (* sqrt(|β^2 - α^2|) *)
b = (Sqrt[p] + Sqrt[x])*(Sqrt[p] + Sqrt[y])*(Sqrt[p] + Sqrt[z]); (* β + α *)
(* Circular case *)
If[δ > 0, Assert[CarlsonRC[α^2, β^2] == 2*ArcTan[b, a]/a]];
(* Hyperbolic case *)
If[δ < 0, Assert[CarlsonRC[α^2, β^2] == Log[1 + a*(1 + a/b)/β]/a]];
(* Degenerate case *)
If[δ == 0, Assert[CarlsonRC[α^2, β^2] == 1/α == 1/β]];Rhombic -> Rectangular:
Module[{a0, a1, a2, x, y, z, pr, rc, rd, rf, rg, rj},
{rc, rd, rf, rg, rj} = Function[{f}, Module[{args, res},
args = {##};
res = f @@ SetAccuracy[args, Accuracy@args + 5];
SetAccuracy[res, Accuracy@res - 5]
] &] /@ {CarlsonRC, CarlsonRD, CarlsonRF, CarlsonRG, CarlsonRJ};
pr = 50;
a0 = RandomReal[{0, 1}, WorkingPrecision -> pr];
a1 = RandomComplex[{0, 1 + I}, WorkingPrecision -> pr];
a2 = Conjugate[a1];
x = (1/2)*(Abs[a1] - Abs[a0 - a1] + a0);
y = (1/2)*(Abs[a1] + Re[a1]);
z = (1/2)*(Abs[a1] + Abs[a0 - a1] + a0);
Assert[x <= y <= z];
Assert[rf[a0, a1, a2] == rf[x, y, z]];
] & // Block[{$AssertFunction = Automatic}, #[]] &;Relations to elliptic divisibility sequence:
Module[{z, w1, w2, g2, g3, a, a2, a3, a4, b, σ},
w1 = 1;
w2 = I*Sqrt[Tan[RandomReal[WorkingPrecision -> 50]*Pi/2]];
{g2, g3} = WeierstrassInvariants[{w1, w2}];
z = w1*(2*RandomReal[WorkingPrecision -> 50] - 1);
σ[z__] := Times @@ WeierstrassSigma[{z}, {g2, g3}];
a[n_] := σ[n*z]/σ[z]^(n^2);
x = WeierstrassP[z, {g2, g3}];
y = WeierstrassPPrime[z, {g2, g3}];
{a2, a3, a4} = a /@ {2, 3, 4};
b = (a2^5 + a4)/(a2*a3); (* ℘''(z) *)
Assert[x == (b^2/4 + a3)/(3*a2^2)];
Assert[y == -a2];
Assert[g2 == 12*x^2 - 2*b];
Assert[g3 == 4*x^3 - g2*x - a2^2];
] & // 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}];
(* from https://dlmf.nist.gov/23.10#E8 *)
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[I*Sqrt[3]*(q^(1/9)*EllipticTheta[1, I*Log[q]/3, q]/EllipticTheta[1, Pi/3, q])^3 == 1 - a/(a^3 - 1)^(1/3)];
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
}
BranchPoints[InverseDixonSM[z^(1/3), α], z] == {
0,
(Sqrt[1 + α^3] + α^(3/2))^2,
(Sqrt[1 + α^3] - α^(3/2))^2,
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}, #[]] &;Specific values for equianharmonic (α=0) case: