/*	Copyright (c) 1990 UNIX System Laboratories, Inc.	*/
/*	Copyright (c) 1984, 1986, 1987, 1988, 1989, 1990 AT&T	*/
/*	  All Rights Reserved  	*/

/*	THIS IS UNPUBLISHED PROPRIETARY SOURCE CODE OF     	*/
/*	UNIX System Laboratories, Inc.                     	*/
/*	The copyright notice above does not evidence any   	*/
/*	actual or intended publication of such source code.	*/


#ident	"@(#)ucblibc:port/gen/double_decim.c	1.1.3.1"

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

		PROPRIETARY NOTICE (Combined)

This source code is unpublished proprietary information
constituting, or derived under license from AT&T's UNIX(r) System V.
In addition, portions of such source code were derived from Berkeley
4.3 BSD under license from the Regents of the University of
California.



		Copyright Notice 

Notice of copyright on this source code product does not indicate 
publication.

	(c) 1986,1987,1988,1989  Sun Microsystems, Inc
	(c) 1983,1984,1985,1986,1987,1988,1989  AT&T.
	          All rights reserved.
********************************************************************/ 

/* Conversion between binary and decimal floating point. */

#include "base_conv.h"

/* PRIVATE FUNCTIONS */

void
decimal_round(pm, pd, ps, round, sticky)
	decimal_mode   *pm;
	decimal_record *pd;
	fp_exception_field_type *ps;
	char            round;
	unsigned        sticky;

/*
 * Rounds decimal record *pd according to modes in *pm, recording exceptions
 * for inexact or overflow in *ps.  round is the round digit and sticky is 0
 * or non-zero to indicate exact or inexact. pd->ndigits is expected to be
 * correctly set.
 */

{
	int             lsd, i;

	if ((round == '0') && (sticky == 0)) {	/* Exact. */
		goto done;
	}
	*ps |= 1 << fp_inexact;

	switch (pm->rd) {
	case fp_nearest:
		if (round < '5')
			goto done;
		if (round > '5')
			goto roundup;
		if (sticky != 0)
			goto roundup;
		/* Now in ambiguous case; round up if lsd is odd. */
		if (pd->ndigits <= 0)
			goto done;	/* Presumed 0. */
		lsd = pd->ds[pd->ndigits - 1] - '0';
		if ((lsd % 2) == 0)
			goto done;
		goto roundup;
	case fp_positive:
		if (pd->sign != 0)
			goto done;
		goto roundup;
	case fp_negative:
		if (pd->sign == 0)
			goto done;
		goto roundup;
	case fp_tozero:
		goto done;
	}
roundup:
	for (i = (pd->ndigits - 1); (pd->ds[i] == '9') && (i >= 0); i--)
		pd->ds[i] = '0';
	if (i >= 0)
		(pd->ds[i])++;
	else {			/* Rounding carry out has occurred. */
		pd->ds[0] = '1';
		if (pm->df == floating_form) {	/* For E format, simply
						 * adjust exponent. */
			pd->exponent++;
		} else {	/* For F format, increase length of string. */
			if (pd->ndigits > 0)
				pd->ds[pd->ndigits] = '0';
			pd->ndigits++;
		}
	}
	goto ret;
done:
	if (pd->ndigits <= 0) {	/* Create zero string. */
		pd->ds[0] = '0';
		pd->ndigits = 1;
	}
ret:
	pd->ds[pd->ndigits] = 0;/* Terminate string. */
	return;
}

void
binary_to_decimal_integer(pu, nsig, ds, nzeros, ndigs)
	unpacked       *pu;	/* Input unpacked integer value input. */
	unsigned        nsig;	/* Input number of significant digits
				 * required. */
	char            ds[];	/* Output decimal integer string output -
				 * must be large enough. */
unsigned       *nzeros;		/* Output number of implicit trailing zeros
				 * produced. */
unsigned       *ndigs;		/* Output number of explicit digits produced
				 * in ds. */

/*
 * Converts an unpacked integer value *pu into a decimal string in *ds, of
 * length returned in *ndigs. /* Inexactness is indicated by setting
 * ds[ndigs-1] odd.
 */

{

	_big_float     *pd, b, d;
	int             e, i, is;
	_BIG_FLOAT_DIGIT stickyshift;
	char            s[4];

	b.bsize = _BIG_FLOAT_SIZE;	/* Initialize sizes of big floats. */
	d.bsize = _BIG_FLOAT_SIZE;
	_unpacked_to_big_float(pu, &b, &e);
	if (e < 0) {
		_right_shift_base_two(&b, (short unsigned) -e, &stickyshift);
#ifdef DEBUG
		assert(stickyshift == 0);
#endif
	}
	_big_binary_to_big_decimal(&b, &d);
	if (e <= 0)
		pd = &d;
	else {
		int	ipd = (int)pd;	/* kludge: some compiler does not allow
				 	 * pd in switch */
		_big_float_times_power(&d, 2, e, (int) nsig, &pd);
		switch (ipd) {
		case BIG_FLOAT_TIMES_TOOBIG:
			{
				char            bcastring[80];

				(void) sprintf(bcastring, " binary exponent %d ", e);
				_base_conversion_abort(ERANGE, bcastring);
				break;
			}
		case BIG_FLOAT_TIMES_NOMEM:
			{
				char            bcastring[80];

				(void) sprintf(bcastring, " binary exponent %d ", e);
				_base_conversion_abort(ENOMEM, bcastring);
				break;
			}
		default:
#ifdef DEBUG
			if (pd != &d)
				(void) printf(" large binary exponent %d needs heap buffer \n", e);
			printf(" product ");
			_display_big_float(pd, 10);
#endif
			break;
		}
	}
	_fourdigitsquick((short unsigned) pd->bsignificand[pd->blength - 1], s);
	for (i = 0; s[i] == '0'; i++);	/* Find first non-zero digit. */
	for (is = 0; i <= 3;)
		ds[is++] = s[i++];	/* Copy subsequent digits. */

	for (i = (pd->blength - 2); i >= 0; i--) {	/* Convert powers of
							 * 10**4 to decimal
							 * digits. */
		_fourdigitsquick((short unsigned) pd->bsignificand[i], &(ds[is]));
		is += 4;
	}

	ds[is] = 0;
	*ndigs = is;
	*nzeros = pd->bexponent;
	if (pd != &d)
		_free_big_float(pd);

#ifdef DEBUG
	printf(" binary to decimal integer result %s * 10**%d \n", ds, *nzeros);
#endif
}

void
binary_to_decimal_fraction(pu, nsig, nfrac, ds, nzeros, ndigs)
	unpacked       *pu;	/* Input unpacked fraction value output < 1
				 * in magnitude. */
	unsigned        nsig;	/* Input number of significant digits
				 * required. */
	unsigned        nfrac;	/* Input number of digits after point
				 * required. */
	char            ds[];	/* Output decimal integer string output -
				 * must be large enough. */
int            *nzeros;		/* Output number of implicit leading zeros
				 * produced. */
int            *ndigs;		/* Output number of explicit digits produced
				 * in ds. */

/*
 * Converts an unpacked fraction value *pu into a decimal string consisting
 * of a) an implicit '.' b) *nzeros implicit leading zeros c) *ndigs explicit
 * digits in ds ds contains at least nsig significant digits. nzeros + *
 * *ndigs is at least nfrac digits after the point. Inexactness is indicated
 * by sticking to the lsb.
 */

{
	_big_float     *pb, b, d;
	int		ipb = (int)pb;		/* kludge: some compiler does no
						 * allow 'pb' to be used in switch */

	int             e, i, j, is, excess;
	char            s[4];
	int             tensig, tenpower;
	_BIG_FLOAT_DIGIT stickyshift;

	*nzeros = 0;
	if (pu->fpclass == fp_zero) {	/* Exact zero. */
		for (i = 0; i <= nfrac; i++)
			ds[i] = '0';
		for (; i <= nsig; i++)
			ds[i] = '0';
		*ndigs = i;
		return;
	}
	b.bsize = _BIG_FLOAT_SIZE;	/* Initialize sizes of big floats. */
	d.bsize = _BIG_FLOAT_SIZE;
	_unpacked_to_big_float(pu, &b, &e);
	/*
	 * e < 0 always
	 */
	b.bexponent = e;
	tenpower = nsig + (int) (((17 - e - 16 * b.blength) * (unsigned long) 19729) >> 16);
	if (tenpower < nfrac)
		tenpower = nfrac;
	tensig = nfrac;
	if (nsig > tensig)
		tensig = nsig;
	tensig = 1 + (((tensig + 2) * 217706) >> 16);
	tensig = -tensig;

#ifdef DEBUG
	printf(" binary to decimal fraction exponent 2**%d \n", e);
	printf(" binary to decimal fraction nsig %d nfrac %d tenpower %d tensig %d \n", nsig, nfrac, tenpower, tensig);
#endif
	_big_float_times_power(&b, 10, tenpower, tensig, &pb);
	switch (ipb) {
	case BIG_FLOAT_TIMES_TOOBIG:
		{
			char            bcastring[80];

			(void) sprintf(bcastring, " decimal exponent %d ", tenpower);
			_base_conversion_abort(ERANGE, bcastring);
			break;
		}
	case BIG_FLOAT_TIMES_NOMEM:
		{
			char            bcastring[80];

			(void) sprintf(bcastring, " decimal exponent %d ", tenpower);
			_base_conversion_abort(ENOMEM, bcastring);
			break;
		}
	default:
#ifdef DEBUG
		if (pb != &b)
			printf(" large decimal exponent %d needs heap buffer \n", tenpower);
		printf(" product ");
		_display_big_float(pb, 2);
#endif
		break;
	}

	if (pb->bexponent <= -16) {
		/* Have computed appropriate decimal part; now toss fraction. */
		excess = (-pb->bexponent) / 16;
#ifdef DEBUG
		printf(" discard %d excess fraction bits \n", 16 * excess);
#endif
		for (i = 0; (i < excess) && (pb->bsignificand[i] == 0); i++);
		if (i < excess)
			pb->bsignificand[excess] |= 1;	/* Sticky bit for
							 * discarded fraction. */
		for (i = excess; i < pb->blength; i++)
			pb->bsignificand[i - excess] = pb->bsignificand[i];

		pb->blength -= excess;
		pb->bexponent += 16 * excess;
	}
	if (pb->bexponent < 0) {
		_right_shift_base_two(pb, (short unsigned) -pb->bexponent, &stickyshift);
		if (stickyshift != 0)
			pb->bsignificand[0] |= 1;	/* Stick to lsb. */
	}
	_big_binary_to_big_decimal(pb, &d);

	i = d.blength - 1;
	while (d.bsignificand[i] == 0)
		i--;
	_fourdigitsquick((short unsigned) d.bsignificand[i], s);
	for (j = 0; s[j] == '0'; j++);	/* Find first non-zero digit. */
	for (is = 0; j <= 3;)
		ds[is++] = s[j++];	/* Copy subsequent digits. */

	for (i--; i >= 0; i--) {/* Convert powers of 10**4 to decimal digits. */
		_fourdigitsquick((short unsigned) d.bsignificand[i], &(ds[is]));
		is += 4;
	}

	ds[is] = 0;
	*ndigs = is;
#ifdef DEBUG
	assert(tenpower >= is);
#endif
	*nzeros = tenpower - is;/* There were supposed to be tenpower leading
				 * digits, and is were found. */

	if (pb != &b)
		_free_big_float(pb);

#ifdef DEBUG
	printf(" binary to decimal fraction result .%s * 10**%d \n", ds, -(*nzeros));
#endif

}

void
_unpacked_to_decimal(px, pm, pd, ps)
	unpacked       *px;
	decimal_mode   *pm;
	decimal_record *pd;
	fp_exception_field_type *ps;
{
	unpacked        fx, ix;
	unsigned        fmask, imask;
	int             i, intdigs, fracdigs, fraczeros, fracsigs, ids, idsbound, lzbound;
	unsigned        nsig, nfrac, intzeros, intsigs;
	char            is[_INTEGER_SIZE], fs[DECIMAL_STRING_LENGTH];
	char            round = '0';
	unsigned        sticky = 0;

	pd->sign = px->sign;
	pd->fpclass = px->fpclass;
	if ((px->fpclass != fp_normal) && (px->fpclass != fp_subnormal))
		return;
	if ((pm->ndigits >= DECIMAL_STRING_LENGTH) ||
	    ((pm->df == floating_form) && (pm->ndigits < 1))) {	/* Gross overflow or bad
								 * spec. */
overflow:
		*ps |= 1 << fp_overflow;
		return;
	}
	/* Divide x up into integer part ix and fraction part fx.	 */

	ix = *px;
	fx = ix;
	if (ix.exponent <= -1) {/* All fraction. */
		ix.fpclass = fp_zero;
	} else if (ix.exponent >= 159) {	/* All integer. */
		fx.fpclass = fp_zero;
	} else if ((ix.exponent % 32) == 31) {	/* Integer/fraction boundary
						 * is conveniently on a word
						 * boundary. */
		imask = (ix.exponent + 1) / 32;	/* Words 0..imask-1 are
						 * integer; imask..SIZE are
						 * fraction. */
		for (i = 0; i < imask; i++)
			fx.significand[i] = 0;
		for (; i < UNPACKED_SIZE; i++)
			ix.significand[i] = 0;
		_fp_normalize(&fx);
	} else {		/* Integer/fraction boundary falls in the
				 * middle of a word. */
		imask = (ix.exponent + 1) / 32;	/* Words 0..imask-1 are
						 * integer; imask is integer
						 * and fraction ;
						 * imask+1..SIZE are
						 * fraction. */
		for (i = 0; i < imask; i++)
			fx.significand[i] = 0;
		fmask = (1 << (31 - (ix.exponent % 32))) - 1;
		fx.significand[imask] &= fmask;
		ix.significand[imask] &= ~fmask;
		for (i = (imask + 1); i < UNPACKED_SIZE; i++)
			ix.significand[i] = 0;
		_fp_normalize(&fx);
	}
	if (ix.fpclass != fp_zero) {	/* Compute integer part of result. */
		if (pm->df == floating_form)
			nsig = pm->ndigits + 1;	/* Significant digits wanted
						 * for E format, plus one for
						 * rounding. */
		else
			nsig = _INTEGER_SIZE;	/* Significant digits wanted
						 * for F format == all. */

		binary_to_decimal_integer(&ix, nsig, is, &intzeros, &intsigs);
	} else {
		intsigs = 0;
		intzeros = 0;
	}
	intdigs = intsigs + intzeros;
	fracdigs = 0;
	if (((pm->df == fixed_form) && (pm->ndigits >= 0)) ||
	    ((pm->df == floating_form) && ((pm->ndigits + 1) > intdigs))) {	/* Need to compute
										 * fraction part. */
		if (pm->df == floating_form) {	/* Need more significant
						 * digits. */
			nsig = pm->ndigits + 2 - intdigs;	/* Add two for rounding,
								 * sticky. */
			if (nsig > DECIMAL_STRING_LENGTH)
				nsig = DECIMAL_STRING_LENGTH;
			nfrac = 1;
		} else {	/* Need fraction digits. */
			nsig = 0;
			nfrac = pm->ndigits + 2;	/* Add two for rounding,
							 * sticky. */
			if (nfrac > DECIMAL_STRING_LENGTH)
				nfrac = DECIMAL_STRING_LENGTH;
		}
		binary_to_decimal_fraction(&fx, nsig, nfrac, fs, &fraczeros, &fracsigs);
		fracdigs = fraczeros + fracsigs;
	}
	if (pm->df == floating_form) {	/* Combine integer and fraction for E
					 * format. */
		idsbound = intsigs;
		if (idsbound > pm->ndigits)
			idsbound = pm->ndigits;
		for (ids = 0; ids < idsbound; ids++)
			pd->ds[ids] = is[ids];
		/* Put integer into output string. */
		idsbound = intsigs + intzeros;
		if (idsbound > pm->ndigits)
			idsbound = pm->ndigits;
		for (; ids < idsbound; ids++)
			pd->ds[ids] = '0';
		if (ids == pm->ndigits) {	/* Integer part had enough
						 * significant digits. */
			pd->ndigits = ids;
			pd->exponent = intdigs - ids;
			if (ids < intdigs) {	/* Gather rounding info. */
				if (ids < intsigs)
					round = is[ids++];
				else
					round = '0';
				for (; (is[ids] == '0') && (ids < intsigs); ids++);
				if (ids < intsigs)
					sticky = 1;
				if (fx.fpclass != fp_zero)
					sticky = 1;
			} else {/* Integer part is exact - round from
				 * fraction. */
				if (fx.fpclass != fp_zero) {
					int             stickystart;
					/* Fraction non-zero. */
					if (fraczeros > 0) {	/* Round digit is zero. */
						round = '0';
						stickystart = 0;	/* Stickies start with
									 * fs[0]. */
					} else {	/* Round digit is fs[0]. */
						round = fs[0];
						stickystart = 1;	/* Stickies start with
									 * fs[1]. */
					}
					if (sticky == 0) {	/* Search for sticky
								 * bits. */
						for (ids = stickystart; (fs[ids] == '0') && (ids < fracdigs); ids++);
						if (ids < fracdigs)
							sticky = 1;
					}
				}
			}
		} else {	/* Need more significant digits from fraction
				 * part. */
			idsbound = pm->ndigits - ids;
			if (ids == 0) {	/* No integer part - find first
					 * significant digit. */
				for (i = 0; fs[i] == '0'; i++);
				idsbound = i + idsbound + fraczeros;
				i += fraczeros;	/* Point i at first
						 * significant digit. */
			} else
				i = 0;
			if (idsbound > fracdigs)
				idsbound = fracdigs;
			pd->exponent = -idsbound;

			if (fraczeros < idsbound)	/* Compute number of
							 * leading zeros
							 * required. */
				lzbound = fraczeros;
			else
				lzbound = idsbound;
			for (; (i < lzbound); i++)
				pd->ds[ids++] = '0';
			for (; (i < idsbound); i++)
				pd->ds[ids++] = fs[i - fraczeros];
			i -= fraczeros;	/* Don't worry about leading zeros
					 * from now on, we're just rounding */
			if (i < fracsigs) {	/* Gather rounding info.  */
				if (i < 0)
					round = '0';
				else
					round = fs[i];
				i++;
				if (sticky == 0) {	/* Find out if remainder
							 * is exact. */
					if (i < 0)
						i = 0;
					for (; (fs[i] == '0') && (i < fracsigs); i++);
					if (i < fracsigs)
						sticky = 1;
				}
			} else {/* Fraction part is exact - add zero digits
				 * if required. */
				for (; ids < pm->ndigits; ids++)
					pd->ds[ids] = '0';
			}
			pd->ndigits = ids;
		}
		decimal_round(pm, pd, ps, round, sticky);
	} else {		/* Combine integer and fraction for F format. */
		if (pm->ndigits >= 0) {	/* Normal F format. */
			if ((intdigs + pm->ndigits) >= DECIMAL_STRING_LENGTH)
				goto overflow;
			for (ids = 0; ids < intsigs; ids++)
				pd->ds[ids] = is[ids];
			for (; ids < intdigs; ids++)
				pd->ds[ids] = '0';
			/* Copy integer digits. */
			idsbound = fracdigs;
			if (idsbound > pm->ndigits)
				idsbound = pm->ndigits;
			if (fraczeros < idsbound)	/* Compute number of
							 * leading zeros
							 * required. */
				lzbound = fraczeros;
			else
				lzbound = idsbound;
			for (i = 0; (i < lzbound); i++)
				pd->ds[ids++] = '0';
			for (; (i < idsbound); i++)
				pd->ds[ids++] = fs[i - fraczeros];	/* Copy fraction digits. */
			for (; i < pm->ndigits; i++)
				pd->ds[ids++] = '0';
			/* Copy trailing zeros if necessary. */
			pd->ndigits = ids;
			pd->exponent = intdigs - ids;
			i -= fraczeros;	/* Don't worry about leading zeros
					 * from now on, we're just rounding */
			if (i < fracsigs) {	/* Gather rounding info.  */
				if (i < 0)
					round = '0';
				else
					round = fs[i];
				i++;
				if (sticky == 0) {	/* Find out if remainder
							 * is exact. */
					if (i < 0)
						i = 0;
					for (; (fs[i] == '0') && (i < fracsigs); i++);
					if (i < fracsigs)
						sticky = 1;
				}
			}
			decimal_round(pm, pd, ps, round, sticky);
		} else {	/* Bizarre F format - round to left of point. */
			int             roundpos = -pm->ndigits;

			if (intdigs >= DECIMAL_STRING_LENGTH)
				goto overflow;
			if (roundpos >= DECIMAL_STRING_LENGTH)
				goto overflow;
			if (intdigs <= roundpos) {	/* Not enough integer
							 * digits. */
				if (intdigs == roundpos) {
					round = is[0];
					i = 1;
				} else {
					round = '0';
					i = 0;
				}
				for (; (is[i] == '0') && (i < intsigs); i++);
				/* Search for sticky bits. */
				if (i < intsigs)
					sticky = 1;
				pd->ndigits = 0;
			} else {/* Some integer digits do not get rounded
				 * away. */
				for (ids = 0; ids < (intsigs - roundpos); ids++)
					pd->ds[ids] = is[ids];
				for (ids = 0; ids < (intdigs - roundpos); ids++)
					pd->ds[ids] = '0';
				/* Copy integer digits. */
				pd->ndigits = ids;
				if (ids < intsigs) {	/* Inexact. */
					round = is[ids++];
					for (; (is[ids] == '0') && (ids < intsigs); ids++);
					/* Search for non-zero digits. */
					if (ids < intsigs)
						sticky = 1;
				}
			}
			if (fx.fpclass != fp_zero)
				sticky = 1;
			decimal_round(pm, pd, ps, round, sticky);
			for (i = pd->ndigits; i < (pd->ndigits + roundpos); i++)
				pd->ds[i] = '0';	/* Blank out rounded
							 * away digits. */
			pd->exponent = 0;
			pd->ndigits = i;
			pd->ds[i] = 0;	/* Terminate string. */
		}
	}
}

void
double_to_decimal(px, pm, pd, ps)
	double         *px;
	decimal_mode   *pm;
	decimal_record *pd;
	fp_exception_field_type *ps;
{
	double_equivalence kluge;
	unpacked        u;

	*ps = 0;		/* Initialize *ps. */
	kluge.x = *px;
	pd->sign = kluge.f.msw.sign;
	pd->fpclass = _class_double(px);
	switch (pd->fpclass) {
	case fp_zero:
		break;
	case fp_infinity:
		break;
	case fp_quiet:
		break;
	case fp_signaling:
		break;
	default:
		_unpack_double(&u, &kluge.x);
		_unpacked_to_decimal(&u, pm, pd, ps);
	}
}

void
quadruple_to_decimal(px, pm, pd, ps)
	quadruple      *px;
	decimal_mode   *pm;
	decimal_record *pd;
	fp_exception_field_type *ps;
{
	quadruple_equivalence kluge;
	unpacked        u;
	int             i;

	*ps = 0;		/* Initialize *ps - no exceptions. */
	for (i = 0; i < 4; i++)
		kluge.x.u[i] = px->u[i];
	pd->sign = kluge.f.msw.sign;
	pd->fpclass = _class_quadruple(px);
	switch (pd->fpclass) {
	case fp_zero:
		break;
	case fp_infinity:
		break;
	case fp_quiet:
		break;
	case fp_signaling:
		break;
	default:
		_unpack_quadruple(&u, px);
		_unpacked_to_decimal(&u, pm, pd, ps);
	}
}
