diff options
-rw-r--r-- | src/lib/math/Makefile | 16 | ||||
-rw-r--r-- | src/lib/math/all_test.go | 35 | ||||
-rw-r--r-- | src/lib/math/asin.go | 9 | ||||
-rw-r--r-- | src/lib/math/atan.go | 47 | ||||
-rw-r--r-- | src/lib/math/atan2.go | 15 | ||||
-rw-r--r-- | src/lib/math/const.go | 17 | ||||
-rw-r--r-- | src/lib/math/exp.go | 31 | ||||
-rw-r--r-- | src/lib/math/log.go | 53 | ||||
-rw-r--r-- | src/lib/math/pow10.go | 15 | ||||
-rw-r--r-- | src/lib/math/sin.go | 43 | ||||
-rw-r--r-- | src/lib/math/sinh.go | 32 | ||||
-rw-r--r-- | src/lib/math/tan.go | 35 |
12 files changed, 168 insertions, 180 deletions
diff --git a/src/lib/math/Makefile b/src/lib/math/Makefile index ef31b174a..d0debdbbc 100644 --- a/src/lib/math/Makefile +++ b/src/lib/math/Makefile @@ -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 diff --git a/src/lib/math/all_test.go b/src/lib/math/all_test.go index ddcb1e3ac..fa0cc7c86 100644 --- a/src/lib/math/all_test.go +++ b/src/lib/math/all_test.go @@ -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); } } diff --git a/src/lib/math/asin.go b/src/lib/math/asin.go index effaba85e..e866d95c6 100644 --- a/src/lib/math/asin.go +++ b/src/lib/math/asin.go @@ -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); } diff --git a/src/lib/math/atan.go b/src/lib/math/atan.go index f2fd7ed33..730d50a51 100644 --- a/src/lib/math/atan.go +++ b/src/lib/math/atan.go @@ -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)); } /* diff --git a/src/lib/math/atan2.go b/src/lib/math/atan2.go index f0b384258..f8c00aa24 100644 --- a/src/lib/math/atan2.go +++ b/src/lib/math/atan2.go @@ -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; } diff --git a/src/lib/math/const.go b/src/lib/math/const.go index a1c5f8e1a..927d6daf2 100644 --- a/src/lib/math/const.go +++ b/src/lib/math/const.go @@ -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; ) diff --git a/src/lib/math/exp.go b/src/lib/math/exp.go index e1402f02a..efc17aec0 100644 --- a/src/lib/math/exp.go +++ b/src/lib/math/exp.go @@ -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): diff --git a/src/lib/math/log.go b/src/lib/math/log.go index 9c54858b4..5eee5daec 100644 --- a/src/lib/math/log.go +++ b/src/lib/math/log.go @@ -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); } diff --git a/src/lib/math/pow10.go b/src/lib/math/pow10.go index e1e9c2e05..1d0a0e0cc 100644 --- a/src/lib/math/pow10.go +++ b/src/lib/math/pow10.go @@ -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]; } } diff --git a/src/lib/math/sin.go b/src/lib/math/sin.go index 077506b30..6e3fbb28f 100644 --- a/src/lib/math/sin.go +++ b/src/lib/math/sin.go @@ -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 { diff --git a/src/lib/math/sinh.go b/src/lib/math/sinh.go index 622467ed7..e0201eef3 100644 --- a/src/lib/math/sinh.go +++ b/src/lib/math/sinh.go @@ -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 { diff --git a/src/lib/math/tan.go b/src/lib/math/tan.go index fff36c67b..ddfe80d0f 100644 --- a/src/lib/math/tan.go +++ b/src/lib/math/tan.go @@ -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) { |