From 1484d274dfbb159005640dcd9f2aa5dc71c45c39 Mon Sep 17 00:00:00 2001 From: Zack Date: Sun, 7 May 2023 09:59:33 -0700 Subject: [PATCH 1/3] add sin.go This version of sin.go comes from pre Go1. Futures versions rely on unsafe which is not yet implemented in gnovm --- gnovm/stdlibs/math/sin.gno | 227 +++++++++++++++++++++++++++++++++++++ 1 file changed, 227 insertions(+) create mode 100644 gnovm/stdlibs/math/sin.gno diff --git a/gnovm/stdlibs/math/sin.gno b/gnovm/stdlibs/math/sin.gno new file mode 100644 index 00000000000..40a611edf70 --- /dev/null +++ b/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 +// moshier@na-net.ornl.gov + +// 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 +} From 934231785534c4ec0b90e0f8f6d9181dd5e80771 Mon Sep 17 00:00:00 2001 From: Zack Scholl Date: Sun, 14 May 2023 14:54:42 -0700 Subject: [PATCH 2/3] add initial tests and signbit --- gnovm/stdlibs/math/all_test.gno | 146 ++++++++++++++++++++++++++++++++ gnovm/stdlibs/math/signbit.gno | 10 +++ 2 files changed, 156 insertions(+) create mode 100644 gnovm/stdlibs/math/all_test.gno create mode 100644 gnovm/stdlibs/math/signbit.gno diff --git a/gnovm/stdlibs/math/all_test.gno b/gnovm/stdlibs/math/all_test.gno new file mode 100644 index 00000000000..ba7b1528d23 --- /dev/null +++ b/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]) + } + } +} diff --git a/gnovm/stdlibs/math/signbit.gno b/gnovm/stdlibs/math/signbit.gno new file mode 100644 index 00000000000..38c2178932e --- /dev/null +++ b/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 +} From 25ffd8abab40f195e51284f69467e9dc410d564e Mon Sep 17 00:00:00 2001 From: Zack Date: Mon, 22 May 2023 14:59:54 -0700 Subject: [PATCH 3/3] partial math addition --- gnovm/docs/go-gno-compatibility.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/gnovm/docs/go-gno-compatibility.md b/gnovm/docs/go-gno-compatibility.md index a579a96b6cc..2da379f9dda 100644 --- a/gnovm/docs/go-gno-compatibility.md +++ b/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 |