ldtoa.c

Newlib 嵌入式 C库 标准实现代码

#ifdef NANS
if( eisnan(a) )
	{
	emov (a, c);
	return;
	}
if( eisnan(b) )
	{
	emov(b,c);
	return;
	}
/* Infinity minus infinity is a NaN.
 * Test for subtracting infinities of the same sign.
 */
if( eisinf(a) && eisinf(b) && ((eisneg (a) ^ eisneg (b)) == 0))
	{
	mtherr( "esub", DOMAIN );
	enan( c, NBITS );
	return;
	}
#endif
eadd1( a, b, c, 1, ldp );
}



static void eadd1(short unsigned int *a, short unsigned int *b, short unsigned int *c, int subflg, LDPARMS *ldp)
{
unsigned short ai[NI], bi[NI], ci[NI];
int i, lost, j, k;
long lt, lta, ltb;

#ifdef USE_INFINITY
if( eisinf(a) )
	{
	emov(a,c);
	if( subflg )
		eneg(c);
	return;
	}
if( eisinf(b) )
	{
	emov(b,c);
	return;
	}
#endif
emovi( a, ai );
emovi( b, bi );
if( subflg )
	ai[0] = ~ai[0];

/* compare exponents */
lta = ai[E];
ltb = bi[E];
lt = lta - ltb;
if( lt > 0L )
	{	/* put the larger number in bi */
	emovz( bi, ci );
	emovz( ai, bi );
	emovz( ci, ai );
	ltb = bi[E];
	lt = -lt;
	}
lost = 0;
if( lt != 0L )
	{
	if( lt < (long )(-NBITS-1) )
		goto done;	/* answer same as larger addend */
	k = (int )lt;
	lost = eshift( ai, k ); /* shift the smaller number down */
	}
else
	{
/* exponents were the same, so must compare significands */
	i = ecmpm( ai, bi );
	if( i == 0 )
		{ /* the numbers are identical in magnitude */
		/* if different signs, result is zero */
		if( ai[0] != bi[0] )
			{
			eclear(c);
			return;
			}
		/* if same sign, result is double */
		/* double denomalized tiny number */
		if( (bi[E] == 0) && ((bi[3] & 0x8000) == 0) )
			{
			eshup1( bi );
			goto done;
			}
		/* add 1 to exponent unless both are zero! */
		for( j=1; j<NI-1; j++ )
			{
			if( bi[j] != 0 )
				{
/* This could overflow, but let emovo take care of that. */
				ltb += 1;
				break;
				}
			}
		bi[E] = (unsigned short )ltb;
		goto done;
		}
	if( i > 0 )
		{	/* put the larger number in bi */
		emovz( bi, ci );
		emovz( ai, bi );
		emovz( ci, ai );
		}
	}
if( ai[0] == bi[0] )
	{
	eaddm( ai, bi );
	subflg = 0;
	}
else
	{
	esubm( ai, bi );
	subflg = 1;
	}
emdnorm( bi, lost, subflg, ltb, 64, ldp );

done:
emovo( bi, c, ldp );
}



/*
;	Divide.
;
;	unsigned short a[NE], b[NE], c[NE];
;       LDPARMS *ldp;
;	ediv( a, b, c, ldp );	c = b / a
*/
static void ediv(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp)
{
unsigned short ai[NI], bi[NI];
int i;
long lt, lta, ltb;

#ifdef NANS
/* Return any NaN input. */
if( eisnan(a) )
	{
	emov(a,c);
	return;
	}
if( eisnan(b) )
	{
	emov(b,c);
	return;
	}
/* Zero over zero, or infinity over infinity, is a NaN. */
if( ((ecmp(a,ezero) == 0) && (ecmp(b,ezero) == 0))
	|| (eisinf (a) && eisinf (b)) )
	{
	mtherr( "ediv", DOMAIN );
	enan( c, NBITS );
	return;
	}
#endif
/* Infinity over anything else is infinity. */
#ifdef USE_INFINITY
if( eisinf(b) )
	{
	if( eisneg(a) ^ eisneg(b) )
		*(c+(NE-1)) = 0x8000;
	else
		*(c+(NE-1)) = 0;
	einfin(c, ldp);
	return;
	}
if( eisinf(a) )
	{
	eclear(c);
	return;
	}
#endif
emovi( a, ai );
emovi( b, bi );
lta = ai[E];
ltb = bi[E];
if( bi[E] == 0 )
	{ /* See if numerator is zero. */
	for( i=1; i<NI-1; i++ )
		{
		if( bi[i] != 0 )
			{
			ltb -= enormlz( bi );
			goto dnzro1;
			}
		}
	eclear(c);
	return;
	}
dnzro1:

if( ai[E] == 0 )
	{	/* possible divide by zero */
	for( i=1; i<NI-1; i++ )
		{
		if( ai[i] != 0 )
			{
			lta -= enormlz( ai );
			goto dnzro2;
			}
		}
	if( ai[0] == bi[0] )
		*(c+(NE-1)) = 0;
	else
		*(c+(NE-1)) = 0x8000;
	einfin(c, ldp);
	mtherr( "ediv", SING );
	return;
	}
dnzro2:

i = edivm( ai, bi, ldp );
/* calculate exponent */
lt = ltb - lta + EXONE;
emdnorm( bi, i, 0, lt, 64, ldp );
/* set the sign */
if( ai[0] == bi[0] )
	bi[0] = 0;
else
	bi[0] = 0Xffff;
emovo( bi, c, ldp );
}



/*
;	Multiply.
;
;	unsigned short a[NE], b[NE], c[NE];
;       LDPARMS *ldp
;	emul( a, b, c, ldp );	c = b * a
*/
static void emul(short unsigned int *a, short unsigned int *b, short unsigned int *c, LDPARMS *ldp)
{
unsigned short ai[NI], bi[NI];
int i, j;
long lt, lta, ltb;

#ifdef NANS
/* NaN times anything is the same NaN. */
if( eisnan(a) )
	{
	emov(a,c);
	return;
	}
if( eisnan(b) )
	{
	emov(b,c);
	return;
	}
/* Zero times infinity is a NaN. */
if( (eisinf(a) && (ecmp(b,ezero) == 0))
	|| (eisinf(b) && (ecmp(a,ezero) == 0)) )
	{
	mtherr( "emul", DOMAIN );
	enan( c, NBITS );
	return;
	}
#endif
/* Infinity times anything else is infinity. */
#ifdef USE_INFINITY
if( eisinf(a) || eisinf(b) )
	{
	if( eisneg(a) ^ eisneg(b) )
		*(c+(NE-1)) = 0x8000;
	else
		*(c+(NE-1)) = 0;
	einfin(c, ldp);
	return;
	}
#endif
emovi( a, ai );
emovi( b, bi );
lta = ai[E];
ltb = bi[E];
if( ai[E] == 0 )
	{
	for( i=1; i<NI-1; i++ )
		{
		if( ai[i] != 0 )
			{
			lta -= enormlz( ai );
			goto mnzer1;
			}
		}
	eclear(c);
	return;
	}
mnzer1:

if( bi[E] == 0 )
	{
	for( i=1; i<NI-1; i++ )
		{
		if( bi[i] != 0 )
			{
			ltb -= enormlz( bi );
			goto mnzer2;
			}
		}
	eclear(c);
	return;
	}
mnzer2:

/* Multiply significands */
j = emulm( ai, bi, ldp );
/* calculate exponent */
lt = lta + ltb - (EXONE - 1);
emdnorm( bi, j, 0, lt, 64, ldp );
/* calculate sign of product */
if( ai[0] == bi[0] )
	bi[0] = 0;
else
	bi[0] = 0xffff;
emovo( bi, c, ldp );
}



#if LDBL_MANT_DIG > 64
static void e113toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp)
{
register unsigned short r;
unsigned short *e, *p;
unsigned short yy[NI];
int denorm, i;

e = pe;
denorm = 0;
ecleaz(yy);
#ifdef IBMPC
e += 7;
#endif
r = *e;
yy[0] = 0;
if( r & 0x8000 )
	yy[0] = 0xffff;
r &= 0x7fff;
#ifdef USE_INFINITY
if( r == 0x7fff )
	{
#ifdef NANS
#ifdef IBMPC
	for( i=0; i<7; i++ )
		{
		if( pe[i] != 0 )
			{
			enan( y, NBITS );
			return;
			}
		}
#else  /* !IBMPC */
	for( i=1; i<8; i++ )
		{
		if( pe[i] != 0 )
			{
			enan( y, NBITS );
			return;
			}
		}
#endif /* !IBMPC */
#endif /* NANS */
	eclear( y );
	einfin( y, ldp );
	if( *e & 0x8000 )
		eneg(y);
	return;
	}
#endif  /* INFINITY */
yy[E] = r;
p = &yy[M + 1];
#ifdef IBMPC
for( i=0; i<7; i++ )
	*p++ = *(--e);
#else  /* IBMPC */
++e;
for( i=0; i<7; i++ )
	*p++ = *e++;
#endif /* IBMPC */ 
/* If denormal, remove the implied bit; else shift down 1. */
if( r == 0 )
	{
	yy[M] = 0;
	}
else
	{
	yy[M] = 1;
	eshift( yy, -1 );
	}
emovo(yy,y,ldp);
}

/* move out internal format to ieee long double */
static void toe113(short unsigned int *a, short unsigned int *b)
{
register unsigned short *p, *q;
unsigned short i;

#ifdef NANS
if( eiisnan(a) )
	{
	enan( b, 113 );
	return;
	}
#endif
p = a;
#ifdef MIEEE
q = b;
#else
q = b + 7;			/* point to output exponent */
#endif

/* If not denormal, delete the implied bit. */
if( a[E] != 0 )
	{
	eshup1 (a);
	}
/* combine sign and exponent */
i = *p++;
#ifdef MIEEE
if( i )
	*q++ = *p++ | 0x8000;
else
	*q++ = *p++;
#else
if( i )
	*q-- = *p++ | 0x8000;
else
	*q-- = *p++;
#endif
/* skip over guard word */
++p;
/* move the significand */
#ifdef MIEEE
for (i = 0; i < 7; i++)
	*q++ = *p++;
#else
for (i = 0; i < 7; i++)
	*q-- = *p++;
#endif
}
#endif /* LDBL_MANT_DIG > 64 */


#if LDBL_MANT_DIG == 64
static void e64toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp)
{
unsigned short yy[NI];
unsigned short *p, *q, *e;
int i;

e = pe;
p = yy;

for( i=0; i<NE-5; i++ )
	*p++ = 0;
#ifdef IBMPC
for( i=0; i<5; i++ )
	*p++ = *e++;
#endif
#ifdef DEC
for( i=0; i<5; i++ )
	*p++ = *e++;
#endif
#ifdef MIEEE
p = &yy[0] + (NE-1);
*p-- = *e++;
++e;  /* MIEEE skips over 2nd short */
for( i=0; i<4; i++ )
	*p-- = *e++;
#endif

#ifdef IBMPC
/* For Intel long double, shift denormal significand up 1
   -- but only if the top significand bit is zero.  */
if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
  {
    unsigned short temp[NI+1];
    emovi(yy, temp);
    eshup1(temp);
    emovo(temp,y,ldp);
    return;
  }
#endif
#ifdef USE_INFINITY
/* Point to the exponent field.  */
p = &yy[NE-1];
if( (*p & 0x7fff) == 0x7fff )
	{
#ifdef NANS
#ifdef IBMPC
	for( i=0; i<4; i++ )
		{
		if((i != 3 && pe[i] != 0)
		   /* Check for Intel long double infinity pattern.  */
		   || (i == 3 && pe[i] != 0x8000))
			{
			enan( y, NBITS );
			return;
			}
		}
#endif
#ifdef MIEEE
	for( i=2; i<=5; i++ )
		{
		if( pe[i] != 0 )
			{
			enan( y, NBITS );
			return;
			}
		}
#endif
#endif /* NANS */
	eclear( y );
	einfin( y, ldp );
	if( *p & 0x8000 )
		eneg(y);
	return;
	}
#endif /* USE_INFINITY */
p = yy;
q = y;
for( i=0; i<NE; i++ )
	*q++ = *p++;
}

/* move out internal format to ieee long double */
static void toe64(short unsigned int *a, short unsigned int *b)
{
register unsigned short *p, *q;
unsigned short i;

#ifdef NANS
if( eiisnan(a) )
	{
	enan( b, 64 );
	return;
	}
#endif
#ifdef IBMPC
/* Shift Intel denormal significand down 1.  */
if( a[E] == 0 )
  eshdn1(a);
#endif
p = a;
#ifdef MIEEE
q = b;
#else
q = b + 4; /* point to output exponent */
/* NOTE: Intel data type is 96 bits wide, clear the last word here. */
*(q+1)= 0;
#endif

/* combine sign and exponent */
i = *p++;
#ifdef MIEEE
if( i )
	*q++ = *p++ | 0x8000;
else
	*q++ = *p++;
*q++ = 0; /* leave 2nd short blank */
#else
if( i )
	*q-- = *p++ | 0x8000;
else
	*q-- = *p++;
#endif
/* skip over guard word */
++p;
/* move the significand */
#ifdef MIEEE
for( i=0; i<4; i++ )
	*q++ = *p++;
#else
#ifdef USE_INFINITY
#ifdef IBMPC
if (eiisinf (a))
        {
	/* Intel long double infinity.  */
	*q-- = 0x8000;
	*q-- = 0;
	*q-- = 0;
	*q = 0;
	return;
	}
#endif /* IBMPC */
#endif /* USE_INFINITY */
for( i=0; i<4; i++ )
	*q-- = *p++;
#endif
}

#endif /* LDBL_MANT_DIG == 64 */

#if LDBL_MANT_DIG == 53
/*
; Convert IEEE double precision to e type
;	double d;
;	unsigned short x[N+2];
;	e53toe( &d, x );
*/
static void e53toe(short unsigned int *pe, short unsigned int *y, LDPARMS *ldp)
{
#ifdef DEC

dectoe( pe, y ); /* see etodec.c */

#else

register unsigned short r;
register unsigned short *p, *e;
unsigned short yy[NI];
int denorm, k;

e = pe;
denorm = 0;	/* flag if denormalized number */
ecleaz(yy);
#ifdef IBMPC
e += 3;
#endif
#ifdef DEC
e += 3;
#endif 
r = *e;
yy[0] = 0;
if( r & 0x8000 )
	yy[0] = 0xffff;
yy[M] = (r & 0x0f) | 0x10;
r &= ~0x800f;	/* strip sign and 4 significand bits */
#ifdef USE_INFINITY
if( r == 0x7ff0 )
	{
#ifdef NANS
#ifdef IBMPC
	if( ((pe[3] & 0xf) != 0) || (pe[2] != 0)
		|| (pe[1] != 0) || (pe[0] != 0) )
		{
		enan( y, NBITS );
		return;
		}
#else  /* !IBMPC */
	if( ((pe[0] & 0xf) != 0) || (pe[1] != 0)
		 || (pe[2] != 0) || (pe[3] != 0) )
		{
		enan( y, NBITS );
		return;
		}
#endif /* !IBMPC */
#endif  /* NANS */
	eclear( y );
	einfin( y, ldp );
	if( yy[0] )
		eneg(y);
	return;
	}
#endif
r >>= 4;
/* If zero exponent, then the significand is denormalized.
 * So, take back the understood high significand bit. */ 
if( r == 0 )
	{
	denorm = 1;
	yy[M] &= ~0x10;
	}
r += EXONE - 01777;
yy[E] = r;
p = &yy[M+1];
#ifdef IBMPC
*p++ = *(--e);
*p++ = *(--e);
*p++ = *(--e);
#else  /* !IBMPC */
++e;
*p++ = *e++;
*p++ = *e++;
*p++ = *e++;
#endif /* !IBMPC */
(void )eshift( yy, -5 );
if( denorm )