tensor/internal/core/expint_test.go

// Copyright (c) 2026 Petr Balvín <opensource@petrbalvin.org> (https://petrbalvin.org)
// SPDX-License-Identifier: MIT

package core

import (
	"math"
	"testing"
)

// TestExpIntegralE1 checks E1 against tabulated values on both the
// series branch (x ≤ 2) and the continued-fraction branch (x > 2).
func TestExpIntegralE1(t *testing.T) {
	x := mustFloats(t, []float64{0.1, 0.5, 1, 2, 5, 10})
	got, err := ExpIntegralE1(x)
	if err != nil {
		t.Fatalf("ExpIntegralE1: %v", err)
	}
	want := []float64{
		1.8229239584193906,
		0.5597735947761609,
		0.21938393439552027,
		0.04890051070806112,
		0.0011482955912753257,
		4.156968929685324e-06,
	}
	for i := range want {
		if math.Abs(got.FloatAt(i)-want[i]) > 1e-13*(1+math.Abs(want[i])) {
			t.Fatalf("E1(%v) = %.16g, want %.16g", x.FloatAt(i), got.FloatAt(i), want[i])
		}
	}
	neg, nerr := ExpIntegralE1(mustFloats(t, []float64{-1}))
	if nerr != nil {
		t.Fatalf("ExpIntegralE1(−1): %v", nerr)
	}
	if !math.IsNaN(neg.FloatAt(0)) {
		t.Fatalf("E1(−1) = %v, want NaN outside the domain", neg.FloatAt(0))
	}
}

// TestExpIntegralEi checks Ei on both branches, including the
// reflection Ei(−x) = −E1(x) for x < 0.
func TestExpIntegralEi(t *testing.T) {
	x := mustFloats(t, []float64{-1, -0.5, 0.5, 1, 5})
	got, err := ExpIntegralEi(x)
	if err != nil {
		t.Fatalf("ExpIntegralEi: %v", err)
	}
	want := []float64{
		-0.21938393439552027,
		-0.5597735947761609,
		0.4542199048631725,
		1.8951178163559368,
		40.18527535580318,
	}
	for i := range want {
		if math.Abs(got.FloatAt(i)-want[i]) > 1e-13*(1+math.Abs(want[i])) {
			t.Fatalf("Ei(%v) = %.16g, want %.16g", x.FloatAt(i), got.FloatAt(i), want[i])
		}
	}
}

// TestDigamma checks the polygamma family against exact values:
// ψ(1) = −γ, ψ(2) = 1 − γ, ψ(½) = −γ − 2 ln 2, and the same for ψ′.
func TestDigamma(t *testing.T) {
	x := mustFloats(t, []float64{0.5, 1, 2, 5})
	got, err := Digamma(x)
	if err != nil {
		t.Fatalf("Digamma: %v", err)
	}
	want := []float64{
		-eulerGamma - 2*math.Log(2),
		-eulerGamma,
		1 - eulerGamma,
		1.5061176684318005,
	}
	for i := range want {
		if math.Abs(got.FloatAt(i)-want[i]) > 1e-11*(1+math.Abs(want[i])) {
			t.Fatalf("psi(%v) = %.16g, want %.16g", x.FloatAt(i), got.FloatAt(i), want[i])
		}
	}
	// Negative non-integer via the reflection formula: ψ(−½) =
	// ψ(1.5) + π·cot(π/2) = ψ(1.5) = 0.03648997397857652.
	neg := mustFloats(t, []float64{-0.5})
	got, err = Digamma(neg)
	if err != nil {
		t.Fatalf("Digamma(−½): %v", err)
	}
	if want := 0.03648997397857652; math.Abs(got.FloatAt(0)-want) > 1e-11 {
		t.Fatalf("psi(−½) = %.16g, want %.16g", got.FloatAt(0), want)
	}
	// Poles at non-positive integers return NaN, the same IEEE
	// convention the gamma family applies.
	pole, perr := Digamma(mustFloats(t, []float64{-2}))
	if perr != nil {
		t.Fatalf("Digamma(−2): %v", perr)
	}
	if !math.IsNaN(pole.FloatAt(0)) {
		t.Fatalf("psi(−2) = %v, want NaN at the pole", pole.FloatAt(0))
	}
}

// TestTrigamma checks ψ′ against exact values: ψ′(1) = π²/6,
// ψ′(½) = π²/2, ψ′(2) = 1 − π²/6.
func TestTrigamma(t *testing.T) {
	x := mustFloats(t, []float64{0.5, 1, 2, 5})
	got, err := Trigamma(x)
	if err != nil {
		t.Fatalf("Trigamma: %v", err)
	}
	want := []float64{
		math.Pi * math.Pi / 2,
		math.Pi * math.Pi / 6,
		math.Pi*math.Pi/6 - 1,
		0.22132295573718011, // π²/6 − (1 + ¼ + ¹⁄₉ + ¹⁄₁₆)
	}
	for i := range want {
		if math.Abs(got.FloatAt(i)-want[i]) > 1e-11*(1+math.Abs(want[i])) {
			t.Fatalf("psi'(%v) = %.16g, want %.16g", x.FloatAt(i), got.FloatAt(i), want[i])
		}
	}
}

// TestFresnel checks C and S against high-precision series values,
// the odd symmetry and the ½ limits at large argument.
func TestFresnel(t *testing.T) {
	x := mustFloats(t, []float64{0.5, 1, 2, 4, 6})
	c, err := FresnelC(x)
	if err != nil {
		t.Fatalf("FresnelC: %v", err)
	}
	s, err := FresnelS(x)
	if err != nil {
		t.Fatalf("FresnelS: %v", err)
	}
	// References at 60-digit precision, all pinned at 1e-12: the plain
	// power series covers the first three (x = 4 sits at its boundary)
	// and x = 6 sums the same series in extended precision.
	wantC := []float64{0.4923442258714464, 0.7798934003768228, 0.4882534060753408, 0.4984260330381776, 0.4995314678555011}
	wantS := []float64{0.06473243286000028, 0.4382591473903548, 0.3434156783636982, 0.4205157542469284, 0.4469607612369303}
	for i := range wantC {
		if math.Abs(c.FloatAt(i)-wantC[i]) > 1e-12 {
			t.Fatalf("C(%v) = %.16g, want %.16g", x.FloatAt(i), c.FloatAt(i), wantC[i])
		}
		if math.Abs(s.FloatAt(i)-wantS[i]) > 1e-12 {
			t.Fatalf("S(%v) = %.16g, want %.16g", x.FloatAt(i), s.FloatAt(i), wantS[i])
		}
	}
	// Odd symmetry.
	nx := mustFloats(t, []float64{-1})
	nc, _ := FresnelC(nx)
	if math.Abs(nc.FloatAt(0)+0.7798934003768228) > 1e-12 {
		t.Fatalf("C(−1) = %v, want −C(1)", nc.FloatAt(0))
	}
}