ref: 75c92428225428c8fde2d015f010e608a0b12f1d
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--;
}
}