home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / qfloat / qflt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  29.6 KB  |  1,778 lines
  1. /*                            qflt.c
  2.  *            QFLOAT
  3.  *
  4.  *    Extended precision floating point routines
  5.  *
  6.  *    asctoq( string, q )    ascii string to q type
  7.  *    dtoq( &d, q )        DEC double precision to q type
  8.  *    etoq( &d, q )        IEEE double precision to q type
  9.  *    ltoq( &l, q )        long integer to q type
  10.  *    qabs(q)            absolute value
  11.  *    qadd( a, b, c )        c = b + a
  12.  *    qclear(q)        q = 0
  13.  *    qcmp( a, b )        compare a to b
  14.  *    qdiv( a, b, c )        c = b / a
  15.  *    qifrac( x, &l, frac )   x to integer part l and q type fraction
  16.  *    qinfin( x )        set x to infinity, leaving its sign alone
  17.  *    qmov( a, b )        b = a
  18.  *    qmul( a, b, c )        c = b * a
  19.  *    qmuli( a, b, c )    c = b * a, a has only 16 significant bits
  20.  *    qneg(q)            q = -q
  21.  *    qnrmlz(q)        adjust exponent and mantissa
  22.  *    qsub( a, b, c )        c = b - a
  23.  *    qtoasc( a, s, n )    q to ASCII string, n digits after decimal
  24.  *    qtod( q, &d )        convert q type to DEC double precision
  25.  *    qtoe( q, &d )        convert q type to IEEE double precision
  26.  *
  27.  * Data structure of the number (a "word" is 16 bits)
  28.  *
  29.  *    sign word        (0 for positive, -1 for negative)
  30.  *    exponent        (0x4001 for 1.0)
  31.  *    high guard word        (always zero after normalization)
  32.  *    N-1 mantissa words    (most significant word first,
  33.  *                 most significant bit is set)
  34.  *
  35.  * Numbers are stored in C language as arrays.  All routines
  36.  * use pointers to the arrays as arguments.
  37.  *
  38.  * The result is always normalized after each arithmetic operation.
  39.  * All arithmetic results are chopped. No rounding is performed except
  40.  * on conversion to double precision.
  41.  */
  42.  
  43. /*
  44.  * Revision history:
  45.  *
  46.  * SLM,  5 Jan 84    PDP-11 assembly language version
  47.  * SLM,  2 Mar 86    fixed bug in asctoq()
  48.  * SLM,  6 Dec 86    C language version
  49.  *
  50.  */
  51.  
  52. #include "mconf.h"
  53. #include "qhead.h"
  54.  
  55.     /* number of 16 bit words in mantissa area */
  56. /*define N 10*/
  57.     /* max number of decimal digits in conversion */
  58. /*define NDEC 46*/
  59.     /* number of bits of precision */
  60. /*define NBITS (OMG-1)*16*/
  61.     /* array index of the exponent */
  62. #define E 1
  63.     /* array index of first mantissa word */
  64. #define M 2
  65.  
  66. /* Define 1 for sticky bit rounding */
  67. #define STICKY 0
  68.  
  69.     /* accumulators */
  70. static short ac1[NQ+1] = {0};
  71. static short ac2[NQ+1] = {0};
  72. static short ac3[NQ+1] = {0};
  73. static short ac4[NQ+1] = {0};
  74.     /* shift count register */
  75. static int SC = 0;
  76.  
  77. extern short qzero[], qone[];
  78.  
  79.  
  80. /*
  81. ;    Absolute value
  82. ;
  83. ;    short q[NQ];
  84. ;    qabs( q );
  85. */
  86.  
  87. qabs(x)
  88. short *x;    /* x is the memory address of a short */
  89. {
  90.  
  91. *x = 0;        /* sign is first word of array */
  92. }
  93.  
  94.  
  95.  
  96.  
  97. /*
  98. ;    Negate
  99. ;
  100. ;    short q[NQ];
  101. ;    qneg( q );
  102. */
  103.  
  104. qneg(x)
  105. short *x;
  106. {
  107.  
  108. *x = ~(*x);    /* complement the sign */
  109. }
  110.  
  111. /*
  112. ; convert long (32-bit) integer to q type
  113. ;
  114. ;    long l;
  115. ;    short q[NQ];
  116. ;    ltoq( &l, q );
  117. ; note &l is the memory address of l
  118. */
  119.  
  120. ltoq( lp, y )
  121. long *lp;    /* lp is the memory address of a long integer */
  122. short *y;    /* y is the address of a short */
  123. {
  124. register short *p, *q;    /* use processor registers if possible */
  125. long ll;
  126.  
  127. qclear( ac1 );
  128. ll = *lp;        /* local copy of lp */
  129. if( ll < 0 )
  130.     {
  131.     ll = -ll;    /* make it positive */
  132.     ac1[0] = -1;    /* put correct sign in the q type number */
  133.     }
  134.  
  135. /* move the long integer to ac1 mantissa area */
  136. p = &ac1[M];
  137. *p++ = ll >> 16;
  138. *p = ll;
  139. /*
  140.  * q = (short *)≪
  141.  * #if DEC
  142.  * p = &ac1[M];
  143.  * *p++ = *q++;
  144.  * #endif
  145.  * #if IBMPC
  146.  * p = &ac1[M+1];
  147.  * *p-- = *q++;
  148.  * #endif
  149.  * *p = *q;
  150.  */
  151. ac1[E] = 040020;    /* exponent if normalize shift count were 0 */
  152. ac1[NQ] = 0;
  153. if( normlz( ac1 ) )    /* normalize the mantissa */
  154.     qclear( ac1 );    /* it was zero */
  155. else
  156.     ac1[E] -= SC;    /* else subtract shift count from exponent */
  157. qmov( ac1, y );        /* output the answer */
  158. }
  159.  
  160. /*
  161. ;    convert DEC double precision to q type
  162. ;    double d;
  163. ;    short q[NQ];
  164. ;    dtoq( &d, q );
  165. */
  166. dtoq( d, y )
  167. short *d;
  168. short *y;
  169. {
  170. register short r, *p;
  171.  
  172. qclear(y);        /* start with a zero */
  173. p = y;            /* point to our number */
  174. r = *d;            /* get DEC exponent word */
  175. if( *d < 0 )
  176.     *p = -1;    /* fill in our sign */
  177. ++p;            /* bump pointer to our exponent word */
  178. r &= 0x7fff;        /* strip the sign bit */
  179. if( r == 0 )        /* answer = 0 if high order DEC word = 0 */
  180.     return;
  181.  
  182.  
  183. r >>= 7;    /* shift exponent word down 7 bits */
  184. r -= 0200;    /* subtract DEC exponent offset */
  185. r += 040000;    /* add our q type exponent offset */
  186. *p++ = r;    /* to form our exponent */
  187.  
  188. r = *d++;    /* now do the high order mantissa */
  189. r &= 0177;    /* strip off the DEC exponent and sign bits */
  190. r |= 0200;    /* the DEC understood high order mantissa bit */
  191. *p++ = r;    /* put result in our high guard word */
  192.  
  193. *p++ = *d++;    /* fill in the rest of our mantissa */
  194. *p++ = *d++;
  195. *p = *d;
  196.  
  197. shdn8(y);    /* shift our mantissa down 8 bits */
  198. }
  199.  
  200. /*
  201. ;    convert q type to DEC double precision
  202. ;    double d;
  203. ;    short q[NQ];
  204. ;    qtod( q, &d );
  205. */
  206.  
  207. qtod( x, d )
  208. short *x, *d;
  209. {
  210. register short r;
  211. int i, j;
  212.  
  213. *d = 0;
  214. if( *x < 0 )
  215.     *d = 0100000;
  216. qmovz( x, ac1 );
  217. r = ac1[E];
  218. if( r < 037600 )
  219.     goto zout;
  220. i = ac1[M+4];
  221. if( (i & 0200) != 0 )
  222.     {
  223.     if( (i & 0377) == 0200 )
  224.         {
  225.         if( (i & 0400) != 0 )
  226.             {
  227.         /* check all less significant bits */
  228.             for( j=M+5; j<=NQ; j++ )
  229.                 {
  230.                 if( ac1[j] != 0 )
  231.                     goto yesrnd;
  232.                 }
  233.             }
  234.         goto nornd;
  235.         }
  236. yesrnd:
  237.     qclear( ac2 );
  238.     ac2[ M+4 ] = 0200;
  239.     ac2[NQ] = 0;
  240.     addm( ac2, ac1 );
  241.     normlz(ac1);
  242.     r -= SC;
  243.     }
  244.  
  245. nornd:
  246.  
  247. r -= 040000;
  248. r += 0200;
  249. if( r < 0 )
  250.     {
  251. zout:
  252.     *d++ = 0;
  253.     *d++ = 0;
  254.     *d++ = 0;
  255.     *d++ = 0;
  256.     return;
  257.     }
  258. if( r >= 0377 )
  259.     {
  260.     *d++ = 077777;
  261.     *d++ = -1;
  262.     *d++ = -1;
  263.     *d++ = -1;
  264.     return;
  265.     }
  266. r &= 0377;
  267. r <<= 7;
  268. shup8( ac1 );
  269. ac1[M] &= 0177;
  270. r |= ac1[M];
  271. *d++ |= r;
  272. *d++ = ac1[M+1];
  273. *d++ = ac1[M+2];
  274. *d++ = ac1[M+3];
  275. }
  276.  
  277. /*
  278. ;    Find integer and fractional parts
  279.  
  280. ;    long i;
  281. ;    short x[NQ], frac[NQ];
  282. ;    qifrac( x, &i, frac );
  283. */
  284.  
  285. qifrac( x, i, frac )
  286. short *x;
  287. long *i;
  288. short *frac;
  289. {
  290. register short *p;
  291.  
  292. qmovz( x, ac1 );
  293. SC = ac1[E] - 040000;
  294. if( SC <= 0 )
  295.     {
  296. /* if exponent <= 0, integer = 0 and argument is fraction */
  297.     *i = 0L;
  298.     qmov( ac1, frac );
  299.     return;
  300.     }
  301. if( SC > 31 )
  302.     {
  303. /*
  304. ;    long integer overflow: output large integer
  305. ;    and correct fraction
  306. */
  307.     *i = 0x7fffffff;
  308.     shift( ac1 );
  309.     goto lab10;
  310.     }
  311.  
  312. if( SC > 16 )
  313.     {
  314. /*
  315. ; shift more than 16 bits: shift up SC-16, output the integer,
  316. ; then complete the shift to get the fraction.
  317. */
  318.     SC -= 16;
  319.     shift( ac1 );
  320.     *i = ((unsigned long )ac1[M] << 16) | (unsigned short )ac1[M+1];
  321. /*
  322.  *     p = (short *)i;
  323.  * #ifdef DEC
  324.  *     *p++ = ac1[M];
  325.  *     *p++ = ac1[M+1];
  326.  * #else
  327.  *     *p++ = ac1[M+1];
  328.  *     *p++ = ac1[M];
  329.  * #endif
  330.  */
  331.     shup16( ac1 );
  332.     goto lab10;
  333.     }
  334.  
  335. /* shift not more than 16 bits */
  336. shift( ac1 );
  337. *i = ac1[M] & 0xffff;
  338.  
  339. lab10:
  340. if( x[0] )
  341.     *i = -(*i);
  342. ac1[0] = 0;
  343. ac1[E] = 040000;
  344. ac1[M] = 0;
  345. if( normlz(ac1) )
  346.     qclear( ac1 );
  347. else
  348.     ac1[E] -= SC;
  349.  
  350. qmov( ac1, frac );
  351. }
  352.  
  353.  
  354. /*
  355. ;    subtract
  356. ;
  357. ;    short a[NQ], b[NQ], c[NQ];
  358. ;    qsub( a, b, c );     c = b - a
  359. */
  360.  
  361. static short subflg = 0;
  362.  
  363. qsub( a, b, c )
  364. short *a, *b, *c;
  365. {
  366.  
  367. subflg = 1;
  368. qadd1( a, b, c );
  369. }
  370.  
  371.  
  372. /*
  373. ;    add
  374. ;
  375. ;    short a[NQ], b[NQ], c[NQ];
  376. ;    qadd( a, b, c );     c = b + a
  377. */
  378.  
  379. qadd( a, b, c )
  380. short *a, *b, *c;
  381. {
  382.  
  383. subflg = 0;
  384. qadd1( a, b, c );
  385. }
  386.  
  387. qadd1( a, b, c )
  388. short *a, *b, *c;
  389. {
  390. long lt;
  391. int i;
  392. #if STICKY
  393. int lost;
  394. #endif
  395.  
  396. qmovz( a, ac1 );
  397. qmovz( b, ac2 );
  398. if( subflg )
  399.     ac1[0] = ~ac1[0];
  400.  
  401. /* compare exponents */
  402. SC = ac1[E] - ac2[E]; 
  403. if( SC > 0 )
  404.     {    /* put the larger number in ac2 */
  405.     qmovz( ac2, ac3 );
  406.     qmov( ac1, ac2 );
  407.     qmov( ac3, ac1 );
  408.     SC = -SC;
  409.     }
  410.  
  411. #if STICKY
  412. lost = 0;
  413. #endif
  414. if( SC != 0 )
  415.     {
  416.     if( SC < -NBITS-1 )
  417.         goto done;    /* answer same as larger addend */
  418. #if STICKY
  419.     lost = shift( ac1, SC ); /* shift the smaller number down */
  420. #else
  421.     shift( ac1, SC ); /* shift the smaller number down */
  422. #endif
  423.     }
  424. else
  425.     {
  426. /* exponents were the same, so must compare mantissae */
  427.     i = cmpm( ac1, ac2 );
  428.     if( i == 0 )
  429.         {    /* the numbers are identical */
  430.         /* if different signs, result is zero */
  431.         if( ac1[0] != ac2[0] )
  432.             goto underf;
  433.         /* if exponents zero, result is zero */
  434.         if( ac1[E] == 0 )
  435.             goto underf;
  436.         /* if same sign, result is double */
  437.         ac2[E] += 1;
  438.         goto done;
  439.         }
  440.     if( i > 0 )
  441.         {    /* put the larger number in ac2 */
  442.         qmovz( ac2, ac3 );
  443.         qmov( ac1, ac2 );
  444.         qmov( ac3, ac1 );
  445.         }
  446.     }
  447.  
  448. if( ac1[0] == ac2[0] )
  449.     {
  450.     addm( ac1, ac2 );
  451.     subflg = 0;
  452.     }
  453. else
  454.     {
  455.     subm( ac1, ac2 );
  456.     subflg = 1;
  457.     }
  458. if( normlz(ac2) )
  459.     goto underf;
  460.  
  461. lt = (long )ac2[E] - SC;
  462. if( lt > 32767 )
  463.     goto overf;
  464. if( lt < 0 )
  465.     {
  466. /*    mtherr( "qadd", UNDERFLOW );*/
  467.     goto underf;    
  468.     }
  469. ac2[E] = lt;
  470.  
  471. /* round off */
  472. i = ac2[NQ] & 0xffff;
  473. if( i & 0x8000 )
  474.     {
  475. #if STICKY
  476.     if( i == 0x8000 )
  477.         {
  478.         if( lost == 0 )
  479.             {
  480. /* Critical case, round to even */
  481.             if( (ac2[NQ-1] & 1) == 0 )
  482.                 goto done;
  483.             }
  484.         else
  485.             {
  486.             if( subflg != 0 )
  487.                 goto done;
  488.             }
  489.         }
  490. #else
  491.     if( subflg != 0 )
  492.         goto done;
  493. #endif
  494.     
  495.     qclear( ac1 );
  496.     ac1[NQ] = 0;
  497.     ac1[NQ-1] = 1;
  498.     addm( ac1, ac2 );
  499.     normlz(ac2);
  500.     if( SC )
  501.         {
  502.         lt = (long )ac2[E] - SC;
  503.         if( lt > 32767 )
  504.         goto overf;
  505.         ac2[E] = lt;
  506.         }
  507.     }
  508. done:
  509. qmov( ac2, c );
  510. return;
  511.  
  512. underf:
  513. qclear(c);
  514. return;
  515.  
  516. overf:
  517. mtherr( "qadd", OVERFLOW );
  518. qinfin(c);
  519. return;
  520. }
  521.  
  522. /*
  523. ;    divide
  524. ;
  525. ;    short a[NQ], b[NQ], c[NQ];
  526. ;    qdiv( a, b, c );    c = b / a
  527. */
  528. /* for Newton iteration version:
  529.  * extern short qtwo[];
  530.  * static short qt[NQ] = {0};
  531.  * static short qu[NQ] = {0};
  532.  */
  533. qdiv( a, b, c )
  534. short *a, *b, *c;
  535. {
  536. int ctr, i;
  537. long lt;
  538. double da;
  539.  
  540. if( b[E] == 0 )
  541.     {
  542.     qclear(c);    /* numerator is zero */
  543.     return;
  544.     }
  545.  
  546. if( a[E] == 0 )
  547.     {    /* divide by zero */
  548.     qinfin(c);
  549.     mtherr( "qdiv", SING );
  550.     return;
  551.     }
  552. qmovz( b, ac3 );
  553. divm( a, ac3 );
  554. /* calculate exponent */
  555. lt = (long )ac3[E] - (long )a[E];
  556. ac3[E] = lt;
  557. ac3[NQ] = 0;
  558. normlz(ac3);
  559. /* Newton iteration version */
  560. lt = lt - SC + 16386L;
  561. /* Taylor series version */
  562. /*lt = lt - SC + 16385L;*/
  563. ac3[E] = lt;
  564.  
  565. if( a[0] == b[0] )
  566.     ac3[0] = 0;
  567. else
  568.     ac3[0] = -1;
  569.  
  570. qmov( ac3, c );
  571. }
  572.  
  573. /*
  574. ;    multiply
  575. ;
  576. ;    short a[NQ], b[NQ], c[NQ];
  577. ;    qmul( a, b, c );    c = b * a
  578. */
  579. qmul( a, b, c )
  580. short *a, *b, *c;
  581. {
  582. register short *p;
  583. register int ctr;
  584. long lt;
  585.  
  586. if( (a[E] == 0) || (b[E] == 0) )
  587.     {
  588.     qclear(c);
  589.     return;
  590.     }
  591. /* detect multiplication by small integer a */
  592. if( a[4] == 0 )
  593.     {
  594.     p = &a[5];
  595.     for( ctr=5; ctr<NQ; ctr++ )
  596.         {
  597.         if( *p++ != 0 )
  598.             goto nota;
  599.         }
  600.     qmov( b, ac3 );
  601.     mulin( a, ac3 );
  602.     lt = (long)a[E] + (long )ac3[E] - 32768L;
  603.     goto mulcon;
  604.     }
  605.  
  606. nota:
  607. /* detect multiplication by small integer b */
  608. if( b[4] == 0 )
  609.     {
  610.     p = &b[5];
  611.     for( ctr=5; ctr<NQ; ctr++ )
  612.         {
  613.         if( *p++ != 0 )
  614.             goto notb;
  615.         }
  616.     qmov( a, ac3 );
  617.     mulin( b, ac3 );
  618.     lt = (long)b[E] + (long )ac3[E] - 32768L;
  619.     goto mulcon;
  620.     }
  621.  
  622. notb:
  623.  
  624. qmov( a, ac3 );
  625. mulm( b, ac3 );
  626. lt = (long)b[E] + (long )ac3[E] - 32768L;
  627.  
  628. mulcon:
  629. /* calculate sign of product */
  630. if( b[0] == a[0] )
  631.     ac3[0] = 0;
  632. else
  633.     ac3[0] = -1;
  634.  
  635. if( lt > 32767 )
  636.     goto overf;
  637. ac3[E] = lt;
  638. ac3[NQ] = 0;
  639. if( normlz(ac3) )
  640.     goto underf;
  641. lt = lt - SC + 16384L;
  642. if( lt > 32767 )
  643.     goto overf;
  644. if( lt < 0 )
  645.     goto underf;
  646. ac3[E] = lt;
  647. qmov( ac3, c );
  648. return;
  649.  
  650. underf:
  651. qclear(c);
  652. return;
  653.  
  654. overf:
  655. qinfin(c);
  656. mtherr( "qmul", OVERFLOW );
  657. return;
  658. }
  659.  
  660.  
  661.  
  662.  
  663. /* Multiply, a has at most 16 significant bits */
  664.  
  665. qmuli( a, b, c )
  666. short *a, *b, *c;
  667. {
  668. int ctr;
  669. long lt;
  670.  
  671. if( (a[E] == 0) || (b[E] == 0) )
  672.     {
  673.     qclear(c);
  674.     return;
  675.     }
  676.  
  677. qmov( b, ac3 );
  678. mulin( a, ac3 );
  679.  
  680. /* calculate sign of product */
  681. if( b[0] == a[0] )
  682.     ac3[0] = 0;
  683. else
  684.     ac3[0] = -1;
  685.  
  686. /* calculate exponent */
  687. lt = (long)ac3[E] + (long )a[E] - 32768L;
  688. if( lt > 32767 )
  689.     goto overf;
  690. ac3[E] = lt;
  691. ac3[NQ] = 0;
  692. if( normlz(ac3) )
  693.     goto underf;
  694. lt = lt - SC + 16384L;
  695. if( lt > 32767 )
  696.     goto overf;
  697. if( lt < 0 )
  698.     goto underf;
  699. ac3[E] = lt;
  700. qmov( ac3, c );
  701. return;
  702.  
  703. underf:
  704. qclear(c);
  705. return;
  706.  
  707. overf:
  708. qinfin(c);
  709. mtherr( "qmuli", OVERFLOW );
  710. return;
  711. }
  712.  
  713.  
  714.  
  715.  
  716. /*
  717. ;    Compare mantissas
  718. ;
  719. ;    short a[NQ], b[NQ];
  720. ;    cmpm( a, b );
  721. ;
  722. ;    for the mantissas:
  723. ;    returns    +1 if a > b
  724. ;         0 if a == b
  725. ;        -1 if a < b
  726. */
  727.  
  728. cmpm( a, b )
  729. register unsigned short *a, *b;
  730. {
  731. int i;
  732.  
  733. a += 2; /* skip up to mantissa area */
  734. b += 2;
  735. for( i=0; i<OMG; i++ )
  736.     {
  737.     if( *a++ != *b++ )
  738.         goto difrnt;
  739.     }
  740. return(0);
  741.  
  742. difrnt:
  743. if( *(--a) > *(--b) )
  744.     return(1);
  745. else
  746.     return(-1);
  747. }
  748.  
  749. /*
  750. ;    shift mantissa
  751. ;
  752. ;    Shifts mantissa area up or down by the number of bits
  753. ;    given by the variable SC.
  754. */
  755.  
  756. shift( x )
  757. short *x;
  758. {
  759. short *p;
  760. #if STICKY
  761. int lost;
  762. #endif
  763.  
  764. if( SC == 0 )
  765.     {
  766. #if STICKY
  767.     return(0);
  768. #else
  769.     return;
  770. #endif
  771.     }
  772.  
  773. #if STICKY
  774. lost = 0;
  775. #endif
  776. if( SC < 0 )
  777.     {
  778.     p = x + NQ;
  779.     SC = -SC;
  780.     while( SC >= 16 )
  781.         {
  782. #if STICKY
  783.         lost |= *p;
  784. #endif
  785.         shdn16(x);
  786.         SC -= 16;
  787.         }
  788.  
  789.     while( SC >= 8 )
  790.         {
  791. #if STICKY
  792.         lost |= *p & 0xff;
  793. #endif
  794.         shdn8(x);
  795.         SC -= 8;
  796.         }
  797.  
  798.     while( SC > 0 )
  799.         {
  800. #if STICKY
  801.         lost |= *p & 1;
  802. #endif
  803.         shdn1(x);
  804.         SC -= 1;
  805.         }
  806.     }
  807. else
  808.     {
  809.     while( SC >= 16 )
  810.         {
  811.         shup16(x);
  812.         SC -= 16;
  813.         }
  814.  
  815.     while( SC >= 8 )
  816.         {
  817.         shup8(x);
  818.         SC -= 8;
  819.         }
  820.  
  821.     while( SC > 0 )
  822.         {
  823.         shup1(x);
  824.         SC -= 1;
  825.         }
  826.     }
  827. #if STICKY
  828. return( lost );
  829. #endif
  830. }
  831.  
  832. /*
  833. ;    normalize
  834. ;
  835. ; shift normalizes the mantissa area pointed to by R1
  836. ; shift count (up = positive) returned in SC
  837. */
  838.  
  839. normlz(x)
  840. short x[];
  841. {
  842. register short *p;
  843.  
  844. SC = 0;
  845. p = &x[M];
  846. if( *p != 0 )
  847.     goto normdn;
  848. ++p;
  849. if( *p < 0 )
  850.     return(0);    /* already normalized */
  851.  
  852. while( *p == 0 )
  853.     {
  854.     shup16(x);
  855.     SC += 16;
  856. /* With guard word, there are NBITS+16 bits available.
  857.  * return true if all are zero.
  858.  */
  859.     if( SC > NBITS )
  860.         return(1);
  861.     }
  862.  
  863. /* see if high byte is zero */
  864. while( (*p & 0xff00) == 0 )
  865.     {
  866.     shup8(x);
  867.     SC += 8;
  868.     }
  869. /* now shift 1 bit at a time */
  870.  
  871. while( *p > 0)
  872.     {
  873.     shup1(x);
  874.     SC += 1;
  875. /*
  876.     if( SC > NBITS )
  877.         {
  878.         printf( "normlz error\n");
  879.         return(0);
  880.         }
  881. */
  882.     }
  883. return(0);
  884.  
  885. /* normalize by shifting down out of the high guard word
  886.    of the mantissa */
  887.  
  888. normdn:
  889.  
  890. if( *p & 0xff00 )
  891.     {
  892.     shdn8(x);
  893.     SC -= 8;
  894.     }
  895. while( *p != 0 )
  896.     {
  897.     shdn1(x);
  898.     SC -= 1;
  899. /*
  900.     if( SC < -NBITS )
  901.         {
  902.         printf( "low normlz error\n");
  903.         return(0);
  904.         }
  905. */
  906.     }
  907. return(0);
  908. }
  909.  
  910. /*
  911. ; Clear out entire number, including sign and exponent, pointed
  912. ; to by x
  913. ;
  914. ; short x[];
  915. ; qclear( x );
  916. */
  917.  
  918. qclear( x )
  919. register short *x;
  920. {
  921. register int i;
  922.  
  923. for( i=0; i<NQ; i++ )
  924.     *x++ = 0;
  925. }
  926.  
  927. /*
  928. ; Fill entire number, including exponent and mantissa, with
  929. ; largest possible number.
  930. */
  931.  
  932. qinfin(x)
  933. register short *x;
  934. {
  935. register int i;
  936.  
  937. ++x;    /* skip over the sign */
  938. *x++ = 0x7fff;
  939. *x++ = 0;
  940. for( i=0; i<OMG-1; i++ )
  941.     *x++ = 0xffff;
  942. }
  943.  
  944. /* normalization program */
  945. qnrmlz(x)
  946. short *x;    /* x is the memory address of a short */
  947. {
  948.  
  949. qmovz( x, ac1 );
  950. normlz( ac1 );    /* shift normalize the mantissa */
  951. ac1[E] -= SC;    /* subtract the shift count from the exponent */
  952. qmov( ac1, x );
  953. }
  954.  
  955.  
  956.  
  957.  
  958. /*
  959. ; Convert IEEE double precision to Q type
  960. ;    double d;
  961. ;    short q[N+2];
  962. ;    etoq( &d, q );
  963. */
  964.  
  965. etoq( e, y )
  966. short *e, *y;
  967. {
  968. #ifdef DEC
  969. dtoq(e,y);
  970. #else
  971. register short r;
  972. register short *p;
  973. short yy[NQ+1];
  974. int denorm;
  975.  
  976.  
  977. denorm = 0;    /* flag if denormalized number */
  978. qclear(yy);
  979. yy[NQ] = 0;
  980.  
  981. #ifdef MIEEE
  982.  
  983. #endif
  984.  
  985.  
  986. #ifdef IBMPC
  987. e += 3;
  988. #endif
  989.  
  990. /*
  991. r = *e & 0x7fff;
  992. if( r == 0 )
  993.     return;
  994. */
  995.  
  996. r = *e;
  997. yy[0] = 0;
  998. if( r < 0 )
  999.     yy[0] = -1;
  1000.  
  1001. yy[M] = (r & 0x0f) | 0x10;
  1002. r &= ~0x800f;    /* strip sign and 4 mantissa bits */
  1003. r >>= 4;
  1004. /* If zero exponent, then the mantissa is denormalized.
  1005.  * So, take back the understood high mantissa bit. */ 
  1006. if( r == 0 )
  1007.     {
  1008.     denorm = 1;
  1009.     yy[M] &= ~0x10;
  1010.     }
  1011. r += 040001 - 01777;
  1012. yy[E] = r;
  1013. p = &yy[M+1];
  1014. #ifdef MIEEE
  1015. ++e;
  1016. *p++ = *e++;
  1017. *p++ = *e++;
  1018. *p++ = *e++;
  1019. #endif
  1020. #ifdef IBMPC
  1021. *p++ = *(--e);
  1022. *p++ = *(--e);
  1023. *p++ = *(--e);
  1024. #endif
  1025. SC = -5;
  1026. shift(yy);
  1027. if( denorm )
  1028.     { /* if zero exponent, then normalize the mantissa */
  1029.     if( normlz( yy ) )
  1030.         qclear(yy);
  1031.     else
  1032.         yy[E] -= SC-1;
  1033.     }
  1034. qmov( yy, y );
  1035. /* not DEC */
  1036. #endif
  1037. }
  1038.  
  1039. /*
  1040. ; Q type to IEEE double precision
  1041. ;    double d;
  1042. ;    short q[N+2];
  1043. ;    qtoe( q, &d );
  1044. */
  1045. qtoe( x, e )
  1046. short *x, *e;
  1047. {
  1048. #ifdef DEC
  1049. qtod(x,e);
  1050. #else
  1051. short i, j, k;
  1052. register short *p;
  1053.  
  1054.  
  1055. #ifdef MIEEE
  1056.  
  1057. #endif
  1058. #ifdef IBMPC
  1059. e += 3;
  1060. #endif
  1061.  
  1062.  
  1063. *e = 0;    /* output high order */
  1064. p = &ac1[0];
  1065. qmovz( x, p );
  1066. if( *p++ < 0 )
  1067.     *e = 0x8000;    /* output sign bit */
  1068.  
  1069. normlz(ac1);
  1070. *p -= SC;
  1071.  
  1072. i = *p++ -(040001 - 01777);    /* adjust exponent for offsets */
  1073.  
  1074. /* Handle denormalized small numbers.
  1075.  * The tiniest number represented appears to be 2**-1074.
  1076.  */
  1077. if( i <= 0 )
  1078.     {
  1079.     if( i > -53 )
  1080.         {
  1081.         SC = i - 1;
  1082.         shift( ac1 );
  1083.         i = 0;
  1084.         }
  1085.     else
  1086.         {
  1087. /*ozero:*/
  1088. #ifdef MIEEE
  1089.         ++e;
  1090.         *e++ = 0;
  1091.         *e++ = 0;
  1092.         *e++ = 0;
  1093. #endif
  1094. #ifdef IBMPC
  1095.         *(--e) = 0;
  1096.         *(--e) = 0;
  1097.         *(--e) = 0;
  1098. #endif
  1099.         return;
  1100.         }
  1101.     }
  1102.  
  1103. /* round off to nearest or even */
  1104. k = ac1[M+4];
  1105. if( (k & 0x400) != 0 )
  1106.     {
  1107.     if( (k & 0x07ff) == 0x400 )
  1108.         {
  1109.         if( (k & 0x800) != 0 )
  1110.             {
  1111.         /* check all less significant bits */
  1112.             for( j=M+5; j<=NQ; j++ )
  1113.                 {
  1114.                 if( ac1[j] != 0 )
  1115.                     goto yesrnd;
  1116.                 }
  1117.             }
  1118.         goto nornd;
  1119.         }
  1120. yesrnd:
  1121.     qclear( ac2 );
  1122.     ac2[NQ] = 0;
  1123.     ac2[M+4] = 0x800;
  1124.     addm( ac2, ac1 );
  1125.     if( ac1[2] )
  1126.         {
  1127.         shdn1(ac1);
  1128.         i += 1;
  1129.         }
  1130.     if( (i == 0) && (ac1[3] & 0x8000) )
  1131.         i += 1;
  1132.     }
  1133.  
  1134. nornd:
  1135.  
  1136. if( i > 2047 )
  1137.     {    /* Saturate at largest number less than infinity. */
  1138.     mtherr( "qtoe", OVERFLOW );
  1139.     *e |= 0x7fef;
  1140. #ifdef MIEEE
  1141.     ++e;
  1142.     *e++ = 0xffff;
  1143.     *e++ = 0xffff;
  1144.     *e++ = 0xffff;
  1145. #endif
  1146. #ifdef IBMPC
  1147.     *(--e) = 0xffff;
  1148.     *(--e) = 0xffff;
  1149.     *(--e) = 0xffff;
  1150. #endif
  1151.     return;
  1152.     }
  1153.  
  1154.  
  1155. i <<= 4;
  1156. SC = 5;
  1157. shift( ac1 );
  1158. i |= *p++ & 0x0f;    /* *p = ac1[M] */
  1159. *e |= i;    /* high order output already has sign bit set */
  1160. #ifdef MIEEE
  1161. ++e;
  1162. *e++ = *p++;
  1163. *e++ = *p++;
  1164. *e++ = *p++;
  1165. #endif
  1166. #ifdef IBMPC
  1167. *(--e) = *p++;
  1168. *(--e) = *p++;
  1169. *(--e) = *p;
  1170. #endif
  1171. /* not DEC */
  1172. #endif
  1173. }
  1174.  
  1175.  
  1176.  
  1177. /*                        qtoasc.c    */
  1178. /* Convert q type number to ASCII string */
  1179.  
  1180. /*include "qhead.h"*/
  1181.  
  1182. #define NTEN 12
  1183.  
  1184. #if NQ < 13
  1185. #define NTT 12
  1186. static short qtens[NTEN+1][NTT] = {
  1187. /* 4096 */
  1188. {0x0000,0x7527,0x0000,0xc460,0x5202,0x8a20,0x979a,0xc94c,
  1189. 0x153f,0x804a,0x4a92,0x6576,},
  1190. /* 2048 */
  1191. {0x0000,0x5a94,0x0000,0x9e8b,0x3b5d,0xc53d,0x5de4,0xa74d,
  1192. 0x28ce,0x329a,0xce52,0x6a32,},
  1193. /* 1024 */
  1194. {0x0000,0x4d4a,0x0000,0xc976,0x7586,0x8175,0x0c17,0x650d,
  1195. 0x3d28,0xf18b,0x50ce,0x526c,},
  1196. /* 512 */
  1197. {0x0000,0x46a5,0x0000,0xe319,0xa0ae,0xa60e,0x91c6,0xcc65,
  1198. 0x5c54,0xbc50,0x58f8,0x9c66,},
  1199. /* 256 */
  1200. {0x0000,0x4353,0x0000,0xaa7e,0xebfb,0x9df9,0xde8d,0xddbb,
  1201. 0x901b,0x98fe,0xeab7,0x851e,},
  1202. /* 128 */
  1203. {0x0000,0x41aa,0x0000,0x93ba,0x47c9,0x80e9,0x8cdf,0xc66f,
  1204. 0x336c,0x36b1,0x0137,0x0235,},
  1205. /* 64 */
  1206. {0x0000,0x40d5,0x0000,0xc278,0x1f49,0xffcf,0xa6d5,0x3cbf,
  1207. 0x6b71,0xc76b,0x25fb,0x50f8,},
  1208. /* 32 */
  1209. {0000000,0040153,0000000,0116705,0126650,0025560,
  1210. 0132635,0170040,0000000,0000000,0000000,0000000,},
  1211. /* 16 */
  1212. {0000000,0040066,0000000,0107033,0144677,0002000,
  1213. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1214. /* 8 */
  1215. {0000000,0040033,0000000,0137274,0020000,0000000,
  1216. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1217. /* 4 */
  1218. {0000000,0040016,0000000,0116100,0000000,0000000,
  1219. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1220. /* 2 */
  1221. {0000000,0040007,0000000,0144000,0000000,0000000,
  1222. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1223. /* 1 */
  1224. {0000000,0040004,0000000,0120000,0000000,0000000,
  1225. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1226. };
  1227.  
  1228. static short qmtens[NTEN+1][NTT] = {
  1229. /* -4096 */
  1230. {0x0000,0x0ada,0x0000,0xa6dd,0x04c8,0xd2ce,0x9fde,0x2de3,
  1231. 0x8123,0xa1c3,0xcffc,0x2030,},
  1232. /* -2048 */
  1233. {0x0000,0x256d,0x0000,0xceae,0x534f,0x3436,0x2de4,0x4925,
  1234. 0x12d4,0xf2ea,0xd2cb,0x8264,},
  1235. /* -1024 */
  1236. {0x0000,0x32b7,0x0000,0xa2a6,0x82a5,0xda57,0xc0bd,0x87a6,
  1237. 0x0158,0x6bd3,0xf698,0xf53f,},
  1238. /* -512 */
  1239. {0x0000,0x395c,0x0000,0x9049,0xee32,0xdb23,0xd21c,0x7132,
  1240. 0xd332,0xe3f2,0x04d4,0xe731,},
  1241. /* -256 */
  1242. {0x0000,0x3cae,0x0000,0xc031,0x4325,0x637a,0x1939,0xfa91,
  1243. 0x1155,0xfefb,0x5308,0xa23e,},
  1244. /* -128 */
  1245. {0x0000,0x3e57,0x0000,0xddd0,0x467c,0x64bc,0xe4a0,0xac7c,
  1246. 0xb3f6,0xd05d,0xdbde,0xe26d,},
  1247. /* -64 */
  1248. {0x0000,0x3f2c,0x0000,0xa87f,0xea27,0xa539,0xe9a5,0x3f23,
  1249. 0x98d7,0x47b3,0x6224,0x2a20,},
  1250. /* -32 */
  1251. {0x0000,0x3f96,0x0000,0xcfb1,0x1ead,0x4539,0x94ba,0x67de,
  1252. 0x18ed,0xa581,0x4af2,0x0b5b,},
  1253. /* -16 */
  1254. {0x0000,0x3fcb,0x0000,0xe695,0x94be,0xc44d,0xe15b,0x4c2e,
  1255. 0xbe68,0x7989,0xa9b3,0xbf71,},
  1256. /* -8 */
  1257. {0x0000,0x3fe6,0x0000,0xabcc,0x7711,0x8461,0xcefc,0xfdc2,
  1258. 0x0d2b,0x36ba,0x7c3d,0x3d4d,},
  1259. /* -4 */
  1260. {0x0000,0x3ff3,0x0000,0xd1b7,0x1758,0xe219,0x652b,0xd3c3,
  1261. 0x6113,0x404e,0xa4a8,0xc155,},
  1262. /* -2 */
  1263. {0x0000,0x3ffa,0x0000,0xa3d7,0x0a3d,0x70a3,0xd70a,0x3d70,
  1264. 0xa3d7,0x0a3d,0x70a3,0xd70a,},
  1265. /* -1 */
  1266. {0x0000,0x3ffd,0x0000,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,
  1267. 0xcccc,0xcccc,0xcccc,0xcccd,},
  1268. };
  1269.  
  1270. #else
  1271.  
  1272. #define NTT 24
  1273. static short qtens[NTEN+1][NTT] = {
  1274. /* 4096 */
  1275. {0x0000,0x7527,0x0000,0xc460,0x5202,0x8a20,0x979a,0xc94c,
  1276. 0x153f,0x804a,0x4a92,0x6576,0x1fb2,0x444e,0x2267,0xdd5c,
  1277. 0xf7c9,0x45f2,0x2a3f,0xfff0,0xd1d2,0xe184,0x69b1,0xea39,},
  1278. /* 2048 */
  1279. {0x0000,0x5a94,0x0000,0x9e8b,0x3b5d,0xc53d,0x5de4,0xa74d,
  1280. 0x28ce,0x329a,0xce52,0x6a31,0x97bb,0xebe3,0x034f,0x7715,
  1281. 0x4ce2,0xbcba,0x1964,0x8b21,0xc11e,0xb962,0xb1b6,0x1b94,},
  1282. /* 1024 */
  1283. {0x0000,0x4d4a,0x0000,0xc976,0x7586,0x8175,0x0c17,0x650d,
  1284. 0x3d28,0xf18b,0x50ce,0x526b,0x9882,0x7524,0x9b0f,0xd6f4,
  1285. 0xb6d2,0x7bd1,0xc61c,0x2530,0x69a5,0xc329,0xf9af,0x4a9c,},
  1286. /* 512 */
  1287. {0x0000,0x46a5,0x0000,0xe319,0xa0ae,0xa60e,0x91c6,0xcc65,
  1288. 0x5c54,0xbc50,0x58f8,0x9c65,0x8398,0x1d13,0x4cba,0x422d,
  1289. 0x38ea,0x3584,0xcde4,0x0b9a,0x1d7d,0x634d,0xf2d8,0x74a2,},
  1290. /* 256 */
  1291. {0x0000,0x4353,0x0000,0xaa7e,0xebfb,0x9df9,0xde8d,0xddbb,
  1292. 0x901b,0x98fe,0xeab7,0x851e,0x4cbf,0x3de2,0xf98a,0xae78,
  1293. 0x0c7f,0xea81,0xc788,0x5a6b,0x43a2,0xa6c4,0x95b8,0xccd6,},
  1294. /* 128 */
  1295. {0x0000,0x41aa,0x0000,0x93ba,0x47c9,0x80e9,0x8cdf,0xc66f,
  1296. 0x336c,0x36b1,0x0137,0x0234,0xf3fd,0x7b08,0xdd39,0x0bc3,
  1297. 0xc54e,0x3f40,0xf7e6,0x424b,0xa54f,0x8040,0x0000,0x0000,},
  1298. /* 64 */
  1299. {0x0000,0x40d5,0x0000,0xc278,0x1f49,0xffcf,0xa6d5,0x3cbf,
  1300. 0x6b71,0xc76b,0x25fb,0x50f8,0x0800,0x0000,0x0000,0x0000,
  1301. 0x0000,0x0000,0x0000,0x0000,0x0000,0x0000,0x0000,0x0000,},
  1302. /* 32 */
  1303. {0000000,0040153,0000000,0116705,0126650,0025560,
  1304. 0132635,0170040,0000000,0000000,0000000,0000000,
  1305. 0000000,0000000,0000000,0000000,0000000,0000000,
  1306. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1307. /* 16 */
  1308. {0000000,0040066,0000000,0107033,0144677,0002000,
  1309. 0000000,0000000,0000000,0000000,0000000,0000000,
  1310. 0000000,0000000,0000000,0000000,0000000,0000000,
  1311. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1312. /* 8 */
  1313. {0000000,0040033,0000000,0137274,0020000,0000000,
  1314. 0000000,0000000,0000000,0000000,0000000,0000000,
  1315. 0000000,0000000,0000000,0000000,0000000,0000000,
  1316. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1317. /* 4 */
  1318. {0000000,0040016,0000000,0116100,0000000,0000000,
  1319. 0000000,0000000,0000000,0000000,0000000,0000000,
  1320. 0000000,0000000,0000000,0000000,0000000,0000000,
  1321. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1322. /* 2 */
  1323. {0000000,0040007,0000000,0144000,0000000,0000000,
  1324. 0000000,0000000,0000000,0000000,0000000,0000000,
  1325. 0000000,0000000,0000000,0000000,0000000,0000000,
  1326. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1327. /* 1 */
  1328. {0000000,0040004,0000000,0120000,0000000,0000000,
  1329. 0000000,0000000,0000000,0000000,0000000,0000000,
  1330. 0000000,0000000,0000000,0000000,0000000,0000000,
  1331. 0000000,0000000,0000000,0000000,0000000,0000000,},
  1332. };
  1333.  
  1334. static short qmtens[NTEN+1][NTT] = {
  1335. /* -4096 */
  1336. {0x0000,0x0ada,0x0000,0xa6dd,0x04c8,0xd2ce,0x9fde,0x2de3,
  1337. 0x8123,0xa1c3,0xcffc,0x2030,0x5d02,0x44e0,0x91ba,0x5e2d,
  1338. 0x7403,0x972f,0x6f2b,0x04e6,0x63ac,0x16d5,0x1215,0x9834,},
  1339. /* -2048 */
  1340. {0x0000,0x256d,0x0000,0xceae,0x534f,0x3436,0x2de4,0x4925,
  1341. 0x12d4,0xf2ea,0xd2cb,0x8263,0xca5c,0xbc77,0x4bd9,0x71aa,
  1342. 0xd590,0x46c7,0x4249,0xdddd,0xc0df,0xbb0d,0x6f2c,0x0584,},
  1343. /* -1024 */
  1344. {0x0000,0x32b7,0x0000,0xa2a6,0x82a5,0xda57,0xc0bd,0x87a6,
  1345. 0x0158,0x6bd3,0xf698,0xf53e,0x94d1,0xb235,0x7c32,0xc0ea,
  1346. 0xff37,0x55a2,0xddcd,0x5191,0xa70f,0x9e33,0x9aa9,0x6270,},
  1347. /* -512 */
  1348. {0x0000,0x395c,0x0000,0x9049,0xee32,0xdb23,0xd21c,0x7132,
  1349. 0xd332,0xe3f2,0x04d4,0xe731,0x7d62,0x209b,0x6a93,0xd4c9,
  1350. 0x4a9d,0xa069,0x3e0c,0xfefd,0xe89b,0x01e5,0x1068,0xe0e4,},
  1351. /* -256 */
  1352. {0x0000,0x3cae,0x0000,0xc031,0x4325,0x637a,0x1939,0xfa91,
  1353. 0x1155,0xfefb,0x5308,0xa23e,0x2ed2,0x7766,0xe8cc,0x9b03,
  1354. 0x5377,0x08b1,0x648f,0x7a73,0x0e52,0xe6e7,0x4360,0x55f4,},
  1355. /* -128 */
  1356. {0x0000,0x3e57,0x0000,0xddd0,0x467c,0x64bc,0xe4a0,0xac7c,
  1357. 0xb3f6,0xd05d,0xdbde,0xe26c,0xa606,0x3461,0xfffa,0x4ed7,
  1358. 0x75fc,0x49f2,0x7952,0xf2f3,0xa45f,0x9009,0xf3c9,0xee49,},
  1359. /* -64 */
  1360. {0x0000,0x3f2c,0x0000,0xa87f,0xea27,0xa539,0xe9a5,0x3f23,
  1361. 0x98d7,0x47b3,0x6224,0x2a1f,0xee40,0xd90a,0xab31,0x0e12,
  1362. 0x8b5d,0x938c,0xfb3f,0x6f20,0x3147,0x025d,0x1129,0xe492,},
  1363. /* -32 */
  1364. {0x0000,0x3f96,0x0000,0xcfb1,0x1ead,0x4539,0x94ba,0x67de,
  1365. 0x18ed,0xa581,0x4af2,0x0b5b,0x1aa0,0x28cc,0xd99e,0x59e3,
  1366. 0x38e3,0x87ad,0x8e28,0x5676,0x0892,0x1621,0x96af,0x3c7f,},
  1367. /* -16 */
  1368. {0x0000,0x3fcb,0x0000,0xe695,0x94be,0xc44d,0xe15b,0x4c2e,
  1369. 0xbe68,0x7989,0xa9b3,0xbf71,0x6c1a,0xdd27,0xf085,0x23cc,
  1370. 0xd348,0x4db6,0x70aa,0xd6e0,0x12cf,0x7fa7,0xcf42,0x8583,},
  1371. /* -8 */
  1372. {0x0000,0x3fe6,0x0000,0xabcc,0x7711,0x8461,0xcefc,0xfdc2,
  1373. 0x0d2b,0x36ba,0x7c3d,0x3d4d,0x3d75,0x8161,0x697c,0x7068,
  1374. 0xf3b4,0x6d2f,0x8350,0x5705,0x755f,0xd37a,0x3b04,0x8dd9,},
  1375. /* -4 */
  1376. {0x0000,0x3ff3,0x0000,0xd1b7,0x1758,0xe219,0x652b,0xd3c3,
  1377. 0x6113,0x404e,0xa4a8,0xc154,0xc985,0xf06f,0x6944,0x6738,
  1378. 0x1d7d,0xbf48,0x7fcb,0x923a,0x29c7,0x79a6,0xb50b,0x0f28,},
  1379. /* -2 */
  1380. {0x0000,0x3ffa,0x0000,0xa3d7,0x0a3d,0x70a3,0xd70a,0x3d70,
  1381. 0xa3d7,0x0a3d,0x70a3,0xd70a,0x3d70,0xa3d7,0x0a3d,0x70a3,
  1382. 0xd70a,0x3d70,0xa3d7,0x0a3d,0x70a3,0xd70a,0x3d70,0xa3d7,},
  1383. /* -1 */
  1384. {0x0000,0x3ffd,0x0000,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,
  1385. 0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,
  1386. 0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccc,0xcccd,},
  1387. };
  1388. #endif
  1389.  
  1390.  
  1391.  
  1392. qtoasc( q, string, ndigs )
  1393. short q[];
  1394. char *string;
  1395. int ndigs;
  1396. {
  1397. long digit;
  1398. short x[NTT], xt[NTT];
  1399. short *p, *r, *ten, *tenth;
  1400. short sign;
  1401. int i, k, expon;
  1402. char *s, *ss;
  1403.  
  1404. qmov( q, x );
  1405. sign = x[0];
  1406. x[0] = 0;
  1407. expon = 0;
  1408. ten = &qtens[NTEN][0];
  1409.  
  1410. i = qcmp( qone, x );
  1411. if( i == 0 )
  1412.     goto isone;
  1413. if( x[1] == 0 )
  1414.     {
  1415.     qclear( x );
  1416.     goto isone;
  1417.     }
  1418.  
  1419. if( i < 0 )
  1420.     {
  1421.     k = 4096;
  1422.     p = &qtens[0][0];
  1423.     qmov( qone, ac4 );
  1424.     qmov( x, xt );
  1425.     while( qcmp( ten, x ) <= 0 )
  1426.         {
  1427.         if( qcmp( p, xt ) <= 0 )
  1428.             {
  1429.             qdiv( p, xt, xt );
  1430.             qmul( p, ac4, ac4 );
  1431.             expon += k;
  1432.             }
  1433.         k >>= 1;
  1434.         if( k == 0 )
  1435.             break;
  1436.         p += NTT;
  1437.         }
  1438.     qdiv( ac4, x, x );
  1439.     }
  1440. else
  1441.     {
  1442.     k = -4096;
  1443.     p = &qmtens[0][0];
  1444.     r = &qtens[0][0];
  1445.     tenth = &qmtens[NTEN][0];
  1446.     while( qcmp( tenth, x ) > 0 )
  1447.         {
  1448.         if( qcmp( p, x ) >= 0 )
  1449.             {
  1450.             qmul( r, x, x );
  1451.             expon += k;
  1452.             }
  1453.         k /= 2;
  1454. /* Prevent infinite loop due to arithmetic error: */
  1455.         if( k == 0 )
  1456.             break;
  1457.         p += NTT;
  1458.         r += NTT;
  1459.         }
  1460.     qmuli( ten, x, x );
  1461.     expon -= 1;
  1462.     }
  1463.  
  1464. isone:
  1465. qifrac( x, &digit, x );
  1466. /* The following check handles numbers very close to 10**(2**n)
  1467.  * when there is a mistake due to arithmetic error.
  1468.  */
  1469. if( digit >= 10 )
  1470.     {
  1471.     qdiv( ten, x, x );
  1472.     expon += 1;
  1473.     digit = 1;
  1474.     }
  1475. s = string;
  1476. if( sign < 0 )
  1477.     *s++ = '-';
  1478. else
  1479.     *s++ = ' ';
  1480. *s++ = (char )digit + 060;
  1481. *s++ = '.';
  1482. if( ndigs < 0 )
  1483.     ndigs = 0;
  1484. if( ndigs > NDEC )
  1485.     ndigs = NDEC;
  1486. for( k=0; k<ndigs; k++ )
  1487.     {
  1488.     qmuli( ten, x, x );
  1489.     qifrac( x, &digit, x );
  1490.     *s++ = (char )digit + 060;
  1491.     }
  1492.  
  1493. *s = '\0';    /* mark end of string */
  1494. ss = s;
  1495.  
  1496. /* round off the ASCII string */
  1497.  
  1498. qmuli( ten, x, x );
  1499. qifrac( x, &digit, x );
  1500. if( digit > 4 )
  1501.     {
  1502. /* Check for critical rounding case */
  1503.     if( digit == 5 )
  1504.         {
  1505.         if( qcmp( x, qzero ) != 0 )
  1506.             goto roun;    /* round to nearest */
  1507.         if( (*(s-1) & 1) == 0 )
  1508.             goto doexp;    /* round to even */
  1509.         }
  1510. roun:
  1511.     --s;
  1512.     k = *s & 0x7f;
  1513. /* Carry out to most significant digit? */
  1514.     if( k == '.' )
  1515.         {
  1516.         --s;
  1517.         k = *s & 0x7f;
  1518.         k += 1;
  1519.         *s = k;
  1520. /* Most significant digit rounds to 10? */
  1521.         if( k > '9' )
  1522.             {
  1523.             *s = '1';
  1524.             expon += 1;
  1525.             }
  1526.         goto doexp;
  1527.         }
  1528. /* Round up and carry out from less significant digits. */
  1529.     k += 1;
  1530.     *s = k;
  1531.     if( k > '9' )
  1532.         {
  1533.         *s = '0';
  1534.         goto roun;
  1535.         }
  1536.     }
  1537.  
  1538. doexp:
  1539.  
  1540. sprintf( ss, "E%d\0", expon );
  1541. }
  1542.  
  1543.  
  1544. /*  short a[NQ], b[NQ];
  1545.  * qcmp( a, b )
  1546.  *
  1547.  *  returns +1 if a > b
  1548.  *           0 if a == b
  1549.  *          -1 if a < b
  1550.  */
  1551.  
  1552. qcmp( p, q )
  1553. register unsigned short *p, *q;
  1554. {
  1555. unsigned short r[NQ];
  1556. register int i;
  1557. int msign;
  1558.  
  1559. if( ( *(p+1) <= NBITS)  && ( *(q+1) <= NBITS ) )
  1560.     {
  1561.     qsub( q, p, r );
  1562.     if( r[E] == 0 )
  1563.         return( 0 );
  1564.     if( r[0] == 0 )
  1565.         return( 1 );
  1566.     return( -1 );
  1567.     }
  1568.  
  1569. if( *p != *q )
  1570.     { /* the signs are different */
  1571.     if( *p == 0 )
  1572.         return( 1 );
  1573.     else
  1574.         return( -1 );
  1575.     }
  1576.  
  1577. /* both are the same sign */
  1578. if( *p == 0 )
  1579.     msign = 1;
  1580. else
  1581.     msign = -1;
  1582.     
  1583. i = NQ;
  1584. do
  1585.     {
  1586.     if( *p++ != *q++ )
  1587.         {
  1588.         goto diff;
  1589.         }
  1590.     }
  1591. while( --i > 0 );
  1592.  
  1593. return(0);    /* equality */
  1594.  
  1595.  
  1596.  
  1597. diff:
  1598.  
  1599. if( *(--p) > *(--q) )
  1600.     return( msign );        /* p is bigger */
  1601. else
  1602.     return( -msign );    /* p is littler */
  1603. }
  1604.  
  1605.  
  1606.  
  1607. /*
  1608. ;                                ASCTOQ
  1609. ;        ASCTOQ.MAC        LATEST REV: 11 JAN 84
  1610. ;                    SLM, 3 JAN 78
  1611. ;
  1612. ;    Convert ASCII string to quadruple precision floating point
  1613. ;
  1614. ;        Numeric input is free field decimal number
  1615. ;        with max of 15 digits with or without 
  1616. ;        decimal point entered as ASCII from teletype.
  1617. ;    Entering E after the number followed by a second
  1618. ;    number causes the second number to be interpreted
  1619. ;    as a power of 10 to be multiplied by the first number
  1620. ;    (i.e., "scientific" notation).
  1621. ;
  1622. ;    Usage:
  1623. ;        asctoq( string, q );
  1624. */
  1625.  
  1626. /*if NQ < 13*/
  1627. /*static short qten[NTT] = {0,040004,0,0120000,0,0,0,0,0,0,0,0};*/
  1628.  
  1629. static short qten[24] = {
  1630.  0,040004,0,0120000,0,0,0,0,0,0,0,0,
  1631.  0,0,0,0,0,0,0,0,0,0,0,0
  1632. };
  1633.  
  1634.  
  1635. asctoq( s, y )
  1636. char *s;
  1637. short *y;
  1638. {
  1639. short yy[NQ+1], qt[NQ];
  1640. int esign, nsign, decflg, sgnflg, nexp, exp, prec;
  1641. short *p;
  1642.  
  1643.  
  1644. nsign = 0;
  1645. decflg = 0;
  1646. sgnflg = 0;
  1647. nexp = 0;
  1648. exp = 0;
  1649. prec = 0;
  1650. qclear( yy );
  1651. yy[NQ] = 0;
  1652.  
  1653. nxtcom:
  1654. if( (*s >= '0') && (*s <= '9') )
  1655. {
  1656. if( (prec == 0) && (decflg == 0) && (*s == '0') )
  1657.     goto donchr;
  1658. if( prec < NDEC )
  1659.     {
  1660.     if( decflg )
  1661.         nexp += 1;    /* count digits after decimal point */
  1662.     shup1( yy );    /* multiply current number by 10 */
  1663.     qmovz( yy, ac2 );
  1664.     shup1( ac2 );
  1665.     shup1( ac2 );
  1666.     addm( ac2, yy );
  1667.     qclear( ac2 );
  1668.     ac2[OMG+1] = *s - '0';
  1669.     addm( ac2, yy );
  1670.     }
  1671. prec += 1;
  1672. goto donchr;
  1673. }
  1674.  
  1675. switch( *s )
  1676. {
  1677. case ' ':
  1678.     break;
  1679. case 'E':
  1680. case 'e':
  1681.     goto expnt;
  1682. case '.':    /* decimal point */
  1683.     if( decflg )
  1684.         goto error;
  1685.     ++decflg;
  1686.     break;
  1687. case '-':
  1688.     nsign = -1;
  1689. case '+':
  1690.     if( sgnflg )
  1691.         goto error;
  1692.     ++sgnflg;
  1693.     break;
  1694. case '\0':
  1695. /* For Microware OS-9 operating system: */
  1696. #ifndef OSK
  1697. case '\n':
  1698. #endif
  1699. case '\r':
  1700.     goto daldone;
  1701. default:
  1702. error:
  1703.     printf( "asctoq conversion error\n" );
  1704.     qclear(y);
  1705.     return;
  1706. }
  1707. donchr:
  1708. ++s;
  1709. goto nxtcom;
  1710.  
  1711.  
  1712. /*        EXPONENT INTERPRETATION */
  1713. expnt:
  1714.  
  1715. esign = 1;
  1716. exp = 0;
  1717. ++s;
  1718. /* check for + or - */
  1719. if( *s == '-' )
  1720.     {
  1721.     esign = -1;
  1722.     ++s;
  1723.     }
  1724. if( *s == '+' )
  1725.     ++s;
  1726. while( (*s >= '0') && (*s <= '9') )
  1727.     {
  1728.     exp *= 10;
  1729.     exp += *s++ - '0';
  1730.     }
  1731. if( esign < 0 )
  1732.     exp = -exp;
  1733.  
  1734. daldone:
  1735. nexp = exp - nexp;
  1736.  
  1737. if( normlz(yy) )
  1738.     {
  1739.     qclear(y);
  1740.     return;
  1741.     }
  1742.  
  1743. yy[E] = 040000 + NBITS - SC;
  1744. qmov( yy, y );
  1745. y[0] = nsign;
  1746.  
  1747. /* multiply or divide by 10**NEXP */
  1748. if( nexp == 0 )
  1749.     goto aexit;
  1750. esign = 0;
  1751. if( nexp < 0 )
  1752.     {
  1753.     esign = -1;
  1754.     nexp = -nexp;
  1755.     }
  1756.  
  1757. p = &qtens[NTEN][0];
  1758. exp = 1;
  1759. qmov( qone, qt );
  1760.  
  1761. do
  1762.     {
  1763.     if( exp & nexp )
  1764.         qmul( p, qt, qt );
  1765.     exp <<= 1;
  1766.     p -= NQ;
  1767.     }
  1768. while( exp < 8192 );
  1769.  
  1770. if( esign < 0 )
  1771.     qdiv( qt, y, y );
  1772. else
  1773.     qmul( qt, y, y );
  1774. aexit:
  1775. return;
  1776. }
  1777.  
  1778.