aboutsummaryrefslogtreecommitdiff
path: root/src/base/math
diff options
context:
space:
mode:
Diffstat (limited to 'src/base/math')
-rw-r--r--src/base/math/abs.c18
-rw-r--r--src/base/math/asin.c40
-rw-r--r--src/base/math/atan.c76
-rw-r--r--src/base/math/atan2.c24
-rw-r--r--src/base/math/copysign.c12
-rw-r--r--src/base/math/exp.c40
-rw-r--r--src/base/math/fabs.c10
-rw-r--r--src/base/math/floor.c27
-rw-r--r--src/base/math/fmod.c31
-rw-r--r--src/base/math/frexp.c120
-rw-r--r--src/base/math/nan.c54
-rw-r--r--src/base/math/pow.c69
-rw-r--r--src/base/math/pow10.c48
-rw-r--r--src/base/math/rules.mk1
-rw-r--r--src/base/math/sin.c67
-rw-r--r--src/base/math/sinh.c62
-rw-r--r--src/base/math/sqrt.c54
-rw-r--r--src/base/math/tan.c67
-rw-r--r--src/base/math/tanh.c25
19 files changed, 845 insertions, 0 deletions
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 <u.h>
+#include <base.h>
+
+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 <u.h>
+#include <base.h>
+
+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 <u.h>
+#include <base.h>
+
+#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 <u.h>
+#include <base.h>
+/*
+ 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 <u.h>
+#include <base.h>
+
+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 <u.h>
+#include <base.h>
+
+#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 <u.h>
+#include <base.h>
+
+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 <u.h>
+#include <base.h>
+/*
+ * 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 <u.h>
+#include <base.h>
+
+/*
+ * 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 <u.h>
+#include <base.h>
+
+/*
+ * 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<<SIG) * d;
+ *ep = -SIG;
+ }
+ *ep += ((x.hi >> 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<<SHIFT)-1);
+ x.hi &= ~((1<<SHIFT)-1);
+ x.hi |= (1<<(SHIFT-1)) | (z>>1);
+ }
+ while(deltae > 0){ /* shift bits down */
+ bits = deltae;
+ if(bits > SHIFT)
+ bits = SHIFT;
+ x.lo >>= bits;
+ x.lo |= (x.hi&((1<<bits)-1)) << (32-bits);
+ z = x.hi & ((1<<SHIFT)-1);
+ x.hi &= ~((1<<SHIFT)-1);
+ x.hi |= z>>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 <u.h>
+#include <base.h>
+
+#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 <u.h>
+#include <base.h>
+
+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 <u.h>
+#include <base.h>
+
+/*
+ * 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 <u.h>
+#include <base.h>
+
+#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 <u.h>
+#include <base.h>
+
+/*
+ * 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 <u.h>
+#include <base.h>
+
+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 <u.h>
+#include <base.h>
+
+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 <u.h>
+#include <base.h>
+
+/*
+ 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);
+}
generated by cgit v1.2.1 (git 2.18.0) at 2024-07-06 00:10:45 +0000