Subversion Repositories planix.SVN

/* derived from /netlib/fp/dtoa.c assuming IEEE, Standard C */

/* kudos to dmg@bell-labs.com, gripes to ehg@bell-labs.com */

/* Let x be the exact mathematical number defined by a decimal

 *	string s.  Then atof(s) is the round-nearest-even IEEE

 *	floating point value.

 * Let y be an IEEE floating point value and let s be the string

 *	printed as %.17g.  Then atof(s) is exactly y.

 */

#include <u.h>

#include <libc.h>

static Lock _dtoalk[2];

#define ACQUIRE_DTOA_LOCK(n)	lock(&_dtoalk[n])

#define FREE_DTOA_LOCK(n)	unlock(&_dtoalk[n])

#define PRIVATE_mem ((2000+sizeof(double)-1)/sizeof(double))

static double private_mem[PRIVATE_mem], *pmem_next = private_mem;

#define FLT_ROUNDS	1

#define DBL_DIG		15

#define DBL_MAX_10_EXP	308

#define DBL_MAX_EXP	1024

#define FLT_RADIX	2

#define Storeinc(a,b,c) (*a++ = b << 16 | c & 0xffff)

/* Ten_pmax = floor(P*log(2)/log(5)) */

/* Bletch = (highest power of 2 < DBL_MAX_10_EXP) / 16 */

/* Quick_max = floor((P-1)*log(FLT_RADIX)/log(10) - 1) */

/* Int_max = floor(P*log(FLT_RADIX)/log(10) - 1) */

#define Exp_shift  20

#define Exp_shift1 20

#define Exp_msk1    0x100000

#define Exp_msk11   0x100000

#define Exp_mask  0x7ff00000

#define P 53

#define Bias 1023

#define Emin (-1022)

#define Exp_1  0x3ff00000

#define Exp_11 0x3ff00000

#define Ebits 11

#define Frac_mask  0xfffff

#define Frac_mask1 0xfffff

#define Ten_pmax 22

#define Bletch 0x10

#define Bndry_mask  0xfffff

#define Bndry_mask1 0xfffff

#define LSB 1

#define Sign_bit 0x80000000

#define Log2P 1

#define Tiny0 0

#define Tiny1 1

#define Quick_max 14

#define Int_max 14

#define Avoid_Underflow

#define rounded_product(a,b) a *= b

#define rounded_quotient(a,b) a /= b

#define Big0 (Frac_mask1 | Exp_msk1*(DBL_MAX_EXP+Bias-1))

#define Big1 0xffffffff

#define FFFFFFFF 0xffffffffUL

#define Kmax 15

typedef struct Bigint Bigint;

typedef struct Ulongs Ulongs;

struct Bigint {

	Bigint *next;

	int	k, maxwds, sign, wds;

	unsigned x[1];

};

struct Ulongs {

	ulong	hi;

	ulong	lo;

};

static Bigint *freelist[Kmax+1];

Ulongs

double2ulongs(double d)

{

	FPdbleword dw;

	Ulongs uls;

	dw.x = d;

	uls.hi = dw.hi;

	uls.lo = dw.lo;

	return uls;

}

double

ulongs2double(Ulongs uls)

{

	FPdbleword dw;

	dw.hi = uls.hi;

	dw.lo = uls.lo;

	return dw.x;

}

static Bigint *

Balloc(int k)

{

	int	x;

	Bigint * rv;

	unsigned int	len;

	ACQUIRE_DTOA_LOCK(0);

	if (rv = freelist[k]) {

		freelist[k] = rv->next;

	} else {

		x = 1 << k;

		len = (sizeof(Bigint) + (x - 1) * sizeof(unsigned int) + sizeof(double) -1)

		 / sizeof(double);

		if (pmem_next - private_mem + len <= PRIVATE_mem) {

			rv = (Bigint * )pmem_next;

			pmem_next += len;

		} else

			rv = (Bigint * )malloc(len * sizeof(double));

		rv->k = k;

		rv->maxwds = x;

	}

	FREE_DTOA_LOCK(0);

	rv->sign = rv->wds = 0;

	return rv;

}

static void	

Bfree(Bigint *v)

{

	if (v) {

		ACQUIRE_DTOA_LOCK(0);

		v->next = freelist[v->k];

		freelist[v->k] = v;

		FREE_DTOA_LOCK(0);

	}

}

#define Bcopy(x,y) memcpy((char *)&x->sign, (char *)&y->sign, \

y->wds*sizeof(int) + 2*sizeof(int))

static Bigint *

multadd(Bigint *b, int m, int a)	/* multiply by m and add a */

{

	int	i, wds;

	unsigned int carry, *x, y;

	unsigned int xi, z;

	Bigint * b1;

	wds = b->wds;

	x = b->x;

	i = 0;

	carry = a;

	do {

		xi = *x;

		y = (xi & 0xffff) * m + carry;

		z = (xi >> 16) * m + (y >> 16);

		carry = z >> 16;

		*x++ = (z << 16) + (y & 0xffff);

	} while (++i < wds);

	if (carry) {

		if (wds >= b->maxwds) {

			b1 = Balloc(b->k + 1);

			Bcopy(b1, b);

			Bfree(b);

			b = b1;

		}

		b->x[wds++] = carry;

		b->wds = wds;

	}

	return b;

}

static Bigint *

s2b(const char *s, int nd0, int nd, unsigned int y9)

{

	Bigint * b;

	int	i, k;

	int x, y;

	x = (nd + 8) / 9;

	for (k = 0, y = 1; x > y; y <<= 1, k++) 

		;

	b = Balloc(k);

	b->x[0] = y9;

	b->wds = 1;

	i = 9;

	if (9 < nd0) {

		s += 9;

		do 

			b = multadd(b, 10, *s++ - '0');

		while (++i < nd0);

		s++;

	} else

		s += 10;

	for (; i < nd; i++)

		b = multadd(b, 10, *s++ - '0');

	return b;

}

static int	

hi0bits(register unsigned int x)

{

	register int	k = 0;

	if (!(x & 0xffff0000)) {

		k = 16;

		x <<= 16;

	}

	if (!(x & 0xff000000)) {

		k += 8;

		x <<= 8;

	}

	if (!(x & 0xf0000000)) {

		k += 4;

		x <<= 4;

	}

	if (!(x & 0xc0000000)) {

		k += 2;

		x <<= 2;

	}

	if (!(x & 0x80000000)) {

		k++;

		if (!(x & 0x40000000))

			return 32;

	}

	return k;

}

static int	

lo0bits(unsigned int *y)

{

	register int	k;

	register unsigned int x = *y;

	if (x & 7) {

		if (x & 1)

			return 0;

		if (x & 2) {

			*y = x >> 1;

			return 1;

		}

		*y = x >> 2;

		return 2;

	}

	k = 0;

	if (!(x & 0xffff)) {

		k = 16;

		x >>= 16;

	}

	if (!(x & 0xff)) {

		k += 8;

		x >>= 8;

	}

	if (!(x & 0xf)) {

		k += 4;

		x >>= 4;

	}

	if (!(x & 0x3)) {

		k += 2;

		x >>= 2;

	}

	if (!(x & 1)) {

		k++;

		x >>= 1;

		if (!x & 1)

			return 32;

	}

	*y = x;

	return k;

}

static Bigint *

i2b(int i)

{

	Bigint * b;

	b = Balloc(1);

	b->x[0] = i;

	b->wds = 1;

	return b;

}

static Bigint *

mult(Bigint *a, Bigint *b)

{

	Bigint * c;

	int	k, wa, wb, wc;

	unsigned int * x, *xa, *xae, *xb, *xbe, *xc, *xc0;

	unsigned int y;

	unsigned int carry, z;

	unsigned int z2;

	if (a->wds < b->wds) {

		c = a;

		a = b;

		b = c;

	}

	k = a->k;

	wa = a->wds;

	wb = b->wds;

	wc = wa + wb;

	if (wc > a->maxwds)

		k++;

	c = Balloc(k);

	for (x = c->x, xa = x + wc; x < xa; x++)

		*x = 0;

	xa = a->x;

	xae = xa + wa;

	xb = b->x;

	xbe = xb + wb;

	xc0 = c->x;

	for (; xb < xbe; xb++, xc0++) {

		if (y = *xb & 0xffff) {

			x = xa;

			xc = xc0;

			carry = 0;

			do {

				z = (*x & 0xffff) * y + (*xc & 0xffff) + carry;

				carry = z >> 16;

				z2 = (*x++ >> 16) * y + (*xc >> 16) + carry;

				carry = z2 >> 16;

				Storeinc(xc, z2, z);

			} while (x < xae);

			*xc = carry;

		}

		if (y = *xb >> 16) {

			x = xa;

			xc = xc0;

			carry = 0;

			z2 = *xc;

			do {

				z = (*x & 0xffff) * y + (*xc >> 16) + carry;

				carry = z >> 16;

				Storeinc(xc, z, z2);

				z2 = (*x++ >> 16) * y + (*xc & 0xffff) + carry;

				carry = z2 >> 16;

			} while (x < xae);

			*xc = z2;

		}

	}

	for (xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) 

		;

	c->wds = wc;

	return c;

}

static Bigint *p5s;

static Bigint *

pow5mult(Bigint *b, int k)

{

	Bigint * b1, *p5, *p51;

	int	i;

	static int	p05[3] = { 

		5, 25, 125 	};

	if (i = k & 3)

		b = multadd(b, p05[i-1], 0);

	if (!(k >>= 2))

		return b;

	if (!(p5 = p5s)) {

		/* first time */

		ACQUIRE_DTOA_LOCK(1);

		if (!(p5 = p5s)) {

			p5 = p5s = i2b(625);

			p5->next = 0;

		}

		FREE_DTOA_LOCK(1);

	}

	for (; ; ) {

		if (k & 1) {

			b1 = mult(b, p5);

			Bfree(b);

			b = b1;

		}

		if (!(k >>= 1))

			break;

		if (!(p51 = p5->next)) {

			ACQUIRE_DTOA_LOCK(1);

			if (!(p51 = p5->next)) {

				p51 = p5->next = mult(p5, p5);

				p51->next = 0;

			}

			FREE_DTOA_LOCK(1);

		}

		p5 = p51;

	}

	return b;

}

static Bigint *

lshift(Bigint *b, int k)

{

	int	i, k1, n, n1;

	Bigint * b1;

	unsigned int * x, *x1, *xe, z;

	n = k >> 5;

	k1 = b->k;

	n1 = n + b->wds + 1;

	for (i = b->maxwds; n1 > i; i <<= 1)

		k1++;

	b1 = Balloc(k1);

	x1 = b1->x;

	for (i = 0; i < n; i++)

		*x1++ = 0;

	x = b->x;

	xe = x + b->wds;

	if (k &= 0x1f) {

		k1 = 32 - k;

		z = 0;

		do {

			*x1++ = *x << k | z;

			z = *x++ >> k1;

		} while (x < xe);

		if (*x1 = z)

			++n1;

	} else 

		do

			*x1++ = *x++;

		while (x < xe);

	b1->wds = n1 - 1;

	Bfree(b);

	return b1;

}

static int	

cmp(Bigint *a, Bigint *b)

{

	unsigned int * xa, *xa0, *xb, *xb0;

	int	i, j;

	i = a->wds;

	j = b->wds;

	if (i -= j)

		return i;

	xa0 = a->x;

	xa = xa0 + j;

	xb0 = b->x;

	xb = xb0 + j;

	for (; ; ) {

		if (*--xa != *--xb)

			return * xa < *xb ? -1 : 1;

		if (xa <= xa0)

			break;

	}

	return 0;

}

static Bigint *

diff(Bigint *a, Bigint *b)

{

	Bigint * c;

	int	i, wa, wb;

	unsigned int * xa, *xae, *xb, *xbe, *xc;

	unsigned int borrow, y;

	unsigned int z;

	i = cmp(a, b);

	if (!i) {

		c = Balloc(0);

		c->wds = 1;

		c->x[0] = 0;

		return c;

	}

	if (i < 0) {

		c = a;

		a = b;

		b = c;

		i = 1;

	} else

		i = 0;

	c = Balloc(a->k);

	c->sign = i;

	wa = a->wds;

	xa = a->x;

	xae = xa + wa;

	wb = b->wds;

	xb = b->x;

	xbe = xb + wb;

	xc = c->x;

	borrow = 0;

	do {

		y = (*xa & 0xffff) - (*xb & 0xffff) - borrow;

		borrow = (y & 0x10000) >> 16;

		z = (*xa++ >> 16) - (*xb++ >> 16) - borrow;

		borrow = (z & 0x10000) >> 16;

		Storeinc(xc, z, y);

	} while (xb < xbe);

	while (xa < xae) {

		y = (*xa & 0xffff) - borrow;

		borrow = (y & 0x10000) >> 16;

		z = (*xa++ >> 16) - borrow;

		borrow = (z & 0x10000) >> 16;

		Storeinc(xc, z, y);

	}

	while (!*--xc)

		wa--;

	c->wds = wa;

	return c;

}

static double	

ulp(double x)

{

	ulong L;

	Ulongs uls;

	uls = double2ulongs(x);

	L = (uls.hi & Exp_mask) - (P - 1) * Exp_msk1;

	return ulongs2double((Ulongs){L, 0});

}

static double	

b2d(Bigint *a, int *e)

{

	unsigned *xa, *xa0, w, y, z;

	int	k;

	ulong d0, d1;

	xa0 = a->x;

	xa = xa0 + a->wds;

	y = *--xa;

	k = hi0bits(y);

	*e = 32 - k;

	if (k < Ebits) {

		w = xa > xa0 ? *--xa : 0;

		d1 = y << (32 - Ebits) + k | w >> Ebits - k;

		return ulongs2double((Ulongs){Exp_1 | y >> Ebits - k, d1});

	}

	z = xa > xa0 ? *--xa : 0;

	if (k -= Ebits) {

		d0 = Exp_1 | y << k | z >> 32 - k;

		y = xa > xa0 ? *--xa : 0;

		d1 = z << k | y >> 32 - k;

	} else {

		d0 = Exp_1 | y;

		d1 = z;

	}

	return ulongs2double((Ulongs){d0, d1});

}

static Bigint *

d2b(double d, int *e, int *bits)

{

	Bigint * b;

	int	de, i, k;

	unsigned *x, y, z;

	Ulongs uls;

	b = Balloc(1);

	x = b->x;

	uls = double2ulongs(d);

	z = uls.hi & Frac_mask;

	uls.hi &= 0x7fffffff;		/* clear sign bit, which we ignore */

	de = (int)(uls.hi >> Exp_shift);

	z |= Exp_msk11;

	if (y = uls.lo) {		/* assignment = */

		if (k = lo0bits(&y)) {	/* assignment = */

			x[0] = y | z << 32 - k;

			z >>= k;

		} else

			x[0] = y;

		i = b->wds = (x[1] = z) ? 2 : 1;

	} else {

		k = lo0bits(&z);

		x[0] = z;

		i = b->wds = 1;

		k += 32;

	}

	USED(i);

	*e = de - Bias - (P - 1) + k;

	*bits = P - k;

	return b;

}

static double	

ratio(Bigint *a, Bigint *b)

{

	double	da, db;

	int	k, ka, kb;

	Ulongs uls;

	da = b2d(a, &ka);

	db = b2d(b, &kb);

	k = ka - kb + 32 * (a->wds - b->wds);

	if (k > 0) {

		uls = double2ulongs(da);

		uls.hi += k * Exp_msk1;

		da = ulongs2double(uls);

	} else {

		k = -k;

		uls = double2ulongs(db);

		uls.hi += k * Exp_msk1;

		db = ulongs2double(uls);

	}

	return da / db;

}

static const double

tens[] = {

	1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9,

	1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19,

	1e20, 1e21, 1e22

};

static const double

bigtens[] = { 

	1e16, 1e32, 1e64, 1e128, 1e256 };

static const double tinytens[] = { 

	1e-16, 1e-32, 1e-64, 1e-128,

	9007199254740992.e-256

};

/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */

/* flag unnecessarily.  It leads to a song and dance at the end of strtod. */

#define Scale_Bit 0x10

#define n_bigtens 5

#define NAN_WORD0 0x7ff80000

#define NAN_WORD1 0

static int	

match(const char **sp, char *t)

{

	int	c, d;

	const char * s = *sp;

	while (d = *t++) {

		if ((c = *++s) >= 'A' && c <= 'Z')

			c += 'a' - 'A';

		if (c != d)

			return 0;

	}

	*sp = s + 1;

	return 1;

}

static void	

gethex(double *rvp, const char **sp)

{

	unsigned int c, x[2];

	const char * s;

	int	havedig, udx0, xshift;

	x[0] = x[1] = 0;

	havedig = xshift = 0;

	udx0 = 1;

	s = *sp;

	while (c = *(const unsigned char * )++s) {

		if (c >= '0' && c <= '9')

			c -= '0';

		else if (c >= 'a' && c <= 'f')

			c += 10 - 'a';

		else if (c >= 'A' && c <= 'F')

			c += 10 - 'A';

		else if (c <= ' ') {

			if (udx0 && havedig) {

				udx0 = 0;

				xshift = 1;

			}

			continue;

		} else if (/*(*/ c == ')') {

			*sp = s + 1;

			break;

		} else

			return;	/* invalid form: don't change *sp */

		havedig = 1;

		if (xshift) {

			xshift = 0;

			x[0] = x[1];

			x[1] = 0;

		}

		if (udx0)

			x[0] = (x[0] << 4) | (x[1] >> 28);

		x[1] = (x[1] << 4) | c;

	}

	if ((x[0] &= 0xfffff) || x[1])

		*rvp = ulongs2double((Ulongs){Exp_mask | x[0], x[1]});

}

static int	

quorem(Bigint *b, Bigint *S)

{

	int	n;

	unsigned int * bx, *bxe, q, *sx, *sxe;

	unsigned int borrow, carry, y, ys;

	unsigned int si, z, zs;

	n = S->wds;

	if (b->wds < n)

		return 0;

	sx = S->x;

	sxe = sx + --n;

	bx = b->x;

	bxe = bx + n;

	q = *bxe / (*sxe + 1);	/* ensure q <= true quotient */

	if (q) {

		borrow = 0;

		carry = 0;

		do {

			si = *sx++;

			ys = (si & 0xffff) * q + carry;

			zs = (si >> 16) * q + (ys >> 16);

			carry = zs >> 16;

			y = (*bx & 0xffff) - (ys & 0xffff) - borrow;

			borrow = (y & 0x10000) >> 16;

			z = (*bx >> 16) - (zs & 0xffff) - borrow;

			borrow = (z & 0x10000) >> 16;

			Storeinc(bx, z, y);

		} while (sx <= sxe);

		if (!*bxe) {

			bx = b->x;

			while (--bxe > bx && !*bxe)

				--n;

			b->wds = n;

		}

	}

	if (cmp(b, S) >= 0) {

		q++;

		borrow = 0;

		carry = 0;

		bx = b->x;

		sx = S->x;

		do {

			si = *sx++;

			ys = (si & 0xffff) + carry;

			zs = (si >> 16) + (ys >> 16);

			carry = zs >> 16;

			y = (*bx & 0xffff) - (ys & 0xffff) - borrow;

			borrow = (y & 0x10000) >> 16;

			z = (*bx >> 16) - (zs & 0xffff) - borrow;

			borrow = (z & 0x10000) >> 16;

			Storeinc(bx, z, y);

		} while (sx <= sxe);

		bx = b->x;

		bxe = bx + n;

		if (!*bxe) {

			while (--bxe > bx && !*bxe)

				--n;

			b->wds = n;

		}

	}

	return q;

}

static char	*

rv_alloc(int i)

{

	int	j, k, *r;

	j = sizeof(unsigned int);

	for (k = 0; 

	    sizeof(Bigint) - sizeof(unsigned int) - sizeof(int) + j <= i; 

	    j <<= 1)

		k++;

	r = (int * )Balloc(k);

	*r = k;

	return

	    (char *)(r + 1);

}

static char	*

nrv_alloc(char *s, char **rve, int n)

{

	char	*rv, *t;

	t = rv = rv_alloc(n);

	while (*t = *s++) 

		t++;

	if (rve)

		*rve = t;

	return rv;

}

/* freedtoa(s) must be used to free values s returned by dtoa

 * when MULTIPLE_THREADS is #defined.  It should be used in all cases,

 * but for consistency with earlier versions of dtoa, it is optional

 * when MULTIPLE_THREADS is not defined.

 */

void

freedtoa(char *s)

{

	Bigint * b = (Bigint * )((int *)s - 1);

	b->maxwds = 1 << (b->k = *(int * )b);

	Bfree(b);

}

/* dtoa for IEEE arithmetic (dmg): convert double to ASCII string.

 *

 * Inspired by "How to Print Floating-Point Numbers Accurately" by

 * Guy L. Steele, Jr. and Jon L. White [Proc. ACM SIGPLAN '90, pp. 92-101].

 *

 * Modifications:

 *	1. Rather than iterating, we use a simple numeric overestimate

 *	   to determine k = floor(log10(d)).  We scale relevant

 *	   quantities using O(log2(k)) rather than O(k) multiplications.

 *	2. For some modes > 2 (corresponding to ecvt and fcvt), we don't

 *	   try to generate digits strictly left to right.  Instead, we

 *	   compute with fewer bits and propagate the carry if necessary

 *	   when rounding the final digit up.  This is often faster.

 *	3. Under the assumption that input will be rounded nearest,

 *	   mode 0 renders 1e23 as 1e23 rather than 9.999999999999999e22.

 *	   That is, we allow equality in stopping tests when the

 *	   round-nearest rule will give the same floating-point value

 *	   as would satisfaction of the stopping test with strict

 *	   inequality.

 *	4. We remove common factors of powers of 2 from relevant

 *	   quantities.

 *	5. When converting floating-point integers less than 1e16,

 *	   we use floating-point arithmetic rather than resorting

 *	   to multiple-precision integers.

 *	6. When asked to produce fewer than 15 digits, we first try

 *	   to get by with floating-point arithmetic; we resort to

 *	   multiple-precision integer arithmetic only if we cannot

 *	   guarantee that the floating-point calculation has given

 *	   the correctly rounded result.  For k requested digits and

 *	   "uniformly" distributed input, the probability is

 *	   something like 10^(k-15) that we must resort to the int

 *	   calculation.

 */

char	*

dtoa(double d, int mode, int ndigits, int *decpt, int *sign, char **rve)

{

	/*	Arguments ndigits, decpt, sign are similar to those

	of ecvt and fcvt; trailing zeros are suppressed from

	the returned string.  If not null, *rve is set to point

	to the end of the return value.  If d is +-Infinity or NaN,

	then *decpt is set to 9999.

	mode:

			and rounded to nearest.

		1 ==> like 0, but with Steele & White stopping rule;

			e.g. with IEEE P754 arithmetic , mode 0 gives

			1e23 whereas mode 1 gives 9.999999999999999e22.

		2 ==> max(1,ndigits) significant digits.  This gives a

			return value similar to that of ecvt, except

			that trailing zeros are suppressed.

		3 ==> through ndigits past the decimal point.  This

			gives a return value similar to that from fcvt,

			except that trailing zeros are suppressed, and

			ndigits can be negative.

		4-9 should give the same return values as 2-3, i.e.,

			4 <= mode <= 9 ==> same return as mode

			2 + (mode & 1).  These modes are mainly for

			debugging; often they run slower but sometimes

			faster than modes 2-3.

		4,5,8,9 ==> left-to-right digit generation.