feat: add Sine + Cosine + Signbit

gnovm/docs/go-gno-compatibility.md

@@ -233,7 +233,7 @@ Legend: full, partial, missing, TBD.
 | log/slog/internal                           | TBD      |
 | log/syslog                                  | TBD      |
 | maps                                        | TBD      |
-| math                                        | TBD      |
+| math                                        | partial      |
 | math/big                                    | TBD      |
 | math/bits                                   | TBD      |
 | math/cmplx                                  | TBD      |

gnovm/stdlibs/math/all_test.gno

@@ -0,0 +1,146 @@
+// Copyright 2009 The Go 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 math
+
+import (
+	"testing"
+)
+
+// The expected results below were computed by the high precision calculators
+// at https://keisan.casio.com/.  More exact input values (array vf[], above)
+// were obtained by printing them with "%.26f".  The answers were calculated
+// to 26 digits (by using the "Digit number" drop-down control of each
+// calculator).
+var vf = []float64{
+	4.9790119248836735e+00,
+	7.7388724745781045e+00,
+	-2.7688005719200159e-01,
+	-5.0106036182710749e+00,
+	9.6362937071984173e+00,
+	2.9263772392439646e+00,
+	5.2290834314593066e+00,
+	2.7279399104360102e+00,
+	1.8253080916808550e+00,
+	-8.6859247685756013e+00,
+}
+
+var signbit = []bool{
+	false,
+	false,
+	true,
+	true,
+	false,
+	false,
+	false,
+	false,
+	false,
+	true,
+}
+
+var sin = []float64{
+	-9.6466616586009283766724726e-01,
+	9.9338225271646545763467022e-01,
+	-2.7335587039794393342449301e-01,
+	9.5586257685042792878173752e-01,
+	-2.099421066779969164496634e-01,
+	2.135578780799860532750616e-01,
+	-8.694568971167362743327708e-01,
+	4.019566681155577786649878e-01,
+	9.6778633541687993721617774e-01,
+	-6.734405869050344734943028e-01,
+}
+
+var vfsignbitSC = []float64{
+	Inf(-1),
+	Copysign(0, -1),
+	0,
+	Inf(1),
+	NaN(),
+}
+
+var signbitSC = []bool{
+	true,
+	true,
+	false,
+	false,
+	false,
+}
+
+var vfsinSC = []float64{
+	Inf(-1),
+	Copysign(0, -1),
+	0,
+	Inf(1),
+	NaN(),
+}
+
+var sinSC = []float64{
+	NaN(),
+	Copysign(0, -1),
+	0,
+	NaN(),
+	NaN(),
+}
+
+func tolerance(a, b, e float64) bool {
+	// Multiplying by e here can underflow denormal values to zero.
+	// Check a==b so that at least if a and b are small and identical
+	// we say they match.
+	if a == b {
+		return true
+	}
+	d := a - b
+	if d < 0 {
+		d = -d
+	}
+
+	// note: b is correct (expected) value, a is actual value.
+	// make error tolerance a fraction of b, not a.
+	if b != 0 {
+		e = e * b
+		if e < 0 {
+			e = -e
+		}
+	}
+	return d < e
+}
+func close(a, b float64) bool      { return tolerance(a, b, 1e-14) }
+func veryclose(a, b float64) bool  { return tolerance(a, b, 4e-16) }
+func soclose(a, b, e float64) bool { return tolerance(a, b, e) }
+func alike(a, b float64) bool {
+	switch {
+	case IsNaN(a) && IsNaN(b):
+		return true
+	case a == b:
+		return Signbit(a) == Signbit(b)
+	}
+	return false
+}
+
+func TestSin(t *testing.T) {
+	for i := 0; i < len(vf); i++ {
+		if f := Sin(vf[i]); !veryclose(sin[i], f) {
+			t.Errorf("Sin(%g) = %g, want %g", vf[i], f, sin[i])
+		}
+	}
+	for i := 0; i < len(vfsinSC); i++ {
+		if f := Sin(vfsinSC[i]); !alike(sinSC[i], f) {
+			t.Errorf("Sin(%g) = %g, want %g", vfsinSC[i], f, sinSC[i])
+		}
+	}
+}
+
+func TestSignbit(t *testing.T) {
+	for i := 0; i < len(vf); i++ {
+		if f := Signbit(vf[i]); signbit[i] != f {
+			t.Errorf("Signbit(%g) = %t, want %t", vf[i], f, signbit[i])
+		}
+	}
+	for i := 0; i < len(vfsignbitSC); i++ {
+		if f := Signbit(vfsignbitSC[i]); signbitSC[i] != f {
+			t.Errorf("Signbit(%g) = %t, want %t", vfsignbitSC[i], f, signbitSC[i])
+		}
+	}
+}

gnovm/stdlibs/math/signbit.gno

@@ -0,0 +1,10 @@
+package math
+
+import (
+	imath "internal/math"
+)
+
+// Signbit reports whether x is negative or negative zero.
+func Signbit(x float64) bool {
+	return imath.Float64bits(x)&(1<<63) != 0
+}

gnovm/stdlibs/math/sin.gno

@@ -0,0 +1,227 @@
+// Copyright 2011 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+// from https://github.com/golang/go/blob/e4de2e7fd04c92d4035cd268d5043f2380aef437/src/pkg/math/sin.go
+
+package math
+
+/*
+	Floating-point sine and cosine.
+*/
+
+// The original C code, the long comment, and the constants
+// below were from http://netlib.sandia.gov/cephes/cmath/sin.c,
+// available from http://www.netlib.org/cephes/cmath.tgz.
+// The go code is a simplified version of the original C.
+//
+//      sin.c
+//
+//      Circular sine
+//
+// SYNOPSIS:
+//
+// double x, y, sin();
+// y = sin( x );
+//
+// DESCRIPTION:
+//
+// Range reduction is into intervals of pi/4.  The reduction error is nearly
+// eliminated by contriving an extended precision modular arithmetic.
+//
+// Two polynomial approximating functions are employed.
+// Between 0 and pi/4 the sine is approximated by
+//      x  +  x**3 P(x**2).
+// Between pi/4 and pi/2 the cosine is represented as
+//      1  -  x**2 Q(x**2).
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain      # trials      peak         rms
+//    DEC       0, 10       150000       3.0e-17     7.8e-18
+//    IEEE -1.07e9,+1.07e9  130000       2.1e-16     5.4e-17
+//
+// Partial loss of accuracy begins to occur at x = 2**30 = 1.074e9.  The loss
+// is not gradual, but jumps suddenly to about 1 part in 10e7.  Results may
+// be meaningless for x > 2**49 = 5.6e14.
+//
+//      cos.c
+//
+//      Circular cosine
+//
+// SYNOPSIS:
+//
+// double x, y, cos();
+// y = cos( x );
+//
+// DESCRIPTION:
+//
+// Range reduction is into intervals of pi/4.  The reduction error is nearly
+// eliminated by contriving an extended precision modular arithmetic.
+//
+// Two polynomial approximating functions are employed.
+// Between 0 and pi/4 the cosine is approximated by
+//      1  -  x**2 Q(x**2).
+// Between pi/4 and pi/2 the sine is represented as
+//      x  +  x**3 P(x**2).
+//
+// ACCURACY:
+//
+//                      Relative error:
+// arithmetic   domain      # trials      peak         rms
+//    IEEE -1.07e9,+1.07e9  130000       2.1e-16     5.4e-17
+//    DEC        0,+1.07e9   17000       3.0e-17     7.2e-18
+//
+// Cephes Math Library Release 2.8:  June, 2000
+// Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
+//
+// The readme file at http://netlib.sandia.gov/cephes/ says:
+//    Some software in this archive may be from the book _Methods and
+// Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
+// International, 1989) or from the Cephes Mathematical Library, a
+// commercial product. In either event, it is copyrighted by the author.
+// What you see here may be used freely but it comes with no support or
+// guarantee.
+//
+//   The two known misprints in the book are repaired here in the
+// source listings for the gamma function and the incomplete beta
+// integral.
+//
+//   Stephen L. Moshier
+//   [email protected]
+
+// sin coefficients
+var _sin = [...]float64{
+	1.58962301576546568060e-10, // 0x3de5d8fd1fd19ccd
+	-2.50507477628578072866e-8, // 0xbe5ae5e5a9291f5d
+	2.75573136213857245213e-6,  // 0x3ec71de3567d48a1
+	-1.98412698295895385996e-4, // 0xbf2a01a019bfdf03
+	8.33333333332211858878e-3,  // 0x3f8111111110f7d0
+	-1.66666666666666307295e-1, // 0xbfc5555555555548
+}
+
+// cos coefficients
+var _cos = [...]float64{
+	-1.13585365213876817300e-11, // 0xbda8fa49a0861a9b
+	2.08757008419747316778e-9,   // 0x3e21ee9d7b4e3f05
+	-2.75573141792967388112e-7,  // 0xbe927e4f7eac4bc6
+	2.48015872888517045348e-5,   // 0x3efa01a019c844f5
+	-1.38888888888730564116e-3,  // 0xbf56c16c16c14f91
+	4.16666666666665929218e-2,   // 0x3fa555555555554b
+}
+
+// Cos returns the cosine of x.
+//
+// Special cases are:
+//
+//	Cos(±Inf) = NaN
+//	Cos(NaN) = NaN
+func Cos(x float64) float64 {
+	const (
+		PI4A = 7.85398125648498535156e-1                             // 0x3fe921fb40000000, Pi/4 split into three parts
+		PI4B = 3.77489470793079817668e-8                             // 0x3e64442d00000000,
+		PI4C = 2.69515142907905952645e-15                            // 0x3ce8469898cc5170,
+		M4PI = 1.273239544735162542821171882678754627704620361328125 // 4/pi
+	)
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
+	// special cases
+	switch {
+	case x != x || x < -MaxFloat64 || x > MaxFloat64: // IsNaN(x) || IsInf(x, 0):
+		return NaN()
+	}
+
+	// make argument positive
+	sign := false
+	if x < 0 {
+		x = -x
+	}
+
+	j := int64(x * M4PI) // integer part of x/(Pi/4), as integer for tests on the phase angle
+	y := float64(j)      // integer part of x/(Pi/4), as float
+
+	// map zeros to origin
+	if j&1 == 1 {
+		j += 1
+		y += 1
+	}
+	j &= 7 // octant modulo 2Pi radians (360 degrees)
+	if j > 3 {
+		j -= 4
+		sign = !sign
+	}
+	if j > 1 {
+		sign = !sign
+	}
+
+	z := ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
+	zz := z * z
+	if j == 1 || j == 2 {
+		y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
+	} else {
+		y = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5])
+	}
+	if sign {
+		y = -y
+	}
+	return y
+}
+
+// Sin returns the sine of x.
+//
+// Special cases are:
+//
+//	Sin(±0) = ±0
+//	Sin(±Inf) = NaN
+//	Sin(NaN) = NaN
+func Sin(x float64) float64 {
+	const (
+		PI4A = 7.85398125648498535156e-1                             // 0x3fe921fb40000000, Pi/4 split into three parts
+		PI4B = 3.77489470793079817668e-8                             // 0x3e64442d00000000,
+		PI4C = 2.69515142907905952645e-15                            // 0x3ce8469898cc5170,
+		M4PI = 1.273239544735162542821171882678754627704620361328125 // 4/pi
+	)
+	// TODO(rsc): Remove manual inlining of IsNaN, IsInf
+	// when compiler does it for us
+	// special cases
+	switch {
+	case x == 0 || x != x: // x == 0 || IsNaN():
+		return x // return ±0 || NaN()
+	case x < -MaxFloat64 || x > MaxFloat64: // IsInf(x, 0):
+		return NaN()
+	}
+
+	// make argument positive but save the sign
+	sign := false
+	if x < 0 {
+		x = -x
+		sign = true
+	}
+
+	j := int64(x * M4PI) // integer part of x/(Pi/4), as integer for tests on the phase angle
+	y := float64(j)      // integer part of x/(Pi/4), as float
+
+	// map zeros to origin
+	if j&1 == 1 {
+		j += 1
+		y += 1
+	}
+	j &= 7 // octant modulo 2Pi radians (360 degrees)
+	// reflect in x axis
+	if j > 3 {
+		sign = !sign
+		j -= 4
+	}
+
+	z := ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic
+	zz := z * z
+	if j == 1 || j == 2 {
+		y = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5])
+	} else {
+		y = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5])
+	}
+	if sign {
+		y = -y
+	}
+	return y
+}