Added special condition and benchmarks for Cbrt, Sincos. Took Frexp and Ldexp out of bits.go.
R=rsc
CC=golang-dev
https://golang.org/cl/206084
exp2_386.$O\
fabs_386.$O\
floor_386.$O\
+ frexp_386.$O\
fmod_386.$O\
hypot_386.$O\
+ ldexp_386.$O\
log_386.$O\
log1p_386.$O\
modf_386.$O\
sin_386.$O\
+ sincos_386.$O\
sqrt_386.$O\
tan_386.$O\
atanh.go\
atan2.go\
bits.go\
+ cbrt.go\
const.go\
copysign.go\
erf.go\
fdim.go\
floor.go\
fmod.go\
+ frexp.go\
hypot.go\
+ ldexp.go\
log.go\
log1p.go\
modf.go\
pow.go\
pow10.go\
sin.go\
+ sincos.go\
sinh.go\
sqrt.go\
sqrt_port.go\
1.3902530903455392306872261e+00,
2.2859857424479142655411058e+00,
}
+var cbrt = []float64{
+ 1.7075799841925094446722675e+00,
+ 1.9779982212970353936691498e+00,
+ -6.5177429017779910853339447e-01,
+ -1.7111838886544019873338113e+00,
+ 2.1279920909827937423960472e+00,
+ 1.4303536770460741452312367e+00,
+ 1.7357021059106154902341052e+00,
+ 1.3972633462554328350552916e+00,
+ 1.2221149580905388454977636e+00,
+ -2.0556003730500069110343596e+00,
+}
var ceil = []float64{
5.0000000000000000e+00,
8.0000000000000000e+00,
NaN(),
}
+var vfcbrtSC = []float64{
+ Inf(-1),
+ Inf(1),
+ NaN(),
+}
+var cbrtSC = []float64{
+ Inf(-1),
+ Inf(1),
+ NaN(),
+}
+
var vfceilSC = []float64{
Inf(-1),
Inf(1),
}
}
+func TestCbrt(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ if f := Cbrt(vf[i]); !veryclose(cbrt[i], f) {
+ t.Errorf("Cbrt(%g) = %g, want %g\n", vf[i], f, cbrt[i])
+ }
+ }
+ for i := 0; i < len(vfcbrtSC); i++ {
+ if f := Cbrt(vfcbrtSC[i]); !alike(cbrtSC[i], f) {
+ t.Errorf("Cbrt(%g) = %g, want %g\n", vfcbrtSC[i], f, cbrtSC[i])
+ }
+ }
+}
+
func TestCeil(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := Ceil(vf[i]); ceil[i] != f {
}
}
+func TestSincos(t *testing.T) {
+ for i := 0; i < len(vf); i++ {
+ if s, c := Sincos(vf[i]); !close(sin[i], s) || !close(cos[i], c) {
+ t.Errorf("Sincos(%g) = %g, %g want %g, %g\n", vf[i], s, c, sin[i], cos[i])
+ }
+ }
+}
+
func TestSinh(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := Sinh(vf[i]); !close(sinh[i], f) {
// Check that math functions of high angle values
// return similar results to low angle values
+func TestLargeCos(t *testing.T) {
+ large := float64(100000 * Pi)
+ for i := 0; i < len(vf); i++ {
+ f1 := Cos(vf[i])
+ f2 := Cos(vf[i] + large)
+ if !kindaclose(f1, f2) {
+ t.Errorf("Cos(%g) = %g, want %g\n", vf[i]+large, f2, f1)
+ }
+ }
+}
+
func TestLargeSin(t *testing.T) {
large := float64(100000 * Pi)
for i := 0; i < len(vf); i++ {
}
}
-func TestLargeCos(t *testing.T) {
+func TestLargeSincos(t *testing.T) {
large := float64(100000 * Pi)
for i := 0; i < len(vf); i++ {
- f1 := Cos(vf[i])
- f2 := Cos(vf[i] + large)
- if !kindaclose(f1, f2) {
- t.Errorf("Cos(%g) = %g, want %g\n", vf[i]+large, f2, f1)
+ f1, g1 := Sincos(vf[i])
+ f2, g2 := Sincos(vf[i] + large)
+ if !kindaclose(f1, f2) || !kindaclose(g1, g2) {
+ t.Errorf("Sincos(%g) = %g, %g, want %g, %g\n", vf[i]+large, f2, g2, f1, g1)
}
}
}
}
}
+func BenchmarkCbrt(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Cbrt(10)
+ }
+}
+
func BenchmarkCeil(b *testing.B) {
for i := 0; i < b.N; i++ {
Ceil(.5)
}
}
+func BenchmarkSincos(b *testing.B) {
+ for i := 0; i < b.N; i++ {
+ Sincos(.5)
+ }
+}
+
func BenchmarkSinh(b *testing.B) {
for i := 0; i < b.N; i++ {
Sinh(2.5)
// return sign >= 0 && x == uvinf || sign <= 0 && x == uvneginf;
return sign >= 0 && f > MaxFloat64 || sign <= 0 && f < -MaxFloat64
}
-
-// Frexp breaks f into a normalized fraction
-// and an integral power of two.
-// It returns frac and exp satisfying f == frac × 2<sup>exp</sup>,
-// with the absolute value of frac in the interval [½, 1).
-func Frexp(f float64) (frac float64, exp int) {
- // TODO(rsc): Remove manual inlining of IsNaN, IsInf
- // when compiler does it for us
- // special cases
- switch {
- case f == 0:
- return
- case f < -MaxFloat64 || f > MaxFloat64 || f != f: // IsInf(f, 0) || IsNaN(f):
- frac = f
- return
- }
- x := Float64bits(f)
- exp = int((x>>shift)&mask) - bias
- x &^= mask << shift
- x |= bias << shift
- frac = Float64frombits(x)
- return
-}
-
-// Ldexp is the inverse of Frexp.
-// It returns frac × 2<sup>exp</sup>.
-func Ldexp(frac float64, exp int) float64 {
- // TODO(rsc): Remove manual inlining of IsNaN, IsInf
- // when compiler does it for us
- // special cases
- if frac != frac { // IsNaN(frac)
- return NaN()
- }
- x := Float64bits(frac)
- exp += int(x>>shift) & mask
- if exp <= 0 {
- return 0 // underflow
- }
- if exp >= mask { // overflow
- if frac < 0 {
- return Inf(-1)
- }
- return Inf(1)
- }
- x &^= mask << shift
- x |= uint64(exp) << shift
- return Float64frombits(x)
-}
--- /dev/null
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+/*
+ The algorithm is based in part on "Optimal Partitioning of
+ Newton's Method for Calculating Roots", by Gunter Meinardus
+ and G. D. Taylor, Mathematics of Computation © 1980 American
+ Mathematical Society.
+ (http://www.jstor.org/stable/2006387?seq=9, accessed 11-Feb-2010)
+*/
+
+// Cbrt returns the cube root of its argument.
+//
+// Special cases are:
+// Exp(+Inf) = +Inf
+// Exp(-Inf) = -Inf
+// Exp(NaN) = NaN
+func Cbrt(x float64) float64 {
+ const (
+ A1 = 1.662848358e-01
+ A2 = 1.096040958e+00
+ A3 = 4.105032829e-01
+ A4 = 5.649335816e-01
+ B1 = 2.639607233e-01
+ B2 = 8.699282849e-01
+ B3 = 1.629083358e-01
+ B4 = 2.824667908e-01
+ C1 = 4.190115298e-01
+ C2 = 6.904625373e-01
+ C3 = 6.46502159e-02
+ C4 = 1.412333954e-01
+ )
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
+ // special cases
+ switch {
+ case x != x || x < -MaxFloat64 || x > MaxFloat64: // IsNaN(x) || IsInf(x, 0):
+ return x
+ }
+ sign := false
+ if x < 0 {
+ x = -x
+ sign = true
+ }
+ // Reduce argument
+ f, e := Frexp(x)
+ m := e % 3
+ if m > 0 {
+ m -= 3
+ e -= m // e is multiple of 3
+ }
+ f = Ldexp(f, m) // 0.125 <= f < 1.0
+
+ // Estimate cube root
+ switch m {
+ case 0: // 0.5 <= f < 1.0
+ f = A1*f + A2 - A3/(A4+f)
+ case -1: // 0.25 <= f < 0.5
+ f = B1*f + B2 - B3/(B4+f)
+ default: // 0.125 <= f < 0.25
+ f = C1*f + C2 - C3/(C4+f)
+ }
+ y := Ldexp(f, e/3) // e/3 = exponent of cube root
+
+ // Iterate
+ s := y * y * y
+ t := s + x
+ y *= (t + x) / (s + t)
+ // Reiterate
+ s = (y*y*y - x) / x
+ y -= y * (((14.0/81.0)*s-(2.0/9.0))*s + (1.0 / 3.0)) * s
+ if sign {
+ y = -y
+ }
+ return y
+}
--- /dev/null
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+// Frexp breaks f into a normalized fraction
+// and an integral power of two.
+// It returns frac and exp satisfying f == frac × 2<sup>exp</sup>,
+// with the absolute value of frac in the interval [½, 1).
+func Frexp(f float64) (frac float64, exp int) {
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
+ // special cases
+ switch {
+ case f == 0:
+ return
+ case f < -MaxFloat64 || f > MaxFloat64 || f != f: // IsInf(f, 0) || IsNaN(f):
+ frac = f
+ return
+ }
+ x := Float64bits(f)
+ exp = int((x>>shift)&mask) - bias
+ x &^= mask << shift
+ x |= bias << shift
+ frac = Float64frombits(x)
+ return
+}
--- /dev/null
+// Copyright 2010 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.
+
+// func Frexp(x float64) (f float64, e int)
+TEXT ·Frexp(SB),7,$0
+ FMOVD x+0(FP), F0 // F0=x
+ FXAM
+ FSTSW AX
+ SAHF
+ JNP nan_zero_inf
+ JCS nan_zero_inf
+ FXTRACT // F0=f (0<=f<1), F1=e
+ FMULD $(0.5), F0 // F0=f (0.5<=f<1), F1=e
+ FMOVDP F0, f+8(FP) // F0=e
+ FLD1 // F0=1, F1=e
+ FADDDP F0, F1 // F0=e+1
+ FMOVLP F0, e+16(FP) // (int=int32)
+ RET
+nan_zero_inf:
+ FMOVDP F0, f+8(FP) // F0=e
+ MOVL $0, e+16(FP) // e=0
+ RET
--- /dev/null
+// Copyright 2010 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
+
+func Frexp(x float64) (f float64, e int)
--- /dev/null
+// Copyright 2009 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+// Ldexp is the inverse of Frexp.
+// It returns frac × 2<sup>exp</sup>.
+func Ldexp(frac float64, exp int) float64 {
+ // TODO(rsc): Remove manual inlining of IsNaN, IsInf
+ // when compiler does it for us
+ // special cases
+ if frac != frac { // IsNaN(frac)
+ return NaN()
+ }
+ x := Float64bits(frac)
+ exp += int(x>>shift) & mask
+ if exp <= 0 {
+ return 0 // underflow
+ }
+ if exp >= mask { // overflow
+ if frac < 0 {
+ return Inf(-1)
+ }
+ return Inf(1)
+ }
+ x &^= mask << shift
+ x |= uint64(exp) << shift
+ return Float64frombits(x)
+}
--- /dev/null
+// Copyright 2010 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.
+
+// func Ldexp(f float64, e int) float64
+TEXT ·Ldexp(SB),7,$0
+ FMOVL e+8(FP), F0 // F0=e
+ FMOVD x+0(FP), F0 // F0=x, F1=e
+ FSCALE // F0=x*2**e, F1=e
+ FMOVDP F0, F1 // F0=x*2**e
+ FMOVDP F0, r+12(FP)
+ RET
--- /dev/null
+// Copyright 2010 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
+
+func Ldexp(f float64, e int) float64
--- /dev/null
+// Copyright 2010 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
+
+// Sincos(x) returns Sin(x), Cos(x).
+//
+// Special conditions are:
+// Sincos(+Inf) = NaN, NaN
+// Sincos(-Inf) = NaN, NaN
+// Sincos(NaN) = NaN, NaN
+func Sincos(x float64) (sin, cos float64) { return Sin(x), Cos(x) }
--- /dev/null
+// Copyright 2010 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.
+
+// func Sincos(x float64) (sin, cos float64)
+TEXT ·Sincos(SB),7,$0
+ FMOVD x+0(FP), F0 // F0=x
+ FSINCOS // F0=cos(x), F1=sin(x) if -2**63 < x < 2**63
+ FSTSW AX // AX=status word
+ ANDW $0x0400, AX
+ JNE 4(PC) // jump if x outside range
+ FMOVDP F0, c+16(FP) // F0=sin(x)
+ FMOVDP F0, s+8(FP)
+ RET
+ FLDPI // F0=Pi, F1=x
+ FADDD F0, F0 // F0=2*Pi, F1=x
+ FXCHD F0, F1 // F0=x, F1=2*Pi
+ FPREM1 // F0=reduced_x, F1=2*Pi
+ FSTSW AX // AX=status word
+ ANDW $0x0400, AX
+ JNE -3(PC) // jump if reduction incomplete
+ FMOVDP F0, F1 // F0=reduced_x
+ FSINCOS // F0=cos(reduced_x), F1=sin(reduced_x)
+ FMOVDP F0, c+16(FP) // F0=sin(reduced_x)
+ FMOVDP F0, s+8(FP)
+ RET
--- /dev/null
+// Copyright 2010 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
+
+func Sincos(x float64) (sin, cos float64)