/* 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 INFINITY 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.
*/
/* Check if e-type number is not a number.
*/ int eisnan(x) unsignedshort 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.
*/
/* Move in external format number, * converting it to internal format.
*/ void emovi( a, b ) unsignedshort *a, *b;
{ registerunsignedshort *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 INFINITY 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 ) unsignedshort *a, *b;
{ registerunsignedshort *p, *q; unsignedshort 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 INFINITY 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 ) registerunsignedshort *xi;
{ registerint i;
for( i=0; i<NI; i++ )
*xi++ = 0;
}
/* same, but don't touch the sign. */
void ecleazs( xi ) registerunsignedshort *xi;
{ registerint i;
++xi; for(i=0; i<NI-1; i++)
*xi++ = 0;
}
/* Move internal format number from a to b.
*/ void emovz( a, b ) registerunsignedshort *a, *b;
{ registerint i;
#ifdef INFINITY /* Return nonzero if internal format number is infinite. */
staticint
eiisinf (x) unsignedshort 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 ) registerunsignedshort *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);
p = &a[NI-2];
k = NBITS; while( *p == 0 ) /* significand is not supposed to be all zero */
{
eshdn6(a);
k -= 16;
} if( (*p & 0xff) == 0 )
{
eshdn8(a);
k -= 8;
}
q = &equot[NI-1];
j = 0; for( i=0; i<k; i++ )
{ if( *p & 1 )
eaddm(b, equot); /* remember if there were any nonzero bits shifted out */ if( *q & 1 )
j |= 1;
eshdn1(a);
eshdn1(equot);
}
for( i=0; i<NI; i++ )
b[i] = equot[i];
/* return flag for lost nonzero bits */ return(j);
}
#else
/* Multiply significand of e-type number b
by 16-bit quantity a, e-type result to c. */
void m16m( a, b, c ) unsignedshort a; unsignedshort b[], c[];
{ registerunsignedshort *pp; registerunsignedlong carry; unsignedshort *ps; unsignedshort p[NI]; unsignedlong aa, m; int i;
/* Do not execute the divide instruction if it will overflow. */ if( (tdenm * ((unsignedlong)0xffffL)) < tnum )
tquot = 0xffff; else
tquot = tnum / tdenm;
/* Prove that the divide worked. */ /* tcheck = (unsigned long )tquot * tdenm; if( tnum - tcheck > tdenm ) tquot = 0xffff;
*/ /* Multiply denominator by trial quotient digit. */
m16m( tquot, den, tprod ); /* The quotient digit may have been overestimated. */ if( ecmpm( tprod, num ) > 0 )
{
tquot -= 1;
esubm( den, tprod ); if( ecmpm( tprod, num ) > 0 )
{
tquot -= 1;
esubm( den, tprod );
}
} /* if( ecmpm( tprod, num ) > 0 ) { eshow( "tprod", tprod ); eshow( "num ", num ); printf( "tnum = %08lx, tden = %04x, tquot = %04x\n", tnum, den[M+1], tquot ); }
*/
esubm( tprod, num ); /* if( ecmpm( num, den ) >= 0 ) { eshow( "num ", num ); eshow( "den ", den ); printf( "tnum = %08lx, tden = %04x, tquot = %04x\n", tnum, den[M+1], tquot ); }
*/
equot[i] = tquot;
eshup6(num);
} /* test for nonzero remainder after roundoff bit */
p = &num[M];
j = 0; for( i=M; i<NI; i++ )
{
j |= *p++;
} if( j )
j = 1;
for( i=0; i<NI; i++ )
num[i] = equot[i];
return( (int )j );
}
/* Multiply significands */ int emulm( a, b ) unsignedshort a[], b[];
{ unsignedshort *p, *q; unsignedshort pprod[NI]; unsignedshort j; int i;
/* * Normalize and round off. * * The internal format number to be rounded is "s". * Input "lost" indicates whether the number is exact. * This is the so-called sticky bit. * * Input "subflg" indicates whether the number was obtained * by a subtraction operation. In that case if lost is nonzero * then the number is slightly smaller than indicated. * * Input "exp" is the biased exponent, which may be negative. * the exponent field of "s" is ignored but is replaced by * "exp" as adjusted by normalization and rounding. * * Input "rcntrl" is the rounding control.
*/
/* 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] = (unsignedshort )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 );
done:
emovo( bi, c );
}
/* ; Divide. ; ; unsigned short a[NE], b[NE], c[NE]; ; ediv( a, b, c ); c = b / a
*/ void ediv( a, b, c ) unsignedshort *a, *b, *c;
{ unsignedshort ai[NI], bi[NI]; int i, 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 /* 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 INFINITY if( eisinf(b) )
{
einfin(c); goto divsign;
} if( eisinf(a) )
{
eclear(c); goto divsign;
} #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); goto divsign;
}
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;
}
}
einfin(c);
mtherr( "ediv", SING ); goto divsign;
}
dnzro2:
i = edivm( ai, bi ); /* calculate exponent */
lt = ltb - lta + EXONE;
emdnorm( bi, i, 0, lt, 64 );
emovo( bi, c );
/* ; Multiply. ; ; unsigned short a[NE], b[NE], c[NE]; ; emul( a, b, c ); c = b * a
*/ void emul( a, b, c ) unsignedshort *a, *b, *c;
{ unsignedshort 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 INFINITY 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:
/* 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 ) unsignedshort *pe, *y;
{ #ifdef DEC
dectoe( pe, y ); /* see etodec.c */
#else
registerunsignedshort r; registerunsignedshort *p, *e; unsignedshort 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] -= (unsignedshort )(k-1);
}
emovo( yy, y ); #endif/* not DEC */
}
void e64toe( pe, y ) unsignedshort *pe, *y;
{ unsignedshort yy[NI]; unsignedshort *p, *q, *e; int i;
/* ; e type to IEEE single precision ; float d; ; unsigned short x[N+2]; ; xtod( x, &d );
*/ void etoe24( x, e ) unsignedshort *x, *e;
{ long exp; unsignedshort 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 INFINITY 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 );
}
#ifdef NANS if( eiisnan(x) )
{
enan( y, 24 ); return;
} #endif
p = &x[0]; #ifdef IBMPC
y += 1; #endif #ifdef DEC
y += 1; #endif
*y = 0; /* output high order */ if( *p++ )
*y = 0x8000; /* output sign bit */
i = *p++; if( i >= 255 )
{ /* Saturate at largest number less than infinity. */ #ifdef INFINITY
*y |= (unsignedshort )0x7f80; #ifdef IBMPC
*(--y) = 0; #endif #ifdef DEC
*(--y) = 0; #endif #ifdef MIEEE
++y;
*y = 0; #endif #else
*y |= (unsignedshort )0x7f7f; #ifdef IBMPC
*(--y) = 0xffff; #endif #ifdef DEC
*(--y) = 0xffff; #endif #ifdef MIEEE
++y;
*y = 0xffff; #endif #endif return;
} if( i == 0 )
{
(void )eshift( x, 7 );
} else
{
i <<= 7;
(void )eshift( x, 8 );
}
i |= *p++ & (unsignedshort )0x7f; /* *p = xi[M] */
*y |= i; /* high order output already has sign bit set */ #ifdef IBMPC
*(--y) = *p; #endif #ifdef DEC
*(--y) = *p; #endif #ifdef MIEEE
++y;
*y = *p; #endif
}
/* Compare two e type numbers. * * unsigned short a[NE], b[NE]; * 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.
*/ int ecmp( a, b ) unsignedshort *a, *b;
{ unsignedshort ai[NI], bi[NI]; registerunsignedshort *p, *q; registerint i; int msign;
#ifdef NANS if (eisnan (a) || eisnan (b)) return( -2 ); #endif
emovi( a, ai );
p = ai;
emovi( b, bi );
q = bi;
if( *(--p) > *(--q) ) return( msign ); /* p is bigger */ else return( -msign ); /* p is littler */
}
/* Find nearest integer to x = floor( x + 0.5 ) * * unsigned short x[NE], y[NE] * eround( x, y );
*/ void eround( x, y ) unsignedshort *x, *y;
{
eadd( ehalf, x, y );
efloor( y, y );
}
/* ; convert long (32-bit) integer to e type ; ; long l; ; unsigned short x[NE]; ; ltoe( &l, x ); ; note &l is the memory address of l
*/ void ltoe( lp, y ) long *lp; /* lp is the memory address of a long integer */ unsignedshort *y; /* y is the address of a short */
{ unsignedshort yi[NI]; unsignedlong ll; int k;
ecleaz( yi ); if( *lp < 0 )
{
ll = (unsignedlong )( -(*lp) ); /* make it positive */
yi[0] = 0xffff; /* put correct sign in the e type number */
} else
{
ll = (unsignedlong )( *lp );
} /* move the long integer to yi significand area */ if( sizeof(long) == 8 )
{
yi[M] = (unsignedshort) (ll >> (LONGBITS - 16));
yi[M + 1] = (unsignedshort) (ll >> (LONGBITS - 32));
yi[M + 2] = (unsignedshort) (ll >> 16);
yi[M + 3] = (unsignedshort) ll;
yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
} else
{
yi[M] = (unsignedshort )(ll >> 16);
yi[M+1] = (unsignedshort )ll;
yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
} if( (k = enormlz( yi )) > NBITS ) /* normalize the significand */
ecleaz( yi ); /* it was zero */ else
yi[E] -= (unsignedshort )k; /* subtract shift count from exponent */
emovo( yi, y ); /* output the answer */
}
/* ; convert unsigned long (32-bit) integer to e type ; ; unsigned long l; ; unsigned short x[NE]; ; ltox( &l, x ); ; note &l is the memory address of l
*/ void ultoe( lp, y ) unsignedlong *lp; /* lp is the memory address of a long integer */ unsignedshort *y; /* y is the address of a short */
{ unsignedshort yi[NI]; unsignedlong ll; int k;
ecleaz( yi );
ll = *lp;
/* move the long integer to ayi significand area */ if( sizeof(long) == 8 )
{
yi[M] = (unsignedshort) (ll >> (LONGBITS - 16));
yi[M + 1] = (unsignedshort) (ll >> (LONGBITS - 32));
yi[M + 2] = (unsignedshort) (ll >> 16);
yi[M + 3] = (unsignedshort) ll;
yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
} else
{
yi[M] = (unsignedshort )(ll >> 16);
yi[M+1] = (unsignedshort )ll;
yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
} if( (k = enormlz( yi )) > NBITS ) /* normalize the significand */
ecleaz( yi ); /* it was zero */ else
yi[E] -= (unsignedshort )k; /* subtract shift count from exponent */
emovo( yi, y ); /* output the answer */
}
/* ; Find long integer and fractional parts
; long i; ; unsigned short x[NE], frac[NE]; ; xifrac( x, &i, frac ); The integer output has the sign of the input. The fraction is the positive fractional part of abs(x).
*/ void eifrac( x, i, frac ) unsignedshort *x; long *i; unsignedshort *frac;
{ unsignedshort xi[NI]; int j, k; unsignedlong ll;
emovi( x, xi );
k = (int )xi[E] - (EXONE - 1); if( k <= 0 )
{ /* if exponent <= 0, integer = 0 and real output is fraction */
*i = 0L;
emovo( xi, frac ); return;
} if( k > (8 * sizeof(long) - 1) )
{ /* ; long integer overflow: output large integer ; and correct fraction
*/
j = 8 * sizeof(long) - 1; if( xi[0] )
*i = (long) ((unsignedlong) 1) << j; else
*i = (long) (((unsignedlong) (~(0L))) >> 1);
(void )eshift( xi, k );
} if( k > 16 )
{ /* Shift more than 16 bits: shift up k-16 mod 16 then shift by 16's.
*/
j = k - ((k >> 4) << 4);
eshift (xi, j);
ll = xi[M];
k -= j; do
{
eshup6 (xi);
ll = (ll << 16) | xi[M];
} while ((k -= 16) > 0);
*i = ll; if (xi[0])
*i = -(*i);
} else
{ /* shift not more than 16 bits */
eshift( xi, k );
*i = (long )xi[M] & 0xffff; if( xi[0] )
*i = -(*i);
}
xi[0] = 0;
xi[E] = EXONE - 1;
xi[M] = 0; if( (k = enormlz( xi )) > NBITS )
ecleaz( xi ); else
xi[E] -= (unsignedshort )k;
emovo( xi, frac );
}
/* ; Find unsigned long integer and fractional parts
; unsigned long i; ; unsigned short x[NE], frac[NE]; ; xifrac( x, &i, frac );
A negative e type input yields integer output = 0 but correct fraction.
*/ void euifrac( x, i, frac ) unsignedshort *x; unsignedlong *i; unsignedshort *frac;
{ unsignedshort xi[NI]; int j, k; unsignedlong ll;
emovi( x, xi );
k = (int )xi[E] - (EXONE - 1); if( k <= 0 )
{ /* if exponent <= 0, integer = 0 and argument is fraction */
*i = 0L;
emovo( xi, frac ); return;
} if( k > (8 * sizeof(long)) )
{ /* ; long integer overflow: output large integer ; and correct fraction
*/
*i = ~(0L);
(void )eshift( xi, k );
} elseif( k > 16 )
{ /* Shift more than 16 bits: shift up k-16 mod 16 then shift up by 16's.
*/
j = k - ((k >> 4) << 4);
eshift (xi, j);
ll = xi[M];
k -= j; do
{
eshup6 (xi);
ll = (ll << 16) | xi[M];
} while ((k -= 16) > 0);
*i = ll;
} else
{ /* shift not more than 16 bits */
eshift( xi, k );
*i = (long )xi[M] & 0xffff;
}
if( xi[0] ) /* A negative value yields unsigned integer 0. */
*i = 0L;
/* ; Shift significand ; ; Shifts significand area up or down by the number of bits ; given by the variable sc.
*/ int eshift( x, sc ) unsignedshort *x; int sc;
{ unsignedshort lost; unsignedshort *p;
/* ; normalize ; ; Shift normalizes the significand area pointed to by argument ; shift count (up = positive) is returned.
*/ int enormlz(x) unsignedshort x[];
{ registerunsignedshort *p; int sc;
sc = 0;
p = &x[M]; if( *p != 0 ) goto normdn;
++p; if( *p & 0x8000 ) return( 0 ); /* already normalized */ while( *p == 0 )
{
eshup6(x);
sc += 16; /* With guard word, there are NBITS+16 bits available. * return true if all are zero.
*/ if( sc > NBITS ) return( sc );
} /* see if high byte is zero */ while( (*p & 0xff00) == 0 )
{
eshup8(x);
sc += 8;
} /* now shift 1 bit at a time */ while( (*p & 0x8000) == 0)
{
eshup1(x);
sc += 1; if( sc > (NBITS+16) )
{
mtherr( "enormlz", UNDERFLOW ); return( sc );
}
} return( sc );
/* Normalize by shifting down out of the high guard word
of the significand */
normdn:
/* Test for exponent nonzero but significand denormalized. * This is an error condition.
*/ if( (y[NE-1] != 0) && ((y[NE-2] & 0x8000) == 0) )
{
mtherr( "etoasc", DOMAIN );
sprintf( string, "NaN" ); goto bxit;
}
/* Compare to 1.0 */
i = ecmp( eone, y ); if( i == 0 ) goto isone;
if( i < 0 )
{ /* Number is greater than 1 */ /* Convert significand to an integer and strip trailing decimal zeros. */
emov( y, u );
u[NE-1] = EXONE + NBITS - 1;
p = &etens[NTEN-4][0];
m = 16; do
{
ediv( p, u, t );
efloor( t, w ); 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] - (unsignedint )(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; unsignedshort *y;
{
asctoeg( s, y, 24 );
}
/* ASCII to double */ void asctoe53( s, y ) char *s; unsignedshort *y;
{ #ifdef DEC
asctoeg( s, y, 56 ); #else
asctoeg( s, y, 53 ); #endif
}
/* ASCII to long double */ void asctoe64( s, y ) char *s; unsignedshort *y;
{
asctoeg( s, y, 64 );
}
/* ASCII to 128-bit long double */ void asctoe113 (s, y) char *s; unsignedshort *y;
{
asctoeg( s, y, 113 );
}
/* ASCII to super double */ void asctoe( s, y ) char *s; unsignedshort *y;
{
asctoeg( s, y, NBITS );
}
/* Space to make a copy of the input string: */ staticchar lstr[82] = {0};
void asctoeg( ss, y, oprec ) char *ss; unsignedshort *y; int oprec;
{ unsignedshort yy[NI], xt[NI], tt[NI]; int esign, decflg, sgnflg, nexp, exp, prec, lost; int k, trail, c, rndsav; long lexp; unsignedshort nsign, *p; char *sp, *s;
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 );
*/ staticunsignedshort bmask[] = {
0xffff,
0xfffe,
0xfffc,
0xfff8,
0xfff0,
0xffe0,
0xffc0,
0xff80,
0xff00,
0xfe00,
0xfc00,
0xf800,
0xf000,
0xe000,
0xc000,
0x8000,
0x0000,
};
void efloor( x, y ) unsignedshort x[], y[];
{ registerunsignedshort *p; int e, expon, i; unsignedshort 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;
--> --------------------
--> maximum size reached
--> --------------------
Messung V0.5
¤ Dauer der Verarbeitung: 0.39 Sekunden
(vorverarbeitet)
¤
Die Informationen auf dieser Webseite wurden
nach bestem Wissen sorgfältig zusammengestellt. Es wird jedoch weder Vollständigkeit, noch Richtigkeit,
noch Qualität der bereit gestellten Informationen zugesichert.
Bemerkung:
Die farbliche Syntaxdarstellung und die Messung sind noch experimentell.
