$(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
-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;
}
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]);
}
}
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]);
}
}
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]);
}
}
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]);
}
}
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]);
}
}
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]);
}
}
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]);
}
}
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]);
}
}
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]);
}
}
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);
}
}
* Arctan is called after appropriate range reduction.
*/
-const
-(
- pio2 = .15707963267948966192313216e1
-)
-
export func Asin(arg float64) float64 {
var temp, x float64;
var sign bool;
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);
}
if arg > 1 || arg < -1 {
return sys.NaN();
}
- return pio2 - Asin(arg);
+ return Pi/2 - Asin(arg);
}
package math
+import "math"
+
/*
* floating-point arctangent
*
* 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;
}
* 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));
}
/*
* 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;
}
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;
)
// 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):
// 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
// 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) ;
// 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):
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);
}
* 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];
}
}
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);
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 {
* 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;
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 {
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;
x = -x;
sign = true;
}
- x = x * piu4; /* overflow? */
+ x = x * (4/Pi); /* overflow? */
var e float64;
e, x = sys.modf(x);
i := int32(e);
}
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) {