Math Is Fun Forum

lanxiyu

Member

Registered: 2022-05-10

Posts: 64

Mathematica package for Dedekind eta function

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 - 1

Last edited by lanxiyu (2025-06-25 17:42:58)