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

package core

import (
	"math"
	"math/big"
	"testing"
)

// cosm1Big evaluates cos(x) − 1 by the Taylor series in 256-bit
// arithmetic, the exact referent the float64 implementation is held
// against across the crossover.
func cosm1Big(x float64) *big.Float {
	const prec = 256
	xf := new(big.Float).SetPrec(prec).SetFloat64(x)
	sq := new(big.Float).SetPrec(prec).Mul(xf, xf)
	sum := new(big.Float).SetPrec(prec).SetInt64(1)
	term := new(big.Float).SetPrec(prec).SetInt64(1)
	for k := 1; k <= 60; k++ {
		term.Mul(term, sq)
		term.Quo(term, new(big.Float).SetPrec(prec).SetInt64(int64(2*k-1)))
		term.Quo(term, new(big.Float).SetPrec(prec).SetInt64(int64(2*k)))
		term.Neg(term)
		sum.Add(sum, term)
	}
	return sum.Sub(sum, new(big.Float).SetPrec(prec).SetInt64(1))
}

// TestCosm1AgainstSeries holds the implementation against the exact
// referent on both sides of the crossover, from arguments whose answer
// is −x²/2 as far as the format can see up to ones where the direct
// subtraction carries it alone.
func TestCosm1AgainstSeries(t *testing.T) {
	points := []float64{
		1e-300, 1e-200, 5e-12, 1e-10, 1e-9, 1e-8, 1e-6, 1e-4, 0.01, 0.1,
		0.3, 0.5, 0.78, math.Pi / 4, 0.79, 1, 2, -0.3, -0.78, -2,
	}
	for _, x := range points {
		got, err := Cosm1(mustFloats(t, []float64{x}))
		if err != nil {
			t.Fatalf("Cosm1(%v): %v", x, err)
		}
		want, _ := cosm1Big(x).Float64()
		v := got.FloatAt(0)
		if want == 0 {
			if v != 0 {
				t.Fatalf("Cosm1(%v) = %v, want 0", x, v)
			}
			continue
		}
		if d := math.Abs(v-want) / math.Abs(want); d > 6e-16 {
			t.Fatalf("Cosm1(%v) = %.17g, want %.17g (relative %.3g)", x, v, want, d)
		}
	}
}

// TestCosm1Pins pins the behaviour the referent cannot speak for: the
// exact zero at the origin, the underflowed answer keeping the minus
// sign the function's range promises, the bit-equality with the direct
// subtraction above the crossover, and the integer promotion.
func TestCosm1Pins(t *testing.T) {
	zero, err := Cosm1(mustFloats(t, []float64{0}))
	if err != nil {
		t.Fatal(err)
	}
	if zero.FloatAt(0) != 0 {
		t.Fatalf("Cosm1(0) = %v, want 0", zero.FloatAt(0))
	}
	tiny, err := Cosm1(mustFloats(t, []float64{1e-300}))
	if err != nil {
		t.Fatal(err)
	}
	if tiny.FloatAt(0) != 0 || !math.Signbit(tiny.FloatAt(0)) {
		t.Fatalf("Cosm1(1e-300) = %v, want the negative zero the true answer underflows to", tiny.FloatAt(0))
	}
	for _, x := range []float64{0.79, 1, 2, 10, 100, 1e6, -3.5} {
		got, err := Cosm1(mustFloats(t, []float64{x}))
		if err != nil {
			t.Fatalf("Cosm1(%v): %v", x, err)
		}
		if want := math.Cos(x) - 1; got.FloatAt(0) != want {
			t.Fatalf("Cosm1(%v) = %.17g, want the direct %.17g", x, got.FloatAt(0), want)
		}
	}
	ints, err := FromInts([]int64{0, 1}, 2)
	if err != nil {
		t.Fatal(err)
	}
	promoted, err := Cosm1(ints)
	if err != nil {
		t.Fatalf("Cosm1 over ints: %v", err)
	}
	if promoted.FloatAt(1) != math.Cos(1)-1 {
		t.Fatalf("Cosm1(int 1) = %.17g, want %.17g", promoted.FloatAt(1), math.Cos(1)-1)
	}
	cx, _ := FromComplexes([]complex128{1}, 1)
	if _, err := Cosm1(cx); err == nil {
		t.Fatal("Cosm1 over complex: want an error")
	}
}