home *** CD-ROM | disk | FTP | other *** search
/ Simtel MSDOS - Coast to Coast / simteldosarchivecoasttocoast2.iso / astrnomy / de118i.zip / IEEE.C < prev    next >
C/C++ Source or Header  |  1993-02-28  |  53KB  |  3,118 lines
  1. /*                            ieee.c
  2.  *
  3.  *    Extended precision IEEE binary floating point arithmetic routines
  4.  *
  5.  * Numbers are stored in C language as arrays of 16-bit unsigned
  6.  * short integers.  The arguments of the routines are pointers to
  7.  * the arrays.
  8.  *
  9.  *
  10.  * External e type data structure, simulates Intel 8087 chip
  11.  * temporary real format but possibly with a larger significand:
  12.  *
  13.  *    NE-1 significand words    (least significant word first,
  14.  *                 most significant bit is normally set)
  15.  *    exponent        (value = EXONE for 1.0,
  16.  *                top bit is the sign)
  17.  *
  18.  *
  19.  * Internal data structure of a number (a "word" is 16 bits):
  20.  *
  21.  * ei[0]    sign word    (0 for positive, 0xffff for negative)
  22.  * ei[1]    biased exponent    (value = EXONE for the number 1.0)
  23.  * ei[2]    high guard word    (always zero after normalization)
  24.  * ei[3]
  25.  * to ei[NI-2]    significand    (NI-4 significand words,
  26.  *                 most significant word first,
  27.  *                 most significant bit is set)
  28.  * ei[NI-1]    low guard word    (0x8000 bit is rounding place)
  29.  *
  30.  *
  31.  *
  32.  *        Routines for external format numbers
  33.  *
  34.  *    asctoe( string, e )    ASCII string to extended double e type
  35.  *    asctoe64( string, &d )    ASCII string to long double
  36.  *    asctoe53( string, &d )    ASCII string to double
  37.  *    asctoe24( string, &f )    ASCII string to single
  38.  *    asctoeg( string, e, prec ) ASCII string to specified precision
  39.  *    e24toe( &f, e )        IEEE single precision to e type
  40.  *    e53toe( &d, e )        IEEE double precision to e type
  41.  *    e64toe( &d, e )        IEEE long double precision to e type
  42.  *    eabs(e)            absolute value
  43.  *    eadd( a, b, c )        c = b + a
  44.  *    eclear(e)        e = 0
  45.  *    ecmp( a, b )        compare a to b, return 1, 0, or -1
  46.  *    ediv( a, b, c )        c = b / a
  47.  *    efloor( a, b )        truncate to integer, toward -infinity
  48.  *    efrexp( a, exp, s )    extract exponent and significand
  49.  *    eifrac( e, &l, frac )   e to long integer and e type fraction
  50.  *    euifrac( e, &l, frac )  e to unsigned long integer and e type fraction
  51.  *    einfin( e )        set e to infinity, leaving its sign alone
  52.  *    eldexp( a, n, b )    multiply by 2**n
  53.  *    emov( a, b )        b = a
  54.  *    emul( a, b, c )        c = b * a
  55.  *    eneg(e)            e = -e
  56.  *    eround( a, b )        b = nearest integer value to a
  57.  *    esub( a, b, c )        c = b - a
  58.  *    e24toasc( &f, str, n )    single to ASCII string, n digits after decimal
  59.  *    e53toasc( &d, str, n )    double to ASCII string, n digits after decimal
  60.  *    e64toasc( &d, str, n )    long double to ASCII string
  61.  *    etoasc( e, str, n )    e to ASCII string, n digits after decimal
  62.  *    etoe24( e, &f )        convert e type to IEEE single precision
  63.  *    etoe53( e, &d )        convert e type to IEEE double precision
  64.  *    etoe64( e, &d )        convert e type to IEEE long double precision
  65.  *    ltoe( &l, e )        long (32 bit) integer to e type
  66.  *    ultoe( &l, e )        unsigned long (32 bit) integer to e type
  67.  *      eisneg( e )             1 if sign bit of e != 0, else 0
  68.  *      eisinf( e )             1 if e has maximum exponent
  69.  *
  70.  *
  71.  *        Routines for internal format numbers
  72.  *
  73.  *    eaddm( ai, bi )        add significands, bi = bi + ai
  74.  *    ecleaz(ei)        ei = 0
  75.  *    ecleazs(ei)        set ei = 0 but leave its sign alone
  76.  *    ecmpm( ai, bi )        compare significands, return 1, 0, or -1
  77.  *    edivm( ai, bi )        divide  significands, bi = bi / ai
  78.  *    emdnorm(ai,l,s,exp)    normalize and round off
  79.  *    emovi( a, ai )        convert external a to internal ai
  80.  *    emovo( ai, a )        convert internal ai to external a
  81.  *    emovz( ai, bi )        bi = ai, low guard word of bi = 0
  82.  *    emulm( ai, bi )        multiply significands, bi = bi * ai
  83.  *    enormlz(ei)        left-justify the significand
  84.  *    eshdn1( ai )        shift significand and guards down 1 bit
  85.  *    eshdn8( ai )        shift down 8 bits
  86.  *    eshdn6( ai )        shift down 16 bits
  87.  *    eshift( ai, n )        shift ai n bits up (or down if n < 0)
  88.  *    eshup1( ai )        shift significand and guards up 1 bit
  89.  *    eshup8( ai )        shift up 8 bits
  90.  *    eshup6( ai )        shift up 16 bits
  91.  *    esubm( ai, bi )        subtract significands, bi = bi - ai
  92.  *
  93.  *
  94.  * The result is always normalized and rounded to NI-4 word precision
  95.  * after each arithmetic operation.
  96.  *
  97.  * Exception flags and NaNs are NOT fully supported.
  98.  * This arithmetic should never produce a NaN output, but it might
  99.  * be confused by a NaN input.
  100.  * Define INFINITY in mconf.h for support of infinity; otherwise a
  101.  * saturation arithmetic is implemented.
  102.  * Denormals are always supported here where appropriate (e.g., not
  103.  * for conversion to DEC numbers).
  104.  */
  105.  
  106. /*
  107.  * Revision history:
  108.  *
  109.  *  5 Jan 84    PDP-11 assembly language version
  110.  *  2 Mar 86    fixed bug in asctoq()
  111.  *  6 Dec 86    C language version
  112.  * 30 Aug 88    100 digit version, improved rounding
  113.  * 15 May 92    80-bit long double support
  114.  *
  115.  * Author:  S. L. Moshier.
  116.  */
  117.  
  118. #include <stdio.h>
  119. #include "ehead.h"
  120. #include "mconf.h"
  121.  
  122.  
  123. /* Control register for rounding precision.
  124.  * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
  125.  */
  126. int rndprc = NBITS;
  127. extern int rndprc;
  128.  
  129. void eaddm(), esubm(), emdnorm(), asctoeg();
  130. static void toe24(), toe53(), toe64();
  131. void eremain(), einit(), eiremain();
  132. int ecmpm(), edivm(), emulm(), eisneg(), eisinf();
  133. void emovi(), emovo(), emovz(), ecleaz(), eadd1();
  134. void etodec(), todec(), dectoe();
  135.  
  136.  
  137.  
  138.  
  139. void einit()
  140. {
  141. }
  142.  
  143. /*
  144. ; Clear out entire external format number.
  145. ;
  146. ; unsigned short x[];
  147. ; eclear( x );
  148. */
  149.  
  150. void eclear( x )
  151. register unsigned short *x;
  152. {
  153. register int i;
  154.  
  155. for( i=0; i<NE; i++ )
  156.     *x++ = 0;
  157. }
  158.  
  159.  
  160.  
  161. /* Move external format number from a to b.
  162.  *
  163.  * emov( a, b );
  164.  */
  165.  
  166. void emov( a, b )
  167. register unsigned short *a, *b;
  168. {
  169. register int i;
  170.  
  171. for( i=0; i<NE; i++ )
  172.     *b++ = *a++;
  173. }
  174.  
  175.  
  176. /*
  177. ;    Absolute value of external format number
  178. ;
  179. ;    short x[NE];
  180. ;    eabs( x );
  181. */
  182.  
  183. void eabs(x)
  184. unsigned short x[];    /* x is the memory address of a short */
  185. {
  186.  
  187. x[NE-1] &= 0x7fff; /* sign is top bit of last word of external format */
  188. }
  189.  
  190.  
  191.  
  192.  
  193. /*
  194. ;    Negate external format number
  195. ;
  196. ;    unsigned short x[NE];
  197. ;    eneg( x );
  198. */
  199.  
  200. void eneg(x)
  201. unsigned short x[];
  202. {
  203.  
  204. x[NE-1] ^= 0x8000; /* Toggle the sign bit */
  205. }
  206.  
  207.  
  208.  
  209. /* Return 1 if external format number is negative,
  210.  * else return zero.
  211.  */
  212. int eisneg(x)
  213. unsigned short x[];
  214. {
  215.  
  216. if( x[NE-1] & 0x8000 )
  217.     return( 1 );
  218. else
  219.     return( 0 );
  220. }
  221.  
  222.  
  223. /* Return 1 if external format number has maximum possible exponent,
  224.  * else return zero.
  225.  */
  226. int eisinf(x)
  227. unsigned short x[];
  228. {
  229.  
  230. if( (x[NE-1] & 0x7fff) == 0x7fff )
  231.     return( 1 );
  232. else
  233.     return( 0 );
  234. }
  235.  
  236.  
  237. /*
  238. ; Fill entire number, including exponent and significand, with
  239. ; largest possible number.  These programs implement a saturation
  240. ; value that is an ordinary, legal number.  A special value
  241. ; "infinity" may also be implemented; this would require tests
  242. ; for that value and implementation of special rules for arithmetic
  243. ; operations involving inifinity.
  244. */
  245.  
  246. void einfin(x)
  247. register unsigned short *x;
  248. {
  249. register int i;
  250.  
  251. #ifdef INFINITY
  252. for( i=0; i<NE-1; i++ )
  253.     *x++ = 0;
  254. *x |= 32767;
  255. #else
  256. for( i=0; i<NE-1; i++ )
  257.     *x++ = 0xffff;
  258. *x |= 32766;
  259. if( rndprc < NBITS )
  260.     {
  261.     if( rndprc == 64 )
  262.         {
  263.         *(x-5) = 0;
  264.         }
  265.     if( rndprc == 53 )
  266.         {
  267.         *(x-4) = 0xf800;
  268.         }
  269.     else
  270.         {
  271.         *(x-4) = 0;
  272.         *(x-3) = 0;
  273.         *(x-2) = 0xff00;
  274.         }
  275.     }
  276. #endif
  277. }
  278.  
  279.  
  280.  
  281. /* Move in external format number,
  282.  * converting it to internal format.
  283.  */
  284. void emovi( a, b )
  285. unsigned short *a, *b;
  286. {
  287. register unsigned short *p, *q;
  288. int i;
  289.  
  290. q = b;
  291. p = a + (NE-1);    /* point to last word of external number */
  292. /* get the sign bit */
  293. if( *p & 0x8000 )
  294.     *q++ = 0xffff;
  295. else
  296.     *q++ = 0;
  297. /* get the exponent */
  298. *q = *p--;
  299. *q++ &= 0x7fff;    /* delete the sign bit */
  300. #ifdef INFINITY
  301. if( (*(q-1) & 0x7fff) == 0x7fff )
  302.     {
  303.     for( i=2; i<NI; i++ )
  304.         *q++ = 0;
  305.     return;
  306.     }
  307. #endif
  308. /* clear high guard word */
  309. *q++ = 0;
  310. /* move in the significand */
  311. for( i=0; i<NE-1; i++ )
  312.     *q++ = *p--;
  313. /* clear low guard word */
  314. *q = 0;
  315. }
  316.  
  317.  
  318. /* Move internal format number out,
  319.  * converting it to external format.
  320.  */
  321. void emovo( a, b )
  322. unsigned short *a, *b;
  323. {
  324. register unsigned short *p, *q;
  325. unsigned short i;
  326.  
  327. p = a;
  328. q = b + (NE-1); /* point to output exponent */
  329. /* combine sign and exponent */
  330. i = *p++;
  331. if( i )
  332.     *q-- = *p++ | 0x8000;
  333. else
  334.     *q-- = *p++;
  335. #ifdef INFINITY
  336. if( *(p-1) == 0x7fff )
  337.     {
  338.     einfin(b);
  339.     return;
  340.     }
  341. #endif
  342. /* skip over guard word */
  343. ++p;
  344. /* move the significand */
  345. for( i=0; i<NE-1; i++ )
  346.     *q-- = *p++;
  347. }
  348.  
  349.  
  350.  
  351.  
  352. /* Clear out internal format number.
  353.  */
  354.  
  355. void ecleaz( xi )
  356. register unsigned short *xi;
  357. {
  358. register int i;
  359.  
  360. for( i=0; i<NI; i++ )
  361.     *xi++ = 0;
  362. }
  363.  
  364. /* same, but don't touch the sign. */
  365.  
  366. void ecleazs( xi )
  367. register unsigned short *xi;
  368. {
  369. register int i;
  370.  
  371. ++xi;
  372. for(i=0; i<NI-1; i++)
  373.     *xi++ = 0;
  374. }
  375.  
  376.  
  377.  
  378.  
  379. /* Move internal format number from a to b.
  380.  */
  381. void emovz( a, b )
  382. register unsigned short *a, *b;
  383. {
  384. register int i;
  385.  
  386. for( i=0; i<NI-1; i++ )
  387.     *b++ = *a++;
  388. /* clear low guard word */
  389. *b = 0;
  390. }
  391.  
  392.  
  393. /*
  394. ;    Compare significands of numbers in internal format.
  395. ;    Guard words are included in the comparison.
  396. ;
  397. ;    unsigned short a[NI], b[NI];
  398. ;    cmpm( a, b );
  399. ;
  400. ;    for the significands:
  401. ;    returns    +1 if a > b
  402. ;         0 if a == b
  403. ;        -1 if a < b
  404. */
  405. int ecmpm( a, b )
  406. register unsigned short *a, *b;
  407. {
  408. int i;
  409.  
  410. a += M; /* skip up to significand area */
  411. b += M;
  412. for( i=M; i<NI; i++ )
  413.     {
  414.     if( *a++ != *b++ )
  415.         goto difrnt;
  416.     }
  417. return(0);
  418.  
  419. difrnt:
  420. if( *(--a) > *(--b) )
  421.     return(1);
  422. else
  423.     return(-1);
  424. }
  425.  
  426.  
  427. /*
  428. ;    Shift significand down by 1 bit
  429. */
  430.  
  431. void eshdn1(x)
  432. register unsigned short *x;
  433. {
  434. register unsigned short bits;
  435. int i;
  436.  
  437. x += M;    /* point to significand area */
  438.  
  439. bits = 0;
  440. for( i=M; i<NI; i++ )
  441.     {
  442.     if( *x & 1 )
  443.         bits |= 1;
  444.     *x >>= 1;
  445.     if( bits & 2 )
  446.         *x |= 0x8000;
  447.     bits <<= 1;
  448.     ++x;
  449.     }    
  450. }
  451.  
  452.  
  453.  
  454. /*
  455. ;    Shift significand up by 1 bit
  456. */
  457.  
  458. void eshup1(x)
  459. register unsigned short *x;
  460. {
  461. register unsigned short bits;
  462. int i;
  463.  
  464. x += NI-1;
  465. bits = 0;
  466.  
  467. for( i=M; i<NI; i++ )
  468.     {
  469.     if( *x & 0x8000 )
  470.         bits |= 1;
  471.     *x <<= 1;
  472.     if( bits & 2 )
  473.         *x |= 1;
  474.     bits <<= 1;
  475.     --x;
  476.     }
  477. }
  478.  
  479.  
  480.  
  481. /*
  482. ;    Shift significand down by 8 bits
  483. */
  484.  
  485. void eshdn8(x)
  486. register unsigned short *x;
  487. {
  488. register unsigned short newbyt, oldbyt;
  489. int i;
  490.  
  491. x += M;
  492. oldbyt = 0;
  493. for( i=M; i<NI; i++ )
  494.     {
  495.     newbyt = *x << 8;
  496.     *x >>= 8;
  497.     *x |= oldbyt;
  498.     oldbyt = newbyt;
  499.     ++x;
  500.     }
  501. }
  502.  
  503. /*
  504. ;    Shift significand up by 8 bits
  505. */
  506.  
  507. void eshup8(x)
  508. register unsigned short *x;
  509. {
  510. int i;
  511. register unsigned short newbyt, oldbyt;
  512.  
  513. x += NI-1;
  514. oldbyt = 0;
  515.  
  516. for( i=M; i<NI; i++ )
  517.     {
  518.     newbyt = *x >> 8;
  519.     *x <<= 8;
  520.     *x |= oldbyt;
  521.     oldbyt = newbyt;
  522.     --x;
  523.     }
  524. }
  525.  
  526. /*
  527. ;    Shift significand up by 16 bits
  528. */
  529.  
  530. void eshup6(x)
  531. register unsigned short *x;
  532. {
  533. int i;
  534. register unsigned short *p;
  535.  
  536. p = x + M;
  537. x += M + 1;
  538.  
  539. for( i=M; i<NI-1; i++ )
  540.     *p++ = *x++;
  541.  
  542. *p = 0;
  543. }
  544.  
  545. /*
  546. ;    Shift significand down by 16 bits
  547. */
  548.  
  549. void eshdn6(x)
  550. register unsigned short *x;
  551. {
  552. int i;
  553. register unsigned short *p;
  554.  
  555. x += NI-1;
  556. p = x + 1;
  557.  
  558. for( i=M; i<NI-1; i++ )
  559.     *(--p) = *(--x);
  560.  
  561. *(--p) = 0;
  562. }
  563.  
  564. /*
  565. ;    Add significands
  566. ;    x + y replaces y
  567. */
  568.  
  569. void eaddm( x, y )
  570. unsigned short *x, *y;
  571. {
  572. register unsigned long a;
  573. int i;
  574. unsigned int carry;
  575.  
  576. x += NI-1;
  577. y += NI-1;
  578. carry = 0;
  579. for( i=M; i<NI; i++ )
  580.     {
  581.     a = (unsigned long )(*x) + (unsigned long )(*y) + carry;
  582.     if( a & 0x10000 )
  583.         carry = 1;
  584.     else
  585.         carry = 0;
  586.     *y = (unsigned short )a;
  587.     --x;
  588.     --y;
  589.     }
  590. }
  591.  
  592. /*
  593. ;    Subtract significands
  594. ;    y - x replaces y
  595. */
  596.  
  597. void esubm( x, y )
  598. unsigned short *x, *y;
  599. {
  600. unsigned long a;
  601. int i;
  602. unsigned int carry;
  603.  
  604. x += NI-1;
  605. y += NI-1;
  606. carry = 0;
  607. for( i=M; i<NI; i++ )
  608.     {
  609.     a = (unsigned long )(*y) - (unsigned long )(*x) - carry;
  610.     if( a & 0x10000 )
  611.         carry = 1;
  612.     else
  613.         carry = 0;
  614.     *y = (unsigned short )a;
  615.     --x;
  616.     --y;
  617.     }
  618. }
  619.  
  620.  
  621. /* Divide significands */
  622.  
  623. static unsigned short equot[NI] = {0};
  624.  
  625. int edivm( den, num )
  626. unsigned short den[], num[];
  627. {
  628. int i;
  629. register unsigned short *p, *q;
  630. unsigned short j;
  631.  
  632. p = &equot[0];
  633. *p++ = num[0];
  634. *p++ = num[1];
  635.  
  636. for( i=M; i<NI; i++ )
  637.     {
  638.     *p++ = 0;
  639.     }
  640.  
  641. /* Use faster compare and subtraction if denominator
  642.  * has only 15 bits of significance.
  643.  */
  644. p = &den[M+2];
  645. if( *p++ == 0 )
  646.     {
  647.     for( i=M+3; i<NI; i++ )
  648.         {
  649.         if( *p++ != 0 )
  650.             goto fulldiv;
  651.         }
  652.     if( (den[M+1] & 1) != 0 )
  653.         goto fulldiv;
  654.     eshdn1(num);
  655.     eshdn1(den);
  656.  
  657.     p = &den[M+1];
  658.     q = &num[M+1];
  659.  
  660.     for( i=0; i<NBITS+2; i++ )
  661.         {
  662.         if( *p <= *q )
  663.             {
  664.             *q -= *p;
  665.             j = 1;
  666.             }
  667.         else
  668.             {
  669.             j = 0;
  670.             }
  671.         eshup1(equot);
  672.         equot[NI-2] |= j;
  673.         eshup1(num);
  674.         }
  675.     goto divdon;
  676.     }
  677.  
  678. /* The number of quotient bits to calculate is
  679.  * NBITS + 1 scaling guard bit + 1 roundoff bit.
  680.  */
  681. fulldiv:
  682.  
  683. p = &equot[NI-2];
  684. for( i=0; i<NBITS+2; i++ )
  685.     {
  686.     if( ecmpm(den,num) <= 0 )
  687.         {
  688.         esubm(den, num);
  689.         j = 1;    /* quotient bit = 1 */
  690.         }
  691.     else
  692.         j = 0;
  693.     eshup1(equot);
  694.     *p |= j;
  695.     eshup1(num);
  696.     }
  697.  
  698. divdon:
  699.  
  700. eshdn1( equot );
  701. eshdn1( equot );
  702.  
  703. /* test for nonzero remainder after roundoff bit */
  704. p = &num[M];
  705. j = 0;
  706. for( i=M; i<NI; i++ )
  707.     {
  708.     j |= *p++;
  709.     }
  710. if( j )
  711.     j = 1;
  712.  
  713.  
  714. for( i=0; i<NI; i++ )
  715.     num[i] = equot[i];
  716. return( (int )j );
  717. }
  718.  
  719.  
  720. /* Multiply significands */
  721. int emulm( a, b )
  722. unsigned short a[], b[];
  723. {
  724. unsigned short *p, *q;
  725. int i, j, k;
  726.  
  727. equot[0] = b[0];
  728. equot[1] = b[1];
  729. for( i=M; i<NI; i++ )
  730.     equot[i] = 0;
  731.  
  732. p = &a[NI-2];
  733. k = NBITS;
  734. while( *p == 0 ) /* significand is not supposed to be all zero */
  735.     {
  736.     eshdn6(a);
  737.     k -= 16;
  738.     }
  739. if( (*p & 0xff) == 0 )
  740.     {
  741.     eshdn8(a);
  742.     k -= 8;
  743.     }
  744.  
  745. q = &equot[NI-1];
  746. j = 0;
  747. for( i=0; i<k; i++ )
  748.     {
  749.     if( *p & 1 )
  750.         eaddm(b, equot);
  751. /* remember if there were any nonzero bits shifted out */
  752.     if( *q & 1 )
  753.         j |= 1;
  754.     eshdn1(a);
  755.     eshdn1(equot);
  756.     }
  757.  
  758. for( i=0; i<NI; i++ )
  759.     b[i] = equot[i];
  760.  
  761. /* return flag for lost nonzero bits */
  762. return(j);
  763. }
  764.  
  765.  
  766.  
  767. /*
  768.  * Normalize and round off.
  769.  *
  770.  * The internal format number to be rounded is "s".
  771.  * Input "lost" indicates whether or not the number is exact.
  772.  * This is the so-called sticky bit.
  773.  *
  774.  * Input "subflg" indicates whether the number was obtained
  775.  * by a subtraction operation.  In that case if lost is nonzero
  776.  * then the number is slightly smaller than indicated.
  777.  *
  778.  * Input "exp" is the biased exponent, which may be negative.
  779.  * the exponent field of "s" is ignored but is replaced by
  780.  * "exp" as adjusted by normalization and rounding.
  781.  *
  782.  * Input "rcntrl" is the rounding control.
  783.  */
  784.  
  785. static int rlast = -1;
  786. static int rw = 0;
  787. static unsigned short rmsk = 0;
  788. static unsigned short rmbit = 0;
  789. static unsigned short rebit = 0;
  790. static int re = 0;
  791. static unsigned short rbit[NI] = {0,0,0,0,0,0,0,0};
  792.  
  793. void emdnorm( s, lost, subflg, exp, rcntrl )
  794. unsigned short s[];
  795. int lost;
  796. int subflg;
  797. long exp;
  798. int rcntrl;
  799. {
  800. int i, j;
  801. unsigned short r;
  802.  
  803. /* Normalize */
  804. j = enormlz( s );
  805.  
  806. /* a blank significand could mean either zero or infinity. */
  807. #ifndef INFINITY
  808. if( j > NBITS )
  809.     {
  810.     ecleazs( s );
  811.     return;
  812.     }
  813. #endif
  814. exp -= j;
  815. #ifndef INFINITY
  816. if( exp >= 32767L )
  817.     goto overf;
  818. #endif
  819. if( exp < 0L )
  820.     {
  821.     if( exp > (long )(-NBITS-1) )
  822.         {
  823.         j = (int )exp;
  824.         i = eshift( s, j );
  825.         if( i )
  826.             lost = 1;
  827.         }
  828.     else
  829.         {
  830.         ecleazs( s );
  831.         return;
  832.         }
  833.     }
  834. /* Round off, unless told not to by rcntrl. */
  835. if( rcntrl == 0 )
  836.     goto mdfin;
  837. /* Set up rounding parameters if the control register changed. */
  838. if( rndprc != rlast )
  839.     {
  840.     ecleaz( rbit );
  841.     switch( rndprc )
  842.         {
  843.         default:
  844.         case NBITS:
  845.             rw = NI-1; /* low guard word */
  846.             rmsk = 0xffff;
  847.             rmbit = 0x8000;
  848.             rbit[rw-1] = 1;
  849.             re = NI-2;
  850.             rebit = 1;
  851.             break;
  852.         case 64:
  853.             rw = 7;
  854.             rmsk = 0xffff;
  855.             rmbit = 0x8000;
  856.             rbit[rw-1] = 1;
  857.             re = rw-1;
  858.             rebit = 1;
  859.             break;
  860. /* For DEC arithmetic */
  861.         case 56:
  862.             rw = 6;
  863.             rmsk = 0xff;
  864.             rmbit = 0x80;
  865.             rbit[rw] = 0x100;
  866.             re = rw;
  867.             rebit = 0x100;
  868.             break;
  869.         case 53:
  870.             rw = 6;
  871.             rmsk = 0x7ff;
  872.             rmbit = 0x0400;
  873.             rbit[rw] = 0x800;
  874.             re = rw;
  875.             rebit = 0x800;
  876.             break;
  877.         case 24:
  878.             rw = 4;
  879.             rmsk = 0xff;
  880.             rmbit = 0x80;
  881.             rbit[rw] = 0x100;
  882.             re = rw;
  883.             rebit = 0x100;
  884.             break;
  885.         }
  886.     rlast = rndprc;
  887.     }
  888.  
  889. if( rndprc >= 64 )
  890.     {
  891.     r = s[rw] & rmsk;
  892.     if( rndprc == 64 )
  893.         {
  894.         i = rw+1;
  895.         while( i < NI )
  896.             {
  897.             if( s[i] )
  898.                 r |= 1;
  899.             s[i] = 0;
  900.             ++i;
  901.             }
  902.         }
  903.     }
  904. else
  905.     {
  906.     if( exp <= 0 )
  907.         eshdn1(s);
  908.     r = s[rw] & rmsk;
  909. /* These tests assume NI = 8 */
  910.     i = rw+1;
  911.     while( i < NI )
  912.         {
  913.         if( s[i] )
  914.             r |= 1;
  915.         s[i] = 0;
  916.         ++i;
  917.         }
  918. /*
  919.     if( rndprc == 24 )
  920.         {
  921.         if( s[5] || s[6] )
  922.             r |= 1;
  923.         s[5] = 0;
  924.         s[6] = 0;
  925.         }
  926. */
  927.     s[rw] &= ~rmsk;
  928.     }
  929. if( (r & rmbit) != 0 )
  930.     {
  931.     if( r == rmbit)
  932.         {
  933.         if( lost == 0 )
  934.             { /* round to even */
  935.             if( (s[re] & rebit) == 0 )
  936.                 goto mddone;
  937.             }
  938.         else
  939.             {
  940.             if( subflg != 0 )
  941.                 goto mddone;
  942.             }
  943.         }
  944.     eaddm( rbit, s );
  945.     }
  946. mddone:
  947. if( (rndprc < 64) && (exp <= 0) )
  948.     {
  949.     eshup1(s);
  950.     }
  951. if( s[2] != 0 )
  952.     { /* overflow on roundoff */
  953.     eshdn1(s);
  954.     exp += 1;
  955.     }
  956. mdfin:
  957. s[NI-1] = 0;
  958. if( exp >= 32767L )
  959.     {
  960. #ifndef INFINITY
  961. overf:
  962. #endif
  963. #ifdef INFINITY
  964.     s[1] = 32767;
  965.     for( i=2; i<NI-1; i++ )
  966.         s[i] = 0;
  967. #else
  968.     s[1] = 32766;
  969.     s[2] = 0;
  970.     for( i=M+1; i<NI-1; i++ )
  971.         s[i] = 0xffff;
  972.     s[NI-1] = 0;
  973.     if( rndprc < 64 )
  974.         {
  975.         s[rw] &= ~rmsk;
  976.         if( rndprc == 24 )
  977.             {
  978.             s[5] = 0;
  979.             s[6] = 0;
  980.             }
  981.         }
  982. #endif
  983.     return;
  984.     }
  985. if( exp < 0 )
  986.     s[1] = 0;
  987. else
  988.     s[1] = (unsigned short )exp;
  989. }
  990.  
  991.  
  992.  
  993. /*
  994. ;    Subtract external format numbers.
  995. ;
  996. ;    unsigned short a[NE], b[NE], c[NE];
  997. ;    esub( a, b, c );     c = b - a
  998. */
  999.  
  1000. static int subflg = 0;
  1001.  
  1002. void esub( a, b, c )
  1003. unsigned short *a, *b, *c;
  1004. {
  1005.  
  1006. subflg = 1;
  1007. eadd1( a, b, c );
  1008. }
  1009.  
  1010.  
  1011. /*
  1012. ;    Add.
  1013. ;
  1014. ;    unsigned short a[NE], b[NE], c[NE];
  1015. ;    eadd( a, b, c );     c = b + a
  1016. */
  1017. void eadd( a, b, c )
  1018. unsigned short *a, *b, *c;
  1019. {
  1020.  
  1021. subflg = 0;
  1022. eadd1( a, b, c );
  1023. }
  1024.  
  1025. void eadd1( a, b, c )
  1026. unsigned short *a, *b, *c;
  1027. {
  1028. unsigned short ai[NI], bi[NI], ci[NI];
  1029. int i, lost, j, k;
  1030. long lt, lta, ltb;
  1031.  
  1032. #ifdef INFINITY
  1033. if( eisinf(a) )
  1034.     {
  1035.     emov(a,c);
  1036.     if( subflg )
  1037.         eneg(c);
  1038.     return;
  1039.     }
  1040. if( eisinf(b) )
  1041.     {
  1042.     emov(b,c);
  1043.     return;
  1044.     }
  1045. #endif
  1046. emovi( a, ai );
  1047. emovi( b, bi );
  1048. if( subflg )
  1049.     ai[0] = ~ai[0];
  1050.  
  1051. /* compare exponents */
  1052. lta = ai[E];
  1053. ltb = bi[E];
  1054. lt = lta - ltb;
  1055. if( lt > 0L )
  1056.     {    /* put the larger number in bi */
  1057.     emovz( bi, ci );
  1058.     emovz( ai, bi );
  1059.     emovz( ci, ai );
  1060.     ltb = bi[E];
  1061.     lt = -lt;
  1062.     }
  1063. lost = 0;
  1064. if( lt != 0L )
  1065.     {
  1066.     if( lt < (long )(-NBITS-1) )
  1067.         goto done;    /* answer same as larger addend */
  1068.     k = (int )lt;
  1069.     lost = eshift( ai, k ); /* shift the smaller number down */
  1070.     }
  1071. else
  1072.     {
  1073. /* exponents were the same, so must compare significands */
  1074.     i = ecmpm( ai, bi );
  1075.     if( i == 0 )
  1076.         { /* the numbers are identical in magnitude */
  1077.         /* if different signs, result is zero */
  1078.         if( ai[0] != bi[0] )
  1079.             {
  1080.             eclear(c);
  1081.             return;
  1082.             }
  1083.         /* if same sign, result is double */
  1084.         /* double denomalized tiny number */
  1085.         if( (bi[E] == 0) && ((bi[3] & 0x8000) == 0) )
  1086.             {
  1087.             eshup1( bi );
  1088.             goto done;
  1089.             }
  1090.         /* add 1 to exponent unless both are zero! */
  1091.         for( j=1; j<NI-1; j++ )
  1092.             {
  1093.             if( bi[j] != 0 )
  1094.                 {
  1095. /* This could overflow, but let emovo take care of that. */
  1096.                 ltb += 1;
  1097.                 break;
  1098.                 }
  1099.             }
  1100.         bi[E] = (unsigned short )ltb;
  1101.         goto done;
  1102.         }
  1103.     if( i > 0 )
  1104.         {    /* put the larger number in bi */
  1105.         emovz( bi, ci );
  1106.         emovz( ai, bi );
  1107.         emovz( ci, ai );
  1108.         }
  1109.     }
  1110. if( ai[0] == bi[0] )
  1111.     {
  1112.     eaddm( ai, bi );
  1113.     subflg = 0;
  1114.     }
  1115. else
  1116.     {
  1117.     esubm( ai, bi );
  1118.     subflg = 1;
  1119.     }
  1120. emdnorm( bi, lost, subflg, ltb, 64 );
  1121.  
  1122. done:
  1123. emovo( bi, c );
  1124. }
  1125.  
  1126.  
  1127.  
  1128. /*
  1129. ;    Divide.
  1130. ;
  1131. ;    unsigned short a[NE], b[NE], c[NE];
  1132. ;    ediv( a, b, c );    c = b / a
  1133. */
  1134. void ediv( a, b, c )
  1135. unsigned short *a, *b, *c;
  1136. {
  1137. unsigned short ai[NI], bi[NI];
  1138. int i;
  1139. long lt, lta, ltb;
  1140.  
  1141. #ifdef INFINITY
  1142. if( eisinf(b) )
  1143.     {
  1144.     if( eisneg(a) ^ eisneg(b) )
  1145.         *(c+(NE-1)) = 0x8000;
  1146.     else
  1147.         *(c+(NE-1)) = 0;
  1148.     einfin(c);
  1149.     return;
  1150.     }
  1151. if( eisinf(a) )
  1152.     {
  1153.     eclear(c);
  1154.     return;
  1155.     }
  1156. #endif
  1157. emovi( a, ai );
  1158. emovi( b, bi );
  1159. lta = ai[E];
  1160. ltb = bi[E];
  1161. if( bi[E] == 0 )
  1162.     { /* See if numerator is zero. */
  1163.     for( i=1; i<NI-1; i++ )
  1164.         {
  1165.         if( bi[i] != 0 )
  1166.             {
  1167.             ltb -= enormlz( bi );
  1168.             goto dnzro1;
  1169.             }
  1170.         }
  1171.     eclear(c);
  1172.     return;
  1173.     }
  1174. dnzro1:
  1175.  
  1176. if( ai[E] == 0 )
  1177.     {    /* possible divide by zero */
  1178.     for( i=1; i<NI-1; i++ )
  1179.         {
  1180.         if( ai[i] != 0 )
  1181.             {
  1182.             lta -= enormlz( ai );
  1183.             goto dnzro2;
  1184.             }
  1185.         }
  1186.     if( ai[0] == bi[0] )
  1187.         *(c+(NE-1)) = 0;
  1188.     else
  1189.         *(c+(NE-1)) = 0x8000;
  1190.     einfin(c);
  1191.     mtherr( "ediv", SING );
  1192.     return;
  1193.     }
  1194. dnzro2:
  1195.  
  1196. i = edivm( ai, bi );
  1197. /* calculate exponent */
  1198. lt = ltb - lta + EXONE;
  1199. emdnorm( bi, i, 0, lt, 64 );
  1200. /* set the sign */
  1201. if( ai[0] == bi[0] )
  1202.     bi[0] = 0;
  1203. else
  1204.     bi[0] = 0Xffff;
  1205. emovo( bi, c );
  1206. }
  1207.  
  1208.  
  1209.  
  1210. /*
  1211. ;    Multiply.
  1212. ;
  1213. ;    unsigned short a[NE], b[NE], c[NE];
  1214. ;    emul( a, b, c );    c = b * a
  1215. */
  1216. void emul( a, b, c )
  1217. unsigned short *a, *b, *c;
  1218. {
  1219. unsigned short ai[NI], bi[NI];
  1220. int i, j;
  1221. long lt, lta, ltb;
  1222.  
  1223. #ifdef INFINITY
  1224. if( eisinf(a) || eisinf(b) )
  1225.     {
  1226.     if( eisneg(a) ^ eisneg(b) )
  1227.         *(c+(NE-1)) = 0x8000;
  1228.     else
  1229.         *(c+(NE-1)) = 0;
  1230.     einfin(c);
  1231.     return;
  1232.     }
  1233. #endif
  1234. emovi( a, ai );
  1235. emovi( b, bi );
  1236. lta = ai[E];
  1237. ltb = bi[E];
  1238. if( ai[E] == 0 )
  1239.     {
  1240.     for( i=1; i<NI-1; i++ )
  1241.         {
  1242.         if( ai[i] != 0 )
  1243.             {
  1244.             lta -= enormlz( ai );
  1245.             goto mnzer1;
  1246.             }
  1247.         }
  1248.     eclear(c);
  1249.     return;
  1250.     }
  1251. mnzer1:
  1252.  
  1253. if( bi[E] == 0 )
  1254.     {
  1255.     for( i=1; i<NI-1; i++ )
  1256.         {
  1257.         if( bi[i] != 0 )
  1258.             {
  1259.             ltb -= enormlz( bi );
  1260.             goto mnzer2;
  1261.             }
  1262.         }
  1263.     eclear(c);
  1264.     return;
  1265.     }
  1266. mnzer2:
  1267.  
  1268. /* Multiply significands */
  1269. j = emulm( ai, bi );
  1270. /* calculate exponent */
  1271. lt = lta + ltb - (EXONE - 1);
  1272. emdnorm( bi, j, 0, lt, 64 );
  1273. /* calculate sign of product */
  1274. if( ai[0] == bi[0] )
  1275.     bi[0] = 0;
  1276. else
  1277.     bi[0] = 0xffff;
  1278. emovo( bi, c );
  1279. }
  1280.  
  1281.  
  1282.  
  1283.  
  1284. /*
  1285. ; Convert IEEE double precision to e type
  1286. ;    double d;
  1287. ;    unsigned short x[N+2];
  1288. ;    e53toe( &d, x );
  1289. */
  1290. void e53toe( e, y )
  1291. unsigned short *e, *y;
  1292. {
  1293. #ifdef DEC
  1294.  
  1295. dectoe( e, y ); /* see etodec.c */
  1296.  
  1297. #else
  1298.  
  1299. register unsigned short r;
  1300. register unsigned short *p;
  1301. unsigned short yy[NI];
  1302. int denorm, k;
  1303.  
  1304. denorm = 0;    /* flag if denormalized number */
  1305. ecleaz(yy);
  1306. #ifdef IBMPC
  1307. e += 3;
  1308. #endif
  1309. r = *e;
  1310. yy[0] = 0;
  1311. if( r & 0x8000 )
  1312.     yy[0] = 0xffff;
  1313. yy[M] = (r & 0x0f) | 0x10;
  1314. r &= ~0x800f;    /* strip sign and 4 significand bits */
  1315. #ifdef INFINITY
  1316. if( r == 0x7ff0 )
  1317.     {
  1318.     einfin( y );
  1319.     if( r & 0x8000 )
  1320.         eneg(y);
  1321.     return;
  1322.     }
  1323. #endif
  1324. r >>= 4;
  1325. /* If zero exponent, then the significand is denormalized.
  1326.  * So, take back the understood high significand bit. */ 
  1327. if( r == 0 )
  1328.     {
  1329.     denorm = 1;
  1330.     yy[M] &= ~0x10;
  1331.     }
  1332. r += EXONE - 01777;
  1333. yy[E] = r;
  1334. p = &yy[M+1];
  1335. #ifdef IBMPC
  1336. *p++ = *(--e);
  1337. *p++ = *(--e);
  1338. *p++ = *(--e);
  1339. #endif
  1340. #ifdef MIEEE
  1341. ++e;
  1342. *p++ = *e++;
  1343. *p++ = *e++;
  1344. *p++ = *e++;
  1345. #endif
  1346. (void )eshift( yy, -5 );
  1347. if( denorm )
  1348.     { /* if zero exponent, then normalize the significand */
  1349.     if( (k = enormlz(yy)) > NBITS )
  1350.         ecleazs(yy);
  1351.     else
  1352.         yy[E] -= (unsigned short )(k-1);
  1353.     }
  1354. emovo( yy, y );
  1355. #endif /* not DEC */
  1356. }
  1357.  
  1358. void e64toe( e, y )
  1359. unsigned short *e, *y;
  1360. {
  1361. unsigned short yy[NI];
  1362. unsigned short *p, *q;
  1363. int i;
  1364.  
  1365. p = yy;
  1366. for( i=0; i<NE-5; i++ )
  1367.     *p++ = 0;
  1368. #ifdef IBMPC
  1369. for( i=0; i<5; i++ )
  1370.     *p++ = *e++;
  1371. #endif
  1372. #ifdef MIEEE
  1373. p = &yy[0] + (NE-1);
  1374. *p-- = *e++;
  1375. ++e;
  1376. for( i=0; i<4; i++ )
  1377.     *p-- = *e++;
  1378. #endif
  1379. p = yy;
  1380. q = y;
  1381. #ifdef INFINITY
  1382. if( *p == 0x7fff )
  1383.     {
  1384.     einfin( y );
  1385.     if( *p & 0x8000 )
  1386.         eneg(y);
  1387.     return;
  1388.     }
  1389. #endif
  1390. for( i=0; i<NE; i++ )
  1391.     *q++ = *p++;
  1392. }
  1393.  
  1394.  
  1395. /*
  1396. ; Convert IEEE single precision to e type
  1397. ;    float d;
  1398. ;    unsigned short x[N+2];
  1399. ;    dtox( &d, x );
  1400. */
  1401. void e24toe( e, y )
  1402. unsigned short *e, *y;
  1403. {
  1404. register unsigned short r;
  1405. register unsigned short *p;
  1406. unsigned short yy[NI];
  1407. int denorm, k;
  1408.  
  1409. denorm = 0;    /* flag if denormalized number */
  1410. ecleaz(yy);
  1411. #ifdef IBMPC
  1412. e += 1;
  1413. #endif
  1414. #ifdef DEC
  1415. e += 1;
  1416. #endif
  1417. r = *e;
  1418. yy[0] = 0;
  1419. if( r & 0x8000 )
  1420.     yy[0] = 0xffff;
  1421. yy[M] = (r & 0x7f) | 0200;
  1422. r &= ~0x807f;    /* strip sign and 7 significand bits */
  1423. #ifdef INFINITY
  1424. if( r == 0x7f80 )
  1425.     {
  1426.     einfin( y );
  1427.     if( r & 0x8000 )
  1428.         eneg(y);
  1429.     return;
  1430.     }
  1431. #endif
  1432. r >>= 7;
  1433. /* If zero exponent, then the significand is denormalized.
  1434.  * So, take back the understood high significand bit. */ 
  1435. if( r == 0 )
  1436.     {
  1437.     denorm = 1;
  1438.     yy[M] &= ~0200;
  1439.     }
  1440. r += EXONE - 0177;
  1441. yy[E] = r;
  1442. p = &yy[M+1];
  1443. #ifdef IBMPC
  1444. *p++ = *(--e);
  1445. #endif
  1446. #ifdef DEC
  1447. *p++ = *(--e);
  1448. #endif
  1449. #ifdef MIEEE
  1450. ++e;
  1451. *p++ = *e++;
  1452. #endif
  1453. (void )eshift( yy, -8 );
  1454. if( denorm )
  1455.     { /* if zero exponent, then normalize the significand */
  1456.     if( (k = enormlz(yy)) > NBITS )
  1457.         ecleazs(yy);
  1458.     else
  1459.         yy[E] -= (unsigned short )(k-1);
  1460.     }
  1461. emovo( yy, y );
  1462. }
  1463.  
  1464.  
  1465. void etoe64( x, e )
  1466. unsigned short *x, *e;
  1467. {
  1468. unsigned short xi[NI];
  1469. long exp;
  1470. int rndsav;
  1471.  
  1472. emovi( x, xi );
  1473. exp = (long )xi[E]; /* adjust exponent for offset */
  1474. #ifdef INFINITY
  1475. if( eisinf(x) )
  1476.     goto nonorm;
  1477. #endif
  1478. /* round off to nearest or even */
  1479. rndsav = rndprc;
  1480. rndprc = 64;
  1481. emdnorm( xi, 0, 0, exp, 64 );
  1482. rndprc = rndsav;
  1483. nonorm:
  1484. toe64( xi, e );
  1485. }
  1486.  
  1487. /* move out internal format to ieee long double */
  1488. static void toe64( a, b )
  1489. unsigned short *a, *b;
  1490. {
  1491. register unsigned short *p, *q;
  1492. unsigned short i;
  1493.  
  1494. p = a;
  1495. #ifdef MIEEE
  1496. q = b;
  1497. #else
  1498. q = b + 4; /* point to output exponent */
  1499. #endif
  1500.  
  1501. /* combine sign and exponent */
  1502. i = *p++;
  1503. #ifdef MIEEE
  1504. if( i )
  1505.     *q++ = *p++ | 0x8000;
  1506. else
  1507.     *q++ = *p++;
  1508. *q++ = 0;
  1509. #else
  1510. if( i )
  1511.     *q-- = *p++ | 0x8000;
  1512. else
  1513.     *q-- = *p++;
  1514. #endif
  1515. /* skip over guard word */
  1516. ++p;
  1517. /* move the significand */
  1518. #ifdef MIEEE
  1519. for( i=0; i<4; i++ )
  1520.     *q++ = *p++;
  1521. #else
  1522. for( i=0; i<4; i++ )
  1523.     *q-- = *p++;
  1524. #endif
  1525. }
  1526.  
  1527.  
  1528. /*
  1529. ; e type to IEEE double precision
  1530. ;    double d;
  1531. ;    unsigned short x[NE];
  1532. ;    etoe53( x, &d );
  1533. */
  1534.  
  1535. #ifdef DEC
  1536.  
  1537. void etoe53( x, e )
  1538. unsigned short *x, *e;
  1539. {
  1540. etodec( x, e ); /* see etodec.c */
  1541. }
  1542.  
  1543. static void toe53( x, y )
  1544. unsigned short *x, *y;
  1545. {
  1546. todec( x, y );
  1547. }
  1548.  
  1549. #else
  1550.  
  1551. void etoe53( x, e )
  1552. unsigned short *x, *e;
  1553. {
  1554. unsigned short xi[NI];
  1555. long exp;
  1556. int rndsav;
  1557.  
  1558. emovi( x, xi );
  1559. exp = (long )xi[E] - (EXONE - 0x3ff); /* adjust exponent for offsets */
  1560. #ifdef INFINITY
  1561. if( eisinf(x) )
  1562.     goto nonorm;
  1563. #endif
  1564. /* round off to nearest or even */
  1565. rndsav = rndprc;
  1566. rndprc = 53;
  1567. emdnorm( xi, 0, 0, exp, 64 );
  1568. rndprc = rndsav;
  1569. nonorm:
  1570. toe53( xi, e );
  1571. }
  1572.  
  1573.  
  1574. static void toe53( x, y )
  1575. unsigned short *x, *y;
  1576. {
  1577. unsigned short i;
  1578. unsigned short *p;
  1579.  
  1580.  
  1581. p = &x[0];
  1582. #ifdef IBMPC
  1583. y += 3;
  1584. #endif
  1585. *y = 0;    /* output high order */
  1586. if( *p++ )
  1587.     *y = 0x8000;    /* output sign bit */
  1588.  
  1589. i = *p++;
  1590. if( i >= (unsigned int )2047 )
  1591.     {    /* Saturate at largest number less than infinity. */
  1592. #ifdef INFINITY
  1593.     *y |= 0x7ff0;
  1594. #ifdef IBMPC
  1595.     *(--y) = 0;
  1596.     *(--y) = 0;
  1597.     *(--y) = 0;
  1598. #endif
  1599. #ifdef MIEEE
  1600.     ++y;
  1601.     *y++ = 0;
  1602.     *y++ = 0;
  1603.     *y++ = 0;
  1604. #endif
  1605. #else
  1606.     *y |= (unsigned short )0x7fef;
  1607. #ifdef IBMPC
  1608.     *(--y) = 0xffff;
  1609.     *(--y) = 0xffff;
  1610.     *(--y) = 0xffff;
  1611. #endif
  1612. #ifdef MIEEE
  1613.     ++y;
  1614.     *y++ = 0xffff;
  1615.     *y++ = 0xffff;
  1616.     *y++ = 0xffff;
  1617. #endif
  1618. #endif
  1619.     return;
  1620.     }
  1621. if( i == 0 )
  1622.     {
  1623.     (void )eshift( x, 4 );
  1624.     }
  1625. else
  1626.     {
  1627.     i <<= 4;
  1628.     (void )eshift( x, 5 );
  1629.     }
  1630. i |= *p++ & (unsigned short )0x0f;    /* *p = xi[M] */
  1631. *y |= (unsigned short )i; /* high order output already has sign bit set */
  1632. #ifdef IBMPC
  1633. *(--y) = *p++;
  1634. *(--y) = *p++;
  1635. *(--y) = *p;
  1636. #endif
  1637. #ifdef MIEEE
  1638. ++y;
  1639. *y++ = *p++;
  1640. *y++ = *p++;
  1641. *y++ = *p++;
  1642. #endif
  1643. }
  1644.  
  1645. #endif /* not DEC */
  1646.  
  1647.  
  1648.  
  1649. /*
  1650. ; e type to IEEE single precision
  1651. ;    float d;
  1652. ;    unsigned short x[N+2];
  1653. ;    xtod( x, &d );
  1654. */
  1655. void etoe24( x, e )
  1656. unsigned short *x, *e;
  1657. {
  1658. long exp;
  1659. unsigned short xi[NI];
  1660. int rndsav;
  1661.  
  1662. emovi( x, xi );
  1663. exp = (long )xi[E] - (EXONE - 0177); /* adjust exponent for offsets */
  1664. #ifdef INFINITY
  1665. if( eisinf(x) )
  1666.     goto nonorm;
  1667. #endif
  1668. /* round off to nearest or even */
  1669. rndsav = rndprc;
  1670. rndprc = 24;
  1671. emdnorm( xi, 0, 0, exp, 64 );
  1672. rndprc = rndsav;
  1673. nonorm:
  1674. toe24( xi, e );
  1675. }
  1676.  
  1677. static void toe24( x, y )
  1678. unsigned short *x, *y;
  1679. {
  1680. unsigned short i;
  1681. unsigned short *p;
  1682.  
  1683. p = &x[0];
  1684. #ifdef IBMPC
  1685. y += 1;
  1686. #endif
  1687. #ifdef DEC
  1688. y += 1;
  1689. #endif
  1690. *y = 0;    /* output high order */
  1691. if( *p++ )
  1692.     *y = 0x8000;    /* output sign bit */
  1693.  
  1694. i = *p++;
  1695. if( i >= 255 )
  1696.     {    /* Saturate at largest number less than infinity. */
  1697. #ifdef INFINITY
  1698.     *y |= (unsigned short )0x7f80;
  1699. #ifdef IBMPC
  1700.     *(--y) = 0;
  1701. #endif
  1702. #ifdef DEC
  1703.     *(--y) = 0;
  1704. #endif
  1705. #ifdef MIEEE
  1706.     ++y;
  1707.     *y = 0;
  1708. #endif
  1709. #else
  1710.     *y |= (unsigned short )0x7f7f;
  1711. #ifdef IBMPC
  1712.     *(--y) = 0xffff;
  1713. #endif
  1714. #ifdef DEC
  1715.     *(--y) = 0xffff;
  1716. #endif
  1717. #ifdef MIEEE
  1718.     ++y;
  1719.     *y = 0xffff;
  1720. #endif
  1721. #endif
  1722.     return;
  1723.     }
  1724. if( i == 0 )
  1725.     {
  1726.     (void )eshift( x, 7 );
  1727.     }
  1728. else
  1729.     {
  1730.     i <<= 7;
  1731.     (void )eshift( x, 8 );
  1732.     }
  1733. i |= *p++ & (unsigned short )0x7f;    /* *p = xi[M] */
  1734. *y |= i;    /* high order output already has sign bit set */
  1735. #ifdef IBMPC
  1736. *(--y) = *p;
  1737. #endif
  1738. #ifdef DEC
  1739. *(--y) = *p;
  1740. #endif
  1741. #ifdef MIEEE
  1742. ++y;
  1743. *y = *p;
  1744. #endif
  1745. }
  1746.  
  1747.  
  1748. /* Compare two e type numbers.
  1749.  *
  1750.  * unsigned short a[NE], b[NE];
  1751.  * ecmp( a, b );
  1752.  *
  1753.  *  returns +1 if a > b
  1754.  *           0 if a == b
  1755.  *          -1 if a < b
  1756.  */
  1757. int ecmp( a, b )
  1758. unsigned short *a, *b;
  1759. {
  1760. unsigned short ai[NI], bi[NI];
  1761. register unsigned short *p, *q;
  1762. register int i;
  1763. int msign;
  1764.  
  1765. emovi( a, ai );
  1766. p = ai;
  1767. emovi( b, bi );
  1768. q = bi;
  1769.  
  1770. if( *p != *q )
  1771.     { /* the signs are different */
  1772. /* -0 equals + 0 */
  1773.     for( i=1; i<NI-1; i++ )
  1774.         {
  1775.         if( ai[i] != 0 )
  1776.             goto nzro;
  1777.         if( bi[i] != 0 )
  1778.             goto nzro;
  1779.         }
  1780.     return(0);
  1781. nzro:
  1782.     if( *p == 0 )
  1783.         return( 1 );
  1784.     else
  1785.         return( -1 );
  1786.     }
  1787. /* both are the same sign */
  1788. if( *p == 0 )
  1789.     msign = 1;
  1790. else
  1791.     msign = -1;
  1792. i = NI-1;
  1793. do
  1794.     {
  1795.     if( *p++ != *q++ )
  1796.         {
  1797.         goto diff;
  1798.         }
  1799.     }
  1800. while( --i > 0 );
  1801.  
  1802. return(0);    /* equality */
  1803.  
  1804.  
  1805.  
  1806. diff:
  1807.  
  1808. if( *(--p) > *(--q) )
  1809.     return( msign );        /* p is bigger */
  1810. else
  1811.     return( -msign );    /* p is littler */
  1812. }
  1813.  
  1814.  
  1815.  
  1816.  
  1817. /* Find nearest integer to x = floor( x + 0.5 )
  1818.  *
  1819.  * unsigned short x[NE], y[NE]
  1820.  * eround( x, y );
  1821.  */
  1822. void eround( x, y )
  1823. unsigned short *x, *y;
  1824. {
  1825.  
  1826. eadd( ehalf, x, y );
  1827. efloor( y, y );
  1828. }
  1829.  
  1830.  
  1831.  
  1832.  
  1833. /*
  1834. ; convert long (32-bit) integer to e type
  1835. ;
  1836. ;    long l;
  1837. ;    unsigned short x[NE];
  1838. ;    ltoe( &l, x );
  1839. ; note &l is the memory address of l
  1840. */
  1841. void ltoe( lp, y )
  1842. long *lp;    /* lp is the memory address of a long integer */
  1843. unsigned short *y;    /* y is the address of a short */
  1844. {
  1845. unsigned short yi[NI];
  1846. unsigned long ll;
  1847. int k;
  1848.  
  1849. ecleaz( yi );
  1850. if( *lp < 0 )
  1851.     {
  1852.     ll =  (unsigned long )( -(*lp) ); /* make it positive */
  1853.     yi[0] = 0xffff; /* put correct sign in the e type number */
  1854.     }
  1855. else
  1856.     {
  1857.     ll = (unsigned long )( *lp );
  1858.     }
  1859. /* move the long integer to yi significand area */
  1860. yi[M] = (unsigned short )(ll >> 16); 
  1861. yi[M+1] = (unsigned short )ll;
  1862.  
  1863. yi[E] = EXONE + 15;    /* exponent if normalize shift count were 0 */
  1864. if( (k = enormlz( yi )) > NBITS ) /* normalize the significand */
  1865.     ecleaz( yi );    /* it was zero */
  1866. else
  1867.     yi[E] -= (unsigned short )k; /* subtract shift count from exponent */
  1868. emovo( yi, y );    /* output the answer */
  1869. }
  1870.  
  1871. /*
  1872. ; convert unsigned long (32-bit) integer to e type
  1873. ;
  1874. ;    unsigned long l;
  1875. ;    unsigned short x[NE];
  1876. ;    ltox( &l, x );
  1877. ; note &l is the memory address of l
  1878. */
  1879. void ultoe( lp, y )
  1880. unsigned long *lp; /* lp is the memory address of a long integer */
  1881. unsigned short *y;    /* y is the address of a short */
  1882. {
  1883. unsigned short yi[NI];
  1884. unsigned long ll;
  1885. int k;
  1886.  
  1887. ecleaz( yi );
  1888. ll = *lp;
  1889.  
  1890. /* move the long integer to ayi significand area */
  1891. yi[M] = (unsigned short )(ll >> 16); 
  1892. yi[M+1] = (unsigned short )ll;
  1893.  
  1894. yi[E] = EXONE + 15;    /* exponent if normalize shift count were 0 */
  1895. if( (k = enormlz( yi )) > NBITS ) /* normalize the significand */
  1896.     ecleaz( yi );    /* it was zero */
  1897. else
  1898.     yi[E] -= (unsigned short )k; /* subtract shift count from exponent */
  1899. emovo( yi, y );    /* output the answer */
  1900. }
  1901.  
  1902.  
  1903. /*
  1904. ;    Find long integer and fractional parts
  1905.  
  1906. ;    long i;
  1907. ;    unsigned short x[NE], frac[NE];
  1908. ;    xifrac( x, &i, frac );
  1909.  
  1910.   The integer output has the sign of the input.  The fraction is
  1911.   the positive fractional part of abs(x).
  1912. */
  1913. void eifrac( x, i, frac )
  1914. unsigned short *x;
  1915. long *i;
  1916. unsigned short *frac;
  1917. {
  1918. unsigned short xi[NI];
  1919. int k;
  1920.  
  1921. emovi( x, xi );
  1922. k = (int )xi[E] - (EXONE - 1);
  1923. if( k <= 0 )
  1924.     {
  1925. /* if exponent <= 0, integer = 0 and real output is fraction */
  1926.     *i = 0L;
  1927.     emovo( xi, frac );
  1928.     return;
  1929.     }
  1930. if( k > 31 )
  1931.     {
  1932. /*
  1933. ;    long integer overflow: output large integer
  1934. ;    and correct fraction
  1935. */
  1936.     if( xi[0] )
  1937.         *i = 0x80000000;
  1938.     else
  1939.         *i = 0x7fffffff;
  1940.     (void )eshift( xi, k );
  1941.     goto lab11;
  1942.     }
  1943.  
  1944. if( k > 16 )
  1945.     {
  1946. /*
  1947. ; shift more than 16 bits: shift up k-16, output the integer,
  1948. ; then complete the shift to get the fraction.
  1949. */
  1950.     k -= 16;
  1951.     (void )eshift( xi, k );
  1952.  
  1953.     *i = (long )(((unsigned long )xi[M] << 16) | xi[M+1]);
  1954.     eshup6( xi );
  1955.     goto lab10;
  1956.     }
  1957.  
  1958. /* shift not more than 16 bits */
  1959. (void )eshift( xi, k );
  1960. *i = (long )xi[M] & 0xffff;
  1961.  
  1962. lab10:
  1963.  
  1964. if( xi[0] )
  1965.     *i = -(*i);
  1966. lab11:
  1967.  
  1968. xi[0] = 0;
  1969. xi[E] = EXONE - 1;
  1970. xi[M] = 0;
  1971. if( (k = enormlz( xi )) > NBITS )
  1972.     ecleaz( xi );
  1973. else
  1974.     xi[E] -= (unsigned short )k;
  1975.  
  1976. emovo( xi, frac );
  1977. }
  1978.  
  1979.  
  1980. /*
  1981. ;    Find unsigned long integer and fractional parts
  1982.  
  1983. ;    unsigned long i;
  1984. ;    unsigned short x[NE], frac[NE];
  1985. ;    xifrac( x, &i, frac );
  1986.  
  1987.   A negative e type input yields integer output = 0
  1988.   but correct fraction.
  1989. */
  1990. void euifrac( x, i, frac )
  1991. unsigned short *x;
  1992. long *i;
  1993. unsigned short *frac;
  1994. {
  1995. unsigned short xi[NI];
  1996. int k;
  1997.  
  1998. emovi( x, xi );
  1999. k = (int )xi[E] - (EXONE - 1);
  2000. if( k <= 0 )
  2001.     {
  2002. /* if exponent <= 0, integer = 0 and argument is fraction */
  2003.     *i = 0L;
  2004.     emovo( xi, frac );
  2005.     return;
  2006.     }
  2007. if( k > 32 )
  2008.     {
  2009. /*
  2010. ;    long integer overflow: output large integer
  2011. ;    and correct fraction
  2012. */
  2013.     *i = 0xffffffff;
  2014.     (void )eshift( xi, k );
  2015.     goto lab10;
  2016.     }
  2017.  
  2018. if( k > 16 )
  2019.     {
  2020. /*
  2021. ; shift more than 16 bits: shift up k-16, output the integer,
  2022. ; then complete the shift to get the fraction.
  2023. */
  2024.     k -= 16;
  2025.     (void )eshift( xi, k );
  2026.  
  2027.     *i = (long )(((unsigned long )xi[M] << 16) | xi[M+1]);
  2028.     eshup6( xi );
  2029.     goto lab10;
  2030.     }
  2031.  
  2032. /* shift not more than 16 bits */
  2033. (void )eshift( xi, k );
  2034. *i = (long )xi[M] & 0xffff;
  2035.  
  2036. lab10:
  2037.  
  2038. if( xi[0] )
  2039.     *i = 0L;
  2040.  
  2041. xi[0] = 0;
  2042. xi[E] = EXONE - 1;
  2043. xi[M] = 0;
  2044. if( (k = enormlz( xi )) > NBITS )
  2045.     ecleaz( xi );
  2046. else
  2047.     xi[E] -= (unsigned short )k;
  2048.  
  2049. emovo( xi, frac );
  2050. }
  2051.  
  2052.  
  2053.  
  2054. /*
  2055. ;    Shift significand
  2056. ;
  2057. ;    Shifts significand area up or down by the number of bits
  2058. ;    given by the variable sc.
  2059. */
  2060. int eshift( x, sc )
  2061. unsigned short *x;
  2062. int sc;
  2063. {
  2064. unsigned short lost;
  2065. unsigned short *p;
  2066.  
  2067. if( sc == 0 )
  2068.     return( 0 );
  2069.  
  2070. lost = 0;
  2071. p = x + NI-1;
  2072.  
  2073. if( sc < 0 )
  2074.     {
  2075.     sc = -sc;
  2076.     while( sc >= 16 )
  2077.         {
  2078.         lost |= *p;    /* remember lost bits */
  2079.         eshdn6(x);
  2080.         sc -= 16;
  2081.         }
  2082.  
  2083.     while( sc >= 8 )
  2084.         {
  2085.         lost |= *p & 0xff;
  2086.         eshdn8(x);
  2087.         sc -= 8;
  2088.         }
  2089.  
  2090.     while( sc > 0 )
  2091.         {
  2092.         lost |= *p & 1;
  2093.         eshdn1(x);
  2094.         sc -= 1;
  2095.         }
  2096.     }
  2097. else
  2098.     {
  2099.     while( sc >= 16 )
  2100.         {
  2101.         eshup6(x);
  2102.         sc -= 16;
  2103.         }
  2104.  
  2105.     while( sc >= 8 )
  2106.         {
  2107.         eshup8(x);
  2108.         sc -= 8;
  2109.         }
  2110.  
  2111.     while( sc > 0 )
  2112.         {
  2113.         eshup1(x);
  2114.         sc -= 1;
  2115.         }
  2116.     }
  2117. if( lost )
  2118.     lost = 1;
  2119. return( (int )lost );
  2120. }
  2121.  
  2122.  
  2123.  
  2124. /*
  2125. ;    normalize
  2126. ;
  2127. ; Shift normalizes the significand area pointed to by argument
  2128. ; shift count (up = positive) is returned.
  2129. */
  2130. int enormlz(x)
  2131. unsigned short x[];
  2132. {
  2133. register unsigned short *p;
  2134. int sc;
  2135.  
  2136. sc = 0;
  2137. p = &x[M];
  2138. if( *p != 0 )
  2139.     goto normdn;
  2140. ++p;
  2141. if( *p & 0x8000 )
  2142.     return( 0 );    /* already normalized */
  2143. while( *p == 0 )
  2144.     {
  2145.     eshup6(x);
  2146.     sc += 16;
  2147. /* With guard word, there are NBITS+16 bits available.
  2148.  * return true if all are zero.
  2149.  */
  2150.     if( sc > NBITS )
  2151.         return( sc );
  2152.     }
  2153. /* see if high byte is zero */
  2154. while( (*p & 0xff00) == 0 )
  2155.     {
  2156.     eshup8(x);
  2157.     sc += 8;
  2158.     }
  2159. /* now shift 1 bit at a time */
  2160. while( (*p  & 0x8000) == 0)
  2161.     {
  2162.     eshup1(x);
  2163.     sc += 1;
  2164.     if( sc > NBITS )
  2165.         {
  2166.         mtherr( "enormlz", UNDERFLOW );
  2167.         return( sc );
  2168.         }
  2169.     }
  2170. return( sc );
  2171.  
  2172. /* Normalize by shifting down out of the high guard word
  2173.    of the significand */
  2174. normdn:
  2175.  
  2176. if( *p & 0xff00 )
  2177.     {
  2178.     eshdn8(x);
  2179.     sc -= 8;
  2180.     }
  2181. while( *p != 0 )
  2182.     {
  2183.     eshdn1(x);
  2184.     sc -= 1;
  2185.  
  2186.     if( sc < -NBITS )
  2187.         {
  2188.         mtherr( "enormlz", OVERFLOW );
  2189.         return( sc );
  2190.         }
  2191.     }
  2192. return( sc );
  2193. }
  2194.  
  2195.  
  2196.  
  2197.  
  2198. /* Convert e type number to decimal format ASCII string.
  2199.  * The constants are for 64 bit precision.
  2200.  */
  2201.  
  2202. #define NTEN 12
  2203. #define MAXP 4096
  2204.  
  2205. static unsigned short etens[NTEN+1][NE] = {
  2206. {0xc94c,0x979a,0x8a20,0x5202,0xc460,0x7525,},/* 10**4096 */
  2207. {0xa74d,0x5de4,0xc53d,0x3b5d,0x9e8b,0x5a92,},/* 10**2048 */
  2208. {0x650d,0x0c17,0x8175,0x7586,0xc976,0x4d48,},
  2209. {0xcc65,0x91c6,0xa60e,0xa0ae,0xe319,0x46a3,},
  2210. {0xddbc,0xde8d,0x9df9,0xebfb,0xaa7e,0x4351,},
  2211. {0xc66f,0x8cdf,0x80e9,0x47c9,0x93ba,0x41a8,},
  2212. {0x3cbf,0xa6d5,0xffcf,0x1f49,0xc278,0x40d3,},
  2213. {0xf020,0xb59d,0x2b70,0xada8,0x9dc5,0x4069,},
  2214. {0x0000,0x0000,0x0400,0xc9bf,0x8e1b,0x4034,},
  2215. {0x0000,0x0000,0x0000,0x2000,0xbebc,0x4019,},
  2216. {0x0000,0x0000,0x0000,0x0000,0x9c40,0x400c,},
  2217. {0x0000,0x0000,0x0000,0x0000,0xc800,0x4005,},
  2218. {0x0000,0x0000,0x0000,0x0000,0xa000,0x4002,}, /* 10**1 */
  2219. };
  2220.  
  2221. static unsigned short emtens[NTEN+1][NE] = {
  2222. {0x2de4,0x9fde,0xd2ce,0x04c8,0xa6dd,0x0ad8,}, /* 10**-4096 */
  2223. {0x4925,0x2de4,0x3436,0x534f,0xceae,0x256b,}, /* 10**-2048 */
  2224. {0x87a6,0xc0bd,0xda57,0x82a5,0xa2a6,0x32b5,},
  2225. {0x7133,0xd21c,0xdb23,0xee32,0x9049,0x395a,},
  2226. {0xfa91,0x1939,0x637a,0x4325,0xc031,0x3cac,},
  2227. {0xac7d,0xe4a0,0x64bc,0x467c,0xddd0,0x3e55,},
  2228. {0x3f24,0xe9a5,0xa539,0xea27,0xa87f,0x3f2a,},
  2229. {0x67de,0x94ba,0x4539,0x1ead,0xcfb1,0x3f94,},
  2230. {0x4c2f,0xe15b,0xc44d,0x94be,0xe695,0x3fc9,},
  2231. {0xfdc2,0xcefc,0x8461,0x7711,0xabcc,0x3fe4,},
  2232. {0xd3c3,0x652b,0xe219,0x1758,0xd1b7,0x3ff1,},
  2233. {0x3d71,0xd70a,0x70a3,0x0a3d,0xa3d7,0x3ff8,},
  2234. {0xcccd,0xcccc,0xcccc,0xcccc,0xcccc,0x3ffb,}, /* 10**-1 */
  2235. };
  2236.  
  2237. void e24toasc( x, string, ndigs )
  2238. unsigned short x[];
  2239. char *string;
  2240. int ndigs;
  2241. {
  2242. unsigned short w[NI];
  2243.  
  2244. #ifdef INFINITY
  2245. #ifdef IBMPC
  2246. if( (x[1] & 0x7f80) == 0x7f80 )
  2247. #else
  2248. if( (x[0] & 0x7f80) == 0x7f80 )
  2249. #endif
  2250.     {
  2251. #ifdef IBMPC
  2252.     if( x[1] & 0x8000 )
  2253. #else
  2254.     if( x[0] & 0x8000 )
  2255. #endif
  2256.         sprintf( string, " -Infinity " );
  2257.     else
  2258.         sprintf( string, " Infinity " );
  2259.     return;
  2260.     }
  2261. #endif
  2262. e24toe( x, w );
  2263. etoasc( w, string, ndigs );
  2264. }
  2265.  
  2266.  
  2267. void e53toasc( x, string, ndigs )
  2268. unsigned short x[];
  2269. char *string;
  2270. int ndigs;
  2271. {
  2272. unsigned short w[NI];
  2273.  
  2274. #ifdef INFINITY
  2275. #ifdef IBMPC
  2276. if( (x[3] & 0x7ff0) == 0x7ff0 )
  2277. #else
  2278. if( (x[0] & 0x7ff0) == 0x7ff0 )
  2279. #endif
  2280.     {
  2281. #ifdef IBMPC
  2282.     if( x[3] & 0x8000 )
  2283. #else
  2284.     if( x[0] & 0x8000 )
  2285. #endif
  2286.         sprintf( string, " -Infinity " );
  2287.     else
  2288.         sprintf( string, " Infinity " );
  2289.     return;
  2290.     }
  2291. #endif
  2292. e53toe( x, w );
  2293. etoasc( w, string, ndigs );
  2294. }
  2295.  
  2296.  
  2297. void e64toasc( x, string, ndigs )
  2298. unsigned short x[];
  2299. char *string;
  2300. int ndigs;
  2301. {
  2302. unsigned short w[NI];
  2303.  
  2304. #ifdef INFINITY
  2305. #ifdef IBMPC
  2306. if( (x[4] & 0x7fff) == 0x7fff )
  2307. #else
  2308. if( (x[0] & 0x7fff) == 0x7fff )
  2309. #endif
  2310.     {
  2311. #ifdef IBMPC
  2312.     if( x[4] & 0x8000 )
  2313. #else
  2314.     if( x[0] & 0x8000 )
  2315. #endif
  2316.         sprintf( string, " -Infinity " );
  2317.     else
  2318.         sprintf( string, " Infinity " );
  2319.     return;
  2320.     }
  2321. #endif
  2322. e64toe( x, w );
  2323. etoasc( w, string, ndigs );
  2324. }
  2325.  
  2326.  
  2327.  
  2328. void etoasc( x, string, ndigs )
  2329. unsigned short x[];
  2330. char *string;
  2331. int ndigs;
  2332. {
  2333. long digit;
  2334. unsigned short y[NI], t[NI], u[NI], w[NI];
  2335. unsigned short *p, *r, *ten;
  2336. unsigned short sign;
  2337. int i, j, k, expon, rndsav;
  2338. char *s, *ss;
  2339. unsigned short m;
  2340.  
  2341. rndsav = rndprc;
  2342. rndprc = NBITS;        /* set to full precision */
  2343. emov( x, y ); /* retain external format */
  2344. if( y[NE-1] & 0x8000 )
  2345.     {
  2346.     sign = 0xffff;
  2347.     y[NE-1] &= 0x7fff;
  2348.     }
  2349. else
  2350.     {
  2351.     sign = 0;
  2352.     }
  2353. expon = 0;
  2354. ten = &etens[NTEN][0];
  2355. emov( eone, t );
  2356. /* Test for zero exponent */
  2357. if( y[NE-1] == 0 )
  2358.     {
  2359.     for( k=0; k<NE-1; k++ )
  2360.         {
  2361.         if( y[k] != 0 )
  2362.             goto tnzro; /* denormalized number */
  2363.         }
  2364.     goto isone; /* legal all zeros */
  2365.     }
  2366. tnzro:
  2367.  
  2368. /* Test for infinity.  Don't bother with illegal infinities.
  2369.  */
  2370. if( y[NE-1] == 0x7fff )
  2371.     {
  2372.     if( sign )
  2373.         sprintf( string, " -Infinity " );
  2374.     else
  2375.         sprintf( string, " Infinity " );
  2376.     goto bxit;
  2377.     }
  2378.  
  2379. /* Test for exponent nonzero but significand denormalized.
  2380.  * This is an error condition.
  2381.  */
  2382. if( (y[NE-1] != 0) && ((y[NE-2] & 0x8000) == 0) )
  2383.     {
  2384.     mtherr( "etoasc", DOMAIN );
  2385.     sprintf( string, "NaN" );
  2386.     goto bxit;
  2387.     }
  2388.  
  2389. /* Compare to 1.0 */
  2390. i = ecmp( eone, y );
  2391. if( i == 0 )
  2392.     goto isone;
  2393.  
  2394. if( i < 0 )
  2395.     { /* Number is greater than 1 */
  2396. /* Convert significand to an integer and strip trailing decimal zeros. */
  2397.     emov( y, u );
  2398.     u[NE-1] = EXONE + NBITS - 1;
  2399.  
  2400.     p = &etens[NTEN-4][0];
  2401.     m = 16;
  2402. do
  2403.     {
  2404.     ediv( p, u, t );
  2405.     efloor( t, w );
  2406.     for( j=0; j<NE-1; j++ )
  2407.         {
  2408.         if( t[j] != w[j] )
  2409.             goto noint;
  2410.         }
  2411.     emov( t, u );
  2412.     expon += (int )m;
  2413. noint:
  2414.     p += NE;
  2415.     m >>= 1;
  2416.     }
  2417. while( m != 0 );
  2418.  
  2419. /* Rescale from integer significand */
  2420.     u[NE-1] += y[NE-1] - (unsigned int )(EXONE + NBITS - 1);
  2421.     emov( u, y );
  2422. /* Find power of 10 */
  2423.     emov( eone, t );
  2424.     m = MAXP;
  2425.     p = &etens[0][0];
  2426.     while( ecmp( ten, u ) <= 0 )
  2427.         {
  2428.         if( ecmp( p, u ) <= 0 )
  2429.             {
  2430.             ediv( p, u, u );
  2431.             emul( p, t, t );
  2432.             expon += (int )m;
  2433.             }
  2434.         m >>= 1;
  2435.         if( m == 0 )
  2436.             break;
  2437.         p += NE;
  2438.         }
  2439.     }
  2440. else
  2441.     { /* Number is less than 1.0 */
  2442. /* Pad significand with trailing decimal zeros. */
  2443.     if( y[NE-1] == 0 )
  2444.         {
  2445.         while( (y[NE-2] & 0x8000) == 0 )
  2446.             {
  2447.             emul( ten, y, y );
  2448.             expon -= 1;
  2449.             }
  2450.         }
  2451.     else
  2452.         {
  2453.         emovi( y, w );
  2454.         for( i=0; i<NDEC+1; i++ )
  2455.             {
  2456.             if( (w[NI-1] & 0x7) != 0 )
  2457.                 break;
  2458. /* multiply by 10 */
  2459.             emovz( w, u );
  2460.             eshdn1( u );
  2461.             eshdn1( u );
  2462.             eaddm( w, u );
  2463.             u[1] += 3;
  2464.             while( u[2] != 0 )
  2465.                 {
  2466.                 eshdn1(u);
  2467.                 u[1] += 1;
  2468.                 }
  2469.             if( u[NI-1] != 0 )
  2470.                 break;
  2471.             if( eone[NE-1] <= u[1] )
  2472.                 break;
  2473.             emovz( u, w );
  2474.             expon -= 1;
  2475.             }
  2476.         emovo( w, y );
  2477.         }
  2478.     k = -MAXP;
  2479.     p = &emtens[0][0];
  2480.     r = &etens[0][0];
  2481.     emov( y, w );
  2482.     emov( eone, t );
  2483.     while( ecmp( eone, w ) > 0 )
  2484.         {
  2485.         if( ecmp( p, w ) >= 0 )
  2486.             {
  2487.             emul( r, w, w );
  2488.             emul( r, t, t );
  2489.             expon += k;
  2490.             }
  2491.         k /= 2;
  2492.         if( k == 0 )
  2493.             break;
  2494.         p += NE;
  2495.         r += NE;
  2496.         }
  2497.     ediv( t, eone, t );
  2498.     }
  2499. isone:
  2500. /* Find the first (leading) digit. */
  2501. emovi( t, w );
  2502. emovz( w, t );
  2503. emovi( y, w );
  2504. emovz( w, y );
  2505. eiremain( t, y );
  2506. digit = equot[NI-1];
  2507. while( (digit == 0) && (ecmp(y,ezero) != 0) )
  2508.     {
  2509.     eshup1( y );
  2510.     emovz( y, u );
  2511.     eshup1( u );
  2512.     eshup1( u );
  2513.     eaddm( u, y );
  2514.     eiremain( t, y );
  2515.     digit = equot[NI-1];
  2516.     expon -= 1;
  2517.     }
  2518. s = string;
  2519. if( sign )
  2520.     *s++ = '-';
  2521. else
  2522.     *s++ = ' ';
  2523. *s++ = (char )digit + '0';
  2524. *s++ = '.';
  2525. /* Examine number of digits requested by caller. */
  2526. if( ndigs < 0 )
  2527.     ndigs = 0;
  2528. if( ndigs > NDEC )
  2529.     ndigs = NDEC;
  2530. /* Generate digits after the decimal point. */
  2531. for( k=0; k<=ndigs; k++ )
  2532.     {
  2533. /* multiply current number by 10, without normalizing */
  2534.     eshup1( y );
  2535.     emovz( y, u );
  2536.     eshup1( u );
  2537.     eshup1( u );
  2538.     eaddm( u, y );
  2539.     eiremain( t, y );
  2540.     *s++ = (char )equot[NI-1] + '0';
  2541.     }
  2542. digit = equot[NI-1];
  2543. --s;
  2544. ss = s;
  2545. /* round off the ASCII string */
  2546. if( digit > 4 )
  2547.     {
  2548. /* Test for critical rounding case in ASCII output. */
  2549.     if( digit == 5 )
  2550.         {
  2551.         emovo( y, t );
  2552.         if( ecmp(t,ezero) != 0 )
  2553.             goto roun;    /* round to nearest */
  2554.         if( (*(s-1) & 1) == 0 )
  2555.             goto doexp;    /* round to even */
  2556.         }
  2557. /* Round up and propagate carry-outs */
  2558. roun:
  2559.     --s;
  2560.     k = *s & 0x7f;
  2561. /* Carry out to most significant digit? */
  2562.     if( k == '.' )
  2563.         {
  2564.         --s;
  2565.         k = *s;
  2566.         k += 1;
  2567.         *s = (char )k;
  2568. /* Most significant digit carries to 10? */
  2569.         if( k > '9' )
  2570.             {
  2571.             expon += 1;
  2572.             *s = '1';
  2573.             }
  2574.         goto doexp;
  2575.         }
  2576. /* Round up and carry out from less significant digits */
  2577.     k += 1;
  2578.     *s = (char )k;
  2579.     if( k > '9' )
  2580.         {
  2581.         *s = '0';
  2582.         goto roun;
  2583.         }
  2584.     }
  2585. doexp:
  2586. /*
  2587. if( expon >= 0 )
  2588.     sprintf( ss, "e+%d", expon );
  2589. else
  2590.     sprintf( ss, "e%d", expon );
  2591. */
  2592.     sprintf( ss, "E%d", expon );
  2593. bxit:
  2594. rndprc = rndsav;
  2595. }
  2596.  
  2597.  
  2598.  
  2599.  
  2600. /*
  2601. ;                                ASCTOQ
  2602. ;        ASCTOQ.MAC        LATEST REV: 11 JAN 84
  2603. ;                    SLM, 3 JAN 78
  2604. ;
  2605. ;    Convert ASCII string to quadruple precision floating point
  2606. ;
  2607. ;        Numeric input is free field decimal number
  2608. ;        with max of 15 digits with or without 
  2609. ;        decimal point entered as ASCII from teletype.
  2610. ;    Entering E after the number followed by a second
  2611. ;    number causes the second number to be interpreted
  2612. ;    as a power of 10 to be multiplied by the first number
  2613. ;    (i.e., "scientific" notation).
  2614. ;
  2615. ;    Usage:
  2616. ;        asctoq( string, q );
  2617. */
  2618.  
  2619. /* ASCII to single */
  2620. void asctoe24( s, y )
  2621. char *s;
  2622. unsigned short *y;
  2623. {
  2624. asctoeg( s, y, 24 );
  2625. }
  2626.  
  2627.  
  2628. /* ASCII to double */
  2629. void asctoe53( s, y )
  2630. char *s;
  2631. unsigned short *y;
  2632. {
  2633. #ifdef DEC
  2634. asctoeg( s, y, 56 );
  2635. #else
  2636. asctoeg( s, y, 53 );
  2637. #endif
  2638. }
  2639.  
  2640.  
  2641. /* ASCII to long double */
  2642. void asctoe64( s, y )
  2643. char *s;
  2644. unsigned short *y;
  2645. {
  2646. asctoeg( s, y, 64 );
  2647. }
  2648.  
  2649. /* ASCII to super double */
  2650. void asctoe( s, y )
  2651. char *s;
  2652. unsigned short *y;
  2653. {
  2654. asctoeg( s, y, NBITS );
  2655. }
  2656.  
  2657. /* Space to make a copy of the input string: */
  2658. static char lstr[82] = {0};
  2659.  
  2660. void asctoeg( ss, y, oprec )
  2661. char *ss;
  2662. unsigned short *y;
  2663. int oprec;
  2664. {
  2665. unsigned short yy[NI], xt[NI], tt[NI];
  2666. int esign, decflg, sgnflg, nexp, exp, prec, lost;
  2667. int k, trail, c, rndsav;
  2668. long lexp;
  2669. unsigned short nsign, *p;
  2670. char *sp, *s;
  2671.  
  2672. /* Copy the input string. */
  2673. s = ss;
  2674. while( *s == ' ' ) /* skip leading spaces */
  2675.     ++s;
  2676. sp = lstr;
  2677. for( k=0; k<79; k++ )
  2678.     {
  2679.     if( (*sp++ = *s++) == '\0' )
  2680.         break;
  2681.     }
  2682. *sp = '\0';
  2683. s = lstr;
  2684.  
  2685. rndsav = rndprc;
  2686. rndprc = NBITS; /* Set to full precision */
  2687. lost = 0;
  2688. nsign = 0;
  2689. decflg = 0;
  2690. sgnflg = 0;
  2691. nexp = 0;
  2692. exp = 0;
  2693. prec = 0;
  2694. ecleaz( yy );
  2695. trail = 0;
  2696.  
  2697. nxtcom:
  2698. k = *s - '0';
  2699. if( (k >= 0) && (k <= 9) )
  2700.     {
  2701. /* Ignore leading zeros */
  2702.     if( (prec == 0) && (decflg == 0) && (k == 0) )
  2703.         goto donchr;
  2704. /* Identify and strip trailing zeros after the decimal point. */
  2705.     if( (trail == 0) && (decflg != 0) )
  2706.         {
  2707.         sp = s;
  2708.         while( (*sp >= '0') && (*sp <= '9') )
  2709.             ++sp;
  2710. /* Check for syntax error */
  2711.         c = *sp & 0x7f;
  2712.         if( (c != 'e') && (c != 'E') && (c != '\0')
  2713.             && (c != '\n') && (c != '\r') && (c != ' ')
  2714.             && (c != ',') )
  2715.             goto error;
  2716.         --sp;
  2717.         while( *sp == '0' )
  2718.             *sp-- = 'z';
  2719.         trail = 1;
  2720.         if( *s == 'z' )
  2721.             goto donchr;
  2722.         }
  2723. /* If enough digits were given to more than fill up the yy register,
  2724.  * continuing until overflow into the high guard word yy[2]
  2725.  * guarantees that there will be a roundoff bit at the top
  2726.  * of the low guard word after normalization.
  2727.  */
  2728.     if( yy[2] == 0 )
  2729.         {
  2730.         if( decflg )
  2731.             nexp += 1; /* count digits after decimal point */
  2732.         eshup1( yy );    /* multiply current number by 10 */
  2733.         emovz( yy, xt );
  2734.         eshup1( xt );
  2735.         eshup1( xt );
  2736.         eaddm( xt, yy );
  2737.         ecleaz( xt );
  2738.         xt[NI-2] = (unsigned short )k;
  2739.         eaddm( xt, yy );
  2740.         }
  2741.     else
  2742.         {
  2743.         lost |= k;
  2744.         }
  2745.     prec += 1;
  2746.     goto donchr;
  2747.     }
  2748.  
  2749. switch( *s )
  2750.     {
  2751.     case 'z':
  2752.         break;
  2753.     case 'E':
  2754.     case 'e':
  2755.         goto expnt;
  2756.     case '.':    /* decimal point */
  2757.         if( decflg )
  2758.             goto error;
  2759.         ++decflg;
  2760.         break;
  2761.     case '-':
  2762.         nsign = 0xffff;
  2763.         if( sgnflg )
  2764.             goto error;
  2765.         ++sgnflg;
  2766.         break;
  2767.     case '+':
  2768.         if( sgnflg )
  2769.             goto error;
  2770.         ++sgnflg;
  2771.         break;
  2772.     case ',':
  2773.     case ' ':
  2774.     case '\0':
  2775.     case '\n':
  2776.     case '\r':
  2777.         goto daldone;
  2778.     case 'i':
  2779.     case 'I':
  2780.         ecleaz(yy);
  2781.         yy[E] = 0x7fff;  /* infinity */
  2782.         goto aexit;
  2783.     default:
  2784.     error:
  2785.         mtherr( "asctoe", DOMAIN );
  2786.         eclear(y);
  2787.         goto aexit;
  2788.     }
  2789. donchr:
  2790. ++s;
  2791. goto nxtcom;
  2792.  
  2793. /* Exponent interpretation */
  2794. expnt:
  2795.  
  2796. esign = 1;
  2797. exp = 0;
  2798. ++s;
  2799. /* check for + or - */
  2800. if( *s == '-' )
  2801.     {
  2802.     esign = -1;
  2803.     ++s;
  2804.     }
  2805. if( *s == '+' )
  2806.     ++s;
  2807. while( (*s >= '0') && (*s <= '9') )
  2808.     {
  2809.     exp *= 10;
  2810.     exp += *s++ - '0';
  2811.     }
  2812. if( esign < 0 )
  2813.     exp = -exp;
  2814.  
  2815. daldone:
  2816. nexp = exp - nexp;
  2817. /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
  2818. while( (nexp > 0) && (yy[2] == 0) )
  2819.     {
  2820.     emovz( yy, xt );
  2821.     eshup1( xt );
  2822.     eshup1( xt );
  2823.     eaddm( yy, xt );
  2824.     eshup1( xt );
  2825.     if( xt[2] != 0 )
  2826.         break;
  2827.     nexp -= 1;
  2828.     emovz( xt, yy );
  2829.     }
  2830. if( (k = enormlz(yy)) > NBITS )
  2831.     {
  2832.     ecleaz(yy);
  2833.     goto aexit;
  2834.     }
  2835. lexp = (EXONE - 1 + NBITS) - k;
  2836. emdnorm( yy, lost, 0, lexp, 64 );
  2837. /* convert to external format */
  2838.  
  2839.  
  2840. /* Multiply by 10**nexp.  If precision is 64 bits,
  2841.  * the maximum relative error incurred in forming 10**n
  2842.  * for 0 <= n <= 324 is 8.2e-20, at 10**180.
  2843.  * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
  2844.  * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
  2845.  */
  2846. lexp = yy[E];
  2847. if( nexp == 0 )
  2848.     {
  2849.     k = 0;
  2850.     goto expdon;
  2851.     }
  2852. esign = 1;
  2853. if( nexp < 0 )
  2854.     {
  2855.     nexp = -nexp;
  2856.     esign = -1;
  2857.     if( nexp > 4096 )
  2858.         { /* Punt.  Can't handle this without 2 divides. */
  2859.         emovi( etens[0], tt );
  2860.         lexp -= tt[E];
  2861.         k = edivm( tt, yy );
  2862.         lexp += EXONE;
  2863.         nexp -= 4096;
  2864.         }
  2865.     }
  2866. p = &etens[NTEN][0];
  2867. emov( eone, xt );
  2868. exp = 1;
  2869. do
  2870.     {
  2871.     if( exp & nexp )
  2872.         emul( p, xt, xt );
  2873.     p -= NE;
  2874.     exp = exp + exp;
  2875.     }
  2876. while( exp <= MAXP );
  2877.  
  2878. emovi( xt, tt );
  2879. if( esign < 0 )
  2880.     {
  2881.     lexp -= tt[E];
  2882.     k = edivm( tt, yy );
  2883.     lexp += EXONE;
  2884.     }
  2885. else
  2886.     {
  2887.     lexp += tt[E];
  2888.     k = emulm( tt, yy );
  2889.     lexp -= EXONE - 1;
  2890.     }
  2891.  
  2892. expdon:
  2893.  
  2894. /* Round and convert directly to the destination type */
  2895. if( oprec == 53 )
  2896.     lexp -= EXONE - 0x3ff;
  2897. else if( oprec == 24 )
  2898.     lexp -= EXONE - 0177;
  2899. #ifdef DEC
  2900. else if( oprec == 56 )
  2901.     lexp -= EXONE - 0201;
  2902. #endif
  2903. rndprc = oprec;
  2904. emdnorm( yy, k, 0, lexp, 64 );
  2905.  
  2906. aexit:
  2907.  
  2908. rndprc = rndsav;
  2909. yy[0] = nsign;
  2910. switch( oprec )
  2911.     {
  2912. #ifdef DEC
  2913.     case 56:
  2914.         todec( yy, y ); /* see etodec.c */
  2915.         break;
  2916. #endif
  2917.     case 53:
  2918.         toe53( yy, y );
  2919.         break;
  2920.     case 24:
  2921.         toe24( yy, y );
  2922.         break;
  2923.     case 64:
  2924.         toe64( yy, y );
  2925.         break;
  2926.     case NBITS:
  2927.         emovo( yy, y );
  2928.         break;
  2929.     }
  2930. }
  2931.  
  2932.  
  2933.  
  2934. /* y = largest integer not greater than x
  2935.  * (truncated toward minus infinity)
  2936.  *
  2937.  * unsigned short x[NE], y[NE]
  2938.  *
  2939.  * efloor( x, y );
  2940.  */
  2941. static unsigned short bmask[] = {
  2942. 0xffff,
  2943. 0xfffe,
  2944. 0xfffc,
  2945. 0xfff8,
  2946. 0xfff0,
  2947. 0xffe0,
  2948. 0xffc0,
  2949. 0xff80,
  2950. 0xff00,
  2951. 0xfe00,
  2952. 0xfc00,
  2953. 0xf800,
  2954. 0xf000,
  2955. 0xe000,
  2956. 0xc000,
  2957. 0x8000,
  2958. 0x0000,
  2959. };
  2960.  
  2961. void efloor( x, y )
  2962. unsigned short x[], y[];
  2963. {
  2964. register unsigned short *p;
  2965. int e, expon, i;
  2966. unsigned short f[NE];
  2967.  
  2968. emov( x, f ); /* leave in external format */
  2969. expon = (int )f[NE-1];
  2970. e = (expon & 0x7fff) - (EXONE - 1);
  2971. if( e <= 0 )
  2972.     {
  2973.     eclear(y);
  2974.     goto isitneg;
  2975.     }
  2976. /* number of bits to clear out */
  2977. e = NBITS - e;
  2978. emov( f, y );
  2979. if( e <= 0 )
  2980.     return;
  2981.  
  2982. p = &y[0];
  2983. while( e >= 16 )
  2984.     {
  2985.     *p++ = 0;
  2986.     e -= 16;
  2987.     }
  2988. /* clear the remaining bits */
  2989. *p &= bmask[e];
  2990. /* truncate negatives toward minus infinity */
  2991. isitneg:
  2992.  
  2993. if( (unsigned short )expon & (unsigned short )0x8000 )
  2994.     {
  2995.     for( i=0; i<NE-1; i++ )
  2996.         {
  2997.         if( f[i] != y[i] )
  2998.             {
  2999.             esub( eone, y, y );
  3000.             break;
  3001.             }
  3002.         }
  3003.     }
  3004. }
  3005.  
  3006.  
  3007. /* unsigned short x[], s[];
  3008.  * long *exp;
  3009.  *
  3010.  * efrexp( x, exp, s );
  3011.  *
  3012.  * Returns s and exp such that  s * 2**exp = x and .5 <= s < 1.
  3013.  * For example, 1.1 = 0.55 * 2**1
  3014.  * Handles denormalized numbers properly using long integer exp.
  3015.  */
  3016. void efrexp( x, exp, s )
  3017. unsigned short x[];
  3018. long *exp;
  3019. unsigned short s[];
  3020. {
  3021. unsigned short xi[NI];
  3022. long li;
  3023.  
  3024. emovi( x, xi );
  3025. li = (long )((short )xi[1]);
  3026.  
  3027. if( li == 0 )
  3028.     {
  3029.     li -= enormlz( xi );
  3030.     }
  3031. xi[1] = 0x3ffe;
  3032. emovo( xi, s );
  3033. *exp = li - 0x3ffe;
  3034. }
  3035.  
  3036.  
  3037.  
  3038. /* unsigned short x[], y[];
  3039.  * long pwr2;
  3040.  *
  3041.  * eldexp( x, pwr2, y );
  3042.  *
  3043.  * Returns y = x * 2**pwr2.
  3044.  */
  3045. void eldexp( x, pwr2, y )
  3046. unsigned short x[];
  3047. long pwr2;
  3048. unsigned short y[];
  3049. {
  3050. unsigned short xi[NI];
  3051. long li;
  3052. int i;
  3053.  
  3054. emovi( x, xi );
  3055. li = xi[1];
  3056. li += pwr2;
  3057. i = 0;
  3058. emdnorm( xi, i, i, li, 64 );
  3059. emovo( xi, y );
  3060. }
  3061.  
  3062.  
  3063. /* c = remainder after dividing b by a
  3064.  * Least significant integer quotient bits left in equot[].
  3065.  */
  3066. void eremain( a, b, c )
  3067. unsigned short a[], b[], c[];
  3068. {
  3069. unsigned short den[NI], num[NI];
  3070.  
  3071. if( ecmp(a,ezero) == 0 )
  3072.     {
  3073.     mtherr( "eremain", SING );
  3074.     eclear( c );
  3075.     return;
  3076.     }
  3077. emovi( a, den );
  3078. emovi( b, num );
  3079. eiremain( den, num );
  3080. /* Sign of remainder = sign of quotient */
  3081. if( a[0] == b[0] )
  3082.     num[0] = 0;
  3083. else
  3084.     num[0] = 0xffff;
  3085. emovo( num, c );
  3086. }
  3087.  
  3088. void eiremain( den, num )
  3089. unsigned short den[], num[];
  3090. {
  3091. long ld, ln;
  3092. unsigned short j;
  3093.  
  3094. ld = den[E];
  3095. ld -= enormlz( den );
  3096. ln = num[E];
  3097. ln -= enormlz( num );
  3098. ecleaz( equot );
  3099. while( ln >= ld )
  3100.     {
  3101.     if( ecmpm(den,num) <= 0 )
  3102.         {
  3103.         esubm(den, num);
  3104.         j = 1;
  3105.         }
  3106.     else
  3107.         {
  3108.         j = 0;
  3109.         }
  3110.     eshup1(equot);
  3111.     equot[NI-1] |= j;
  3112.     eshup1(num);
  3113.     ln -= 1;
  3114.     }
  3115. emdnorm( num, 0, 0, ln, 0 );
  3116. }
  3117.  
  3118.