ieee.c

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

/*							ieee.c
 *
 *    Extended precision IEEE binary floating point arithmetic routines
 *
 * Numbers are stored in C language as arrays of 16-bit unsigned
 * short integers.  The arguments of the routines are pointers to
 * the arrays.
 *
 *
 * External e type data structure, simulates Intel 8087 chip
 * temporary real format but possibly with a larger significand:
 *
 *	NE-1 significand words	(least significant word first,
 *				 most significant bit is normally set)
 *	exponent		(value = EXONE for 1.0,
 *				top bit is the sign)
 *
 *
 * Internal data structure of a number (a "word" is 16 bits):
 *
 * ei[0]	sign word	(0 for positive, 0xffff for negative)
 * ei[1]	biased exponent	(value = EXONE for the number 1.0)
 * ei[2]	high guard word	(always zero after normalization)
 * ei[3]
 * to ei[NI-2]	significand	(NI-4 significand words,
 *				 most significant word first,
 *				 most significant bit is set)
 * ei[NI-1]	low guard word	(0x8000 bit is rounding place)
 *
 *
 *
 *		Routines for external format numbers
 *
 *	asctoe( string, e )	ASCII string to extended double e type
 *	asctoe64( string, &d )	ASCII string to long double
 *	asctoe53( string, &d )	ASCII string to double
 *	asctoe24( string, &f )	ASCII string to single
 *	asctoeg( string, e, prec ) ASCII string to specified precision
 *	e24toe( &f, e )		IEEE single precision to e type
 *	e53toe( &d, e )		IEEE double precision to e type
 *	e64toe( &d, e )		IEEE long double precision to e type
 *	eabs(e)			absolute value
 *	eadd( a, b, c )		c = b + a
 *	eclear(e)		e = 0
 *	ecmp (a, b)		Returns 1 if a > b, 0 if a == b,
 *				-1 if a < b, -2 if either a or b is a NaN.
 *	ediv( a, b, c )		c = b / a
 *	efloor( a, b )		truncate to integer, toward -infinity
 *	efrexp( a, exp, s )	extract exponent and significand
 *	eifrac( e, &l, frac )   e to long integer and e type fraction
 *	euifrac( e, &l, frac )  e to unsigned long integer and e type fraction
 *	einfin( e )		set e to infinity, leaving its sign alone
 *	eldexp( a, n, b )	multiply by 2**n
 *	emov( a, b )		b = a
 *	emul( a, b, c )		c = b * a
 *	eneg(e)			e = -e
 *	eround( a, b )		b = nearest integer value to a
 *	esub( a, b, c )		c = b - a
 *	e24toasc( &f, str, n )	single to ASCII string, n digits after decimal
 *	e53toasc( &d, str, n )	double to ASCII string, n digits after decimal
 *	e64toasc( &d, str, n )	long double to ASCII string
 *	etoasc( e, str, n )	e to ASCII string, n digits after decimal
 *	etoe24( e, &f )		convert e type to IEEE single precision
 *	etoe53( e, &d )		convert e type to IEEE double precision
 *	etoe64( e, &d )		convert e type to IEEE long double precision
 *	ltoe( &l, e )		long (32 bit) integer to e type
 *	ultoe( &l, e )		unsigned long (32 bit) integer to e type
 *      eisneg( e )             1 if sign bit of e != 0, else 0
 *      eisinf( e )             1 if e has maximum exponent (non-IEEE)
 *				or is infinite (IEEE)
 *      eisnan( e )             1 if e is a NaN
 *	esqrt( a, b )		b = square root of a
 *
 *
 *		Routines for internal format numbers
 *
 *	eaddm( ai, bi )		add significands, bi = bi + ai
 *	ecleaz(ei)		ei = 0
 *	ecleazs(ei)		set ei = 0 but leave its sign alone
 *	ecmpm( ai, bi )		compare significands, return 1, 0, or -1
 *	edivm( ai, bi )		divide  significands, bi = bi / ai
 *	emdnorm(ai,l,s,exp)	normalize and round off
 *	emovi( a, ai )		convert external a to internal ai
 *	emovo( ai, a )		convert internal ai to external a
 *	emovz( ai, bi )		bi = ai, low guard word of bi = 0
 *	emulm( ai, bi )		multiply significands, bi = bi * ai
 *	enormlz(ei)		left-justify the significand
 *	eshdn1( ai )		shift significand and guards down 1 bit
 *	eshdn8( ai )		shift down 8 bits
 *	eshdn6( ai )		shift down 16 bits
 *	eshift( ai, n )		shift ai n bits up (or down if n < 0)
 *	eshup1( ai )		shift significand and guards up 1 bit
 *	eshup8( ai )		shift up 8 bits
 *	eshup6( ai )		shift up 16 bits
 *	esubm( ai, bi )		subtract significands, bi = bi - ai
 *
 *
 * The result is always normalized and rounded to NI-4 word precision
 * after each arithmetic operation.
 *
 * Exception flags are NOT fully supported.
 *
 * Define INFINITIES in mconf.h for support of infinity; otherwise a
 * saturation arithmetic is implemented.
 *
 * Define NANS for support of Not-a-Number items; otherwise the
 * arithmetic will never produce a NaN output, and might be confused
 * by a NaN input.
 * If NaN's are supported, the output of ecmp(a,b) is -2 if
 * either a or b is a NaN. This means asking if(ecmp(a,b) < 0)
 * may not be legitimate. Use if(ecmp(a,b) == -1) for less-than
 * if in doubt.
 * Signaling NaN's are NOT supported; they are treated the same
 * as quiet NaN's.
 *
 * Denormals are always supported here where appropriate (e.g., not
 * for conversion to DEC numbers).
 */

/*
 * Revision history:
 *
 *  5 Jan 84	PDP-11 assembly language version
 *  2 Mar 86	fixed bug in asctoq()
 *  6 Dec 86	C language version
 * 30 Aug 88	100 digit version, improved rounding
 * 15 May 92    80-bit long double support
 *
 * Author:  S. L. Moshier.
 */

#include <stdio.h>
#include "mconf.h"
#include "ehead.h"

/* Change UNK into something else. */
#ifdef UNK
#undef UNK
#if BIGENDIAN
#define MIEEE 1
#else
#define IBMPC 1
#endif
#endif

/* NaN's require infinity support. */
#ifdef NANS
#ifndef INFINITIES
#define INFINITIES
#endif
#endif

/* This handles 64-bit long ints. */
#define LONGBITS (8 * sizeof(long))

/* Control register for rounding precision.
 * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
 */
int rndprc = NBITS;
extern int rndprc;

void eaddm(), esubm(), emdnorm(), asctoeg(), enan();
static void toe24(), toe53(), toe64(), toe113();
void eremain(), einit(), eiremain();
int ecmpm(), edivm(), emulm(), eisneg(), eisinf();
void emovi(), emovo(), emovz(), ecleaz(), eadd1();
void etodec(), todec(), dectoe();
int eisnan(), eiisnan();



void einit()
{
}

/*
; Clear out entire external format number.
;
; unsigned short x[];
; eclear( x );
*/

void eclear( x )
register unsigned short *x;
{
register int i;

for( i=0; i<NE; i++ )
	*x++ = 0;
}



/* Move external format number from a to b.
 *
 * emov( a, b );
 */

void emov( a, b )
register unsigned short *a, *b;
{
register int i;

for( i=0; i<NE; i++ )
	*b++ = *a++;
}


/*
;	Absolute value of external format number
;
;	short x[NE];
;	eabs( x );
*/

void eabs(x)
unsigned short x[];	/* x is the memory address of a short */
{

x[NE-1] &= 0x7fff; /* sign is top bit of last word of external format */
}




/*
;	Negate external format number
;
;	unsigned short x[NE];
;	eneg( x );
*/

void eneg(x)
unsigned short x[];
{

#ifdef NANS
if( eisnan(x) )
	return;
#endif
x[NE-1] ^= 0x8000; /* Toggle the sign bit */
}



/* Return 1 if external format number is negative,
 * else return zero.
 */
int eisneg(x)
unsigned short x[];
{

#ifdef NANS
if( eisnan(x) )
	return( 0 );
#endif
if( x[NE-1] & 0x8000 )
	return( 1 );
else
	return( 0 );
}


/* Return 1 if external format number has maximum possible exponent,
 * else return zero.
 */
int eisinf(x)
unsigned short x[];
{

if( (x[NE-1] & 0x7fff) == 0x7fff )
	{
#ifdef NANS
	if( eisnan(x) )
		return( 0 );
#endif
	return( 1 );
	}
else
	return( 0 );
}

/* Check if e-type number is not a number.
 */
int eisnan(x)
unsigned short x[];
{

#ifdef NANS
int i;
/* NaN has maximum exponent */
if( (x[NE-1] & 0x7fff) != 0x7fff )
	return (0);
/* ... and non-zero significand field. */
for( i=0; i<NE-1; i++ )
	{
	if( *x++ != 0 )
		return (1);
	}
#endif
return (0);
}

/*
; Fill entire number, including exponent and significand, with
; largest possible number.  These programs implement a saturation
; value that is an ordinary, legal number.  A special value
; "infinity" may also be implemented; this would require tests
; for that value and implementation of special rules for arithmetic
; operations involving inifinity.
*/

void einfin(x)
register unsigned short *x;
{
register int i;

#ifdef INFINITIES
for( i=0; i<NE-1; i++ )
	*x++ = 0;
*x |= 32767;
#else
for( i=0; i<NE-1; i++ )
	*x++ = 0xffff;
*x |= 32766;
if( rndprc < NBITS )
	{
	if (rndprc == 113)
		{
		*(x - 9) = 0;
		*(x - 8) = 0;
		}
	if( rndprc == 64 )
		{
		*(x-5) = 0;
		}
	if( rndprc == 53 )
		{
		*(x-4) = 0xf800;
		}
	else
		{
		*(x-4) = 0;
		*(x-3) = 0;
		*(x-2) = 0xff00;
		}
	}
#endif
}



/* Move in external format number,
 * converting it to internal format.
 */
void emovi( a, b )
unsigned short *a, *b;
{
register unsigned short *p, *q;
int i;

q = b;
p = a + (NE-1);	/* point to last word of external number */
/* get the sign bit */
if( *p & 0x8000 )
	*q++ = 0xffff;
else
	*q++ = 0;
/* get the exponent */
*q = *p--;
*q++ &= 0x7fff;	/* delete the sign bit */
#ifdef INFINITIES
if( (*(q-1) & 0x7fff) == 0x7fff )
	{
#ifdef NANS
	if( eisnan(a) )
		{
		*q++ = 0;
		for( i=3; i<NI; i++ )
			*q++ = *p--;
		return;
		}
#endif
	for( i=2; i<NI; i++ )
		*q++ = 0;
	return;
	}
#endif
/* clear high guard word */
*q++ = 0;
/* move in the significand */
for( i=0; i<NE-1; i++ )
	*q++ = *p--;
/* clear low guard word */
*q = 0;
}


/* Move internal format number out,
 * converting it to external format.
 */
void emovo( a, b )
unsigned short *a, *b;
{
register unsigned short *p, *q;
unsigned short i;

p = a;
q = b + (NE-1); /* point to output exponent */
/* combine sign and exponent */
i = *p++;
if( i )
	*q-- = *p++ | 0x8000;
else
	*q-- = *p++;
#ifdef INFINITIES
if( *(p-1) == 0x7fff )
	{
#ifdef NANS
	if( eiisnan(a) )
		{
		enan( b, NBITS );
		return;
		}
#endif
	einfin(b);
	return;
	}
#endif
/* skip over guard word */
++p;
/* move the significand */
for( i=0; i<NE-1; i++ )
	*q-- = *p++;
}




/* Clear out internal format number.
 */

void ecleaz( xi )
register unsigned short *xi;
{
register int i;

for( i=0; i<NI; i++ )
	*xi++ = 0;
}

/* same, but don't touch the sign. */

void ecleazs( xi )
register unsigned short *xi;
{
register int i;

++xi;
for(i=0; i<NI-1; i++)
	*xi++ = 0;
}




/* Move internal format number from a to b.
 */
void emovz( a, b )
register unsigned short *a, *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.
 */

int eiisnan (x)
unsigned short x[];
{
int i;

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

#ifdef INFINITIES
/* 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

/*
;	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
*/
int ecmpm( a, b )
register unsigned short *a, *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
*/

void eshdn1(x)
register unsigned short *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
*/

void eshup1(x)
register unsigned short *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
*/

void eshdn8(x)
register unsigned short *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
*/

void eshup8(x)
register unsigned short *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
*/

void eshup6(x)
register unsigned short *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
*/

void eshdn6(x)
register unsigned short *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
*/

void eaddm( x, y )
unsigned short *x, *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
*/

void esubm( x, y )
unsigned short *x, *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 */

static unsigned short equot[NI] = {0}; /* was static */

#if 0
int edivm( den, num )
unsigned short den[], num[];
{
int i;
register unsigned short *p, *q;
unsigned short j;

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

for( i=M; i<NI; i++ )
	{
	*p++ = 0;
	}

/* Use faster compare and subtraction if denominator
 * has only 15 bits of significance.
 */
p = &den[M+2];
if( *p++ == 0 )
	{
	for( i=M+3; i<NI; i++ )
		{
		if( *p++ != 0 )
			goto fulldiv;
		}
	if( (den[M+1] & 1) != 0 )
		goto fulldiv;
	eshdn1(num);
	eshdn1(den);

	p = &den[M+1];
	q = &num[M+1];

	for( i=0; i<NBITS+2; i++ )
		{
		if( *p <= *q )
			{
			*q -= *p;
			j = 1;
			}
		else
			{
			j = 0;
			}