ref: 12e226ac6798412271020d42645f69138715c519
dir: /utils/cc/mpatof.c/
#include "cc.h" enum { Mpscale = 29, /* safely smaller than bits in a long */ Mpprec = 36, /* Mpscale*Mpprec sb > largest fp exp */ Mpbase = 1L<<Mpscale, }; typedef struct { long a[Mpprec]; char ovf; } Mp; int mpatof(char*, double*); int mpatov(char *s, vlong *v); void mpint(Mp*, int); void mppow(Mp*, int, int); void mpmul(Mp*, int); void mpadd(Mp*, Mp*); int mptof(Mp*, double*); /* * convert a string, s, to floating in *d * return conversion overflow. * required syntax is [+-]d*[.]d*[e[+-]d*] */ int mpatof(char *s, double *d) { Mp a, b; int dp, c, f, ef, ex, zer; double d1, d2; dp = 0; /* digits after decimal point */ f = 0; /* sign */ ex = 0; /* exponent */ zer = 1; /* zero */ memset(&a, 0, sizeof(a)); for(;;) { switch(c = *s++) { default: goto bad; case '-': f = 1; case ' ': case '\t': case '+': continue; case '.': dp = 1; continue; case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9': zer = 0; case '0': mpint(&b, c-'0'); mpmul(&a, 10); mpadd(&a, &b); if(dp) dp++; continue; case 'E': case 'e': ex = 0; ef = 0; for(;;) { c = *s++; if(c == '+' || c == ' ' || c == '\t') continue; if(c == '-') { ef = 1; continue; } if(c >= '0' && c <= '9') { ex = ex*10 + (c-'0'); continue; } break; } if(ef) ex = -ex; case 0: break; } break; } if(a.ovf) goto bad; if(zer) { *d = 0; return 0; } if(dp) dp--; dp -= ex; if(dp > 0) { /* * must divide by 10**dp */ if(mptof(&a, &d1)) goto bad; /* * trial exponent of 8**dp * 8 (being between 5 and 10) * should pick up all underflows * in the division of 5**dp. */ d2 = frexp(d1, &ex); d2 = ldexp(d2, ex-3*dp); if(d2 == 0) goto bad; /* * decompose each 10 into 5*2. * create 5**dp in fixed point * and then play with the exponent * for the remaining 2**dp. * note that 5**dp will overflow * with as few as 134 input digits. */ mpint(&a, 1); mppow(&a, 5, dp); if(mptof(&a, &d2)) goto bad; d1 = frexp(d1/d2, &ex); d1 = ldexp(d1, ex-dp); if(d1 == 0) goto bad; } else { /* * must multiply by 10**|dp| -- * just do it in fixed point. */ mppow(&a, 10, -dp); if(mptof(&a, &d1)) goto bad; } if(f) d1 = -d1; *d = d1; return 0; bad: return 1; } /* * convert a to floating in *d * return conversion overflow */ int mptof(Mp *a, double *d) { double f, g; long x, *a1; int i; if(a->ovf) return 1; a1 = a->a; f = ldexp(*a1++, 0); for(i=Mpscale; i<Mpprec*Mpscale; i+=Mpscale) if(x = *a1++) { g = ldexp(x, i); /* * NOTE: the test (g==0) is plan9 * specific. ansi compliant overflow * is signaled by HUGE and errno==ERANGE. * change this for your particular ldexp. */ if(g == 0) return 1; f += g; /* this could bomb! */ } *d = f; return 0; } /* * return a += b */ void mpadd(Mp *a, Mp *b) { int i, c; long x, *a1, *b1; if(b->ovf) a->ovf = 1; if(a->ovf) return; c = 0; a1 = a->a; b1 = b->a; for(i=0; i<Mpprec; i++) { x = *a1 + *b1++ + c; c = 0; if(x >= Mpbase) { x -= Mpbase; c = 1; } *a1++ = x; } a->ovf = c; } /* * return a = c */ void mpint(Mp *a, int c) { memset(a, 0, sizeof(*a)); a->a[0] = c; } /* * return a *= c */ void mpmul(Mp *a, int c) { Mp p; int b; memmove(&p, a, sizeof(p)); if(!(c & 1)) memset(a, 0, sizeof(*a)); c &= ~1; for(b=2; c; b<<=1) { mpadd(&p, &p); if(c & b) { mpadd(a, &p); c &= ~b; } } } /* * return a *= b**e */ void mppow(Mp *a, int b, int e) { int b1; b1 = b*b; b1 = b1*b1; while(e >= 4) { mpmul(a, b1); e -= 4; if(a->ovf) return; } while(e > 0) { mpmul(a, b); e--; } }