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. --- include/base.h | 1 + include/base/math.h | 39 ++++ include/base/string.h | 8 +- src/base/fmt/float.c | 507 +---------------------------------------- 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 ++ src/base/rng/exponential.c | 2 +- src/base/rng/normal.c | 4 +- src/base/rng/poisson.c | 21 +- src/base/rules.mk | 1 + src/base/string/raw/asfloat.c | 474 ++++++++++++++++++++++++++++++++++++++ src/base/string/raw/atof.c | 8 + src/base/string/raw/atoi.c | 4 +- src/libbio/newick.c | 2 +- src/libbio/phylo.c | 4 +- sys/linux/amd64/arch/types.h | 11 + sys/linux/arm/arch/types.h | 10 + sys/linux/arm64/arch/types.h | 12 + sys/linux/i386/arch/types.h | 11 + sys/linux/riscv64/arch/types.h | 12 + 37 files changed, 1452 insertions(+), 524 deletions(-) create mode 100644 include/base/math.h 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 create mode 100644 src/base/string/raw/asfloat.c create mode 100644 src/base/string/raw/atof.c diff --git a/include/base.h b/include/base.h index 801b06c..99cd40b 100644 --- a/include/base.h +++ b/include/base.h @@ -46,6 +46,7 @@ typedef wchar_t wchar; #include #include #include +#include #include // ----------------------------------------------------------------------------- diff --git a/include/base/math.h b/include/base/math.h new file mode 100644 index 0000000..69d344e --- /dev/null +++ b/include/base/math.h @@ -0,0 +1,39 @@ +#pragma once + +double math·NaN(void); +double math·Inf(int); +int math·isNaN(double); +int math·isInf(double, int); +ulong math·umuldiv(ulong, ulong, ulong); +long math·muldiv(long, long, long); + +double math·pow(double, double); +double math·atan2(double, double); +double math·fabs(double); +double math·atan(double); +int math·abs(int); +long math·labs(long); +double math·ldexp(double, int); +double math·log(double); +double math·log10(double); +double math·exp(double); +double math·floor(double); +double math·ceil(double); +double math·sin(double); +double math·cos(double); +double math·tan(double); +double math·asin(double); +double math·acos(double); +double math·sinh(double); +double math·cosh(double); +double math·tanh(double); +double math·sqrt(double); +double math·fmod(double, double); +double math·modf(double, double*); +double math·frexp(double, int*); + +double math·copysign(double, double); + +#define math·HUGE 3.4028234e38 +#define math·PIO2 1.570796326794896619231e0 +#define math·PI (math·PIO2+math·PIO2) diff --git a/include/base/string.h b/include/base/string.h index 238ebf9..a59c553 100644 --- a/include/base/string.h +++ b/include/base/string.h @@ -29,9 +29,15 @@ int str·compare(char *, char *); int str·ncompare(char *, intptr len, char *); int str·ecompare(char *, char *end, char *); -int str·atoi(char *s); +/* old school interfaces (no error checking) */ +long str·atoi(char *s); +double str·atof(char *s); char *str·itoa(char *s, int x); +/* nicer */ +vlong str·asint(char *s, char **end); +double str·asfloat(char *s, char **end); + /* augmented string functions */ string string·makecap(char *s, vlong len, vlong cap); string string·makelen(char *s, vlong len); diff --git a/src/base/fmt/float.c b/src/base/fmt/float.c index ba2637e..c99630c 100644 --- a/src/base/fmt/float.c +++ b/src/base/fmt/float.c @@ -157,507 +157,6 @@ sub1(char *a, int n) rt·exit(1); } -// ----------------------------------------------------------------------- -// strtod - -#define Nbits 28 -#define Nmant 53 -#define Prec ((Nmant+Nbits+1)/Nbits) - -#define Sigbit (1<<(Prec*Nbits-Nmant)) /* first significant bit of Prec-th word */ -#define Ndig 1500 -#define One (ulong)(1<>1) -#define Maxe 310 - -#define Fsign (1<<0) /* found - */ -#define Fesign (1<<1) /* found e- */ -#define Fdpoint (1<<2) /* found . */ - -#define S0 0 /* _ _S0 +S1 #S2 .S3 */ -#define S1 1 /* _+ #S2 .S3 */ -#define S2 2 /* _+# #S2 .S4 eS5 */ -#define S3 3 /* _+. #S4 */ -#define S4 4 /* _+#.# #S4 eS5 */ -#define S5 5 /* _+#.#e +S6 #S7 */ -#define S6 6 /* _+#.#e+ #S7 */ -#define S7 7 /* _+#.#e+# #S7 */ - -typedef struct Tab Tab; -struct Tab -{ - int bp; - int siz; - char *cmp; -}; - -static ulong -umuldiv(ulong a, ulong b, ulong c) -{ - double d; - - d = ((double)a * (double)b) / (double)c; - if(d >= 4294967295.) - d = 4294967295.; - return (ulong)d; -} - -static void -frnorm(ulong *f) -{ - int i, c; - - c = 0; - for(i=Prec-1; i>0; i--) { - f[i] += c; - c = f[i] >> Nbits; - f[i] &= One-1; - } - f[0] += c; -} - -static int -fpcmp(char *a, ulong* f) -{ - ulong tf[Prec]; - int i, d, c; - - for(i=0; i> Nbits) + '0'; - tf[0] &= One-1; - - /* compare next digit */ - c = *a; - if(c == 0) { - if('0' < d) - return -1; - if(tf[0] != 0) - goto cont; - for(i=1; i d) - return +1; - if(c < d) - return -1; - a++; - cont:; -} -} - -static void -divby(char *a, int *na, int b) -{ - int n, c; - char *p; - - p = a; - n = 0; - while(n>>b == 0){ - c = *a++; - if(c == 0) { - while(n) { - c = n*10; - if(c>>b) - break; - n = c; - } - goto xx; - } - n = n*10 + c-'0'; - (*na)--; - } - for(;;){ - c = n>>b; - n -= c<>b; - n -= c<= (int)(arrlen(tab1))) - d = (int)(arrlen(tab1))-1; - t = tab1 + d; - b = t->bp; - if(mem·compare(a, t->siz, t->cmp) > 0) - d--; - *dp -= d; - *bp += b; - divby(a, na, b); -} - -static void -mulby(char *a, char *p, char *q, int b) -{ - int n, c; - - n = 0; - *p = 0; - for(;;) { - q--; - if(q < a) - break; - c = *q - '0'; - c = (c<= (int)(arrlen(tab2))) - d = (int)(arrlen(tab2))-1; - t = tab2 + d; - b = t->bp; - if(mem·compare(a, t->siz, t->cmp) < 0) - d--; - p = a + *na; - *bp -= b; - *dp += d; - *na += d; - mulby(a, p+d, p, b); -} - -static int -cmp(char *a, char *b) -{ - int c1, c2; - - while((c1 = *b++) != '\0') { - c2 = *a++; - if(utf8·isupper(c2)) - c2 = utf8·tolower(c2); - if(c1 != c2) - return 1; - } - return 0; -} - -double -fmtstrtod(char *as, char **aas) -{ - int na, ex, dp, bp, c, i, flag, state; - ulong low[Prec], hig[Prec], mid[Prec]; - double d; - char *s, a[Ndig]; - - flag = 0; /* Fsign, Fesign, Fdpoint */ - na = 0; /* number of digits of a[] */ - dp = 0; /* na of decimal point */ - ex = 0; /* exonent */ - - state = S0; - for(s=as;;s++){ - c = *s; - if('0' <= c && c <= '9'){ - switch(state){ - case S0: case S1: case S2: - state = S2; - break; - case S3: case S4: - state = S4; - break; - case S5: case S6: case S7: - state = S7; - ex = ex*10 + (c-'0'); - continue; - } - - if(na == 0 && c == '0'){ - dp--; - continue; - } - if(na < Ndig-50) - a[na++] = c; - continue; - } - switch(c){ - case '\t': case '\n': case '\v': case '\f': case '\r': case ' ': - if(state == S0) - continue; - break; - case '-': - if(state == S0) - flag |= Fsign; - else - flag |= Fesign; - case '+': - if(state == S0) - state = S1; - else - if(state == S5) - state = S6; - else - break; /* syntax */ - continue; - case '.': - flag |= Fdpoint; - dp = na; - if(state == S0 || state == S1){ - state = S3; - continue; - } - if(state == S2){ - state = S4; - continue; - } - break; - case 'e': case 'E': - if(state == S2 || state == S4){ - state = S5; - continue; - } - break; - } - break; - } - - /* clean up return char-pointer */ - switch(state) { - case S0: - if(cmp(s, "nan") == 0){ - if(aas != nil) - *aas = s+3; - goto retnan; - } - case S1: - if(cmp(s, "infinity") == 0){ - if(aas != nil) - *aas = s+8; - goto retinf; - } - if(cmp(s, "inf") == 0){ - if(aas != nil) - *aas = s+3; - goto retinf; - } - case S3: - if(aas != nil) - *aas = as; - goto ret0; /* no digits found */ - case S6: - s--; /* back over +- */ - case S5: - s--; /* back over e */ - break; - } - if(aas != nil) - *aas = s; - - if(flag & Fdpoint) - while(na > 0 && a[na-1] == '0') - na--; - if(na == 0) - goto ret0; /* zero */ - a[na] = 0; - if(!(flag & Fdpoint)) - dp = na; - if(flag & Fesign) - ex = -ex; - dp += ex; - if(dp < -Maxe){ - /* errno = ERANGE; */ - goto ret0; /* underflow by exp */ - } else - if(dp > +Maxe) - goto retinf; /* overflow by exp */ - - /* - * normalize the decimal ascii number - * to range .[5-9][0-9]* e0 - */ - bp = 0; /* binary exponent */ - while(dp > 0) - divascii(a, &na, &dp, &bp); - while(dp < 0 || a[0] < '5') - mulascii(a, &na, &dp, &bp); - - /* close approx by naive conversion */ - mid[0] = 0; - mid[1] = 1; - for(i=0; (c=a[i]) != '\0'; i++) { - mid[0] = mid[0]*10 + (c-'0'); - mid[1] = mid[1]*10; - if(i >= 8) - break; - } - low[0] = umuldiv(mid[0], One, mid[1]); - hig[0] = umuldiv(mid[0]+1, One, mid[1]); - for(i=1; i>= 1; - } - frnorm(mid); - - /* compare */ - c = fpcmp(a, mid); - if(c > 0) { - c = 1; - for(i=0; i= Sigbit/2) { - mid[Prec-1] += Sigbit; - frnorm(mid); - } - goto out; - -ret0: - return 0; - -retnan: - return NaN(); - -retinf: - /* Unix strtod requires these. Plan 9 would return Inf(0) or Inf(-1). */ - /* errno = ERANGE; */ - if(flag & Fsign) - return -HUGE_VAL; - return HUGE_VAL; - -out: - d = 0; - for(i=0; i g) { if(add1(s, NSIGNIF)){ /* gained a digit */ @@ -763,7 +262,7 @@ dtoa(double f, char *s, int *exp, int *neg, int *len) c = s[i]; if(c != '0'){ s[i] = '0'; - g=fmtstrtod(s, nil); + g=str·asfloat(s, nil); if(g != f){ s[i] = c; break; 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); +} diff --git a/src/base/rng/exponential.c b/src/base/rng/exponential.c index c07e007..9aa0d0c 100644 --- a/src/base/rng/exponential.c +++ b/src/base/rng/exponential.c @@ -7,5 +7,5 @@ rng·exponential(double lambda) double f; f = rng·random(); - return -log(1 - f)/lambda; + return -math·log(1 - f)/lambda; } diff --git a/src/base/rng/normal.c b/src/base/rng/normal.c index aab5731..8e3e5d4 100644 --- a/src/base/rng/normal.c +++ b/src/base/rng/normal.c @@ -53,7 +53,7 @@ erfinv(double x) z2 = ((((((b7*r+b6)*r+b5)*r+b4)*r+b3)*r+b2)*r+b1)*r + b0; return s*(x*z1) / z2; } - r = sqrt(Ln2 - log(1.0-x)); + r = math·sqrt(Ln2 - math·log(1.0-x)); if(r <= 5.0) { r -= 1.6; z1 = ((((((c7*r+c6)*r+c5)*r+c4)*r+c3)*r+c2)*r+c1)*r + c0; @@ -73,5 +73,5 @@ rng·normal(void) double f; f = rng·random(); - return sqrt(2)*erfinv(2*f-1); + return math·sqrt(2)*erfinv(2*f-1); } diff --git a/src/base/rng/poisson.c b/src/base/rng/poisson.c index 3ec15c9..776070c 100644 --- a/src/base/rng/poisson.c +++ b/src/base/rng/poisson.c @@ -28,7 +28,7 @@ log1pmx(double x, double off) return x*x*r; } - return log(1+x) - off; + return math·log(1+x) - off; } static inline @@ -49,14 +49,14 @@ procf(double mu, double s, int64 K, double *px, double *py, double *fx, double * if(K < 10) { *px = -mu; - *py = pow(mu,K) / factorial[K]; + *py = math·pow(mu,K) / factorial[K]; }else{ d = 0.08333333333333333/K; d = d - 4.8*d*d*d; V = (mu - K) / K; *px = K*log1pmx(V,mu-K) - d; - *py = 0.3989422804014327/sqrt(K); + *py = 0.3989422804014327/math·sqrt(K); } X = (K - mu + 0.5)/s; @@ -66,21 +66,20 @@ procf(double mu, double s, int64 K, double *px, double *py, double *fx, double * return c; } -static inline -uint64 +static inline uint64 bigpoisson(double mu) { int64 L,K; double G,s,d,U,E,T; double px,py,fx,fy,c; - s = sqrt(mu); + s = math·sqrt(mu); d = 6*mu*mu; - L = floor(mu - 1.1484); + L = math·floor(mu - 1.1484); stepN: G = mu + s*rng·normal(); - K = floor(G); + K = math·floor(G); if(K<0) goto stepP; stepI: @@ -99,13 +98,13 @@ stepE: E = rng·exponential(1.0); U = rng·random(); U = U + U - 1; - T = 1.8 + copysign(E,U); + T = 1.8 + math·copysign(E,U); if(T < 0.6744) goto stepE; - K = floor(mu + s*T); + K = math·floor(mu + s*T); c = procf(mu, s, K, &px, &py, &fx, &fy); stepH: - if(c*fabs(U) > (py*exp(px + E) - fy*exp(fx + E))) + if(c*math·fabs(U) > (py*math·exp(px + E) - fy*math·exp(fx + E))) goto stepE; return K; } diff --git a/src/base/rules.mk b/src/base/rules.mk index fe7bb13..8199f25 100644 --- a/src/base/rules.mk +++ b/src/base/rules.mk @@ -15,6 +15,7 @@ include $(d)/fs/rules.mk include $(d)/gz/rules.mk include $(d)/mem/rules.mk include $(d)/mmap/rules.mk +include $(d)/math/rules.mk include $(d)/rng/rules.mk include $(d)/sort/rules.mk include $(d)/string/rules.mk diff --git a/src/base/string/raw/asfloat.c b/src/base/string/raw/asfloat.c new file mode 100644 index 0000000..7d35626 --- /dev/null +++ b/src/base/string/raw/asfloat.c @@ -0,0 +1,474 @@ +#include +#include + +#define Nbits 28 +#define Nmant 53 +#define Prec ((Nmant+Nbits+1)/Nbits) + +#define Sigbit (1<<(Prec*Nbits-Nmant)) /* first significant bit of Prec-th word */ +#define Ndig 1500 +#define One (ulong)(1<>1) +#define Maxe 310 + +#define Fsign (1<<0) /* found - */ +#define Fesign (1<<1) /* found e- */ +#define Fdpoint (1<<2) /* found . */ + +#define S0 0 /* _ _S0 +S1 #S2 .S3 */ +#define S1 1 /* _+ #S2 .S3 */ +#define S2 2 /* _+# #S2 .S4 eS5 */ +#define S3 3 /* _+. #S4 */ +#define S4 4 /* _+#.# #S4 eS5 */ +#define S5 5 /* _+#.#e +S6 #S7 */ +#define S6 6 /* _+#.#e+ #S7 */ +#define S7 7 /* _+#.#e+# #S7 */ + +typedef struct Tab Tab; +struct Tab +{ + int bp; + int siz; + char *cmp; +}; + +static ulong +umuldiv(ulong a, ulong b, ulong c) +{ + double d; + + d = ((double)a * (double)b) / (double)c; + if(d >= 4294967295.) + d = 4294967295.; + return (ulong)d; +} + +static void +frnorm(ulong *f) +{ + int i, c; + + c = 0; + for(i=Prec-1; i>0; i--) { + f[i] += c; + c = f[i] >> Nbits; + f[i] &= One-1; + } + f[0] += c; +} + +static int +fpcmp(char *a, ulong* f) +{ + ulong tf[Prec]; + int i, d, c; + + for(i=0; i> Nbits) + '0'; + tf[0] &= One-1; + + /* compare next digit */ + c = *a; + if(c == 0) { + if('0' < d) + return -1; + if(tf[0] != 0) + goto cont; + for(i=1; i d) + return +1; + if(c < d) + return -1; + a++; + cont:; +} +} + +static void +divby(char *a, int *na, int b) +{ + int n, c; + char *p; + + p = a; + n = 0; + while(n>>b == 0){ + c = *a++; + if(c == 0) { + while(n) { + c = n*10; + if(c>>b) + break; + n = c; + } + goto xx; + } + n = n*10 + c-'0'; + (*na)--; + } + for(;;){ + c = n>>b; + n -= c<>b; + n -= c<= (int)(arrlen(tab1))) + d = (int)(arrlen(tab1))-1; + t = tab1 + d; + b = t->bp; + if(mem·compare(a, t->siz, t->cmp) > 0) + d--; + *dp -= d; + *bp += b; + divby(a, na, b); +} + +static void +mulby(char *a, char *p, char *q, int b) +{ + int n, c; + + n = 0; + *p = 0; + for(;;) { + q--; + if(q < a) + break; + c = *q - '0'; + c = (c<= (int)(arrlen(tab2))) + d = (int)(arrlen(tab2))-1; + t = tab2 + d; + b = t->bp; + if(mem·compare(a, t->siz, t->cmp) < 0) + d--; + p = a + *na; + *bp -= b; + *dp += d; + *na += d; + mulby(a, p+d, p, b); +} + +static int +cmp(char *a, char *b) +{ + int c1, c2; + + while((c1 = *b++) != '\0') { + c2 = *a++; + if(utf8·isupper(c2)) + c2 = utf8·tolower(c2); + if(c1 != c2) + return 1; + } + return 0; +} + +double +str·asfloat(char *as, char **aas) +{ + int na, ex, dp, bp, c, i, flag, state; + ulong low[Prec], hig[Prec], mid[Prec]; + double d; + char *s, a[Ndig]; + + flag = 0; /* Fsign, Fesign, Fdpoint */ + na = 0; /* number of digits of a[] */ + dp = 0; /* na of decimal point */ + ex = 0; /* exonent */ + + state = S0; + for(s=as;;s++){ + c = *s; + if('0' <= c && c <= '9'){ + switch(state){ + case S0: case S1: case S2: + state = S2; + break; + case S3: case S4: + state = S4; + break; + case S5: case S6: case S7: + state = S7; + ex = ex*10 + (c-'0'); + continue; + } + + if(na == 0 && c == '0'){ + dp--; + continue; + } + if(na < Ndig-50) + a[na++] = c; + continue; + } + switch(c){ + case '\t': case '\n': case '\v': case '\f': case '\r': case ' ': + if(state == S0) + continue; + break; + case '-': + if(state == S0) + flag |= Fsign; + else + flag |= Fesign; + case '+': + if(state == S0) + state = S1; + else + if(state == S5) + state = S6; + else + break; /* syntax */ + continue; + case '.': + flag |= Fdpoint; + dp = na; + if(state == S0 || state == S1){ + state = S3; + continue; + } + if(state == S2){ + state = S4; + continue; + } + break; + case 'e': case 'E': + if(state == S2 || state == S4){ + state = S5; + continue; + } + break; + } + break; + } + + /* clean up return char-pointer */ + switch(state) { + case S0: + if(cmp(s, "nan") == 0){ + if(aas != nil) + *aas = s+3; + goto retnan; + } + case S1: + if(cmp(s, "infinity") == 0){ + if(aas != nil) + *aas = s+8; + goto retinf; + } + if(cmp(s, "inf") == 0){ + if(aas != nil) + *aas = s+3; + goto retinf; + } + case S3: + if(aas != nil) + *aas = as; + goto ret0; /* no digits found */ + case S6: + s--; /* back over +- */ + case S5: + s--; /* back over e */ + break; + } + if(aas != nil) + *aas = s; + + if(flag & Fdpoint) + while(na > 0 && a[na-1] == '0') + na--; + if(na == 0) + goto ret0; /* zero */ + a[na] = 0; + if(!(flag & Fdpoint)) + dp = na; + if(flag & Fesign) + ex = -ex; + dp += ex; + if(dp < -Maxe){ + /* errno = ERANGE; */ + goto ret0; /* underflow by exp */ + } else + if(dp > +Maxe) + goto retinf; /* overflow by exp */ + + /* + * normalize the decimal ascii number + * to range .[5-9][0-9]* e0 + */ + bp = 0; /* binary exponent */ + while(dp > 0) + divascii(a, &na, &dp, &bp); + while(dp < 0 || a[0] < '5') + mulascii(a, &na, &dp, &bp); + + /* close approx by naive conversion */ + mid[0] = 0; + mid[1] = 1; + for(i=0; (c=a[i]) != '\0'; i++) { + mid[0] = mid[0]*10 + (c-'0'); + mid[1] = mid[1]*10; + if(i >= 8) + break; + } + low[0] = umuldiv(mid[0], One, mid[1]); + hig[0] = umuldiv(mid[0]+1, One, mid[1]); + for(i=1; i>= 1; + } + frnorm(mid); + + /* compare */ + c = fpcmp(a, mid); + if(c > 0) { + c = 1; + for(i=0; i= Sigbit/2) { + mid[Prec-1] += Sigbit; + frnorm(mid); + } + goto out; + +ret0: + return 0; + +retnan: + return math·NaN(); + +retinf: + if(flag & Fsign) + return math·Inf(-1); + return math·Inf(+1); + +out: + d = 0; + for(i=0; i +#include + +double +str·atof(char *s) +{ + return str·asfloat(s, nil); +} diff --git a/src/base/string/raw/atoi.c b/src/base/string/raw/atoi.c index 8084e3e..e62ec62 100644 --- a/src/base/string/raw/atoi.c +++ b/src/base/string/raw/atoi.c @@ -1,7 +1,7 @@ -int +long str·atoi(char *s) { - int n = 0; + long n = 0; while(*s) n = 10*n + (*s++ - '0'); diff --git a/src/libbio/newick.c b/src/libbio/newick.c index a555379..6071c1e 100644 --- a/src/libbio/newick.c +++ b/src/libbio/newick.c @@ -118,7 +118,7 @@ lex(io·Peeker s, void* fp) *c = 0; tok.kind = tok·number; - tok.lit.x = atof(b); + tok.lit.x = str·atof(b); return tok; case '\"': diff --git a/src/libbio/phylo.c b/src/libbio/phylo.c index 1b0ae59..2343e23 100644 --- a/src/libbio/phylo.c +++ b/src/libbio/phylo.c @@ -393,10 +393,10 @@ phylo·reroot(bio·Tree *tree, bio·Node *node, double d) // TODO: should check that node is part of this tree? // TODO: should we check if node->parent != nil? - if(fabs(d) < PREC){ + if(math·fabs(d) < PREC){ new = node; rotateparent(node->parent, node); - }else if(fabs(d-node->dist) < PREC){ + }else if(math·fabs(d-node->dist) < PREC){ new = node->parent; if (new->parent->parent) { rotateparent(new->parent->parent, new->parent); diff --git a/sys/linux/amd64/arch/types.h b/sys/linux/amd64/arch/types.h index d22169c..b6bcd3d 100644 --- a/sys/linux/amd64/arch/types.h +++ b/sys/linux/amd64/arch/types.h @@ -1,5 +1,16 @@ #pragma once +typedef union sys·DoubleWord sys·DoubleWord; + +/* little-endian */ +union sys·DoubleWord +{ + double x; + struct{ + uint32 lo; + uint32 hi; + }; +}; #if 0 /* * copied from musl: diff --git a/sys/linux/arm/arch/types.h b/sys/linux/arm/arch/types.h index 536ca5f..8a59abb 100644 --- a/sys/linux/arm/arch/types.h +++ b/sys/linux/arm/arch/types.h @@ -1,5 +1,15 @@ #pragma once +typedef union sys·DoubleWord sys·DoubleWord; + +union sys·DoubleWord +{ + double x; + struct{ + uint32 lo; + uint32 hi; + }; +}; #if 0 /* * copied from musl: diff --git a/sys/linux/arm64/arch/types.h b/sys/linux/arm64/arch/types.h index 65d05f6..3cb2277 100644 --- a/sys/linux/arm64/arch/types.h +++ b/sys/linux/arm64/arch/types.h @@ -1,5 +1,17 @@ #pragma once +typedef union sys·DoubleWord sys·DoubleWord; + +/* little-endian */ +union sys·DoubleWord +{ + double x; + struct{ + uint32 lo; + uint32 hi; + }; +}; + #if 0 struct sys·Info { diff --git a/sys/linux/i386/arch/types.h b/sys/linux/i386/arch/types.h index be2421a..0e39df6 100644 --- a/sys/linux/i386/arch/types.h +++ b/sys/linux/i386/arch/types.h @@ -1,5 +1,16 @@ #pragma once +typedef union sys·DoubleWord sys·DoubleWord; + +union sys·DoubleWord +{ + double x; + struct{ + uint32 lo; + uint32 hi; + }; +}; + #if 0 /* * copied from musl: diff --git a/sys/linux/riscv64/arch/types.h b/sys/linux/riscv64/arch/types.h index af908d2..5109adf 100644 --- a/sys/linux/riscv64/arch/types.h +++ b/sys/linux/riscv64/arch/types.h @@ -1,5 +1,17 @@ #pragma once +typedef union sys·DoubleWord sys·DoubleWord; + +/* little-endian */ +union sys·DoubleWord +{ + double x; + struct{ + uint32 lo; + uint32 hi; + }; +}; + #if 0 struct sys·Info { -- cgit v1.2.1