/* _Ldtob function */
#include <float.h>
#include <string.h>
#include "xmath.h"
#include "xstdio.h"
#include <stdlib.h>

/* macros */
#define NDIG	8

/* static data */
static const ldouble pows[] = {
	1e1L, 1e2L, 1e4L, 1e8L, 1e16L, 1e32L,
#if 0x100 < _LBIAS	/* assume IEEE 754 8- or 10-byte */
	1e64L, 1e128L, 1e256L,
#if _DLONG	/* assume IEEE 754 10-byte */
	1e512L, 1e1024L, 1e2048L, 1e4096L,
#endif
#endif
};

static short _Ldunscale(short *, ldouble *);
static void _Genld(_Pft *, char, char *, short, short);

void _Ldtob(_Pft *px, char code)
{	/* convert long double to text */
	char ac[32];
	char *p = ac;
	ldouble ldval = px->v.ld;
	short errx, nsig, xexp;

	if (px->prec < 0)
		px->prec = 6;
	else if (px->prec == 0 && (code == 'g' || code == 'G'))
		px->prec = 1;
	if (0 < (errx = _Ldunscale(&xexp, &px->v.ld))) {
		/* x == Nan, x == INF */
		memcpy(px->s, errx == NAN ? "NaN" : "Inf", px->n1 = 3);
		return;
	} else if (0 == errx)	/*x == 0 */
		nsig = 0, xexp = 0;
	else {	/* 0 < |x|, convert it */
		{	/* scale ldval to ~~10^(NDIG/2) */
			int i, n;

			if (ldval < 0.0)
				ldval = -ldval;
			if ((xexp = xexp * 30103L / 100000L - NDIG/2) < 0) {
				/* scale up */
				n = (-xexp + (NDIG/2-1)) & ~(NDIG/2-1), xexp = -n;
				for (i = 0; 0 < n; n >>= 1, ++i)
					if (n & 1)
						ldval *= pows[i];
			} else if (0 < xexp) {	/* scale down */
				ldouble factor = 1.0;

				xexp &= ~(NDIG/2-1);
				for (n = xexp, i = 0; 0 < n; n >>= 1, ++i)
					if (n & 1)
						factor *= pows[i];
				ldval /= factor;
			}
		}
		{	/* convert significant digits */
			int gen = px->prec +
				(code == 'f' ? xexp + 2 + NDIG : 2 + NDIG / 2);

			if (LDBL_DIG + NDIG / 2 < gen)
				gen = LDBL_DIG + NDIG / 2;
			for (*p++ = '0'; 0 < gen && 0.0 < ldval; p += NDIG) {
				/* convert NDIG at a time */
				int j;
				long lo = (long) ldval;

				if (0 < (gen -= NDIG))
					ldval = (ldval - (ldouble) lo) * 1e8L;
				for (p += NDIG, j = NDIG; 0 < lo && 0 <= --j; ) {
					/* convert NDIG digits */
					ldiv_t qr;

					qr = ldiv(lo, 10);
					*--p = qr.rem + '0', lo = qr.quot;
				}
				while (0 <= --j)
					*--p = '0';
			}
			gen = p - &ac[1];
			for (p = &ac[1], xexp += NDIG - 1; *p == '0'; ++p)
				--gen, --xexp;	/* correct xexp */
			nsig = px->prec + (code == 'f' ? xexp + 1 :
				code == 'e' || code == 'E' ? 1 : 0);
			if (gen < nsig)
				nsig = gen;
			if (0 < nsig) {	/* round and strip trailing zeros */
				const char drop = nsig < gen && '5' <= p[nsig] ? '9' : '0';
				int n;

				for (n = nsig; p[--n] == drop; )
					--nsig;
				if (drop == '9')
					++p[n];
				if (n < 0)
					--p, ++nsig, ++xexp;
			}
		}
	}
	_Genld(px, code, p, nsig, xexp);
}

#if _DLONG	/* 10-byte IEEE format */
#define _LMASK	0x7fff
#define _LMAX	0x7fff
#define _LSIGN	0x8000
#if _D0==3	/* little-endian order */
#define _L0	4
#define _L1	3
#define _L2	2
#define _L3	1
#define _L4	0
#else
#define _L0	0
#define _L1	1
#define _L2	2
#define _L3	3
#define _L4	4
#endif

static short dnorm(unsigned short *ps)
{	/* normalize long double fraction */
	short xchar;

	for (xchar = 0; ps[_L1] == 0; xchar -= 16) {	/* shift left by 16 */
		ps[_L1] = ps[_L2], ps[_L2] = ps[_L3];
		ps[_L3] = ps[_L4], ps[_L4] = 0;
	}
	for (; ps[_L1] < 1U<<_LOFF; --xchar) {	/* shift left by 1 */
		ps[_L1] = ps[_L1] << 1 | ps[_L2] >> 15;
		ps[_L2] = ps[_L2] << 1 | ps[_L3] >> 15;
		ps[_L3] = ps[_L3] << 1 | ps[_L4] >> 15;
		ps[_L4] <<= 1;
	}
	return (xchar);
}

static short _Ldunscale(short *pex, ldouble *px)
{	/* separate *px to |frac| < 1/2 and 2^*pex */
	unsigned short *ps = (unsigned short *) px;
	short xchar = ps[_L0] & _LMAX;

	if (xchar == _LMAX) {	/* NaN or INF */
		*pex = 0;
		return (ps[_L1] & 0x7fff || ps[_L2] || ps[_L3] || ps[_L4] ?
			NAN : INF);
	} else if (ps[_L1] == 0 && ps[_L2] == 0 && ps[_L3] == 0 && ps[_L4] == 0) {
		/* zero */
		*pex = 0;
		return (0);
	} else {	/* finite, reduce to [1/2, 1) */
		xchar += dnorm(ps);
		ps[_L0] = ps[_L0] & _LSIGN | _LBIAS;
		*pex = xchar - _LBIAS;
		return (FINITE);
	}
}

#else	/* long double same as double */
/*
static short _Dnorm(unsigned short *ps)
{
	short xchar;
	unsigned short sign = ps[_D0] & _DSIGN;

	xchar = 0;
	if ((ps[_D0] &= _DFRAC) != 0 || ps[_D1] || ps[_D2] || ps[_D3]) {
		for (; ps[_D0] == 0; xchar -= 16) {
			ps[_D0] = ps[_D1], ps[_D1] = ps[_D2];
			ps[_D2] = ps[_D3], ps[_D3] = 0;
		}
		for (; ps[_D0] < 1<<_DOFF; --xchar) {
			ps[_D0] = ps[_D0] << 1 | ps[_D1] >> 15;
			ps[_D1] = ps[_D1] << 1 | ps[_D2] >> 15;
			ps[_D2] = ps[_D2] << 1 | ps[_D3] >> 15;
			ps[_D3] <<= 1;
		}
		for (; 1<<_DOFF+1 <= ps[_D0]; ++xchar) {
			ps[_D3] = ps[_D3] >> 1 | ps[_D2] << 15;
			ps[_D2] = ps[_D2] >> 1 | ps[_D1] << 15;
			ps[_D1] = ps[_D1] >> 1 | ps[_D0] << 15;
			ps[_D0] >>= 1;
		}
		ps[_D0] &= _DFRAC;
	}
	ps[_D0] |= sign;
	return (xchar);
}
*/
static short _Ldunscale(short *pex, ldouble *px)
{	/* separate *px to |frac| < 1/2 and 2^*pex */
	unsigned short *ps = (unsigned short *) px;
	short xchar = (ps[_D0] & _DMASK) >> _DOFF;

	if (xchar == _DMAX) {	/* NaN or INF */
		*pex = 0;
		return (ps[_D0] & _DFRAC || ps[_D1] || ps[_D2] || ps[_D3] ?
			NAN : INF);
	} else if (0 < xchar /* || (xchar = _Dnorm(ps)) != 0 */) {
		/* finite, reduce to [1/2, 1) */
		ps[_D0] = ps[_D0] & ~_DMASK | _DBIAS << _DOFF;
#if _LONG_DOUBLE
		*pex = xchar - _DBIAS;
#else
		*pex = xchar - _DBIAS + 1;	/* for SGI */
#endif
		return (FINITE);
	} else if (xchar < 0) {	/* error! */
		return (NAN);
	} else {	/* zero */
		*pex = 0;
		return (0);
	}
}

#endif

static void _Genld(_Pft *px, char code, char *p, short nsig, short xexp)
{	/* generate long double text */
	const char point = '.';

	if (nsig <= 0)
		nsig = 1, p = "0";
	if (code == 'f' || (code == 'g' || code == 'G') &&
		-4 <= xexp && xexp < px->prec) {	/* 'f' format */
		++xexp;		/* change to leading digit count */
		if (code != 'f') {	/* fixup for 'g' */
			if (!(px->flags & _FNO) && nsig < px->prec)
				px->prec = nsig;
			if ((px->prec -= xexp) < 0)
				px->prec = 0;
		}
		if (xexp <= 0) {	/* digits only to right of point */
			px->s[px->n1++] = '0';
			if (0 < px->prec || px->flags & _FNO)
				px->s[px->n1++] = point;
			if (px->prec < -xexp)
				xexp = -px->prec;
			px->nz1 = -xexp;
			px->prec += xexp;
			if (px->prec < nsig)
				nsig = px->prec;
			memcpy(&px->s[px->n1], p, px->n2 = nsig);
			px->nz2 = px->prec - nsig;
		} else if (nsig < xexp) {	/* zeros before point */
			memcpy(&px->s[px->n1], p, nsig);
			px->n1 += nsig;
			px->nz1 = xexp - nsig;
			if (0 < px->prec || px->flags & _FNO)
				px->s[px->n1] = point, ++px->n2;
			px->nz2 = px->prec;
		} else {	/* enough digits before point */
			memcpy(&px->s[px->n1], p, xexp);
			px->n1 += xexp;
			nsig -= xexp;
			if (0 < px->prec || px->flags & _FNO)
				px->s[px->n1++] = point;
			if (px->prec < nsig)
				nsig = px->prec;
			memcpy(&px->s[px->n1], p + xexp, nsig);
			px->n1 += nsig;
			px->nz1 = px->prec - nsig;
		}
	} else {	/* 'e' format */
		if (code == 'g' || code == 'G') {	/* fixup for 'g' */
			if (nsig < px->prec)
				px->prec = nsig;
			if (--px->prec < 0)
				px->prec = 0;
			code = code == 'g' ? 'e' : 'E';
		}
		px->s[px->n1++] = *p++;
		if (0 < px->prec || px->flags & _FNO)
			px->s[px->n1++] = point;
		if (0 < px->prec) {	/* put fraction digits */
			if (px->prec < --nsig)
				nsig = px->prec;
			memcpy(&px->s[px->n1], p, nsig);
			px->n1 += nsig;
			px->nz1 = px->prec - nsig;
		}
		p = &px->s[px->n1];	/* put exponent */
		*p++ = code;
		if (0 <= xexp)
			*p++ = '+';
		else {	/* negative exponent */
			*p++ = '-';
			xexp = -xexp;
		}
		if (100 <= xexp) {	/* put oversize exponent */
			if (1000 <= xexp)
				*p++ = xexp / 1000 + '0', xexp %= 1000;
			*p++ = xexp / 100 + '0', xexp %= 100;
		}
		*p++ = xexp / 10 + '0', xexp %= 10;
		*p++ = xexp + '0';
		px->n2 = p - &px->s[px->n1];
	}
	if ((px->flags & (_FMI | _FZE)) == _FZE) {	/* pad with leading zeros */
		int n = px->n0 + px->n1 + px->nz1 + px->n2 + px->nz2;

		if (n < px->width)
			px->nz0 = px->width - n;
	}
}