/* qflt.c * QFLOAT * * Extended precision floating point routines * * asctoq( string, q ) ascii string to q type * dtoq( &d, q ) DEC double precision to q type * etoq( &d, q ) IEEE double precision to q type * e24toq( &d, q ) IEEE single precision to q type * e113toq( &d, q ) 128-bit long double precision to q type * ltoq( &l, q ) long integer to q type * qabs(q) absolute value * qadd( a, b, c ) c = b + a * qclear(q) q = 0 * qcmp( a, b ) compare a to b * qdiv( a, b, c ) c = b / a * qifrac( x, &l, frac ) x to integer part l and q type fraction * qfrexp( x, l, y ) find exponent l and fraction y between .5 and 1 * qldexp( x, l, y ) multiply x by 2^l * qinfin( x ) set x to infinity, leaving its sign alone * qmov( a, b ) b = a * qmul( a, b, c ) c = b * a * qmuli( a, b, c ) c = b * a, a has only 16 significant bits * qisneg(q) returns sign of q * qneg(q) q = -q * qnrmlz(q) adjust exponent and mantissa * qsub( a, b, c ) c = b - a * qtoasc( a, s, n ) q to ASCII string, n digits after decimal * qtod( q, &d ) convert q type to DEC double precision * qtoe( q, &d ) convert q type to IEEE double precision * qtoe24( q, &d ) convert q type to IEEE single precision * qtoe113( q, &d ) convert q type to 128-bit long double precision * * Data structure of the number (a "word" is 16 bits) * * sign word (0 for positive, -1 for negative) * exponent (EXPONE for 1.0) * high guard word (always zero after normalization) * N-1 mantissa words (most significant word first, * most significant bit is set) * * Numbers are stored in C language as arrays. All routines * use pointers to the arrays as arguments. * * The result is always normalized after each arithmetic operation. * All arithmetic results are chopped. No rounding is performed except * on conversion to double precision. */ /* * Revision history: * * SLM, 5 Jan 84 PDP-11 assembly language version * SLM, 2 Mar 86 fixed bug in asctoq() * SLM, 6 Dec 86 C language version * */ #include #include "mconf.h" #include "qhead.h" #ifdef UNK #if BIGENDIAN #define MIEEE #else #define IBMPC #endif #undef UNK #endif /* number of 16 bit words in mantissa area */ /*define N 10*/ /* max number of decimal digits in conversion */ /*define NDEC 46*/ /* number of bits of precision */ /*define NBITS (OMG-1)*16*/ /* array index of the sign (a whole int wasted) */ #define SIGNWORD 0 /* array index of the exponent */ #define E 1 /* array index of first mantissa word */ #define M 2 /* Define 1 for sticky bit rounding */ #ifndef STICKY #define STICKY 0 #endif /* accumulators */ static QELT ac1[NQ+1]; static QELT ac2[NQ+1]; static QELT ac3[NQ+1]; static QELT ac4[NQ+1]; /* shift count register */ static int SC; extern QELT qzero[], qone[]; #ifdef ANSIPROT int shdn1( QELT * ), shup1( QELT * ), shdn8( QELT * ), shup8( QELT * ); int shdn16( QELT * ), shup16( QELT * ); int mulm( QELT *, QELT * ), mulin( QELT *, QELT * ), divn( QELT *, QELT * ); int qmovz( QELT *, QELT * ), normlz( QELT * ), cmpm( QELT *, QELT * ); int addm( QELT *, QELT * ), subm( QELT *, QELT * ), shift( QELT * ); int qadd1( QELT *, QELT *, QELT * ); #else int shdn1(), shup1(), shdn8(), shup8(), shdn16(), shup16(); int mulm(), mulin(), divm(), qmovz(), normlz(); int cmpm(), addm(), subm(), qadd1(), shift(); #endif /* ; Absolute value ; ; QELT q[NQ]; ; qabs( q ); */ int qabs(x) QELT x[]; { x[SIGNWORD] = 0; return 0; } /* Test if negative */ int qisneg(x) QELT *x; { return (x[SIGNWORD] != 0); } /* ; Negate ; ; short q[NQ]; ; qneg( q ); */ int qneg(x) QELT *x; { *x = ~(*x); /* complement the sign */ return 0; } /* ; convert long (32-bit) integer to q type ; ; long l; ; QELT q[NQ]; ; ltoq( &l, q ); ; note &l is the memory address of l */ int ltoq( lp, y ) long *lp; QELT y[]; { register QELT *p; /* use processor registers if possible */ long ll; qclear( ac1 ); ll = *lp; /* local copy of lp */ if( ll < 0 ) { ll = -ll; /* make it positive */ ac1[SIGNWORD] = -1; /* put correct sign in the q type number */ } /* move the long integer to ac1 mantissa area */ p = &ac1[M]; *p++ = ll >> 16; *p = ll; /* * q = (short *)≪ * #if DEC * p = &ac1[M]; * *p++ = *q++; * #endif * #if IBMPC * p = &ac1[M+1]; * *p-- = *q++; * #endif * *p = *q; */ ac1[E] = EXPONE+15; /* exponent if normalize shift count were 0 */ ac1[NQ] = 0; if( normlz( ac1 ) ) /* normalize the mantissa */ qclear( ac1 ); /* it was zero */ else ac1[E] -= SC; /* else subtract shift count from exponent */ qmov( ac1, y ); /* output the answer */ return 0; } /* ; convert DEC double precision to q type ; double d; ; QELT q[NQ]; ; dtoq( &d, q ); */ int dtoq( d, y ) unsigned short *d; QELT *y; { register QELT r, *p; qclear(y); /* start with a zero */ p = y; /* point to our number */ r = *d; /* get DEC exponent word */ if( r & 0x8000 ) *p = -1; /* fill in our sign */ ++p; /* bump pointer to our exponent word */ r &= 0x7fff; /* strip the sign bit */ if( r == 0 ) /* answer = 0 if high order DEC word = 0 */ return 0; r >>= 7; /* shift exponent word down 7 bits */ r -= 0200; /* subtract DEC exponent offset */ r += EXPONE - 1; /* add our q type exponent offset */ *p++ = r; /* to form our exponent */ r = *d++; /* now do the high order mantissa */ r &= 0177; /* strip off the DEC exponent and sign bits */ r |= 0200; /* the DEC understood high order mantissa bit */ *p++ = r; /* put result in our high guard word */ *p++ = *d++; /* fill in the rest of our mantissa */ *p++ = *d++; *p = *d; shdn8(y); /* shift our mantissa down 8 bits */ return 0; } /* ; convert q type to DEC double precision ; double d; ; QELT q[NQ]; ; qtod( q, &d ); */ int qtod( x, d ) QELT *x; unsigned short *d; { register int r; int i, j; *d = 0; if( *x != 0 ) *d = 0100000; qmovz( x, ac1 ); r = ac1[E]; if( r < (EXPONE - 1 - 0200) ) goto zout; i = ac1[M+4]; if( (i & 0200) != 0 ) { if( (i & 0377) == 0200 ) { if( (i & 0400) != 0 ) { /* check all less significant bits */ for( j=M+5; j<=NQ; j++ ) { if( ac1[j] != 0 ) goto yesrnd; } } goto nornd; } yesrnd: qclear( ac2 ); ac2[ M+4 ] = 0200; ac2[NQ] = 0; addm( ac2, ac1 ); normlz(ac1); r -= SC; } nornd: r -= EXPONE - 1; r += 0200; if( r < 0 ) { zout: *d++ = 0; *d++ = 0; *d++ = 0; *d++ = 0; return 0; } if( r >= 0377 ) { *d++ = 077777; *d++ = -1; *d++ = -1; *d++ = -1; return 0; } r &= 0377; r <<= 7; shup8( ac1 ); ac1[M] &= 0177; r |= ac1[M]; *d++ |= r; *d++ = ac1[M+1]; *d++ = ac1[M+2]; *d++ = ac1[M+3]; return 0; } /* ; Find integer and fractional parts ; long i; ; QELT x[NQ], frac[NQ]; ; qifrac( x, &i, frac ); */ int qifrac( x, i, frac ) QELT x[]; long *i; QELT frac[]; { qmovz( x, ac1 ); SC = ac1[E] - (EXPONE - 1); if( SC <= 0 ) { /* if exponent <= 0, integer = 0 and argument is fraction */ *i = 0L; qmov( ac1, frac ); return 0; } if( SC > 31 ) { /* ; long integer overflow: output large integer ; and correct fraction */ *i = 0x7fffffff; shift( ac1 ); goto lab10; } if( SC > 16 ) { /* ; shift more than 16 bits: shift up SC-16, output the integer, ; then complete the shift to get the fraction. */ SC -= 16; shift( ac1 ); *i = ((unsigned long )ac1[M] << 16) | (unsigned short )ac1[M+1]; /* * p = (short *)i; * #ifdef DEC * *p++ = ac1[M]; * *p++ = ac1[M+1]; * #else * *p++ = ac1[M+1]; * *p++ = ac1[M]; * #endif */ shup16( ac1 ); goto lab10; } /* shift not more than 16 bits */ shift( ac1 ); *i = ac1[M] & 0xffff; lab10: if( x[SIGNWORD] ) *i = -(*i); ac1[SIGNWORD] = 0; ac1[E] = EXPONE - 1; ac1[M] = 0; if( normlz(ac1) ) qclear( ac1 ); else ac1[E] -= SC; qmov( ac1, frac ); return 0; } int qldexp( x, n, y) QELT *x; long n; QELT *y; { long k; k = (long) x[E] + n; qmov( x, y ); y[E] = k; if( (k > MAXEXP) || (n > (2 * (long)MAXEXP)) ) qinfin(y); if( (k <= 0) || (n < (-2 * (long)MAXEXP)) ) qclear(y); return(0); } int qfrexp (x, e, y) QELT *x; long *e; QELT *y; { if( x[E] == 0 ) { *e = 0; qclear(y); } else { *e = (long) x[E] - (long) EXPONE + 1; qmov( x, y ); y[1] = EXPONE - 1; } return 0; } /* ; subtract ; ; QELT a[NQ], b[NQ], c[NQ]; ; qsub( a, b, c ); c = b - a */ static short subflg = 0; int qsub( a, b, c ) QELT *a, *b, *c; { subflg = 1; qadd1( a, b, c ); return 0; } /* ; add ; ; QELT a[NQ], b[NQ], c[NQ]; ; qadd( a, b, c ); c = b + a */ int qadd( a, b, c ) QELT *a, *b, *c; { subflg = 0; qadd1( a, b, c ); return 0; } int qadd1( a, b, c ) QELT *a, *b, *c; { long lt; int i; #if STICKY int lost; #endif qmovz( a, ac1 ); qmovz( b, ac2 ); if( subflg ) ac1[SIGNWORD] = ~ac1[SIGNWORD]; /* compare exponents */ lt = (long) ac1[E] - (long) ac2[E]; if( lt > 0 ) { /* put the larger number in ac2 */ qmovz( ac2, ac3 ); qmov( ac1, ac2 ); qmov( ac3, ac1 ); lt = -lt; } SC = lt; #if STICKY lost = 0; #endif if( lt != 0 ) { if( lt < -NBITS-1 ) goto done; /* answer same as larger addend */ #if STICKY lost = shift( ac1 ); /* shift the smaller number down */ #else shift( ac1 ); /* shift the smaller number down */ #endif } else { /* exponents were the same, so must compare mantissae */ i = cmpm( ac1, ac2 ); if( i == 0 ) { /* the numbers are identical */ /* if different signs, result is zero */ if( ac1[SIGNWORD] != ac2[SIGNWORD] ) goto underf; /* if exponents zero, result is zero */ if( ac1[E] == 0 ) goto underf; /* if same sign, result is double */ if( ac1[E] >= MAXEXP ) { qclear(c); if( ac1[SIGNWORD] != 0 ) qneg(c); goto overf; } ac2[E] += 1; goto done; } if( i > 0 ) { /* put the larger number in ac2 */ qmovz( ac2, ac3 ); qmov( ac1, ac2 ); qmov( ac3, ac1 ); } } if( ac1[SIGNWORD] == ac2[SIGNWORD] ) { addm( ac1, ac2 ); subflg = 0; } else { subm( ac1, ac2 ); subflg = 1; } if( normlz(ac2) ) goto underf; lt = (long )ac2[E] - SC; if( lt > MAXEXP ) goto overf; if( lt < 0 ) { /* mtherr( "qadd", UNDERFLOW );*/ goto underf; } ac2[E] = lt; /* round off */ i = ac2[NQ] & 0xffff; if( i & SIGNBIT ) { #if STICKY if( i == SIGNBIT ) { if( lost == 0 ) { /* Critical case, round to even */ if( (ac2[NQ-1] & 1) == 0 ) goto done; } else { if( subflg != 0 ) goto done; } } #else if( subflg != 0 ) goto done; #endif qclear( ac1 ); ac1[NQ] = 0; ac1[NQ-1] = 1; addm( ac1, ac2 ); normlz(ac2); if( SC ) { lt = (long )ac2[E] - SC; if( lt > MAXEXP ) goto overf; ac2[E] = lt; } } done: qmov( ac2, c ); return 0; underf: qclear(c); return 0; overf: mtherr( "qadd", OVERFLOW ); qinfin(c); return 0; } /* ; divide ; ; QELT a[NQ], b[NQ], c[NQ]; ; qdiv( a, b, c ); c = b / a */ /* for Newton iteration version: * extern short qtwo[]; * static short qt[NQ] = {0}; * static short qu[NQ] = {0}; */ int qdiv( a, b, c ) QELT *a, *b, *c; { long lt; if( b[E] == 0 ) { divunderf: qclear(c); /* numerator is zero */ return 0; } if( a[E] == 0 ) { /* divide by zero */ qinfin(c); mtherr( "qdiv", SING ); return 0; } qmovz( b, ac3 ); /* Avoid exponent underflow in mdnorm. */ lt = (long) ac3[E]; ac3[E] = 4; divm( a, ac3 ); /* calculate exponent */ lt = lt + (long )ac3[E] -4L - (long )a[E]; ac3[E] = lt; ac3[NQ] = 0; normlz(ac3); lt = lt - SC + EXPONE + 1; if( lt > MAXEXP ) { qinfin(ac3); mtherr( "qdiv", OVERFLOW ); } else if( lt <= 0 ) goto divunderf; else ac3[E] = lt; if( a[SIGNWORD] == b[SIGNWORD] ) ac3[SIGNWORD] = 0; else ac3[SIGNWORD] = -1; qmov( ac3, c ); return 0; } /* ; multiply ; ; QELT a[NQ], b[NQ], c[NQ]; ; qmul( a, b, c ); c = b * a */ int qmul( a, b, c ) QELT *a, *b, *c; { QELT *p; register int ctr; long lt; if( (a[E] == 0) || (b[E] == 0) ) { qclear(c); return 0; } /* detect multiplication by small integer a */ if( a[M+2] == 0 ) { p = &a[M+3]; for( ctr=M+3; ctr MAXEXP ) goto overf; ac3[E] = lt; ac3[NQ] = 0; if( normlz(ac3) ) goto underf; lt = lt - SC + EXPONE -1; if( lt > MAXEXP ) goto overf; if( lt < 0 ) goto underf; ac3[E] = lt; qmov( ac3, c ); return 0; underf: qclear(c); return 0; overf: qinfin(c); mtherr( "qmul", OVERFLOW ); return 0; } /* Multiply, a has at most WORDSIZE significant bits */ int qmuli( a, b, c ) QELT *a, *b, *c; { long lt; if( (a[E] == 0) || (b[E] == 0) ) { qclear(c); return 0; } qmov( b, ac3 ); mulin( a, ac3 ); /* calculate sign of product */ if( b[SIGNWORD] == a[SIGNWORD] ) ac3[SIGNWORD] = 0; else ac3[SIGNWORD] = -1; /* calculate exponent */ lt = ((long)ac3[E] - (EXPONE-1)) + ((long )a[E] - (EXPONE - 1)); if( lt > MAXEXP ) goto overf; ac3[E] = lt; ac3[NQ] = 0; if( normlz(ac3) ) goto underf; lt = lt - SC + EXPONE - 1; if( lt > MAXEXP ) goto overf; if( lt < 0 ) goto underf; ac3[E] = lt; qmov( ac3, c ); return 0; underf: qclear(c); return 0; overf: qinfin(c); mtherr( "qmuli", OVERFLOW ); return 0; } /* ; Compare mantissas ; ; QELT a[NQ], b[NQ]; ; cmpm( a, b ); ; ; for the mantissas: ; returns +1 if a > b ; 0 if a == b ; -1 if a < b */ int cmpm( a, b ) register QELT *a, *b; { int i; a += M; /* skip up to mantissa area */ b += M; for( i=0; i (unsigned int) *(--b) ) return(1); else return(-1); } /* ; shift mantissa ; ; Shifts mantissa area up or down by the number of bits ; given by the variable SC. */ int shift( x ) QELT *x; { QELT *p; #if STICKY int lost; #endif if( SC == 0 ) return(0); #if STICKY lost = 0; #endif if( SC < 0 ) { p = x + NQ; SC = -SC; while( SC >= 16 ) { #if STICKY lost |= *p; #endif shdn16(x); SC -= 16; } while( SC >= 8 ) { #if STICKY lost |= *p & 0xff; #endif shdn8(x); SC -= 8; } while( SC > 0 ) { #if STICKY lost |= *p & 1; #endif shdn1(x); SC -= 1; } } else { while( SC >= 16 ) { shup16(x); SC -= 16; } while( SC >= 8 ) { shup8(x); SC -= 8; } while( SC > 0 ) { shup1(x); SC -= 1; } } #if STICKY return( lost ); #else return(0); #endif } /* ; normalize ; ; shift normalizes the mantissa area pointed to by R1 ; shift count (up = positive) returned in SC */ int normlz(x) QELT x[]; { register QELT *p; SC = 0; p = &x[M]; if( *p != 0 ) goto normdn; ++p; if( *p & SIGNBIT ) return(0); /* already normalized */ while( *p == 0 ) { shup16(x); SC += 16; /* With guard word, there are NBITS+WORDSIZE bits available. * return true if all are zero. */ if( SC > NBITS ) return(1); } /* see if high byte is zero */ while( (*p & 0xff00) == 0 ) { shup8(x); SC += 8; } /* now shift 1 bit at a time */ while( (*p & SIGNBIT) == 0) { shup1(x); SC += 1; /* if( SC > NBITS ) { printf( "normlz error\n"); return(0); } */ } return(0); /* normalize by shifting down out of the high guard word of the mantissa */ normdn: if( *p & 0xff00 ) { shdn8(x); SC -= 8; } while( *p != 0 ) { shdn1(x); SC -= 1; /* if( SC < -NBITS ) { printf( "low normlz error\n"); return(0); } */ } return(0); } /* ; Clear out entire number, including sign and exponent, pointed ; to by x ; ; QELT x[]; ; qclear( x ); */ /* Moved to qfltb.c */ #if 0 int qclear( x ) register QELT *x; { register int i; for( i=0; i>= 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 += EXPONE - 0177; yy[E] = r; p = &yy[M+1]; #ifdef IBMPC m = *(--e); #endif #ifdef DEC m = *(--e); #endif #ifdef MIEEE ++e; m = *e++; #endif #if WORDSIZE == 32 m <<= 16; #endif *p++ = m; shdn8( yy ); if( denorm ) { /* if zero exponent, then normalize the significand */ normlz(yy); if( SC > NBITS ) qclear(yy); else yy[E] -= (QELT) (SC-1); } qmov( yy, y ); return 0; } /* ; e type to IEEE single precision ; float d; ; QELT x[NQ]; ; xtod( x, &d ); */ int qtoe24( x, e ) QELT *x; unsigned short *e; { long exp; QELT xi[NQ+1]; QELT *p; QELT k; unsigned long f; int j, sign; qmov( x, xi ); xi[NQ] = 0; p = &xi[0]; f = 0L; if( *p++ ) sign = 0x8000; /* output sign bit */ else sign = 0; exp = (long )*p++ - (EXPONE - 0177); /* adjust exponent for offsets */ if( exp <= 0 ) { if( exp > -24 ) { SC = exp - 1; shift( xi ); exp = 0; } else { f = 0; goto fout; } } /* round off to nearest or even */ #if WORDSIZE == 32 k = xi[M+1]; #else k = xi[M+2]; #endif if( (k & 0x80) != 0 ) { if( (k & 0x0ff) == 0x80 ) { if( (k & 0x100) != 0 ) { /* check all less significant bits */ #if WORDSIZE == 32 for( j=M+2; j<=NQ; j++ ) #else for( j=M+3; j<=NQ; j++ ) #endif { if( xi[j] != 0 ) goto yesrnd; } } goto nornd; } yesrnd: qclear( ac2 ); ac2[NQ] = 0; #if WORDSIZE == 32 ac2[M+1] = 0x80; #else ac2[M+2] = 0x80; #endif addm( ac2, xi ); if( xi[2] ) { shdn1(xi); exp += 1; } if( (exp == 0) && (xi[M+1] & SIGNBIT) ) exp += 1; } nornd: if( exp >= 255 ) { /* Saturate at largest number less than infinity. */ #ifdef INFINITY f = (unsigned long )0x7f800000; #else f = (unsigned long )0x7f7fffff; #endif goto fout; } f = *p++; /* Skip over high guard word. */ f = *p++; #if WORDSIZE == 16 f = (f << 16) | *p; #endif f &= 0x7fffffff; /* f += 0x80; */ /* Retire old DEC rounding. */ f >>= 8; f |= (exp << 23) & 0x7f800000; fout: if(sign) f |= 0x80000000L; #ifdef IBMPC *e++ = f; *e = f >> 16; #endif #ifdef DEC *e++ = f; *e = f >> 16; #endif #ifdef MIEEE *e++ = f >> 16; *e = f; #endif return 0; } /* ; Convert IEEE double precision to Q type ; double d; ; QELT q[NQ]; ; etoq( &d, q ); */ int etoq( e, y ) unsigned short *e; QELT *y; { #ifdef DEC dtoq(e,y); #else register int r; register QELT *p; QELT yy[NQ+1]; int denorm; denorm = 0; /* flag if denormalized number */ qclear(yy); yy[NQ] = 0; #ifdef MIEEE #endif #ifdef IBMPC e += M+1; #endif /* r = *e & 0x7fff; if( r == 0 ) return 0; */ r = *e; yy[SIGNWORD] = 0; if( r & 0x8000 ) yy[SIGNWORD] = -1; yy[M] = (r & 0x0f) | 0x10; r &= ~0x800f; /* strip sign and 4 mantissa bits */ r >>= 4; /* If zero exponent, then the mantissa is denormalized. * So take back the understood high mantissa bit. */ if( r == 0 ) { denorm = 1; yy[M] &= ~0x10; } r += EXPONE - 01777; yy[E] = r; p = &yy[M+1]; #ifdef MIEEE ++e; *p++ = *e++; *p++ = *e++; *p++ = *e++; #endif #ifdef IBMPC *p++ = *(--e); *p++ = *(--e); *p++ = *(--e); #endif SC = -5; shift(yy); if( denorm ) { /* if zero exponent, then normalize the mantissa */ if( normlz( yy ) ) qclear(yy); else yy[E] -= SC-1; } qmov( yy, y ); /* not DEC */ #endif return 0; } /* ; Q type to IEEE double precision ; double d; ; short q[NQ]; ; qtoe( q, &d ); */ int qtoe( x, e ) QELT *x; unsigned short *e; { #ifdef DEC qtod(x,e); #else int j, k; long i; register QELT *p; #ifdef MIEEE #endif #ifdef IBMPC e += M+1; #endif *e = 0; /* output high order */ p = &ac1[SIGNWORD]; qmovz( x, ac1 ); if( *p++ != 0 ) *e = 0x8000; /* output sign bit */ if( normlz(ac1) ) goto ozero; *p -= SC; i = (long) *p++ - (EXPONE - 1023); /* adjust exponent for offsets */ /* Handle denormalized small numbers. */ if( i <= 0 ) { if( i > -53 ) { SC = i - 1; shift( ac1 ); i = 0; } else { ozero: #ifdef MIEEE ++e; *e++ = 0; *e++ = 0; *e++ = 0; #endif #ifdef IBMPC *(--e) = 0; *(--e) = 0; *(--e) = 0; #endif return 0; } } /* round off to nearest or even */ k = ac1[M+4]; if( (k & 0x400) != 0 ) { if( (k & 0x07ff) == 0x400 ) { if( (k & 0x800) != 0 ) { /* check all less significant bits */ for( j=M+5; j<=NQ; j++ ) { if( ac1[j] != 0 ) goto yesrnd; } } goto nornd; } yesrnd: qclear( ac2 ); ac2[NQ] = 0; ac2[M+4] = 0x800; addm( ac2, ac1 ); if( ac1[2] ) { shdn1(ac1); i += 1; } if( (i == 0) && (ac1[3] & SIGNBIT) ) i += 1; } nornd: if( i > 2047 ) { /* Saturate at largest number less than infinity. */ mtherr( "qtoe", OVERFLOW ); *e |= 0x7fef; #ifdef MIEEE ++e; *e++ = 0xffff; *e++ = 0xffff; *e++ = 0xffff; #endif #ifdef IBMPC *(--e) = 0xffff; *(--e) = 0xffff; *(--e) = 0xffff; #endif return 0; } i <<= 4; SC = 5; shift( ac1 ); i |= *p++ & 0x0f; /* *p = ac1[M] */ *e |= i; /* high order output already has sign bit set */ #ifdef MIEEE ++e; *e++ = *p++; *e++ = *p++; *e++ = *p++; #endif #ifdef IBMPC *(--e) = *p++; *(--e) = *p++; *(--e) = *p; #endif /* not DEC */ #endif return 0; } /* ; Convert 80-bit IEEE long double precision to Q type ; long double d; ; QELT q[N+2]; ; etoq( &d, q ); */ int e64toq( e, y ) unsigned short *e; QELT *y; { register int r; register QELT *p; QELT yy[NQ+1]; int denorm; denorm = 0; /* flag if denormalized number */ qclear(yy); yy[NQ] = 0; #ifdef MIEEE #endif #ifdef IBMPC e += 4; #endif r = *e; yy[SIGNWORD] = 0; if( r & 0x8000 ) yy[SIGNWORD] = -1; r &= 0x7fff; /* If zero exponent, then the mantissa is denormalized. */ if( r == 0 ) { denorm = 1; } r += (EXPONE-0x3fff); yy[E] = r; p = &yy[M+1]; #if WORDSIZE== 32 #ifdef MIEEE ++e; *p = ((unsigned int) *e++) << 16; *p++ |= *e++; *p = ((unsigned int) *e++) << 16; *p++ |= *e++; #endif #ifdef IBMPC *p = ((unsigned int) *(--e)) << 16; *p++ |= *(--e); *p = ((unsigned int) *(--e)) << 16; *p++ |= *(--e); #endif #else /* 16-bit wordsize */ #ifdef MIEEE ++e; *p++ = *e++; *p++ = *e++; *p++ = *e++; *p++ = *e++; #endif #ifdef IBMPC *p++ = *(--e); *p++ = *(--e); *p++ = *(--e); *p++ = *(--e); #endif #endif if( denorm ) { /* if zero exponent, then normalize the mantissa */ #ifdef IBMPC /* For Intel long double, shift denormal significand up 1 -- but only if the top significand bit is zero. */ if( (yy[M+1] & SIGNBIT) == 0 ) shup1( yy ); #endif if( normlz( yy ) ) qclear(yy); else yy[E] -= SC; } qmov( yy, y ); return 0; } /* ; Q type to 80-bit IEEE long double precision ; long double d; ; QELT q[N+2]; ; qtoe( q, &d ); */ int qtoe64( x, e ) QELT *x; unsigned short *e; { QELT k; int j; long i; register QELT *p; #ifdef MIEEE #endif #ifdef IBMPC e += 4; #if 0 /* NOTE: if data type is 96 bits wide, clear the last word here. */ *(e+1)= 0; #endif #endif *e = 0; /* output high order */ p = &ac1[0]; qmovz( x, p ); if( ac1[SIGNWORD] ) *e = 0x8000; /* output sign bit */ ++p; if( normlz(ac1) ) goto o64zero; *p -= SC; i = (long) *p++ - (EXPONE-0x3fff); /* adjust exponent for offsets */ if( i >= 0x7fff) { #ifdef MIEEE *e++ |= 0x7ffe; *e++ = 0xffff; *e++ = 0xffff; *e++ = 0xffff; *e++ = 0xffff; #endif #ifdef IBMPC *e-- |= 0x7ffe; *e-- = 0xffff; *e-- = 0xffff; *e-- = 0xffff; *e-- = 0xffff; #endif return 0; } /* We can't handle denormal numbers if the q-type exponenent is 0. */ if( i <= 0 ) { if( i > -65 ) { #if BIGENDIAN SC = i; #else /* Intel 80x87 loses a bit. */ SC = i - 1; #endif shift( ac1 ); i = 0; } else { o64zero: #ifdef MIEEE ++e; *e++ = 0; *e++ = 0; *e++ = 0; *e++ = 0; #endif #ifdef IBMPC *(--e) = 0; *(--e) = 0; *(--e) = 0; *(--e) = 0; #endif return 0; } } /* round off to nearest or even */ #if WORDSIZE == 32 k = ac1[M+3]; #else k = ac1[M+5]; #endif if( (k & SIGNBIT) != 0 ) { if( (k & ((QELT) SIGNBIT - 1)) == 0 ) { /* check all less significant bits */ #if WORDSIZE == 32 /* 0 1 2 3 4 5 6 */ /* S E M 1,2 3,4 5,6 7,8 */ for( j=M+4; j<=NQ; j++ ) #else /* 0 1 2 3 4 5 6 7 8 9 10 */ /* S E M 1 2 3 4 5 6 7 8 */ for( j=M+6; j<=NQ; j++ ) #endif { if( ac1[j] != 0 ) goto yesrnd; } /* round to even */ #if WORDSIZE == 32 if( (ac1[4] & 1) == 0 ) #else if( (ac1[6] & 1) == 0 ) #endif goto nornd; } yesrnd: qclear( ac2 ); ac2[NQ] = 0; #if WORDSIZE == 32 ac2[M+3] = SIGNBIT; #else ac2[M+5] = SIGNBIT; #endif addm( ac2, ac1 ); if( ac1[M] ) { shdn1(ac1); i += 1; } if( (i == 0) && (ac1[3] & SIGNBIT) ) i += 1; } nornd: *e |= i; /* high order output already has sign bit set */ p = &ac1[M+1]; #if WORDSIZE == 32 #ifdef MIEEE ++e; *e++ = *p >> 16; *e++ = *p++; *e++ = *p >> 16; *e++ = *p++; #endif #ifdef IBMPC *(--e) = *p >> 16; *(--e) = *p++; *(--e) = *p >> 16; *(--e) = *p++; #endif #else /* 16-bit words */ #ifdef MIEEE ++e; *e++ = *p++; *e++ = *p++; *e++ = *p++; *e++ = *p; #endif #ifdef IBMPC *(--e) = *p++; *(--e) = *p++; *(--e) = *p++; *(--e) = *p; #endif #endif return 0; } /* ; Convert 128-bit IEEE long double precision to Q type ; long double d; ; QELT q[NQ]; ; e113toq( &d, q ); */ int e113toq( e, y ) unsigned short *e; QELT *y; { register short r; register QELT *p; short yy[NQ+1]; int denorm; denorm = 0; /* flag if denormalized number */ qclear(yy); yy[NQ] = 0; #ifdef MIEEE #endif #ifdef IBMPC e += 7; #endif r = *e; yy[0] = 0; if( r & 0x8000 ) yy[0] = -1; r &= 0x7fff; yy[M] = 1; /* If zero exponent, then the mantissa is denormalized. * So, take back the understood high mantissa bit. */ if( r == 0 ) { denorm = 1; yy[M] = 0; } r += (EXPONE-0x3fff); yy[E] = r; p = &yy[M+1]; #ifdef MIEEE ++e; *p++ = *e++; *p++ = *e++; *p++ = *e++; *p++ = *e++; *p++ = *e++; *p++ = *e++; *p++ = *e++; #endif #ifdef IBMPC *p++ = *(--e); *p++ = *(--e); *p++ = *(--e); *p++ = *(--e); *p++ = *(--e); *p++ = *(--e); *p++ = *(--e); #endif SC = -1; shift(yy); if( denorm ) { /* if zero exponent, then normalize the mantissa */ if( normlz( yy ) ) qclear(yy); else yy[E] -= SC-1; } qmov( yy, y ); return 0; } /* ; Q type to 128-bit IEEE long double precision ; long double d; ; QELT q[NQ]; ; qtoe113( q, &d ); */ int qtoe113( x, e ) QELT *x; unsigned short *e; { QELT k; short i, j; register QELT *p; #ifdef MIEEE #endif #ifdef IBMPC e += 7; #endif *e = 0; /* output high order */ p = &ac1[0]; qmovz( x, p ); if( *p++ != 0 ) *e = 0x8000; /* output sign bit */ normlz(ac1); *p -= SC; i = *p++ - (EXPONE-0x3fff); /* adjust exponent for offsets */ /* We can't handle denormal numbers. */ if( i <= 0 ) { if( i > -2 ) { SC = i - 1; shift( ac1 ); i = 0; } else { /*ozero:*/ #ifdef MIEEE ++e; *e++ = 0; *e++ = 0; *e++ = 0; *e++ = 0; *e++ = 0; *e++ = 0; *e++ = 0; #endif #ifdef IBMPC *(--e) = 0; *(--e) = 0; *(--e) = 0; *(--e) = 0; *(--e) = 0; *(--e) = 0; *(--e) = 0; #endif return 0; } } /* round off to nearest or even */ /* 0 1 2 3 4 5 6 7 8 9 10 */ /* S E M 1 2 3 4 5 6 7 8 */ k = ac1[M+8]; if( (k & 0x4000) != 0 ) { if( (k & 0x7fff) == 0x4000 ) { /* check all less significant bits */ for( j=M+9; j<=NQ; j++ ) { if( ac1[j] != 0 ) goto yesrnd; } /* round to even */ if( (k & 0x8000) == 0 ) goto nornd; } yesrnd: qclear( ac2 ); ac2[NQ] = 0; ac2[M+8] = 0x8000; addm( ac2, ac1 ); if( ac1[2] ) { shdn1(ac1); i += 1; } if( (i == 0) && (ac1[3] & 0x8000) ) i += 1; } nornd: SC = 1; shift( ac1 ); *e |= i; /* high order output already has sign bit set */ p = &ac1[M+1]; #ifdef MIEEE ++e; *e++ = *p++; *e++ = *p++; *e++ = *p++; *e++ = *p++; *e++ = *p++; *e++ = *p++; *e++ = *p++; #endif #ifdef IBMPC *(--e) = *p++; *(--e) = *p++; *(--e) = *p++; *(--e) = *p++; *(--e) = *p++; *(--e) = *p++; *(--e) = *p; #endif return 0; } /* qtoasc.c */ /* Convert q type number to ASCII string */ /* Get values for powers of ten. */ #include "qtens.h" int qtoasc( q, string, ndigs ) QELT q[]; char *string; int ndigs; { long digit; QELT x[NTT], xt[NTT]; QELT *p, *r, *ten, *tenth; QELT sign; int i, k, expon; char *s, *ss; qmov( q, x ); sign = x[SIGNWORD]; x[SIGNWORD] = 0; expon = 0; ten = &qtens[NTEN][0]; i = qcmp( qone, x ); if( i == 0 ) goto isone; if( x[1] == 0 ) { qclear( x ); goto isone; } if( i < 0 ) { k = MAXNTEN; p = &qtens[0][0]; qmov( qone, ac4 ); qmov( x, xt ); while( qcmp( ten, x ) <= 0 ) { if( qcmp( p, xt ) <= 0 ) { qdiv( p, xt, xt ); qmul( p, ac4, ac4 ); expon += k; } k >>= 1; if( k == 0 ) break; p += NTT; } qdiv( ac4, x, x ); } else { k = MINNTEN; p = &qmtens[0][0]; r = &qtens[0][0]; tenth = &qmtens[NTEN][0]; while( qcmp( tenth, x ) > 0 ) { if( qcmp( p, x ) >= 0 ) { qmul( r, x, x ); expon += k; } k /= 2; /* Prevent infinite loop due to arithmetic error: */ if( k == 0 ) break; p += NTT; r += NTT; } qmuli( ten, x, x ); expon -= 1; } isone: qifrac( x, &digit, x ); /* The following check handles numbers very close to 10**(2**n) * when there is a mistake due to arithmetic error. */ if( digit >= 10 ) { qdiv( ten, x, x ); expon += 1; digit = 1; } s = string; if( sign != 0 ) *s++ = '-'; else *s++ = ' '; *s++ = (char )digit + 060; *s++ = '.'; if( ndigs < 0 ) ndigs = 0; if( ndigs > NDEC ) ndigs = NDEC; for( k=0; k 4 ) { /* Check for critical rounding case */ if( digit == 5 ) { if( qcmp( x, qzero ) != 0 ) goto roun; /* round to nearest */ if( (*(s-1) & 1) == 0 ) goto doexp; /* round to even */ } roun: --s; k = *s & 0x7f; /* Carry out to most significant digit? */ if( k == '.' ) { --s; k = *s & 0x7f; k += 1; *s = k; /* Most significant digit rounds to 10? */ if( k > '9' ) { *s = '1'; expon += 1; } goto doexp; } /* Round up and carry out from less significant digits. */ k += 1; *s = k; if( k > '9' ) { *s = '0'; goto roun; } } doexp: sprintf( ss, "E%d", expon ); return 0; } /* QELT a[NQ], b[NQ]; * qcmp( a, b ) * * returns +1 if a > b * 0 if a == b * -1 if a < b */ int qcmp( p, q ) register QELT *p, *q; { QELT r[NQ]; register int i; int msign; if( ( *(p+1) <= (QELT)NBITS) && ( *(q+1) <= (QELT)NBITS ) ) { qsub( q, p, r ); if( r[E] == 0 ) return( 0 ); if( r[SIGNWORD] == 0 ) return( 1 ); return( -1 ); } if( *p != *q ) { /* the signs are different */ if( p[SIGNWORD] == 0 ) return( 1 ); else return( -1 ); } /* both are the same sign */ if( *p == 0 ) msign = 1; else msign = -1; i = NQ; do { if( *p++ != *q++ ) { goto diff; } } while( --i > 0 ); return(0); /* equality */ diff: if( (unsigned int) *(--p) > (unsigned int) *(--q) ) return( msign ); /* p is bigger */ else return( -msign ); /* p is littler */ } /* ; 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 ); */ int asctoq( s, y ) char *s; QELT *y; { QELT yy[NQ+1], qt[NQ]; int esign, nsign, decflg, sgnflg, nexp, exp, prec; QELT *p; nsign = 0; esign = 1; decflg = 0; sgnflg = 0; nexp = 0; exp = 0; prec = 0; qclear( yy ); yy[NQ] = 0; nxtcom: if( (*s >= '0') && (*s <= '9') ) { if( (prec == 0) && (decflg == 0) && (*s == '0') ) goto donchr; if( prec < NDEC ) { if( decflg ) nexp += 1; /* count digits after decimal point */ shup1( yy ); /* multiply current number by 10 */ qmovz( yy, ac2 ); shup1( ac2 ); shup1( ac2 ); addm( ac2, yy ); qclear( ac2 ); ac2[OMG+1] = *s - '0'; addm( ac2, yy ); } prec += 1; goto donchr; } switch( *s ) { case ' ': break; case 'E': case 'e': goto expnt; case '.': /* decimal point */ if( decflg ) goto error; ++decflg; break; case '-': nsign = -1; case '+': if( sgnflg ) goto error; ++sgnflg; break; case '\0': /* For Microware OS-9 operating system: */ #ifndef OSK case '\n': #endif case '\r': goto daldone; default: error: printf( "asctoq conversion error\n" ); qclear(y); return 0; } donchr: ++s; goto nxtcom; /* EXPONENT INTERPRETATION */ expnt: /* 0.0eXXX is zero, regardless of XXX. Check for the 0.0. */ for( exp = 0; exp < NQ + 1; exp++ ) { if( yy[exp] != 0 ) goto read_expnt; } goto zero; read_expnt: exp = 0; ++s; /* check for + or - */ if( *s == '-' ) { esign = -1; ++s; } if( *s == '+' ) ++s; while( (*s >= '0') && (*s <= '9') ) { /* Check for oversize decimal exponent. */ if( exp >= 3276 || exp < 0 ) { if( esign < 0 ) goto zero; else goto infinite; } exp *= 10; exp += *s++ - '0'; } if( esign < 0 ) exp = -exp; daldone: nexp = exp - nexp; if( normlz(yy) ) { qclear(y); return 0; } yy[E] = EXPONE - 1 + NBITS - SC; qmov( yy, y ); y[SIGNWORD] = nsign; /* Escape from excessively large exponent. */ if( nexp >= 2 * MAXNTEN ) { infinite: y[SIGNWORD] = nsign; qinfin(y); mtherr( "asctoq", OVERFLOW ); goto aexit; } if( nexp <= -2 * MAXNTEN ) { zero: qclear(y); mtherr( "asctoq", UNDERFLOW ); goto aexit; } /* multiply or divide by 10**NEXP */ if( nexp == 0 ) goto aexit; esign = 0; if( nexp < 0 ) { esign = -1; nexp = -nexp; } p = &qtens[NTEN][0]; exp = 1; qmov( qone, qt ); do { if( exp & nexp ) qmul( p, qt, qt ); exp <<= 1; p -= NQ; } while( exp <= MAXNTEN ); if( esign < 0 ) qdiv( qt, y, y ); else qmul( qt, y, y ); aexit: return 0; }