ieee.c

128位长双精度型数字运算包

	for( j=0; j<NE-1; j++ )
		{
		if( t[j] != w[j] )
			goto noint;
		}
	emov( t, u );
	expon += (int )m;
noint:
	p += NE;
	m >>= 1;
	}
while( m != 0 );

/* Rescale from integer significand */
	u[NE-1] += y[NE-1] - (unsigned int )(EXONE + NBITS - 1);
	emov( u, y );
/* Find power of 10 */
	emov( eone, t );
	m = MAXP;
	p = &etens[0][0];
	while( ecmp( ten, u ) <= 0 )
		{
		if( ecmp( p, u ) <= 0 )
			{
			ediv( p, u, u );
			emul( p, t, t );
			expon += (int )m;
			}
		m >>= 1;
		if( m == 0 )
			break;
		p += NE;
		}
	}
else
	{ /* Number is less than 1.0 */
/* Pad significand with trailing decimal zeros. */
	if( y[NE-1] == 0 )
		{
		while( (y[NE-2] & 0x8000) == 0 )
			{
			emul( ten, y, y );
			expon -= 1;
			}
		}
	else
		{
		emovi( y, w );
		for( i=0; i<NDEC+1; i++ )
			{
			if( (w[NI-1] & 0x7) != 0 )
				break;
/* multiply by 10 */
			emovz( w, u );
			eshdn1( u );
			eshdn1( u );
			eaddm( w, u );
			u[1] += 3;
			while( u[2] != 0 )
				{
				eshdn1(u);
				u[1] += 1;
				}
			if( u[NI-1] != 0 )
				break;
			if( eone[NE-1] <= u[1] )
				break;
			emovz( u, w );
			expon -= 1;
			}
		emovo( w, y );
		}
	k = -MAXP;
	p = &emtens[0][0];
	r = &etens[0][0];
	emov( y, w );
	emov( eone, t );
	while( ecmp( eone, w ) > 0 )
		{
		if( ecmp( p, w ) >= 0 )
			{
			emul( r, w, w );
			emul( r, t, t );
			expon += k;
			}
		k /= 2;
		if( k == 0 )
			break;
		p += NE;
		r += NE;
		}
	ediv( t, eone, t );
	}
isone:
/* Find the first (leading) digit. */
emovi( t, w );
emovz( w, t );
emovi( y, w );
emovz( w, y );
eiremain( t, y );
digit = equot[NI-1];
while( (digit == 0) && (ecmp(y,ezero) != 0) )
	{
	eshup1( y );
	emovz( y, u );
	eshup1( u );
	eshup1( u );
	eaddm( u, y );
	eiremain( t, y );
	digit = equot[NI-1];
	expon -= 1;
	}
s = string;
if( sign )
	*s++ = '-';
else
	*s++ = ' ';
/* Examine number of digits requested by caller. */
if( ndigs < 0 )
	ndigs = 0;
if( ndigs > NDEC )
	ndigs = NDEC;
if( digit == 10 )
	{
	*s++ = '1';
	*s++ = '.';
	if( ndigs > 0 )
		{
		*s++ = '0';
		ndigs -= 1;
		}
	expon += 1;
	}
else
	{
	*s++ = (char )digit + '0';
	*s++ = '.';
	}
/* Generate digits after the decimal point. */
for( k=0; k<=ndigs; k++ )
	{
/* multiply current number by 10, without normalizing */
	eshup1( y );
	emovz( y, u );
	eshup1( u );
	eshup1( u );
	eaddm( u, y );
	eiremain( t, y );
	*s++ = (char )equot[NI-1] + '0';
	}
digit = equot[NI-1];
--s;
ss = s;
/* round off the ASCII string */
if( digit > 4 )
	{
/* Test for critical rounding case in ASCII output. */
	if( digit == 5 )
		{
		emovo( y, t );
		if( ecmp(t,ezero) != 0 )
			goto roun;	/* round to nearest */
		if( (*(s-1) & 1) == 0 )
			goto doexp;	/* round to even */
		}
/* Round up and propagate carry-outs */
roun:
	--s;
	k = *s & 0x7f;
/* Carry out to most significant digit? */
	if( k == '.' )
		{
		--s;
		k = *s;
		k += 1;
		*s = (char )k;
/* Most significant digit carries to 10? */
		if( k > '9' )
			{
			expon += 1;
			*s = '1';
			}
		goto doexp;
		}
/* Round up and carry out from less significant digits */
	k += 1;
	*s = (char )k;
	if( k > '9' )
		{
		*s = '0';
		goto roun;
		}
	}
doexp:
/*
if( expon >= 0 )
	sprintf( ss, "e+%d", expon );
else
	sprintf( ss, "e%d", expon );
*/
	sprintf( ss, "E%d", expon );
bxit:
rndprc = rndsav;
}




/*
;								ASCTOQ
;		ASCTOQ.MAC		LATEST REV: 11 JAN 84
;					SLM, 3 JAN 78
;
;	Convert ASCII string to quadruple precision floating point
;
;		Numeric input is free field decimal number
;		with max of 15 digits with or without 
;		decimal point entered as ASCII from teletype.
;	Entering E after the number followed by a second
;	number causes the second number to be interpreted
;	as a power of 10 to be multiplied by the first number
;	(i.e., "scientific" notation).
;
;	Usage:
;		asctoq( string, q );
*/

/* ASCII to single */
void asctoe24( s, y )
char *s;
unsigned short *y;
{
asctoeg( s, y, 24 );
}


/* ASCII to double */
void asctoe53( s, y )
char *s;
unsigned short *y;
{
#ifdef DEC
asctoeg( s, y, 56 );
#else
asctoeg( s, y, 53 );
#endif
}


/* ASCII to long double */
void asctoe64( s, y )
char *s;
unsigned short *y;
{
asctoeg( s, y, 64 );
}

/* ASCII to 128-bit long double */
void asctoe113 (s, y)
char *s;
unsigned short *y;
{
asctoeg( s, y, 113 );
}

/* ASCII to super double */
void asctoe( s, y )
char *s;
unsigned short *y;
{
asctoeg( s, y, NBITS );
}

/* Space to make a copy of the input string: */
static char lstr[82] = {0};

void asctoeg( ss, y, oprec )
char *ss;
unsigned short *y;
int oprec;
{
unsigned short yy[NI], xt[NI], tt[NI];
int esign, decflg, sgnflg, nexp, exp, prec, lost;
int k, trail, c, rndsav;
long lexp;
unsigned short nsign, *p;
char *sp, *s;

/* Copy the input string. */
s = ss;
while( *s == ' ' ) /* skip leading spaces */
	++s;
sp = lstr;
for( k=0; k<79; k++ )
	{
	if( (*sp++ = *s++) == '\0' )
		break;
	}
*sp = '\0';
s = lstr;

rndsav = rndprc;
rndprc = NBITS; /* Set to full precision */
lost = 0;
nsign = 0;
decflg = 0;
sgnflg = 0;
nexp = 0;
exp = 0;
prec = 0;
ecleaz( yy );
trail = 0;

nxtcom:
k = *s - '0';
if( (k >= 0) && (k <= 9) )
	{
/* Ignore leading zeros */
	if( (prec == 0) && (decflg == 0) && (k == 0) )
		goto donchr;
/* Identify and strip trailing zeros after the decimal point. */
	if( (trail == 0) && (decflg != 0) )
		{
		sp = s;
		while( (*sp >= '0') && (*sp <= '9') )
			++sp;
/* Check for syntax error */
		c = *sp & 0x7f;
		if( (c != 'e') && (c != 'E') && (c != '\0')
			&& (c != '\n') && (c != '\r') && (c != ' ')
			&& (c != ',') )
			goto error;
		--sp;
		while( *sp == '0' )
			*sp-- = 'z';
		trail = 1;
		if( *s == 'z' )
			goto donchr;
		}
/* If enough digits were given to more than fill up the yy register,
 * continuing until overflow into the high guard word yy[2]
 * guarantees that there will be a roundoff bit at the top
 * of the low guard word after normalization.
 */
	if( yy[2] == 0 )
		{
		if( decflg )
			nexp += 1; /* count digits after decimal point */
		eshup1( yy );	/* multiply current number by 10 */
		emovz( yy, xt );
		eshup1( xt );
		eshup1( xt );
		eaddm( xt, yy );
		ecleaz( xt );
		xt[NI-2] = (unsigned short )k;
		eaddm( xt, yy );
		}
	else
		{
		/* Mark any lost non-zero digit.  */
		lost |= k;
		/* Count lost digits before the decimal point.  */
		if (decflg == 0)
		        nexp -= 1;
		}
	prec += 1;
	goto donchr;
	}

switch( *s )
	{
	case 'z':
		break;
	case 'E':
	case 'e':
		goto expnt;
	case '.':	/* decimal point */
		if( decflg )
			goto error;
		++decflg;
		break;
	case '-':
		nsign = 0xffff;
		if( sgnflg )
			goto error;
		++sgnflg;
		break;
	case '+':
		if( sgnflg )
			goto error;
		++sgnflg;
		break;
	case ',':
	case ' ':
	case '\0':
	case '\n':
	case '\r':
		goto daldone;
	case 'i':
	case 'I':
		goto infinite;
	default:
	error:
#ifdef NANS
		enan( yy, NI*16 );
#else
		mtherr( "asctoe", DOMAIN );
		ecleaz(yy);
#endif
		goto aexit;
	}
donchr:
++s;
goto nxtcom;

/* Exponent interpretation */
expnt:

esign = 1;
exp = 0;
++s;
/* check for + or - */
if( *s == '-' )
	{
	esign = -1;
	++s;
	}
if( *s == '+' )
	++s;
while( (*s >= '0') && (*s <= '9') )
	{
	exp *= 10;
	exp += *s++ - '0';
	if (exp > 4977)
		{
		if (esign < 0)
			goto zero;
		else
			goto infinite;
		}
	}
if( esign < 0 )
	exp = -exp;
if( exp > 4932 )
	{
infinite:
	ecleaz(yy);
	yy[E] = 0x7fff;  /* infinity */
	goto aexit;
	}
if( exp < -4977 )
	{
zero:
	ecleaz(yy);
	goto aexit;
	}

daldone:
nexp = exp - nexp;
/* Pad trailing zeros to minimize power of 10, per IEEE spec. */
while( (nexp > 0) && (yy[2] == 0) )
	{
	emovz( yy, xt );
	eshup1( xt );
	eshup1( xt );
	eaddm( yy, xt );
	eshup1( xt );
	if( xt[2] != 0 )
		break;
	nexp -= 1;
	emovz( xt, yy );
	}
if( (k = enormlz(yy)) > NBITS )
	{
	ecleaz(yy);
	goto aexit;
	}
lexp = (EXONE - 1 + NBITS) - k;
emdnorm( yy, lost, 0, lexp, 64 );
/* convert to external format */


/* Multiply by 10**nexp.  If precision is 64 bits,
 * the maximum relative error incurred in forming 10**n
 * for 0 <= n <= 324 is 8.2e-20, at 10**180.
 * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
 * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
 */
lexp = yy[E];
if( nexp == 0 )
	{
	k = 0;
	goto expdon;
	}
esign = 1;
if( nexp < 0 )
	{
	nexp = -nexp;
	esign = -1;
	if( nexp > 4096 )
		{ /* Punt.  Can't handle this without 2 divides. */
		emovi( etens[0], tt );
		lexp -= tt[E];
		k = edivm( tt, yy );
		lexp += EXONE;
		nexp -= 4096;
		}
	}
p = &etens[NTEN][0];
emov( eone, xt );
exp = 1;
do
	{
	if( exp & nexp )
		emul( p, xt, xt );
	p -= NE;
	exp = exp + exp;
	}
while( exp <= MAXP );

emovi( xt, tt );
if( esign < 0 )
	{
	lexp -= tt[E];
	k = edivm( tt, yy );
	lexp += EXONE;
	}
else
	{
	lexp += tt[E];
	k = emulm( tt, yy );
	lexp -= EXONE - 1;
	}

expdon:

/* Round and convert directly to the destination type */
if( oprec == 53 )
	lexp -= EXONE - 0x3ff;
else if( oprec == 24 )
	lexp -= EXONE - 0177;
#ifdef DEC
else if( oprec == 56 )
	lexp -= EXONE - 0201;
#endif
rndprc = oprec;
emdnorm( yy, k, 0, lexp, 64 );

aexit:

rndprc = rndsav;
yy[0] = nsign;
switch( oprec )
	{
#ifdef DEC
	case 56:
		todec( yy, y ); /* see etodec.c */
		break;
#endif
	case 53:
		toe53( yy, y );
		break;
	case 24:
		toe24( yy, y );
		break;
	case 64:
		toe64( yy, y );
		break;
	case 113:
		toe113( yy, y );
		break;
	case NBITS:
		emovo( yy, y );
		break;
	}
}


 
/* y = largest integer not greater than x
 * (truncated toward minus infinity)
 *
 * unsigned short x[NE], y[NE]
 *
 * efloor( x, y );
 */
static unsigned short bmask[] = {
0xffff,
0xfffe,
0xfffc,
0xfff8,
0xfff0,
0xffe0,
0xffc0,
0xff80,
0xff00,
0xfe00,
0xfc00,
0xf800,
0xf000,
0xe000,
0xc000,
0x8000,
0x0000,
};

void efloor( x, y )
unsigned short x[], y[];
{
register unsigned short *p;
int e, expon, i;
unsigned short f[NE];

emov( x, f ); /* leave in external format */
expon = (int )f[NE-1];
e = (expon & 0x7fff) - (EXONE - 1);
if( e <= 0 )
	{
	eclear(y);
	goto isitneg;
	}
/* number of bits to clear out */
e = NBITS - e;
emov( f, y );
if( e <= 0 )
	return;

p = &y[0];
while( e >= 16 )
	{
	*p++ = 0;
	e -= 16;
	}
/* clear the remaining bits */
*p &= bmask[e];
/* truncate negatives toward minus infinity */
isitneg:

if( (unsigned short )expon & (unsigned short )0x8000 )
	{
	for( i=0; i<NE-1; i++ )
		{
		if( f[i] != y[i] )
			{
			esub( eone, y, y );
			break;
			}
		}
	}
}


/* unsigned short x[], s[];
 * long *exp;
 *
 * efrexp( x, exp, s );
 *
 * Returns s and exp such that  s * 2**exp = x and .5 <= s < 1.
 * For example, 1.1 = 0.55 * 2**1
 * Handles denormalized numbers properly using long integer exp.
 */
void efrexp( x, exp, s )
unsigned short x[];
long *exp;
unsigned short s[];
{
unsigned short xi[NI];
long li;

emovi( x, xi );
li = (long )((short )xi[1]);

if( li == 0 )
	{
	li -= enormlz( xi );
	}
xi[1] = 0x3ffe;
emovo( xi, s );
*exp = li - 0x3ffe;
}



/* unsigned short x[], y[];
 * long pwr2;
 *
 * eldexp( x, pwr2, y );
 *
 * Returns y = x * 2**pwr2.
 */
void eldexp( x, pwr2, y )
unsigned short x[];
long pwr2;
unsigned short y[];
{
unsigned short xi[NI];
long li;
int i;

emovi( x, xi );
li = xi[1];
li += pwr2;
i = 0;
emdnorm( xi, i, i, li, 64 );
emovo( xi, y );
}


/* c = remainder after dividing b by a
 * Least significant integer quotient bits left in equot[].
 */
void eremain( a, b, c )
unsigned short a[], b[], c[];
{
unsigned short den[NI], num[NI];

#ifdef NANS
if( eisinf(b) || (ecmp(a,ezero) == 0) || eisnan(a) || eisnan(b))
	{
	enan( c, NBITS );
	return;
	}
#endif
if( ecmp(a,ezero) == 0 )
	{
	mtherr( "eremain", SING );
	eclear( c );
	return;
	}
emovi( a, den );
emovi( b, num );
eiremain( den, num );
/* Sign of remainder = sign of quotient */
if( a[0] == b[0] )
	num[0] = 0;
else
	num[0] = 0xffff;
emovo( num, c );
}


void eiremain( den, num )
unsigned short den[], num[];
{
long ld, ln;
unsigned short j;

ld = den[E];
ld -= enormlz( den );
ln = num[E];
ln -= enormlz( num );
ecleaz( equot );
while( ln >= ld )
	{
	if( ecmpm(den,num) <= 0 )
		{
		esubm(den, num);
		j = 1;
		}
	else
		{
		j = 0;
		}
	eshup1(equot);
	equot[NI-1] |= j;
	eshup1(num);
	ln -= 1;
	}
emdnorm( num, 0, 0, ln, 0 );
}

/* NaN bit patterns
 */
#ifdef MIEEE
unsigned short nan113[8] = {
  0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
unsigned short nan64[6] = {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
unsigned short nan53[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
unsigned short nan24[2] = {0x7fff, 0xffff};
#endif

#ifdef IBMPC
unsigned short nan113[8] = {0, 0, 0, 0, 0, 0, 0xc000, 0xffff};
unsigned short nan64[6] = {0, 0, 0, 0xc000, 0xffff, 0};
unsigned short nan53[4] = {0, 0, 0, 0xfff8};
unsigned short nan24[2] = {0, 0xffc0};
#endif


void enan (nan, size)
unsigned short *nan;
int size;
{
i