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