Skip to content

G
go
Commit 2c8d9a56 authored by Russ Cox's avatar Russ Cox
Browse files
Options

redo and clean up math. 

R=r
DELTA=243  (60 added, 72 deleted, 111 changed)
OCL=22909
CL=22912
parent 293c8f8c
Hide whitespace changes
Inline Side-by-side
Showing with 168 additions and 180 deletions
src/lib/math/Makefile
View file @ 2c8d9a56
......@@ -32,40 +32,40 @@ coverage: packages
$(AS) $*.s
O1=\
atan.$O\
exp.$O\
fabs.$O\
floor.$O\
fmod.$O\
hypot.$O\
pow10.$O\
sin.$O\
sqrt.$O\
tan.$O\
const.$O\
O2=\
asin.$O\
atan2.$O\
atan.$O\
log.$O\
sin.$O\
sinh.$O\
tan.$O\
O3=\
asin.$O\
atan2.$O\
pow.$O\
tanh.$O\
math.a: a1 a2 a3
a1: $(O1)
$(AR) grc math.a atan.$O exp.$O fabs.$O floor.$O fmod.$O hypot.$O pow10.$O sin.$O sqrt.$O tan.$O const.$O
$(AR) grc math.a exp.$O fabs.$O floor.$O fmod.$O hypot.$O pow10.$O sqrt.$O const.$O
rm -f $(O1)
a2: $(O2)
$(AR) grc math.a asin.$O atan2.$O log.$O sinh.$O
$(AR) grc math.a atan.$O log.$O sin.$O sinh.$O tan.$O
rm -f $(O2)
a3: $(O3)
$(AR) grc math.a pow.$O tanh.$O
$(AR) grc math.a asin.$O atan2.$O pow.$O tanh.$O
rm -f $(O3)
newpkg: clean
......
src/lib/math/all_test.go
View file @ 2c8d9a56
......@@ -154,7 +154,7 @@ var tanh = []float64 {
-9.9999994291374019e-01,
}
func Tolerance(a,b,e float64) bool {
func tolerance(a,b,e float64) bool {
d := a-b;
if d < 0 {
d = -d;
......@@ -168,16 +168,16 @@ func Tolerance(a,b,e float64) bool {
}
return d < e;
}
func Close(a,b float64) bool {
return Tolerance(a, b, 1e-14);
func close(a,b float64) bool {
return tolerance(a, b, 1e-14);
}
func VeryClose(a,b float64) bool {
return Tolerance(a, b, 4e-16);
func veryclose(a,b float64) bool {
return tolerance(a, b, 4e-16);
}
export func TestAsin(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := math.Asin(vf[i]/10); !VeryClose(asin[i], f) {
if f := math.Asin(vf[i]/10); !veryclose(asin[i], f) {
t.Errorf("math.Asin(%g) = %g, want %g\n", vf[i]/10, f, asin[i]);
}
}
......@@ -185,7 +185,7 @@ export func TestAsin(t *testing.T) {
export func TestAtan(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := math.Atan(vf[i]); !VeryClose(atan[i], f) {
if f := math.Atan(vf[i]); !veryclose(atan[i], f) {
t.Errorf("math.Atan(%g) = %g, want %g\n", vf[i], f, atan[i]);
}
}
......@@ -193,7 +193,7 @@ export func TestAtan(t *testing.T) {
export func TestExp(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := math.Exp(vf[i]); !VeryClose(exp[i], f) {
if f := math.Exp(vf[i]); !veryclose(exp[i], f) {
t.Errorf("math.Exp(%g) = %g, want %g\n", vf[i], f, exp[i]);
}
}
......@@ -214,15 +214,14 @@ export func TestLog(t *testing.T) {
t.Errorf("math.Log(%g) = %g, want %g\n", a, f, log[i]);
}
}
const Ln10 = 2.30258509299404568401799145468436421;
if f := math.Log(10); f != Ln10 {
t.Errorf("math.Log(%g) = %g, want %g\n", 10, f, Ln10);
if f := math.Log(10); f != math.Ln10 {
t.Errorf("math.Log(%g) = %g, want %g\n", 10, f, math.Ln10);
}
}
export func TestPow(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := math.Pow(10, vf[i]); !Close(pow[i], f) {
if f := math.Pow(10, vf[i]); !close(pow[i], f) {
t.Errorf("math.Pow(10, %.17g) = %.17g, want %.17g\n", vf[i], f, pow[i]);
}
}
......@@ -230,7 +229,7 @@ export func TestPow(t *testing.T) {
export func TestSin(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := math.Sin(vf[i]); !Close(sin[i], f) {
if f := math.Sin(vf[i]); !close(sin[i], f) {
t.Errorf("math.Sin(%g) = %g, want %g\n", vf[i], f, sin[i]);
}
}
......@@ -238,7 +237,7 @@ export func TestSin(t *testing.T) {
export func TestSinh(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := math.Sinh(vf[i]); !VeryClose(sinh[i], f) {
if f := math.Sinh(vf[i]); !veryclose(sinh[i], f) {
t.Errorf("math.Sinh(%g) = %g, want %g\n", vf[i], f, sinh[i]);
}
}
......@@ -247,7 +246,7 @@ export func TestSinh(t *testing.T) {
export func TestSqrt(t *testing.T) {
for i := 0; i < len(vf); i++ {
a := math.Fabs(vf[i]);
if f := math.Sqrt(a); !VeryClose(sqrt[i], f) {
if f := math.Sqrt(a); !veryclose(sqrt[i], f) {
t.Errorf("math.Sqrt(%g) = %g, want %g\n", a, f, floor[i]);
}
}
......@@ -255,7 +254,7 @@ export func TestSqrt(t *testing.T) {
export func TestTan(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := math.Tan(vf[i]); !Close(tan[i], f) {
if f := math.Tan(vf[i]); !close(tan[i], f) {
t.Errorf("math.Tan(%g) = %g, want %g\n", vf[i], f, tan[i]);
}
}
......@@ -263,7 +262,7 @@ export func TestTan(t *testing.T) {
export func TestTanh(t *testing.T) {
for i := 0; i < len(vf); i++ {
if f := math.Tanh(vf[i]); !VeryClose(tanh[i], f) {
if f := math.Tanh(vf[i]); !veryclose(tanh[i], f) {
t.Errorf("math.Tanh(%g) = %g, want %g\n", vf[i], f, tanh[i]);
}
}
......@@ -272,7 +271,7 @@ export func TestTanh(t *testing.T) {
export func TestHypot(t *testing.T) {
for i := 0; i < len(vf); i++ {
a := math.Fabs(tanh[i]*math.Sqrt(2));
if f := math.Hypot(tanh[i], tanh[i]); !VeryClose(a, f) {
if f := math.Hypot(tanh[i], tanh[i]); !veryclose(a, f) {
t.Errorf("math.Hypot(%g, %g) = %g, want %g\n", tanh[i], tanh[i], f, a);
}
}
......
src/lib/math/asin.go
View file @ 2c8d9a56
......@@ -13,11 +13,6 @@ import "math"
* Arctan is called after appropriate range reduction.
*/
const
(
pio2 = .15707963267948966192313216e1
)
export func Asin(arg float64) float64 {
var temp, x float64;
var sign bool;
......@@ -34,7 +29,7 @@ export func Asin(arg float64) float64 {
temp = Sqrt(1 - x*x);
if x > 0.7 {
temp = pio2 - Atan(temp/x);
temp = Pi/2 - Atan(temp/x);
} else {
temp = Atan(x/temp);
}
......@@ -49,5 +44,5 @@ export func Acos(arg float64) float64 {
if arg > 1 || arg < -1 {
return sys.NaN();
}
return pio2 - Asin(arg);
return Pi/2 - Asin(arg);
}
src/lib/math/atan.go
View file @ 2c8d9a56
......@@ -4,6 +4,8 @@
package math
import "math"
/*
* floating-point arctangent
*
......@@ -13,32 +15,27 @@ package math
* coefficients are #5077 from Hart & Cheney. (19.56D)
*/
const
(
ap4 = .161536412982230228262e2;
ap3 = .26842548195503973794141e3;
ap2 = .11530293515404850115428136e4;
ap1 = .178040631643319697105464587e4;
ap0 = .89678597403663861959987488e3;
aq4 = .5895697050844462222791e2;
aq3 = .536265374031215315104235e3;
aq2 = .16667838148816337184521798e4;
aq1 = .207933497444540981287275926e4;
aq0 = .89678597403663861962481162e3;
apio2 = .15707963267948966192313216e1;
apio4 = .7853981633974483096156608e0;
asq2p1 = .2414213562373095048802e1; // sqrt(2)+1
asq2m1 = .414213562373095048802e0; // sqrt(2)-1
)
/*
* xatan evaluates a series valid in the
* range [-0.414...,+0.414...]. (tan(pi/8))
*/
func xatan(arg float64) float64 {
argsq := arg*arg;
value := ((((ap4*argsq + ap3)*argsq + ap2)*argsq + ap1)*argsq + ap0);
value = value/(((((argsq + aq4)*argsq + aq3)*argsq + aq2)*argsq + aq1)*argsq + aq0);
const
(
P4 = .161536412982230228262e2;
P3 = .26842548195503973794141e3;
P2 = .11530293515404850115428136e4;
P1 = .178040631643319697105464587e4;
P0 = .89678597403663861959987488e3;
Q4 = .5895697050844462222791e2;
Q3 = .536265374031215315104235e3;
Q2 = .16667838148816337184521798e4;
Q1 = .207933497444540981287275926e4;
Q0 = .89678597403663861962481162e3;
)
sq := arg*arg;
value := ((((P4*sq + P3)*sq + P2)*sq + P1)*sq + P0);
value = value/(((((sq + Q4)*sq + Q3)*sq + Q2)*sq + Q1)*sq + Q0);
return value*arg;
}
......@@ -47,13 +44,13 @@ func xatan(arg float64) float64 {
* to the range [0,0.414...] and calls xatan.
*/
func satan(arg float64) float64 {
if arg < asq2m1 {
if arg < Sqrt2 - 1 {
return xatan(arg);
}
if arg > asq2p1 {
return apio2 - xatan(1/arg);
if arg > Sqrt2 + 1 {
return Pi/2 - xatan(1/arg);
}
return apio4 + xatan((arg-1)/(arg+1));
return Pi/4 + xatan((arg-1)/(arg+1));
}
/*
......
src/lib/math/atan2.go
View file @ 2c8d9a56
......@@ -10,26 +10,19 @@ import "math"
* atan2 discovers what quadrant the angle
* is in and calls atan.
*/
const
(
pio2 = .15707963267948966192313216e1;
pi = .3141592653589793238462643383276e1;
)
export func Atan2(arg1, arg2 float64) float64 {
if arg1+arg2 == arg1 {
if arg1 >= 0 {
return pio2;
return Pi/2;
}
return -pio2;
return -Pi/2;
}
x := Atan(arg1/arg2);
if arg2 < 0 {
if x <= 0 {
return x + pi;
return x + Pi;
}
return x - pi;
return x - Pi;
}
return x;
}
src/lib/math/const.go
View file @ 2c8d9a56
......@@ -5,5 +5,20 @@
package math
export const (
Sqrt2 = 1.41421356237309504880168872420969808;
// Mathematical constants.
// Reference: http://www.research.att.com/~njas/sequences/Axxxxxx
E = 2.71828182845904523536028747135266249775724709369995957496696763; // A001113
Pi = 3.14159265358979323846264338327950288419716939937510582097494459; // A000796
Phi = 1.61803398874989484820458683436563811772030917980576286213544862; // A001622
Sqrt2 = 1.41421356237309504880168872420969807856967187537694807317667974; // A002193
SqrtE = 1.64872127070012814684865078781416357165377610071014801157507931; // A019774
SqrtPi = 1.77245385090551602729816748334114518279754945612238712821380779; // A002161
SqrtPhi = 1.27201964951406896425242246173749149171560804184009624861664038; // A139339
Ln2 = 0.693147180559945309417232121458176568075500134360255254120680009; // A002162
Log2E = 1/Ln2;
Ln10 = 2.30258509299404568401799145468436420760110148862877297603332790; // A002392
Log10E = 1/Ln10;
)
src/lib/math/exp.go
View file @ 2c8d9a56
......@@ -82,26 +82,23 @@ import "math"
// compiler will convert from decimal to binary accurately enough
// to produce the hexadecimal values shown.
export const (
Ln2 = 0.693147180559945309417232121458176568;
HalfLn2 = 0.346573590279972654708616060729088284;
Ln2Hi = 6.93147180369123816490e-01;
Ln2Lo = 1.90821492927058770002e-10;
Log2e = 1.44269504088896338700e+00;
export func Exp(x float64) float64 {
const (
Ln2Hi = 6.93147180369123816490e-01;
Ln2Lo = 1.90821492927058770002e-10;
Log2e = 1.44269504088896338700e+00;
P1 = 1.66666666666666019037e-01; /* 0x3FC55555; 0x5555553E */
P2 = -2.77777777770155933842e-03; /* 0xBF66C16C; 0x16BEBD93 */
P3 = 6.61375632143793436117e-05; /* 0x3F11566A; 0xAF25DE2C */
P4 = -1.65339022054652515390e-06; /* 0xBEBBBD41; 0xC5D26BF1 */
P5 = 4.13813679705723846039e-08; /* 0x3E663769; 0x72BEA4D0 */
P1 = 1.66666666666666019037e-01; /* 0x3FC55555; 0x5555553E */
P2 = -2.77777777770155933842e-03; /* 0xBF66C16C; 0x16BEBD93 */
P3 = 6.61375632143793436117e-05; /* 0x3F11566A; 0xAF25DE2C */
P4 = -1.65339022054652515390e-06; /* 0xBEBBBD41; 0xC5D26BF1 */
P5 = 4.13813679705723846039e-08; /* 0x3E663769; 0x72BEA4D0 */
Overflow = 7.09782712893383973096e+02;
Underflow = -7.45133219101941108420e+02;
NearZero = 1.0/(1<<28); // 2^-28
)
Overflow = 7.09782712893383973096e+02;
Underflow = -7.45133219101941108420e+02;
NearZero = 1.0/(1<<28); // 2^-28
)
export func Exp(x float64) float64 {
// special cases
switch {
case sys.isNaN(x) || sys.isInf(x, 1):
......
src/lib/math/log.go
View file @ 2c8d9a56
......@@ -37,11 +37,11 @@ import "math"
// of this polynomial approximation is bounded by 2**-58.45. In
// other words,
// 2 4 6 8 10 12 14
// R(z) ~ lg1*s +lg2*s +lg3*s +lg4*s +lg5*s +lg6*s +lg7*s
// (the values of lg1 to lg7 are listed in the program)
// R(z) ~ L1*s +L2*s +L3*s +L4*s +L5*s +L6*s +L7*s
// (the values of L1 to L7 are listed in the program)
// and
// | 2 14 | -58.45
// | lg1*s +...+lg7*s - R(z) | <= 2
// | L1*s +...+L7*s - R(z) | <= 2
// | |
// Note that 2s = f - s*f = f - hfsq + s*hfsq, where hfsq = f*f/2.
// In order to guarantee error in log below 1ulp, we compute log
......@@ -49,11 +49,11 @@ import "math"
// log(1+f) = f - s*(f - R) (if f is not too large)
// log(1+f) = f - (hfsq - s*(hfsq+R)). (better accuracy)
//
// 3. Finally, log(x) = k*ln2 + log(1+f).
// = k*ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*ln2_lo)))
// Here ln2 is split into two floating point number:
// ln2_hi + ln2_lo,
// where n*ln2_hi is always exact for |n| < 2000.
// 3. Finally, log(x) = k*Ln2 + log(1+f).
// = k*Ln2_hi+(f-(hfsq-(s*(hfsq+R)+k*Ln2_lo)))
// Here Ln2 is split into two floating point number:
// Ln2_hi + Ln2_lo,
// where n*Ln2_hi is always exact for |n| < 2000.
//
// Special cases:
// log(x) is NaN with signal if x < 0 (including -INF) ;
......@@ -70,19 +70,19 @@ import "math"
// compiler will convert from decimal to binary accurately enough
// to produce the hexadecimal values shown.
const (
ln2Hi = 6.93147180369123816490e-01; /* 3fe62e42 fee00000 */
ln2Lo = 1.90821492927058770002e-10; /* 3dea39ef 35793c76 */
lg1 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
lg2 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
lg3 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
lg4 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
lg5 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
lg6 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
lg7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
)
export func Log(x float64) float64 {
const (
Ln2Hi = 6.93147180369123816490e-01; /* 3fe62e42 fee00000 */
Ln2Lo = 1.90821492927058770002e-10; /* 3dea39ef 35793c76 */
L1 = 6.666666666666735130e-01; /* 3FE55555 55555593 */
L2 = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
L3 = 2.857142874366239149e-01; /* 3FD24924 94229359 */
L4 = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
L5 = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
L6 = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
L7 = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
)
// special cases
switch {
case sys.isNaN(x) || sys.isInf(x, 1):
......@@ -106,23 +106,18 @@ export func Log(x float64) float64 {
s := f/(2+f);
s2 := s*s;
s4 := s2*s2;
t1 := s2*(lg1 + s4*(lg3 + s4*(lg5 + s4*lg7)));
t2 := s4*(lg2 + s4*(lg4 + s4*lg6));
t1 := s2*(L1 + s4*(L3 + s4*(L5 + s4*L7)));
t2 := s4*(L2 + s4*(L4 + s4*L6));
R := t1 + t2;
hfsq := 0.5*f*f;
return k*ln2Hi - ((hfsq-(s*(hfsq+R)+k*ln2Lo)) - f);
return k*Ln2Hi - ((hfsq-(s*(hfsq+R)+k*Ln2Lo)) - f);
}
const
(
ln10u1 = .4342944819032518276511;
)
export func Log10(arg float64) float64 {
if arg <= 0 {
return sys.NaN();
}
return Log(arg) * ln10u1;
return Log(arg) * (1/Ln10);
}
src/lib/math/pow10.go
View file @ 2c8d9a56
......@@ -13,25 +13,24 @@ package math
* than multipication of lower powers of 10.
*/
const tabsize = 70;
var tab[tabsize] float64;
var pow10tab [70]float64;
export func Pow10(e int) float64 {
if e < 0 {
return 1/Pow10(-e);
}
if e < tabsize {
return tab[e];
if e < len(pow10tab) {
return pow10tab[e];
}
m := e/2;
return Pow10(m) * Pow10(e-m);
}
func init() {
tab[0] = 1.0e0;
tab[1] = 1.0e1;
for i:=2; i<tabsize; i++ {
pow10tab[0] = 1.0e0;
pow10tab[1] = 1.0e1;
for i:=2; i<len(pow10tab); i++ {
m := i/2;
tab[i] = tab[m] * tab[i-m];
pow10tab[i] = pow10tab[m] * pow10tab[i-m];
}
}
src/lib/math/sin.go
View file @ 2c8d9a56
......@@ -4,37 +4,34 @@
package math
/*
Coefficients are #3370 from Hart & Cheney (18.80D).
*/
const
(
sp0 = .1357884097877375669092680e8;
sp1 = -.4942908100902844161158627e7;
sp2 = .4401030535375266501944918e6;
sp3 = -.1384727249982452873054457e5;
sp4 = .1459688406665768722226959e3;
sq0 = .8644558652922534429915149e7;
sq1 = .4081792252343299749395779e6;
sq2 = .9463096101538208180571257e4;
sq3 = .1326534908786136358911494e3;
spiu2 = .6366197723675813430755350e0; // 2/pi
)
import "math"
func sinus(arg float64, quad int) float64 {
// Coefficients are #3370 from Hart & Cheney (18.80D).
const
(
P0 = .1357884097877375669092680e8;
P1 = -.4942908100902844161158627e7;
P2 = .4401030535375266501944918e6;
P3 = -.1384727249982452873054457e5;
P4 = .1459688406665768722226959e3;
Q0 = .8644558652922534429915149e7;
Q1 = .4081792252343299749395779e6;
Q2 = .9463096101538208180571257e4;
Q3 = .1326534908786136358911494e3;
)
x := arg;
if(x < 0) {
x = -x;
quad = quad+2;
}
x = x * spiu2; /* underflow? */
x = x * (2/Pi); /* underflow? */
var y float64;
if x > 32764 {
var e float64;
e, y = sys.modf(x);
e = e + float64(quad);
temsp1, f := sys.modf(0.25*e);
temp1, f := sys.modf(0.25*e);
quad = int(e - 4*f);
} else {
k := int32(x);
......@@ -49,10 +46,10 @@ func sinus(arg float64, quad int) float64 {
y = -y;
}
ysq := y*y;
temsp1 := ((((sp4*ysq+sp3)*ysq+sp2)*ysq+sp1)*ysq+sp0)*y;
temsp2 := ((((ysq+sq3)*ysq+sq2)*ysq+sq1)*ysq+sq0);
return temsp1/temsp2;
yy := y*y;
temp1 := ((((P4*yy+P3)*yy+P2)*yy+P1)*yy+P0)*y;
temp2 := ((((yy+Q3)*yy+Q2)*yy+Q1)*yy+Q0);
return temp1/temp2;
}
export func Cos(arg float64) float64 {
......
src/lib/math/sinh.go
View file @ 2c8d9a56
......@@ -14,25 +14,25 @@ import "math"
* greater in magnitude than 0.5.
*
* A series is used for arguments smaller in magnitude than 0.5.
* The coefficients are #2029 from Hart & Cheney. (20.36D)
*
* cosh(arg) is computed from the exponential func for
* all arguments.
*/
const
(
shp0 = -0.6307673640497716991184787251e+6;
shp1 = -0.8991272022039509355398013511e+5;
shp2 = -0.2894211355989563807284660366e+4;
shp3 = -0.2630563213397497062819489e+2;
shq0 = -0.6307673640497716991212077277e+6;
shq1 = 0.1521517378790019070696485176e+5;
shq2 = -0.173678953558233699533450911e+3;
)
export func Sinh(arg float64) float64 {
sign := false;
// The coefficients are #2029 from Hart & Cheney. (20.36D)
const
(
P0 = -0.6307673640497716991184787251e+6;
P1 = -0.8991272022039509355398013511e+5;
P2 = -0.2894211355989563807284660366e+4;
P3 = -0.2630563213397497062819489e+2;
Q0 = -0.6307673640497716991212077277e+6;
Q1 = 0.1521517378790019070696485176e+5;
Q2 = -0.173678953558233699533450911e+3;
)
sign := false;
if arg < 0 {
arg = -arg;
sign = true;
......@@ -47,9 +47,9 @@ export func Sinh(arg float64) float64 {
temp = (Exp(arg) - Exp(-arg))/2;
default:
argsq := arg*arg;
temp = (((shp3*argsq+shp2)*argsq+shp1)*argsq+shp0)*arg;
temp = temp/(((argsq+shq2)*argsq+shq1)*argsq+shq0);
sq := arg*arg;
temp = (((P3*sq+P2)*sq+P1)*sq+P0)*arg;
temp = temp/(((sq+Q2)*sq+Q1)*sq+Q0);
}
if sign {
......
src/lib/math/tan.go
View file @ 2c8d9a56
......@@ -4,25 +4,26 @@
package math
import "math"
/*
* floating point tangent
* Coefficients are #4285 from Hart & Cheney. (19.74D)
*/
const
(
p0 = -.1306820264754825668269611177e+5;
p1 = .1055970901714953193602353981e+4;
p2 = -.1550685653483266376941705728e+2;
p3 = .3422554387241003435328470489e-1;
p4 = .3386638642677172096076369e-4;
q0 = -.1663895238947119001851464661e+5;
q1 = .4765751362916483698926655581e+4;
q2 = -.1555033164031709966900124574e+3;
piu4 = .1273239544735162686151070107e+1; // 4/pi
)
export func Tan(arg float64) float64 {
// Coefficients are #4285 from Hart & Cheney. (19.74D)
const
(
P0 = -.1306820264754825668269611177e+5;
P1 = .1055970901714953193602353981e+4;
P2 = -.1550685653483266376941705728e+2;
P3 = .3422554387241003435328470489e-1;
P4 = .3386638642677172096076369e-4;
Q0 = -.1663895238947119001851464661e+5;
Q1 = .4765751362916483698926655581e+4;
Q2 = -.1555033164031709966900124574e+3;
)
flag := false;
sign := false;
x := arg;
......@@ -30,7 +31,7 @@ export func Tan(arg float64) float64 {
x = -x;
sign = true;
}
x = x * piu4; /* overflow? */
x = x * (4/Pi); /* overflow? */
var e float64;
e, x = sys.modf(x);
i := int32(e);
......@@ -50,8 +51,8 @@ export func Tan(arg float64) float64 {
}
xsq := x*x;
temp := ((((p4*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0)*x;
temp = temp/(((xsq+q2)*xsq+q1)*xsq+q0);
temp := ((((P4*xsq+P3)*xsq+P2)*xsq+P1)*xsq+P0)*x;
temp = temp/(((xsq+Q2)*xsq+Q1)*xsq+Q0);
if flag {
if(temp == 0) {
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment
GitLab Nexedi Edition | About GitLab | About Nexedi | 沪ICP备2021021310号-2 | 沪ICP备2021021310号-7