ieee.c

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

divsign:

if( sign )
	*(c+(NE-1)) |= 0x8000;
else
	*(c+(NE-1)) &= ~0x8000;
}



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

/* IEEE says if result is not a NaN, the sign is "-" if and only if
   operands have opposite signs -- but flush -0 to 0 later if not IEEE.  */
sign = eisneg(a) ^ eisneg(b);

#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 INFINITIES
if( eisinf(a) || eisinf(b) )
	{
	einfin(c);
	goto mulsign;
	}
#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);
	goto mulsign;
	}
mnzer1:

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

/* Multiply significands */
j = emulm( ai, bi );
/* calculate exponent */
lt = lta + ltb - (EXONE - 1);
emdnorm( bi, j, 0, lt, 64 );
emovo( bi, c );
/*  IEEE says sign is "-" if and only if operands have opposite signs.  */
mulsign:
if( sign )
	*(c+(NE-1)) |= 0x8000;
else
	*(c+(NE-1)) &= ~0x8000;
}




/*
; Convert IEEE double precision to e type
;	double d;
;	unsigned short x[N+2];
;	e53toe( &d, x );
*/
void e53toe( pe, y )
unsigned short *pe, *y;
{
#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
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 INFINITIES
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
	if( ((pe[0] & 0xf) != 0) || (pe[1] != 0)
		 || (pe[2] != 0) || (pe[3] != 0) )
		{
		enan( y, NBITS );
		return;
		}
#endif
#endif  /* NANS */
	eclear( y );
	einfin( y );
	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);
#endif
#ifdef MIEEE
++e;
*p++ = *e++;
*p++ = *e++;
*p++ = *e++;
#endif
(void )eshift( yy, -5 );
if( denorm )
	{ /* if zero exponent, then normalize the significand */
	if( (k = enormlz(yy)) > NBITS )
		ecleazs(yy);
	else
		yy[E] -= (unsigned short )(k-1);
	}
emovo( yy, y );
#endif /* not DEC */
}

void e64toe( pe, y )
unsigned short *pe, *y;
{
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;
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);
    return;
  }
#endif
#ifdef INFINITIES
/* 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;
			}
		}
#else
	/* In Motorola extended precision format, the most significant
	   bit of an infinity mantissa could be either 1 or 0.  It is
	   the lower order bits that tell whether the value is a NaN.  */
	if ((pe[2] & 0x7fff) != 0)
		goto bigend_nan;

	for( i=3; i<=5; i++ )
		{
		if( pe[i] != 0 )
			{
bigend_nan:
			enan( y, NBITS );
			return;
			}
		}
#endif
#endif /* NANS */
	eclear( y );
	einfin( y );
	if( *p & 0x8000 )
		eneg(y);
	return;
	}
#endif
p = yy;
q = y;
for( i=0; i<NE; i++ )
	*q++ = *p++;
}

void e113toe(pe,y)
unsigned short *pe, *y;
{
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 INFINITIES
if( r == 0x7fff )
	{
#ifdef NANS
#ifdef IBMPC
	for( i=0; i<7; i++ )
		{
		if( pe[i] != 0 )
			{
			enan( y, NBITS );
			return;
			}
		}
#else
	for( i=1; i<8; i++ )
		{
		if( pe[i] != 0 )
			{
			enan( y, NBITS );
			return;
			}
		}
#endif
#endif /* NANS */
	eclear( y );
	einfin( y );
	if( *e & 0x8000 )
		eneg(y);
	return;
	}
#endif  /* INFINITIES */
yy[E] = r;
p = &yy[M + 1];
#ifdef IBMPC
for( i=0; i<7; i++ )
	*p++ = *(--e);
#endif
#ifdef MIEEE
++e;
for( i=0; i<7; i++ )
	*p++ = *e++;
#endif
/* 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);
}


/*
; Convert IEEE single precision to e type
;	float d;
;	unsigned short x[N+2];
;	dtox( &d, x );
*/
void e24toe( pe, y )
unsigned short *pe, *y;
{
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 += 1;
#endif
#ifdef DEC
e += 1;
#endif
r = *e;
yy[0] = 0;
if( r & 0x8000 )
	yy[0] = 0xffff;
yy[M] = (r & 0x7f) | 0200;
r &= ~0x807f;	/* strip sign and 7 significand bits */
#ifdef INFINITIES
if( r == 0x7f80 )
	{
#ifdef NANS
#ifdef MIEEE
	if( ((pe[0] & 0x7f) != 0) || (pe[1] != 0) )
		{
		enan( y, NBITS );
		return;
		}
#else
	if( ((pe[1] & 0x7f) != 0) || (pe[0] != 0) )
		{
		enan( y, NBITS );
		return;
		}
#endif
#endif  /* NANS */
	eclear( y );
	einfin( y );
	if( yy[0] )
		eneg(y);
	return;
	}
#endif
r >>= 7;
/* If zero exponent, then the significand is denormalized.
 * So, take back the understood high significand bit. */ 
if( r == 0 )
	{
	denorm = 1;
	yy[M] &= ~0200;
	}
r += EXONE - 0177;
yy[E] = r;
p = &yy[M+1];
#ifdef IBMPC
*p++ = *(--e);
#endif
#ifdef DEC
*p++ = *(--e);
#endif
#ifdef MIEEE
++e;
*p++ = *e++;
#endif
(void )eshift( yy, -8 );
if( denorm )
	{ /* if zero exponent, then normalize the significand */
	if( (k = enormlz(yy)) > NBITS )
		ecleazs(yy);
	else
		yy[E] -= (unsigned short )(k-1);
	}
emovo( yy, y );
}

void etoe113(x,e)
unsigned short *x, *e;
{
unsigned short xi[NI];
long exp;
int rndsav;

#ifdef NANS
if( eisnan(x) )
	{
	enan( e, 113 );
	return;
	}
#endif
emovi( x, xi );
exp = (long )xi[E];
#ifdef INFINITIES
if( eisinf(x) )
	goto nonorm;
#endif
/* round off to nearest or even */
rndsav = rndprc;
rndprc = 113;
emdnorm( xi, 0, 0, exp, 64 );
rndprc = rndsav;
nonorm:
toe113 (xi, e);
}

/* move out internal format to ieee long double */
static void toe113(a,b)
unsigned short *a, *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
}


void etoe64( x, e )
unsigned short *x, *e;
{
unsigned short xi[NI];
long exp;
int rndsav;

#ifdef NANS
if( eisnan(x) )
	{
	enan( e, 64 );
	return;
	}
#endif
emovi( x, xi );
exp = (long )xi[E]; /* adjust exponent for offset */
#ifdef INFINITIES
if( eisinf(x) )
	goto nonorm;
#endif
/* round off to nearest or even */
rndsav = rndprc;
rndprc = 64;
emdnorm( xi, 0, 0, exp, 64 );
rndprc = rndsav;
nonorm:
toe64( xi, e );
}

/* move out internal format to ieee long double */
static void toe64( a, b )
unsigned short *a, *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 */
#if 1
/* NOTE: if data type is 96 bits wide, clear the last word here. */
*(q+1)= 0;
#endif
#endif

/* combine sign and exponent */
i = *p++;
#ifdef MIEEE
if( i )
	*q++ = *p++ | 0x8000;
else
	*q++ = *p++;
*q++ = 0;
#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 INFINITIES
if (eiisinf (a))
        {
	/* Intel long double infinity.  */
	*q-- = 0x8000;
	*q-- = 0;
	*q-- = 0;
	*q = 0;
	return;
	}
#endif
for( i=0; i<4; i++ )
	*q-- = *p++;
#endif
}


/*
; e type to IEEE double precision
;	double d;
;	unsigned short x[NE];
;	etoe53( x, &d );
*/

#ifdef DEC

void etoe53( x, e )
unsigned short *x, *e;
{
etodec( x, e ); /* see etodec.c */
}

static void toe53( x, y )
unsigned short *x, *y;
{
todec( x, y );
}

#else

void etoe53( x, e )
unsigned short *x, *e;
{
unsigned short xi[NI];
long exp;
int rndsav;

#ifdef NANS
if( eisnan(x) )
	{
	enan( e, 53 );
	return;
	}
#endif
emovi( x, xi );
exp = (long )xi[E] - (EXONE - 0x3ff); /* adjust exponent for offsets */
#ifdef INFINITIES
if( eisinf(x) )
	goto nonorm;
#endif
/* round off to nearest or even */
rndsav = rndprc;
rndprc = 53;
emdnorm( xi, 0, 0, exp, 64 );
rndprc = rndsav;
nonorm:
toe53( xi, e );
}


static void toe53( x, y )
unsigned short *x, *y;
{
unsigned short i;
unsigned short *p;


#ifdef NANS
if( eiisnan(x) )
	{
	enan( y, 53 );
	return;
	}
#endif
p = &x[0];
#ifdef IBMPC
y += 3;
#endif
*y = 0;	/* output high order */
if( *p++ )
	*y = 0x8000;	/* output sign bit */

i = *p++;
if( i >= (unsigned int )2047 )
	{	/* Saturate at largest number less than infinity. */
#ifdef INFINITIES
	*y |= 0x7ff0;
#ifdef IBMPC
	*(--y) = 0;
	*(--y) = 0;
	*(--y) = 0;
#endif
#ifdef MIEEE
	++y;
	*y++ = 0;
	*y++ = 0;
	*y++ = 0;
#endif
#else
	*y |= (unsigned short )0x7fef;
#ifdef IBMPC
	*(--y) = 0xffff;
	*(--y) = 0xffff;
	*(--y) = 0xffff;
#endif
#ifdef MIEEE
	++y;
	*y++ = 0xffff;
	*y++ = 0xffff;
	*y++ = 0xffff;
#endif
#endif
	return;
	}
if( i == 0 )
	{
	(void )eshift( x, 4 );
	}
else
	{
	i <<= 4;
	(void )eshift( x, 5 );
	}
i |= *p++ & (unsigned short )0x0f;	/* *p = xi[M] */
*y |= (unsigned short )i; /* high order output already has sign bit set */
#ifdef IBMPC
*(--y) = *p++;
*(--y) = *p++;
*(--y) = *p;
#endif
#ifdef MIEEE
++y;
*y++ = *p++;
*y++ = *p++;
*y++ = *p++;
#endif
}

#endif /* not DEC */



/*
; e type to IEEE single precision
;	float d;
;	unsigned short x[N+2];
;	xtod( x, &d );
*/
void etoe24( x, e )
unsigned short *x, *e;
{
long exp;
unsigned short xi[NI];
int rndsav;

#ifdef NANS
if( eisnan(x) )
	{
	enan( e, 24 );
	return;
	}
#endif
emovi( x, xi );
exp = (long )xi[E] - (EXONE - 0177); /* adjust exponent for offsets */
#ifdef INFINITIES
if( eisinf(x) )
	goto nonorm;
#endif
/* round off to nearest or even */
rndsav = rndprc;
rndprc = 24;
emdnorm( xi, 0, 0, exp, 64 );
rndprc = rndsav;
nonorm:
toe24( xi, e );
}