Math Is Fun Forum

BeginPackage["CubicTheta`"];

Unprotect[CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime];

ClearAll[CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime];

Begin["`Private`"];

CubicTheta[u_, v_, q_] /;
(VectorQ[{u, v, q}, NumericQ] && Abs[q] < 1 && Precision[{u, v, q}] < Infinity) :=
Module[{p = q^3, t, a = u + v, b = u - v},
	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))];
	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[{CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime}, {Listable, NumericFunction, ReadProtected}];

Protect[CubicTheta, CubicThetaA, CubicThetaB, CubicThetaC, CubicThetaAPrime, CubicThetaBPrime, CubicThetaCPrime];

EndPackage[];

Module[{a, b, c, q},
q = RandomReal[WorkingPrecision -> 20];
a = CubicThetaA[q];
b = CubicThetaB[q];
c = CubicThetaC[q];
On[Assert];
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))
]]];
]