diff options
Diffstat (limited to 'libgo/go/math')
27 files changed, 1464 insertions, 264 deletions
diff --git a/libgo/go/math/all_test.go b/libgo/go/math/all_test.go index 1ac9d71..3aae037 100644 --- a/libgo/go/math/all_test.go +++ b/libgo/go/math/all_test.go @@ -827,6 +827,8 @@ var vfatan2SC = [][2]float64{ {+Pi, Inf(-1)}, {+Pi, 0}, {+Pi, Inf(1)}, + {1.0, Inf(1)}, + {-1.0, Inf(1)}, {+Pi, NaN()}, {Inf(1), Inf(-1)}, {Inf(1), -Pi}, @@ -864,6 +866,8 @@ var atan2SC = []float64{ Pi, // atan2(+Pi, -Inf) Pi / 2, // atan2(+Pi, +0) 0, // atan2(+Pi, +Inf) + 0, // atan2(+1, +Inf) + Copysign(0, -1), // atan2(-1, +Inf) NaN(), // atan2(+Pi, NaN) 3 * Pi / 4, // atan2(+Inf, -Inf) Pi / 2, // atan2(+Inf, -Pi) diff --git a/libgo/go/math/arith_s390x.go b/libgo/go/math/arith_s390x.go index 8d1fa6a..c28651e 100644 --- a/libgo/go/math/arith_s390x.go +++ b/libgo/go/math/arith_s390x.go @@ -6,6 +6,8 @@ package math +import "internal/cpu" + func log10TrampolineSetup(x float64) float64 func log10Asm(x float64) float64 @@ -72,8 +74,6 @@ func expm1Asm(x float64) float64 func powTrampolineSetup(x, y float64) float64 func powAsm(x, y float64) float64 -// hasVectorFacility reports whether the machine has the z/Architecture +// hasVX reports whether the machine has the z/Architecture // vector facility installed and enabled. -func hasVectorFacility() bool - -var hasVX = hasVectorFacility() +var hasVX = cpu.S390X.HasVX diff --git a/libgo/go/math/big/arith_decl.go b/libgo/go/math/big/arith_decl.go index 0a139f1..61df0df 100644 --- a/libgo/go/math/big/arith_decl.go +++ b/libgo/go/math/big/arith_decl.go @@ -3,7 +3,7 @@ // license that can be found in the LICENSE file. // +build ignore -// +build !math_big_pure_go,!riscv64 +// +build !math_big_pure_go package big diff --git a/libgo/go/math/big/arith_decl_pure.go b/libgo/go/math/big/arith_decl_pure.go index 8853eb6..ee8f922 100644 --- a/libgo/go/math/big/arith_decl_pure.go +++ b/libgo/go/math/big/arith_decl_pure.go @@ -2,7 +2,7 @@ // Use of this source code is governed by a BSD-style // license that can be found in the LICENSE file. -// -build math_big_pure_go riscv64 +// -build math_big_pure_go package big diff --git a/libgo/go/math/big/arith_decl_s390x.go b/libgo/go/math/big/arith_decl_s390x.go index 3936d3f..712f115 100644 --- a/libgo/go/math/big/arith_decl_s390x.go +++ b/libgo/go/math/big/arith_decl_s390x.go @@ -7,18 +7,13 @@ package big +import "internal/cpu" + func addVV_check(z, x, y []Word) (c Word) func addVV_vec(z, x, y []Word) (c Word) func addVV_novec(z, x, y []Word) (c Word) func subVV_check(z, x, y []Word) (c Word) func subVV_vec(z, x, y []Word) (c Word) func subVV_novec(z, x, y []Word) (c Word) -func addVW_check(z, x []Word, y Word) (c Word) -func addVW_vec(z, x []Word, y Word) (c Word) -func addVW_novec(z, x []Word, y Word) (c Word) -func subVW_check(z, x []Word, y Word) (c Word) -func subVW_vec(z, x []Word, y Word) (c Word) -func subVW_novec(z, x []Word, y Word) (c Word) -func hasVectorFacility() bool -var hasVX = hasVectorFacility() +var hasVX = cpu.S390X.HasVX diff --git a/libgo/go/math/big/arith_s390x_test.go b/libgo/go/math/big/arith_s390x_test.go index 0e8ac85..be1ca7d 100644 --- a/libgo/go/math/big/arith_s390x_test.go +++ b/libgo/go/math/big/arith_s390x_test.go @@ -31,15 +31,3 @@ func TestFunVVnovec(t *testing.T) { } } } - -func TestFunVWnovec(t *testing.T) { - if hasVX == true { - for _, a := range sumVW { - arg := a - testFunVW(t, "addVW_novec", addVW_novec, arg) - - arg = argVW{a.x, a.z, a.y, a.c} - testFunVW(t, "subVW_novec", subVW_novec, arg) - } - } -} diff --git a/libgo/go/math/big/float.go b/libgo/go/math/big/float.go index b3c3295..da964ee 100644 --- a/libgo/go/math/big/float.go +++ b/libgo/go/math/big/float.go @@ -224,7 +224,9 @@ func (x *Float) Mode() RoundingMode { return x.mode } -// Acc returns the accuracy of x produced by the most recent operation. +// Acc returns the accuracy of x produced by the most recent +// operation, unless explicitly documented otherwise by that +// operation. func (x *Float) Acc() Accuracy { return x.acc } diff --git a/libgo/go/math/big/floatconv.go b/libgo/go/math/big/floatconv.go index 95e32d3..57b7df3 100644 --- a/libgo/go/math/big/floatconv.go +++ b/libgo/go/math/big/floatconv.go @@ -290,7 +290,7 @@ func ParseFloat(s string, base int, prec uint, mode RoundingMode) (f *Float, b i return new(Float).SetPrec(prec).SetMode(mode).Parse(s, base) } -var _ fmt.Scanner = &floatZero // *Float must implement fmt.Scanner +var _ fmt.Scanner = (*Float)(nil) // *Float must implement fmt.Scanner // Scan is a support routine for fmt.Scanner; it sets z to the value of // the scanned number. It accepts formats whose verbs are supported by diff --git a/libgo/go/math/big/int.go b/libgo/go/math/big/int.go index 18f122e..65f3248 100644 --- a/libgo/go/math/big/int.go +++ b/libgo/go/math/big/int.go @@ -447,11 +447,26 @@ func (z *Int) SetBytes(buf []byte) *Int { } // Bytes returns the absolute value of x as a big-endian byte slice. +// +// To use a fixed length slice, or a preallocated one, use FillBytes. func (x *Int) Bytes() []byte { buf := make([]byte, len(x.abs)*_S) return buf[x.abs.bytes(buf):] } +// FillBytes sets buf to the absolute value of x, storing it as a zero-extended +// big-endian byte slice, and returns buf. +// +// If the absolute value of x doesn't fit in buf, FillBytes will panic. +func (x *Int) FillBytes(buf []byte) []byte { + // Clear whole buffer. (This gets optimized into a memclr.) + for i := range buf { + buf[i] = 0 + } + x.abs.bytes(buf) + return buf +} + // BitLen returns the length of the absolute value of x in bits. // The bit length of 0 is 0. func (x *Int) BitLen() int { @@ -465,8 +480,8 @@ func (x *Int) TrailingZeroBits() uint { } // Exp sets z = x**y mod |m| (i.e. the sign of m is ignored), and returns z. -// If m == nil or m == 0, z = x**y unless y <= 0 then z = 1. If m > 0, y < 0, -// and x and n are not relatively prime, z is unchanged and nil is returned. +// If m == nil or m == 0, z = x**y unless y <= 0 then z = 1. If m != 0, y < 0, +// and x and m are not relatively prime, z is unchanged and nil is returned. // // Modular exponentiation of inputs of a particular size is not a // cryptographically constant-time operation. diff --git a/libgo/go/math/big/int_test.go b/libgo/go/math/big/int_test.go index e3a1587..3c85573 100644 --- a/libgo/go/math/big/int_test.go +++ b/libgo/go/math/big/int_test.go @@ -1840,3 +1840,57 @@ func BenchmarkDiv(b *testing.B) { }) } } + +func TestFillBytes(t *testing.T) { + checkResult := func(t *testing.T, buf []byte, want *Int) { + t.Helper() + got := new(Int).SetBytes(buf) + if got.CmpAbs(want) != 0 { + t.Errorf("got 0x%x, want 0x%x: %x", got, want, buf) + } + } + panics := func(f func()) (panic bool) { + defer func() { panic = recover() != nil }() + f() + return + } + + for _, n := range []string{ + "0", + "1000", + "0xffffffff", + "-0xffffffff", + "0xffffffffffffffff", + "0x10000000000000000", + "0xabababababababababababababababababababababababababa", + "0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff", + } { + t.Run(n, func(t *testing.T) { + t.Logf(n) + x, ok := new(Int).SetString(n, 0) + if !ok { + panic("invalid test entry") + } + + // Perfectly sized buffer. + byteLen := (x.BitLen() + 7) / 8 + buf := make([]byte, byteLen) + checkResult(t, x.FillBytes(buf), x) + + // Way larger, checking all bytes get zeroed. + buf = make([]byte, 100) + for i := range buf { + buf[i] = 0xff + } + checkResult(t, x.FillBytes(buf), x) + + // Too small. + if byteLen > 0 { + buf = make([]byte, byteLen-1) + if !panics(func() { x.FillBytes(buf) }) { + t.Errorf("expected panic for small buffer and value %x", x) + } + } + }) + } +} diff --git a/libgo/go/math/big/link_test.go b/libgo/go/math/big/link_test.go new file mode 100644 index 0000000..2212bd4 --- /dev/null +++ b/libgo/go/math/big/link_test.go @@ -0,0 +1,63 @@ +// Copyright 2020 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 big + +import ( + "bytes" + "internal/testenv" + "io/ioutil" + "os/exec" + "path/filepath" + "testing" +) + +// Tests that the linker is able to remove references to Float, Rat, +// and Int if unused (notably, not used by init). +func TestLinkerGC(t *testing.T) { + if testing.Short() { + t.Skip("skipping in short mode") + } + t.Parallel() + tmp := t.TempDir() + goBin := testenv.GoToolPath(t) + goFile := filepath.Join(tmp, "x.go") + file := []byte(`package main +import _ "math/big" +func main() {} +`) + if err := ioutil.WriteFile(goFile, file, 0644); err != nil { + t.Fatal(err) + } + cmd := exec.Command(goBin, "build", "-o", "x.exe", "x.go") + cmd.Dir = tmp + if out, err := cmd.CombinedOutput(); err != nil { + t.Fatalf("compile: %v, %s", err, out) + } + + cmd = exec.Command(goBin, "tool", "nm", "x.exe") + cmd.Dir = tmp + nm, err := cmd.CombinedOutput() + if err != nil { + t.Fatalf("nm: %v, %s", err, nm) + } + const want = "runtime.(*Frames).Next" + if !bytes.Contains(nm, []byte(want)) { + // Test the test. + t.Errorf("expected symbol %q not found", want) + } + bad := []string{ + "math/big.(*Float)", + "math/big.(*Rat)", + "math/big.(*Int)", + } + for _, sym := range bad { + if bytes.Contains(nm, []byte(sym)) { + t.Errorf("unexpected symbol %q found", sym) + } + } + if t.Failed() { + t.Logf("Got: %s", nm) + } +} diff --git a/libgo/go/math/big/nat.go b/libgo/go/math/big/nat.go index 1b771ca..6a3989b 100644 --- a/libgo/go/math/big/nat.go +++ b/libgo/go/math/big/nat.go @@ -740,7 +740,8 @@ func (z nat) divLarge(u, uIn, vIn nat) (q, r nat) { // The remainder overwrites input u. // // Precondition: -// - len(q) >= len(u)-len(v) +// - q is large enough to hold the quotient u / v +// which has a maximum length of len(u)-len(v)+1. func (q nat) divBasic(u, v nat) { n := len(v) m := len(u) - n @@ -779,6 +780,8 @@ func (q nat) divBasic(u, v nat) { } // D4. + // Compute the remainder u - (q̂*v) << (_W*j). + // The subtraction may overflow if q̂ estimate was off by one. qhatv[n] = mulAddVWW(qhatv[0:n], v, qhat, 0) qhl := len(qhatv) if j+qhl > len(u) && qhatv[n] == 0 { @@ -787,7 +790,11 @@ func (q nat) divBasic(u, v nat) { c := subVV(u[j:j+qhl], u[j:], qhatv) if c != 0 { c := addVV(u[j:j+n], u[j:], v) - u[j+n] += c + // If n == qhl, the carry from subVV and the carry from addVV + // cancel out and don't affect u[j+n]. + if n < qhl { + u[j+n] += c + } qhat-- } @@ -827,6 +834,10 @@ func (z nat) divRecursive(u, v nat) { putNat(tmp) } +// divRecursiveStep computes the division of u by v. +// - z must be large enough to hold the quotient +// - the quotient will overwrite z +// - the remainder will overwrite u func (z nat) divRecursiveStep(u, v nat, depth int, tmp *nat, temps []*nat) { u = u.norm() v = v.norm() @@ -1465,19 +1476,26 @@ func (z nat) expNNMontgomery(x, y, m nat) nat { } // bytes writes the value of z into buf using big-endian encoding. -// len(buf) must be >= len(z)*_S. The value of z is encoded in the -// slice buf[i:]. The number i of unused bytes at the beginning of -// buf is returned as result. +// The value of z is encoded in the slice buf[i:]. If the value of z +// cannot be represented in buf, bytes panics. The number i of unused +// bytes at the beginning of buf is returned as result. func (z nat) bytes(buf []byte) (i int) { i = len(buf) for _, d := range z { for j := 0; j < _S; j++ { i-- - buf[i] = byte(d) + if i >= 0 { + buf[i] = byte(d) + } else if byte(d) != 0 { + panic("math/big: buffer too small to fit value") + } d >>= 8 } } + if i < 0 { + i = 0 + } for i < len(buf) && buf[i] == 0 { i++ } diff --git a/libgo/go/math/big/nat_test.go b/libgo/go/math/big/nat_test.go index 32f29e3..89e913f 100644 --- a/libgo/go/math/big/nat_test.go +++ b/libgo/go/math/big/nat_test.go @@ -786,3 +786,21 @@ func TestNatDiv(t *testing.T) { } } } + +// TestIssue37499 triggers the edge case of divBasic where +// the inaccurate estimate of the first word's quotient +// happens at the very beginning of the loop. +func TestIssue37499(t *testing.T) { + // Choose u and v such that v is slightly larger than u >> N. + // This tricks divBasic into choosing 1 as the first word + // of the quotient. This works in both 32-bit and 64-bit settings. + u := natFromString("0x2b6c385a05be027f5c22005b63c42a1165b79ff510e1706b39f8489c1d28e57bb5ba4ef9fd9387a3e344402c0a453381") + v := natFromString("0x2b6c385a05be027f5c22005b63c42a1165b79ff510e1706c") + + q := nat(nil).make(8) + q.divBasic(u, v) + q = q.norm() + if s := string(q.utoa(16)); s != "fffffffffffffffffffffffffffffffffffffffffffffffb" { + t.Fatalf("incorrect quotient: %s", s) + } +} diff --git a/libgo/go/math/big/sqrt.go b/libgo/go/math/big/sqrt.go index 53403aa..0d50164 100644 --- a/libgo/go/math/big/sqrt.go +++ b/libgo/go/math/big/sqrt.go @@ -4,17 +4,29 @@ package big -import "math" - -var ( - three = NewFloat(3.0) +import ( + "math" + "sync" ) +var threeOnce struct { + sync.Once + v *Float +} + +func three() *Float { + threeOnce.Do(func() { + threeOnce.v = NewFloat(3.0) + }) + return threeOnce.v +} + // Sqrt sets z to the rounded square root of x, and returns it. // // If z's precision is 0, it is changed to x's precision before the // operation. Rounding is performed according to z's precision and -// rounding mode. +// rounding mode, but z's accuracy is not computed. Specifically, the +// result of z.Acc() is undefined. // // The function panics if z < 0. The value of z is undefined in that // case. @@ -61,61 +73,14 @@ func (z *Float) Sqrt(x *Float) *Float { } // 0.25 <= z < 2.0 - // Solving x² - z = 0 directly requires a Quo call, but it's - // faster for small precisions. - // - // Solving 1/x² - z = 0 avoids the Quo call and is much faster for - // high precisions. - // - // 128bit precision is an empirically chosen threshold. - if z.prec <= 128 { - z.sqrtDirect(z) - } else { - z.sqrtInverse(z) - } + // Solving 1/x² - z = 0 avoids Quo calls and is faster, especially + // for high precisions. + z.sqrtInverse(z) // re-attach halved exponent return z.SetMantExp(z, b/2) } -// Compute √x (up to prec 128) by solving -// t² - x = 0 -// for t, starting with a 53 bits precision guess from math.Sqrt and -// then using at most two iterations of Newton's method. -func (z *Float) sqrtDirect(x *Float) { - // let - // f(t) = t² - x - // then - // g(t) = f(t)/f'(t) = ½(t² - x)/t - // and the next guess is given by - // t2 = t - g(t) = ½(t² + x)/t - u := new(Float) - ng := func(t *Float) *Float { - u.prec = t.prec - u.Mul(t, t) // u = t² - u.Add(u, x) // = t² + x - u.exp-- // = ½(t² + x) - return t.Quo(u, t) // = ½(t² + x)/t - } - - xf, _ := x.Float64() - sq := NewFloat(math.Sqrt(xf)) - - switch { - case z.prec > 128: - panic("sqrtDirect: only for z.prec <= 128") - case z.prec > 64: - sq.prec *= 2 - sq = ng(sq) - fallthrough - default: - sq.prec *= 2 - sq = ng(sq) - } - - z.Set(sq) -} - // Compute √x (to z.prec precision) by solving // 1/t² - x = 0 // for t (using Newton's method), and then inverting. @@ -128,6 +93,7 @@ func (z *Float) sqrtInverse(x *Float) { // t2 = t - g(t) = ½t(3 - xt²) u := newFloat(z.prec) v := newFloat(z.prec) + three := three() ng := func(t *Float) *Float { u.prec = t.prec v.prec = t.prec @@ -137,7 +103,6 @@ func (z *Float) sqrtInverse(x *Float) { 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_test.go b/libgo/go/math/bits/bits_test.go index c0f4309..23b4539 100644 --- a/libgo/go/math/bits/bits_test.go +++ b/libgo/go/math/bits/bits_test.go @@ -806,6 +806,130 @@ func TestAddSubUint64(t *testing.T) { } } +func TestAdd64OverflowPanic(t *testing.T) { + // Test that 64-bit overflow panics fire correctly. + // These are designed to improve coverage of compiler intrinsics. + tests := []func(uint64, uint64) uint64{ + func(a, b uint64) uint64 { + x, c := Add64(a, b, 0) + if c > 0 { + panic("overflow") + } + return x + }, + func(a, b uint64) uint64 { + x, c := Add64(a, b, 0) + if c != 0 { + panic("overflow") + } + return x + }, + func(a, b uint64) uint64 { + x, c := Add64(a, b, 0) + if c == 1 { + panic("overflow") + } + return x + }, + func(a, b uint64) uint64 { + x, c := Add64(a, b, 0) + if c != 1 { + return x + } + panic("overflow") + }, + func(a, b uint64) uint64 { + x, c := Add64(a, b, 0) + if c == 0 { + return x + } + panic("overflow") + }, + } + for _, test := range tests { + shouldPanic := func(f func()) { + defer func() { + if err := recover(); err == nil { + t.Fatalf("expected panic") + } + }() + f() + } + + // overflow + shouldPanic(func() { test(_M64, 1) }) + shouldPanic(func() { test(1, _M64) }) + shouldPanic(func() { test(_M64, _M64) }) + + // no overflow + test(_M64, 0) + test(0, 0) + test(1, 1) + } +} + +func TestSub64OverflowPanic(t *testing.T) { + // Test that 64-bit overflow panics fire correctly. + // These are designed to improve coverage of compiler intrinsics. + tests := []func(uint64, uint64) uint64{ + func(a, b uint64) uint64 { + x, c := Sub64(a, b, 0) + if c > 0 { + panic("overflow") + } + return x + }, + func(a, b uint64) uint64 { + x, c := Sub64(a, b, 0) + if c != 0 { + panic("overflow") + } + return x + }, + func(a, b uint64) uint64 { + x, c := Sub64(a, b, 0) + if c == 1 { + panic("overflow") + } + return x + }, + func(a, b uint64) uint64 { + x, c := Sub64(a, b, 0) + if c != 1 { + return x + } + panic("overflow") + }, + func(a, b uint64) uint64 { + x, c := Sub64(a, b, 0) + if c == 0 { + return x + } + panic("overflow") + }, + } + for _, test := range tests { + shouldPanic := func(f func()) { + defer func() { + if err := recover(); err == nil { + t.Fatalf("expected panic") + } + }() + f() + } + + // overflow + shouldPanic(func() { test(0, 1) }) + shouldPanic(func() { test(1, _M64) }) + shouldPanic(func() { test(_M64-1, _M64) }) + + // no overflow + test(_M64, 0) + test(0, 0) + test(1, 1) + } +} + func TestMulDiv(t *testing.T) { testMul := func(msg string, f func(x, y uint) (hi, lo uint), x, y, hi, lo uint) { hi1, lo1 := f(x, y) diff --git a/libgo/go/math/cmplx/abs.go b/libgo/go/math/cmplx/abs.go index f3cd107..2f89d1b 100644 --- a/libgo/go/math/cmplx/abs.go +++ b/libgo/go/math/cmplx/abs.go @@ -3,7 +3,8 @@ // license that can be found in the LICENSE file. // Package cmplx provides basic constants and mathematical functions for -// complex numbers. +// complex numbers. Special case handling conforms to the C99 standard +// Annex G IEC 60559-compatible complex arithmetic. package cmplx import "math" diff --git a/libgo/go/math/cmplx/asin.go b/libgo/go/math/cmplx/asin.go index 062f324..30d019e 100644 --- a/libgo/go/math/cmplx/asin.go +++ b/libgo/go/math/cmplx/asin.go @@ -49,8 +49,31 @@ import "math" // Asin returns the inverse sine of x. func Asin(x complex128) complex128 { - if imag(x) == 0 && math.Abs(real(x)) <= 1 { - return complex(math.Asin(real(x)), imag(x)) + switch re, im := real(x), imag(x); { + case im == 0 && math.Abs(re) <= 1: + return complex(math.Asin(re), im) + case re == 0 && math.Abs(im) <= 1: + return complex(re, math.Asinh(im)) + case math.IsNaN(im): + switch { + case re == 0: + return complex(re, math.NaN()) + case math.IsInf(re, 0): + return complex(math.NaN(), re) + default: + return NaN() + } + case math.IsInf(im, 0): + switch { + case math.IsNaN(re): + return x + case math.IsInf(re, 0): + return complex(math.Copysign(math.Pi/4, re), im) + default: + return complex(math.Copysign(0, re), im) + } + case math.IsInf(re, 0): + return complex(math.Copysign(math.Pi/2, re), math.Copysign(re, im)) } ct := complex(-imag(x), real(x)) // i * x xx := x * x @@ -62,8 +85,31 @@ func Asin(x complex128) complex128 { // Asinh returns the inverse hyperbolic sine of x. func Asinh(x complex128) complex128 { - if imag(x) == 0 && math.Abs(real(x)) <= 1 { - return complex(math.Asinh(real(x)), imag(x)) + switch re, im := real(x), imag(x); { + case im == 0 && math.Abs(re) <= 1: + return complex(math.Asinh(re), im) + case re == 0 && math.Abs(im) <= 1: + return complex(re, math.Asin(im)) + case math.IsInf(re, 0): + switch { + case math.IsInf(im, 0): + return complex(re, math.Copysign(math.Pi/4, im)) + case math.IsNaN(im): + return x + default: + return complex(re, math.Copysign(0.0, im)) + } + case math.IsNaN(re): + switch { + case im == 0: + return x + case math.IsInf(im, 0): + return complex(im, re) + default: + return NaN() + } + case math.IsInf(im, 0): + return complex(math.Copysign(im, re), math.Copysign(math.Pi/2, im)) } xx := x * x x1 := complex(1+real(xx), imag(xx)) // 1 + x*x @@ -91,6 +137,9 @@ func Acos(x complex128) complex128 { // Acosh returns the inverse hyperbolic cosine of x. func Acosh(x complex128) complex128 { + if x == 0 { + return complex(0, math.Copysign(math.Pi/2, imag(x))) + } w := Acos(x) if imag(w) <= 0 { return complex(-imag(w), real(w)) // i * w @@ -133,6 +182,19 @@ func Acosh(x complex128) complex128 { // Atan returns the inverse tangent of x. func Atan(x complex128) complex128 { + switch re, im := real(x), imag(x); { + case im == 0: + return complex(math.Atan(re), im) + case re == 0 && math.Abs(im) <= 1: + return complex(re, math.Atanh(im)) + case math.IsInf(im, 0) || math.IsInf(re, 0): + if math.IsNaN(re) { + return complex(math.NaN(), math.Copysign(0, im)) + } + return complex(math.Copysign(math.Pi/2, re), math.Copysign(0, im)) + case math.IsNaN(re) || math.IsNaN(im): + return NaN() + } x2 := real(x) * real(x) a := 1 - x2 - imag(x)*imag(x) if a == 0 { diff --git a/libgo/go/math/cmplx/cmath_test.go b/libgo/go/math/cmplx/cmath_test.go index 57ba76a..3011e83 100644 --- a/libgo/go/math/cmplx/cmath_test.go +++ b/libgo/go/math/cmplx/cmath_test.go @@ -291,48 +291,219 @@ var tanh = []complex128{ (-1.0000000491604982429364892e+00 - 2.901873195374433112227349e-08i), } -// special cases +// huge values along the real axis for testing reducePi in Tan +var hugeIn = []complex128{ + 1 << 28, + 1 << 29, + 1 << 30, + 1 << 35, + -1 << 120, + 1 << 240, + 1 << 300, + -1 << 480, + 1234567891234567 << 180, + -1234567891234567 << 300, +} + +// Results for tanHuge[i] calculated with https://github.com/robpike/ivy +// using 4096 bits of working precision. +var tanHuge = []complex128{ + 5.95641897939639421, + -0.34551069233430392, + -0.78469661331920043, + 0.84276385870875983, + 0.40806638884180424, + -0.37603456702698076, + 4.60901287677810962, + 3.39135965054779932, + -6.76813854009065030, + -0.76417695016604922, +} + +// special cases conform to C99 standard appendix G.6 Complex arithmetic +var inf, nan = math.Inf(1), math.NaN() + var vcAbsSC = []complex128{ NaN(), } var absSC = []float64{ math.NaN(), } -var vcAcosSC = []complex128{ - NaN(), -} -var acosSC = []complex128{ - NaN(), -} -var vcAcoshSC = []complex128{ - NaN(), -} -var acoshSC = []complex128{ - NaN(), -} -var vcAsinSC = []complex128{ - NaN(), -} -var asinSC = []complex128{ - NaN(), -} -var vcAsinhSC = []complex128{ - NaN(), -} -var asinhSC = []complex128{ - NaN(), -} -var vcAtanSC = []complex128{ - NaN(), -} -var atanSC = []complex128{ - NaN(), -} -var vcAtanhSC = []complex128{ - NaN(), -} -var atanhSC = []complex128{ - NaN(), +var acosSC = []struct { + in, + want complex128 +}{ + // G.6.1.1 + {complex(zero, zero), + complex(math.Pi/2, -zero)}, + {complex(-zero, zero), + complex(math.Pi/2, -zero)}, + {complex(zero, nan), + complex(math.Pi/2, nan)}, + {complex(-zero, nan), + complex(math.Pi/2, nan)}, + {complex(1.0, inf), + complex(math.Pi/2, -inf)}, + {complex(1.0, nan), + NaN()}, + {complex(-inf, 1.0), + complex(math.Pi, -inf)}, + {complex(inf, 1.0), + complex(0.0, -inf)}, + {complex(-inf, inf), + complex(3*math.Pi/4, -inf)}, + {complex(inf, inf), + complex(math.Pi/4, -inf)}, + {complex(inf, nan), + complex(nan, -inf)}, // imaginary sign unspecified + {complex(-inf, nan), + complex(nan, inf)}, // imaginary sign unspecified + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + complex(nan, -inf)}, + {NaN(), + NaN()}, +} +var acoshSC = []struct { + in, + want complex128 +}{ + // G.6.2.1 + {complex(zero, zero), + complex(zero, math.Pi/2)}, + {complex(-zero, zero), + complex(zero, math.Pi/2)}, + {complex(1.0, inf), + complex(inf, math.Pi/2)}, + {complex(1.0, nan), + NaN()}, + {complex(-inf, 1.0), + complex(inf, math.Pi)}, + {complex(inf, 1.0), + complex(inf, zero)}, + {complex(-inf, inf), + complex(inf, 3*math.Pi/4)}, + {complex(inf, inf), + complex(inf, math.Pi/4)}, + {complex(inf, nan), + complex(inf, nan)}, + {complex(-inf, nan), + complex(inf, nan)}, + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + complex(inf, nan)}, + {NaN(), + NaN()}, +} +var asinSC = []struct { + in, + want complex128 +}{ + // Derived from Asin(z) = -i * Asinh(i * z), G.6 #7 + {complex(zero, zero), + complex(zero, zero)}, + {complex(1.0, inf), + complex(0, inf)}, + {complex(1.0, nan), + NaN()}, + {complex(inf, 1), + complex(math.Pi/2, inf)}, + {complex(inf, inf), + complex(math.Pi/4, inf)}, + {complex(inf, nan), + complex(nan, inf)}, // imaginary sign unspecified + {complex(nan, zero), + NaN()}, + {complex(nan, 1), + NaN()}, + {complex(nan, inf), + complex(nan, inf)}, + {NaN(), + NaN()}, +} +var asinhSC = []struct { + in, + want complex128 +}{ + // G.6.2.2 + {complex(zero, zero), + complex(zero, zero)}, + {complex(1.0, inf), + complex(inf, math.Pi/2)}, + {complex(1.0, nan), + NaN()}, + {complex(inf, 1.0), + complex(inf, zero)}, + {complex(inf, inf), + complex(inf, math.Pi/4)}, + {complex(inf, nan), + complex(inf, nan)}, + {complex(nan, zero), + complex(nan, zero)}, + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + complex(inf, nan)}, // sign of real part unspecified + {NaN(), + NaN()}, +} +var atanSC = []struct { + in, + want complex128 +}{ + // Derived from Atan(z) = -i * Atanh(i * z), G.6 #7 + {complex(0, zero), + complex(0, zero)}, + {complex(0, nan), + NaN()}, + {complex(1.0, zero), + complex(math.Pi/4, zero)}, + {complex(1.0, inf), + complex(math.Pi/2, zero)}, + {complex(1.0, nan), + NaN()}, + {complex(inf, 1), + complex(math.Pi/2, zero)}, + {complex(inf, inf), + complex(math.Pi/2, zero)}, + {complex(inf, nan), + complex(math.Pi/2, zero)}, + {complex(nan, 1), + NaN()}, + {complex(nan, inf), + complex(nan, zero)}, + {NaN(), + NaN()}, +} +var atanhSC = []struct { + in, + want complex128 +}{ + // G.6.2.3 + {complex(zero, zero), + complex(zero, zero)}, + {complex(zero, nan), + complex(zero, nan)}, + {complex(1.0, zero), + complex(inf, zero)}, + {complex(1.0, inf), + complex(0, math.Pi/2)}, + {complex(1.0, nan), + NaN()}, + {complex(inf, 1.0), + complex(zero, math.Pi/2)}, + {complex(inf, inf), + complex(zero, math.Pi/2)}, + {complex(inf, nan), + complex(0, nan)}, + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + complex(zero, math.Pi/2)}, // sign of real part not specified. + {NaN(), + NaN()}, } var vcConjSC = []complex128{ NaN(), @@ -340,23 +511,105 @@ var vcConjSC = []complex128{ var conjSC = []complex128{ NaN(), } -var vcCosSC = []complex128{ - NaN(), -} -var cosSC = []complex128{ - NaN(), -} -var vcCoshSC = []complex128{ - NaN(), -} -var coshSC = []complex128{ - NaN(), -} -var vcExpSC = []complex128{ - NaN(), -} -var expSC = []complex128{ - NaN(), +var cosSC = []struct { + in, + want complex128 +}{ + // Derived from Cos(z) = Cosh(i * z), G.6 #7 + {complex(zero, zero), + complex(1.0, -zero)}, + {complex(zero, inf), + complex(inf, -zero)}, + {complex(zero, nan), + complex(nan, zero)}, // imaginary sign unspecified + {complex(1.0, inf), + complex(inf, -inf)}, + {complex(1.0, nan), + NaN()}, + {complex(inf, zero), + complex(nan, -zero)}, + {complex(inf, 1.0), + NaN()}, + {complex(inf, inf), + complex(inf, nan)}, // real sign unspecified + {complex(inf, nan), + NaN()}, + {complex(nan, zero), + complex(nan, -zero)}, // imaginary sign unspecified + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + complex(inf, nan)}, + {NaN(), + NaN()}, +} +var coshSC = []struct { + in, + want complex128 +}{ + // G.6.2.4 + {complex(zero, zero), + complex(1.0, zero)}, + {complex(zero, inf), + complex(nan, zero)}, // imaginary sign unspecified + {complex(zero, nan), + complex(nan, zero)}, // imaginary sign unspecified + {complex(1.0, inf), + NaN()}, + {complex(1.0, nan), + NaN()}, + {complex(inf, zero), + complex(inf, zero)}, + {complex(inf, 1.0), + complex(inf*math.Cos(1.0), inf*math.Sin(1.0))}, // +inf cis(y) + {complex(inf, inf), + complex(inf, nan)}, // real sign unspecified + {complex(inf, nan), + complex(inf, nan)}, + {complex(nan, zero), + complex(nan, zero)}, // imaginary sign unspecified + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + NaN()}, + {NaN(), + NaN()}, +} +var expSC = []struct { + in, + want complex128 +}{ + // G.6.3.1 + {complex(zero, zero), + complex(1.0, zero)}, + {complex(-zero, zero), + complex(1.0, zero)}, + {complex(1.0, inf), + NaN()}, + {complex(1.0, nan), + NaN()}, + {complex(inf, zero), + complex(inf, zero)}, + {complex(-inf, 1.0), + complex(math.Copysign(0.0, math.Cos(1.0)), math.Copysign(0.0, math.Sin(1.0)))}, // +0 cis(y) + {complex(inf, 1.0), + complex(inf*math.Cos(1.0), inf*math.Sin(1.0))}, // +inf cis(y) + {complex(-inf, inf), + complex(zero, zero)}, // real and imaginary sign unspecified + {complex(inf, inf), + complex(inf, nan)}, // real sign unspecified + {complex(-inf, nan), + complex(zero, zero)}, // real and imaginary sign unspecified + {complex(inf, nan), + complex(inf, nan)}, // real sign unspecified + {complex(nan, zero), + complex(nan, zero)}, + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + NaN()}, + {NaN(), + NaN()}, } var vcIsNaNSC = []complex128{ complex(math.Inf(-1), math.Inf(-1)), @@ -380,17 +633,70 @@ var isNaNSC = []bool{ false, true, } -var vcLogSC = []complex128{ - NaN(), -} -var logSC = []complex128{ - NaN(), -} -var vcLog10SC = []complex128{ - NaN(), -} -var log10SC = []complex128{ - NaN(), + +var logSC = []struct { + in, + want complex128 +}{ + // G.6.3.2 + {complex(zero, zero), + complex(-inf, zero)}, + {complex(-zero, zero), + complex(-inf, math.Pi)}, + {complex(1.0, inf), + complex(inf, math.Pi/2)}, + {complex(1.0, nan), + NaN()}, + {complex(-inf, 1.0), + complex(inf, math.Pi)}, + {complex(inf, 1.0), + complex(inf, 0.0)}, + {complex(-inf, inf), + complex(inf, 3*math.Pi/4)}, + {complex(inf, inf), + complex(inf, math.Pi/4)}, + {complex(-inf, nan), + complex(inf, nan)}, + {complex(inf, nan), + complex(inf, nan)}, + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + complex(inf, nan)}, + {NaN(), + NaN()}, +} +var log10SC = []struct { + in, + want complex128 +}{ + // derived from Log special cases via Log10(x) = math.Log10E*Log(x) + {complex(zero, zero), + complex(-inf, zero)}, + {complex(-zero, zero), + complex(-inf, float64(math.Log10E)*float64(math.Pi))}, + {complex(1.0, inf), + complex(inf, float64(math.Log10E)*float64(math.Pi/2))}, + {complex(1.0, nan), + NaN()}, + {complex(-inf, 1.0), + complex(inf, float64(math.Log10E)*float64(math.Pi))}, + {complex(inf, 1.0), + complex(inf, 0.0)}, + {complex(-inf, inf), + complex(inf, float64(math.Log10E)*float64(3*math.Pi/4))}, + {complex(inf, inf), + complex(inf, float64(math.Log10E)*float64(math.Pi/4))}, + {complex(-inf, nan), + complex(inf, nan)}, + {complex(inf, nan), + complex(inf, nan)}, + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + complex(inf, nan)}, + {NaN(), + NaN()}, } var vcPolarSC = []complex128{ NaN(), @@ -406,35 +712,153 @@ var powSC = []complex128{ NaN(), NaN(), } -var vcSinSC = []complex128{ - NaN(), -} -var sinSC = []complex128{ - NaN(), -} -var vcSinhSC = []complex128{ - NaN(), -} -var sinhSC = []complex128{ - NaN(), -} -var vcSqrtSC = []complex128{ - NaN(), -} -var sqrtSC = []complex128{ - NaN(), -} -var vcTanSC = []complex128{ - NaN(), -} -var tanSC = []complex128{ - NaN(), +var sinSC = []struct { + in, + want complex128 +}{ + // Derived from Sin(z) = -i * Sinh(i * z), G.6 #7 + {complex(zero, zero), + complex(zero, zero)}, + {complex(zero, inf), + complex(zero, inf)}, + {complex(zero, nan), + complex(zero, nan)}, + {complex(1.0, inf), + complex(inf, inf)}, + {complex(1.0, nan), + NaN()}, + {complex(inf, zero), + complex(nan, zero)}, + {complex(inf, 1.0), + NaN()}, + {complex(inf, inf), + complex(nan, inf)}, + {complex(inf, nan), + NaN()}, + {complex(nan, zero), + complex(nan, zero)}, + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + complex(nan, inf)}, + {NaN(), + NaN()}, } -var vcTanhSC = []complex128{ - NaN(), + +var sinhSC = []struct { + in, + want complex128 +}{ + // G.6.2.5 + {complex(zero, zero), + complex(zero, zero)}, + {complex(zero, inf), + complex(zero, nan)}, // real sign unspecified + {complex(zero, nan), + complex(zero, nan)}, // real sign unspecified + {complex(1.0, inf), + NaN()}, + {complex(1.0, nan), + NaN()}, + {complex(inf, zero), + complex(inf, zero)}, + {complex(inf, 1.0), + complex(inf*math.Cos(1.0), inf*math.Sin(1.0))}, // +inf cis(y) + {complex(inf, inf), + complex(inf, nan)}, // real sign unspecified + {complex(inf, nan), + complex(inf, nan)}, // real sign unspecified + {complex(nan, zero), + complex(nan, zero)}, + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + NaN()}, + {NaN(), + NaN()}, } -var tanhSC = []complex128{ - NaN(), + +var sqrtSC = []struct { + in, + want complex128 +}{ + // G.6.4.2 + {complex(zero, zero), + complex(zero, zero)}, + {complex(-zero, zero), + complex(zero, zero)}, + {complex(1.0, inf), + complex(inf, inf)}, + {complex(nan, inf), + complex(inf, inf)}, + {complex(1.0, nan), + NaN()}, + {complex(-inf, 1.0), + complex(zero, inf)}, + {complex(inf, 1.0), + complex(inf, zero)}, + {complex(-inf, nan), + complex(nan, inf)}, // imaginary sign unspecified + {complex(inf, nan), + complex(inf, nan)}, + {complex(nan, 1.0), + NaN()}, + {NaN(), + NaN()}, +} +var tanSC = []struct { + in, + want complex128 +}{ + // Derived from Tan(z) = -i * Tanh(i * z), G.6 #7 + {complex(zero, zero), + complex(zero, zero)}, + {complex(zero, nan), + complex(zero, nan)}, + {complex(1.0, inf), + complex(zero, 1.0)}, + {complex(1.0, nan), + NaN()}, + {complex(inf, 1.0), + NaN()}, + {complex(inf, inf), + complex(zero, 1.0)}, + {complex(inf, nan), + NaN()}, + {complex(nan, zero), + NaN()}, + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + complex(zero, 1.0)}, + {NaN(), + NaN()}, +} +var tanhSC = []struct { + in, + want complex128 +}{ + // G.6.2.6 + {complex(zero, zero), + complex(zero, zero)}, + {complex(1.0, inf), + NaN()}, + {complex(1.0, nan), + NaN()}, + {complex(inf, 1.0), + complex(1.0, math.Copysign(0.0, math.Sin(2*1.0)))}, // 1 + i 0 sin(2y) + {complex(inf, inf), + complex(1.0, zero)}, // imaginary sign unspecified + {complex(inf, nan), + complex(1.0, zero)}, // imaginary sign unspecified + {complex(nan, zero), + complex(nan, zero)}, + {complex(nan, 1.0), + NaN()}, + {complex(nan, inf), + NaN()}, + {NaN(), + NaN()}, } // branch cut continuity checks @@ -496,13 +920,25 @@ func cTolerance(a, b complex128, e float64) bool { func cSoclose(a, b complex128, e float64) bool { return cTolerance(a, b, e) } func cVeryclose(a, b complex128) bool { return cTolerance(a, b, 4e-16) } func cAlike(a, b complex128) bool { - switch { - case IsNaN(a) && IsNaN(b): - return true - case a == b: - return math.Signbit(real(a)) == math.Signbit(real(b)) && math.Signbit(imag(a)) == math.Signbit(imag(b)) - } - return false + var realAlike, imagAlike bool + if isExact(real(b)) { + realAlike = alike(real(a), real(b)) + } else { + // Allow non-exact special cases to have errors in ULP. + realAlike = veryclose(real(a), real(b)) + } + if isExact(imag(b)) { + imagAlike = alike(imag(a), imag(b)) + } else { + // Allow non-exact special cases to have errors in ULP. + imagAlike = veryclose(imag(a), imag(b)) + } + return realAlike && imagAlike +} +func isExact(x float64) bool { + // Special cases that should match exactly. Other cases are multiples + // of Pi that may not be last bit identical on all platforms. + return math.IsNaN(x) || math.IsInf(x, 0) || x == 0 || x == 1 || x == -1 } func TestAbs(t *testing.T) { @@ -523,9 +959,17 @@ func TestAcos(t *testing.T) { t.Errorf("Acos(%g) = %g, want %g", vc[i], f, acos[i]) } } - for i := 0; i < len(vcAcosSC); i++ { - if f := Acos(vcAcosSC[i]); !cAlike(acosSC[i], f) { - t.Errorf("Acos(%g) = %g, want %g", vcAcosSC[i], f, acosSC[i]) + for _, v := range acosSC { + if f := Acos(v.in); !cAlike(v.want, f) { + t.Errorf("Acos(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Acos(Conj(z)) == Conj(Acos(z)) + if f := Acos(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Acos(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) } } for _, pt := range branchPoints { @@ -540,10 +984,19 @@ func TestAcosh(t *testing.T) { t.Errorf("Acosh(%g) = %g, want %g", vc[i], f, acosh[i]) } } - for i := 0; i < len(vcAcoshSC); i++ { - if f := Acosh(vcAcoshSC[i]); !cAlike(acoshSC[i], f) { - t.Errorf("Acosh(%g) = %g, want %g", vcAcoshSC[i], f, acoshSC[i]) + for _, v := range acoshSC { + if f := Acosh(v.in); !cAlike(v.want, f) { + t.Errorf("Acosh(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue } + // Acosh(Conj(z)) == Conj(Acosh(z)) + if f := Acosh(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Acosh(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + } for _, pt := range branchPoints { if f0, f1 := Acosh(pt[0]), Acosh(pt[1]); !cVeryclose(f0, f1) { @@ -557,9 +1010,25 @@ func TestAsin(t *testing.T) { t.Errorf("Asin(%g) = %g, want %g", vc[i], f, asin[i]) } } - for i := 0; i < len(vcAsinSC); i++ { - if f := Asin(vcAsinSC[i]); !cAlike(asinSC[i], f) { - t.Errorf("Asin(%g) = %g, want %g", vcAsinSC[i], f, asinSC[i]) + for _, v := range asinSC { + if f := Asin(v.in); !cAlike(v.want, f) { + t.Errorf("Asin(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Asin(Conj(z)) == Asin(Sinh(z)) + if f := Asin(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Asin(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + if math.IsNaN(real(v.in)) || math.IsNaN(real(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Asin(-z) == -Asin(z) + if f := Asin(-v.in); !cAlike(-v.want, f) && !cAlike(v.in, -v.in) { + t.Errorf("Asin(%g) = %g, want %g", -v.in, f, -v.want) } } for _, pt := range branchPoints { @@ -574,9 +1043,25 @@ func TestAsinh(t *testing.T) { t.Errorf("Asinh(%g) = %g, want %g", vc[i], f, asinh[i]) } } - for i := 0; i < len(vcAsinhSC); i++ { - if f := Asinh(vcAsinhSC[i]); !cAlike(asinhSC[i], f) { - t.Errorf("Asinh(%g) = %g, want %g", vcAsinhSC[i], f, asinhSC[i]) + for _, v := range asinhSC { + if f := Asinh(v.in); !cAlike(v.want, f) { + t.Errorf("Asinh(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Asinh(Conj(z)) == Asinh(Sinh(z)) + if f := Asinh(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Asinh(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + if math.IsNaN(real(v.in)) || math.IsNaN(real(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Asinh(-z) == -Asinh(z) + if f := Asinh(-v.in); !cAlike(-v.want, f) && !cAlike(v.in, -v.in) { + t.Errorf("Asinh(%g) = %g, want %g", -v.in, f, -v.want) } } for _, pt := range branchPoints { @@ -591,9 +1076,25 @@ func TestAtan(t *testing.T) { t.Errorf("Atan(%g) = %g, want %g", vc[i], f, atan[i]) } } - for i := 0; i < len(vcAtanSC); i++ { - if f := Atan(vcAtanSC[i]); !cAlike(atanSC[i], f) { - t.Errorf("Atan(%g) = %g, want %g", vcAtanSC[i], f, atanSC[i]) + for _, v := range atanSC { + if f := Atan(v.in); !cAlike(v.want, f) { + t.Errorf("Atan(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Atan(Conj(z)) == Conj(Atan(z)) + if f := Atan(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Atan(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + if math.IsNaN(real(v.in)) || math.IsNaN(real(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Atan(-z) == -Atan(z) + if f := Atan(-v.in); !cAlike(-v.want, f) && !cAlike(v.in, -v.in) { + t.Errorf("Atan(%g) = %g, want %g", -v.in, f, -v.want) } } for _, pt := range branchPoints { @@ -608,9 +1109,25 @@ func TestAtanh(t *testing.T) { t.Errorf("Atanh(%g) = %g, want %g", vc[i], f, atanh[i]) } } - for i := 0; i < len(vcAtanhSC); i++ { - if f := Atanh(vcAtanhSC[i]); !cAlike(atanhSC[i], f) { - t.Errorf("Atanh(%g) = %g, want %g", vcAtanhSC[i], f, atanhSC[i]) + for _, v := range atanhSC { + if f := Atanh(v.in); !cAlike(v.want, f) { + t.Errorf("Atanh(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Atanh(Conj(z)) == Conj(Atanh(z)) + if f := Atanh(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Atanh(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + if math.IsNaN(real(v.in)) || math.IsNaN(real(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Atanh(-z) == -Atanh(z) + if f := Atanh(-v.in); !cAlike(-v.want, f) && !cAlike(v.in, -v.in) { + t.Errorf("Atanh(%g) = %g, want %g", -v.in, f, -v.want) } } for _, pt := range branchPoints { @@ -637,9 +1154,25 @@ func TestCos(t *testing.T) { t.Errorf("Cos(%g) = %g, want %g", vc[i], f, cos[i]) } } - for i := 0; i < len(vcCosSC); i++ { - if f := Cos(vcCosSC[i]); !cAlike(cosSC[i], f) { - t.Errorf("Cos(%g) = %g, want %g", vcCosSC[i], f, cosSC[i]) + for _, v := range cosSC { + if f := Cos(v.in); !cAlike(v.want, f) { + t.Errorf("Cos(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Cos(Conj(z)) == Cos(Cosh(z)) + if f := Cos(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Cos(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + if math.IsNaN(real(v.in)) || math.IsNaN(real(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Cos(-z) == Cos(z) + if f := Cos(-v.in); !cAlike(v.want, f) && !cAlike(v.in, -v.in) { + t.Errorf("Cos(%g) = %g, want %g", -v.in, f, v.want) } } } @@ -649,9 +1182,25 @@ func TestCosh(t *testing.T) { t.Errorf("Cosh(%g) = %g, want %g", vc[i], f, cosh[i]) } } - for i := 0; i < len(vcCoshSC); i++ { - if f := Cosh(vcCoshSC[i]); !cAlike(coshSC[i], f) { - t.Errorf("Cosh(%g) = %g, want %g", vcCoshSC[i], f, coshSC[i]) + for _, v := range coshSC { + if f := Cosh(v.in); !cAlike(v.want, f) { + t.Errorf("Cosh(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Cosh(Conj(z)) == Conj(Cosh(z)) + if f := Cosh(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Cosh(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + if math.IsNaN(real(v.in)) || math.IsNaN(real(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Cosh(-z) == Cosh(z) + if f := Cosh(-v.in); !cAlike(v.want, f) && !cAlike(v.in, -v.in) { + t.Errorf("Cosh(%g) = %g, want %g", -v.in, f, v.want) } } } @@ -661,9 +1210,17 @@ func TestExp(t *testing.T) { t.Errorf("Exp(%g) = %g, want %g", vc[i], f, exp[i]) } } - for i := 0; i < len(vcExpSC); i++ { - if f := Exp(vcExpSC[i]); !cAlike(expSC[i], f) { - t.Errorf("Exp(%g) = %g, want %g", vcExpSC[i], f, expSC[i]) + for _, v := range expSC { + if f := Exp(v.in); !cAlike(v.want, f) { + t.Errorf("Exp(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Exp(Conj(z)) == Exp(Cosh(z)) + if f := Exp(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Exp(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) } } } @@ -680,9 +1237,17 @@ func TestLog(t *testing.T) { t.Errorf("Log(%g) = %g, want %g", vc[i], f, log[i]) } } - for i := 0; i < len(vcLogSC); i++ { - if f := Log(vcLogSC[i]); !cAlike(logSC[i], f) { - t.Errorf("Log(%g) = %g, want %g", vcLogSC[i], f, logSC[i]) + for _, v := range logSC { + if f := Log(v.in); !cAlike(v.want, f) { + t.Errorf("Log(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Log(Conj(z)) == Conj(Log(z)) + if f := Log(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Log(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) } } for _, pt := range branchPoints { @@ -697,9 +1262,17 @@ func TestLog10(t *testing.T) { t.Errorf("Log10(%g) = %g, want %g", vc[i], f, log10[i]) } } - for i := 0; i < len(vcLog10SC); i++ { - if f := Log10(vcLog10SC[i]); !cAlike(log10SC[i], f) { - t.Errorf("Log10(%g) = %g, want %g", vcLog10SC[i], f, log10SC[i]) + for _, v := range log10SC { + if f := Log10(v.in); !cAlike(v.want, f) { + t.Errorf("Log10(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Log10(Conj(z)) == Conj(Log10(z)) + if f := Log10(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Log10(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) } } } @@ -764,9 +1337,25 @@ func TestSin(t *testing.T) { t.Errorf("Sin(%g) = %g, want %g", vc[i], f, sin[i]) } } - for i := 0; i < len(vcSinSC); i++ { - if f := Sin(vcSinSC[i]); !cAlike(sinSC[i], f) { - t.Errorf("Sin(%g) = %g, want %g", vcSinSC[i], f, sinSC[i]) + for _, v := range sinSC { + if f := Sin(v.in); !cAlike(v.want, f) { + t.Errorf("Sin(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Sin(Conj(z)) == Conj(Sin(z)) + if f := Sin(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Sinh(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + if math.IsNaN(real(v.in)) || math.IsNaN(real(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Sin(-z) == -Sin(z) + if f := Sin(-v.in); !cAlike(-v.want, f) && !cAlike(v.in, -v.in) { + t.Errorf("Sinh(%g) = %g, want %g", -v.in, f, -v.want) } } } @@ -776,9 +1365,25 @@ func TestSinh(t *testing.T) { t.Errorf("Sinh(%g) = %g, want %g", vc[i], f, sinh[i]) } } - for i := 0; i < len(vcSinhSC); i++ { - if f := Sinh(vcSinhSC[i]); !cAlike(sinhSC[i], f) { - t.Errorf("Sinh(%g) = %g, want %g", vcSinhSC[i], f, sinhSC[i]) + for _, v := range sinhSC { + if f := Sinh(v.in); !cAlike(v.want, f) { + t.Errorf("Sinh(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Sinh(Conj(z)) == Conj(Sinh(z)) + if f := Sinh(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Sinh(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + if math.IsNaN(real(v.in)) || math.IsNaN(real(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Sinh(-z) == -Sinh(z) + if f := Sinh(-v.in); !cAlike(-v.want, f) && !cAlike(v.in, -v.in) { + t.Errorf("Sinh(%g) = %g, want %g", -v.in, f, -v.want) } } } @@ -788,9 +1393,17 @@ func TestSqrt(t *testing.T) { t.Errorf("Sqrt(%g) = %g, want %g", vc[i], f, sqrt[i]) } } - for i := 0; i < len(vcSqrtSC); i++ { - if f := Sqrt(vcSqrtSC[i]); !cAlike(sqrtSC[i], f) { - t.Errorf("Sqrt(%g) = %g, want %g", vcSqrtSC[i], f, sqrtSC[i]) + for _, v := range sqrtSC { + if f := Sqrt(v.in); !cAlike(v.want, f) { + t.Errorf("Sqrt(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Sqrt(Conj(z)) == Conj(Sqrt(z)) + if f := Sqrt(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Sqrt(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) } } for _, pt := range branchPoints { @@ -805,9 +1418,25 @@ func TestTan(t *testing.T) { t.Errorf("Tan(%g) = %g, want %g", vc[i], f, tan[i]) } } - for i := 0; i < len(vcTanSC); i++ { - if f := Tan(vcTanSC[i]); !cAlike(tanSC[i], f) { - t.Errorf("Tan(%g) = %g, want %g", vcTanSC[i], f, tanSC[i]) + for _, v := range tanSC { + if f := Tan(v.in); !cAlike(v.want, f) { + t.Errorf("Tan(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Tan(Conj(z)) == Conj(Tan(z)) + if f := Tan(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Tan(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + if math.IsNaN(real(v.in)) || math.IsNaN(real(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Tan(-z) == -Tan(z) + if f := Tan(-v.in); !cAlike(-v.want, f) && !cAlike(v.in, -v.in) { + t.Errorf("Tan(%g) = %g, want %g", -v.in, f, -v.want) } } } @@ -817,9 +1446,25 @@ func TestTanh(t *testing.T) { t.Errorf("Tanh(%g) = %g, want %g", vc[i], f, tanh[i]) } } - for i := 0; i < len(vcTanhSC); i++ { - if f := Tanh(vcTanhSC[i]); !cAlike(tanhSC[i], f) { - t.Errorf("Tanh(%g) = %g, want %g", vcTanhSC[i], f, tanhSC[i]) + for _, v := range tanhSC { + if f := Tanh(v.in); !cAlike(v.want, f) { + t.Errorf("Tanh(%g) = %g, want %g", v.in, f, v.want) + } + if math.IsNaN(imag(v.in)) || math.IsNaN(imag(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Tanh(Conj(z)) == Conj(Tanh(z)) + if f := Tanh(Conj(v.in)); !cAlike(Conj(v.want), f) && !cAlike(v.in, Conj(v.in)) { + t.Errorf("Tanh(%g) = %g, want %g", Conj(v.in), f, Conj(v.want)) + } + if math.IsNaN(real(v.in)) || math.IsNaN(real(v.want)) { + // Negating NaN is undefined with regard to the sign bit produced. + continue + } + // Tanh(-z) == -Tanh(z) + if f := Tanh(-v.in); !cAlike(-v.want, f) && !cAlike(v.in, -v.in) { + t.Errorf("Tanh(%g) = %g, want %g", -v.in, f, -v.want) } } } diff --git a/libgo/go/math/cmplx/exp.go b/libgo/go/math/cmplx/exp.go index 485ed2c..d5d0a5d 100644 --- a/libgo/go/math/cmplx/exp.go +++ b/libgo/go/math/cmplx/exp.go @@ -49,6 +49,23 @@ import "math" // Exp returns e**x, the base-e exponential of x. func Exp(x complex128) complex128 { + switch re, im := real(x), imag(x); { + case math.IsInf(re, 0): + switch { + case re > 0 && im == 0: + return x + case math.IsInf(im, 0) || math.IsNaN(im): + if re < 0 { + return complex(0, math.Copysign(0, im)) + } else { + return complex(math.Inf(1.0), math.NaN()) + } + } + case math.IsNaN(re): + if im == 0 { + return complex(math.NaN(), im) + } + } r := math.Exp(real(x)) s, c := math.Sincos(imag(x)) return complex(r*c, r*s) diff --git a/libgo/go/math/cmplx/huge_test.go b/libgo/go/math/cmplx/huge_test.go new file mode 100644 index 0000000..f8e60c2 --- /dev/null +++ b/libgo/go/math/cmplx/huge_test.go @@ -0,0 +1,22 @@ +// Copyright 2020 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. + +// Disabled for s390x because it uses assembly routines that are not +// accurate for huge arguments. + +// +build !s390x + +package cmplx + +import ( + "testing" +) + +func TestTanHuge(t *testing.T) { + for i, x := range hugeIn { + if f := Tan(x); !cSoclose(tanHuge[i], f, 3e-15) { + t.Errorf("Tan(%g) = %g, want %g", x, f, tanHuge[i]) + } + } +} diff --git a/libgo/go/math/cmplx/log.go b/libgo/go/math/cmplx/log.go index 881a064..fd39c76 100644 --- a/libgo/go/math/cmplx/log.go +++ b/libgo/go/math/cmplx/log.go @@ -60,5 +60,6 @@ func Log(x complex128) complex128 { // Log10 returns the decimal logarithm of x. func Log10(x complex128) complex128 { - return math.Log10E * Log(x) + z := Log(x) + return complex(math.Log10E*real(z), math.Log10E*imag(z)) } diff --git a/libgo/go/math/cmplx/sin.go b/libgo/go/math/cmplx/sin.go index 2c57536..febac0e 100644 --- a/libgo/go/math/cmplx/sin.go +++ b/libgo/go/math/cmplx/sin.go @@ -51,6 +51,19 @@ import "math" // Sin returns the sine of x. func Sin(x complex128) complex128 { + switch re, im := real(x), imag(x); { + case im == 0 && (math.IsInf(re, 0) || math.IsNaN(re)): + return complex(math.NaN(), im) + case math.IsInf(im, 0): + switch { + case re == 0: + return x + case math.IsInf(re, 0) || math.IsNaN(re): + return complex(math.NaN(), im) + } + case re == 0 && math.IsNaN(im): + return x + } s, c := math.Sincos(real(x)) sh, ch := sinhcosh(imag(x)) return complex(s*ch, c*sh) @@ -71,6 +84,19 @@ func Sin(x complex128) complex128 { // Sinh returns the hyperbolic sine of x. func Sinh(x complex128) complex128 { + switch re, im := real(x), imag(x); { + case re == 0 && (math.IsInf(im, 0) || math.IsNaN(im)): + return complex(re, math.NaN()) + case math.IsInf(re, 0): + switch { + case im == 0: + return complex(re, im) + case math.IsInf(im, 0) || math.IsNaN(im): + return complex(re, math.NaN()) + } + case im == 0 && math.IsNaN(re): + return complex(math.NaN(), im) + } s, c := math.Sincos(imag(x)) sh, ch := sinhcosh(real(x)) return complex(c*sh, s*ch) @@ -96,6 +122,19 @@ func Sinh(x complex128) complex128 { // Cos returns the cosine of x. func Cos(x complex128) complex128 { + switch re, im := real(x), imag(x); { + case im == 0 && (math.IsInf(re, 0) || math.IsNaN(re)): + return complex(math.NaN(), -im*math.Copysign(0, re)) + case math.IsInf(im, 0): + switch { + case re == 0: + return complex(math.Inf(1), -re*math.Copysign(0, im)) + case math.IsInf(re, 0) || math.IsNaN(re): + return complex(math.Inf(1), math.NaN()) + } + case re == 0 && math.IsNaN(im): + return complex(math.NaN(), 0) + } s, c := math.Sincos(real(x)) sh, ch := sinhcosh(imag(x)) return complex(c*ch, -s*sh) @@ -115,6 +154,19 @@ func Cos(x complex128) complex128 { // Cosh returns the hyperbolic cosine of x. func Cosh(x complex128) complex128 { + switch re, im := real(x), imag(x); { + case re == 0 && (math.IsInf(im, 0) || math.IsNaN(im)): + return complex(math.NaN(), re*math.Copysign(0, im)) + case math.IsInf(re, 0): + switch { + case im == 0: + return complex(math.Inf(1), im*math.Copysign(0, re)) + case math.IsInf(im, 0) || math.IsNaN(im): + return complex(math.Inf(1), math.NaN()) + } + case im == 0 && math.IsNaN(re): + return complex(math.NaN(), im) + } s, c := math.Sincos(imag(x)) sh, ch := sinhcosh(real(x)) return complex(c*ch, s*sh) diff --git a/libgo/go/math/cmplx/sqrt.go b/libgo/go/math/cmplx/sqrt.go index 0fbdcde..eddce2f 100644 --- a/libgo/go/math/cmplx/sqrt.go +++ b/libgo/go/math/cmplx/sqrt.go @@ -40,7 +40,7 @@ import "math" // 1/2 // Im w = [ (r - x)/2 ] . // -// Cancelation error in r-x or r+x is avoided by using the +// Cancellation error in r-x or r+x is avoided by using the // identity 2 Re w Im w = y. // // Note that -w is also a square root of z. The root chosen @@ -65,6 +65,8 @@ func Sqrt(x complex128) complex128 { return complex(0, math.Copysign(math.Sqrt(-real(x)), imag(x))) } return complex(math.Sqrt(real(x)), imag(x)) + } else if math.IsInf(imag(x), 0) { + return complex(math.Inf(1.0), imag(x)) } if real(x) == 0 { if imag(x) < 0 { diff --git a/libgo/go/math/cmplx/tan.go b/libgo/go/math/cmplx/tan.go index 0243ea0..67a1133 100644 --- a/libgo/go/math/cmplx/tan.go +++ b/libgo/go/math/cmplx/tan.go @@ -4,7 +4,10 @@ package cmplx -import "math" +import ( + "math" + "math/bits" +) // The original C code, the long comment, and the constants // below are from http://netlib.sandia.gov/cephes/c9x-complex/clog.c. @@ -42,7 +45,7 @@ import "math" // cos 2x + cosh 2y // // On the real axis the denominator is zero at odd multiples -// of PI/2. The denominator is evaluated by its Taylor +// of PI/2. The denominator is evaluated by its Taylor // series near these points. // // ctan(z) = -i ctanh(iz). @@ -57,6 +60,16 @@ import "math" // Tan returns the tangent of x. func Tan(x complex128) complex128 { + switch re, im := real(x), imag(x); { + case math.IsInf(im, 0): + switch { + case math.IsInf(re, 0) || math.IsNaN(re): + return complex(math.Copysign(0, re), math.Copysign(1, im)) + } + return complex(math.Copysign(0, math.Sin(2*re)), math.Copysign(1, im)) + case re == 0 && math.IsNaN(im): + return x + } d := math.Cos(2*real(x)) + math.Cosh(2*imag(x)) if math.Abs(d) < 0.25 { d = tanSeries(x) @@ -81,6 +94,16 @@ func Tan(x complex128) complex128 { // Tanh returns the hyperbolic tangent of x. func Tanh(x complex128) complex128 { + switch re, im := real(x), imag(x); { + case math.IsInf(re, 0): + switch { + case math.IsInf(im, 0) || math.IsNaN(im): + return complex(math.Copysign(1, re), math.Copysign(0, im)) + } + return complex(math.Copysign(1, re), math.Copysign(0, math.Sin(2*im))) + case im == 0 && math.IsNaN(re): + return x + } d := math.Cosh(2*real(x)) + math.Cos(2*imag(x)) if d == 0 { return Inf() @@ -88,22 +111,110 @@ func Tanh(x complex128) complex128 { return complex(math.Sinh(2*real(x))/d, math.Sin(2*imag(x))/d) } -// Program to subtract nearest integer multiple of PI +// reducePi reduces the input argument x to the range (-Pi/2, Pi/2]. +// x must be greater than or equal to 0. For small arguments it +// uses Cody-Waite reduction in 3 float64 parts based on: +// "Elementary Function Evaluation: Algorithms and Implementation" +// Jean-Michel Muller, 1997. +// For very large arguments it uses Payne-Hanek range reduction based on: +// "ARGUMENT REDUCTION FOR HUGE ARGUMENTS: Good to the Last Bit" +// K. C. Ng et al, March 24, 1992. func reducePi(x float64) float64 { + // reduceThreshold is the maximum value of x where the reduction using + // Cody-Waite reduction still gives accurate results. This threshold + // is set by t*PIn being representable as a float64 without error + // where t is given by t = floor(x * (1 / Pi)) and PIn are the leading partial + // terms of Pi. Since the leading terms, PI1 and PI2 below, have 30 and 32 + // trailing zero bits respectively, t should have less than 30 significant bits. + // t < 1<<30 -> floor(x*(1/Pi)+0.5) < 1<<30 -> x < (1<<30-1) * Pi - 0.5 + // So, conservatively we can take x < 1<<30. + const reduceThreshold float64 = 1 << 30 + if math.Abs(x) < reduceThreshold { + // Use Cody-Waite reduction in three parts. + const ( + // PI1, PI2 and PI3 comprise an extended precision value of PI + // such that PI ~= PI1 + PI2 + PI3. The parts are chosen so + // that PI1 and PI2 have an approximately equal number of trailing + // zero bits. This ensures that t*PI1 and t*PI2 are exact for + // large integer values of t. The full precision PI3 ensures the + // approximation of PI is accurate to 102 bits to handle cancellation + // during subtraction. + PI1 = 3.141592502593994 // 0x400921fb40000000 + PI2 = 1.5099578831723193e-07 // 0x3e84442d00000000 + PI3 = 1.0780605716316238e-14 // 0x3d08469898cc5170 + ) + t := x / math.Pi + t += 0.5 + t = float64(int64(t)) // int64(t) = the multiple + return ((x - t*PI1) - t*PI2) - t*PI3 + } + // Must apply Payne-Hanek range reduction const ( - // extended precision value of PI: - DP1 = 3.14159265160560607910e0 // ?? 0x400921fb54000000 - DP2 = 1.98418714791870343106e-9 // ?? 0x3e210b4610000000 - DP3 = 1.14423774522196636802e-17 // ?? 0x3c6a62633145c06e + mask = 0x7FF + shift = 64 - 11 - 1 + bias = 1023 + fracMask = 1<<shift - 1 ) - t := x / math.Pi - if t >= 0 { - t += 0.5 - } else { - t -= 0.5 + // Extract out the integer and exponent such that, + // x = ix * 2 ** exp. + ix := math.Float64bits(x) + exp := int(ix>>shift&mask) - bias - shift + ix &= fracMask + ix |= 1 << shift + + // mPi is the binary digits of 1/Pi as a uint64 array, + // that is, 1/Pi = Sum mPi[i]*2^(-64*i). + // 19 64-bit digits give 1216 bits of precision + // to handle the largest possible float64 exponent. + var mPi = [...]uint64{ + 0x0000000000000000, + 0x517cc1b727220a94, + 0xfe13abe8fa9a6ee0, + 0x6db14acc9e21c820, + 0xff28b1d5ef5de2b0, + 0xdb92371d2126e970, + 0x0324977504e8c90e, + 0x7f0ef58e5894d39f, + 0x74411afa975da242, + 0x74ce38135a2fbf20, + 0x9cc8eb1cc1a99cfa, + 0x4e422fc5defc941d, + 0x8ffc4bffef02cc07, + 0xf79788c5ad05368f, + 0xb69b3f6793e584db, + 0xa7a31fb34f2ff516, + 0xba93dd63f5f2f8bd, + 0x9e839cfbc5294975, + 0x35fdafd88fc6ae84, + 0x2b0198237e3db5d5, + } + // Use the exponent to extract the 3 appropriate uint64 digits from mPi, + // B ~ (z0, z1, z2), such that the product leading digit has the exponent -64. + // Note, exp >= 50 since x >= reduceThreshold and exp < 971 for maximum float64. + digit, bitshift := uint(exp+64)/64, uint(exp+64)%64 + z0 := (mPi[digit] << bitshift) | (mPi[digit+1] >> (64 - bitshift)) + z1 := (mPi[digit+1] << bitshift) | (mPi[digit+2] >> (64 - bitshift)) + z2 := (mPi[digit+2] << bitshift) | (mPi[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) + // Find the magnitude of the fraction. + 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 + x = math.Float64frombits(hi) + // map to (-Pi/2, Pi/2] + if x > 0.5 { + x-- } - t = float64(int64(t)) // int64(t) = the multiple - return ((x - t*DP1) - t*DP2) - t*DP3 + return math.Pi * x } // Taylor series expansion for cosh(2y) - cos(2x) diff --git a/libgo/go/math/example_test.go b/libgo/go/math/example_test.go index 50c7426..9fc1967 100644 --- a/libgo/go/math/example_test.go +++ b/libgo/go/math/example_test.go @@ -219,3 +219,22 @@ func ExampleTrunc() { // 3.00 // -1.00 } + +func ExampleCbrt() { + fmt.Printf("%.2f\n", math.Cbrt(8)) + fmt.Printf("%.2f\n", math.Cbrt(27)) + // Output: + // 2.00 + // 3.00 +} + +func ExampleModf() { + int, frac := math.Modf(3.14) + fmt.Printf("%.2f, %.2f\n", int, frac) + + int, frac = math.Modf(-2.71) + fmt.Printf("%.2f, %.2f\n", int, frac) + // Output: + // 3.00, 0.14 + // -2.00, -0.71 +} diff --git a/libgo/go/math/huge_test.go b/libgo/go/math/huge_test.go index 0b45dbf..9448edc 100644 --- a/libgo/go/math/huge_test.go +++ b/libgo/go/math/huge_test.go @@ -16,6 +16,10 @@ import ( // Inputs to test trig_reduce var trigHuge = []float64{ + 1 << 28, + 1 << 29, + 1 << 30, + 1 << 35, 1 << 120, 1 << 240, 1 << 480, @@ -29,6 +33,10 @@ var trigHuge = []float64{ // 102 decimal digits (1 << 120, 1 << 240, 1 << 480, 1234567891234567 << 180) // were confirmed via https://keisan.casio.com/ var cosHuge = []float64{ + -0.16556897949057876, + -0.94517382606089662, + 0.78670712294118812, + -0.76466301249635305, -0.92587902285483787, 0.93601042593353793, -0.28282777640193788, @@ -38,6 +46,10 @@ var cosHuge = []float64{ } var sinHuge = []float64{ + -0.98619821183697566, + 0.32656766301856334, + -0.61732641504604217, + -0.64443035102329113, 0.37782010936075202, -0.35197227524865778, 0.95917070894368716, @@ -47,6 +59,10 @@ var sinHuge = []float64{ } var tanHuge = []float64{ + 5.95641897939639421, + -0.34551069233430392, + -0.78469661331920043, + 0.84276385870875983, -0.40806638884180424, -0.37603456702698076, -3.39135965054779932, diff --git a/libgo/go/math/trig_reduce.go b/libgo/go/math/trig_reduce.go index 6f8eaba..5cdf4fa 100644 --- a/libgo/go/math/trig_reduce.go +++ b/libgo/go/math/trig_reduce.go @@ -8,13 +8,19 @@ 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) +// reduceThreshold is the maximum value of x where the reduction using Pi/4 +// in 3 float64 parts still gives accurate results. This threshold +// is set by y*C being representable as a float64 without error +// where y is given by y = floor(x * (4 / Pi)) and C is the leading partial +// terms of 4/Pi. Since the leading terms (PI4A and PI4B in sin.go) have 30 +// and 32 trailing zero bits, y should have less than 30 significant bits. +// y < 1<<30 -> floor(x*4/Pi) < 1<<30 -> x < (1<<30 - 1) * Pi/4 +// So, conservatively we can take x < 1<<29. +// Above this threshold Payne-Hanek range reduction must be used. +const reduceThreshold = 1 << 29 // trigReduce implements Payne-Hanek range reduction by Pi/4 -// for x > 0. It returns the integer part mod 8 (j) and +// 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" |