Math Is Fun Forum

lanxiyu

Member

Registered: 2022-05-10

Posts: 66

Mathematica test for Jacobi elliptic functions

Module[{Δ, Δ1, m, k, kp, u, v, d, c, f, ϕ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];

	f[ϕ1_, ϕ2_, ϕ3_] := Module[{x, y, r, ϕ}, 
		x = 1 + Δ[ϕ1] + Δ[ϕ2] + Δ[ϕ3];
		y = m*Sin[ϕ1]*Sin[ϕ2]*Sin[ϕ3];
		r = Sqrt[2*(1 + Δ[ϕ1])*(1 + Δ[ϕ2])*(1 + Δ[ϕ3])];
		ϕ = ArcTan[x, y];
		Assert[ϕ == 2*ArcTan[x + r, y] == ArcSin[y/r]];
		2*ϕ
	];

	Assert[ϕP == ϕu + ϕv - f[ϕu, ϕv, ϕP]];
	Assert[ϕM == ϕu - ϕv + f[ϕu, ϕv, ϕM]];

	Clear[f];

	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[1 + Δ[ϕP] == ((Δ[ϕu] + Δ[ϕv])^2 + m*Sin[ϕu - ϕv]^2)/(2*d)];
	Assert[1 + Δ[ϕM] == ((Δ[ϕu] + Δ[ϕv])^2 + m*Sin[ϕu + ϕv]^2)/(2*d)];
	Assert[1 - Δ[ϕP] == ((Δ[ϕu] - Δ[ϕv])^2 + m*Sin[ϕu + ϕv]^2)/(2*d)];
	Assert[1 - Δ[ϕM] == ((Δ[ϕu] - Δ[ϕv])^2 + m*Sin[ϕu - ϕv]^2)/(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[Cos[ϕP] == Cos[ϕu]*Cos[ϕv] - Sin[ϕu]*Sin[ϕv]*Δ[ϕP]];
	Assert[Cos[ϕM] == Cos[ϕu]*Cos[ϕv] + Sin[ϕu]*Sin[ϕv]*Δ[ϕM]];
	Assert[Δ[ϕP] == Δ[ϕu]*Δ[ϕv] - m*Sin[ϕu]*Sin[ϕv]*Cos[ϕP]];
	Assert[Δ[ϕM] == Δ[ϕu]*Δ[ϕv] + m*Sin[ϕu]*Sin[ϕv]*Cos[ϕM]];

	f[x_, y_, z_, a_, b_] := a*(1 - x^2)*(1 - y^2)*(1 - z^2) == b*(2*x*y*z - x^2 - y^2 - z^2 + 1);
	Assert@f[Cos[ϕu], Cos[ϕv], Cos[ϕP], m, 1];
	Assert@f[Cos[ϕu], Cos[ϕv], Cos[ϕM], m, 1];
	Assert@f[Δ[ϕu], Δ[ϕv], Δ[ϕP], 1, m];
	Assert@f[Δ[ϕu], Δ[ϕv], Δ[ϕM], 1, m];
	Assert@f[Sec[ϕu], Sec[ϕv], Sec[ϕP], m1, 1];
	Assert@f[Sec[ϕu], Sec[ϕv], Sec[ϕM], m1, 1];
	Clear[f];

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

	f[ϕ1_, ϕ2_, ϕ3_] := Module[{x, y, r}, 
		x = m*sn[c]*cn[c]*dn[c]*Sin[ϕ1]*Sin[ϕ2]*Sin[ϕ3];
		y = dn[c]^2 + m*sn[c]^2*Cos[ϕ1]*Cos[ϕ2]*Cos[ϕ3];
		r = Sqrt[Times @@ ((1 - m*sn[c]^2*Sin[#]^2) & /@ {ϕ1, ϕ2, ϕ3})];
		d == ArcTanh[x/y] == ArcSinh[x/r] == Log[(x + y)/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, ϕP, m] == EllipticPi[m*sn[c]^2, ϕu, m] + EllipticPi[m*sn[c]^2, ϕv, m] + (tn[c]/dn[c])*d];
	Assert@f[ϕu, ϕv, ϕP];

	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];
	Assert@f[ϕu, ϕv, ϕM];

	Clear[f, d, c];

	(* Inequalities for addition formulas *)
	With[{ϕ = am[(2/3)*EllipticK[m]]},
		Assert[Abs[ϕP - (ϕu + ϕv)] <= 3*ϕ - Pi];
		Assert[Abs[ϕM - (ϕu - ϕv)] <= 3*ϕ - Pi];
		Assert[Abs[En[ϕP] - (En[ϕu] + En[ϕv])] <= m*Sin[ϕ]^3];
		Assert[Abs[En[ϕM] - (En[ϕu] - En[ϕv])] <= m*Sin[ϕ]^3];
	];

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

	f = (#1[#3/2]*#1[#4/2]*#1[(#3 + #4)/2])*(#2[#3/2]*#2[#4/2]*#2[(#3 - #4)/2]) &;
	d = (1 + Δ[ϕu])*(1 + Δ[ϕv]) - m*(1 + Cos[ϕu])*(1 + Cos[ϕv]);

	Assert[f[sn, 1&, u, v] == (+Sin[ϕu]*Δ[ϕv] + Sin[ϕv]*Δ[ϕu] - Sin[ϕu + ϕv])/d];
	Assert[f[1&, sn, u, v] == (-Sin[ϕu]*Δ[ϕv] + Sin[ϕv]*Δ[ϕu] + Sin[ϕu - ϕv])/d];
	Assert[f[cn, 1&, u, v] == m1*(1 + Cos[ϕu] + Cos[ϕv] + Cos[ϕu + ϕv])/d];
	Assert[f[1&, cn, u, v] == m1*(1 + Cos[ϕu] + Cos[ϕv] + Cos[ϕu - ϕv])/d];
	Assert[f[dn, 1&, u, v] == m1*((1 + Δ[ϕu])*(1 + Δ[ϕv]) - m*Sin[ϕu]*Sin[ϕv])/d];
	Assert[f[1&, dn, u, v] == m1*((1 + Δ[ϕu])*(1 + Δ[ϕv]) + m*Sin[ϕu]*Sin[ϕv])/d];
	Assert[f[sn, sn, u, v] == ((1 - Cos[ϕu])*(1 - Δ[ϕv]) - (1 - Cos[ϕv])*(1 - Δ[ϕu]))/(m*d)];
	Assert[f[sn, cn, u, v] == m1*(+Sin[ϕu]*(1 - Δ[ϕv]) + Sin[ϕv]*(1 - Δ[ϕu]))/(m*d)];
	Assert[f[sn, dn, u, v] == m1*(+Sin[ϕu] + Sin[ϕv] - Sin[ϕu + ϕv])/d];
	Assert[f[cn, sn, u, v] == m1*(-Sin[ϕu]*(1 - Δ[ϕv]) + Sin[ϕv]*(1 - Δ[ϕu]))/(m*d)];
	Assert[f[cn, cn, u, v] == m1*((Δ[ϕu] + Cos[ϕu])*(Δ[ϕv] + Cos[ϕv]) - m1*(1 + Cos[ϕu])*(1 + Cos[ϕv]))/(m*d)];
	Assert[f[cn, dn, u, v] == m1*((Δ[ϕu] + Cos[ϕu])*(Δ[ϕv] + Cos[ϕv]) - m1*Sin[ϕu]*Sin[ϕv])/d];
	Assert[f[dn, sn, u, v] == m1*(-Sin[ϕu] + Sin[ϕv] + Sin[ϕu - ϕv])/d];
	Assert[f[dn, cn, u, v] == m1*((Δ[ϕu] + Cos[ϕu])*(Δ[ϕv] + Cos[ϕv]) + m1*Sin[ϕu]*Sin[ϕv])/d];
	Assert[f[dn, dn, u, v] == m1*(m1*(1 + Δ[ϕu])*(1 + Δ[ϕv]) + m*(Δ[ϕu] + Cos[ϕu])*(Δ[ϕv] + Cos[ϕv]))/d];

	Clear[f, d];

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

	(* Theta function notation *)
	a = RandomReal[{0, 1}, WorkingPrecision -> pr];
	Assert[f[a, 0, m] == (2 - a)*Pi/2 + 2*Arg@EllipticTheta[1, (I*Log@EllipticNomeQ[m] + a*Pi)/2, EllipticNomeQ[m]^2]];

	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

Last edited by lanxiyu (2025-11-26 12:15:43)