From 69ef2f751c71508075700db7cd12eee86d0eeb67 Mon Sep 17 00:00:00 2001 From: n3wbie Date: Tue, 7 Feb 2023 17:42:52 +0900 Subject: [PATCH] feat: additional math libraries --- stdlibs/math/frexp.gno | 43 +++++++++++ stdlibs/math/log.gno | 129 +++++++++++++++++++++++++++++++++ stdlibs/math/log10.gno | 37 ++++++++++ stdlibs/math/pow.gno | 161 +++++++++++++++++++++++++++++++++++++++++ stdlibs/math/sqrt.gno | 149 ++++++++++++++++++++++++++++++++++++++ 5 files changed, 519 insertions(+) create mode 100644 stdlibs/math/frexp.gno create mode 100644 stdlibs/math/log.gno create mode 100644 stdlibs/math/log10.gno create mode 100644 stdlibs/math/pow.gno create mode 100644 stdlibs/math/sqrt.gno diff --git a/stdlibs/math/frexp.gno b/stdlibs/math/frexp.gno new file mode 100644 index 00000000000..ab4c2809722 --- /dev/null +++ b/stdlibs/math/frexp.gno @@ -0,0 +1,43 @@ +// 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 ( + imath "internal/math" +) + +// Frexp breaks f into a normalized fraction +// and an integral power of two. +// It returns frac and exp satisfying f == frac × 2**exp, +// with the absolute value of frac in the interval [½, 1). +// +// Special cases are: +// +// Frexp(±0) = ±0, 0 +// Frexp(±Inf) = ±Inf, 0 +// Frexp(NaN) = NaN, 0 +func Frexp(f float64) (frac float64, exp int) { + // if haveArchFrexp { + // return archFrexp(f) + // } + return frexp(f) +} + +func frexp(f float64) (frac float64, exp int) { + // special cases + switch { + case f == 0: + return f, 0 // correctly return -0 + case IsInf(f, 0) || IsNaN(f): + return f, 0 + } + f, exp = normalize(f) + x := imath.Float64bits(f) + exp += int((x>>shift)&mask) - bias + 1 + x &^= mask << shift + x |= (-1 + bias) << shift + frac = imath.Float64frombits(x) + return +} \ No newline at end of file diff --git a/stdlibs/math/log.gno b/stdlibs/math/log.gno new file mode 100644 index 00000000000..51309645656 --- /dev/null +++ b/stdlibs/math/log.gno @@ -0,0 +1,129 @@ +// 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 + +/* + Floating-point logarithm. +*/ + +// The original C code, the long comment, and the constants +// below are from FreeBSD's /usr/src/lib/msun/src/e_log.c +// and came with this notice. The go code is a simpler +// version of the original C. +// +// ==================================================== +// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +// +// Developed at SunPro, a Sun Microsystems, Inc. business. +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== +// +// __ieee754_log(x) +// Return the logarithm of x +// +// Method : +// 1. Argument Reduction: find k and f such that +// x = 2**k * (1+f), +// where sqrt(2)/2 < 1+f < sqrt(2) . +// +// 2. Approximation of log(1+f). +// Let s = f/(2+f) ; based on log(1+f) = log(1+s) - log(1-s) +// = 2s + 2/3 s**3 + 2/5 s**5 + ....., +// = 2s + s*R +// We use a special Reme algorithm on [0,0.1716] to generate +// a polynomial of degree 14 to approximate R. The maximum error +// of this polynomial approximation is bounded by 2**-58.45. In +// other words, +// 2 4 6 8 10 12 14 +// R(z) ~ L1*s +L2*s +L3*s +L4*s +L5*s +L6*s +L7*s +// (the values of L1 to L7 are listed in the program) and +// | 2 14 | -58.45 +// | L1*s +...+L7*s - R(z) | <= 2 +// | | +// Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2. +// In order to guarantee error in log below 1ulp, we compute log by +// log(1+f) = f - s*(f - R) (if f is not too large) +// log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy) +// +// 3. Finally, log(x) = k*Ln2 + log(1+f). +// = k*Ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*Ln2_lo))) +// Here Ln2 is split into two floating point number: +// Ln2_hi + Ln2_lo, +// where n*Ln2_hi is always exact for |n| < 2000. +// +// Special cases: +// log(x) is NaN with signal if x < 0 (including -INF) ; +// log(+INF) is +INF; log(0) is -INF with signal; +// log(NaN) is that NaN with no signal. +// +// Accuracy: +// according to an error analysis, the error is always less than +// 1 ulp (unit in the last place). +// +// Constants: +// The hexadecimal values are the intended ones for the following +// constants. The decimal values may be used, provided that the +// compiler will convert from decimal to binary accurately enough +// to produce the hexadecimal values shown. + +// Log returns the natural logarithm of x. +// +// Special cases are: +// +// Log(+Inf) = +Inf +// Log(0) = -Inf +// Log(x < 0) = NaN +// Log(NaN) = NaN +func Log(x float64) float64 { + // if haveArchLog { + // return archLog(x) + // } + return log(x) +} + +func log(x float64) float64 { + const ( + Ln2Hi = 6.93147180369123816490e-01 /* 3fe62e42 fee00000 */ + Ln2Lo = 1.90821492927058770002e-10 /* 3dea39ef 35793c76 */ + L1 = 6.666666666666735130e-01 /* 3FE55555 55555593 */ + L2 = 3.999999999940941908e-01 /* 3FD99999 9997FA04 */ + L3 = 2.857142874366239149e-01 /* 3FD24924 94229359 */ + L4 = 2.222219843214978396e-01 /* 3FCC71C5 1D8E78AF */ + L5 = 1.818357216161805012e-01 /* 3FC74664 96CB03DE */ + L6 = 1.531383769920937332e-01 /* 3FC39A09 D078C69F */ + L7 = 1.479819860511658591e-01 /* 3FC2F112 DF3E5244 */ + ) + + // special cases + switch { + case IsNaN(x) || IsInf(x, 1): + return x + case x < 0: + return NaN() + case x == 0: + return Inf(-1) + } + + // reduce + f1, ki := Frexp(x) + if f1 < Sqrt2/2 { + f1 *= 2 + ki-- + } + f := f1 - 1 + k := float64(ki) + + // compute + s := f / (2 + f) + s2 := s * s + s4 := s2 * s2 + t1 := s2 * (L1 + s4*(L3+s4*(L5+s4*L7))) + t2 := s4 * (L2 + s4*(L4+s4*L6)) + R := t1 + t2 + hfsq := 0.5 * f * f + return k*Ln2Hi - ((hfsq - (s*(hfsq+R) + k*Ln2Lo)) - f) +} \ No newline at end of file diff --git a/stdlibs/math/log10.gno b/stdlibs/math/log10.gno new file mode 100644 index 00000000000..1b1d8db9054 --- /dev/null +++ b/stdlibs/math/log10.gno @@ -0,0 +1,37 @@ +// 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 + +// Log10 returns the decimal logarithm of x. +// The special cases are the same as for Log. +func Log10(x float64) float64 { + // if haveArchLog10 { + // return archLog10(x) + // } + return log10(x) +} + +func log10(x float64) float64 { + return Log(x) * (1 / Ln10) +} + +// Log2 returns the binary logarithm of x. +// The special cases are the same as for Log. +func Log2(x float64) float64 { + // if haveArchLog2 { + // return archLog2(x) + // } + return log2(x) +} + +func log2(x float64) float64 { + frac, exp := Frexp(x) + // Make sure exact powers of two give an exact answer. + // Don't depend on Log(0.5)*(1/Ln2)+exp being exactly exp-1. + if frac == 0.5 { + return float64(exp - 1) + } + return Log(frac)*(1/Ln2) + float64(exp) +} \ No newline at end of file diff --git a/stdlibs/math/pow.gno b/stdlibs/math/pow.gno new file mode 100644 index 00000000000..00549408bb0 --- /dev/null +++ b/stdlibs/math/pow.gno @@ -0,0 +1,161 @@ +// 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 ( + imath "internal/math" +) + +func isOddInt(x float64) bool { + xi, xf := Modf(x) + return xf == 0 && int64(xi)&1 == 1 +} + +// Special cases taken from FreeBSD's /usr/src/lib/msun/src/e_pow.c +// updated by IEEE Std. 754-2008 "Section 9.2.1 Special values". + +// Pow returns x**y, the base-x exponential of y. +// +// Special cases are (in order): +// +// Pow(x, ±0) = 1 for any x +// Pow(1, y) = 1 for any y +// Pow(x, 1) = x for any x +// Pow(NaN, y) = NaN +// Pow(x, NaN) = NaN +// Pow(±0, y) = ±Inf for y an odd integer < 0 +// Pow(±0, -Inf) = +Inf +// Pow(±0, +Inf) = +0 +// Pow(±0, y) = +Inf for finite y < 0 and not an odd integer +// Pow(±0, y) = ±0 for y an odd integer > 0 +// Pow(±0, y) = +0 for finite y > 0 and not an odd integer +// Pow(-1, ±Inf) = 1 +// Pow(x, +Inf) = +Inf for |x| > 1 +// Pow(x, -Inf) = +0 for |x| > 1 +// Pow(x, +Inf) = +0 for |x| < 1 +// Pow(x, -Inf) = +Inf for |x| < 1 +// Pow(+Inf, y) = +Inf for y > 0 +// Pow(+Inf, y) = +0 for y < 0 +// Pow(-Inf, y) = Pow(-0, -y) +// Pow(x, y) = NaN for finite x < 0 and finite non-integer y +func Pow(x, y float64) float64 { + // if haveArchPow { + // return archPow(x, y) + // } + return pow(x, y) +} + +func pow(x, y float64) float64 { + switch { + case y == 0 || x == 1: + return 1 + case y == 1: + return x + case IsNaN(x) || IsNaN(y): + return NaN() + case x == 0: + switch { + case y < 0: + if isOddInt(y) { + return Copysign(Inf(1), x) + } + return Inf(1) + case y > 0: + if isOddInt(y) { + return x + } + return 0 + } + case IsInf(y, 0): + switch { + case x == -1: + return 1 + case (Abs(x) < 1) == IsInf(y, 1): + return 0 + default: + return Inf(1) + } + case IsInf(x, 0): + if IsInf(x, -1) { + return Pow(1/x, -y) // Pow(-0, -y) + } + switch { + case y < 0: + return 0 + case y > 0: + return Inf(1) + } + case y == 0.5: + return Sqrt(x) + case y == -0.5: + return 1 / Sqrt(x) + } + + yi, yf := Modf(Abs(y)) + if yf != 0 && x < 0 { + return NaN() + } + if yi >= 1<<63 { + // yi is a large even int that will lead to overflow (or underflow to 0) + // for all x except -1 (x == 1 was handled earlier) + switch { + case x == -1: + return 1 + case (Abs(x) < 1) == (y > 0): + return 0 + default: + return Inf(1) + } + } + + // ans = a1 * 2**ae (= 1 for now). + a1 := 1.0 + ae := 0 + + // ans *= x**yf + if yf != 0 { + if yf > 0.5 { + yf-- + yi++ + } + a1 = Exp(yf * Log(x)) + } + + // ans *= x**yi + // by multiplying in successive squarings + // of x according to bits of yi. + // accumulate powers of two into exp. + x1, xe := Frexp(x) + for i := int64(yi); i != 0; i >>= 1 { + if xe < -1<<12 || 1<<12 < xe { + // catch xe before it overflows the left shift below + // Since i !=0 it has at least one bit still set, so ae will accumulate xe + // on at least one more iteration, ae += xe is a lower bound on ae + // the lower bound on ae exceeds the size of a float64 exp + // so the final call to Ldexp will produce under/overflow (0/Inf) + ae += xe + break + } + if i&1 == 1 { + a1 *= x1 + ae += xe + } + x1 *= x1 + xe <<= 1 + if x1 < .5 { + x1 += x1 + xe-- + } + } + + // ans = a1*2**ae + // if y < 0 { ans = 1 / ans } + // but in the opposite order + if y < 0 { + a1 = 1 / a1 + ae = -ae + } + return Ldexp(a1, ae) +} \ No newline at end of file diff --git a/stdlibs/math/sqrt.gno b/stdlibs/math/sqrt.gno new file mode 100644 index 00000000000..4731337bed9 --- /dev/null +++ b/stdlibs/math/sqrt.gno @@ -0,0 +1,149 @@ +// 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 ( + imath "internal/math" +) + +// The original C code and the long comment below are +// from FreeBSD's /usr/src/lib/msun/src/e_sqrt.c and +// came with this notice. The go code is a simplified +// version of the original C. +// +// ==================================================== +// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. +// +// Developed at SunPro, a Sun Microsystems, Inc. business. +// Permission to use, copy, modify, and distribute this +// software is freely granted, provided that this notice +// is preserved. +// ==================================================== +// +// __ieee754_sqrt(x) +// Return correctly rounded sqrt. +// ----------------------------------------- +// | Use the hardware sqrt if you have one | +// ----------------------------------------- +// Method: +// Bit by bit method using integer arithmetic. (Slow, but portable) +// 1. Normalization +// Scale x to y in [1,4) with even powers of 2: +// find an integer k such that 1 <= (y=x*2**(2k)) < 4, then +// sqrt(x) = 2**k * sqrt(y) +// 2. Bit by bit computation +// Let q = sqrt(y) truncated to i bit after binary point (q = 1), +// i 0 +// i+1 2 +// s = 2*q , and y = 2 * ( y - q ). (1) +// i i i i +// +// To compute q from q , one checks whether +// i+1 i +// +// -(i+1) 2 +// (q + 2 ) <= y. (2) +// i +// -(i+1) +// If (2) is false, then q = q ; otherwise q = q + 2 . +// i+1 i i+1 i +// +// With some algebraic manipulation, it is not difficult to see +// that (2) is equivalent to +// -(i+1) +// s + 2 <= y (3) +// i i +// +// The advantage of (3) is that s and y can be computed by +// i i +// the following recurrence formula: +// if (3) is false +// +// s = s , y = y ; (4) +// i+1 i i+1 i +// +// otherwise, +// -i -(i+1) +// s = s + 2 , y = y - s - 2 (5) +// i+1 i i+1 i i +// +// One may easily use induction to prove (4) and (5). +// Note. Since the left hand side of (3) contain only i+2 bits, +// it is not necessary to do a full (53-bit) comparison +// in (3). +// 3. Final rounding +// After generating the 53 bits result, we compute one more bit. +// Together with the remainder, we can decide whether the +// result is exact, bigger than 1/2ulp, or less than 1/2ulp +// (it will never equal to 1/2ulp). +// The rounding mode can be detected by checking whether +// huge + tiny is equal to huge, and whether huge - tiny is +// equal to huge for some floating point number "huge" and "tiny". +// +// +// Notes: Rounding mode detection omitted. The constants "mask", "shift", +// and "bias" are found in src/math/bits.go + +// Sqrt returns the square root of x. +// +// Special cases are: +// +// Sqrt(+Inf) = +Inf +// Sqrt(±0) = ±0 +// Sqrt(x < 0) = NaN +// Sqrt(NaN) = NaN +func Sqrt(x float64) float64 { + return sqrt(x) +} + +// Note: On systems where Sqrt is a single instruction, the compiler +// may turn a direct call into a direct use of that instruction instead. + +func sqrt(x float64) float64 { + // special cases + switch { + case x == 0 || IsNaN(x) || IsInf(x, 1): + return x + case x < 0: + return NaN() + } + ix := imath.Float64bits(x) + // normalize x + exp := int((ix >> shift) & mask) + if exp == 0 { // subnormal x + for ix&(1<>= 1 // exp = exp/2, exponent of square root + // generate sqrt(x) bit by bit + ix <<= 1 + var q, s uint64 // q = sqrt(x) + r := uint64(1 << (shift + 1)) // r = moving bit from MSB to LSB + for r != 0 { + t := s + r + if t <= ix { + s = t + r + ix -= t + q += r + } + ix <<= 1 + r >>= 1 + } + // final rounding + if ix != 0 { // remainder, result not exact + q += q & 1 // round according to extra bit + } + ix = q>>1 + uint64(exp-1+bias)<