diff --git a/src/pkg/math/Makefile b/src/pkg/math/Makefile index 5472dc1ace..ad012b64cb 100644 --- a/src/pkg/math/Makefile +++ b/src/pkg/math/Makefile @@ -15,7 +15,6 @@ OFILES_amd64=\ exp_amd64.$O\ hypot_amd64.$O\ log_amd64.$O\ - sincos_amd64.$O\ sqrt_amd64.$O\ OFILES_386=\ diff --git a/src/pkg/math/sincos.go b/src/pkg/math/sincos.go index e8261bca71..f5412fd726 100644 --- a/src/pkg/math/sincos.go +++ b/src/pkg/math/sincos.go @@ -4,9 +4,66 @@ package math +// Coefficients _sin[] and _cos[] are found in pkg/math/sin.go. + // Sincos(x) returns Sin(x), Cos(x). // // Special conditions are: +// Sincos(±0) = ±0, 1 // Sincos(±Inf) = NaN, NaN // Sincos(NaN) = NaN, NaN -func Sincos(x float64) (sin, cos float64) { return Sin(x), Cos(x) } +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 + ) + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + // special cases + switch { + case x == 0: + return x, 1 // return ±0.0, 1.0 + case x != x || x < -MaxFloat64 || x > MaxFloat64: // IsNaN(x) || IsInf(x, 0): + return NaN(), NaN() + } + + // make argument positive + sinSign, cosSign := false, false + if x < 0 { + x = -x + 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 + + if j&1 == 1 { // map zeros to origin + j += 1 + y += 1 + } + j &= 7 // octant modulo 2Pi radians (360 degrees) + if j > 3 { // reflect in x axis + j -= 4 + sinSign, cosSign = !sinSign, !cosSign + } + if j > 1 { + 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]) + if j == 1 || j == 2 { + sin, cos = cos, sin + } + if cosSign { + cos = -cos + } + if sinSign { + sin = -sin + } + return +}