Category : C Source Code
Archive   : CEPHES22.ZIP
Filename : QFLT.C

/* 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 {
if( *p++ != 0 )
goto nota;
}
qmov( b, ac3 );
mulin( a, ac3 );
lt = (long)a[E] + (long )ac3[E] - 32768L;
goto mulcon;
}

nota:
/* detect multiplication by small integer b */
if( b[4] == 0 )
{
p = &b[5];
for( ctr=5; ctr {
if( *p++ != 0 )
goto notb;
}
qmov( a, ac3 );
mulin( b, ac3 );
lt = (long)b[E] + (long )ac3[E] - 32768L;
goto mulcon;
}

notb:

qmov( a, ac3 );
mulm( b, ac3 );
lt = (long)b[E] + (long )ac3[E] - 32768L;

mulcon:
/* calculate sign of product */
if( b[0] == a[0] )
ac3[0] = 0;
else
ac3[0] = -1;

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( "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 {
if( *a++ != *b++ )
goto difrnt;
}
return(0);

difrnt:
if( *(--a) > *(--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 *x++ = 0;
}

/*
; Fill entire number, including exponent and mantissa, with
; largest possible number.
*/

qinfin(x)
register short *x;
{
register int i;

++x; /* skip over the sign */
*x++ = 0x7fff;
*x++ = 0;
for( i=0; i *x++ = 0xffff;
}

/* normalization program */
qnrmlz(x)
short *x; /* x is the memory address of a short */
{

qmovz( x, ac1 );
normlz( ac1 ); /* shift normalize the mantissa */
ac1[E] -= SC; /* subtract the shift count from the exponent */
qmov( ac1, x );
}




/*
; Convert IEEE double precision to Q type
; double d;
; short q[N+2];
; etoq( &d, q );
*/

etoq( e, y )
short *e, *y;
{
#ifdef DEC
dtoq(e,y);
#else
register short r;
register short *p;
short yy[NQ+1];
int denorm;


denorm = 0; /* flag if denormalized number */
qclear(yy);
yy[NQ] = 0;

#ifdef MIEEE

#endif


#ifdef IBMPC
e += 3;
#endif

/*
r = *e & 0x7fff;
if( r == 0 )
return;
*/

r = *e;
yy[0] = 0;
if( r < 0 )
yy[0] = -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 += 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 {
qmuli( ten, x, x );
qifrac( x, &digit, x );
*s++ = (char )digit + 060;
}

*s = '\0'; /* mark end of string */
ss = s;

/* round off the ASCII string */

qmuli( ten, x, x );
qifrac( x, &digit, x );
if( digit > 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;
}