103 lines
3.2 KiB
Go
103 lines
3.2 KiB
Go
// 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")
|
||
}
|
||
}
|