mingw-w64-crt/gdtoa/misc.c - mingw-w64 - Git at Google

 /****************************************************************

 The author of this software is David M. Gay.

 Copyright (C) 1998, 1999 by Lucent Technologies
 All Rights Reserved

 Permission to use, copy, modify, and distribute this software and
 its documentation for any purpose and without fee is hereby
 granted, provided that the above copyright notice appear in all
 copies and that both that the copyright notice and this
 permission notice and warranty disclaimer appear in supporting
 documentation, and that the name of Lucent or any of its entities
 not be used in advertising or publicity pertaining to
 distribution of the software without specific, written prior
 permission.

 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
 THIS SOFTWARE.

 ****************************************************************/

 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  * with " at " changed at "@" and " dot " changed to ".").	*/


 #if defined(__MINGW32__) || defined(__MINGW64__)
 /* we have to include windows.h before gdtoa
    headers, otherwise defines cause conflicts. */
 #ifndef WIN32_LEAN_AND_MEAN
 #define WIN32_LEAN_AND_MEAN
 #endif
 #include <windows.h>

 #define NLOCKS 2

 #ifdef USE_WIN32_SL
 /* Use spin locks. */
 static long dtoa_sl[NLOCKS];

 #define ACQUIRE_DTOA_LOCK(n) \
   while (InterlockedCompareExchange (&dtoa_sl[n], 1, 0) != 0) \
      Sleep (0);
 #define FREE_DTOA_LOCK(n) InterlockedExchange (&dtoa_sl[n], 0);

 #else	/* USE_WIN32_SL */

 #include <stdlib.h>
 static CRITICAL_SECTION dtoa_CritSec[NLOCKS];
 static long dtoa_CS_init = 0;
 /*
    1 = initializing
    2 = initialized
    3 = deleted
 */
 static void dtoa_lock_cleanup (void)
 {
 	long last_CS_init = InterlockedExchange (&dtoa_CS_init,3);
 	if (2 == last_CS_init) {
 		int i;
 		for (i = 0; i < NLOCKS; i++)
 			DeleteCriticalSection (&dtoa_CritSec[i]);
 	}
 }

 static void dtoa_lock (unsigned int n)
 {
 	if (2 == dtoa_CS_init) {
 		EnterCriticalSection (&dtoa_CritSec[n]);
 		return;
 	}
 	else if (0 == dtoa_CS_init) {
 		long last_CS_init = InterlockedExchange (&dtoa_CS_init, 1);
 		if (0 == last_CS_init) {
 			int i;
 			for (i = 0; i < NLOCKS;  i++)
 				InitializeCriticalSection (&dtoa_CritSec[i]);
 			atexit (dtoa_lock_cleanup);
 			dtoa_CS_init = 2;
 		}
 		else if (2 == last_CS_init)
 			dtoa_CS_init = 2;
 	}
 	/*  Another thread is initializing. Wait. */
 	while (1 == dtoa_CS_init)
 		Sleep (1);

 	/* It had better be initialized now. */
 	if (2 == dtoa_CS_init)
 		EnterCriticalSection(&dtoa_CritSec[n]);
 }

 static void dtoa_unlock (unsigned int n)
 {
 	if (2 == dtoa_CS_init)
 		LeaveCriticalSection (&dtoa_CritSec[n]);
 }

 #define ACQUIRE_DTOA_LOCK(n) dtoa_lock(n)
 #define FREE_DTOA_LOCK(n) dtoa_unlock(n)
 #endif	/* USE_WIN32_SL */

 #endif	/* __MINGW32__ / __MINGW64__ */

 #include "gdtoaimp.h"

 static Bigint *freelist[Kmax+1];
 #ifndef Omit_Private_Memory
 #ifndef PRIVATE_MEM
 #define PRIVATE_MEM 2304
 #endif
 #define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
 static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
 #endif

 Bigint *Balloc (int k)
 {
 	int x;
 	Bigint *rv;
 #ifndef Omit_Private_Memory
 	unsigned int len;
 #endif

 	ACQUIRE_DTOA_LOCK(0);
 	/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
 	/* but this case seems very unlikely. */
 	if (k <= Kmax && (rv = freelist[k]) !=0) {
 		freelist[k] = rv->next;
 	}
 	else {
 		x = 1 << k;
 #ifdef Omit_Private_Memory
 		rv = (Bigint *)MALLOC(sizeof(Bigint) + (x-1)*sizeof(ULong));
 		if (rv == NULL)
 			return NULL;
 #else
 		len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
 			/sizeof(double);
 		if (k <= Kmax
 		    && (size_t) (pmem_next - private_mem + len) <= PRIVATE_mem)
 		{
 			rv = (Bigint*)pmem_next;
 			pmem_next += len;
 		}
 		else
 		{
 			rv = (Bigint*)MALLOC(len*sizeof(double));
 			if (rv == NULL)
 				return NULL;
 		}
 #endif
 		rv->k = k;
 		rv->maxwds = x;
 	}
 	FREE_DTOA_LOCK(0);
 	rv->sign = rv->wds = 0;
 	return rv;
 }

 void Bfree (Bigint *v)
 {
 	if (v) {
 		if (v->k > Kmax)
 			free((void*)v);
 		else {
 			ACQUIRE_DTOA_LOCK(0);
 			v->next = freelist[v->k];
 			freelist[v->k] = v;
 			FREE_DTOA_LOCK(0);
 		}
 	}
 }

 /* lo0bits():  Shift y so lowest bit is 1 and return the
  *		 number of bits y was shifted.
  * With GCC, we use an inline wrapper for __builtin_clz()
  */
 #ifndef __GNUC__
 int lo0bits (ULong *y)
 {
 	int k;
 	ULong 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)
 			return 32;
 	}
 	*y = x;
 	return k;
 }
 #endif	/* __GNUC__ */

 Bigint *multadd (Bigint *b, int m, int a)	/* multiply by m and add a */
 {
 	int i, wds;
 #ifdef ULLong
 	ULong *x;
 	ULLong carry, y;
 #else
 	ULong carry, *x, y;
 #ifdef Pack_32
 	ULong xi, z;
 #endif
 #endif
 	Bigint *b1;

 	wds = b->wds;
 	x = b->x;
 	i = 0;
 	carry = a;
 	do {
 #ifdef ULLong
 		y = *x * (ULLong)m + carry;
 		carry = y >> 32;
 		*x++ = y & 0xffffffffUL;
 #else
 #ifdef Pack_32
 		xi = *x;
 		y = (xi & 0xffff) * m + carry;
 		z = (xi >> 16) * m + (y >> 16);
 		carry = z >> 16;
 		*x++ = (z << 16) + (y & 0xffff);
 #else
 		y = *x * m + carry;
 		carry = y >> 16;
 		*x++ = y & 0xffff;
 #endif
 #endif
 	} while(++i < wds);
 	if (carry) {
 		if (wds >= b->maxwds) {
 			b1 = Balloc(b->k+1);
 			if (b1 == NULL)
 				return NULL;
 			Bcopy(b1, b);
 			Bfree(b);
 			b = b1;
 		}
 		b->x[wds++] = carry;
 		b->wds = wds;
 	}
 	return b;
 }

 /* hi0bits();
  * With GCC, we use an inline wrapper for __builtin_clz()
  */
 #ifndef __GNUC__
 int hi0bits_D2A (ULong x)
 {
 	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;
 }
 #endif	/* __GNUC__ */

 Bigint *i2b (int i)
 {
 	Bigint *b;

 	b = Balloc(1);
 	if (b == NULL)
 		return NULL;
 	b->x[0] = i;
 	b->wds = 1;
 	return b;
 }

 Bigint *mult (Bigint *a, Bigint *b)
 {
 	Bigint *c;
 	int k, wa, wb, wc;
 	ULong *x, *xa, *xae, *xb, *xbe, *xc, *xc0;
 	ULong y;
 #ifdef ULLong
 	ULLong carry, z;
 #else
 	ULong carry, z;
 #ifdef Pack_32
 	ULong z2;
 #endif
 #endif

 	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);
 	if (c == NULL)
 		return NULL;
 	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;
 #ifdef ULLong
 	for(; xb < xbe; xc0++) {
 		if ( (y = *xb++) !=0) {
 			x = xa;
 			xc = xc0;
 			carry = 0;
 			do {
 				z = *x++ * (ULLong)y + *xc + carry;
 				carry = z >> 32;
 				*xc++ = z & 0xffffffffUL;
 			} while(x < xae);
 			*xc = carry;
 		}
 	}
 #else
 #ifdef Pack_32
 	for(; xb < xbe; xb++, xc0++) {
 		if ( (y = *xb & 0xffff) !=0) {
 			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) !=0) {
 			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;
 		}
 	}
 #else
 	for(; xb < xbe; xc0++) {
 		if ( (y = *xb++) !=0) {
 			x = xa;
 			xc = xc0;
 			carry = 0;
 			do {
 				z = *x++ * y + *xc + carry;
 				carry = z >> 16;
 				*xc++ = z & 0xffff;
 			} while(x < xae);
 			*xc = carry;
 		}
 	}
 #endif
 #endif
 	for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
 	c->wds = wc;
 	return c;
 }

 static Bigint *p5s;

 Bigint *pow5mult (Bigint *b, int k)
 {
 	Bigint *b1, *p5, *p51;
 	int i;
 	static int p05[3] = { 5, 25, 125 };

 	if ( (i = k & 3) !=0)
 	{
 		b = multadd(b, p05[i-1], 0);
 		if (b == NULL)
 		return NULL;
 	}
 	if (!(k >>= 2))
 		return b;
 	if ((p5 = p5s) == 0) {
 		/* first time */
 #ifdef MULTIPLE_THREADS
 		ACQUIRE_DTOA_LOCK(1);
 		if (!(p5 = p5s)) {
 			p5 = p5s = i2b(625);
 			if (p5 == NULL)
 				return NULL;
 			p5->next = 0;
 		}
 		FREE_DTOA_LOCK(1);
 #else
 		p5 = p5s = i2b(625);
 		if (p5 == NULL)
 			return NULL;
 		p5->next = 0;
 #endif
 	}
 	for(;;) {
 		if (k & 1) {
 			b1 = mult(b, p5);
 			if (b1 == NULL)
 				return NULL;
 			Bfree(b);
 			b = b1;
 		}
 		if (!(k >>= 1))
 			break;
 		if ((p51 = p5->next) == 0) {
 #ifdef MULTIPLE_THREADS
 			ACQUIRE_DTOA_LOCK(1);
 			if (!(p51 = p5->next)) {
 				p51 = p5->next = mult(p5,p5);
 				if (p51 == NULL)
 					return NULL;
 				p51->next = 0;
 			}
 			FREE_DTOA_LOCK(1);
 #else
 			p51 = p5->next = mult(p5,p5);
 			if (p51 == NULL)
 				return NULL;
 			p51->next = 0;
 #endif
 		}
 		p5 = p51;
 	}
 	return b;
 }

 Bigint *lshift (Bigint *b, int k)
 {
 	int i, k1, n, n1;
 	Bigint *b1;
 	ULong *x, *x1, *xe, z;

 	n = k >> kshift;
 	k1 = b->k;
 	n1 = n + b->wds + 1;
 	for(i = b->maxwds; n1 > i; i <<= 1)
 		k1++;
 	b1 = Balloc(k1);
 	if (b1 == NULL)
 		return NULL;
 	x1 = b1->x;
 	for(i = 0; i < n; i++)
 		*x1++ = 0;
 	x = b->x;
 	xe = x + b->wds;
 	if (k &= kmask) {
 #ifdef Pack_32
 		k1 = 32 - k;
 		z = 0;
 		do {
 			*x1++ = *x << k | z;
 			z = *x++ >> k1;
 		} while(x < xe);
 		if ((*x1 = z) !=0)
 			++n1;
 #else
 		k1 = 16 - k;
 		z = 0;
 		do {
 			*x1++ = *x << k  & 0xffff | z;
 			z = *x++ >> k1;
 		} while(x < xe);
 		if (*x1 = z)
 			++n1;
 #endif
 	}
 	else do
 		*x1++ = *x++;
 		while(x < xe);
 	b1->wds = n1 - 1;
 	Bfree(b);
 	return b1;
 }

 int cmp (Bigint *a, Bigint *b)
 {
 	ULong *xa, *xa0, *xb, *xb0;
 	int i, j;

 	i = a->wds;
 	j = b->wds;
 #ifdef DEBUG
 	if (i > 1 && !a->x[i-1])
 		Bug("cmp called with a->x[a->wds-1] == 0");
 	if (j > 1 && !b->x[j-1])
 		Bug("cmp called with b->x[b->wds-1] == 0");
 #endif
 	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;
 }

 Bigint *diff (Bigint *a, Bigint *b)
 {
 	Bigint *c;
 	int i, wa, wb;
 	ULong *xa, *xae, *xb, *xbe, *xc;
 #ifdef ULLong
 	ULLong borrow, y;
 #else
 	ULong borrow, y;
 #ifdef Pack_32
 	ULong z;
 #endif
 #endif

 	i = cmp(a,b);
 	if (!i) {
 		c = Balloc(0);
 		if (c == NULL)
 			return NULL;
 		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);
 	if (c == NULL)
 		return NULL;
 	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;
 #ifdef ULLong
 	do {
 		y = (ULLong)*xa++ - *xb++ - borrow;
 		borrow = y >> 32 & 1UL;
 		*xc++ = y & 0xffffffffUL;
 	} while(xb < xbe);
 	while(xa < xae) {
 		y = *xa++ - borrow;
 		borrow = y >> 32 & 1UL;
 		*xc++ = y & 0xffffffffUL;
 	}
 #else
 #ifdef Pack_32
 	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);
 	}
 #else
 	do {
 		y = *xa++ - *xb++ - borrow;
 		borrow = (y & 0x10000) >> 16;
 		*xc++ = y & 0xffff;
 	} while(xb < xbe);
 	while(xa < xae) {
 		y = *xa++ - borrow;
 		borrow = (y & 0x10000) >> 16;
 		*xc++ = y & 0xffff;
 	}
 #endif
 #endif
 	while(!*--xc)
 		wa--;
 	c->wds = wa;
 	return c;
 }

 double b2d (Bigint *a, int *e)
 {
 	ULong *xa, *xa0, w, y, z;
 	int k;
 	union _dbl_union d;
 #define d0 word0(&d)
 #define d1 word1(&d)

 	xa0 = a->x;
 	xa = xa0 + a->wds;
 	y = *--xa;
 #ifdef DEBUG
 	if (!y) Bug("zero y in b2d");
 #endif
 	k = hi0bits(y);
 	*e = 32 - k;
 #ifdef Pack_32
 	if (k < Ebits) {
 		d0 = Exp_1 | y >> (Ebits - k);
 		w = xa > xa0 ? *--xa : 0;
 		d1 = y << ((32-Ebits) + k) | w >> (Ebits - k);
 		goto ret_d;
 	}
 	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;
 	}
 #else
 	if (k < Ebits + 16) {
 		z = xa > xa0 ? *--xa : 0;
 		d0 = Exp_1 | y << k - Ebits | z >> Ebits + 16 - k;
 		w = xa > xa0 ? *--xa : 0;
 		y = xa > xa0 ? *--xa : 0;
 		d1 = z << k + 16 - Ebits | w << k - Ebits | y >> 16 + Ebits - k;
 		goto ret_d;
 	}
 	z = xa > xa0 ? *--xa : 0;
 	w = xa > xa0 ? *--xa : 0;
 	k -= Ebits + 16;
 	d0 = Exp_1 | y << k + 16 | z << k | w >> 16 - k;
 	y = xa > xa0 ? *--xa : 0;
 	d1 = w << k + 16 | y << k;
 #endif
  ret_d:
 	return dval(&d);
 #undef d0
 #undef d1
 }

 Bigint *d2b (double dd, int *e, int *bits)
 {
 	Bigint *b;
 	union _dbl_union d;
 #ifndef Sudden_Underflow
 	int i;
 #endif
 	int de, k;
 	ULong *x, y, z;
 #define d0 word0(&d)
 #define d1 word1(&d)
 	d.d = dd;

 #ifdef Pack_32
 	b = Balloc(1);
 #else
 	b = Balloc(2);
 #endif
 	if (b == NULL)
 		return NULL;
 	x = b->x;

 	z = d0 & Frac_mask;
 	d0 &= 0x7fffffff;	/* clear sign bit, which we ignore */
 #ifdef Sudden_Underflow
 	de = (int)(d0 >> Exp_shift);
 	z |= Exp_msk11;
 #else
 	if ( (de = (int)(d0 >> Exp_shift)) !=0)
 		z |= Exp_msk1;
 #endif
 #ifdef Pack_32
 	if ( (y = d1) !=0) {
 		if ( (k = lo0bits(&y)) !=0) {
 			x[0] = y | z << (32 - k);
 			z >>= k;
 		}
 		else
 			x[0] = y;
 #ifndef Sudden_Underflow
 		i =
 #endif
 		     b->wds = (x[1] = z) !=0 ? 2 : 1;
 	}
 	else {
 		k = lo0bits(&z);
 		x[0] = z;
 #ifndef Sudden_Underflow
 		i =
 #endif
 		    b->wds = 1;
 		k += 32;
 	}
 #else
 	if ( (y = d1) !=0) {
 		if ( (k = lo0bits(&y)) !=0)
 			if (k >= 16) {
 				x[0] = y | z << 32 - k & 0xffff;
 				x[1] = z >> k - 16 & 0xffff;
 				x[2] = z >> k;
 				i = 2;
 			}
 			else {
 				x[0] = y & 0xffff;
 				x[1] = y >> 16 | z << 16 - k & 0xffff;
 				x[2] = z >> k & 0xffff;
 				x[3] = z >> k+16;
 				i = 3;
 			}
 		else {
 			x[0] = y & 0xffff;
 			x[1] = y >> 16;
 			x[2] = z & 0xffff;
 			x[3] = z >> 16;
 			i = 3;
 		}
 	}
 	else {
 #ifdef DEBUG
 		if (!z)
 			Bug("Zero passed to d2b");
 #endif
 		k = lo0bits(&z);
 		if (k >= 16) {
 			x[0] = z;
 			i = 0;
 		}
 		else {
 			x[0] = z & 0xffff;
 			x[1] = z >> 16;
 			i = 1;
 		}
 		k += 32;
 	}
 	while(!x[i])
 		--i;
 	b->wds = i + 1;
 #endif
 #ifndef Sudden_Underflow
 	if (de) {
 #endif
 		*e = de - Bias - (P-1) + k;
 		*bits = P - k;
 #ifndef Sudden_Underflow
 	}
 	else {
 		*e = de - Bias - (P-1) + 1 + k;
 #ifdef Pack_32
 		*bits = 32*i - hi0bits(x[i-1]);
 #else
 		*bits = (i+2)*16 - hi0bits(x[i]);
 #endif
 	}
 #endif
 	return b;
 #undef d0
 #undef d1
 }

 const double
 bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
 const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };

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

 char *strcp_D2A (char *a, const char *b)
 {
 	while((*a = *b++))
 		a++;
 	return a;
 }

 #ifdef NO_STRING_H
 void *memcpy_D2A (void *a1, void *b1, size_t len)
 {
 	char *a = (char*)a1, *ae = a + len;
 	char *b = (char*)b1, *a0 = a;
 	while(a < ae)
 		*a++ = *b++;
 	return a0;
 }
 #endif /* NO_STRING_H */
	/****************************************************************

	The author of this software is David M. Gay.

	Copyright (C) 1998, 1999 by Lucent Technologies
	All Rights Reserved

	Permission to use, copy, modify, and distribute this software and
	its documentation for any purpose and without fee is hereby
	granted, provided that the above copyright notice appear in all
	copies and that both that the copyright notice and this
	permission notice and warranty disclaimer appear in supporting
	documentation, and that the name of Lucent or any of its entities
	not be used in advertising or publicity pertaining to
	distribution of the software without specific, written prior
	permission.

	LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
	INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
	IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
	SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
	WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
	IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
	ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
	THIS SOFTWARE.

	****************************************************************/

	/* Please send bug reports to David M. Gay (dmg at acm dot org,
	* with " at " changed at "@" and " dot " changed to "."). */


	#if defined(__MINGW32__) \|\| defined(__MINGW64__)
	/* we have to include windows.h before gdtoa
	headers, otherwise defines cause conflicts. */
	#ifndef WIN32_LEAN_AND_MEAN
	#define WIN32_LEAN_AND_MEAN
	#endif
	#include <windows.h>

	#define NLOCKS 2

	#ifdef USE_WIN32_SL
	/* Use spin locks. */
	static long dtoa_sl[NLOCKS];

	#define ACQUIRE_DTOA_LOCK(n) \
	while (InterlockedCompareExchange (&dtoa_sl[n], 1, 0) != 0) \
	Sleep (0);
	#define FREE_DTOA_LOCK(n) InterlockedExchange (&dtoa_sl[n], 0);

	#else /* USE_WIN32_SL */

	#include <stdlib.h>
	static CRITICAL_SECTION dtoa_CritSec[NLOCKS];
	static long dtoa_CS_init = 0;
	/*
	1 = initializing
	2 = initialized
	3 = deleted
	*/
	static void dtoa_lock_cleanup (void)
	{
	long last_CS_init = InterlockedExchange (&dtoa_CS_init,3);
	if (2 == last_CS_init) {
	int i;
	for (i = 0; i < NLOCKS; i++)
	DeleteCriticalSection (&dtoa_CritSec[i]);
	}
	}

	static void dtoa_lock (unsigned int n)
	{
	if (2 == dtoa_CS_init) {
	EnterCriticalSection (&dtoa_CritSec[n]);
	return;
	}
	else if (0 == dtoa_CS_init) {
	long last_CS_init = InterlockedExchange (&dtoa_CS_init, 1);
	if (0 == last_CS_init) {
	int i;
	for (i = 0; i < NLOCKS; i++)
	InitializeCriticalSection (&dtoa_CritSec[i]);
	atexit (dtoa_lock_cleanup);
	dtoa_CS_init = 2;
	}
	else if (2 == last_CS_init)
	dtoa_CS_init = 2;
	}
	/* Another thread is initializing. Wait. */
	while (1 == dtoa_CS_init)
	Sleep (1);

	/* It had better be initialized now. */
	if (2 == dtoa_CS_init)
	EnterCriticalSection(&dtoa_CritSec[n]);
	}

	static void dtoa_unlock (unsigned int n)
	{
	if (2 == dtoa_CS_init)
	LeaveCriticalSection (&dtoa_CritSec[n]);
	}

	#define ACQUIRE_DTOA_LOCK(n) dtoa_lock(n)
	#define FREE_DTOA_LOCK(n) dtoa_unlock(n)
	#endif /* USE_WIN32_SL */

	#endif /* __MINGW32__ / __MINGW64__ */

	#include "gdtoaimp.h"

	static Bigint *freelist[Kmax+1];
	#ifndef Omit_Private_Memory
	#ifndef PRIVATE_MEM
	#define PRIVATE_MEM 2304
	#endif
	#define PRIVATE_mem ((PRIVATE_MEM+sizeof(double)-1)/sizeof(double))
	static double private_mem[PRIVATE_mem], *pmem_next = private_mem;
	#endif

	Bigint *Balloc (int k)
	{
	int x;
	Bigint *rv;
	#ifndef Omit_Private_Memory
	unsigned int len;
	#endif

	ACQUIRE_DTOA_LOCK(0);
	/* The k > Kmax case does not need ACQUIRE_DTOA_LOCK(0), */
	/* but this case seems very unlikely. */
	if (k <= Kmax && (rv = freelist[k]) !=0) {
	freelist[k] = rv->next;
	}
	else {
	x = 1 << k;
	#ifdef Omit_Private_Memory
	rv = (Bigint )MALLOC(sizeof(Bigint) + (x-1)sizeof(ULong));
	if (rv == NULL)
	return NULL;
	#else
	len = (sizeof(Bigint) + (x-1)*sizeof(ULong) + sizeof(double) - 1)
	/sizeof(double);
	if (k <= Kmax
	&& (size_t) (pmem_next - private_mem + len) <= PRIVATE_mem)
	{
	rv = (Bigint*)pmem_next;
	pmem_next += len;
	}
	else
	{
	rv = (Bigint)MALLOC(lensizeof(double));
	if (rv == NULL)
	return NULL;
	}
	#endif
	rv->k = k;
	rv->maxwds = x;
	}
	FREE_DTOA_LOCK(0);
	rv->sign = rv->wds = 0;
	return rv;
	}

	void Bfree (Bigint *v)
	{
	if (v) {
	if (v->k > Kmax)
	free((void*)v);
	else {
	ACQUIRE_DTOA_LOCK(0);
	v->next = freelist[v->k];
	freelist[v->k] = v;
	FREE_DTOA_LOCK(0);
	}
	}
	}

	/* lo0bits(): Shift y so lowest bit is 1 and return the
	* number of bits y was shifted.
	* With GCC, we use an inline wrapper for __builtin_clz()
	*/
	#ifndef __GNUC__
	int lo0bits (ULong *y)
	{
	int k;
	ULong 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)
	return 32;
	}
	*y = x;
	return k;
	}
	#endif /* __GNUC__ */

	Bigint multadd (Bigint b, int m, int a) /* multiply by m and add a */
	{
	int i, wds;
	#ifdef ULLong
	ULong *x;
	ULLong carry, y;
	#else
	ULong carry, *x, y;
	#ifdef Pack_32
	ULong xi, z;
	#endif
	#endif
	Bigint *b1;

	wds = b->wds;
	x = b->x;
	i = 0;
	carry = a;
	do {
	#ifdef ULLong
	y = x (ULLong)m + carry;
	carry = y >> 32;
	*x++ = y & 0xffffffffUL;
	#else
	#ifdef Pack_32
	xi = *x;
	y = (xi & 0xffff) * m + carry;
	z = (xi >> 16) * m + (y >> 16);
	carry = z >> 16;
	*x++ = (z << 16) + (y & 0xffff);
	#else
	y = x m + carry;
	carry = y >> 16;
	*x++ = y & 0xffff;
	#endif
	#endif
	} while(++i < wds);
	if (carry) {
	if (wds >= b->maxwds) {
	b1 = Balloc(b->k+1);
	if (b1 == NULL)
	return NULL;
	Bcopy(b1, b);
	Bfree(b);
	b = b1;
	}
	b->x[wds++] = carry;
	b->wds = wds;
	}
	return b;
	}

	/* hi0bits();
	* With GCC, we use an inline wrapper for __builtin_clz()
	*/
	#ifndef __GNUC__
	int hi0bits_D2A (ULong x)
	{
	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;
	}
	#endif /* __GNUC__ */

	Bigint *i2b (int i)
	{
	Bigint *b;

	b = Balloc(1);
	if (b == NULL)
	return NULL;
	b->x[0] = i;
	b->wds = 1;
	return b;
	}

	Bigint mult (Bigint a, Bigint *b)
	{
	Bigint *c;
	int k, wa, wb, wc;
	ULong x, xa, xae, xb, xbe, xc, *xc0;
	ULong y;
	#ifdef ULLong
	ULLong carry, z;
	#else
	ULong carry, z;
	#ifdef Pack_32
	ULong z2;
	#endif
	#endif

	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);
	if (c == NULL)
	return NULL;
	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;
	#ifdef ULLong
	for(; xb < xbe; xc0++) {
	if ( (y = *xb++) !=0) {
	x = xa;
	xc = xc0;
	carry = 0;
	do {
	z = x++ (ULLong)y + *xc + carry;
	carry = z >> 32;
	*xc++ = z & 0xffffffffUL;
	} while(x < xae);
	*xc = carry;
	}
	}
	#else
	#ifdef Pack_32
	for(; xb < xbe; xb++, xc0++) {
	if ( (y = *xb & 0xffff) !=0) {
	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) !=0) {
	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;
	}
	}
	#else
	for(; xb < xbe; xc0++) {
	if ( (y = *xb++) !=0) {
	x = xa;
	xc = xc0;
	carry = 0;
	do {
	z = x++ y + *xc + carry;
	carry = z >> 16;
	*xc++ = z & 0xffff;
	} while(x < xae);
	*xc = carry;
	}
	}
	#endif
	#endif
	for(xc0 = c->x, xc = xc0 + wc; wc > 0 && !*--xc; --wc) ;
	c->wds = wc;
	return c;
	}

	static Bigint *p5s;

	Bigint pow5mult (Bigint b, int k)
	{
	Bigint b1, p5, *p51;
	int i;
	static int p05[3] = { 5, 25, 125 };

	if ( (i = k & 3) !=0)
	{
	b = multadd(b, p05[i-1], 0);
	if (b == NULL)
	return NULL;
	}
	if (!(k >>= 2))
	return b;
	if ((p5 = p5s) == 0) {
	/* first time */
	#ifdef MULTIPLE_THREADS
	ACQUIRE_DTOA_LOCK(1);
	if (!(p5 = p5s)) {
	p5 = p5s = i2b(625);
	if (p5 == NULL)
	return NULL;
	p5->next = 0;
	}
	FREE_DTOA_LOCK(1);
	#else
	p5 = p5s = i2b(625);
	if (p5 == NULL)
	return NULL;
	p5->next = 0;
	#endif
	}
	for(;;) {
	if (k & 1) {
	b1 = mult(b, p5);
	if (b1 == NULL)
	return NULL;
	Bfree(b);
	b = b1;
	}
	if (!(k >>= 1))
	break;
	if ((p51 = p5->next) == 0) {
	#ifdef MULTIPLE_THREADS
	ACQUIRE_DTOA_LOCK(1);
	if (!(p51 = p5->next)) {
	p51 = p5->next = mult(p5,p5);
	if (p51 == NULL)
	return NULL;
	p51->next = 0;
	}
	FREE_DTOA_LOCK(1);
	#else
	p51 = p5->next = mult(p5,p5);
	if (p51 == NULL)
	return NULL;
	p51->next = 0;
	#endif
	}
	p5 = p51;
	}
	return b;
	}

	Bigint lshift (Bigint b, int k)
	{
	int i, k1, n, n1;
	Bigint *b1;
	ULong x, x1, *xe, z;

	n = k >> kshift;
	k1 = b->k;
	n1 = n + b->wds + 1;
	for(i = b->maxwds; n1 > i; i <<= 1)
	k1++;
	b1 = Balloc(k1);
	if (b1 == NULL)
	return NULL;
	x1 = b1->x;
	for(i = 0; i < n; i++)
	*x1++ = 0;
	x = b->x;
	xe = x + b->wds;
	if (k &= kmask) {
	#ifdef Pack_32
	k1 = 32 - k;
	z = 0;
	do {
	x1++ = x << k \| z;
	z = *x++ >> k1;
	} while(x < xe);
	if ((*x1 = z) !=0)
	++n1;
	#else
	k1 = 16 - k;
	z = 0;
	do {
	x1++ = x << k & 0xffff \| z;
	z = *x++ >> k1;
	} while(x < xe);
	if (*x1 = z)
	++n1;
	#endif
	}
	else do
	x1++ = x++;
	while(x < xe);
	b1->wds = n1 - 1;
	Bfree(b);
	return b1;
	}

	int cmp (Bigint a, Bigint b)
	{
	ULong xa, xa0, xb, xb0;
	int i, j;

	i = a->wds;
	j = b->wds;
	#ifdef DEBUG
	if (i > 1 && !a->x[i-1])
	Bug("cmp called with a->x[a->wds-1] == 0");
	if (j > 1 && !b->x[j-1])
	Bug("cmp called with b->x[b->wds-1] == 0");
	#endif
	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;
	}

	Bigint diff (Bigint a, Bigint *b)
	{
	Bigint *c;
	int i, wa, wb;
	ULong xa, xae, xb, xbe, *xc;
	#ifdef ULLong
	ULLong borrow, y;
	#else
	ULong borrow, y;
	#ifdef Pack_32
	ULong z;
	#endif
	#endif

	i = cmp(a,b);
	if (!i) {
	c = Balloc(0);
	if (c == NULL)
	return NULL;
	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);
	if (c == NULL)
	return NULL;
	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;
	#ifdef ULLong
	do {
	y = (ULLong)xa++ - xb++ - borrow;
	borrow = y >> 32 & 1UL;
	*xc++ = y & 0xffffffffUL;
	} while(xb < xbe);
	while(xa < xae) {
	y = *xa++ - borrow;
	borrow = y >> 32 & 1UL;
	*xc++ = y & 0xffffffffUL;
	}
	#else
	#ifdef Pack_32
	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);
	}
	#else
	do {
	y = xa++ - xb++ - borrow;
	borrow = (y & 0x10000) >> 16;
	*xc++ = y & 0xffff;
	} while(xb < xbe);
	while(xa < xae) {
	y = *xa++ - borrow;
	borrow = (y & 0x10000) >> 16;
	*xc++ = y & 0xffff;
	}
	#endif
	#endif
	while(!*--xc)
	wa--;
	c->wds = wa;
	return c;
	}

	double b2d (Bigint a, int e)
	{
	ULong xa, xa0, w, y, z;
	int k;
	union _dbl_union d;
	#define d0 word0(&d)
	#define d1 word1(&d)

	xa0 = a->x;
	xa = xa0 + a->wds;
	y = *--xa;
	#ifdef DEBUG
	if (!y) Bug("zero y in b2d");
	#endif
	k = hi0bits(y);
	*e = 32 - k;
	#ifdef Pack_32
	if (k < Ebits) {
	d0 = Exp_1 \| y >> (Ebits - k);
	w = xa > xa0 ? *--xa : 0;
	d1 = y << ((32-Ebits) + k) \| w >> (Ebits - k);
	goto ret_d;
	}
	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;
	}
	#else
	if (k < Ebits + 16) {
	z = xa > xa0 ? *--xa : 0;
	d0 = Exp_1 \| y << k - Ebits \| z >> Ebits + 16 - k;
	w = xa > xa0 ? *--xa : 0;
	y = xa > xa0 ? *--xa : 0;
	d1 = z << k + 16 - Ebits \| w << k - Ebits \| y >> 16 + Ebits - k;
	goto ret_d;
	}
	z = xa > xa0 ? *--xa : 0;
	w = xa > xa0 ? *--xa : 0;
	k -= Ebits + 16;
	d0 = Exp_1 \| y << k + 16 \| z << k \| w >> 16 - k;
	y = xa > xa0 ? *--xa : 0;
	d1 = w << k + 16 \| y << k;
	#endif
	ret_d:
	return dval(&d);
	#undef d0
	#undef d1
	}

	Bigint d2b (double dd, int e, int *bits)
	{
	Bigint *b;
	union _dbl_union d;
	#ifndef Sudden_Underflow
	int i;
	#endif
	int de, k;
	ULong *x, y, z;
	#define d0 word0(&d)
	#define d1 word1(&d)
	d.d = dd;

	#ifdef Pack_32
	b = Balloc(1);
	#else
	b = Balloc(2);
	#endif
	if (b == NULL)
	return NULL;
	x = b->x;

	z = d0 & Frac_mask;
	d0 &= 0x7fffffff; /* clear sign bit, which we ignore */
	#ifdef Sudden_Underflow
	de = (int)(d0 >> Exp_shift);
	z \|= Exp_msk11;
	#else
	if ( (de = (int)(d0 >> Exp_shift)) !=0)
	z \|= Exp_msk1;
	#endif
	#ifdef Pack_32
	if ( (y = d1) !=0) {
	if ( (k = lo0bits(&y)) !=0) {
	x[0] = y \| z << (32 - k);
	z >>= k;
	}
	else
	x[0] = y;
	#ifndef Sudden_Underflow
	i =
	#endif
	b->wds = (x[1] = z) !=0 ? 2 : 1;
	}
	else {
	k = lo0bits(&z);
	x[0] = z;
	#ifndef Sudden_Underflow
	i =
	#endif
	b->wds = 1;
	k += 32;
	}
	#else
	if ( (y = d1) !=0) {
	if ( (k = lo0bits(&y)) !=0)
	if (k >= 16) {
	x[0] = y \| z << 32 - k & 0xffff;
	x[1] = z >> k - 16 & 0xffff;
	x[2] = z >> k;
	i = 2;
	}
	else {
	x[0] = y & 0xffff;
	x[1] = y >> 16 \| z << 16 - k & 0xffff;
	x[2] = z >> k & 0xffff;
	x[3] = z >> k+16;
	i = 3;
	}
	else {
	x[0] = y & 0xffff;
	x[1] = y >> 16;
	x[2] = z & 0xffff;
	x[3] = z >> 16;
	i = 3;
	}
	}
	else {
	#ifdef DEBUG
	if (!z)
	Bug("Zero passed to d2b");
	#endif
	k = lo0bits(&z);
	if (k >= 16) {
	x[0] = z;
	i = 0;
	}
	else {
	x[0] = z & 0xffff;
	x[1] = z >> 16;
	i = 1;
	}
	k += 32;
	}
	while(!x[i])
	--i;
	b->wds = i + 1;
	#endif
	#ifndef Sudden_Underflow
	if (de) {
	#endif
	*e = de - Bias - (P-1) + k;
	*bits = P - k;
	#ifndef Sudden_Underflow
	}
	else {
	*e = de - Bias - (P-1) + 1 + k;
	#ifdef Pack_32
	bits = 32i - hi0bits(x[i-1]);
	#else
	bits = (i+2)16 - hi0bits(x[i]);
	#endif
	}
	#endif
	return b;
	#undef d0
	#undef d1
	}

	const double
	bigtens[] = { 1e16, 1e32, 1e64, 1e128, 1e256 };
	const double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, 1e-256 };

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

	char strcp_D2A (char a, const char *b)
	{
	while((a = b++))
	a++;
	return a;
	}

	#ifdef NO_STRING_H
	void memcpy_D2A (void a1, void *b1, size_t len)
	{
	char a = (char)a1, *ae = a + len;
	char b = (char)b1, *a0 = a;
	while(a < ae)
	a++ = b++;
	return a0;
	}
	#endif /* NO_STRING_H */