mathext: add dilogarithm function Li2

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
+}

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)
+		}
+	}
+}