ldtoa.c

KPIT GNU Tools is a set of GNU development tools for Renesas microcontrollers.

/* Move internal format number from a to b.
 */
static void emovz(register short unsigned int *a, register short unsigned int *b)
{
register int i;

for( i=0; i<NI-1; i++ )
	*b++ = *a++;
/* clear low guard word */
*b = 0;
}

/* Return nonzero if internal format number is a NaN.
 */

static int eiisnan (short unsigned int *x)
{
int i;

if( (x[E] & 0x7fff) == 0x7fff )
	{
	for( i=M+1; i<NI; i++ )
		{
		if( x[i] != 0 )
			return(1);
		}
	}
return(0);
}

#if LDBL_MANT_DIG == 64

/* Return nonzero if internal format number is infinite. */
static int 
eiisinf (x)
     unsigned short x[];
{

#ifdef NANS
  if (eiisnan (x))
    return (0);
#endif
  if ((x[E] & 0x7fff) == 0x7fff)
    return (1);
  return (0);
}
#endif /* LDBL_MANT_DIG == 64 */

/*
;	Compare significands of numbers in internal format.
;	Guard words are included in the comparison.
;
;	unsigned short a[NI], b[NI];
;	cmpm( a, b );
;
;	for the significands:
;	returns	+1 if a > b
;		 0 if a == b
;		-1 if a < b
*/
static int ecmpm(register short unsigned int *a, register short unsigned int *b)
{
int i;

a += M; /* skip up to significand area */
b += M;
for( i=M; i<NI; i++ )
	{
	if( *a++ != *b++ )
		goto difrnt;
	}
return(0);

difrnt:
if( *(--a) > *(--b) )
	return(1);
else
	return(-1);
}


/*
;	Shift significand down by 1 bit
*/

static void eshdn1(register short unsigned int *x)
{
register unsigned short bits;
int i;

x += M;	/* point to significand area */

bits = 0;
for( i=M; i<NI; i++ )
	{
	if( *x & 1 )
		bits |= 1;
	*x >>= 1;
	if( bits & 2 )
		*x |= 0x8000;
	bits <<= 1;
	++x;
	}	
}



/*
;	Shift significand up by 1 bit
*/

static void eshup1(register short unsigned int *x)
{
register unsigned short bits;
int i;

x += NI-1;
bits = 0;

for( i=M; i<NI; i++ )
	{
	if( *x & 0x8000 )
		bits |= 1;
	*x <<= 1;
	if( bits & 2 )
		*x |= 1;
	bits <<= 1;
	--x;
	}
}



/*
;	Shift significand down by 8 bits
*/

static void eshdn8(register short unsigned int *x)
{
register unsigned short newbyt, oldbyt;
int i;

x += M;
oldbyt = 0;
for( i=M; i<NI; i++ )
	{
	newbyt = *x << 8;
	*x >>= 8;
	*x |= oldbyt;
	oldbyt = newbyt;
	++x;
	}
}

/*
;	Shift significand up by 8 bits
*/

static void eshup8(register short unsigned int *x)
{
int i;
register unsigned short newbyt, oldbyt;

x += NI-1;
oldbyt = 0;

for( i=M; i<NI; i++ )
	{
	newbyt = *x >> 8;
	*x <<= 8;
	*x |= oldbyt;
	oldbyt = newbyt;
	--x;
	}
}

/*
;	Shift significand up by 16 bits
*/

static void eshup6(register short unsigned int *x)
{
int i;
register unsigned short *p;

p = x + M;
x += M + 1;

for( i=M; i<NI-1; i++ )
	*p++ = *x++;

*p = 0;
}

/*
;	Shift significand down by 16 bits
*/

static void eshdn6(register short unsigned int *x)
{
int i;
register unsigned short *p;

x += NI-1;
p = x + 1;

for( i=M; i<NI-1; i++ )
	*(--p) = *(--x);

*(--p) = 0;
}

/*
;	Add significands
;	x + y replaces y
*/

static void eaddm(short unsigned int *x, short unsigned int *y)
{
register unsigned long a;
int i;
unsigned int carry;

x += NI-1;
y += NI-1;
carry = 0;
for( i=M; i<NI; i++ )
	{
	a = (unsigned long )(*x) + (unsigned long )(*y) + carry;
	if( a & 0x10000 )
		carry = 1;
	else
		carry = 0;
	*y = (unsigned short )a;
	--x;
	--y;
	}
}

/*
;	Subtract significands
;	y - x replaces y
*/

static void esubm(short unsigned int *x, short unsigned int *y)
{
unsigned long a;
int i;
unsigned int carry;

x += NI-1;
y += NI-1;
carry = 0;
for( i=M; i<NI; i++ )
	{
	a = (unsigned long )(*y) - (unsigned long )(*x) - carry;
	if( a & 0x10000 )
		carry = 1;
	else
		carry = 0;
	*y = (unsigned short )a;
	--x;
	--y;
	}
}


/* Divide significands */


/* Multiply significand of e-type number b
by 16-bit quantity a, e-type result to c. */

static void m16m(short unsigned int a, short unsigned int *b, short unsigned int *c)
{
register unsigned short *pp;
register unsigned long carry;
unsigned short *ps;
unsigned short p[NI];
unsigned long aa, m;
int i;

aa = a;
pp = &p[NI-2];
*pp++ = 0;
*pp = 0;
ps = &b[NI-1];

for( i=M+1; i<NI; i++ )
	{
	if( *ps == 0 )
		{
		--ps;
		--pp;
		*(pp-1) = 0;
		}
	else
		{
		m = (unsigned long) aa * *ps--;
		carry = (m & 0xffff) + *pp;
		*pp-- = (unsigned short )carry;
		carry = (carry >> 16) + (m >> 16) + *pp;
		*pp = (unsigned short )carry;
		*(pp-1) = carry >> 16;
		}
	}
for( i=M; i<NI; i++ )
	c[i] = p[i];
}


/* Divide significands. Neither the numerator nor the denominator
is permitted to have its high guard word nonzero.  */


static int edivm(short unsigned int *den, short unsigned int *num, LDPARMS *ldp)
{
int i;
register unsigned short *p;
unsigned long tnum;
unsigned short j, tdenm, tquot;
unsigned short tprod[NI+1];
unsigned short *equot = ldp->equot;

p = &equot[0];
*p++ = num[0];
*p++ = num[1];

for( i=M; i<NI; i++ )
	{
	*p++ = 0;
	}
eshdn1( num );
tdenm = den[M+1];
for( i=M; i<NI; i++ )
	{
	/* Find trial quotient digit (the radix is 65536). */
	tnum = (((unsigned long) num[M]) << 16) + num[M+1];

	/* Do not execute the divide instruction if it will overflow. */
        if( (tdenm * 0xffffUL) < tnum )
		tquot = 0xffff;
	else
		tquot = tnum / tdenm;

		/* Prove that the divide worked. */
/*
	tcheck = (unsigned long )tquot * tdenm;
	if( tnum - tcheck > tdenm )
		tquot = 0xffff;
*/
	/* Multiply denominator by trial quotient digit. */
	m16m( tquot, den, tprod );
	/* The quotient digit may have been overestimated. */
	if( ecmpm( tprod, num ) > 0 )
		{
		tquot -= 1;
		esubm( den, tprod );
		if( ecmpm( tprod, num ) > 0 )
			{
			tquot -= 1;
			esubm( den, tprod );
			}
		}
/*
	if( ecmpm( tprod, num ) > 0 )
		{
		eshow( "tprod", tprod );
		eshow( "num  ", num );
		printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
			 tnum, den[M+1], tquot );
		}
*/
	esubm( tprod, num );
/*
	if( ecmpm( num, den ) >= 0 )
		{
		eshow( "num  ", num );
		eshow( "den  ", den );
		printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
			 tnum, den[M+1], tquot );
		}
*/
	equot[i] = tquot;
	eshup6(num);
	}
/* test for nonzero remainder after roundoff bit */
p = &num[M];
j = 0;
for( i=M; i<NI; i++ )
	{
	j |= *p++;
	}
if( j )
	j = 1;

for( i=0; i<NI; i++ )
	num[i] = equot[i];

return( (int )j );
}



/* Multiply significands */
static int emulm(short unsigned int *a, short unsigned int *b, LDPARMS *ldp) 
{
unsigned short *p, *q;
unsigned short pprod[NI];
unsigned short j;
int i;
unsigned short *equot = ldp->equot;

equot[0] = b[0];
equot[1] = b[1];
for( i=M; i<NI; i++ )
	equot[i] = 0;

j = 0;
p = &a[NI-1];
q = &equot[NI-1];
for( i=M+1; i<NI; i++ )
	{
	if( *p == 0 )
		{
		--p;
		}
	else
		{
		m16m( *p--, b, pprod );
		eaddm(pprod, equot);
		}
	j |= *q;
	eshdn6(equot);
	}

for( i=0; i<NI; i++ )
	b[i] = equot[i];

/* return flag for lost nonzero bits */
return( (int)j );
}


/*
static void eshow(str, x)
char *str;
unsigned short *x;
{
int i;

printf( "%s ", str );
for( i=0; i<NI; i++ )
	printf( "%04x ", *x++ );
printf( "\n" );
}
*/


/*
 * Normalize and round off.
 *
 * The internal format number to be rounded is "s".
 * Input "lost" indicates whether the number is exact.
 * This is the so-called sticky bit.
 *
 * Input "subflg" indicates whether the number was obtained
 * by a subtraction operation.  In that case if lost is nonzero
 * then the number is slightly smaller than indicated.
 *
 * Input "exp" is the biased exponent, which may be negative.
 * the exponent field of "s" is ignored but is replaced by
 * "exp" as adjusted by normalization and rounding.
 *
 * Input "rcntrl" is the rounding control.
 */


static void emdnorm(short unsigned int *s, int lost, int subflg, long int exp, int rcntrl, LDPARMS *ldp)
{
int i, j;
unsigned short r;

/* Normalize */
j = enormlz( s );

/* a blank significand could mean either zero or infinity. */
#ifndef INFINITY
if( j > NBITS )
	{
	ecleazs( s );
	return;
	}
#endif
exp -= j;
#ifndef INFINITY
if( exp >= 32767L )
	goto overf;
#else
if( (j > NBITS) && (exp < 32767L) )
	{
	ecleazs( s );
	return;
	}
#endif
if( exp < 0L )
	{
	if( exp > (long )(-NBITS-1) )
		{
		j = (int )exp;
		i = eshift( s, j );
		if( i )
			lost = 1;
		}
	else
		{
		ecleazs( s );
		return;
		}
	}
/* Round off, unless told not to by rcntrl. */
if( rcntrl == 0 )
	goto mdfin;
/* Set up rounding parameters if the control register changed. */
if( ldp->rndprc != ldp->rlast )
	{
	ecleaz( ldp->rbit );
	switch( ldp->rndprc )
		{
		default:
		case NBITS:
			ldp->rw = NI-1; /* low guard word */
			ldp->rmsk = 0xffff;
			ldp->rmbit = 0x8000;
			ldp->rebit = 1;
			ldp->re = ldp->rw - 1;
			break;
		case 113:
			ldp->rw = 10;
			ldp->rmsk = 0x7fff;
			ldp->rmbit = 0x4000;
			ldp->rebit = 0x8000;
			ldp->re = ldp->rw;
			break;
		case 64:
			ldp->rw = 7;
			ldp->rmsk = 0xffff;
			ldp->rmbit = 0x8000;
			ldp->rebit = 1;
			ldp->re = ldp->rw-1;
			break;
/* For DEC arithmetic */
		case 56:
			ldp->rw = 6;
			ldp->rmsk = 0xff;
			ldp->rmbit = 0x80;
			ldp->rebit = 0x100;
			ldp->re = ldp->rw;
			break;
		case 53:
			ldp->rw = 6;
			ldp->rmsk = 0x7ff;
			ldp->rmbit = 0x0400;
			ldp->rebit = 0x800;
			ldp->re = ldp->rw;
			break;
		case 24:
			ldp->rw = 4;
			ldp->rmsk = 0xff;
			ldp->rmbit = 0x80;
			ldp->rebit = 0x100;
			ldp->re = ldp->rw;
			break;
		}
	ldp->rbit[ldp->re] = ldp->rebit;
	ldp->rlast = ldp->rndprc;
	}

/* Shift down 1 temporarily if the data structure has an implied
 * most significant bit and the number is denormal.
 * For rndprc = 64 or NBITS, there is no implied bit.
 * But Intel long double denormals lose one bit of significance even so.
 */
#if IBMPC
if( (exp <= 0) && (ldp->rndprc != NBITS) )
#else
if( (exp <= 0) && (ldp->rndprc != 64) && (ldp->rndprc != NBITS) )
#endif
	{
	lost |= s[NI-1] & 1;
	eshdn1(s);
	}
/* Clear out all bits below the rounding bit,
 * remembering in r if any were nonzero.
 */
r = s[ldp->rw] & ldp->rmsk;
if( ldp->rndprc < NBITS )
	{
	i = ldp->rw + 1;
	while( i < NI )
		{
		if( s[i] )
			r |= 1;
		s[i] = 0;
		++i;
		}
	}
s[ldp->rw] &= ~ldp->rmsk;
if( (r & ldp->rmbit) != 0 )
	{
	if( r == ldp->rmbit )
		{
		if( lost == 0 )
			{ /* round to even */
			if( (s[ldp->re] & ldp->rebit) == 0 )
				goto mddone;
			}
		else
			{
			if( subflg != 0 )
				goto mddone;
			}
		}
	eaddm( ldp->rbit, s );
	}
mddone:
#if IBMPC
if( (exp <= 0) && (ldp->rndprc != NBITS) )
#else
if( (exp <= 0) && (ldp->rndprc != 64) && (ldp->rndprc != NBITS) )
#endif
	{
	eshup1(s);
	}
if( s[2] != 0 )
	{ /* overflow on roundoff */
	eshdn1(s);
	exp += 1;
	}
mdfin:
s[NI-1] = 0;
if( exp >= 32767L )
	{
#ifndef INFINITY
overf:
#endif
#ifdef INFINITY
	s[1] = 32767;
	for( i=2; i<NI-1; i++ )
		s[i] = 0;
#else
	s[1] = 32766;
	s[2] = 0;
	for( i=M+1; i<NI-1; i++ )
		s[i] = 0xffff;
	s[NI-1] = 0;
	if( (ldp->rndprc < 64) || (ldp->rndprc == 113) )
		{
		s[ldp->rw] &= ~ldp->rmsk;
		if( ldp->rndprc == 24 )
			{
			s[5] = 0;
			s[6] = 0;
			}
		}
#endif
	return;
	}
if( exp < 0 )
	s[1] = 0;
else
	s[1] = (unsigned short )exp;
}



/*
;	Subtract external format numbers.
;
;	unsigned short a[NE], b[NE], c[NE];
;       LDPARMS *ldp;
;	esub( a, b, c, ldp );	 c = b - a
*/

static void esub(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp)
{

#ifdef NANS