1.8253080916808550e+00,
-8.6859247685756013e+00,
}
+var acos = []float64{
+ 1.0496193546107222e+00,
+ 6.858401281366443e-01,
+ 1.598487871457716e+00,
+ 2.095619936147586e+00,
+ 2.7053008467824158e-01,
+ 1.2738121680361776e+00,
+ 1.0205369421140630e+00,
+ 1.2945003481781246e+00,
+ 1.3872364345374451e+00,
+ 2.6231510803970464e+00,
+}
var asin = []float64{
5.2117697218417440e-01,
8.8495619865825236e-01,
9.4936501296239700e-01,
-9.9999994291374019e-01,
}
-var vfsin = []float64{
- NaN(),
- Inf(-1),
- 0,
- Inf(1),
-}
-var vfasin = []float64{
+
+// arguments and expected results for special cases
+var vfasinSC = []float64{
NaN(),
-Pi,
- 0,
Pi,
}
-var vf1 = []float64{
+var asinSC = []float64{
+ NaN(),
+ NaN(),
NaN(),
- Inf(-1),
- -Pi,
- -1,
- 0,
- 1,
- Pi,
- Inf(1),
}
-var vfhypot = [][2]float64{
- [2]float64{Inf(-1), 1},
- [2]float64{Inf(1), 1},
- [2]float64{1, Inf(-1)},
- [2]float64{1, Inf(1)},
- [2]float64{NaN(), Inf(-1)},
- [2]float64{NaN(), Inf(1)},
- [2]float64{1, NaN()},
- [2]float64{NaN(), 1},
+
+var vfatanSC = []float64{
+ NaN(),
+}
+var atanSC = []float64{
+ NaN(),
}
-var vf2 = [][2]float64{
+var vfpowSC = [][2]float64{
[2]float64{-Pi, Pi},
[2]float64{-Pi, -Pi},
[2]float64{Inf(-1), 3},
[2]float64{Inf(1), 0},
[2]float64{NaN(), 0},
}
-var pow2 = []float64{
+var powSC = []float64{
NaN(),
NaN(),
Inf(-1),
return false
}
+func TestAcos(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ // if f := Acos(vf[i] / 10); !veryclose(acos[i], f) {
+ if f := Acos(vf[i] / 10); !close(acos[i], f) {
+ t.Errorf("Acos(%g) = %g, want %g\n", vf[i]/10, f, acos[i])
+ }
+ }
+ for i := 0; i < len(vfasinSC); i++ {
+ if f := Acos(vfasinSC[i]); !alike(asinSC[i], f) {
+ t.Errorf("Acos(%g) = %g, want %g\n", vfasinSC[i], f, asinSC[i])
+ }
+ }
+}
+
func TestAsin(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := Asin(vf[i] / 10); !veryclose(asin[i], f) {
t.Errorf("Asin(%g) = %g, want %g\n", vf[i]/10, f, asin[i])
}
}
+ for i := 0; i < len(vfasinSC); i++ {
+ if f := Asin(vfasinSC[i]); !alike(asinSC[i], f) {
+ t.Errorf("Asin(%g) = %g, want %g\n", vfasinSC[i], f, asinSC[i])
+ }
+ }
}
func TestAtan(t *testing.T) {
t.Errorf("Atan(%g) = %g, want %g\n", vf[i], f, atan[i])
}
}
+ for i := 0; i < len(vfatanSC); i++ {
+ if f := Atan(vfatanSC[i]); !alike(atanSC[i], f) {
+ t.Errorf("Atan(%g) = %g, want %g\n", vfatanSC[i], f, atanSC[i])
+ }
+ }
}
func TestExp(t *testing.T) {
t.Errorf("Pow(10, %.17g) = %.17g, want %.17g\n", vf[i], f, pow[i])
}
}
- for i := 0; i < len(vf2); i++ {
- if f := Pow(vf2[i][0], vf2[i][1]); !alike(pow2[i], f) {
- t.Errorf("Pow(%.17g, %.17g) = %.17g, want %.17g\n", vf2[i][0], vf2[i][1], f, pow2[i])
+ for i := 0; i < len(vfpowSC); i++ {
+ if f := Pow(vfpowSC[i][0], vfpowSC[i][1]); !alike(powSC[i], f) {
+ t.Errorf("Pow(%.17g, %.17g) = %.17g, want %.17g\n", vfpowSC[i][0], vfpowSC[i][1], f, powSC[i])
}
}
}
f1 := Sin(vf[i])
f2 := Sin(vf[i] + large)
if !kindaclose(f1, f2) {
- t.Errorf("Sin(%g) = %g, want %g\n", vf[i]+large, f1, f2)
+ t.Errorf("Sin(%g) = %g, want %g\n", vf[i]+large, f2, f1)
}
}
}
f1 := Cos(vf[i])
f2 := Cos(vf[i] + large)
if !kindaclose(f1, f2) {
- t.Errorf("Cos(%g) = %g, want %g\n", vf[i]+large, f1, f2)
+ t.Errorf("Cos(%g) = %g, want %g\n", vf[i]+large, f2, f1)
}
}
}
f1 := Tan(vf[i])
f2 := Tan(vf[i] + large)
if !kindaclose(f1, f2) {
- t.Errorf("Tan(%g) = %g, want %g\n", vf[i]+large, f1, f2)
+ t.Errorf("Tan(%g) = %g, want %g\n", vf[i]+large, f2, f1)
}
}
}
Pow(2.5, 1.5)
}
}
+
+func BenchmarkAtan(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Atan(.5)
+ }
+}
+
+func BenchmarkAsin(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Asin(.5)
+ }
+}
+
+func BenchmarkAcos(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Acos(.5)
+ }
+}
package math
-
-/*
- * sqrt returns the square root of its floating
- * point argument. Newton's method.
- *
- * calls frexp
- */
+// The original C code and the long comment below are
+// from FreeBSD's /usr/src/lib/msun/src/e_sqrt.c and
+// came with this notice. The go code is a simplified
+// version of the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+// __ieee754_sqrt(x)
+// Return correctly rounded sqrt.
+// -----------------------------------------
+// | Use the hardware sqrt if you have one |
+// -----------------------------------------
+// Method:
+// Bit by bit method using integer arithmetic. (Slow, but portable)
+// 1. Normalization
+// Scale x to y in [1,4) with even powers of 2:
+// find an integer k such that 1 <= (y=x*2^(2k)) < 4, then
+// sqrt(x) = 2^k * sqrt(y)
+// 2. Bit by bit computation
+// Let q = sqrt(y) truncated to i bit after binary point (q = 1),
+// i 0
+// i+1 2
+// s = 2*q , and y = 2 * ( y - q ). (1)
+// i i i i
+//
+// To compute q from q , one checks whether
+// i+1 i
+//
+// -(i+1) 2
+// (q + 2 ) <= y. (2)
+// i
+// -(i+1)
+// If (2) is false, then q = q ; otherwise q = q + 2 .
+// i+1 i i+1 i
+//
+// With some algebric manipulation, it is not difficult to see
+// that (2) is equivalent to
+// -(i+1)
+// s + 2 <= y (3)
+// i i
+//
+// The advantage of (3) is that s and y can be computed by
+// i i
+// the following recurrence formula:
+// if (3) is false
+//
+// s = s , y = y ; (4)
+// i+1 i i+1 i
+//
+// otherwise,
+// -i -(i+1)
+// s = s + 2 , y = y - s - 2 (5)
+// i+1 i i+1 i i
+//
+// One may easily use induction to prove (4) and (5).
+// Note. Since the left hand side of (3) contain only i+2 bits,
+// it does not necessary to do a full (53-bit) comparison
+// in (3).
+// 3. Final rounding
+// After generating the 53 bits result, we compute one more bit.
+// Together with the remainder, we can decide whether the
+// result is exact, bigger than 1/2ulp, or less than 1/2ulp
+// (it will never equal to 1/2ulp).
+// The rounding mode can be detected by checking whether
+// huge + tiny is equal to huge, and whether huge - tiny is
+// equal to huge for some floating point number "huge" and "tiny".
+//
+//
+// Notes: Rounding mode detection omitted. The constants "mask", "shift",
+// and "bias" are found in src/pkg/math/bits.go
// Sqrt returns the square root of x.
//
// Sqrt(+Inf) = +Inf
// Sqrt(0) = 0
// Sqrt(x < 0) = NaN
+// Sqrt(NaN) = NaN
func Sqrt(x float64) float64 {
- if IsInf(x, 1) {
+ // special cases
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
+ switch {
+ case x != x || x > MaxFloat64: // IsNaN(x) || IsInf(x, 1):
return x
- }
-
- if x <= 0 {
- if x < 0 {
- return NaN()
- }
+ case x == 0:
return 0
+ case x < 0:
+ return NaN()
}
-
- y, exp := Frexp(x)
- for y < 0.5 {
- y = y * 2
- exp = exp - 1
- }
-
- if exp&1 != 0 {
- y = y * 2
- exp = exp - 1
- }
- temp := 0.5 * (1 + y)
-
- for exp > 60 {
- temp = temp * float64(1<<30)
- exp = exp - 60
+ ix := Float64bits(x)
+ // normalize x
+ exp := int((ix >> shift) & mask)
+ if exp == 0 { // subnormal x
+ for ix&1<<shift == 0 {
+ ix <<= 1
+ exp--
+ }
+ exp++
}
- for exp < -60 {
- temp = temp / float64(1<<30)
- exp = exp + 60
+ exp -= bias + 1 // unbias exponent
+ ix &^= mask << shift
+ ix |= 1 << shift
+ if exp&1 == 1 { // odd exp, double x to make it even
+ ix <<= 1
}
- if exp >= 0 {
- exp = 1 << uint(exp/2)
- temp = temp * float64(exp)
- } else {
- exp = 1 << uint(-exp/2)
- temp = temp / float64(exp)
+ exp >>= 1 // exp = exp/2, exponent of square root
+ // generate sqrt(x) bit by bit
+ ix <<= 1
+ var q, s uint64 // q = sqrt(x)
+ r := uint64(1 << (shift + 1)) // r = moving bit from MSB to LSB
+ for r != 0 {
+ t := s + r
+ if t <= ix {
+ s = t + r
+ ix -= t
+ q += r
+ }
+ ix <<= 1
+ r >>= 1
}
-
- for i := 0; i <= 4; i++ {
- temp = 0.5 * (temp + x/temp)
+ // final rounding
+ if ix != 0 { // remainder, result not exact
+ q += q & 1 // round according to extra bit
}
- return temp
+ ix = q>>1 + 0x3fe0000000000000 // q/2 + 0.5
+ ix += uint64(exp) << shift
+ return Float64frombits(ix)
}