Math Is Fun Forum

lanxiyu

Member

Registered: 2022-05-10

Posts: 67

Mathematica package for Weierstrass half-periods

BeginPackage["WeierstrassHalfPeriods`"];

Unprotect[WeierstrassHalfPeriodRe, WeierstrassHalfPeriodIm, WeierstrassHalfPeriodReIm];

ClearAll[WeierstrassHalfPeriodRe, WeierstrassHalfPeriodIm, WeierstrassHalfPeriodReIm];

Begin["`Private`"];

w[g2_, g3_, f_] /; VectorQ[{g2, g3}, NumericQ] :=
Module[{x, pr = Precision[{g2, g3}], e1, e2, e3, disc, sign, res},
	If[pr == Infinity, Return[$Failed, Module]];
	disc = g2^3 - 27*g3^2;
	pr = Max[pr, Precision[disc]];
	{e1, e2, e3} = x /. {ToRules@NRoots[4*x^3 - g2*x - g3, x, PrecisionGoal -> pr]};
	sign = 0;
	Which[
		!VectorQ[{g2, g3}, PossibleZeroQ @* Im],
			Null,
		Re[disc] < 0,
			{e3, e1, e2} = Sort[{e1, e2, e3}, Im[#1] < Im[#2] &];
			sign = -1,
		True,
			{e3, e2, e1} = Sort[{e1, e2, e3}, Re[#1] < Re[#2] &];
			sign = 1
	];
	x = f[e1, e2, e3, sign];
	If[!MatrixQ[x, NumberQ], Return[$Failed, Module]];
	Return[Total[Pi/(2*ArithmeticGeometricMean @@ #) & /@ x], Module]
];

w[___] := $Failed;

WeierstrassHalfPeriodRe[{g2_, g3_}] :=
Module[{res = $Failed},
	res = w[g2, g3, Which[
		#4 > 0, {{Sqrt[#1 - #3], Sqrt[#1 - #2]}},
		#4 < 0, With[{t = Sqrt[#1 - #2]}, {{Abs@t, Re@t}}]
	] &];
	res /; NumberQ[res]
];

WeierstrassHalfPeriodIm[{g2_, g3_}] :=
Module[{res = $Failed},
	res = w[g2, g3, Which[
		#4 > 0, {-I*{Sqrt[#1 - #3], Sqrt[#2 - #3]}},
		#4 < 0, With[{t = Sqrt[#1 - #2]}, 2*{{Abs@t, Re@t}, -I*{Abs@t, Im@t}}]
	] &];
	res /; NumberQ[res]
];

WeierstrassHalfPeriodReIm[{g2_, g3_}] := {
	WeierstrassHalfPeriodRe[{g2, g3}],
	WeierstrassHalfPeriodIm[{g2, g3}]
};

End[];

SetAttributes[{WeierstrassHalfPeriodRe, WeierstrassHalfPeriodIm, WeierstrassHalfPeriodReIm}, {ReadProtected}];

Protect[{WeierstrassHalfPeriodRe, WeierstrassHalfPeriodIm, WeierstrassHalfPeriodReIm}];

EndPackage[];

Last edited by lanxiyu (2025-12-02 11:27:14)