Feat: added math library

Used Plan9's libc as starting point. This cleans up dangling references
due to loss of libc.

37 files changed, 1452 insertions, 524 deletions

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 <base/flate.h>
 #include <base/gz.h>
 #include <base/sort.h>
+#include <base/math.h>
 #include <base/rng.h>
 
 // -----------------------------------------------------------------------------
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<<Nbits)
-#define Half (ulong)(One>>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<Prec; i++)
-        tf[i] = f[i];
-
-    for(;;) {
-        /* tf *= 10 */
-        for(i=0; i<Prec; i++)
-            tf[i] = tf[i]*10;
-        frnorm(tf);
-        d = (tf[0] >> 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<Prec; i++)
-                if(tf[i] != 0)
-                    goto cont;
-            return 0;
-        }
-        if(c > 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;
-        *p++ = c + '0';
-        c = *a++;
-        if(c == 0)
-            break;
-        n = n*10 + c-'0';
-    }
-    (*na)++;
-    xx:
-    while(n){
-        n = n*10;
-        c = n>>b;
-        n -= c<<b;
-        *p++ = c + '0';
-        (*na)++;
-    }
-    *p = 0;
-}
-
-static Tab tab1[] =
-{
-     1,  0, "",
-     3,  1, "7",
-     6,  2, "63",
-     9,  3, "511",
-    13,  4, "8191",
-    16,  5, "65535",
-    19,  6, "524287",
-    23,  7, "8388607",
-    26,  8, "67108863",
-    27,  9, "134217727",
-};
-
-static void
-divascii(char *a, int *na, int *dp, int *bp)
-{
-    int b, d;
-    Tab *t;
-
-    d = *dp;
-    if(d >= (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<<b) + n;
-        n = c/10;
-        c -= n*10;
-        p--;
-        *p = c + '0';
-    }
-    while(n) {
-        c = n;
-        n = c/10;
-        c -= n*10;
-        p--;
-        *p = c + '0';
-    }
-}
-
-static Tab tab2[] =
-{
-     1,  1, "",				/* dp = 0-0 */
-     3,  3, "125",
-     6,  5, "15625",
-     9,  7, "1953125",
-    13, 10, "1220703125",
-    16, 12, "152587890625",
-    19, 14, "19073486328125",
-    23, 17, "11920928955078125",
-    26, 19, "1490116119384765625",
-    27, 19, "7450580596923828125",		/* dp 8-9 */
-};
-
-static void
-mulascii(char *a, int *na, int *dp, int *bp)
-{
-    char *p;
-    int d, b;
-    Tab *t;
-
-    d = -*dp;
-    if(d >= (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<Prec; i++) {
-        low[i] = 0;
-        hig[i] = One-1;
-    }
-
-    /* binary search for closest mantissa */
-    for(;;) {
-        /* mid = (hig + low) / 2 */
-        c = 0;
-        for(i=0; i<Prec; i++) {
-            mid[i] = hig[i] + low[i];
-            if(c)
-                mid[i] += One;
-            c = mid[i] & 1;
-            mid[i] >>= 1;
-        }
-        frnorm(mid);
-
-        /* compare */
-        c = fpcmp(a, mid);
-        if(c > 0) {
-            c = 1;
-            for(i=0; i<Prec; i++)
-                if(low[i] != mid[i]) {
-                    c = 0;
-                    low[i] = mid[i];
-                }
-            if(c)
-                break;    /* between mid and hig */
-            continue;
-        }
-        if(c < 0) {
-            for(i=0; i<Prec; i++)
-                hig[i] = mid[i];
-            continue;
-        }
-
-        /* only hard part is if even/odd roundings wants to go up */
-        c = mid[Prec-1] & (Sigbit-1);
-        if(c == Sigbit/2 && (mid[Prec-1]&Sigbit) == 0)
-            mid[Prec-1] -= c;
-        break;    /* exactly mid */
-    }
-
-    /* normal rounding applies */
-    c = mid[Prec-1] & (Sigbit-1);
-    mid[Prec-1] -= c;
-    if(c >= 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<Prec; i++)
-        d = d*One + mid[i];
-    if(flag & Fsign)
-        d = -d;
-    d = ldexp(d, bp - Prec*Nbits);
-    if(d == 0)  /* underflow */
-        /* errno = ERANGE; */
-
-    return d;
-}
-
-#undef Nbits
-#undef Nmant
-#undef Prec
-
-#undef Sigbit
-#undef Ndig
-#undef One
-#undef Half
-#undef Maxe
-
-#undef Fsign
-#undef Fesign
-#undef Fdpoint
-
-#undef S0
-#undef S1
-#undef S2
-#undef S3
-#undef S4
-#undef S5
-#undef S6
-#undef S7
-
 static void
 fmtexp(char *p, int e, int ucase)
 {
@@ -712,7 +211,7 @@ dtoa(double f, char *s, int *exp, int *neg, int *len)
         return;
     }
 
-    frexp(f, &e2);
+    math·frexp(f, &e2);
     e = (int)(e2 * .301029995664);
     g = f * fpow10(-e);
     while(g < 1) {
@@ -736,7 +235,7 @@ dtoa(double f, char *s, int *exp, int *neg, int *len)
     fmtexp(s+NSIGNIF, e, 0);
 
     for(i=0; i<10; i++) {
-        g=fmtstrtod(s, nil);
+        g=str·asfloat(s, nil);
         if(f > 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 <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);
+}
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 <u.h>
+#include <base.h>
+
+#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<<Nbits)
+#define Half (ulong)(One>>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<Prec; i++)
+        tf[i] = f[i];
+
+    for(;;) {
+        /* tf *= 10 */
+        for(i=0; i<Prec; i++)
+            tf[i] = tf[i]*10;
+        frnorm(tf);
+        d = (tf[0] >> 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<Prec; i++)
+                if(tf[i] != 0)
+                    goto cont;
+            return 0;
+        }
+        if(c > 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;
+        *p++ = c + '0';
+        c = *a++;
+        if(c == 0)
+            break;
+        n = n*10 + c-'0';
+    }
+    (*na)++;
+    xx:
+    while(n){
+        n = n*10;
+        c = n>>b;
+        n -= c<<b;
+        *p++ = c + '0';
+        (*na)++;
+    }
+    *p = 0;
+}
+
+static Tab tab1[] =
+{
+     1,  0, "",
+     3,  1, "7",
+     6,  2, "63",
+     9,  3, "511",
+    13,  4, "8191",
+    16,  5, "65535",
+    19,  6, "524287",
+    23,  7, "8388607",
+    26,  8, "67108863",
+    27,  9, "134217727",
+};
+
+static void
+divascii(char *a, int *na, int *dp, int *bp)
+{
+    int b, d;
+    Tab *t;
+
+    d = *dp;
+    if(d >= (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<<b) + n;
+        n = c/10;
+        c -= n*10;
+        p--;
+        *p = c + '0';
+    }
+    while(n) {
+        c = n;
+        n = c/10;
+        c -= n*10;
+        p--;
+        *p = c + '0';
+    }
+}
+
+static Tab tab2[] =
+{
+     1,  1, "",				/* dp = 0-0 */
+     3,  3, "125",
+     6,  5, "15625",
+     9,  7, "1953125",
+    13, 10, "1220703125",
+    16, 12, "152587890625",
+    19, 14, "19073486328125",
+    23, 17, "11920928955078125",
+    26, 19, "1490116119384765625",
+    27, 19, "7450580596923828125",		/* dp 8-9 */
+};
+
+static void
+mulascii(char *a, int *na, int *dp, int *bp)
+{
+    char *p;
+    int d, b;
+    Tab *t;
+
+    d = -*dp;
+    if(d >= (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<Prec; i++) {
+        low[i] = 0;
+        hig[i] = One-1;
+    }
+
+    /* binary search for closest mantissa */
+    for(;;){
+        /* mid = (hig + low) / 2 */
+        c = 0;
+        for(i=0; i<Prec; i++) {
+            mid[i] = hig[i] + low[i];
+            if(c)
+                mid[i] += One;
+            c = mid[i] & 1;
+            mid[i] >>= 1;
+        }
+        frnorm(mid);
+
+        /* compare */
+        c = fpcmp(a, mid);
+        if(c > 0) {
+            c = 1;
+            for(i=0; i<Prec; i++)
+                if(low[i] != mid[i]) {
+                    c = 0;
+                    low[i] = mid[i];
+                }
+            if(c)
+                break;    /* between mid and hig */
+            continue;
+        }
+        if(c < 0) {
+            for(i=0; i<Prec; i++)
+                hig[i] = mid[i];
+            continue;
+        }
+
+        /* only hard part is if even/odd roundings wants to go up */
+        c = mid[Prec-1] & (Sigbit-1);
+        if(c == Sigbit/2 && (mid[Prec-1]&Sigbit) == 0)
+            mid[Prec-1] -= c;
+        break;    /* exactly mid */
+    }
+
+    /* normal rounding applies */
+    c = mid[Prec-1] & (Sigbit-1);
+    mid[Prec-1] -= c;
+    if(c >= 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<Prec; i++)
+        d = d*One + mid[i];
+    if(flag & Fsign)
+        d = -d;
+    d = math·ldexp(d, bp - Prec*Nbits);
+    /* if(d == 0)  /1* underflow *1/ */
+        /* errno = ERANGE; */
+    return d;
+}
diff --git a/src/base/string/raw/atof.c b/src/base/string/raw/atof.c
new file mode 100644
index 0000000..1383985
--- /dev/null
+++ b/src/base/string/raw/atof.c
@@ -0,0 +1,8 @@
+#include <u.h>
+#include <base.h>
+
+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
 {