From 9c251ca02972205ba15bd868a57c53380dd468ed Mon Sep 17 00:00:00 2001 From: Schwarf <32012898+Schwarf@users.noreply.github.com> Date: Wed, 24 Sep 2025 21:44:24 +0200 Subject: [PATCH] mathext: add dilogarithm function Li2 --- mathext/dilog.go | 66 ++++++++++++++++++++++++++++++ mathext/dilog_test.go | 93 +++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 159 insertions(+) create mode 100644 mathext/dilog.go create mode 100644 mathext/dilog_test.go diff --git a/mathext/dilog.go b/mathext/dilog.go new file mode 100644 index 00000000..048b54d7 --- /dev/null +++ b/mathext/dilog.go @@ -0,0 +1,66 @@ +// Copyright ©2025 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package mathext + +import ( + "math" + "math/cmplx" +) + +// Li2 returns the dilogarithm Li2(z) on the principal branch. +// +// For |z| < 1, Li2(z) is defined by the power series +// +// Li2(z) = SUM_{k=1}^{infinity} z^k / k^2 +// +// and is analytically continued to the rest of the complex plane. +// The implementation uses reflection and inversion identities to map z into +// a region where the series converges rapidly, then evaluates the series. +// +// Branch cut: Li2 has a logarithmic branch point at z=1 with the standard +// cut on the real axis for z in the interval (1, infinity). The principal value is taken with +// Arg(z) element of (−Pi, Pi]. +func Li2(z complex128) complex128 { + // Special cases + if z == 0 { + return 0 + } + if z == 1 { + return complex(math.Pi*math.Pi/6, 0) + } + if z == -1 { + return complex(-math.Pi*math.Pi/12, 0) + } + + // Reflection: map Re(z) > 0.5 into left half-plane for better convergence + // This formula is applied before inversion on the principal branch and gives very accurate results + // for real z > 1 because Li2(-1) and similar values are known exactly. + // Inversion first would also be valid, but is more sensitive to the + // Arg(-z)= +/-Pi branch choice in log(-z) and can yield wrong imaginary signs + // if care is not taken. + if real(z) > 0.5 { + return complex(math.Pi*math.Pi/6, 0) - complex128(cmplx.Log(z)*cmplx.Log(1-z)) - Li2(1-z) + } + + // Inversion: map |z| > 1 into unit disk + if cmplx.Abs(z) > 1 { + logmz := cmplx.Log(-z) + return -complex(math.Pi*math.Pi/6, 0) - complex128(0.5*logmz*logmz) - Li2(1/z) + } + // Direct series for |z| <= 1 and Re(z) <= 0.5 + return li2Series(z) +} + +func li2Series(z complex128) complex128 { + const tol = 1e-15 + var sum complex128 + zk := z // zk = z^k + for k := 1.0; cmplx.Abs(zk)/(k*k) > tol; k++ { + k2 := k * k + sum += complex(real(zk)/k2, imag(zk)/k2) + zk *= z + } + return sum +} diff --git a/mathext/dilog_test.go b/mathext/dilog_test.go new file mode 100644 index 00000000..e80861b9 --- /dev/null +++ b/mathext/dilog_test.go @@ -0,0 +1,93 @@ +// Copyright ©2025 The Gonum Authors. All rights reserved. +// Use of this source code is governed by a BSD-style +// license that can be found in the LICENSE file. + +package mathext + +import ( + "math" + "math/cmplx" + "testing" + + "gonum.org/v1/gonum/cmplxs/cscalar" +) + +func TestDiLogValues(t *testing.T) { + t.Parallel() + const tol = 1e-12 + + for _, test := range []struct { + in, want complex128 + }{ + // Reference values were generated using Python's scipy.special.spence + // (where Li2(z) = spence(1 - z)) and verified with a computer algebra system. + + // well known values + {in: 0, want: 0}, + {in: 1, want: math.Pi * math.Pi / 6}, + {in: -1, want: -math.Pi * math.Pi / 12}, + {in: 0.5, want: math.Pi*math.Pi/12 - math.Ln2*math.Ln2/2}, + {in: 2, want: complex(math.Pi*math.Pi/4, -math.Pi*math.Ln2)}, + + // abs(z) < 0.5 + {in: 0.1 + 0.1i, want: 0.099751196798713 + 0.105220383905600i}, + {in: -0.3 + 0.39i, want: -0.306174046754339 + 0.337550668621259i}, + {in: 0.001 - 0.49i, want: -0.055831393285929 - 0.478156430372353i}, + + // 0.5 < abs(z) < 1 + {in: -0.9999 + 0.001i, want: -0.8223978143206974 + 0.0006931664732411i}, + {in: 0.5 + 0.7i, want: 0.359364350522653 + 0.856767712327590i}, + {in: complex(math.Pi/4, -math.Pi/7), want: 0.828086461155377 - 0.739345385309341i}, + {in: -0.8 - 0.0001i, want: -0.679781588954542 - 0.000073473333084i}, + + // abs(z) > 1 + {in: -1.1 + 0.1i, want: -0.8917388814454027 + 0.0674285967726009i}, + {in: 5 + 0i, want: 1.783719161266631 - 5.056198322111862i}, + {in: -10 + 0i, want: -4.198277886858104 + 0i}, + {in: 1000 + 10000i, want: -42.71073756884990 + 15.39396088869304i}, + {in: -1791.91931 + 0.5i, want: -29.70223568904652 + 0.00209038439188i}, + } { + got := Li2(test.in) + diff := cmplx.Abs(got - test.want) + if cmplx.Abs(test.want) != 0 { + if !cscalar.EqualWithinRel(got, test.want, tol) { + t.Errorf("Li2(%g) relative error %g exceeds tol %g", test.in, diff/cmplx.Abs(test.want), tol) + } + } else if !cscalar.EqualWithinAbs(got, test.want, tol) { + t.Errorf("Li2(%g) abs error %g exceeds tol %g", test.in, diff, tol) + } + } +} + +func TestDiLogProperties(t *testing.T) { + t.Parallel() + const tol = 1e-12 + + // Duplication formula: Li2(z^2) = 2*(Li2(z) + Li2(-z)) + for i, z := range []complex128{ + 0.1 + 0.1i, + -0.3 + 0.4i, + 0.5, + 200, + } { + lhs := Li2(z * z) + rhs := 2 * (Li2(z) + Li2(-z)) + if !cscalar.EqualWithinAbs(lhs, rhs, tol) { + t.Errorf("duplication formula failed for case %d, z=%v: got %v want %v", i, z, lhs, rhs) + } + } + + // Conjugation symmetry: Li2(conj(z)) == conj(Li2(z)) + // Valid for all z not on the branch cut (real z > 1). + for i, z := range []complex128{ + 0.3 + 0.5i, + -0.8 + 0.1i, + 2 + 0.7i, + } { + lhs := Li2(cmplx.Conj(z)) + rhs := cmplx.Conj(Li2(z)) + if !cscalar.EqualWithinAbs(lhs, rhs, tol) { + t.Errorf("conjugation symmetry failed for case %d, z=%v: got %v want %v", i, z, lhs, rhs) + } + } +}