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