/* 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 * 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 * 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 * 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 * * Data structure of the number (a "word" is 16 bits) * * sign word (0 for positive, -1 for negative) * exponent (0x4001 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 "mconf.h" #include "qhead.h" /* 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 exponent */ #define E 1 /* array index of first mantissa word */ #define M 2 /* Define 1 for sticky bit rounding */ #define STICKY 0 /* accumulators */ static short ac1[NQ+1] = {0}; static short ac2[NQ+1] = {0}; static short ac3[NQ+1] = {0}; static short ac4[NQ+1] = {0}; /* shift count register */ static int SC = 0; extern short qzero[], qone[]; /* ; Absolute value ; ; short q[NQ]; ; qabs( q ); */ qabs(x) short *x; /* x is the memory address of a short */ { *x = 0; /* sign is first word of array */ } /* ; Negate ; ; short q[NQ]; ; qneg( q ); */ qneg(x) short *x; { *x = ~(*x); /* complement the sign */ } /* ; convert long (32-bit) integer to q type ; ; long l; ; short q[NQ]; ; ltoq( &l, q ); ; note &l is the memory address of l */ ltoq( lp, y ) long *lp; /* lp is the memory address of a long integer */ short *y; /* y is the address of a short */ { register short *p, *q; /* use processor registers if possible */ long ll; qclear( ac1 ); ll = *lp; /* local copy of lp */ if( ll < 0 ) { ll = -ll; /* make it positive */ ac1[0] = -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] = 040020; /* 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 */ } /* ; convert DEC double precision to q type ; double d; ; short q[NQ]; ; dtoq( &d, q ); */ dtoq( d, y ) short *d; short *y; { register short r, *p; qclear(y); /* start with a zero */ p = y; /* point to our number */ r = *d; /* get DEC exponent word */ if( *d < 0 ) *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; r >>= 7; /* shift exponent word down 7 bits */ r -= 0200; /* subtract DEC exponent offset */ r += 040000; /* 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 */ } /* ; convert q type to DEC double precision ; double d; ; short q[NQ]; ; qtod( q, &d ); */ qtod( x, d ) short *x, *d; { register short r; int i, j; *d = 0; if( *x < 0 ) *d = 0100000; qmovz( x, ac1 ); r = ac1[E]; if( r < 037600 ) 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 -= 040000; r += 0200; if( r < 0 ) { zout: *d++ = 0; *d++ = 0; *d++ = 0; *d++ = 0; return; } if( r >= 0377 ) { *d++ = 077777; *d++ = -1; *d++ = -1; *d++ = -1; return; } 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]; } /* ; Find integer and fractional parts ; long i; ; short x[NQ], frac[NQ]; ; qifrac( x, &i, frac ); */ qifrac( x, i, frac ) short *x; long *i; short *frac; { register short *p; qmovz( x, ac1 ); SC = ac1[E] - 040000; if( SC <= 0 ) { /* if exponent <= 0, integer = 0 and argument is fraction */ *i = 0L; qmov( ac1, frac ); return; } 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[0] ) *i = -(*i); ac1[0] = 0; ac1[E] = 040000; ac1[M] = 0; if( normlz(ac1) ) qclear( ac1 ); else ac1[E] -= SC; qmov( ac1, frac ); } /* ; subtract ; ; short a[NQ], b[NQ], c[NQ]; ; qsub( a, b, c ); c = b - a */ static short subflg = 0; qsub( a, b, c ) short *a, *b, *c; { subflg = 1; qadd1( a, b, c ); } /* ; add ; ; short a[NQ], b[NQ], c[NQ]; ; qadd( a, b, c ); c = b + a */ qadd( a, b, c ) short *a, *b, *c; { subflg = 0; qadd1( a, b, c ); } qadd1( a, b, c ) short *a, *b, *c; { long lt; int i; #if STICKY int lost; #endif qmovz( a, ac1 ); qmovz( b, ac2 ); if( subflg ) ac1[0] = ~ac1[0]; /* compare exponents */ SC = ac1[E] - ac2[E]; if( SC > 0 ) { /* put the larger number in ac2 */ qmovz( ac2, ac3 ); qmov( ac1, ac2 ); qmov( ac3, ac1 ); SC = -SC; } #if STICKY lost = 0; #endif if( SC != 0 ) { if( SC < -NBITS-1 ) goto done; /* answer same as larger addend */ #if STICKY lost = shift( ac1, SC ); /* shift the smaller number down */ #else shift( ac1, SC ); /* 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[0] != ac2[0] ) goto underf; /* if exponents zero, result is zero */ if( ac1[E] == 0 ) goto underf; /* if same sign, result is double */ 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[0] == ac2[0] ) { addm( ac1, ac2 ); subflg = 0; } else { subm( ac1, ac2 ); subflg = 1; } if( normlz(ac2) ) goto underf; lt = (long )ac2[E] - SC; if( lt > 32767 ) goto overf; if( lt < 0 ) { /* mtherr( "qadd", UNDERFLOW );*/ goto underf; } ac2[E] = lt; /* round off */ i = ac2[NQ] & 0xffff; if( i & 0x8000 ) { #if STICKY if( i == 0x8000 ) { 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 > 32767 ) goto overf; ac2[E] = lt; } } done: qmov( ac2, c ); return; underf: qclear(c); return; overf: mtherr( "qadd", OVERFLOW ); qinfin(c); return; } /* ; divide ; ; short 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}; */ qdiv( a, b, c ) short *a, *b, *c; { int ctr, i; long lt; double da; if( b[E] == 0 ) { qclear(c); /* numerator is zero */ return; } if( a[E] == 0 ) { /* divide by zero */ qinfin(c); mtherr( "qdiv", SING ); return; } qmovz( b, ac3 ); divm( a, ac3 ); /* calculate exponent */ lt = (long )ac3[E] - (long )a[E]; ac3[E] = lt; ac3[NQ] = 0; normlz(ac3); /* Newton iteration version */ lt = lt - SC + 16386L; /* Taylor series version */ /*lt = lt - SC + 16385L;*/ ac3[E] = lt; if( a[0] == b[0] ) ac3[0] = 0; else ac3[0] = -1; qmov( ac3, c ); } /* ; multiply ; ; short a[NQ], b[NQ], c[NQ]; ; qmul( a, b, c ); c = b * a */ qmul( a, b, c ) short *a, *b, *c; { register short *p; register int ctr; long lt; if( (a[E] == 0) || (b[E] == 0) ) { qclear(c); return; } /* detect multiplication by small integer a */ if( a[4] == 0 ) { p = &a[5]; for( ctr=5; ctr 32767 ) goto overf; ac3[E] = lt; ac3[NQ] = 0; if( normlz(ac3) ) goto underf; lt = lt - SC + 16384L; if( lt > 32767 ) goto overf; if( lt < 0 ) goto underf; ac3[E] = lt; qmov( ac3, c ); return; underf: qclear(c); return; overf: qinfin(c); mtherr( "qmul", OVERFLOW ); return; } /* Multiply, a has at most 16 significant bits */ qmuli( a, b, c ) short *a, *b, *c; { int ctr; long lt; if( (a[E] == 0) || (b[E] == 0) ) { qclear(c); return; } qmov( b, ac3 ); mulin( a, ac3 ); /* calculate sign of product */ if( b[0] == a[0] ) ac3[0] = 0; else ac3[0] = -1; /* calculate exponent */ lt = (long)ac3[E] + (long )a[E] - 32768L; if( lt > 32767 ) goto overf; ac3[E] = lt; ac3[NQ] = 0; if( normlz(ac3) ) goto underf; lt = lt - SC + 16384L; if( lt > 32767 ) goto overf; if( lt < 0 ) goto underf; ac3[E] = lt; qmov( ac3, c ); return; underf: qclear(c); return; overf: qinfin(c); mtherr( "qmuli", OVERFLOW ); return; } /* ; Compare mantissas ; ; short a[NQ], b[NQ]; ; cmpm( a, b ); ; ; for the mantissas: ; returns +1 if a > b ; 0 if a == b ; -1 if a < b */ cmpm( a, b ) register unsigned short *a, *b; { int i; a += 2; /* skip up to mantissa area */ b += 2; for( i=0; i *(--b) ) return(1); else return(-1); } /* ; shift mantissa ; ; Shifts mantissa area up or down by the number of bits ; given by the variable SC. */ shift( x ) short *x; { short *p; #if STICKY int lost; #endif if( SC == 0 ) { #if STICKY return(0); #else return; #endif } #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 ); #endif } /* ; normalize ; ; shift normalizes the mantissa area pointed to by R1 ; shift count (up = positive) returned in SC */ normlz(x) short x[]; { register short *p; SC = 0; p = &x[M]; if( *p != 0 ) goto normdn; ++p; if( *p < 0 ) return(0); /* already normalized */ while( *p == 0 ) { shup16(x); SC += 16; /* With guard word, there are NBITS+16 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 > 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 ; ; short x[]; ; qclear( x ); */ qclear( x ) register short *x; { register int i; for( i=0; i>= 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 += 040001 - 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 } /* ; Q type to IEEE double precision ; double d; ; short q[N+2]; ; qtoe( q, &d ); */ qtoe( x, e ) short *x, *e; { #ifdef DEC qtod(x,e); #else short i, j, k; register short *p; #ifdef MIEEE #endif #ifdef IBMPC e += 3; #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++ -(040001 - 01777); /* adjust exponent for offsets */ /* Handle denormalized small numbers. * The tiniest number represented appears to be 2**-1074. */ 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; } } /* 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] & 0x8000) ) 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; } 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 } /* qtoasc.c */ /* Convert q type number to ASCII string */ /*include "qhead.h"*/ #define NTEN 12 #if NQ < 13 #define NTT 12 static short qtens[NTEN+1][NTT] = { /* 4096 */ {0x0000,0x7527,0x0000,0xc460,0x5202,0x8a20,0x979a,0xc94c, 0x153f,0x804a,0x4a92,0x6576,}, /* 2048 */ {0x0000,0x5a94,0x0000,0x9e8b,0x3b5d,0xc53d,0x5de4,0xa74d, 0x28ce,0x329a,0xce52,0x6a32,}, /* 1024 */ {0x0000,0x4d4a,0x0000,0xc976,0x7586,0x8175,0x0c17,0x650d, 0x3d28,0xf18b,0x50ce,0x526c,}, /* 512 */ {0x0000,0x46a5,0x0000,0xe319,0xa0ae,0xa60e,0x91c6,0xcc65, 0x5c54,0xbc50,0x58f8,0x9c66,}, /* 256 */ {0x0000,0x4353,0x0000,0xaa7e,0xebfb,0x9df9,0xde8d,0xddbb, 0x901b,0x98fe,0xeab7,0x851e,}, /* 128 */ {0x0000,0x41aa,0x0000,0x93ba,0x47c9,0x80e9,0x8cdf,0xc66f, 0x336c,0x36b1,0x0137,0x0235,}, /* 64 */ {0x0000,0x40d5,0x0000,0xc278,0x1f49,0xffcf,0xa6d5,0x3cbf, 0x6b71,0xc76b,0x25fb,0x50f8,}, /* 32 */ {0000000,0040153,0000000,0116705,0126650,0025560, 0132635,0170040,0000000,0000000,0000000,0000000,}, /* 16 */ {0000000,0040066,0000000,0107033,0144677,0002000, 0000000,0000000,0000000,0000000,0000000,0000000,}, /* 8 */ {0000000,0040033,0000000,0137274,0020000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000,}, /* 4 */ {0000000,0040016,0000000,0116100,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000,}, /* 2 */ {0000000,0040007,0000000,0144000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000,}, /* 1 */ {0000000,0040004,0000000,0120000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000,}, }; static short qmtens[NTEN+1][NTT] = { /* -4096 */ {0x0000,0x0ada,0x0000,0xa6dd,0x04c8,0xd2ce,0x9fde,0x2de3, 0x8123,0xa1c3,0xcffc,0x2030,}, /* -2048 */ {0x0000,0x256d,0x0000,0xceae,0x534f,0x3436,0x2de4,0x4925, 0x12d4,0xf2ea,0xd2cb,0x8264,}, /* -1024 */ {0x0000,0x32b7,0x0000,0xa2a6,0x82a5,0xda57,0xc0bd,0x87a6, 0x0158,0x6bd3,0xf698,0xf53f,}, /* -512 */ {0x0000,0x395c,0x0000,0x9049,0xee32,0xdb23,0xd21c,0x7132, 0xd332,0xe3f2,0x04d4,0xe731,}, /* -256 */ {0x0000,0x3cae,0x0000,0xc031,0x4325,0x637a,0x1939,0xfa91, 0x1155,0xfefb,0x5308,0xa23e,}, /* -128 */ {0x0000,0x3e57,0x0000,0xddd0,0x467c,0x64bc,0xe4a0,0xac7c, 0xb3f6,0xd05d,0xdbde,0xe26d,}, /* -64 */ {0x0000,0x3f2c,0x0000,0xa87f,0xea27,0xa539,0xe9a5,0x3f23, 0x98d7,0x47b3,0x6224,0x2a20,}, /* -32 */ {0x0000,0x3f96,0x0000,0xcfb1,0x1ead,0x4539,0x94ba,0x67de, 0x18ed,0xa581,0x4af2,0x0b5b,}, /* -16 */ {0x0000,0x3fcb,0x0000,0xe695,0x94be,0xc44d,0xe15b,0x4c2e, 0xbe68,0x7989,0xa9b3,0xbf71,}, /* -8 */ {0x0000,0x3fe6,0x0000,0xabcc,0x7711,0x8461,0xcefc,0xfdc2, 0x0d2b,0x36ba,0x7c3d,0x3d4d,}, /* -4 */ {0x0000,0x3ff3,0x0000,0xd1b7,0x1758,0xe219,0x652b,0xd3c3, 0x6113,0x404e,0xa4a8,0xc155,}, /* -2 */ {0x0000,0x3ffa,0x0000,0xa3d7,0x0a3d,0x70a3,0xd70a,0x3d70, 0xa3d7,0x0a3d,0x70a3,0xd70a,}, /* -1 */ {0x0000,0x3ffd,0x0000,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc, 0xcccc,0xcccc,0xcccc,0xcccd,}, }; #else #define NTT 24 static short qtens[NTEN+1][NTT] = { /* 4096 */ {0x0000,0x7527,0x0000,0xc460,0x5202,0x8a20,0x979a,0xc94c, 0x153f,0x804a,0x4a92,0x6576,0x1fb2,0x444e,0x2267,0xdd5c, 0xf7c9,0x45f2,0x2a3f,0xfff0,0xd1d2,0xe184,0x69b1,0xea39,}, /* 2048 */ {0x0000,0x5a94,0x0000,0x9e8b,0x3b5d,0xc53d,0x5de4,0xa74d, 0x28ce,0x329a,0xce52,0x6a31,0x97bb,0xebe3,0x034f,0x7715, 0x4ce2,0xbcba,0x1964,0x8b21,0xc11e,0xb962,0xb1b6,0x1b94,}, /* 1024 */ {0x0000,0x4d4a,0x0000,0xc976,0x7586,0x8175,0x0c17,0x650d, 0x3d28,0xf18b,0x50ce,0x526b,0x9882,0x7524,0x9b0f,0xd6f4, 0xb6d2,0x7bd1,0xc61c,0x2530,0x69a5,0xc329,0xf9af,0x4a9c,}, /* 512 */ {0x0000,0x46a5,0x0000,0xe319,0xa0ae,0xa60e,0x91c6,0xcc65, 0x5c54,0xbc50,0x58f8,0x9c65,0x8398,0x1d13,0x4cba,0x422d, 0x38ea,0x3584,0xcde4,0x0b9a,0x1d7d,0x634d,0xf2d8,0x74a2,}, /* 256 */ {0x0000,0x4353,0x0000,0xaa7e,0xebfb,0x9df9,0xde8d,0xddbb, 0x901b,0x98fe,0xeab7,0x851e,0x4cbf,0x3de2,0xf98a,0xae78, 0x0c7f,0xea81,0xc788,0x5a6b,0x43a2,0xa6c4,0x95b8,0xccd6,}, /* 128 */ {0x0000,0x41aa,0x0000,0x93ba,0x47c9,0x80e9,0x8cdf,0xc66f, 0x336c,0x36b1,0x0137,0x0234,0xf3fd,0x7b08,0xdd39,0x0bc3, 0xc54e,0x3f40,0xf7e6,0x424b,0xa54f,0x8040,0x0000,0x0000,}, /* 64 */ {0x0000,0x40d5,0x0000,0xc278,0x1f49,0xffcf,0xa6d5,0x3cbf, 0x6b71,0xc76b,0x25fb,0x50f8,0x0800,0x0000,0x0000,0x0000, 0x0000,0x0000,0x0000,0x0000,0x0000,0x0000,0x0000,0x0000,}, /* 32 */ {0000000,0040153,0000000,0116705,0126650,0025560, 0132635,0170040,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000,}, /* 16 */ {0000000,0040066,0000000,0107033,0144677,0002000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000,}, /* 8 */ {0000000,0040033,0000000,0137274,0020000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000,}, /* 4 */ {0000000,0040016,0000000,0116100,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000,}, /* 2 */ {0000000,0040007,0000000,0144000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000,}, /* 1 */ {0000000,0040004,0000000,0120000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000, 0000000,0000000,0000000,0000000,0000000,0000000,}, }; static short qmtens[NTEN+1][NTT] = { /* -4096 */ {0x0000,0x0ada,0x0000,0xa6dd,0x04c8,0xd2ce,0x9fde,0x2de3, 0x8123,0xa1c3,0xcffc,0x2030,0x5d02,0x44e0,0x91ba,0x5e2d, 0x7403,0x972f,0x6f2b,0x04e6,0x63ac,0x16d5,0x1215,0x9834,}, /* -2048 */ {0x0000,0x256d,0x0000,0xceae,0x534f,0x3436,0x2de4,0x4925, 0x12d4,0xf2ea,0xd2cb,0x8263,0xca5c,0xbc77,0x4bd9,0x71aa, 0xd590,0x46c7,0x4249,0xdddd,0xc0df,0xbb0d,0x6f2c,0x0584,}, /* -1024 */ {0x0000,0x32b7,0x0000,0xa2a6,0x82a5,0xda57,0xc0bd,0x87a6, 0x0158,0x6bd3,0xf698,0xf53e,0x94d1,0xb235,0x7c32,0xc0ea, 0xff37,0x55a2,0xddcd,0x5191,0xa70f,0x9e33,0x9aa9,0x6270,}, /* -512 */ {0x0000,0x395c,0x0000,0x9049,0xee32,0xdb23,0xd21c,0x7132, 0xd332,0xe3f2,0x04d4,0xe731,0x7d62,0x209b,0x6a93,0xd4c9, 0x4a9d,0xa069,0x3e0c,0xfefd,0xe89b,0x01e5,0x1068,0xe0e4,}, /* -256 */ {0x0000,0x3cae,0x0000,0xc031,0x4325,0x637a,0x1939,0xfa91, 0x1155,0xfefb,0x5308,0xa23e,0x2ed2,0x7766,0xe8cc,0x9b03, 0x5377,0x08b1,0x648f,0x7a73,0x0e52,0xe6e7,0x4360,0x55f4,}, /* -128 */ {0x0000,0x3e57,0x0000,0xddd0,0x467c,0x64bc,0xe4a0,0xac7c, 0xb3f6,0xd05d,0xdbde,0xe26c,0xa606,0x3461,0xfffa,0x4ed7, 0x75fc,0x49f2,0x7952,0xf2f3,0xa45f,0x9009,0xf3c9,0xee49,}, /* -64 */ {0x0000,0x3f2c,0x0000,0xa87f,0xea27,0xa539,0xe9a5,0x3f23, 0x98d7,0x47b3,0x6224,0x2a1f,0xee40,0xd90a,0xab31,0x0e12, 0x8b5d,0x938c,0xfb3f,0x6f20,0x3147,0x025d,0x1129,0xe492,}, /* -32 */ {0x0000,0x3f96,0x0000,0xcfb1,0x1ead,0x4539,0x94ba,0x67de, 0x18ed,0xa581,0x4af2,0x0b5b,0x1aa0,0x28cc,0xd99e,0x59e3, 0x38e3,0x87ad,0x8e28,0x5676,0x0892,0x1621,0x96af,0x3c7f,}, /* -16 */ {0x0000,0x3fcb,0x0000,0xe695,0x94be,0xc44d,0xe15b,0x4c2e, 0xbe68,0x7989,0xa9b3,0xbf71,0x6c1a,0xdd27,0xf085,0x23cc, 0xd348,0x4db6,0x70aa,0xd6e0,0x12cf,0x7fa7,0xcf42,0x8583,}, /* -8 */ {0x0000,0x3fe6,0x0000,0xabcc,0x7711,0x8461,0xcefc,0xfdc2, 0x0d2b,0x36ba,0x7c3d,0x3d4d,0x3d75,0x8161,0x697c,0x7068, 0xf3b4,0x6d2f,0x8350,0x5705,0x755f,0xd37a,0x3b04,0x8dd9,}, /* -4 */ {0x0000,0x3ff3,0x0000,0xd1b7,0x1758,0xe219,0x652b,0xd3c3, 0x6113,0x404e,0xa4a8,0xc154,0xc985,0xf06f,0x6944,0x6738, 0x1d7d,0xbf48,0x7fcb,0x923a,0x29c7,0x79a6,0xb50b,0x0f28,}, /* -2 */ {0x0000,0x3ffa,0x0000,0xa3d7,0x0a3d,0x70a3,0xd70a,0x3d70, 0xa3d7,0x0a3d,0x70a3,0xd70a,0x3d70,0xa3d7,0x0a3d,0x70a3, 0xd70a,0x3d70,0xa3d7,0x0a3d,0x70a3,0xd70a,0x3d70,0xa3d7,}, /* -1 */ {0x0000,0x3ffd,0x0000,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc, 0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc, 0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccd,}, }; #endif qtoasc( q, string, ndigs ) short q[]; char *string; int ndigs; { long digit; short x[NTT], xt[NTT]; short *p, *r, *ten, *tenth; short sign; int i, k, expon; char *s, *ss; qmov( q, x ); sign = x[0]; x[0] = 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 = 4096; 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 = -4096; 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\0", expon ); } /* short a[NQ], b[NQ]; * qcmp( a, b ) * * returns +1 if a > b * 0 if a == b * -1 if a < b */ qcmp( p, q ) register unsigned short *p, *q; { unsigned short r[NQ]; register int i; int msign; if( ( *(p+1) <= NBITS) && ( *(q+1) <= NBITS ) ) { qsub( q, p, r ); if( r[E] == 0 ) return( 0 ); if( r[0] == 0 ) return( 1 ); return( -1 ); } if( *p != *q ) { /* the signs are different */ if( *p == 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( *(--p) > *(--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 ); */ /*if NQ < 13*/ /*static short qten[NTT] = {0,040004,0,0120000,0,0,0,0,0,0,0,0};*/ static short qten[24] = { 0,040004,0,0120000,0,0,0,0,0,0,0,0, 0,0,0,0,0,0,0,0,0,0,0,0 }; asctoq( s, y ) char *s; short *y; { short yy[NQ+1], qt[NQ]; int esign, nsign, decflg, sgnflg, nexp, exp, prec; short *p; nsign = 0; 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; } donchr: ++s; goto nxtcom; /* EXPONENT INTERPRETATION */ expnt: esign = 1; exp = 0; ++s; /* check for + or - */ if( *s == '-' ) { esign = -1; ++s; } if( *s == '+' ) ++s; while( (*s >= '0') && (*s <= '9') ) { exp *= 10; exp += *s++ - '0'; } if( esign < 0 ) exp = -exp; daldone: nexp = exp - nexp; if( normlz(yy) ) { qclear(y); return; } yy[E] = 040000 + NBITS - SC; qmov( yy, y ); y[0] = nsign; /* 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 < 8192 ); if( esign < 0 ) qdiv( qt, y, y ); else qmul( qt, y, y ); aexit: return; }