#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;
}