Math Is Fun Forum

BeginPackage["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["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[];

(* https://resources.wolframcloud.com/PacletRepository/resources/JanM/Dixon/ *)
<< "JanM`Dixon`";
Hold[{f, u, v, α, s, c, x, y, z, b, γ, w1, w2, λ, λ0, λ1, λ2, κ, κ0, κ1, κ2, ϕ, 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];
	λ0 = Gamma[1/3]^3/(2*Sqrt[3]*Pi);
	κ0 = 4*Pi^2/(3*Gamma[1/3]^3);
	Assert[N[DixonF[λ0, 0] == κ0, pr]];

	x = λ0*Hypergeometric2F1[1/3, 1/3, 2/3, -α^3];
	y = κ0*α*Hypergeometric2F1[2/3, 2/3, 4/3, -α^3];
	λ = Quiet@Re[x + y];
	λ1 = w1*x + w2*y;
	λ2 = w2*x + w1*y;

	x = κ0*Hypergeometric2F1[-1/3, 2/3, 1/3, -α^3];
	y = λ0*α^2*Hypergeometric2F1[1/3, 4/3, 5/3, -α^3]/2;
	κ = Quiet@Re[x + y];
	κ1 = α + w2*x + w1*y;
	κ2 = α + w1*x + w2*y;

	If[α < 0 && !VectorQ[{λ, κ}, NumericQ], With[{b = 1 + α^(-3)},
		{λ, κ} = 2*Pi*3^(-3/2)*{
			-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)];
	With[{v = λ1 - λ2}, If[NumericQ[v],
		Assert[DixonSM[u + v, α] == w2*s];
		Assert[DixonCM[u + v, α] == w1*c];
	]];

	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];
	];

	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*α]},
		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[α > -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[α]},
		Assert[DixonSM[s*u, α] == s*E^((1/3)*ArcSinh[s^3])];
		Assert[DixonCM[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];
		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[w1*u, w2*α] == {w1*s, c}];
		Assert[f[u*(1 + α)/(w2 - w1), (2 - α)/(1 + α)] == -{s + c - 1, s + w2*c - w1}/(s + w1*c - w2)];
		Assert[f[3*u*γ/(w2 - w1), -α/γ] == -{3*γ*s*c, s^3 + w2*c^3 - w1}/(s^3 + w1*c^3 - w2)];
		Assert[f[u*γ/(w2 - w1), -α/γ]^3 == -{α*s - c + 1, α*s - w2*c + w1}/(α*s - w1*c + w2)];
		Assert[f[b*u, (α - 2)/b] == {b*s*c, -(s - c^2 + c*s^2)}/(c - s^2 + s*c^2)];
		Assert[f[b*u/3, (α - 2)/b]^3 == {α*s - c + 1, α*c - s + 1}/(s + c + α)];
		Assert[f[(γ - α)*u, (2*γ + α)/(γ - α)] == {s^3 + (γ - α)*s*c, c^3 + (γ - α)*s*c}/(1 - (γ - α)*s*c)];
	];

	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)];
		];
	];
] /.
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]
	]
];

Module[{u, q, a, b},
	u = RandomReal[WorkingPrecision -> 50];
	q = RandomReal[WorkingPrecision -> 50];
	a = ((EllipticTheta[4, q]/EllipticTheta[4, Pi/3, q])^2 + 2*EllipticTheta[4, Pi/3, q]/EllipticTheta[4, q])/3;
	b = EllipticTheta[2, q]*EllipticTheta[3, q]*EllipticTheta[4, q]/(2*Sqrt[3]*EllipticTheta[1, Pi/3, q]);
	Assert[EllipticTheta[1, u, q^3]^3 == b*(a*EllipticTheta[1, u, q] - EllipticTheta[1, u + Pi/3, q] + EllipticTheta[1, u + 2*Pi/3, q])];
	Assert[EllipticTheta[2, u, q^3]^3 == b*(a*EllipticTheta[2, u, q] - EllipticTheta[2, u + Pi/3, q] + EllipticTheta[2, u + 2*Pi/3, q])];
	Assert[EllipticTheta[3, u, q^3]^3 == b*(a*EllipticTheta[3, u, q] + EllipticTheta[3, u + Pi/3, q] + EllipticTheta[3, u + 2*Pi/3, q])];
	Assert[EllipticTheta[4, u, q^3]^3 == b*(a*EllipticTheta[4, u, q] + EllipticTheta[4, u + Pi/3, q] + EllipticTheta[4, u + 2*Pi/3, q])];
	Assert[EllipticTheta[1, u, q^3] ==  EllipticTheta[1, u/3, q]*EllipticTheta[1, (u + Pi)/3, q]*EllipticTheta[1, (u + 2*Pi)/3, q]/(3*b)];
	Assert[EllipticTheta[2, u, q^3] == -EllipticTheta[2, u/3, q]*EllipticTheta[2, (u + Pi)/3, q]*EllipticTheta[2, (u + 2*Pi)/3, q]/(3*b)];
	Assert[EllipticTheta[3, u, q^3] ==  EllipticTheta[3, u/3, q]*EllipticTheta[3, (u + Pi)/3, q]*EllipticTheta[3, (u + 2*Pi)/3, q]/(3*b)];
	Assert[EllipticTheta[4, u, q^3] ==  EllipticTheta[4, u/3, q]*EllipticTheta[4, (u + Pi)/3, q]*EllipticTheta[4, (u + 2*Pi)/3, q]/(3*b)];
	Assert[EllipticTheta[1, u, q^9] == -(EllipticTheta[1, u/3, q] - EllipticTheta[1, (u + Pi)/3, q] + EllipticTheta[1, (u + 2*Pi)/3, q])/3];
	Assert[EllipticTheta[2, u, q^9] ==  (EllipticTheta[2, u/3, q] - EllipticTheta[2, (u + Pi)/3, q] + EllipticTheta[2, (u + 2*Pi)/3, q])/3];
	Assert[EllipticTheta[3, u, q^9] ==  (EllipticTheta[3, u/3, q] + EllipticTheta[3, (u + Pi)/3, q] + EllipticTheta[3, (u + 2*Pi)/3, q])/3];
	Assert[EllipticTheta[4, u, q^9] ==  (EllipticTheta[4, u/3, q] + EllipticTheta[4, (u + Pi)/3, q] + EllipticTheta[4, (u + 2*Pi)/3, q])/3];
	With[{s = EllipticTheta[1, u, q], m = EllipticTheta[1, u + Pi/3, q], c = EllipticTheta[1, u + 2*Pi/3, q]},
		Assert[EllipticTheta[1, u+Pi/9, q]^3 == ( EllipticTheta[1, Pi/9, q]^3*c*m^2 + EllipticTheta[1, 2*Pi/9, q]^3*s*c^2 + EllipticTheta[1, 4*Pi/9, q]^3*m*s^2)/EllipticTheta[1, Pi/3, q]^3]
		Assert[EllipticTheta[1, u-Pi/9, q]^3 == (-EllipticTheta[1, Pi/9, q]^3*m*c^2 + EllipticTheta[1, 2*Pi/9, q]^3*s*m^2 - EllipticTheta[1, 4*Pi/9, q]^3*c*s^2)/EllipticTheta[1, Pi/3, q]^3];
	];
] & // Block[{$AssertFunction = Automatic}, #[]] &;

(* 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);
}