From 4f4a855d82a889cebcfca150a7a43909bcb6a346 Mon Sep 17 00:00:00 2001 From: Ian Lance Taylor Date: Fri, 18 Jan 2019 19:04:36 +0000 Subject: libgo: update to Go1.12beta2 Reviewed-on: https://go-review.googlesource.com/c/158019 gotools/: * Makefile.am (go_cmd_vet_files): Update for Go1.12beta2 release. (GOTOOLS_TEST_TIMEOUT): Increase to 600. (check-runtime): Export LD_LIBRARY_PATH before computing GOARCH and GOOS. (check-vet): Copy golang.org/x/tools into check-vet-dir. * Makefile.in: Regenerate. gcc/testsuite/: * go.go-torture/execute/names-1.go: Stop using debug/xcoff, which is no longer externally visible. From-SVN: r268084 --- libgo/go/math/all_test.go | 76 ++++++++ libgo/go/math/big/arith.go | 2 +- libgo/go/math/big/float.go | 12 +- libgo/go/math/big/int.go | 7 + libgo/go/math/big/int_test.go | 26 +++ libgo/go/math/big/nat.go | 56 +++--- libgo/go/math/big/prime.go | 2 +- libgo/go/math/big/rat.go | 7 + libgo/go/math/big/sqrt.go | 20 +-- libgo/go/math/bits/bits.go | 207 +++++++++++++++++++++- libgo/go/math/bits/bits_test.go | 380 ++++++++++++++++++++++++++++++++++++++++ libgo/go/math/cmplx/isinf.go | 2 +- libgo/go/math/cmplx/isnan.go | 2 +- libgo/go/math/example_test.go | 22 +++ libgo/go/math/export_test.go | 3 + libgo/go/math/huge_test.go | 99 +++++++++++ libgo/go/math/log1p.go | 5 +- libgo/go/math/mod.go | 8 +- libgo/go/math/pow.go | 12 +- libgo/go/math/signbit.go | 2 +- libgo/go/math/sin.go | 64 ++++--- libgo/go/math/sincos.go | 29 +-- libgo/go/math/sincos_386.go | 15 -- libgo/go/math/sinh.go | 2 +- libgo/go/math/tan.go | 28 +-- libgo/go/math/trig_reduce.go | 94 ++++++++++ libgo/go/math/unsafe.go | 20 ++- 27 files changed, 1062 insertions(+), 140 deletions(-) create mode 100644 libgo/go/math/huge_test.go delete mode 100644 libgo/go/math/sincos_386.go create mode 100644 libgo/go/math/trig_reduce.go (limited to 'libgo/go/math') diff --git a/libgo/go/math/all_test.go b/libgo/go/math/all_test.go index 00f2058..ed42941 100644 --- a/libgo/go/math/all_test.go +++ b/libgo/go/math/all_test.go @@ -175,6 +175,7 @@ var cosLarge = []float64{ -2.51772931436786954751e-01, -7.3924135157173099849e-01, } + var cosh = []float64{ 7.2668796942212842775517446e+01, 1.1479413465659254502011135e+03, @@ -1527,6 +1528,7 @@ var vflog1pSC = []float64{ 0, Inf(1), NaN(), + 4503599627370496.5, // Issue #29488 } var log1pSC = []float64{ NaN(), @@ -1536,6 +1538,7 @@ var log1pSC = []float64{ 0, Inf(1), NaN(), + 36.04365338911715, // Issue #29488 } var vfmodfSC = []float64{ @@ -3026,6 +3029,41 @@ func TestLargeTan(t *testing.T) { } } +// Check that trigReduce matches the standard reduction results for input values +// below reduceThreshold. +func TestTrigReduce(t *testing.T) { + inputs := make([]float64, len(vf)) + // all of the standard inputs + copy(inputs, vf) + // all of the large inputs + large := float64(100000 * Pi) + for _, v := range vf { + inputs = append(inputs, v+large) + } + // Also test some special inputs, Pi and right below the reduceThreshold + inputs = append(inputs, Pi, Nextafter(ReduceThreshold, 0)) + for _, x := range inputs { + // reduce the value to compare + j, z := TrigReduce(x) + xred := float64(j)*(Pi/4) + z + + if f, fred := Sin(x), Sin(xred); !close(f, fred) { + t.Errorf("Sin(trigReduce(%g)) != Sin(%g), got %g, want %g", x, x, fred, f) + } + if f, fred := Cos(x), Cos(xred); !close(f, fred) { + t.Errorf("Cos(trigReduce(%g)) != Cos(%g), got %g, want %g", x, x, fred, f) + } + if f, fred := Tan(x), Tan(xred); !close(f, fred) { + t.Errorf(" Tan(trigReduce(%g)) != Tan(%g), got %g, want %g", x, x, fred, f) + } + f, g := Sincos(x) + fred, gred := Sincos(xred) + if !close(f, fred) || !close(g, gred) { + t.Errorf(" Sincos(trigReduce(%g)) != Sincos(%g), got %g, %g, want %g, %g", x, x, fred, gred, f, g) + } + } +} + // Check that math constants are accepted by compiler // and have right value (assumes strconv.ParseFloat works). // https://golang.org/issue/201 @@ -3635,3 +3673,41 @@ func BenchmarkYn(b *testing.B) { } GlobalF = x } + +func BenchmarkFloat64bits(b *testing.B) { + y := uint64(0) + for i := 0; i < b.N; i++ { + y = Float64bits(roundNeg) + } + GlobalI = int(y) +} + +var roundUint64 = uint64(5) + +func BenchmarkFloat64frombits(b *testing.B) { + x := 0.0 + for i := 0; i < b.N; i++ { + x = Float64frombits(roundUint64) + } + GlobalF = x +} + +var roundFloat32 = float32(-2.5) + +func BenchmarkFloat32bits(b *testing.B) { + y := uint32(0) + for i := 0; i < b.N; i++ { + y = Float32bits(roundFloat32) + } + GlobalI = int(y) +} + +var roundUint32 = uint32(5) + +func BenchmarkFloat32frombits(b *testing.B) { + x := float32(0.0) + for i := 0; i < b.N; i++ { + x = Float32frombits(roundUint32) + } + GlobalF = float64(x) +} diff --git a/libgo/go/math/big/arith.go b/libgo/go/math/big/arith.go index ad35240..f9db911 100644 --- a/libgo/go/math/big/arith.go +++ b/libgo/go/math/big/arith.go @@ -82,7 +82,7 @@ func nlz(x Word) uint { return uint(bits.LeadingZeros(uint(x))) } -// q = (u1<<_W + u0 - r)/y +// q = (u1<<_W + u0 - r)/v // Adapted from Warren, Hacker's Delight, p. 152. func divWW_g(u1, u0, v Word) (q, r Word) { if u1 >= v { diff --git a/libgo/go/math/big/float.go b/libgo/go/math/big/float.go index 55b93c8..b3c3295 100644 --- a/libgo/go/math/big/float.go +++ b/libgo/go/math/big/float.go @@ -43,7 +43,7 @@ const debugFloat = false // enable for debugging // precision of the argument with the largest precision value before any // rounding takes place, and the rounding mode remains unchanged. Thus, // uninitialized Floats provided as result arguments will have their -// precision set to a reasonable value determined by the operands and +// precision set to a reasonable value determined by the operands, and // their mode is the zero value for RoundingMode (ToNearestEven). // // By setting the desired precision to 24 or 53 and using matching rounding @@ -56,6 +56,12 @@ const debugFloat = false // enable for debugging // The zero (uninitialized) value for a Float is ready to use and represents // the number +0.0 exactly, with precision 0 and rounding mode ToNearestEven. // +// Operations always take pointer arguments (*Float) rather +// than Float values, and each unique Float value requires +// its own unique *Float pointer. To "copy" a Float value, +// an existing (or newly allocated) Float must be set to +// a new value using the Float.Set method; shallow copies +// of Floats are not supported and may lead to errors. type Float struct { prec uint32 mode RoundingMode @@ -293,7 +299,7 @@ func (z *Float) setExpAndRound(exp int64, sbit uint) { z.round(sbit) } -// SetMantExp sets z to mant × 2**exp and and returns z. +// SetMantExp sets z to mant × 2**exp and returns z. // The result z has the same precision and rounding mode // as mant. SetMantExp is an inverse of MantExp but does // not require 0.5 <= |mant| < 1.0. Specifically: @@ -321,7 +327,7 @@ func (z *Float) SetMantExp(mant *Float, exp int) *Float { return z } -// Signbit returns true if x is negative or negative zero. +// Signbit reports whether x is negative or negative zero. func (x *Float) Signbit() bool { return x.neg } diff --git a/libgo/go/math/big/int.go b/libgo/go/math/big/int.go index 47a288a..dab9a5c 100644 --- a/libgo/go/math/big/int.go +++ b/libgo/go/math/big/int.go @@ -15,6 +15,13 @@ import ( // An Int represents a signed multi-precision integer. // The zero value for an Int represents the value 0. +// +// Operations always take pointer arguments (*Int) rather +// than Int values, and each unique Int value requires +// its own unique *Int pointer. To "copy" an Int value, +// an existing (or newly allocated) Int must be set to +// a new value using the Int.Set method; shallow copies +// of Ints are not supported and may lead to errors. type Int struct { neg bool // sign abs nat // absolute value of the integer diff --git a/libgo/go/math/big/int_test.go b/libgo/go/math/big/int_test.go index 9930ed0..7ef2b39 100644 --- a/libgo/go/math/big/int_test.go +++ b/libgo/go/math/big/int_test.go @@ -1727,3 +1727,29 @@ func BenchmarkIntSqr(b *testing.B) { }) } } + +func benchmarkDiv(b *testing.B, aSize, bSize int) { + var r = rand.New(rand.NewSource(1234)) + aa := randInt(r, uint(aSize)) + bb := randInt(r, uint(bSize)) + if aa.Cmp(bb) < 0 { + aa, bb = bb, aa + } + x := new(Int) + y := new(Int) + + b.ResetTimer() + for i := 0; i < b.N; i++ { + x.DivMod(aa, bb, y) + } +} + +func BenchmarkDiv(b *testing.B) { + min, max, step := 10, 100000, 10 + for i := min; i <= max; i *= step { + j := 2 * i + b.Run(fmt.Sprintf("%d/%d", j, i), func(b *testing.B) { + benchmarkDiv(b, j, i) + }) + } +} diff --git a/libgo/go/math/big/nat.go b/libgo/go/math/big/nat.go index a6f79ed..1e4a3b0 100644 --- a/libgo/go/math/big/nat.go +++ b/libgo/go/math/big/nat.go @@ -58,6 +58,10 @@ func (z nat) make(n int) nat { if n <= cap(z) { return z[:n] // reuse z } + if n == 1 { + // Most nats start small and stay that way; don't over-allocate. + return make(nat, 1) + } // Choosing a good value for e has significant performance impact // because it increases the chance that a value can be reused. const e = 4 // extra capacity @@ -680,43 +684,36 @@ func putNat(x *nat) { var natPool sync.Pool -// q = (uIn-r)/v, with 0 <= r < y +// q = (uIn-r)/vIn, with 0 <= r < y // Uses z as storage for q, and u as storage for r if possible. // See Knuth, Volume 2, section 4.3.1, Algorithm D. // Preconditions: -// len(v) >= 2 -// len(uIn) >= len(v) -func (z nat) divLarge(u, uIn, v nat) (q, r nat) { - n := len(v) +// len(vIn) >= 2 +// len(uIn) >= len(vIn) +// u must not alias z +func (z nat) divLarge(u, uIn, vIn nat) (q, r nat) { + n := len(vIn) m := len(uIn) - n - // determine if z can be reused - // TODO(gri) should find a better solution - this if statement - // is very costly (see e.g. time pidigits -s -n 10000) - if alias(z, u) || alias(z, uIn) || alias(z, v) { - z = nil // z is an alias for u or uIn or v - cannot reuse + // D1. + shift := nlz(vIn[n-1]) + // do not modify vIn, it may be used by another goroutine simultaneously + vp := getNat(n) + v := *vp + shlVU(v, vIn, shift) + + // u may safely alias uIn or vIn, the value of uIn is used to set u and vIn was already used + u = u.make(len(uIn) + 1) + u[len(uIn)] = shlVU(u[0:len(uIn)], uIn, shift) + + // z may safely alias uIn or vIn, both values were used already + if alias(z, u) { + z = nil // z is an alias for u - cannot reuse } q = z.make(m + 1) qhatvp := getNat(n + 1) qhatv := *qhatvp - if alias(u, uIn) || alias(u, v) { - u = nil // u is an alias for uIn or v - cannot reuse - } - u = u.make(len(uIn) + 1) - u.clear() // TODO(gri) no need to clear if we allocated a new u - - // D1. - var v1p *nat - shift := nlz(v[n-1]) - if shift > 0 { - // do not modify v, it may be used by another goroutine simultaneously - v1p = getNat(n) - v1 := *v1p - shlVU(v1, v, shift) - v = v1 - } - u[len(uIn)] = shlVU(u[0:len(uIn)], uIn, shift) // D2. vn1 := v[n-1] @@ -756,9 +753,8 @@ func (z nat) divLarge(u, uIn, v nat) (q, r nat) { q[j] = qhat } - if v1p != nil { - putNat(v1p) - } + + putNat(vp) putNat(qhatvp) q = q.norm() diff --git a/libgo/go/math/big/prime.go b/libgo/go/math/big/prime.go index 4c2c152..d9a5f1e 100644 --- a/libgo/go/math/big/prime.go +++ b/libgo/go/math/big/prime.go @@ -51,7 +51,7 @@ func (x *Int) ProbablyPrime(n int) bool { } if w&1 == 0 { - return false // n is even + return false // x is even } const primesA = 3 * 5 * 7 * 11 * 13 * 17 * 19 * 23 * 37 diff --git a/libgo/go/math/big/rat.go b/libgo/go/math/big/rat.go index 46d58fc..5d0800c 100644 --- a/libgo/go/math/big/rat.go +++ b/libgo/go/math/big/rat.go @@ -13,6 +13,13 @@ import ( // A Rat represents a quotient a/b of arbitrary precision. // The zero value for a Rat represents the value 0. +// +// Operations always take pointer arguments (*Rat) rather +// than Rat values, and each unique Rat value requires +// its own unique *Rat pointer. To "copy" a Rat value, +// an existing (or newly allocated) Rat must be set to +// a new value using the Rat.Set method; shallow copies +// of Rats are not supported and may lead to errors. type Rat struct { // To make zero values for Rat work w/o initialization, // a zero value of b (len(b) == 0) acts like b == 1. diff --git a/libgo/go/math/big/sqrt.go b/libgo/go/math/big/sqrt.go index b989649..53403aa 100644 --- a/libgo/go/math/big/sqrt.go +++ b/libgo/go/math/big/sqrt.go @@ -7,8 +7,6 @@ package big import "math" var ( - half = NewFloat(0.5) - two = NewFloat(2.0) three = NewFloat(3.0) ) @@ -57,9 +55,9 @@ func (z *Float) Sqrt(x *Float) *Float { case 0: // nothing to do case 1: - z.Mul(two, z) + z.exp++ case -1: - z.Mul(half, z) + z.exp-- } // 0.25 <= z < 2.0 @@ -96,7 +94,7 @@ func (z *Float) sqrtDirect(x *Float) { u.prec = t.prec u.Mul(t, t) // u = t² u.Add(u, x) // = t² + x - u.Mul(half, u) // = ½(t² + x) + u.exp-- // = ½(t² + x) return t.Quo(u, t) // = ½(t² + x)/t } @@ -133,11 +131,13 @@ func (z *Float) sqrtInverse(x *Float) { ng := func(t *Float) *Float { u.prec = t.prec v.prec = t.prec - u.Mul(t, t) // u = t² - u.Mul(x, u) // = xt² - v.Sub(three, u) // v = 3 - xt² - u.Mul(t, v) // u = t(3 - xt²) - return t.Mul(half, u) // = ½t(3 - xt²) + u.Mul(t, t) // u = t² + u.Mul(x, u) // = xt² + v.Sub(three, u) // v = 3 - xt² + u.Mul(t, v) // u = t(3 - xt²) + u.exp-- // = ½t(3 - xt²) + return t.Set(u) + } xf, _ := x.Float64() diff --git a/libgo/go/math/bits/bits.go b/libgo/go/math/bits/bits.go index 989baac..306fa76 100644 --- a/libgo/go/math/bits/bits.go +++ b/libgo/go/math/bits/bits.go @@ -8,6 +8,8 @@ // functions for the predeclared unsigned integer types. package bits +import _ "unsafe" // for go:linkname + const uintSize = 32 << (^uint(0) >> 32 & 1) // 32 or 64 // UintSize is the size of a uint in bits. @@ -63,7 +65,7 @@ func TrailingZeros8(x uint8) int { } // TrailingZeros16 returns the number of trailing zero bits in x; the result is 16 for x == 0. -func TrailingZeros16(x uint16) (n int) { +func TrailingZeros16(x uint16) int { if x == 0 { return 16 } @@ -328,3 +330,206 @@ func Len64(x uint64) (n int) { } return n + int(len8tab[x]) } + +// --- Add with carry --- + +// Add returns the sum with carry of x, y and carry: sum = x + y + carry. +// The carry input must be 0 or 1; otherwise the behavior is undefined. +// The carryOut output is guaranteed to be 0 or 1. +func Add(x, y, carry uint) (sum, carryOut uint) { + yc := y + carry + sum = x + yc + if sum < x || yc < y { + carryOut = 1 + } + return +} + +// Add32 returns the sum with carry of x, y and carry: sum = x + y + carry. +// The carry input must be 0 or 1; otherwise the behavior is undefined. +// The carryOut output is guaranteed to be 0 or 1. +func Add32(x, y, carry uint32) (sum, carryOut uint32) { + yc := y + carry + sum = x + yc + if sum < x || yc < y { + carryOut = 1 + } + return +} + +// Add64 returns the sum with carry of x, y and carry: sum = x + y + carry. +// The carry input must be 0 or 1; otherwise the behavior is undefined. +// The carryOut output is guaranteed to be 0 or 1. +func Add64(x, y, carry uint64) (sum, carryOut uint64) { + yc := y + carry + sum = x + yc + if sum < x || yc < y { + carryOut = 1 + } + return +} + +// --- Subtract with borrow --- + +// Sub returns the difference of x, y and borrow: diff = x - y - borrow. +// The borrow input must be 0 or 1; otherwise the behavior is undefined. +// The borrowOut output is guaranteed to be 0 or 1. +func Sub(x, y, borrow uint) (diff, borrowOut uint) { + yb := y + borrow + diff = x - yb + if diff > x || yb < y { + borrowOut = 1 + } + return +} + +// Sub32 returns the difference of x, y and borrow, diff = x - y - borrow. +// The borrow input must be 0 or 1; otherwise the behavior is undefined. +// The borrowOut output is guaranteed to be 0 or 1. +func Sub32(x, y, borrow uint32) (diff, borrowOut uint32) { + yb := y + borrow + diff = x - yb + if diff > x || yb < y { + borrowOut = 1 + } + return +} + +// Sub64 returns the difference of x, y and borrow: diff = x - y - borrow. +// The borrow input must be 0 or 1; otherwise the behavior is undefined. +// The borrowOut output is guaranteed to be 0 or 1. +func Sub64(x, y, borrow uint64) (diff, borrowOut uint64) { + yb := y + borrow + diff = x - yb + if diff > x || yb < y { + borrowOut = 1 + } + return +} + +// --- Full-width multiply --- + +// Mul returns the full-width product of x and y: (hi, lo) = x * y +// with the product bits' upper half returned in hi and the lower +// half returned in lo. +func Mul(x, y uint) (hi, lo uint) { + if UintSize == 32 { + h, l := Mul32(uint32(x), uint32(y)) + return uint(h), uint(l) + } + h, l := Mul64(uint64(x), uint64(y)) + return uint(h), uint(l) +} + +// Mul32 returns the 64-bit product of x and y: (hi, lo) = x * y +// with the product bits' upper half returned in hi and the lower +// half returned in lo. +func Mul32(x, y uint32) (hi, lo uint32) { + tmp := uint64(x) * uint64(y) + hi, lo = uint32(tmp>>32), uint32(tmp) + return +} + +// Mul64 returns the 128-bit product of x and y: (hi, lo) = x * y +// with the product bits' upper half returned in hi and the lower +// half returned in lo. +func Mul64(x, y uint64) (hi, lo uint64) { + const mask32 = 1<<32 - 1 + x0 := x & mask32 + x1 := x >> 32 + y0 := y & mask32 + y1 := y >> 32 + w0 := x0 * y0 + t := x1*y0 + w0>>32 + w1 := t & mask32 + w2 := t >> 32 + w1 += x0 * y1 + hi = x1*y1 + w2 + w1>>32 + lo = x * y + return +} + +// --- Full-width divide --- + +// Div returns the quotient and remainder of (hi, lo) divided by y: +// quo = (hi, lo)/y, rem = (hi, lo)%y with the dividend bits' upper +// half in parameter hi and the lower half in parameter lo. +// Div panics for y == 0 (division by zero) or y <= hi (quotient overflow). +func Div(hi, lo, y uint) (quo, rem uint) { + if UintSize == 32 { + q, r := Div32(uint32(hi), uint32(lo), uint32(y)) + return uint(q), uint(r) + } + q, r := Div64(uint64(hi), uint64(lo), uint64(y)) + return uint(q), uint(r) +} + +// Div32 returns the quotient and remainder of (hi, lo) divided by y: +// quo = (hi, lo)/y, rem = (hi, lo)%y with the dividend bits' upper +// half in parameter hi and the lower half in parameter lo. +// Div32 panics for y == 0 (division by zero) or y <= hi (quotient overflow). +func Div32(hi, lo, y uint32) (quo, rem uint32) { + if y != 0 && y <= hi { + panic(getOverflowError()) + } + z := uint64(hi)<<32 | uint64(lo) + quo, rem = uint32(z/uint64(y)), uint32(z%uint64(y)) + return +} + +// Div64 returns the quotient and remainder of (hi, lo) divided by y: +// quo = (hi, lo)/y, rem = (hi, lo)%y with the dividend bits' upper +// half in parameter hi and the lower half in parameter lo. +// Div64 panics for y == 0 (division by zero) or y <= hi (quotient overflow). +func Div64(hi, lo, y uint64) (quo, rem uint64) { + const ( + two32 = 1 << 32 + mask32 = two32 - 1 + ) + if y == 0 { + panic(getDivideError()) + } + if y <= hi { + panic(getOverflowError()) + } + + s := uint(LeadingZeros64(y)) + y <<= s + + yn1 := y >> 32 + yn0 := y & mask32 + un32 := hi<>(64-s) + un10 := lo << s + un1 := un10 >> 32 + un0 := un10 & mask32 + q1 := un32 / yn1 + rhat := un32 - q1*yn1 + + for q1 >= two32 || q1*yn0 > two32*rhat+un1 { + q1-- + rhat += yn1 + if rhat >= two32 { + break + } + } + + un21 := un32*two32 + un1 - q1*y + q0 := un21 / yn1 + rhat = un21 - q0*yn1 + + for q0 >= two32 || q0*yn0 > two32*rhat+un0 { + q0-- + rhat += yn1 + if rhat >= two32 { + break + } + } + + return q1*two32 + q0, (un21*two32 + un0 - q0*y) >> s +} + +//go:linkname getOverflowError runtime.getOverflowError +func getOverflowError() error + +//go:linkname getDivideError runtime.getDivideError +func getDivideError() error diff --git a/libgo/go/math/bits/bits_test.go b/libgo/go/math/bits/bits_test.go index 5c34f6d..1ec5107 100644 --- a/libgo/go/math/bits/bits_test.go +++ b/libgo/go/math/bits/bits_test.go @@ -6,6 +6,7 @@ package bits_test import ( . "math/bits" + "runtime" "testing" "unsafe" ) @@ -705,6 +706,385 @@ func TestLen(t *testing.T) { } } +const ( + _M = 1<> 52) - 1023) + // correction term if k > 0 { c = 1.0 - (u - x) } else { - c = x - (u - 1.0) // correction term - c /= u + c = x - (u - 1.0) } + c /= u } else { u = x iu = Float64bits(u) diff --git a/libgo/go/math/mod.go b/libgo/go/math/mod.go index 0b208f4..436788f 100644 --- a/libgo/go/math/mod.go +++ b/libgo/go/math/mod.go @@ -30,16 +30,12 @@ func mod(x, y float64) float64 { if y == 0 || IsInf(x, 0) || IsNaN(x) || IsNaN(y) { return NaN() } - if y < 0 { - y = -y - } + y = Abs(y) yfr, yexp := Frexp(y) - sign := false r := x if x < 0 { r = -x - sign = true } for r >= y { @@ -49,7 +45,7 @@ func mod(x, y float64) float64 { } r = r - Ldexp(y, rexp-yexp) } - if sign { + if x < 0 { r = -r } return r diff --git a/libgo/go/math/pow.go b/libgo/go/math/pow.go index 314ff90..b812179 100644 --- a/libgo/go/math/pow.go +++ b/libgo/go/math/pow.go @@ -88,13 +88,7 @@ func pow(x, y float64) float64 { return 1 / Sqrt(x) } - absy := y - flip := false - if absy < 0 { - absy = -absy - flip = true - } - yi, yf := Modf(absy) + yi, yf := Modf(Abs(y)) if yf != 0 && x < 0 { return NaN() } @@ -152,9 +146,9 @@ func pow(x, y float64) float64 { } // ans = a1*2**ae - // if flip { ans = 1 / ans } + // if y < 0 { ans = 1 / ans } // but in the opposite order - if flip { + if y < 0 { a1 = 1 / a1 ae = -ae } diff --git a/libgo/go/math/signbit.go b/libgo/go/math/signbit.go index 670cc1a..f6e61d6 100644 --- a/libgo/go/math/signbit.go +++ b/libgo/go/math/signbit.go @@ -4,7 +4,7 @@ package math -// Signbit returns true if x is negative or negative zero. +// Signbit reports whether x is negative or negative zero. func Signbit(x float64) bool { return Float64bits(x)&(1<<63) != 0 } diff --git a/libgo/go/math/sin.go b/libgo/go/math/sin.go index e3166e2..871e8c6 100644 --- a/libgo/go/math/sin.go +++ b/libgo/go/math/sin.go @@ -124,10 +124,9 @@ func Cos(x float64) float64 { 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 + PI4A = 7.85398125648498535156E-1 // 0x3fe921fb40000000, Pi/4 split into three parts + PI4B = 3.77489470793079817668E-8 // 0x3e64442d00000000, + PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170, ) // special cases switch { @@ -139,15 +138,23 @@ func cos(x float64) float64 { sign := false x = Abs(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++ - y++ + var j uint64 + var y, z float64 + if x >= reduceThreshold { + j, z = trigReduce(x) + } else { + j = uint64(x * (4 / Pi)) // 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++ + y++ + } + j &= 7 // octant modulo 2Pi radians (360 degrees) + z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic } - j &= 7 // octant modulo 2Pi radians (360 degrees) + if j > 3 { j -= 4 sign = !sign @@ -156,7 +163,6 @@ func cos(x float64) float64 { 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]) @@ -185,10 +191,9 @@ func Sin(x float64) float64 { 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 + PI4A = 7.85398125648498535156E-1 // 0x3fe921fb40000000, Pi/4 split into three parts + PI4B = 3.77489470793079817668E-8 // 0x3e64442d00000000, + PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170, ) // special cases switch { @@ -205,22 +210,27 @@ func sin(x float64) float64 { 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++ - y++ + var j uint64 + var y, z float64 + if x >= reduceThreshold { + j, z = trigReduce(x) + } else { + j = uint64(x * (4 / Pi)) // 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++ + y++ + } + j &= 7 // octant modulo 2Pi radians (360 degrees) + z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic } - 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]) diff --git a/libgo/go/math/sincos.go b/libgo/go/math/sincos.go index 65f3790..c002db6 100644 --- a/libgo/go/math/sincos.go +++ b/libgo/go/math/sincos.go @@ -2,8 +2,6 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -// -build !386 - package math // Coefficients _sin[] and _cos[] are found in pkg/math/sin.go. @@ -16,10 +14,9 @@ package math // Sincos(NaN) = NaN, NaN func Sincos(x float64) (sin, cos 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 + PI4A = 7.85398125648498535156E-1 // 0x3fe921fb40000000, Pi/4 split into three parts + PI4B = 3.77489470793079817668E-8 // 0x3e64442d00000000, + PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170, ) // special cases switch { @@ -36,14 +33,21 @@ func Sincos(x float64) (sin, cos float64) { sinSign = 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 + var j uint64 + var y, z float64 + if x >= reduceThreshold { + j, z = trigReduce(x) + } else { + j = uint64(x * (4 / Pi)) // 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 - if j&1 == 1 { // map zeros to origin - j++ - y++ + if j&1 == 1 { // map zeros to origin + j++ + y++ + } + j &= 7 // octant modulo 2Pi radians (360 degrees) + z = ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic } - j &= 7 // octant modulo 2Pi radians (360 degrees) if j > 3 { // reflect in x axis j -= 4 sinSign, cosSign = !sinSign, !cosSign @@ -52,7 +56,6 @@ func Sincos(x float64) (sin, cos float64) { cosSign = !cosSign } - z := ((x - y*PI4A) - y*PI4B) - y*PI4C // Extended precision modular arithmetic zz := z * z cos = 1.0 - 0.5*zz + zz*zz*((((((_cos[0]*zz)+_cos[1])*zz+_cos[2])*zz+_cos[3])*zz+_cos[4])*zz+_cos[5]) sin = z + z*zz*((((((_sin[0]*zz)+_sin[1])*zz+_sin[2])*zz+_sin[3])*zz+_sin[4])*zz+_sin[5]) diff --git a/libgo/go/math/sincos_386.go b/libgo/go/math/sincos_386.go deleted file mode 100644 index d5c2457..0000000 --- a/libgo/go/math/sincos_386.go +++ /dev/null @@ -1,15 +0,0 @@ -// Copyright 2017 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. - -// +build ignore - -package math - -// Sincos returns Sin(x), Cos(x). -// -// Special cases are: -// Sincos(±0) = ±0, 1 -// Sincos(±Inf) = NaN, NaN -// Sincos(NaN) = NaN, NaN -func Sincos(x float64) (sin, cos float64) diff --git a/libgo/go/math/sinh.go b/libgo/go/math/sinh.go index 9dff71e..993424d 100644 --- a/libgo/go/math/sinh.go +++ b/libgo/go/math/sinh.go @@ -41,7 +41,7 @@ func Sinh(x float64) float64 { } var temp float64 - switch true { + switch { case x > 21: temp = Exp(x) * 0.5 diff --git a/libgo/go/math/tan.go b/libgo/go/math/tan.go index f5230d3..929a0a9 100644 --- a/libgo/go/math/tan.go +++ b/libgo/go/math/tan.go @@ -89,10 +89,9 @@ func Tan(x float64) float64 { func tan(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 + PI4A = 7.85398125648498535156E-1 // 0x3fe921fb40000000, Pi/4 split into three parts + PI4B = 3.77489470793079817668E-8 // 0x3e64442d00000000, + PI4C = 2.69515142907905952645E-15 // 0x3ce8469898cc5170, ) // special cases switch { @@ -108,17 +107,22 @@ func tan(x float64) float64 { x = -x sign = true } + var j uint64 + var y, z float64 + if x >= reduceThreshold { + j, z = trigReduce(x) + } else { + j = uint64(x * (4 / Pi)) // 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 - 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 and singularities to origin */ + if j&1 == 1 { + j++ + y++ + } - /* map zeros and singularities to origin */ - if j&1 == 1 { - j++ - y++ + z = ((x - y*PI4A) - y*PI4B) - y*PI4C } - - z := ((x - y*PI4A) - y*PI4B) - y*PI4C zz := z * z if zz > 1e-14 { diff --git a/libgo/go/math/trig_reduce.go b/libgo/go/math/trig_reduce.go new file mode 100644 index 0000000..6f8eaba --- /dev/null +++ b/libgo/go/math/trig_reduce.go @@ -0,0 +1,94 @@ +// Copyright 2018 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 ( + "math/bits" +) + +// reduceThreshold is the maximum value where the reduction using Pi/4 +// in 3 float64 parts still gives accurate results. Above this +// threshold Payne-Hanek range reduction must be used. +const reduceThreshold = (1 << 52) / (4 / Pi) + +// trigReduce implements Payne-Hanek range reduction by Pi/4 +// for x > 0. It returns the integer part mod 8 (j) and +// the fractional part (z) of x / (Pi/4). +// The implementation is based on: +// "ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit" +// K. C. Ng et al, March 24, 1992 +// The simulated multi-precision calculation of x*B uses 64-bit integer arithmetic. +func trigReduce(x float64) (j uint64, z float64) { + const PI4 = Pi / 4 + if x < PI4 { + return 0, x + } + // Extract out the integer and exponent such that, + // x = ix * 2 ** exp. + ix := Float64bits(x) + exp := int(ix>>shift&mask) - bias - shift + ix &^= mask << shift + ix |= 1 << shift + // Use the exponent to extract the 3 appropriate uint64 digits from mPi4, + // B ~ (z0, z1, z2), such that the product leading digit has the exponent -61. + // Note, exp >= -53 since x >= PI4 and exp < 971 for maximum float64. + digit, bitshift := uint(exp+61)/64, uint(exp+61)%64 + z0 := (mPi4[digit] << bitshift) | (mPi4[digit+1] >> (64 - bitshift)) + z1 := (mPi4[digit+1] << bitshift) | (mPi4[digit+2] >> (64 - bitshift)) + z2 := (mPi4[digit+2] << bitshift) | (mPi4[digit+3] >> (64 - bitshift)) + // Multiply mantissa by the digits and extract the upper two digits (hi, lo). + z2hi, _ := bits.Mul64(z2, ix) + z1hi, z1lo := bits.Mul64(z1, ix) + z0lo := z0 * ix + lo, c := bits.Add64(z1lo, z2hi, 0) + hi, _ := bits.Add64(z0lo, z1hi, c) + // The top 3 bits are j. + j = hi >> 61 + // Extract the fraction and find its magnitude. + hi = hi<<3 | lo>>61 + lz := uint(bits.LeadingZeros64(hi)) + e := uint64(bias - (lz + 1)) + // Clear implicit mantissa bit and shift into place. + hi = (hi << (lz + 1)) | (lo >> (64 - (lz + 1))) + hi >>= 64 - shift + // Include the exponent and convert to a float. + hi |= e << shift + z = Float64frombits(hi) + // Map zeros to origin. + if j&1 == 1 { + j++ + j &= 7 + z-- + } + // Multiply the fractional part by pi/4. + return j, z * PI4 +} + +// mPi4 is the binary digits of 4/pi as a uint64 array, +// that is, 4/pi = Sum mPi4[i]*2^(-64*i) +// 19 64-bit digits and the leading one bit give 1217 bits +// of precision to handle the largest possible float64 exponent. +var mPi4 = [...]uint64{ + 0x0000000000000001, + 0x45f306dc9c882a53, + 0xf84eafa3ea69bb81, + 0xb6c52b3278872083, + 0xfca2c757bd778ac3, + 0x6e48dc74849ba5c0, + 0x0c925dd413a32439, + 0xfc3bd63962534e7d, + 0xd1046bea5d768909, + 0xd338e04d68befc82, + 0x7323ac7306a673e9, + 0x3908bf177bf25076, + 0x3ff12fffbc0b301f, + 0xde5e2316b414da3e, + 0xda6cfd9e4f96136e, + 0x9e8c7ecd3cbfd45a, + 0xea4f758fd7cbe2f6, + 0x7a0e73ef14a525d4, + 0xd7f6bf623f1aba10, + 0xac06608df8f6d757, +} diff --git a/libgo/go/math/unsafe.go b/libgo/go/math/unsafe.go index 5ae6742..e59f50c 100644 --- a/libgo/go/math/unsafe.go +++ b/libgo/go/math/unsafe.go @@ -6,16 +6,24 @@ package math import "unsafe" -// Float32bits returns the IEEE 754 binary representation of f. +// Float32bits returns the IEEE 754 binary representation of f, +// with the sign bit of f and the result in the same bit position. +// Float32bits(Float32frombits(x)) == x. func Float32bits(f float32) uint32 { return *(*uint32)(unsafe.Pointer(&f)) } -// Float32frombits returns the floating point number corresponding -// to the IEEE 754 binary representation b. +// Float32frombits returns the floating-point number corresponding +// to the IEEE 754 binary representation b, with the sign bit of b +// and the result in the same bit position. +// Float32frombits(Float32bits(x)) == x. func Float32frombits(b uint32) float32 { return *(*float32)(unsafe.Pointer(&b)) } -// Float64bits returns the IEEE 754 binary representation of f. +// Float64bits returns the IEEE 754 binary representation of f, +// with the sign bit of f and the result in the same bit position, +// and Float64bits(Float64frombits(x)) == x. func Float64bits(f float64) uint64 { return *(*uint64)(unsafe.Pointer(&f)) } -// Float64frombits returns the floating point number corresponding -// the IEEE 754 binary representation b. +// Float64frombits returns the floating-point number corresponding +// to the IEEE 754 binary representation b, with the sign bit of b +// and the result in the same bit position. +// Float64frombits(Float64bits(x)) == x. func Float64frombits(b uint64) float64 { return *(*float64)(unsafe.Pointer(&b)) } -- cgit v1.1