📄 ieee.c
字号:
{
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 INFINITY
if( eisinf(b) )
{
if( eisneg(a) ^ eisneg(b) )
*(c+(NE-1)) = 0x8000;
else
*(c+(NE-1)) = 0;
einfin(c);
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);
mtherr( "ediv", SING );
return;
}
dnzro2:
i = edivm( ai, bi );
/* calculate exponent */
lt = ltb - lta + EXONE;
emdnorm( bi, i, 0, lt, 64 );
/* set the sign */
if( ai[0] == bi[0] )
bi[0] = 0;
else
bi[0] = 0Xffff;
emovo( bi, c );
}
/*
; 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;
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 INFINITY
if( eisinf(a) || eisinf(b) )
{
if( eisneg(a) ^ eisneg(b) )
*(c+(NE-1)) = 0x8000;
else
*(c+(NE-1)) = 0;
einfin(c);
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 );
/* calculate exponent */
lt = lta + ltb - (EXONE - 1);
emdnorm( bi, j, 0, lt, 64 );
/* calculate sign of product */
if( ai[0] == bi[0] )
bi[0] = 0;
else
bi[0] = 0xffff;
emovo( bi, c );
}
/*
; 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 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
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 INFINITY
/* Point to the exponent field. */
p = &yy[NE-1];
if( *p == 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
for( i=1; i<=4; i++ )
{
if( pe[i] != 0 )
{
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 INFINITY
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 /* INFINITY */
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 INFINITY
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 INFINITY
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 INFINITY
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 INFINITY
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 INFINITY
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;
{
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -