From 07e77936d535e58b0aeb4f2a11400c1050556739 Mon Sep 17 00:00:00 2001 From: Nicholas Noll Date: Sun, 5 Dec 2021 16:16:21 -0800 Subject: Feat: added math library Used Plan9's libc as starting point. This cleans up dangling references due to loss of libc. --- src/base/math/abs.c | 18 +++++++ src/base/math/asin.c | 40 ++++++++++++++++ src/base/math/atan.c | 76 ++++++++++++++++++++++++++++++ src/base/math/atan2.c | 24 ++++++++++ src/base/math/copysign.c | 12 +++++ src/base/math/exp.c | 40 ++++++++++++++++ src/base/math/fabs.c | 10 ++++ src/base/math/floor.c | 27 +++++++++++ src/base/math/fmod.c | 31 ++++++++++++ src/base/math/frexp.c | 120 +++++++++++++++++++++++++++++++++++++++++++++++ src/base/math/nan.c | 54 +++++++++++++++++++++ src/base/math/pow.c | 69 +++++++++++++++++++++++++++ src/base/math/pow10.c | 48 +++++++++++++++++++ src/base/math/rules.mk | 1 + src/base/math/sin.c | 67 ++++++++++++++++++++++++++ src/base/math/sinh.c | 62 ++++++++++++++++++++++++ src/base/math/sqrt.c | 54 +++++++++++++++++++++ src/base/math/tan.c | 67 ++++++++++++++++++++++++++ src/base/math/tanh.c | 25 ++++++++++ 19 files changed, 845 insertions(+) create mode 100644 src/base/math/abs.c create mode 100644 src/base/math/asin.c create mode 100644 src/base/math/atan.c create mode 100644 src/base/math/atan2.c create mode 100644 src/base/math/copysign.c create mode 100644 src/base/math/exp.c create mode 100644 src/base/math/fabs.c create mode 100644 src/base/math/floor.c create mode 100644 src/base/math/fmod.c create mode 100644 src/base/math/frexp.c create mode 100644 src/base/math/nan.c create mode 100644 src/base/math/pow.c create mode 100644 src/base/math/pow10.c create mode 100644 src/base/math/rules.mk create mode 100644 src/base/math/sin.c create mode 100644 src/base/math/sinh.c create mode 100644 src/base/math/sqrt.c create mode 100644 src/base/math/tan.c create mode 100644 src/base/math/tanh.c (limited to 'src/base/math') diff --git a/src/base/math/abs.c b/src/base/math/abs.c new file mode 100644 index 0000000..b758653 --- /dev/null +++ b/src/base/math/abs.c @@ -0,0 +1,18 @@ +#include +#include + +int +math·abs(int a) +{ + if(a < 0) + return -a; + return a; +} + +long +math·labs(long a) +{ + if(a < 0) + return -a; + return a; +} diff --git a/src/base/math/asin.c b/src/base/math/asin.c new file mode 100644 index 0000000..4687715 --- /dev/null +++ b/src/base/math/asin.c @@ -0,0 +1,40 @@ +/* + * asin(arg) and acos(arg) return the arcsin, arccos, + * respectively of their arguments. + * + * Arctan is called after appropriate range reduction. + */ + +#include +#include + +double +math·asin(double arg) +{ + double temp; + int sign; + + sign = 0; + if(arg < 0) { + arg = -arg; + sign++; + } + if(arg > 1) + return math·NaN(); + temp = math·sqrt(1 - arg*arg); + if(arg > 0.7) + temp = math·PIO2 - math·atan(temp/arg); + else + temp = math·atan(arg/temp); + if(sign) + temp = -temp; + return temp; +} + +double +math·acos(double arg) +{ + if(arg > 1 || arg < -1) + return math·NaN(); + return math·PIO2 - math·asin(arg); +} diff --git a/src/base/math/atan.c b/src/base/math/atan.c new file mode 100644 index 0000000..3108724 --- /dev/null +++ b/src/base/math/atan.c @@ -0,0 +1,76 @@ +/* + floating-point arctangent + + atan returns the value of the arctangent of its + argument in the range [-pi/2,pi/2]. + + atan2 returns the arctangent of arg1/arg2 + in the range [-pi,pi]. + + there are no error returns. + + coefficients are #5077 from Hart & Cheney. (19.56D) +*/ + +#include +#include + +#define sq2p1 2.414213562373095048802e0 +#define sq2m1 .414213562373095048802e0 +#define p4 .161536412982230228262e2 +#define p3 .26842548195503973794141e3 +#define p2 .11530293515404850115428136e4 +#define p1 .178040631643319697105464587e4 +#define p0 .89678597403663861959987488e3 +#define q4 .5895697050844462222791e2 +#define q3 .536265374031215315104235e3 +#define q2 .16667838148816337184521798e4 +#define q1 .207933497444540981287275926e4 +#define q0 .89678597403663861962481162e3 + + +/* + xatan evaluates a series valid in the + range [-0.414...,+0.414...]. (tan(pi/8)) + */ + +static double +xatan(double arg) +{ + double argsq, value; + + argsq = arg*arg; + value = ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0); + value = value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0); + return value*arg; +} + +/* + satan reduces its argument (known to be positive) + to the range [0,0.414...] and calls xatan. + */ + +static double +satan(double arg) +{ + + if(arg < sq2m1) + return xatan(arg); + if(arg > sq2p1) + return math·PIO2 - xatan(1/arg); + return math·PIO2/2 + xatan((arg-1)/(arg+1)); +} + +/* + atan makes its argument positive and + calls the inner routine satan. + */ + +double +math·atan(double arg) +{ + + if(arg > 0) + return satan(arg); + return -satan(-arg); +} diff --git a/src/base/math/atan2.c b/src/base/math/atan2.c new file mode 100644 index 0000000..a681a80 --- /dev/null +++ b/src/base/math/atan2.c @@ -0,0 +1,24 @@ +#include +#include +/* + atan2 discovers what quadrant the angle + is in and calls atan. +*/ + +double +math·atan2(double arg1, double arg2) +{ + + if(arg1+arg2 == arg1) { + if(arg1 >= 0) + return math·PIO2; + return -math·PIO2; + } + arg1 = math·atan(arg1/arg2); + if(arg2 < 0) { + if(arg1 <= 0) + return arg1 + math·PI; + return arg1 - math·PI; + } + return arg1; +} diff --git a/src/base/math/copysign.c b/src/base/math/copysign.c new file mode 100644 index 0000000..0e7540b --- /dev/null +++ b/src/base/math/copysign.c @@ -0,0 +1,12 @@ +#include +#include + +double +math·copysign(double x, double y) +{ + union {double f; uint64 i;} ux={x}, uy={y}; + ux.i &= -1ULL/2; + ux.i |= uy.i & 1ULL<<63; + + return ux.f; +} diff --git a/src/base/math/exp.c b/src/base/math/exp.c new file mode 100644 index 0000000..c062d3b --- /dev/null +++ b/src/base/math/exp.c @@ -0,0 +1,40 @@ +/* + exp returns the exponential function of its + floating-point argument. + + The coefficients are #1069 from Hart and Cheney. (22.35D) +*/ + +#include +#include + +#define p0 .2080384346694663001443843411e7 +#define p1 .3028697169744036299076048876e5 +#define p2 .6061485330061080841615584556e2 +#define q0 .6002720360238832528230907598e7 +#define q1 .3277251518082914423057964422e6 +#define q2 .1749287689093076403844945335e4 +#define log2e 1.4426950408889634073599247 +#define sqrt2 1.4142135623730950488016887 +#define maxf 10000 + +double +math·exp(double arg) +{ + double fract, temp1, temp2, xsq; + int ent; + + if(arg == 0) + return 1; + if(arg < -maxf) + return 0; + if(arg > maxf) + return math·Inf(1); + arg *= log2e; + ent = math·floor(arg); + fract = (arg-ent) - 0.5; + xsq = fract*fract; + temp1 = ((p2*xsq+p1)*xsq+p0)*fract; + temp2 = ((xsq+q2)*xsq+q1)*xsq + q0; + return math·ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent); +} diff --git a/src/base/math/fabs.c b/src/base/math/fabs.c new file mode 100644 index 0000000..feef84e --- /dev/null +++ b/src/base/math/fabs.c @@ -0,0 +1,10 @@ +#include +#include + +double +math·fabs(double arg) +{ + if(arg < 0) + return -arg; + return arg; +} diff --git a/src/base/math/floor.c b/src/base/math/floor.c new file mode 100644 index 0000000..60c5ad0 --- /dev/null +++ b/src/base/math/floor.c @@ -0,0 +1,27 @@ +#include +#include +/* + * floor and ceil-- greatest integer <= arg + * (resp least >=) + */ + +double +math·floor(double d) +{ + double fract; + + if(d < 0) { + fract = math·modf(-d, &d); + if(fract != 0.0) + d += 1; + d = -d; + } else + math·modf(d, &d); + return d; +} + +double +math·ceil(double d) +{ + return -math·floor(-d); +} diff --git a/src/base/math/fmod.c b/src/base/math/fmod.c new file mode 100644 index 0000000..d724a8c --- /dev/null +++ b/src/base/math/fmod.c @@ -0,0 +1,31 @@ +#include +#include + +/* + * floating-point mod function without infinity or NaN checking + */ +double +math·fmod(double x, double y) +{ + int sign, yexp, rexp; + double r, yfr, rfr; + + if (y == 0) + return x; + if (y < 0) + y = -y; + yfr = math·frexp(y, &yexp); + sign = 0; + if(x < 0) { + r = -x; + sign++; + } else + r = x; + while(r >= y) { + rfr = math·frexp(r, &rexp); + r -= math·ldexp(y, rexp - yexp - (rfr < yfr)); + } + if(sign) + r = -r; + return r; +} diff --git a/src/base/math/frexp.c b/src/base/math/frexp.c new file mode 100644 index 0000000..7d119bb --- /dev/null +++ b/src/base/math/frexp.c @@ -0,0 +1,120 @@ +#include +#include + +/* + * this is big/little endian non-portable + * it gets the endian from the FPdbleword + * union in u.h. + */ +#define MASK 0x7ffL +#define SHIFT 20 +#define BIAS 1022L +#define SIG 52 + +double +math·frexp(double d, int *ep) +{ + sys·DoubleWord x, a; + + *ep = 0; + /* order matters: only isNaN can operate on NaN */ + if(math·isNaN(d) || math·isInf(d, 0) || d == 0) + return d; + x.x = d; + a.x = math·fabs(d); + if((a.hi >> SHIFT) == 0){ /* normalize subnormal numbers */ + x.x = (double)(1ULL<> SHIFT) & MASK) - BIAS; + x.hi &= ~(MASK << SHIFT); + x.hi |= BIAS << SHIFT; + return x.x; +} + +double +math·ldexp(double d, int deltae) +{ + int e, bits; + sys·DoubleWord x; + ulong z; + + if(d == 0) + return 0; + x.x = d; + e = (x.hi >> SHIFT) & MASK; + if(deltae >= 0 || e+deltae >= 1){ /* no underflow */ + e += deltae; + if(e >= MASK){ /* overflow */ + if(d < 0) + return math·Inf(-1); + return math·Inf(1); + } + }else{ /* underflow gracefully */ + deltae = -deltae; + /* need to shift d right deltae */ + if(e > 1){ /* shift e-1 by exponent manipulation */ + deltae -= e-1; + e = 1; + } + if(deltae > 0 && e==1){ /* shift 1 by switch from 1.fff to 0.1ff */ + deltae--; + e = 0; + x.lo >>= 1; + x.lo |= (x.hi&1)<<31; + z = x.hi & ((1<>1); + } + while(deltae > 0){ /* shift bits down */ + bits = deltae; + if(bits > SHIFT) + bits = SHIFT; + x.lo >>= bits; + x.lo |= (x.hi&((1<>bits; + deltae -= bits; + } + } + x.hi &= ~(MASK << SHIFT); + x.hi |= (long)e << SHIFT; + return x.x; +} + +double +math·modf(double d, double *ip) +{ + sys·DoubleWord x; + int e; + + x.x = d; + e = (x.hi >> SHIFT) & MASK; + if(e == MASK){ + *ip = d; + if(x.lo != 0 || (x.hi & 0xfffffL) != 0) /* NaN */ + return d; + /* ±Inf */ + x.hi &= 0x80000000L; + return x.x; + } + if(d < 1) { + if(d < 0) { + x.x = math·modf(-d, ip); + *ip = -*ip; + return -x.x; + } + *ip = 0; + return d; + } + e -= BIAS; + if(e <= SHIFT+1) { + x.hi &= ~(0x1fffffL >> e); + x.lo = 0; + } else + if(e <= SHIFT+33) + x.lo &= ~(0x7fffffffL >> (e-SHIFT-2)); + *ip = x.x; + return d - x.x; +} diff --git a/src/base/math/nan.c b/src/base/math/nan.c new file mode 100644 index 0000000..c1f85c2 --- /dev/null +++ b/src/base/math/nan.c @@ -0,0 +1,54 @@ +#include +#include + +#define NANEXP (2047<<20) +#define NANMASK (2047<<20) +#define NANSIGN (1<<31) + +double +math·NaN(void) +{ + sys·DoubleWord a; + + a.hi = NANEXP; + a.lo = 1; + return a.x; +} + +int +math·isNaN(double d) +{ + sys·DoubleWord a; + + a.x = d; + if((a.hi & NANMASK) != NANEXP) + return 0; + return !math·isInf(d, 0); +} + +double +math·Inf(int sign) +{ + sys·DoubleWord a; + + a.hi = NANEXP; + a.lo = 0; + if(sign < 0) + a.hi |= NANSIGN; + return a.x; +} + +int +math·isInf(double d, int sign) +{ + sys·DoubleWord a; + + a.x = d; + if(a.lo != 0) + return 0; + if(a.hi == NANEXP) + return sign >= 0; + if(a.hi == (NANEXP|NANSIGN)) + return sign <= 0; + return 0; +} diff --git a/src/base/math/pow.c b/src/base/math/pow.c new file mode 100644 index 0000000..912a4e1 --- /dev/null +++ b/src/base/math/pow.c @@ -0,0 +1,69 @@ +#include +#include + +double +math·pow(double x, double y) /* return x ^ y (exponentiation) */ +{ + double xy, y1, ye; + long i; + int ex, ey, flip; + + if(y == 0.0) + return 1.0; + + flip = 0; + if(y < 0.){ + y = -y; + flip = 1; + } + y1 = math·modf(y, &ye); + if(y1 != 0.0){ + if(x <= 0.) + goto zreturn; + if(y1 > 0.5) { + y1 -= 1.; + ye += 1.; + } + xy = math·exp(y1 * math·log(x)); + }else + xy = 1.0; + if(ye > 0x7FFFFFFF){ /* should be ~0UL but compiler can't convert double to ulong */ + if(x <= 0.){ + zreturn: + if(x==0. && !flip) + return 0.; + return math·NaN(); + } + if(flip){ + if(y == .5) + return 1./math·sqrt(x); + y = -y; + }else if(y == .5) + return math·sqrt(x); + return math·exp(y * math·log(x)); + } + x = math·frexp(x, &ex); + ey = 0; + i = ye; + if(i) + for(;;){ + if(i & 1){ + xy *= x; + ey += ex; + } + i >>= 1; + if(i == 0) + break; + x *= x; + ex <<= 1; + if(x < .5){ + x += x; + ex -= 1; + } + } + if(flip){ + xy = 1. / xy; + ey = -ey; + } + return math·ldexp(xy, ey); +} diff --git a/src/base/math/pow10.c b/src/base/math/pow10.c new file mode 100644 index 0000000..a62428c --- /dev/null +++ b/src/base/math/pow10.c @@ -0,0 +1,48 @@ +#include +#include + +/* + * this table might overflow 127-bit exponent representations. + * in that case, truncate it after 1.0e38. + * it is important to get all one can from this + * routine since it is used in atof to scale numbers. + * the presumption is that C converts fp numbers better + * than multipication of lower powers of 10. + */ +static double tab[] = +{ + 1.0e0, 1.0e1, 1.0e2, 1.0e3, 1.0e4, 1.0e5, 1.0e6, 1.0e7, 1.0e8, 1.0e9, + 1.0e10, 1.0e11, 1.0e12, 1.0e13, 1.0e14, 1.0e15, 1.0e16, 1.0e17, 1.0e18, 1.0e19, + 1.0e20, 1.0e21, 1.0e22, 1.0e23, 1.0e24, 1.0e25, 1.0e26, 1.0e27, 1.0e28, 1.0e29, + 1.0e30, 1.0e31, 1.0e32, 1.0e33, 1.0e34, 1.0e35, 1.0e36, 1.0e37, 1.0e38, 1.0e39, + 1.0e40, 1.0e41, 1.0e42, 1.0e43, 1.0e44, 1.0e45, 1.0e46, 1.0e47, 1.0e48, 1.0e49, + 1.0e50, 1.0e51, 1.0e52, 1.0e53, 1.0e54, 1.0e55, 1.0e56, 1.0e57, 1.0e58, 1.0e59, + 1.0e60, 1.0e61, 1.0e62, 1.0e63, 1.0e64, 1.0e65, 1.0e66, 1.0e67, 1.0e68, 1.0e69, + 1.0e70, 1.0e71, 1.0e72, 1.0e73, 1.0e74, 1.0e75, 1.0e76, 1.0e77, 1.0e78, 1.0e79, + 1.0e80, 1.0e81, 1.0e82, 1.0e83, 1.0e84, 1.0e85, 1.0e86, 1.0e87, 1.0e88, 1.0e89, + 1.0e90, 1.0e91, 1.0e92, 1.0e93, 1.0e94, 1.0e95, 1.0e96, 1.0e97, 1.0e98, 1.0e99, + 1.0e100,1.0e101,1.0e102,1.0e103,1.0e104,1.0e105,1.0e106,1.0e107,1.0e108,1.0e109, + 1.0e110,1.0e111,1.0e112,1.0e113,1.0e114,1.0e115,1.0e116,1.0e117,1.0e118,1.0e119, + 1.0e120,1.0e121,1.0e122,1.0e123,1.0e124,1.0e125,1.0e126,1.0e127,1.0e128,1.0e129, + 1.0e130,1.0e131,1.0e132,1.0e133,1.0e134,1.0e135,1.0e136,1.0e137,1.0e138,1.0e139, + 1.0e140,1.0e141,1.0e142,1.0e143,1.0e144,1.0e145,1.0e146,1.0e147,1.0e148,1.0e149, + 1.0e150,1.0e151,1.0e152,1.0e153,1.0e154,1.0e155,1.0e156,1.0e157,1.0e158,1.0e159, +}; + +double +math·pow10(int n) +{ + int m; + + if(n < 0) { + n = -n; + if(n < sizeof(tab)/sizeof(tab[0])) + return 1/tab[n]; + m = n/2; + return 1/(math·pow10(m) * math·pow10(n-m)); + } + if(n < sizeof(tab)/sizeof(tab[0])) + return tab[n]; + m = n/2; + return math·pow10(m) * math·pow10(n-m); +} diff --git a/src/base/math/rules.mk b/src/base/math/rules.mk new file mode 100644 index 0000000..b66bfce --- /dev/null +++ b/src/base/math/rules.mk @@ -0,0 +1 @@ +SRCS_$(d)+=$(wildcard $(d)/math/*.c) diff --git a/src/base/math/sin.c b/src/base/math/sin.c new file mode 100644 index 0000000..486c5fa --- /dev/null +++ b/src/base/math/sin.c @@ -0,0 +1,67 @@ +/* + C program for floating point sin/cos. + Calls modf. + There are no error exits. + Coefficients are #3370 from Hart & Cheney (18.80D). +*/ + +#include +#include + +#define p0 .1357884097877375669092680e8 +#define p1 -.4942908100902844161158627e7 +#define p2 .4401030535375266501944918e6 +#define p3 -.1384727249982452873054457e5 +#define p4 .1459688406665768722226959e3 +#define q0 .8644558652922534429915149e7 +#define q1 .4081792252343299749395779e6 +#define q2 .9463096101538208180571257e4 +#define q3 .1326534908786136358911494e3 + +static double +sinus(double arg, int quad) +{ + double e, f, ysq, x, y, temp1, temp2; + int k; + + x = arg; + if(x < 0) { + x = -x; + quad += 2; + } + x *= 1/math·PIO2; /* underflow? */ + if(x > 32764) { + y = math·modf(x, &e); + e += quad; + math·modf(0.25*e, &f); + quad = e - 4*f; + } else { + k = x; + y = x - k; + quad += k; + quad &= 3; + } + if(quad & 1) + y = 1-y; + if(quad > 1) + y = -y; + + ysq = y*y; + temp1 = ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y; + temp2 = ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0); + return temp1/temp2; +} + +double +math·cos(double arg) +{ + if(arg < 0) + arg = -arg; + return sinus(arg, 1); +} + +double +math·sin(double arg) +{ + return sinus(arg, 0); +} diff --git a/src/base/math/sinh.c b/src/base/math/sinh.c new file mode 100644 index 0000000..ce036ae --- /dev/null +++ b/src/base/math/sinh.c @@ -0,0 +1,62 @@ +#include +#include + +/* + * sinh(arg) returns the hyperbolic sine of its floating- + * point argument. + * + * The exponential function is called for arguments + * 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 function for + * all arguments. + */ + +static double p0 = -0.6307673640497716991184787251e+6; +static double p1 = -0.8991272022039509355398013511e+5; +static double p2 = -0.2894211355989563807284660366e+4; +static double p3 = -0.2630563213397497062819489e+2; +static double q0 = -0.6307673640497716991212077277e+6; +static double q1 = 0.1521517378790019070696485176e+5; +static double q2 = -0.173678953558233699533450911e+3; + +double +math·sinh(double arg) +{ + double temp, argsq; + int sign; + + sign = 0; + if(arg < 0) { + arg = -arg; + sign++; + } + if(arg > 21) { + temp = math·exp(arg)/2; + goto out; + } + if(arg > 0.5) { + temp = (math·exp(arg) - math·exp(-arg))/2; + goto out; + } + argsq = arg*arg; + temp = (((p3*argsq+p2)*argsq+p1)*argsq+p0)*arg; + temp /= (((argsq+q2)*argsq+q1)*argsq+q0); +out: + if(sign) + temp = -temp; + return temp; +} + +double +math·cosh(double arg) +{ + if(arg < 0) + arg = - arg; + if(arg > 21) + return math·exp(arg)/2; + return (math·exp(arg) + math·exp(-arg))/2; +} diff --git a/src/base/math/sqrt.c b/src/base/math/sqrt.c new file mode 100644 index 0000000..a94d69f --- /dev/null +++ b/src/base/math/sqrt.c @@ -0,0 +1,54 @@ +/* + sqrt returns the square root of its floating + point argument. Newton's method. + + calls frexp +*/ + +#include +#include + +double +math·sqrt(double arg) +{ + double x, temp; + int exp, i; + + if(arg <= 0) { + if(arg < 0) + return math·NaN(); + return 0; + } + if(math·isInf(arg, 1)) + return arg; + x = math·frexp(arg, &exp); + while(x < 0.5) { + x *= 2; + exp--; + } + /* + * NOTE + * this wont work on 1's comp + */ + if(exp & 1) { + x *= 2; + exp--; + } + temp = 0.5 * (1.0+x); + + while(exp > 60) { + temp *= (1L<<30); + exp -= 60; + } + while(exp < -60) { + temp /= (1L<<30); + exp += 60; + } + if(exp >= 0) + temp *= 1L << (exp/2); + else + temp /= 1L << (-exp/2); + for(i=0; i<=4; i++) + temp = 0.5*(temp + arg/temp); + return temp; +} diff --git a/src/base/math/tan.c b/src/base/math/tan.c new file mode 100644 index 0000000..08b8806 --- /dev/null +++ b/src/base/math/tan.c @@ -0,0 +1,67 @@ +/* + floating point tangent + + A series is used after range reduction. + Coefficients are #4285 from Hart & Cheney. (19.74D) + */ + +#include +#include + +static double p0 = -0.1306820264754825668269611177e+5; +static double p1 = 0.1055970901714953193602353981e+4; +static double p2 = -0.1550685653483266376941705728e+2; +static double p3 = 0.3422554387241003435328470489e-1; +static double p4 = 0.3386638642677172096076369e-4; +static double q0 = -0.1663895238947119001851464661e+5; +static double q1 = 0.4765751362916483698926655581e+4; +static double q2 = -0.1555033164031709966900124574e+3; + +double +math·tan(double arg) +{ + double temp, e, x, xsq; + int flag, sign, i; + + flag = 0; + sign = 0; + if(arg < 0){ + arg = -arg; + sign++; + } + arg = 2*arg/math·PIO2; /* overflow? */ + x = math·modf(arg, &e); + i = e; + switch(i%4) { + case 1: + x = 1 - x; + flag = 1; + break; + + case 2: + sign = !sign; + flag = 1; + break; + + case 3: + x = 1 - x; + sign = !sign; + break; + + case 0: + break; + } + + xsq = x*x; + temp = ((((p4*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0)*x; + temp = temp/(((xsq+q2)*xsq+q1)*xsq+q0); + + if(flag) { + if(temp == 0) + return math·NaN(); + temp = 1/temp; + } + if(sign) + temp = -temp; + return temp; +} diff --git a/src/base/math/tanh.c b/src/base/math/tanh.c new file mode 100644 index 0000000..81e428a --- /dev/null +++ b/src/base/math/tanh.c @@ -0,0 +1,25 @@ +#include +#include + +/* + tanh(arg) computes the hyperbolic tangent of its floating + point argument. + + sinh and cosh are called except for large arguments, which + would cause overflow improperly. + */ + +double +math·tanh(double arg) +{ + + if(arg < 0) { + arg = -arg; + if(arg > 21) + return -1; + return -math·sinh(arg)/math·cosh(arg); + } + if(arg > 21) + return 1; + return math·sinh(arg)/math·cosh(arg); +} -- cgit v1.2.1