Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  1.  /* Extended precision arithmetic functions for long double I/O.
  2.   * This program has been placed in the public domain.
  3.   */
  4.  
  5. #include <_ansi.h>
  6. #include <reent.h>
  7. #include <string.h>
  8. #include <stdlib.h>
  9. #include "mprec.h"
  10.  
  11. /* These are the externally visible entries. */
  12. /* linux name:  long double _IO_strtold (char *, char **); */
  13. long double _strtold (char *, char **);
  14. char *_ldtoa_r (struct _reent *, long double, int, int, int *, int *,
  15.                 char **);
  16. int _ldcheck (long double *);
  17. #if 0
  18. void _IO_ldtostr (long double *, char *, int, int, char);
  19. #endif
  20.  
  21.  /* Number of 16 bit words in external x type format */
  22. #define NE 10
  23.  
  24.  /* Number of 16 bit words in internal format */
  25. #define NI (NE+3)
  26.  
  27.  /* Array offset to exponent */
  28. #define E 1
  29.  
  30.  /* Array offset to high guard word */
  31. #define M 2
  32.  
  33.  /* Number of bits of precision */
  34. #define NBITS ((NI-4)*16)
  35.  
  36.  /* Maximum number of decimal digits in ASCII conversion
  37.   * = NBITS*log10(2)
  38.   */
  39. #define NDEC (NBITS*8/27)
  40.  
  41.  /* The exponent of 1.0 */
  42. #define EXONE (0x3fff)
  43.  
  44.  /* Maximum exponent digits - base 10 */
  45. #define MAX_EXP_DIGITS 5
  46.  
  47. /* Control structure for long double conversion including rounding precision values.
  48.  * rndprc can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
  49.  */
  50. typedef struct
  51. {
  52.   int rlast;
  53.   int rndprc;
  54.   int rw;
  55.   int re;
  56.   int outexpon;
  57.   unsigned short rmsk;
  58.   unsigned short rmbit;
  59.   unsigned short rebit;
  60.   unsigned short rbit[NI];
  61.   unsigned short equot[NI];
  62. } LDPARMS;
  63.  
  64. static void esub (_CONST short unsigned int *a, _CONST short unsigned int *b,
  65.                   short unsigned int *c, LDPARMS * ldp);
  66. static void emul (_CONST short unsigned int *a, _CONST short unsigned int *b,
  67.                   short unsigned int *c, LDPARMS * ldp);
  68. static void ediv (_CONST short unsigned int *a, _CONST short unsigned int *b,
  69.                   short unsigned int *c, LDPARMS * ldp);
  70. static int ecmp (_CONST short unsigned int *a, _CONST short unsigned int *b);
  71. static int enormlz (short unsigned int *x);
  72. static int eshift (short unsigned int *x, int sc);
  73. static void eshup1 (register short unsigned int *x);
  74. static void eshup8 (register short unsigned int *x);
  75. static void eshup6 (register short unsigned int *x);
  76. static void eshdn1 (register short unsigned int *x);
  77. static void eshdn8 (register short unsigned int *x);
  78. static void eshdn6 (register short unsigned int *x);
  79. static void eneg (short unsigned int *x);
  80. static void emov (register _CONST short unsigned int *a,
  81.                   register short unsigned int *b);
  82. static void eclear (register short unsigned int *x);
  83. static void einfin (register short unsigned int *x, register LDPARMS * ldp);
  84. static void efloor (short unsigned int *x, short unsigned int *y,
  85.                     LDPARMS * ldp);
  86. static void etoasc (short unsigned int *x, char *string, int ndigs,
  87.                     int outformat, LDPARMS * ldp);
  88.  
  89. union uconv
  90. {
  91.   unsigned short pe;
  92.   long double d;
  93. };
  94.  
  95. #if LDBL_MANT_DIG == 24
  96. static void e24toe (short unsigned int *pe, short unsigned int *y,
  97.                     LDPARMS * ldp);
  98. #elif LDBL_MANT_DIG == 53
  99. static void e53toe (short unsigned int *pe, short unsigned int *y,
  100.                     LDPARMS * ldp);
  101. #elif LDBL_MANT_DIG == 64
  102. static void e64toe (short unsigned int *pe, short unsigned int *y,
  103.                     LDPARMS * ldp);
  104. #else
  105. static void e113toe (short unsigned int *pe, short unsigned int *y,
  106.                      LDPARMS * ldp);
  107. #endif
  108.  
  109. /*                                                      econst.c        */
  110. /*  e type constants used by high precision check routines */
  111.  
  112. #if NE == 10
  113. /* 0.0 */
  114. static _CONST unsigned short ezero[NE] = { 0x0000, 0x0000, 0x0000, 0x0000,
  115.   0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,
  116. };
  117.  
  118. /* 1.0E0 */
  119. static _CONST unsigned short eone[NE] = { 0x0000, 0x0000, 0x0000, 0x0000,
  120.   0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,
  121. };
  122.  
  123. #else
  124.  
  125. /* 0.0 */
  126. static _CONST unsigned short ezero[NE] = {
  127.   0, 0000000, 0000000, 0000000, 0000000, 0000000,
  128. };
  129.  
  130. /* 1.0E0 */
  131. static _CONST unsigned short eone[NE] = {
  132.   0, 0000000, 0000000, 0000000, 0100000, 0x3fff,
  133. };
  134.  
  135. #endif
  136.  
  137. /* Debugging routine for displaying errors */
  138. #ifdef DEBUG
  139. /* Notice: the order of appearance of the following
  140.  * messages is bound to the error codes defined
  141.  * in mconf.h.
  142.  */
  143. static _CONST char *_CONST ermsg[7] = {
  144.   "unknown",                    /* error code 0 */
  145.   "domain",                     /* error code 1 */
  146.   "singularity",                /* et seq.      */
  147.   "overflow",
  148.   "underflow",
  149.   "total loss of precision",
  150.   "partial loss of precision"
  151. };
  152.  
  153. #define mtherr(name, code) printf( "\n%s %s error\n", name, ermsg[code] );
  154. #else
  155. #define mtherr(name, code)
  156. #endif
  157.  
  158. /*                                                      ieee.c
  159.  *
  160.  *    Extended precision IEEE binary floating point arithmetic routines
  161.  *
  162.  * Numbers are stored in C language as arrays of 16-bit unsigned
  163.  * short integers.  The arguments of the routines are pointers to
  164.  * the arrays.
  165.  *
  166.  *
  167.  * External e type data structure, simulates Intel 8087 chip
  168.  * temporary real format but possibly with a larger significand:
  169.  *
  170.  *      NE-1 significand words  (least significant word first,
  171.  *                               most significant bit is normally set)
  172.  *      exponent                (value = EXONE for 1.0,
  173.  *                              top bit is the sign)
  174.  *
  175.  *
  176.  * Internal data structure of a number (a "word" is 16 bits):
  177.  *
  178.  * ei[0]        sign word       (0 for positive, 0xffff for negative)
  179.  * ei[1]        biased exponent (value = EXONE for the number 1.0)
  180.  * ei[2]        high guard word (always zero after normalization)
  181.  * ei[3]
  182.  * to ei[NI-2]  significand     (NI-4 significand words,
  183.  *                               most significant word first,
  184.  *                               most significant bit is set)
  185.  * ei[NI-1]     low guard word  (0x8000 bit is rounding place)
  186.  *
  187.  *
  188.  *
  189.  *              Routines for external format numbers
  190.  *
  191.  *      asctoe( string, e )     ASCII string to extended double e type
  192.  *      asctoe64( string, &d )  ASCII string to long double
  193.  *      asctoe53( string, &d )  ASCII string to double
  194.  *      asctoe24( string, &f )  ASCII string to single
  195.  *      asctoeg( string, e, prec, ldp ) ASCII string to specified precision
  196.  *      e24toe( &f, e, ldp )    IEEE single precision to e type
  197.  *      e53toe( &d, e, ldp )    IEEE double precision to e type
  198.  *      e64toe( &d, e, ldp )    IEEE long double precision to e type
  199.  *      e113toe( &d, e, ldp )   IEEE long double precision to e type
  200.  *      eabs(e)                 absolute value
  201.  *      eadd( a, b, c )         c = b + a
  202.  *      eclear(e)               e = 0
  203.  *      ecmp (a, b)             Returns 1 if a > b, 0 if a == b,
  204.  *                              -1 if a < b, -2 if either a or b is a NaN.
  205.  *      ediv( a, b, c, ldp )    c = b / a
  206.  *      efloor( a, b, ldp )     truncate to integer, toward -infinity
  207.  *      efrexp( a, exp, s )     extract exponent and significand
  208.  *      eifrac( e, &l, frac )   e to long integer and e type fraction
  209.  *      euifrac( e, &l, frac )  e to unsigned long integer and e type fraction
  210.  *      einfin( e, ldp )        set e to infinity, leaving its sign alone
  211.  *      eldexp( a, n, b )       multiply by 2**n
  212.  *      emov( a, b )            b = a
  213.  *      emul( a, b, c, ldp )    c = b * a
  214.  *      eneg(e)                 e = -e
  215.  *      eround( a, b )          b = nearest integer value to a
  216.  *      esub( a, b, c, ldp )    c = b - a
  217.  *      e24toasc( &f, str, n )  single to ASCII string, n digits after decimal
  218.  *      e53toasc( &d, str, n )  double to ASCII string, n digits after decimal
  219.  *      e64toasc( &d, str, n )  long double to ASCII string
  220.  *      etoasc(e,str,n,fmt,ldp)e to ASCII string, n digits after decimal
  221.  *      etoe24( e, &f )         convert e type to IEEE single precision
  222.  *      etoe53( e, &d )         convert e type to IEEE double precision
  223.  *      etoe64( e, &d )         convert e type to IEEE long double precision
  224.  *      ltoe( &l, e )           long (32 bit) integer to e type
  225.  *      ultoe( &l, e )          unsigned long (32 bit) integer to e type
  226.  *      eisneg( e )             1 if sign bit of e != 0, else 0
  227.  *      eisinf( e )             1 if e has maximum exponent (non-IEEE)
  228.  *                              or is infinite (IEEE)
  229.  *      eisnan( e )             1 if e is a NaN
  230.  *      esqrt( a, b )           b = square root of a
  231.  *
  232.  *
  233.  *              Routines for internal format numbers
  234.  *
  235.  *      eaddm( ai, bi )         add significands, bi = bi + ai
  236.  *      ecleaz(ei)              ei = 0
  237.  *      ecleazs(ei)             set ei = 0 but leave its sign alone
  238.  *      ecmpm( ai, bi )         compare significands, return 1, 0, or -1
  239.  *      edivm( ai, bi, ldp )    divide  significands, bi = bi / ai
  240.  *      emdnorm(ai,l,s,exp,ldp) normalize and round off
  241.  *      emovi( a, ai )          convert external a to internal ai
  242.  *      emovo( ai, a, ldp )     convert internal ai to external a
  243.  *      emovz( ai, bi )         bi = ai, low guard word of bi = 0
  244.  *      emulm( ai, bi, ldp )    multiply significands, bi = bi * ai
  245.  *      enormlz(ei)             left-justify the significand
  246.  *      eshdn1( ai )            shift significand and guards down 1 bit
  247.  *      eshdn8( ai )            shift down 8 bits
  248.  *      eshdn6( ai )            shift down 16 bits
  249.  *      eshift( ai, n )         shift ai n bits up (or down if n < 0)
  250.  *      eshup1( ai )            shift significand and guards up 1 bit
  251.  *      eshup8( ai )            shift up 8 bits
  252.  *      eshup6( ai )            shift up 16 bits
  253.  *      esubm( ai, bi )         subtract significands, bi = bi - ai
  254.  *
  255.  *
  256.  * The result is always normalized and rounded to NI-4 word precision
  257.  * after each arithmetic operation.
  258.  *
  259.  * Exception flags are NOT fully supported.
  260.  *
  261.  * Define USE_INFINITY in mconf.h for support of infinity; otherwise a
  262.  * saturation arithmetic is implemented.
  263.  *
  264.  * Define NANS for support of Not-a-Number items; otherwise the
  265.  * arithmetic will never produce a NaN output, and might be confused
  266.  * by a NaN input.
  267.  * If NaN's are supported, the output of ecmp(a,b) is -2 if
  268.  * either a or b is a NaN. This means asking if(ecmp(a,b) < 0)
  269.  * may not be legitimate. Use if(ecmp(a,b) == -1) for less-than
  270.  * if in doubt.
  271.  * Signaling NaN's are NOT supported; they are treated the same
  272.  * as quiet NaN's.
  273.  *
  274.  * Denormals are always supported here where appropriate (e.g., not
  275.  * for conversion to DEC numbers).
  276.  */
  277.  
  278. /*
  279.  * Revision history:
  280.  *
  281.  *  5 Jan 84    PDP-11 assembly language version
  282.  *  6 Dec 86    C language version
  283.  * 30 Aug 88    100 digit version, improved rounding
  284.  * 15 May 92    80-bit long double support
  285.  * 22 Nov 00    Revised to fit into newlib by Jeff Johnston <jjohnstn@redhat.com>
  286.  *
  287.  * Author:  S. L. Moshier.
  288.  *
  289.  * Copyright (c) 1984,2000 S.L. Moshier
  290.  *
  291.  * Permission to use, copy, modify, and distribute this software for any
  292.  * purpose without fee is hereby granted, provided that this entire notice
  293.  * is included in all copies of any software which is or includes a copy
  294.  * or modification of this software and in all copies of the supporting
  295.  * documentation for such software.
  296.  *
  297.  * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED
  298.  * WARRANTY.  IN PARTICULAR,  THE AUTHOR MAKES NO REPRESENTATION
  299.  * OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY OF THIS
  300.  * SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.
  301.  *
  302.  */
  303.  
  304. #include <stdio.h>
  305. /* #include "\usr\include\stdio.h" */
  306. /*#include "ehead.h"*/
  307. /*#include "mconf.h"*/
  308. /*                                                      mconf.h
  309.  *
  310.  *      Common include file for math routines
  311.  *
  312.  *
  313.  *
  314.  * SYNOPSIS:
  315.  *
  316.  * #include "mconf.h"
  317.  *
  318.  *
  319.  *
  320.  * DESCRIPTION:
  321.  *
  322.  * This file contains definitions for error codes that are
  323.  * passed to the common error handling routine mtherr()
  324.  * (which see).
  325.  *
  326.  * The file also includes a conditional assembly definition
  327.  * for the type of computer arithmetic (IEEE, DEC, Motorola
  328.  * IEEE, or UNKnown).
  329.  *
  330.  * For Digital Equipment PDP-11 and VAX computers, certain
  331.  * IBM systems, and others that use numbers with a 56-bit
  332.  * significand, the symbol DEC should be defined.  In this
  333.  * mode, most floating point constants are given as arrays
  334.  * of octal integers to eliminate decimal to binary conversion
  335.  * errors that might be introduced by the compiler.
  336.  *
  337.  * For computers, such as IBM PC, that follow the IEEE
  338.  * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
  339.  * Std 754-1985), the symbol IBMPC should be defined.  These
  340.  * numbers have 53-bit significands.  In this mode, constants
  341.  * are provided as arrays of hexadecimal 16 bit integers.
  342.  *
  343.  * To accommodate other types of computer arithmetic, all
  344.  * constants are also provided in a normal decimal radix
  345.  * which one can hope are correctly converted to a suitable
  346.  * format by the available C language compiler.  To invoke
  347.  * this mode, the symbol UNK is defined.
  348.  *
  349.  * An important difference among these modes is a predefined
  350.  * set of machine arithmetic constants for each.  The numbers
  351.  * MACHEP (the machine roundoff error), MAXNUM (largest number
  352.  * represented), and several other parameters are preset by
  353.  * the configuration symbol.  Check the file const.c to
  354.  * ensure that these values are correct for your computer.
  355.  *
  356.  * For ANSI C compatibility, define ANSIC equal to 1.  Currently
  357.  * this affects only the atan2() function and others that use it.
  358.  */
  359.  
  360. /* Constant definitions for math error conditions
  361.  */
  362.  
  363. #define DOMAIN          1       /* argument domain error */
  364. #define SING            2       /* argument singularity */
  365. #define OVERFLOW        3       /* overflow range error */
  366. #define UNDERFLOW       4       /* underflow range error */
  367. #define TLOSS           5       /* total loss of precision */
  368. #define PLOSS           6       /* partial loss of precision */
  369.  
  370. #define EDOM            33
  371. #define ERANGE          34
  372.  
  373. typedef struct
  374. {
  375.   double r;
  376.   double i;
  377. } cmplx;
  378.  
  379. /* Type of computer arithmetic */
  380.  
  381. #ifndef DEC
  382. #ifdef __IEEE_LITTLE_ENDIAN
  383. #define IBMPC 1
  384. #else /* !__IEEE_LITTLE_ENDIAN */
  385. #define MIEEE 1
  386. #endif /* !__IEEE_LITTLE_ENDIAN */
  387. #endif /* !DEC */
  388.  
  389. /* Define 1 for ANSI C atan2() function
  390.  * See atan.c and clog.c.
  391.  */
  392. #define ANSIC 1
  393.  
  394. /*define VOLATILE volatile*/
  395. #define VOLATILE
  396.  
  397. #define NANS
  398. #define USE_INFINITY
  399.  
  400. /* NaN's require infinity support. */
  401. #ifdef NANS
  402. #ifndef INFINITY
  403. #define USE_INFINITY
  404. #endif
  405. #endif
  406.  
  407. /* This handles 64-bit long ints. */
  408. #define LONGBITS (8 * sizeof(long))
  409.  
  410.  
  411. static void eaddm (short unsigned int *x, short unsigned int *y);
  412. static void esubm (short unsigned int *x, short unsigned int *y);
  413. static void emdnorm (short unsigned int *s, int lost, int subflg,
  414.                      long int exp, int rcntrl, LDPARMS * ldp);
  415. static int asctoeg (char *ss, short unsigned int *y, int oprec,
  416.                     LDPARMS * ldp);
  417. static void enan (short unsigned int *nan, int size);
  418. #if LDBL_MANT_DIG == 24
  419. static void toe24 (short unsigned int *x, short unsigned int *y);
  420. #elif LDBL_MANT_DIG == 53
  421. static void toe53 (short unsigned int *x, short unsigned int *y);
  422. #elif LDBL_MANT_DIG == 64
  423. static void toe64 (short unsigned int *a, short unsigned int *b);
  424. #else
  425. static void toe113 (short unsigned int *a, short unsigned int *b);
  426. #endif
  427. static void eiremain (short unsigned int *den, short unsigned int *num,
  428.                       LDPARMS * ldp);
  429. static int ecmpm (register short unsigned int *a,
  430.                   register short unsigned int *b);
  431. static int edivm (short unsigned int *den, short unsigned int *num,
  432.                   LDPARMS * ldp);
  433. static int emulm (short unsigned int *a, short unsigned int *b,
  434.                   LDPARMS * ldp);
  435. static int eisneg (_CONST short unsigned int *x);
  436. static int eisinf (_CONST short unsigned int *x);
  437. static void emovi (_CONST short unsigned int *a, short unsigned int *b);
  438. static void emovo (short unsigned int *a, short unsigned int *b,
  439.                    LDPARMS * ldp);
  440. static void emovz (register short unsigned int *a,
  441.                    register short unsigned int *b);
  442. static void ecleaz (register short unsigned int *xi);
  443. static void eadd1 (_CONST short unsigned int *a, _CONST short unsigned int *b,
  444.                    short unsigned int *c, int subflg, LDPARMS * ldp);
  445. static int eisnan (_CONST short unsigned int *x);
  446. static int eiisnan (short unsigned int *x);
  447.  
  448. #ifdef DEC
  449. static void etodec (), todec (), dectoe ();
  450. #endif
  451.  
  452. /*
  453. ; Clear out entire external format number.
  454. ;
  455. ; unsigned short x[];
  456. ; eclear( x );
  457. */
  458.  
  459. static void
  460. eclear (register short unsigned int *x)
  461. {
  462.   register int i;
  463.  
  464.   for (i = 0; i < NE; i++)
  465.     *x++ = 0;
  466. }
  467.  
  468.  
  469.  
  470. /* Move external format number from a to b.
  471.  *
  472.  * emov( a, b );
  473.  */
  474.  
  475. static void
  476. emov (register _CONST short unsigned int *a, register short unsigned int *b)
  477. {
  478.   register int i;
  479.  
  480.   for (i = 0; i < NE; i++)
  481.     *b++ = *a++;
  482. }
  483.  
  484.  
  485. /*
  486. ;       Negate external format number
  487. ;
  488. ;       unsigned short x[NE];
  489. ;       eneg( x );
  490. */
  491.  
  492. static void
  493. eneg (short unsigned int *x)
  494. {
  495.  
  496. #ifdef NANS
  497.   if (eisnan (x))
  498.     return;
  499. #endif
  500.   x[NE - 1] ^= 0x8000;          /* Toggle the sign bit */
  501. }
  502.  
  503.  
  504.  
  505. /* Return 1 if external format number is negative,
  506.  * else return zero.
  507.  */
  508. static int
  509. eisneg (_CONST short unsigned int *x)
  510. {
  511.  
  512. #ifdef NANS
  513.   if (eisnan (x))
  514.     return (0);
  515. #endif
  516.   if (x[NE - 1] & 0x8000)
  517.     return (1);
  518.   else
  519.     return (0);
  520. }
  521.  
  522.  
  523. /* Return 1 if external format number has maximum possible exponent,
  524.  * else return zero.
  525.  */
  526. static int
  527. eisinf (_CONST short unsigned int *x)
  528. {
  529.  
  530.   if ((x[NE - 1] & 0x7fff) == 0x7fff)
  531.     {
  532. #ifdef NANS
  533.       if (eisnan (x))
  534.         return (0);
  535. #endif
  536.       return (1);
  537.     }
  538.   else
  539.     return (0);
  540. }
  541.  
  542. /* Check if e-type number is not a number.
  543.  */
  544. static int
  545. eisnan (_CONST short unsigned int *x)
  546. {
  547.  
  548. #ifdef NANS
  549.   int i;
  550. /* NaN has maximum exponent */
  551.   if ((x[NE - 1] & 0x7fff) != 0x7fff)
  552.     return (0);
  553. /* ... and non-zero significand field. */
  554.   for (i = 0; i < NE - 1; i++)
  555.     {
  556.       if (*x++ != 0)
  557.         return (1);
  558.     }
  559. #endif
  560.   return (0);
  561. }
  562.  
  563. /*
  564. ; Fill entire number, including exponent and significand, with
  565. ; largest possible number.  These programs implement a saturation
  566. ; value that is an ordinary, legal number.  A special value
  567. ; "infinity" may also be implemented; this would require tests
  568. ; for that value and implementation of special rules for arithmetic
  569. ; operations involving inifinity.
  570. */
  571.  
  572. static void
  573. einfin (register short unsigned int *x, register LDPARMS * ldp)
  574. {
  575.   register int i;
  576.  
  577. #ifdef USE_INFINITY
  578.   for (i = 0; i < NE - 1; i++)
  579.     *x++ = 0;
  580.   *x |= 32767;
  581.   ldp = ldp;
  582. #else
  583.   for (i = 0; i < NE - 1; i++)
  584.     *x++ = 0xffff;
  585.   *x |= 32766;
  586.   if (ldp->rndprc < NBITS)
  587.     {
  588.       if (ldp->rndprc == 113)
  589.         {
  590.           *(x - 9) = 0;
  591.           *(x - 8) = 0;
  592.         }
  593.       if (ldp->rndprc == 64)
  594.         {
  595.           *(x - 5) = 0;
  596.         }
  597.       if (ldp->rndprc == 53)
  598.         {
  599.           *(x - 4) = 0xf800;
  600.         }
  601.       else
  602.         {
  603.           *(x - 4) = 0;
  604.           *(x - 3) = 0;
  605.           *(x - 2) = 0xff00;
  606.         }
  607.     }
  608. #endif
  609. }
  610.  
  611. /* Move in external format number,
  612.  * converting it to internal format.
  613.  */
  614. static void
  615. emovi (_CONST short unsigned int *a, short unsigned int *b)
  616. {
  617.   register _CONST unsigned short *p;
  618.   register unsigned short *q;
  619.   int i;
  620.  
  621.   q = b;
  622.   p = a + (NE - 1);             /* point to last word of external number */
  623. /* get the sign bit */
  624.   if (*p & 0x8000)
  625.     *q++ = 0xffff;
  626.   else
  627.     *q++ = 0;
  628. /* get the exponent */
  629.   *q = *p--;
  630.   *q++ &= 0x7fff;               /* delete the sign bit */
  631. #ifdef USE_INFINITY
  632.   if ((*(q - 1) & 0x7fff) == 0x7fff)
  633.     {
  634. #ifdef NANS
  635.       if (eisnan (a))
  636.         {
  637.           *q++ = 0;
  638.           for (i = 3; i < NI; i++)
  639.             *q++ = *p--;
  640.           return;
  641.         }
  642. #endif
  643.       for (i = 2; i < NI; i++)
  644.         *q++ = 0;
  645.       return;
  646.     }
  647. #endif
  648. /* clear high guard word */
  649.   *q++ = 0;
  650. /* move in the significand */
  651.   for (i = 0; i < NE - 1; i++)
  652.     *q++ = *p--;
  653. /* clear low guard word */
  654.   *q = 0;
  655. }
  656.  
  657.  
  658. /* Move internal format number out,
  659.  * converting it to external format.
  660.  */
  661. static void
  662. emovo (short unsigned int *a, short unsigned int *b, LDPARMS * ldp)
  663. {
  664.   register unsigned short *p, *q;
  665.   unsigned short i;
  666.  
  667.   p = a;
  668.   q = b + (NE - 1);             /* point to output exponent */
  669. /* combine sign and exponent */
  670.   i = *p++;
  671.   if (i)
  672.     *q-- = *p++ | 0x8000;
  673.   else
  674.     *q-- = *p++;
  675. #ifdef USE_INFINITY
  676.   if (*(p - 1) == 0x7fff)
  677.     {
  678. #ifdef NANS
  679.       if (eiisnan (a))
  680.         {
  681.           enan (b, NBITS);
  682.           return;
  683.         }
  684. #endif
  685.       einfin (b, ldp);
  686.       return;
  687.     }
  688. #endif
  689. /* skip over guard word */
  690.   ++p;
  691. /* move the significand */
  692.   for (i = 0; i < NE - 1; i++)
  693.     *q-- = *p++;
  694. }
  695.  
  696.  
  697. /* Clear out internal format number.
  698.  */
  699.  
  700. static void
  701. ecleaz (register short unsigned int *xi)
  702. {
  703.   register int i;
  704.  
  705.   for (i = 0; i < NI; i++)
  706.     *xi++ = 0;
  707. }
  708.  
  709. /* same, but don't touch the sign. */
  710.  
  711. static void
  712. ecleazs (register short unsigned int *xi)
  713. {
  714.   register int i;
  715.  
  716.   ++xi;
  717.   for (i = 0; i < NI - 1; i++)
  718.     *xi++ = 0;
  719. }
  720.  
  721.  
  722.  
  723.  
  724. /* Move internal format number from a to b.
  725.  */
  726. static void
  727. emovz (register short unsigned int *a, register short unsigned int *b)
  728. {
  729.   register int i;
  730.  
  731.   for (i = 0; i < NI - 1; i++)
  732.     *b++ = *a++;
  733. /* clear low guard word */
  734.   *b = 0;
  735. }
  736.  
  737. /* Return nonzero if internal format number is a NaN.
  738.  */
  739.  
  740. static int
  741. eiisnan (short unsigned int *x)
  742. {
  743.   int i;
  744.  
  745.   if ((x[E] & 0x7fff) == 0x7fff)
  746.     {
  747.       for (i = M + 1; i < NI; i++)
  748.         {
  749.           if (x[i] != 0)
  750.             return (1);
  751.         }
  752.     }
  753.   return (0);
  754. }
  755.  
  756. #if LDBL_MANT_DIG == 64
  757.  
  758. /* Return nonzero if internal format number is infinite. */
  759. static int
  760. eiisinf (unsigned short x[])
  761. {
  762.  
  763. #ifdef NANS
  764.   if (eiisnan (x))
  765.     return (0);
  766. #endif
  767.   if ((x[E] & 0x7fff) == 0x7fff)
  768.     return (1);
  769.   return (0);
  770. }
  771. #endif /* LDBL_MANT_DIG == 64 */
  772.  
  773. /*
  774. ;       Compare significands of numbers in internal format.
  775. ;       Guard words are included in the comparison.
  776. ;
  777. ;       unsigned short a[NI], b[NI];
  778. ;       cmpm( a, b );
  779. ;
  780. ;       for the significands:
  781. ;       returns +1 if a > b
  782. ;                0 if a == b
  783. ;               -1 if a < b
  784. */
  785. static int
  786. ecmpm (register short unsigned int *a, register short unsigned int *b)
  787. {
  788.   int i;
  789.  
  790.   a += M;                       /* skip up to significand area */
  791.   b += M;
  792.   for (i = M; i < NI; i++)
  793.     {
  794.       if (*a++ != *b++)
  795.         goto difrnt;
  796.     }
  797.   return (0);
  798.  
  799. difrnt:
  800.   if (*(--a) > *(--b))
  801.     return (1);
  802.   else
  803.     return (-1);
  804. }
  805.  
  806.  
  807. /*
  808. ;       Shift significand down by 1 bit
  809. */
  810.  
  811. static void
  812. eshdn1 (register short unsigned int *x)
  813. {
  814.   register unsigned short bits;
  815.   int i;
  816.  
  817.   x += M;                       /* point to significand area */
  818.  
  819.   bits = 0;
  820.   for (i = M; i < NI; i++)
  821.     {
  822.       if (*x & 1)
  823.         bits |= 1;
  824.       *x >>= 1;
  825.       if (bits & 2)
  826.         *x |= 0x8000;
  827.       bits <<= 1;
  828.       ++x;
  829.     }
  830. }
  831.  
  832.  
  833.  
  834. /*
  835. ;       Shift significand up by 1 bit
  836. */
  837.  
  838. static void
  839. eshup1 (register short unsigned int *x)
  840. {
  841.   register unsigned short bits;
  842.   int i;
  843.  
  844.   x += NI - 1;
  845.   bits = 0;
  846.  
  847.   for (i = M; i < NI; i++)
  848.     {
  849.       if (*x & 0x8000)
  850.         bits |= 1;
  851.       *x <<= 1;
  852.       if (bits & 2)
  853.         *x |= 1;
  854.       bits <<= 1;
  855.       --x;
  856.     }
  857. }
  858.  
  859.  
  860.  
  861. /*
  862. ;       Shift significand down by 8 bits
  863. */
  864.  
  865. static void
  866. eshdn8 (register short unsigned int *x)
  867. {
  868.   register unsigned short newbyt, oldbyt;
  869.   int i;
  870.  
  871.   x += M;
  872.   oldbyt = 0;
  873.   for (i = M; i < NI; i++)
  874.     {
  875.       newbyt = *x << 8;
  876.       *x >>= 8;
  877.       *x |= oldbyt;
  878.       oldbyt = newbyt;
  879.       ++x;
  880.     }
  881. }
  882.  
  883. /*
  884. ;       Shift significand up by 8 bits
  885. */
  886.  
  887. static void
  888. eshup8 (register short unsigned int *x)
  889. {
  890.   int i;
  891.   register unsigned short newbyt, oldbyt;
  892.  
  893.   x += NI - 1;
  894.   oldbyt = 0;
  895.  
  896.   for (i = M; i < NI; i++)
  897.     {
  898.       newbyt = *x >> 8;
  899.       *x <<= 8;
  900.       *x |= oldbyt;
  901.       oldbyt = newbyt;
  902.       --x;
  903.     }
  904. }
  905.  
  906. /*
  907. ;       Shift significand up by 16 bits
  908. */
  909.  
  910. static void
  911. eshup6 (register short unsigned int *x)
  912. {
  913.   int i;
  914.   register unsigned short *p;
  915.  
  916.   p = x + M;
  917.   x += M + 1;
  918.  
  919.   for (i = M; i < NI - 1; i++)
  920.     *p++ = *x++;
  921.  
  922.   *p = 0;
  923. }
  924.  
  925. /*
  926. ;       Shift significand down by 16 bits
  927. */
  928.  
  929. static void
  930. eshdn6 (register short unsigned int *x)
  931. {
  932.   int i;
  933.   register unsigned short *p;
  934.  
  935.   x += NI - 1;
  936.   p = x + 1;
  937.  
  938.   for (i = M; i < NI - 1; i++)
  939.     *(--p) = *(--x);
  940.  
  941.   *(--p) = 0;
  942. }
  944. /*
  945. ;       Add significands
  946. ;       x + y replaces y
  947. */
  948.  
  949. static void
  950. eaddm (short unsigned int *x, short unsigned int *y)
  951. {
  952.   register unsigned long a;
  953.   int i;
  954.   unsigned int carry;
  955.  
  956.   x += NI - 1;
  957.   y += NI - 1;
  958.   carry = 0;
  959.   for (i = M; i < NI; i++)
  960.     {
  961.       a = (unsigned long) (*x) + (unsigned long) (*y) + carry;
  962.       if (a & 0x10000)
  963.         carry = 1;
  964.       else
  965.         carry = 0;
  966.       *y = (unsigned short) a;
  967.       --x;
  968.       --y;
  969.     }
  970. }
  971.  
  972. /*
  973. ;       Subtract significands
  974. ;       y - x replaces y
  975. */
  976.  
  977. static void
  978. esubm (short unsigned int *x, short unsigned int *y)
  979. {
  980.   unsigned long a;
  981.   int i;
  982.   unsigned int carry;
  983.  
  984.   x += NI - 1;
  985.   y += NI - 1;
  986.   carry = 0;
  987.   for (i = M; i < NI; i++)
  988.     {
  989.       a = (unsigned long) (*y) - (unsigned long) (*x) - carry;
  990.       if (a & 0x10000)
  991.         carry = 1;
  992.       else
  993.         carry = 0;
  994.       *y = (unsigned short) a;
  995.       --x;
  996.       --y;
  997.     }
  998. }
  999.  
  1000.  
  1001. /* Divide significands */
  1002.  
  1003.  
  1004. /* Multiply significand of e-type number b
  1005. by 16-bit quantity a, e-type result to c. */
  1006.  
  1007. static void
  1008. m16m (short unsigned int a, short unsigned int *b, short unsigned int *c)
  1009. {
  1010.   register unsigned short *pp;
  1011.   register unsigned long carry;
  1012.   unsigned short *ps;
  1013.   unsigned short p[NI];
  1014.   unsigned long aa, m;
  1015.   int i;
  1016.  
  1017.   aa = a;
  1018.   pp = &p[NI - 2];
  1019.   *pp++ = 0;
  1020.   *pp = 0;
  1021.   ps = &b[NI - 1];
  1022.  
  1023.   for (i = M + 1; i < NI; i++)
  1024.     {
  1025.       if (*ps == 0)
  1026.         {
  1027.           --ps;
  1028.           --pp;
  1029.           *(pp - 1) = 0;
  1030.         }
  1031.       else
  1032.         {
  1033.           m = (unsigned long) aa **ps--;
  1034.           carry = (m & 0xffff) + *pp;
  1035.           *pp-- = (unsigned short) carry;
  1036.           carry = (carry >> 16) + (m >> 16) + *pp;
  1037.           *pp = (unsigned short) carry;
  1038.           *(pp - 1) = carry >> 16;
  1039.         }
  1040.     }
  1041.   for (i = M; i < NI; i++)
  1042.     c[i] = p[i];
  1043. }
  1044.  
  1045.  
  1046. /* Divide significands. Neither the numerator nor the denominator
  1047. is permitted to have its high guard word nonzero.  */
  1048.  
  1049.  
  1050. static int
  1051. edivm (short unsigned int *den, short unsigned int *num, LDPARMS * ldp)
  1052. {
  1053.   int i;
  1054.   register unsigned short *p;
  1055.   unsigned long tnum;
  1056.   unsigned short j, tdenm, tquot;
  1057.   unsigned short tprod[NI + 1];
  1058.   unsigned short *equot = ldp->equot;
  1059.  
  1060.   p = &equot[0];
  1061.   *p++ = num[0];
  1062.   *p++ = num[1];
  1063.  
  1064.   for (i = M; i < NI; i++)
  1065.     {
  1066.       *p++ = 0;
  1067.     }
  1068.   eshdn1 (num);
  1069.   tdenm = den[M + 1];
  1070.   for (i = M; i < NI; i++)
  1071.     {
  1072.       /* Find trial quotient digit (the radix is 65536). */
  1073.       tnum = (((unsigned long) num[M]) << 16) + num[M + 1];
  1074.  
  1075.       /* Do not execute the divide instruction if it will overflow. */
  1076.       if ((tdenm * 0xffffUL) < tnum)
  1077.         tquot = 0xffff;
  1078.       else
  1079.         tquot = tnum / tdenm;
  1080.  
  1081.       /* Prove that the divide worked. */
  1082. /*
  1083.         tcheck = (unsigned long )tquot * tdenm;
  1084.         if( tnum - tcheck > tdenm )
  1085.                 tquot = 0xffff;
  1086. */
  1087.       /* Multiply denominator by trial quotient digit. */
  1088.       m16m (tquot, den, tprod);
  1089.       /* The quotient digit may have been overestimated. */
  1090.       if (ecmpm (tprod, num) > 0)
  1091.         {
  1092.           tquot -= 1;
  1093.           esubm (den, tprod);
  1094.           if (ecmpm (tprod, num) > 0)
  1095.             {
  1096.               tquot -= 1;
  1097.               esubm (den, tprod);
  1098.             }
  1099.         }
  1100. /*
  1101.         if( ecmpm( tprod, num ) > 0 )
  1102.                 {
  1103.                 eshow( "tprod", tprod );
  1104.                 eshow( "num  ", num );
  1105.                 printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
  1106.                          tnum, den[M+1], tquot );
  1107.                 }
  1108. */
  1109.       esubm (tprod, num);
  1110. /*
  1111.         if( ecmpm( num, den ) >= 0 )
  1112.                 {
  1113.                 eshow( "num  ", num );
  1114.                 eshow( "den  ", den );
  1115.                 printf( "tnum = %08lx, tden = %04x, tquot = %04x\n",
  1116.                          tnum, den[M+1], tquot );
  1117.                 }
  1118. */
  1119.       equot[i] = tquot;
  1120.       eshup6 (num);
  1121.     }
  1122. /* test for nonzero remainder after roundoff bit */
  1123.   p = &num[M];
  1124.   j = 0;
  1125.   for (i = M; i < NI; i++)
  1126.     {
  1127.       j |= *p++;
  1128.     }
  1129.   if (j)
  1130.     j = 1;
  1131.  
  1132.   for (i = 0; i < NI; i++)
  1133.     num[i] = equot[i];
  1134.  
  1135.   return ((int) j);
  1136. }
  1137.  
  1138.  
  1139.  
  1140. /* Multiply significands */
  1141. static int
  1142. emulm (short unsigned int *a, short unsigned int *b, LDPARMS * ldp)
  1143. {
  1144.   unsigned short *p, *q;
  1145.   unsigned short pprod[NI];
  1146.   unsigned short j;
  1147.   int i;
  1148.   unsigned short *equot = ldp->equot;
  1149.  
  1150.   equot[0] = b[0];
  1151.   equot[1] = b[1];
  1152.   for (i = M; i < NI; i++)
  1153.     equot[i] = 0;
  1154.  
  1155.   j = 0;
  1156.   p = &a[NI - 1];
  1157.   q = &equot[NI - 1];
  1158.   for (i = M + 1; i < NI; i++)
  1159.     {
  1160.       if (*p == 0)
  1161.         {
  1162.           --p;
  1163.         }
  1164.       else
  1165.         {
  1166.           m16m (*p--, b, pprod);
  1167.           eaddm (pprod, equot);
  1168.         }
  1169.       j |= *q;
  1170.       eshdn6 (equot);
  1171.     }
  1172.  
  1173.   for (i = 0; i < NI; i++)
  1174.     b[i] = equot[i];
  1175.  
  1176. /* return flag for lost nonzero bits */
  1177.   return ((int) j);
  1178. }
  1179.  
  1180.  
  1181. /*
  1182. static void eshow(str, x)
  1183. char *str;
  1184. unsigned short *x;
  1185. {
  1186. int i;
  1187.  
  1188. printf( "%s ", str );
  1189. for( i=0; i<NI; i++ )
  1190.         printf( "%04x ", *x++ );
  1191. printf( "\n" );
  1192. }
  1193. */
  1194.  
  1195.  
  1196. /*
  1197.  * Normalize and round off.
  1198.  *
  1199.  * The internal format number to be rounded is "s".
  1200.  * Input "lost" indicates whether the number is exact.
  1201.  * This is the so-called sticky bit.
  1202.  *
  1203.  * Input "subflg" indicates whether the number was obtained
  1204.  * by a subtraction operation.  In that case if lost is nonzero
  1205.  * then the number is slightly smaller than indicated.
  1206.  *
  1207.  * Input "exp" is the biased exponent, which may be negative.
  1208.  * the exponent field of "s" is ignored but is replaced by
  1209.  * "exp" as adjusted by normalization and rounding.
  1210.  *
  1211.  * Input "rcntrl" is the rounding control.
  1212.  */
  1213.  
  1214.  
  1215. static void
  1216. emdnorm (short unsigned int *s, int lost, int subflg, long int exp,
  1217.          int rcntrl, LDPARMS * ldp)
  1218. {
  1219.   int i, j;
  1220.   unsigned short r;
  1221.  
  1222. /* Normalize */
  1223.   j = enormlz (s);
  1224.  
  1225. /* a blank significand could mean either zero or infinity. */
  1226. #ifndef USE_INFINITY
  1227.   if (j > NBITS)
  1228.     {
  1229.       ecleazs (s);
  1230.       return;
  1231.     }
  1232. #endif
  1233.   exp -= j;
  1234. #ifndef USE_INFINITY
  1235.   if (exp >= 32767L)
  1236.     goto overf;
  1237. #else
  1238.   if ((j > NBITS) && (exp < 32767L))
  1239.     {
  1240.       ecleazs (s);
  1241.       return;
  1242.     }
  1243. #endif
  1244.   if (exp < 0L)
  1245.     {
  1246.       if (exp > (long) (-NBITS - 1))
  1247.         {
  1248.           j = (int) exp;
  1249.           i = eshift (s, j);
  1250.           if (i)
  1251.             lost = 1;
  1252.         }
  1253.       else
  1254.         {
  1255.           ecleazs (s);
  1256.           return;
  1257.         }
  1258.     }
  1259. /* Round off, unless told not to by rcntrl. */
  1260.   if (rcntrl == 0)
  1261.     goto mdfin;
  1262. /* Set up rounding parameters if the control register changed. */
  1263.   if (ldp->rndprc != ldp->rlast)
  1264.     {
  1265.       ecleaz (ldp->rbit);
  1266.       switch (ldp->rndprc)
  1267.         {
  1268.         default:
  1269.         case NBITS:
  1270.           ldp->rw = NI - 1;     /* low guard word */
  1271.           ldp->rmsk = 0xffff;
  1272.           ldp->rmbit = 0x8000;
  1273.           ldp->rebit = 1;
  1274.           ldp->re = ldp->rw - 1;
  1275.           break;
  1276.         case 113:
  1277.           ldp->rw = 10;
  1278.           ldp->rmsk = 0x7fff;
  1279.           ldp->rmbit = 0x4000;
  1280.           ldp->rebit = 0x8000;
  1281.           ldp->re = ldp->rw;
  1282.           break;
  1283.         case 64:
  1284.           ldp->rw = 7;
  1285.           ldp->rmsk = 0xffff;
  1286.           ldp->rmbit = 0x8000;
  1287.           ldp->rebit = 1;
  1288.           ldp->re = ldp->rw - 1;
  1289.           break;
  1290. /* For DEC arithmetic */
  1291.         case 56:
  1292.           ldp->rw = 6;
  1293.           ldp->rmsk = 0xff;
  1294.           ldp->rmbit = 0x80;
  1295.           ldp->rebit = 0x100;
  1296.           ldp->re = ldp->rw;
  1297.           break;
  1298.         case 53:
  1299.           ldp->rw = 6;
  1300.           ldp->rmsk = 0x7ff;
  1301.           ldp->rmbit = 0x0400;
  1302.           ldp->rebit = 0x800;
  1303.           ldp->re = ldp->rw;
  1304.           break;
  1305.         case 24:
  1306.           ldp->rw = 4;
  1307.           ldp->rmsk = 0xff;
  1308.           ldp->rmbit = 0x80;
  1309.           ldp->rebit = 0x100;
  1310.           ldp->re = ldp->rw;
  1311.           break;
  1312.         }
  1313.       ldp->rbit[ldp->re] = ldp->rebit;
  1314.       ldp->rlast = ldp->rndprc;
  1315.     }
  1316.  
  1317. /* Shift down 1 temporarily if the data structure has an implied
  1318.  * most significant bit and the number is denormal.
  1319.  * For rndprc = 64 or NBITS, there is no implied bit.
  1320.  * But Intel long double denormals lose one bit of significance even so.
  1321.  */
  1322. #if IBMPC
  1323.   if ((exp <= 0) && (ldp->rndprc != NBITS))
  1324. #else
  1325.   if ((exp <= 0) && (ldp->rndprc != 64) && (ldp->rndprc != NBITS))
  1326. #endif
  1327.     {
  1328.       lost |= s[NI - 1] & 1;
  1329.       eshdn1 (s);
  1330.     }
  1331. /* Clear out all bits below the rounding bit,
  1332.  * remembering in r if any were nonzero.
  1333.  */
  1334.   r = s[ldp->rw] & ldp->rmsk;
  1335.   if (ldp->rndprc < NBITS)
  1336.     {
  1337.       i = ldp->rw + 1;
  1338.       while (i < NI)
  1339.         {
  1340.           if (s[i])
  1341.             r |= 1;
  1342.           s[i] = 0;
  1343.           ++i;
  1344.         }
  1345.     }
  1346.   s[ldp->rw] &= ~ldp->rmsk;
  1347.   if ((r & ldp->rmbit) != 0)
  1348.     {
  1349.       if (r == ldp->rmbit)
  1350.         {
  1351.           if (lost == 0)
  1352.             {                   /* round to even */
  1353.               if ((s[ldp->re] & ldp->rebit) == 0)
  1354.                 goto mddone;
  1355.             }
  1356.           else
  1357.             {
  1358.               if (subflg != 0)
  1359.                 goto mddone;
  1360.             }
  1361.         }
  1362.       eaddm (ldp->rbit, s);
  1363.     }
  1364. mddone:
  1365. #if IBMPC
  1366.   if ((exp <= 0) && (ldp->rndprc != NBITS))
  1367. #else
  1368.   if ((exp <= 0) && (ldp->rndprc != 64) && (ldp->rndprc != NBITS))
  1369. #endif
  1370.     {
  1371.       eshup1 (s);
  1372.     }
  1373.   if (s[2] != 0)
  1374.     {                           /* overflow on roundoff */
  1375.       eshdn1 (s);
  1376.       exp += 1;
  1377.     }
  1378. mdfin:
  1379.   s[NI - 1] = 0;
  1380.   if (exp >= 32767L)
  1381.     {
  1382. #ifndef USE_INFINITY
  1383.     overf:
  1384. #endif
  1385. #ifdef USE_INFINITY
  1386.       s[1] = 32767;
  1387.       for (i = 2; i < NI - 1; i++)
  1388.         s[i] = 0;
  1389. #else
  1390.       s[1] = 32766;
  1391.       s[2] = 0;
  1392.       for (i = M + 1; i < NI - 1; i++)
  1393.         s[i] = 0xffff;
  1394.       s[NI - 1] = 0;
  1395.       if ((ldp->rndprc < 64) || (ldp->rndprc == 113))
  1396.         {
  1397.           s[ldp->rw] &= ~ldp->rmsk;
  1398.           if (ldp->rndprc == 24)
  1399.             {
  1400.               s[5] = 0;
  1401.               s[6] = 0;
  1402.             }
  1403.         }
  1404. #endif
  1405.       return;
  1406.     }
  1407.   if (exp < 0)
  1408.     s[1] = 0;
  1409.   else
  1410.     s[1] = (unsigned short) exp;
  1411. }
  1412.  
  1413.  
  1414.  
  1415. /*
  1416. ;       Subtract external format numbers.
  1417. ;
  1418. ;       unsigned short a[NE], b[NE], c[NE];
  1419. ;       LDPARMS *ldp;
  1420. ;       esub( a, b, c, ldp );    c = b - a
  1421. */
  1422.  
  1423. static void
  1424. esub (_CONST short unsigned int *a, _CONST short unsigned int *b,
  1425.       short unsigned int *c, LDPARMS * ldp)
  1426. {
  1427.  
  1428. #ifdef NANS
  1429.   if (eisnan (a))
  1430.     {
  1431.       emov (a, c);
  1432.       return;
  1433.     }
  1434.   if (eisnan (b))
  1435.     {
  1436.       emov (b, c);
  1437.       return;
  1438.     }
  1439. /* Infinity minus infinity is a NaN.
  1440.  * Test for subtracting infinities of the same sign.
  1441.  */
  1442.   if (eisinf (a) && eisinf (b) && ((eisneg (a) ^ eisneg (b)) == 0))
  1443.     {
  1444.       mtherr ("esub", DOMAIN);
  1445.       enan (c, NBITS);
  1446.       return;
  1447.     }
  1448. #endif
  1449.   eadd1 (a, b, c, 1, ldp);
  1450. }
  1451.  
  1452.  
  1453.  
  1454. static void
  1455. eadd1 (_CONST short unsigned int *a, _CONST short unsigned int *b,
  1456.        short unsigned int *c, int subflg, LDPARMS * ldp)
  1457. {
  1458.   unsigned short ai[NI], bi[NI], ci[NI];
  1459.   int i, lost, j, k;
  1460.   long lt, lta, ltb;
  1461.  
  1462. #ifdef USE_INFINITY
  1463.   if (eisinf (a))
  1464.     {
  1465.       emov (a, c);
  1466.       if (subflg)
  1467.         eneg (c);
  1468.       return;
  1469.     }
  1470.   if (eisinf (b))
  1471.     {
  1472.       emov (b, c);
  1473.       return;
  1474.     }
  1475. #endif
  1476.   emovi (a, ai);
  1477.   emovi (b, bi);
  1478.   if (subflg)
  1479.     ai[0] = ~ai[0];
  1480.  
  1481. /* compare exponents */
  1482.   lta = ai[E];
  1483.   ltb = bi[E];
  1484.   lt = lta - ltb;
  1485.   if (lt > 0L)
  1486.     {                           /* put the larger number in bi */
  1487.       emovz (bi, ci);
  1488.       emovz (ai, bi);
  1489.       emovz (ci, ai);
  1490.       ltb = bi[E];
  1491.       lt = -lt;
  1492.     }
  1493.   lost = 0;
  1494.   if (lt != 0L)
  1495.     {
  1496.       if (lt < (long) (-NBITS - 1))
  1497.         goto done;              /* answer same as larger addend */
  1498.       k = (int) lt;
  1499.       lost = eshift (ai, k);    /* shift the smaller number down */
  1500.     }
  1501.   else
  1502.     {
  1503. /* exponents were the same, so must compare significands */
  1504.       i = ecmpm (ai, bi);
  1505.       if (i == 0)
  1506.         {                       /* the numbers are identical in magnitude */
  1507.           /* if different signs, result is zero */
  1508.           if (ai[0] != bi[0])
  1509.             {
  1510.               eclear (c);
  1511.               return;
  1512.             }
  1513.           /* if same sign, result is double */
  1514.           /* double denomalized tiny number */
  1515.           if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
  1516.             {
  1517.               eshup1 (bi);
  1518.               goto done;
  1519.             }
  1520.           /* add 1 to exponent unless both are zero! */
  1521.           for (j = 1; j < NI - 1; j++)
  1522.             {
  1523.               if (bi[j] != 0)
  1524.                 {
  1525. /* This could overflow, but let emovo take care of that. */
  1526.                   ltb += 1;
  1527.                   break;
  1528.                 }
  1529.             }
  1530.           bi[E] = (unsigned short) ltb;
  1531.           goto done;
  1532.         }
  1533.       if (i > 0)
  1534.         {                       /* put the larger number in bi */
  1535.           emovz (bi, ci);
  1536.           emovz (ai, bi);
  1537.           emovz (ci, ai);
  1538.         }
  1539.     }
  1540.   if (ai[0] == bi[0])
  1541.     {
  1542.       eaddm (ai, bi);
  1543.       subflg = 0;
  1544.     }
  1545.   else
  1546.     {
  1547.       esubm (ai, bi);
  1548.       subflg = 1;
  1549.     }
  1550.   emdnorm (bi, lost, subflg, ltb, 64, ldp);
  1551.  
  1552. done:
  1553.   emovo (bi, c, ldp);
  1554. }
  1555.  
  1556.  
  1557.  
  1558. /*
  1559. ;       Divide.
  1560. ;
  1561. ;       unsigned short a[NE], b[NE], c[NE];
  1562. ;       LDPARMS *ldp;
  1563. ;       ediv( a, b, c, ldp );   c = b / a
  1564. */
  1565. static void
  1566. ediv (_CONST short unsigned int *a, _CONST short unsigned int *b,
  1567.       short unsigned int *c, LDPARMS * ldp)
  1568. {
  1569.   unsigned short ai[NI], bi[NI];
  1570.   int i;
  1571.   long lt, lta, ltb;
  1572.  
  1573. #ifdef NANS
  1574. /* Return any NaN input. */
  1575.   if (eisnan (a))
  1576.     {
  1577.       emov (a, c);
  1578.       return;
  1579.     }
  1580.   if (eisnan (b))
  1581.     {
  1582.       emov (b, c);
  1583.       return;
  1584.     }
  1585. /* Zero over zero, or infinity over infinity, is a NaN. */
  1586.   if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
  1587.       || (eisinf (a) && eisinf (b)))
  1588.     {
  1589.       mtherr ("ediv", DOMAIN);
  1590.       enan (c, NBITS);
  1591.       return;
  1592.     }
  1593. #endif
  1594. /* Infinity over anything else is infinity. */
  1595. #ifdef USE_INFINITY
  1596.   if (eisinf (b))
  1597.     {
  1598.       if (eisneg (a) ^ eisneg (b))
  1599.         *(c + (NE - 1)) = 0x8000;
  1600.       else
  1601.         *(c + (NE - 1)) = 0;
  1602.       einfin (c, ldp);
  1603.       return;
  1604.     }
  1605.   if (eisinf (a))
  1606.     {
  1607.       eclear (c);
  1608.       return;
  1609.     }
  1610. #endif
  1611.   emovi (a, ai);
  1612.   emovi (b, bi);
  1613.   lta = ai[E];
  1614.   ltb = bi[E];
  1615.   if (bi[E] == 0)
  1616.     {                           /* See if numerator is zero. */
  1617.       for (i = 1; i < NI - 1; i++)
  1618.         {
  1619.           if (bi[i] != 0)
  1620.             {
  1621.               ltb -= enormlz (bi);
  1622.               goto dnzro1;
  1623.             }
  1624.         }
  1625.       eclear (c);
  1626.       return;
  1627.     }
  1628. dnzro1:
  1629.  
  1630.   if (ai[E] == 0)
  1631.     {                           /* possible divide by zero */
  1632.       for (i = 1; i < NI - 1; i++)
  1633.         {
  1634.           if (ai[i] != 0)
  1635.             {
  1636.               lta -= enormlz (ai);
  1637.               goto dnzro2;
  1638.             }
  1639.         }
  1640.       if (ai[0] == bi[0])
  1641.         *(c + (NE - 1)) = 0;
  1642.       else
  1643.         *(c + (NE - 1)) = 0x8000;
  1644.       einfin (c, ldp);
  1645.       mtherr ("ediv", SING);
  1646.       return;
  1647.     }
  1648. dnzro2:
  1649.  
  1650.   i = edivm (ai, bi, ldp);
  1651. /* calculate exponent */
  1652.   lt = ltb - lta + EXONE;
  1653.   emdnorm (bi, i, 0, lt, 64, ldp);
  1654. /* set the sign */
  1655.   if (ai[0] == bi[0])
  1656.     bi[0] = 0;
  1657.   else
  1658.     bi[0] = 0Xffff;
  1659.   emovo (bi, c, ldp);
  1660. }
  1661.  
  1662.  
  1663.  
  1664. /*
  1665. ;       Multiply.
  1666. ;
  1667. ;       unsigned short a[NE], b[NE], c[NE];
  1668. ;       LDPARMS *ldp
  1669. ;       emul( a, b, c, ldp );   c = b * a
  1670. */
  1671. static void
  1672. emul (_CONST short unsigned int *a, _CONST short unsigned int *b,
  1673.       short unsigned int *c, LDPARMS * ldp)
  1674. {
  1675.   unsigned short ai[NI], bi[NI];
  1676.   int i, j;
  1677.   long lt, lta, ltb;
  1678.  
  1679. #ifdef NANS
  1680. /* NaN times anything is the same NaN. */
  1681.   if (eisnan (a))
  1682.     {
  1683.       emov (a, c);
  1684.       return;
  1685.     }
  1686.   if (eisnan (b))
  1687.     {
  1688.       emov (b, c);
  1689.       return;
  1690.     }
  1691. /* Zero times infinity is a NaN. */
  1692.   if ((eisinf (a) && (ecmp (b, ezero) == 0))
  1693.       || (eisinf (b) && (ecmp (a, ezero) == 0)))
  1694.     {
  1695.       mtherr ("emul", DOMAIN);
  1696.       enan (c, NBITS);
  1697.       return;
  1698.     }
  1699. #endif
  1700. /* Infinity times anything else is infinity. */
  1701. #ifdef USE_INFINITY
  1702.   if (eisinf (a) || eisinf (b))
  1703.     {
  1704.       if (eisneg (a) ^ eisneg (b))
  1705.         *(c + (NE - 1)) = 0x8000;
  1706.       else
  1707.         *(c + (NE - 1)) = 0;
  1708.       einfin (c, ldp);
  1709.       return;
  1710.     }
  1711. #endif
  1712.   emovi (a, ai);
  1713.   emovi (b, bi);
  1714.   lta = ai[E];
  1715.   ltb = bi[E];
  1716.   if (ai[E] == 0)
  1717.     {
  1718.       for (i = 1; i < NI - 1; i++)
  1719.         {
  1720.           if (ai[i] != 0)
  1721.             {
  1722.               lta -= enormlz (ai);
  1723.               goto mnzer1;
  1724.             }
  1725.         }
  1726.       eclear (c);
  1727.       return;
  1728.     }
  1729. mnzer1:
  1730.  
  1731.   if (bi[E] == 0)
  1732.     {
  1733.       for (i = 1; i < NI - 1; i++)
  1734.         {
  1735.           if (bi[i] != 0)
  1736.             {
  1737.               ltb -= enormlz (bi);
  1738.               goto mnzer2;
  1739.             }
  1740.         }
  1741.       eclear (c);
  1742.       return;
  1743.     }
  1744. mnzer2:
  1745.  
  1746. /* Multiply significands */
  1747.   j = emulm (ai, bi, ldp);
  1748. /* calculate exponent */
  1749.   lt = lta + ltb - (EXONE - 1);
  1750.   emdnorm (bi, j, 0, lt, 64, ldp);
  1751. /* calculate sign of product */
  1752.   if (ai[0] == bi[0])
  1753.     bi[0] = 0;
  1754.   else
  1755.     bi[0] = 0xffff;
  1756.   emovo (bi, c, ldp);
  1757. }
  1758.  
  1759.  
  1760.  
  1761. #if LDBL_MANT_DIG > 64
  1762. static void
  1763. e113toe (short unsigned int *pe, short unsigned int *y, LDPARMS * ldp)
  1764. {
  1765.   register unsigned short r;
  1766.   unsigned short *e, *p;
  1767.   unsigned short yy[NI];
  1768.   int denorm, i;
  1769.  
  1770.   e = pe;
  1771.   denorm = 0;
  1772.   ecleaz (yy);
  1773. #ifdef IBMPC
  1774.   e += 7;
  1775. #endif
  1776.   r = *e;
  1777.   yy[0] = 0;
  1778.   if (r & 0x8000)
  1779.     yy[0] = 0xffff;
  1780.   r &= 0x7fff;
  1781. #ifdef USE_INFINITY
  1782.   if (r == 0x7fff)
  1783.     {
  1784. #ifdef NANS
  1785. #ifdef IBMPC
  1786.       for (i = 0; i < 7; i++)
  1787.         {
  1788.           if (pe[i] != 0)
  1789.             {
  1790.               enan (y, NBITS);
  1791.               return;
  1792.             }
  1793.         }
  1794. #else /* !IBMPC */
  1795.       for (i = 1; i < 8; i++)
  1796.         {
  1797.           if (pe[i] != 0)
  1798.             {
  1799.               enan (y, NBITS);
  1800.               return;
  1801.             }
  1802.         }
  1803. #endif /* !IBMPC */
  1804. #endif /* NANS */
  1805.       eclear (y);
  1806.       einfin (y, ldp);
  1807.       if (*e & 0x8000)
  1808.         eneg (y);
  1809.       return;
  1810.     }
  1811. #endif /* INFINITY */
  1812.   yy[E] = r;
  1813.   p = &yy[M + 1];
  1814. #ifdef IBMPC
  1815.   for (i = 0; i < 7; i++)
  1816.     *p++ = *(--e);
  1817. #else /* IBMPC */
  1818.   ++e;
  1819.   for (i = 0; i < 7; i++)
  1820.     *p++ = *e++;
  1821. #endif /* IBMPC */
  1822. /* If denormal, remove the implied bit; else shift down 1. */
  1823.   if (r == 0)
  1824.     {
  1825.       yy[M] = 0;
  1826.     }
  1827.   else
  1828.     {
  1829.       yy[M] = 1;
  1830.       eshift (yy, -1);
  1831.     }
  1832.   emovo (yy, y, ldp);
  1833. }
  1834.  
  1835. /* move out internal format to ieee long double */
  1836. static void
  1837. toe113 (short unsigned int *a, short unsigned int *b)
  1838. {
  1839.   register unsigned short *p, *q;
  1840.   unsigned short i;
  1841.  
  1842. #ifdef NANS
  1843.   if (eiisnan (a))
  1844.     {
  1845.       enan (b, 113);
  1846.       return;
  1847.     }
  1848. #endif
  1849.   p = a;
  1850. #ifdef MIEEE
  1851.   q = b;
  1852. #else
  1853.   q = b + 7;                    /* point to output exponent */
  1854. #endif
  1855.  
  1856. /* If not denormal, delete the implied bit. */
  1857.   if (a[E] != 0)
  1858.     {
  1859.       eshup1 (a);
  1860.     }
  1861. /* combine sign and exponent */
  1862.   i = *p++;
  1863. #ifdef MIEEE
  1864.   if (i)
  1865.     *q++ = *p++ | 0x8000;
  1866.   else
  1867.     *q++ = *p++;
  1868. #else
  1869.   if (i)
  1870.     *q-- = *p++ | 0x8000;
  1871.   else
  1872.     *q-- = *p++;
  1873. #endif
  1874. /* skip over guard word */
  1875.   ++p;
  1876. /* move the significand */
  1877. #ifdef MIEEE
  1878.   for (i = 0; i < 7; i++)
  1879.     *q++ = *p++;
  1880. #else
  1881.   for (i = 0; i < 7; i++)
  1882.     *q-- = *p++;
  1883. #endif
  1884. }
  1885. #endif /* LDBL_MANT_DIG > 64 */
  1886.  
  1887.  
  1888. #if LDBL_MANT_DIG == 64
  1889. static void
  1890. e64toe (short unsigned int *pe, short unsigned int *y, LDPARMS * ldp)
  1891. {
  1892.   unsigned short yy[NI];
  1893.   unsigned short *p, *q, *e;
  1894.   int i;
  1895.  
  1896.   e = pe;
  1897.   p = yy;
  1898.  
  1899.   for (i = 0; i < NE - 5; i++)
  1900.     *p++ = 0;
  1901. #ifdef IBMPC
  1902.   for (i = 0; i < 5; i++)
  1903.     *p++ = *e++;
  1904. #endif
  1905. #ifdef DEC
  1906.   for (i = 0; i < 5; i++)
  1907.     *p++ = *e++;
  1908. #endif
  1909. #ifdef MIEEE
  1910.   p = &yy[0] + (NE - 1);
  1911.   *p-- = *e++;
  1912.   ++e;                          /* MIEEE skips over 2nd short */
  1913.   for (i = 0; i < 4; i++)
  1914.     *p-- = *e++;
  1915. #endif
  1916.  
  1917. #ifdef IBMPC
  1918. /* For Intel long double, shift denormal significand up 1
  1919.    -- but only if the top significand bit is zero.  */
  1920.   if ((yy[NE - 1] & 0x7fff) == 0 && (yy[NE - 2] & 0x8000) == 0)
  1921.     {
  1922.       unsigned short temp[NI + 1];
  1923.       emovi (yy, temp);
  1924.       eshup1 (temp);
  1925.       emovo (temp, y, ldp);
  1926.       return;
  1927.     }
  1928. #endif
  1929. #ifdef USE_INFINITY
  1930. /* Point to the exponent field.  */
  1931.   p = &yy[NE - 1];
  1932.   if ((*p & 0x7fff) == 0x7fff)
  1933.     {
  1934. #ifdef NANS
  1935. #ifdef IBMPC
  1936.       for (i = 0; i < 4; i++)
  1937.         {
  1938.           if ((i != 3 && pe[i] != 0)
  1939.               /* Check for Intel long double infinity pattern.  */
  1940.               || (i == 3 && pe[i] != 0x8000))
  1941.             {
  1942.               enan (y, NBITS);
  1943.               return;
  1944.             }
  1945.         }
  1946. #endif
  1947. #ifdef MIEEE
  1948.       for (i = 2; i <= 5; i++)
  1949.         {
  1950.           if (pe[i] != 0)
  1951.             {
  1952.               enan (y, NBITS);
  1953.               return;
  1954.             }
  1955.         }
  1956. #endif
  1957. #endif /* NANS */
  1958.       eclear (y);
  1959.       einfin (y, ldp);
  1960.       if (*p & 0x8000)
  1961.         eneg (y);
  1962.       return;
  1963.     }
  1964. #endif /* USE_INFINITY */
  1965.   p = yy;
  1966.   q = y;
  1967.   for (i = 0; i < NE; i++)
  1968.     *q++ = *p++;
  1969. }
  1970.  
  1971. /* move out internal format to ieee long double */
  1972. static void
  1973. toe64 (short unsigned int *a, short unsigned int *b)
  1974. {
  1975.   register unsigned short *p, *q;
  1976.   unsigned short i;
  1977.  
  1978. #ifdef NANS
  1979.   if (eiisnan (a))
  1980.     {
  1981.       enan (b, 64);
  1982.       return;
  1983.     }
  1984. #endif
  1985. #ifdef IBMPC
  1986. /* Shift Intel denormal significand down 1.  */
  1987.   if (a[E] == 0)
  1988.     eshdn1 (a);
  1989. #endif
  1990.   p = a;
  1991. #ifdef MIEEE
  1992.   q = b;
  1993. #else
  1994.   q = b + 4;                    /* point to output exponent */
  1995. /* NOTE: Intel data type is 96 bits wide, clear the last word here. */
  1996.   *(q + 1) = 0;
  1997. #endif
  1998.  
  1999. /* combine sign and exponent */
  2000.   i = *p++;
  2001. #ifdef MIEEE
  2002.   if (i)
  2003.     *q++ = *p++ | 0x8000;
  2004.   else
  2005.     *q++ = *p++;
  2006.   *q++ = 0;                     /* leave 2nd short blank */
  2007. #else
  2008.   if (i)
  2009.     *q-- = *p++ | 0x8000;
  2010.   else
  2011.     *q-- = *p++;
  2012. #endif
  2013. /* skip over guard word */
  2014.   ++p;
  2015. /* move the significand */
  2016. #ifdef MIEEE
  2017.   for (i = 0; i < 4; i++)
  2018.     *q++ = *p++;
  2019. #else
  2020. #ifdef USE_INFINITY
  2021. #ifdef IBMPC
  2022.   if (eiisinf (a))
  2023.     {
  2024.       /* Intel long double infinity.  */
  2025.       *q-- = 0x8000;
  2026.       *q-- = 0;
  2027.       *q-- = 0;
  2028.       *q = 0;
  2029.       return;
  2030.     }
  2031. #endif /* IBMPC */
  2032. #endif /* USE_INFINITY */
  2033.   for (i = 0; i < 4; i++)
  2034.     *q-- = *p++;
  2035. #endif
  2036. }
  2037.  
  2038. #endif /* LDBL_MANT_DIG == 64 */
  2039.  
  2040. #if LDBL_MANT_DIG == 53
  2041. /*
  2042. ; Convert IEEE double precision to e type
  2043. ;       double d;
  2044. ;       unsigned short x[N+2];
  2045. ;       e53toe( &d, x );
  2046. */
  2047. static void
  2048. e53toe (short unsigned int *pe, short unsigned int *y, LDPARMS * ldp)
  2049. {
  2050. #ifdef DEC
  2051.  
  2052.   dectoe (pe, y);               /* see etodec.c */
  2053.  
  2054. #else
  2055.  
  2056.   register unsigned short r;
  2057.   register unsigned short *p, *e;
  2058.   unsigned short yy[NI];
  2059.   int denorm, k;
  2060.  
  2061.   e = pe;
  2062.   denorm = 0;                   /* flag if denormalized number */
  2063.   ecleaz (yy);
  2064. #ifdef IBMPC
  2065.   e += 3;
  2066. #endif
  2067. #ifdef DEC
  2068.   e += 3;
  2069. #endif
  2070.   r = *e;
  2071.   yy[0] = 0;
  2072.   if (r & 0x8000)
  2073.     yy[0] = 0xffff;
  2074.   yy[M] = (r & 0x0f) | 0x10;
  2075.   r &= ~0x800f;                 /* strip sign and 4 significand bits */
  2076. #ifdef USE_INFINITY
  2077.   if (r == 0x7ff0)
  2078.     {
  2079. #ifdef NANS
  2080. #ifdef IBMPC
  2081.       if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
  2082.           || (pe[1] != 0) || (pe[0] != 0))
  2083.         {
  2084.           enan (y, NBITS);
  2085.           return;
  2086.         }
  2087. #else /* !IBMPC */
  2088.       if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
  2089.           || (pe[2] != 0) || (pe[3] != 0))
  2090.         {
  2091.           enan (y, NBITS);
  2092.           return;
  2093.         }
  2094. #endif /* !IBMPC */
  2095. #endif /* NANS */
  2096.       eclear (y);
  2097.       einfin (y, ldp);
  2098.       if (yy[0])
  2099.         eneg (y);
  2100.       return;
  2101.     }
  2102. #endif
  2103.   r >>= 4;
  2104. /* If zero exponent, then the significand is denormalized.
  2105.  * So, take back the understood high significand bit. */
  2106.   if (r == 0)
  2107.     {
  2108.       denorm = 1;
  2109.       yy[M] &= ~0x10;
  2110.     }
  2111.   r += EXONE - 01777;
  2112.   yy[E] = r;
  2113.   p = &yy[M + 1];
  2114. #ifdef IBMPC
  2115.   *p++ = *(--e);
  2116.   *p++ = *(--e);
  2117.   *p++ = *(--e);
  2118. #else /* !IBMPC */
  2119.   ++e;
  2120.   *p++ = *e++;
  2121.   *p++ = *e++;
  2122.   *p++ = *e++;
  2123. #endif /* !IBMPC */
  2124.   (void) eshift (yy, -5);
  2125.   if (denorm)
  2126.     {                           /* if zero exponent, then normalize the significand */
  2127.       if ((k = enormlz (yy)) > NBITS)
  2128.         ecleazs (yy);
  2129.       else
  2130.         yy[E] -= (unsigned short) (k - 1);
  2131.     }
  2132.   emovo (yy, y, ldp);
  2133. #endif /* !DEC */
  2134. }
  2135.  
  2136. /*
  2137. ; e type to IEEE double precision
  2138. ;       double d;
  2139. ;       unsigned short x[NE];
  2140. ;       etoe53( x, &d );
  2141. */
  2142.  
  2143. #ifdef DEC
  2144.  
  2145. static void
  2146. etoe53 (x, e)
  2147.      unsigned short *x, *e;
  2148. {
  2149.   etodec (x, e);                /* see etodec.c */
  2150. }
  2151.  
  2152. static void
  2153. toe53 (x, y)
  2154.      unsigned short *x, *y;
  2155. {
  2156.   todec (x, y);
  2157. }
  2158.  
  2159. #else
  2160.  
  2161. static void
  2162. toe53 (short unsigned int *x, short unsigned int *y)
  2163. {
  2164.   unsigned short i;
  2165.   unsigned short *p;
  2166.  
  2167.  
  2168. #ifdef NANS
  2169.   if (eiisnan (x))
  2170.     {
  2171.       enan (y, 53);
  2172.       return;
  2173.     }
  2174. #endif
  2175.   p = &x[0];
  2176. #ifdef IBMPC
  2177.   y += 3;
  2178. #endif
  2179. #ifdef DEC
  2180.   y += 3;
  2181. #endif
  2182.   *y = 0;                       /* output high order */
  2183.   if (*p++)
  2184.     *y = 0x8000;                /* output sign bit */
  2185.  
  2186.   i = *p++;
  2187.   if (i >= (unsigned int) 2047)
  2188.     {                           /* Saturate at largest number less than infinity. */
  2189. #ifdef USE_INFINITY
  2190.       *y |= 0x7ff0;
  2191. #ifdef IBMPC
  2192.       *(--y) = 0;
  2193.       *(--y) = 0;
  2194.       *(--y) = 0;
  2195. #else /* !IBMPC */
  2196.       ++y;
  2197.       *y++ = 0;
  2198.       *y++ = 0;
  2199.       *y++ = 0;
  2200. #endif /* IBMPC */
  2201. #else /* !USE_INFINITY */
  2202.       *y |= (unsigned short) 0x7fef;
  2203. #ifdef IBMPC
  2204.       *(--y) = 0xffff;
  2205.       *(--y) = 0xffff;
  2206.       *(--y) = 0xffff;
  2207. #else /* !IBMPC */
  2208.       ++y;
  2209.       *y++ = 0xffff;
  2210.       *y++ = 0xffff;
  2211.       *y++ = 0xffff;
  2212. #endif
  2213. #endif /* !USE_INFINITY */
  2214.       return;
  2215.     }
  2216.   if (i == 0)
  2217.     {
  2218.       (void) eshift (x, 4);
  2219.     }
  2220.   else
  2221.     {
  2222.       i <<= 4;
  2223.       (void) eshift (x, 5);
  2224.     }
  2225.   i |= *p++ & (unsigned short) 0x0f;    /* *p = xi[M] */
  2226.   *y |= (unsigned short) i;     /* high order output already has sign bit set */
  2227. #ifdef IBMPC
  2228.   *(--y) = *p++;
  2229.   *(--y) = *p++;
  2230.   *(--y) = *p;
  2231. #else /* !IBMPC */
  2232.   ++y;
  2233.   *y++ = *p++;
  2234.   *y++ = *p++;
  2235.   *y++ = *p++;
  2236. #endif /* !IBMPC */
  2237. }
  2238.  
  2239. #endif /* not DEC */
  2240. #endif /* LDBL_MANT_DIG == 53 */
  2241.  
  2242. #if LDBL_MANT_DIG == 24
  2243. /*
  2244. ; Convert IEEE single precision to e type
  2245. ;       float d;
  2246. ;       unsigned short x[N+2];
  2247. ;       dtox( &d, x );
  2248. */
  2249. void
  2250. e24toe (short unsigned int *pe, short unsigned int *y, LDPARMS * ldp)
  2251. {
  2252.   register unsigned short r;
  2253.   register unsigned short *p, *e;
  2254.   unsigned short yy[NI];
  2255.   int denorm, k;
  2256.  
  2257.   e = pe;
  2258.   denorm = 0;                   /* flag if denormalized number */
  2259.   ecleaz (yy);
  2260. #ifdef IBMPC
  2261.   e += 1;
  2262. #endif
  2263. #ifdef DEC
  2264.   e += 1;
  2265. #endif
  2266.   r = *e;
  2267.   yy[0] = 0;
  2268.   if (r & 0x8000)
  2269.     yy[0] = 0xffff;
  2270.   yy[M] = (r & 0x7f) | 0200;
  2271.   r &= ~0x807f;                 /* strip sign and 7 significand bits */
  2272. #ifdef USE_INFINITY
  2273.   if (r == 0x7f80)
  2274.     {
  2275. #ifdef NANS
  2276. #ifdef MIEEE
  2277.       if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
  2278.         {
  2279.           enan (y, NBITS);
  2280.           return;
  2281.         }
  2282. #else /* !MIEEE */
  2283.       if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
  2284.         {
  2285.           enan (y, NBITS);
  2286.           return;
  2287.         }
  2288. #endif /* !MIEEE */
  2289. #endif /* NANS */
  2290.       eclear (y);
  2291.       einfin (y, ldp);
  2292.       if (yy[0])
  2293.         eneg (y);
  2294.       return;
  2295.     }
  2296. #endif
  2297.   r >>= 7;
  2298. /* If zero exponent, then the significand is denormalized.
  2299.  * So, take back the understood high significand bit. */
  2300.   if (r == 0)
  2301.     {
  2302.       denorm = 1;
  2303.       yy[M] &= ~0200;
  2304.     }
  2305.   r += EXONE - 0177;
  2306.   yy[E] = r;
  2307.   p = &yy[M + 1];
  2308. #ifdef IBMPC
  2309.   *p++ = *(--e);
  2310. #endif
  2311. #ifdef DEC
  2312.   *p++ = *(--e);
  2313. #endif
  2314. #ifdef MIEEE
  2315.   ++e;
  2316.   *p++ = *e++;
  2317. #endif
  2318.   (void) eshift (yy, -8);
  2319.   if (denorm)
  2320.     {                           /* if zero exponent, then normalize the significand */
  2321.       if ((k = enormlz (yy)) > NBITS)
  2322.         ecleazs (yy);
  2323.       else
  2324.         yy[E] -= (unsigned short) (k - 1);
  2325.     }
  2326.   emovo (yy, y, ldp);
  2327. }
  2328.  
  2329. static void
  2330. toe24 (short unsigned int *x, short unsigned int *y)
  2331. {
  2332.   unsigned short i;
  2333.   unsigned short *p;
  2334.  
  2335. #ifdef NANS
  2336.   if (eiisnan (x))
  2337.     {
  2338.       enan (y, 24);
  2339.       return;
  2340.     }
  2341. #endif
  2342.   p = &x[0];
  2343. #ifdef IBMPC
  2344.   y += 1;
  2345. #endif
  2346. #ifdef DEC
  2347.   y += 1;
  2348. #endif
  2349.   *y = 0;                       /* output high order */
  2350.   if (*p++)
  2351.     *y = 0x8000;                /* output sign bit */
  2352.  
  2353.   i = *p++;
  2354.   if (i >= 255)
  2355.     {                           /* Saturate at largest number less than infinity. */
  2356. #ifdef USE_INFINITY
  2357.       *y |= (unsigned short) 0x7f80;
  2358. #ifdef IBMPC
  2359.       *(--y) = 0;
  2360. #endif
  2361. #ifdef DEC
  2362.       *(--y) = 0;
  2363. #endif
  2364. #ifdef MIEEE
  2365.       ++y;
  2366.       *y = 0;
  2367. #endif
  2368. #else /* !USE_INFINITY */
  2369.       *y |= (unsigned short) 0x7f7f;
  2370. #ifdef IBMPC
  2371.       *(--y) = 0xffff;
  2372. #endif
  2373. #ifdef DEC
  2374.       *(--y) = 0xffff;
  2375. #endif
  2376. #ifdef MIEEE
  2377.       ++y;
  2378.       *y = 0xffff;
  2379. #endif
  2380. #endif /* !USE_INFINITY */
  2381.       return;
  2382.     }
  2383.   if (i == 0)
  2384.     {
  2385.       (void) eshift (x, 7);
  2386.     }
  2387.   else
  2388.     {
  2389.       i <<= 7;
  2390.       (void) eshift (x, 8);
  2391.     }
  2392.   i |= *p++ & (unsigned short) 0x7f;    /* *p = xi[M] */
  2393.   *y |= i;                      /* high order output already has sign bit set */
  2394. #ifdef IBMPC
  2395.   *(--y) = *p;
  2396. #endif
  2397. #ifdef DEC
  2398.   *(--y) = *p;
  2399. #endif
  2400. #ifdef MIEEE
  2401.   ++y;
  2402.   *y = *p;
  2403. #endif
  2404. }
  2405. #endif /* LDBL_MANT_DIG == 24 */
  2406.  
  2407. /* Compare two e type numbers.
  2408.  *
  2409.  * unsigned short a[NE], b[NE];
  2410.  * ecmp( a, b );
  2411.  *
  2412.  *  returns +1 if a > b
  2413.  *           0 if a == b
  2414.  *          -1 if a < b
  2415.  *          -2 if either a or b is a NaN.
  2416.  */
  2417. static int
  2418. ecmp (_CONST short unsigned int *a, _CONST short unsigned int *b)
  2419. {
  2420.   unsigned short ai[NI], bi[NI];
  2421.   register unsigned short *p, *q;
  2422.   register int i;
  2423.   int msign;
  2424.  
  2425. #ifdef NANS
  2426.   if (eisnan (a) || eisnan (b))
  2427.     return (-2);
  2428. #endif
  2429.   emovi (a, ai);
  2430.   p = ai;
  2431.   emovi (b, bi);
  2432.   q = bi;
  2433.  
  2434.   if (*p != *q)
  2435.     {                           /* the signs are different */
  2436. /* -0 equals + 0 */
  2437.       for (i = 1; i < NI - 1; i++)
  2438.         {
  2439.           if (ai[i] != 0)
  2440.             goto nzro;
  2441.           if (bi[i] != 0)
  2442.             goto nzro;
  2443.         }
  2444.       return (0);
  2445.     nzro:
  2446.       if (*p == 0)
  2447.         return (1);
  2448.       else
  2449.         return (-1);
  2450.     }
  2451. /* both are the same sign */
  2452.   if (*p == 0)
  2453.     msign = 1;
  2454.   else
  2455.     msign = -1;
  2456.   i = NI - 1;
  2457.   do
  2458.     {
  2459.       if (*p++ != *q++)
  2460.         {
  2461.           goto diff;
  2462.         }
  2463.     }
  2464.   while (--i > 0);
  2465.  
  2466.   return (0);                   /* equality */
  2467.  
  2468.  
  2469.  
  2470. diff:
  2471.  
  2472.   if (*(--p) > *(--q))
  2473.     return (msign);             /* p is bigger */
  2474.   else
  2475.     return (-msign);            /* p is littler */
  2476. }
  2477.  
  2478.  
  2479. /*
  2480. ;       Shift significand
  2481. ;
  2482. ;       Shifts significand area up or down by the number of bits
  2483. ;       given by the variable sc.
  2484. */
  2485. static int
  2486. eshift (short unsigned int *x, int sc)
  2487. {
  2488.   unsigned short lost;
  2489.   unsigned short *p;
  2490.  
  2491.   if (sc == 0)
  2492.     return (0);
  2493.  
  2494.   lost = 0;
  2495.   p = x + NI - 1;
  2496.  
  2497.   if (sc < 0)
  2498.     {
  2499.       sc = -sc;
  2500.       while (sc >= 16)
  2501.         {
  2502.           lost |= *p;           /* remember lost bits */
  2503.           eshdn6 (x);
  2504.           sc -= 16;
  2505.         }
  2506.  
  2507.       while (sc >= 8)
  2508.         {
  2509.           lost |= *p & 0xff;
  2510.           eshdn8 (x);
  2511.           sc -= 8;
  2512.         }
  2513.  
  2514.       while (sc > 0)
  2515.         {
  2516.           lost |= *p & 1;
  2517.           eshdn1 (x);
  2518.           sc -= 1;
  2519.         }
  2520.     }
  2521.   else
  2522.     {
  2523.       while (sc >= 16)
  2524.         {
  2525.           eshup6 (x);
  2526.           sc -= 16;
  2527.         }
  2528.  
  2529.       while (sc >= 8)
  2530.         {
  2531.           eshup8 (x);
  2532.           sc -= 8;
  2533.         }
  2534.  
  2535.       while (sc > 0)
  2536.         {
  2537.           eshup1 (x);
  2538.           sc -= 1;
  2539.         }
  2540.     }
  2541.   if (lost)
  2542.     lost = 1;
  2543.   return ((int) lost);
  2544. }
  2545.  
  2546.  
  2547.  
  2548. /*
  2549. ;       normalize
  2550. ;
  2551. ; Shift normalizes the significand area pointed to by argument
  2552. ; shift count (up = positive) is returned.
  2553. */
  2554. static int
  2555. enormlz (short unsigned int *x)
  2556. {
  2557.   register unsigned short *p;
  2558.   int sc;
  2559.  
  2560.   sc = 0;
  2561.   p = &x[M];
  2562.   if (*p != 0)
  2563.     goto normdn;
  2564.   ++p;
  2565.   if (*p & 0x8000)
  2566.     return (0);                 /* already normalized */
  2567.   while (*p == 0)
  2568.     {
  2569.       eshup6 (x);
  2570.       sc += 16;
  2571. /* With guard word, there are NBITS+16 bits available.
  2572.  * return true if all are zero.
  2573.  */
  2574.       if (sc > NBITS)
  2575.         return (sc);
  2576.     }
  2577. /* see if high byte is zero */
  2578.   while ((*p & 0xff00) == 0)
  2579.     {
  2580.       eshup8 (x);
  2581.       sc += 8;
  2582.     }
  2583. /* now shift 1 bit at a time */
  2584.   while ((*p & 0x8000) == 0)
  2585.     {
  2586.       eshup1 (x);
  2587.       sc += 1;
  2588.       if (sc > (NBITS + 16))
  2589.         {
  2590.           mtherr ("enormlz", UNDERFLOW);
  2591.           return (sc);
  2592.         }
  2593.     }
  2594.   return (sc);
  2595.  
  2596. /* Normalize by shifting down out of the high guard word
  2597.    of the significand */
  2598. normdn:
  2599.  
  2600.   if (*p & 0xff00)
  2601.     {
  2602.       eshdn8 (x);
  2603.       sc -= 8;
  2604.     }
  2605.   while (*p != 0)
  2606.     {
  2607.       eshdn1 (x);
  2608.       sc -= 1;
  2609.  
  2610.       if (sc < -NBITS)
  2611.         {
  2612.           mtherr ("enormlz", OVERFLOW);
  2613.           return (sc);
  2614.         }
  2615.     }
  2616.   return (sc);
  2617. }
  2618.  
  2619.  
  2620.  
  2621.  
  2622. /* Convert e type number to decimal format ASCII string.
  2623.  * The constants are for 64 bit precision.
  2624.  */
  2625.  
  2626. #define NTEN 12
  2627. #define MAXP 4096
  2628.  
  2629. #if NE == 10
  2630. static _CONST unsigned short etens[NTEN + 1][NE] = {
  2631.   {0x6576, 0x4a92, 0x804a, 0x153f,
  2632.    0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
  2633.   {0x6a32, 0xce52, 0x329a, 0x28ce,
  2634.    0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
  2635.   {0x526c, 0x50ce, 0xf18b, 0x3d28,
  2636.    0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
  2637.   {0x9c66, 0x58f8, 0xbc50, 0x5c54,
  2638.    0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
  2639.   {0x851e, 0xeab7, 0x98fe, 0x901b,
  2640.    0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
  2641.   {0x0235, 0x0137, 0x36b1, 0x336c,
  2642.    0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
  2643.   {0x50f8, 0x25fb, 0xc76b, 0x6b71,
  2644.    0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
  2645.   {0x0000, 0x0000, 0x0000, 0x0000,
  2646.    0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
  2647.   {0x0000, 0x0000, 0x0000, 0x0000,
  2648.    0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
  2649.   {0x0000, 0x0000, 0x0000, 0x0000,
  2650.    0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
  2651.   {0x0000, 0x0000, 0x0000, 0x0000,
  2652.    0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
  2653.   {0x0000, 0x0000, 0x0000, 0x0000,
  2654.    0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
  2655.   {0x0000, 0x0000, 0x0000, 0x0000,
  2656.    0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
  2657. };
  2658.  
  2659. static _CONST unsigned short emtens[NTEN + 1][NE] = {
  2660.   {0x2030, 0xcffc, 0xa1c3, 0x8123,
  2661.    0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
  2662.   {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
  2663.    0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
  2664.   {0xf53f, 0xf698, 0x6bd3, 0x0158,
  2665.    0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
  2666.   {0xe731, 0x04d4, 0xe3f2, 0xd332,
  2667.    0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
  2668.   {0xa23e, 0x5308, 0xfefb, 0x1155,
  2669.    0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
  2670.   {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
  2671.    0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
  2672.   {0x2a20, 0x6224, 0x47b3, 0x98d7,
  2673.    0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
  2674.   {0x0b5b, 0x4af2, 0xa581, 0x18ed,
  2675.    0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
  2676.   {0xbf71, 0xa9b3, 0x7989, 0xbe68,
  2677.    0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
  2678.   {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
  2679.    0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
  2680.   {0xc155, 0xa4a8, 0x404e, 0x6113,
  2681.    0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
  2682.   {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
  2683.    0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
  2684.   {0xcccd, 0xcccc, 0xcccc, 0xcccc,
  2685.    0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
  2686. };
  2687. #else
  2688. static _CONST unsigned short etens[NTEN + 1][NE] = {
  2689.   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},    /* 10**4096 */
  2690.   {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
  2691.   {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
  2692.   {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
  2693.   {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
  2694.   {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
  2695.   {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
  2696.   {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
  2697.   {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
  2698.   {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
  2699.   {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
  2700.   {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
  2701.   {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
  2702. };
  2703.  
  2704. static _CONST unsigned short emtens[NTEN + 1][NE] = {
  2705.   {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
  2706.   {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
  2707.   {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
  2708.   {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
  2709.   {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
  2710.   {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
  2711.   {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
  2712.   {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
  2713.   {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
  2714.   {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
  2715.   {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
  2716.   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
  2717.   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
  2718. };
  2719. #endif
  2720.  
  2721.  
  2722.  
  2723. /* ASCII string outputs for unix */
  2724.  
  2725.  
  2726. #if 0
  2727. void
  2728. _IO_ldtostr (x, string, ndigs, flags, fmt)
  2729.      long double *x;
  2730.      char *string;
  2731.      int ndigs;
  2732.      int flags;
  2733.      char fmt;
  2734. {
  2735.   unsigned short w[NI];
  2736.   char *t, *u;
  2737.   LDPARMS rnd;
  2738.   LDPARMS *ldp = &rnd;
  2739.  
  2740.   rnd.rlast = -1;
  2741.   rnd.rndprc = NBITS;
  2742.  
  2743.   if (sizeof (long double) == 16)
  2744.     e113toe ((unsigned short *) x, w, ldp);
  2745.   else
  2746.     e64toe ((unsigned short *) x, w, ldp);
  2747.  
  2748.   etoasc (w, string, ndigs, -1, ldp);
  2749.   if (ndigs == 0 && flags == 0)
  2750.     {
  2751.       /* Delete the decimal point unless alternate format.  */
  2752.       t = string;
  2753.       while (*t != '.')
  2754.         ++t;
  2755.       u = t + 1;
  2756.       while (*t != '\0')
  2757.         *t++ = *u++;
  2758.     }
  2759.   if (*string == ' ')
  2760.     {
  2761.       t = string;
  2762.       u = t + 1;
  2763.       while (*t != '\0')
  2764.         *t++ = *u++;
  2765.     }
  2766.   if (fmt == 'E')
  2767.     {
  2768.       t = string;
  2769.       while (*t != 'e')
  2770.         ++t;
  2771.       *t = 'E';
  2772.     }
  2773. }
  2774.  
  2775. #endif
  2776.  
  2777. /* This routine will not return more than NDEC+1 digits. */
  2778.  
  2779. char *
  2780. _ldtoa_r (struct _reent *ptr, long double d, int mode, int ndigits,
  2781.           int *decpt, int *sign, char **rve)
  2782. {
  2783.   unsigned short e[NI];
  2784.   char *s, *p;
  2785.   int i, j, k;
  2786.   int orig_ndigits;
  2787.   LDPARMS rnd;
  2788.   LDPARMS *ldp = &rnd;
  2789.   char *outstr;
  2790.   char outbuf[NDEC + MAX_EXP_DIGITS + 10];
  2791.   union uconv du;
  2792.   du.d = d;
  2793.  
  2794.   orig_ndigits = ndigits;
  2795.   rnd.rlast = -1;
  2796.   rnd.rndprc = NBITS;
  2797.  
  2798.   _REENT_CHECK_MP (ptr);
  2799.  
  2800. /* reentrancy addition to use mprec storage pool */
  2801.   if (_REENT_MP_RESULT (ptr))
  2802.     {
  2803.       _REENT_MP_RESULT (ptr)->_k = _REENT_MP_RESULT_K (ptr);
  2804.       _REENT_MP_RESULT (ptr)->_maxwds = 1 << _REENT_MP_RESULT_K (ptr);
  2805.       Bfree (ptr, _REENT_MP_RESULT (ptr));
  2806.       _REENT_MP_RESULT (ptr) = 0;
  2807.     }
  2808.  
  2809. #if LDBL_MANT_DIG == 24
  2810.   e24toe (&du.pe, e, ldp);
  2811. #elif LDBL_MANT_DIG == 53
  2812.   e53toe (&du.pe, e, ldp);
  2813. #elif LDBL_MANT_DIG == 64
  2814.   e64toe (&du.pe, e, ldp);
  2815. #else
  2816.   e113toe (&du.pe, e, ldp);
  2817. #endif
  2818.  
  2819.   if (eisneg (e))
  2820.     *sign = 1;
  2821.   else
  2822.     *sign = 0;
  2823. /* Mode 3 is "f" format.  */
  2824.   if (mode != 3)
  2825.     ndigits -= 1;
  2826. /* Mode 0 is for %.999 format, which is supposed to give a
  2827.    minimum length string that will convert back to the same binary value.
  2828.    For now, just ask for 20 digits which is enough but sometimes too many.  */
  2829.   if (mode == 0)
  2830.     ndigits = 20;
  2831.  
  2832. /* This sanity limit must agree with the corresponding one in etoasc, to
  2833.    keep straight the returned value of outexpon.  */
  2834.   if (ndigits > NDEC)
  2835.     ndigits = NDEC;
  2836.  
  2837.   etoasc (e, outbuf, ndigits, mode, ldp);
  2838.   s = outbuf;
  2839.   if (eisinf (e) || eisnan (e))
  2840.     {
  2841.       *decpt = 9999;
  2842.       goto stripspaces;
  2843.     }
  2844.   *decpt = ldp->outexpon + 1;
  2845.  
  2846. /* Transform the string returned by etoasc into what the caller wants.  */
  2847.  
  2848. /* Look for decimal point and delete it from the string. */
  2849.   s = outbuf;
  2850.   while (*s != '\0')
  2851.     {
  2852.       if (*s == '.')
  2853.         goto yesdecpt;
  2854.       ++s;
  2855.     }
  2856.   goto nodecpt;
  2857.  
  2858. yesdecpt:
  2859.  
  2860. /* Delete the decimal point.  */
  2861.   while (*s != '\0')
  2862.     {
  2863.       *s = *(s + 1);
  2864.       ++s;
  2865.     }
  2866.  
  2867. nodecpt:
  2868.  
  2869. /* Back up over the exponent field. */
  2870.   while (*s != 'E' && s > outbuf)
  2871.     --s;
  2872.   *s = '\0';
  2873.  
  2874. stripspaces:
  2875.  
  2876. /* Strip leading spaces and sign. */
  2877.   p = outbuf;
  2878.   while (*p == ' ' || *p == '-')
  2879.     ++p;
  2880.  
  2881. /* Find new end of string.  */
  2882.   s = outbuf;
  2883.   while ((*s++ = *p++) != '\0')
  2884.     ;
  2885.   --s;
  2886.  
  2887. /* Strip trailing zeros.  */
  2888.   if (mode == 2)
  2889.     k = 1;
  2890.   else if (ndigits > ldp->outexpon)
  2891.     k = ndigits;
  2892.   else
  2893.     k = ldp->outexpon;
  2894.  
  2895.   while (*(s - 1) == '0' && ((s - outbuf) > k))
  2896.     *(--s) = '\0';
  2897.  
  2898. /* In f format, flush small off-scale values to zero.
  2899.    Rounding has been taken care of by etoasc. */
  2900.   if (mode == 3 && ((ndigits + ldp->outexpon) < 0))
  2901.     {
  2902.       s = outbuf;
  2903.       *s = '\0';
  2904.       *decpt = 0;
  2905.     }
  2906.  
  2907. /* reentrancy addition to use mprec storage pool */
  2908. /* we want to have enough space to hold the formatted result */
  2909.  
  2910.   if (mode == 3)                /* f format, account for sign + dec digits + decpt + frac */
  2911.     i = *decpt + orig_ndigits + 3;
  2912.   else                          /* account for sign + max precision digs + E + exp sign + exponent */
  2913.     i = orig_ndigits + MAX_EXP_DIGITS + 4;
  2914.  
  2915.   j = sizeof (__ULong);
  2916.   for (_REENT_MP_RESULT_K (ptr) = 0;
  2917.        sizeof (_Bigint) - sizeof (__ULong) + j <= i; j <<= 1)
  2918.     _REENT_MP_RESULT_K (ptr)++;
  2919.   _REENT_MP_RESULT (ptr) = Balloc (ptr, _REENT_MP_RESULT_K (ptr));
  2920.  
  2921. /* Copy from internal temporary buffer to permanent buffer.  */
  2922.   outstr = (char *) _REENT_MP_RESULT (ptr);
  2923.   strcpy (outstr, outbuf);
  2924.  
  2925.   if (rve)
  2926.     *rve = outstr + (s - outbuf);
  2927.  
  2928.   return outstr;
  2929. }
  2930.  
  2931. /* Routine used to tell if long double is NaN or Infinity or regular number.
  2932.    Returns:  0 = regular number
  2933.              1 = Nan
  2934.              2 = Infinity
  2935. */
  2936. int
  2937. _ldcheck (long double *d)
  2938. {
  2939.   unsigned short e[NI];
  2940.   LDPARMS rnd;
  2941.   LDPARMS *ldp = &rnd;
  2942.  
  2943.   union uconv du;
  2944.  
  2945.   rnd.rlast = -1;
  2946.   rnd.rndprc = NBITS;
  2947.   du.d = *d;
  2948. #if LDBL_MANT_DIG == 24
  2949.   e24toe (&du.pe, e, ldp);
  2950. #elif LDBL_MANT_DIG == 53
  2951.   e53toe (&du.pe, e, ldp);
  2952. #elif LDBL_MANT_DIG == 64
  2953.   e64toe (&du.pe, e, ldp);
  2954. #else
  2955.   e113toe (&du.pe, e, ldp);
  2956. #endif
  2957.  
  2958.   if ((e[NE - 1] & 0x7fff) == 0x7fff)
  2959.     {
  2960. #ifdef NANS
  2961.       if (eisnan (e))
  2962.         return (1);
  2963. #endif
  2964.       return (2);
  2965.     }
  2966.   else
  2967.     return (0);
  2968. }                               /* _ldcheck */
  2969.  
  2970. static void
  2971. etoasc (short unsigned int *x, char *string, int ndigits, int outformat,
  2972.         LDPARMS * ldp)
  2973. {
  2974.   long digit;
  2975.   unsigned short y[NI], t[NI], u[NI], w[NI];
  2976.   _CONST unsigned short *p, *r, *ten;
  2977.   unsigned short sign;
  2978.   int i, j, k, expon, rndsav, ndigs;
  2979.   char *s, *ss;
  2980.   unsigned short m;
  2981.   unsigned short *equot = ldp->equot;
  2982.  
  2983.   ndigs = ndigits;
  2984.   rndsav = ldp->rndprc;
  2985. #ifdef NANS
  2986.   if (eisnan (x))
  2987.     {
  2988.       sprintf (string, " NaN ");
  2989.       expon = 9999;
  2990.       goto bxit;
  2991.     }
  2992. #endif
  2993.   ldp->rndprc = NBITS;          /* set to full precision */
  2994.   emov (x, y);                  /* retain external format */
  2995.   if (y[NE - 1] & 0x8000)
  2996.     {
  2997.       sign = 0xffff;
  2998.       y[NE - 1] &= 0x7fff;
  2999.     }
  3000.   else
  3001.     {
  3002.       sign = 0;
  3003.     }
  3004.   expon = 0;
  3005.   ten = &etens[NTEN][0];
  3006.   emov (eone, t);
  3007. /* Test for zero exponent */
  3008.   if (y[NE - 1] == 0)
  3009.     {
  3010.       for (k = 0; k < NE - 1; k++)
  3011.         {
  3012.           if (y[k] != 0)
  3013.             goto tnzro;         /* denormalized number */
  3014.         }
  3015.       goto isone;               /* legal all zeros */
  3016.     }
  3017. tnzro:
  3018.  
  3019. /* Test for infinity.
  3020.  */
  3021.   if (y[NE - 1] == 0x7fff)
  3022.     {
  3023.       if (sign)
  3024.         sprintf (string, " -Infinity ");
  3025.       else
  3026.         sprintf (string, " Infinity ");
  3027.       expon = 9999;
  3028.       goto bxit;
  3029.     }
  3030.  
  3031. /* Test for exponent nonzero but significand denormalized.
  3032.  * This is an error condition.
  3033.  */
  3034.   if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
  3035.     {
  3036.       mtherr ("etoasc", DOMAIN);
  3037.       sprintf (string, "NaN");
  3038.       expon = 9999;
  3039.       goto bxit;
  3040.     }
  3041.  
  3042. /* Compare to 1.0 */
  3043.   i = ecmp (eone, y);
  3044.   if (i == 0)
  3045.     goto isone;
  3046.  
  3047.   if (i < 0)
  3048.     {                           /* Number is greater than 1 */
  3049. /* Convert significand to an integer and strip trailing decimal zeros. */
  3050.       emov (y, u);
  3051.       u[NE - 1] = EXONE + NBITS - 1;
  3052.  
  3053.       p = &etens[NTEN - 4][0];
  3054.       m = 16;
  3055.       do
  3056.         {
  3057.           ediv (p, u, t, ldp);
  3058.           efloor (t, w, ldp);
  3059.           for (j = 0; j < NE - 1; j++)
  3060.             {
  3061.               if (t[j] != w[j])
  3062.                 goto noint;
  3063.             }
  3064.           emov (t, u);
  3065.           expon += (int) m;
  3066.         noint:
  3067.           p += NE;
  3068.           m >>= 1;
  3069.         }
  3070.       while (m != 0);
  3071.  
  3072. /* Rescale from integer significand */
  3073.       u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
  3074.       emov (u, y);
  3075. /* Find power of 10 */
  3076.       emov (eone, t);
  3077.       m = MAXP;
  3078.       p = &etens[0][0];
  3079.       while (ecmp (ten, u) <= 0)
  3080.         {
  3081.           if (ecmp (p, u) <= 0)
  3082.             {
  3083.               ediv (p, u, u, ldp);
  3084.               emul (p, t, t, ldp);
  3085.               expon += (int) m;
  3086.             }
  3087.           m >>= 1;
  3088.           if (m == 0)
  3089.             break;
  3090.           p += NE;
  3091.         }
  3092.     }
  3093.   else
  3094.     {                           /* Number is less than 1.0 */
  3095. /* Pad significand with trailing decimal zeros. */
  3096.       if (y[NE - 1] == 0)
  3097.         {
  3098.           while ((y[NE - 2] & 0x8000) == 0)
  3099.             {
  3100.               emul (ten, y, y, ldp);
  3101.               expon -= 1;
  3102.             }
  3103.         }
  3104.       else
  3105.         {
  3106.           emovi (y, w);
  3107.           for (i = 0; i < NDEC + 1; i++)
  3108.             {
  3109.               if ((w[NI - 1] & 0x7) != 0)
  3110.                 break;
  3111. /* multiply by 10 */
  3112.               emovz (w, u);
  3113.               eshdn1 (u);
  3114.               eshdn1 (u);
  3115.               eaddm (w, u);
  3116.               u[1] += 3;
  3117.               while (u[2] != 0)
  3118.                 {
  3119.                   eshdn1 (u);
  3120.                   u[1] += 1;
  3121.                 }
  3122.               if (u[NI - 1] != 0)
  3123.                 break;
  3124.               if (eone[NE - 1] <= u[1])
  3125.                 break;
  3126.               emovz (u, w);
  3127.               expon -= 1;
  3128.             }
  3129.           emovo (w, y, ldp);
  3130.         }
  3131.       k = -MAXP;
  3132.       p = &emtens[0][0];
  3133.       r = &etens[0][0];
  3134.       emov (y, w);
  3135.       emov (eone, t);
  3136.       while (ecmp (eone, w) > 0)
  3137.         {
  3138.           if (ecmp (p, w) >= 0)
  3139.             {
  3140.               emul (r, w, w, ldp);
  3141.               emul (r, t, t, ldp);
  3142.               expon += k;
  3143.             }
  3144.           k /= 2;
  3145.           if (k == 0)
  3146.             break;
  3147.           p += NE;
  3148.           r += NE;
  3149.         }
  3150.       ediv (t, eone, t, ldp);
  3151.     }
  3152. isone:
  3153. /* Find the first (leading) digit. */
  3154.   emovi (t, w);
  3155.   emovz (w, t);
  3156.   emovi (y, w);
  3157.   emovz (w, y);
  3158.   eiremain (t, y, ldp);
  3159.   digit = equot[NI - 1];
  3160.   while ((digit == 0) && (ecmp (y, ezero) != 0))
  3161.     {
  3162.       eshup1 (y);
  3163.       emovz (y, u);
  3164.       eshup1 (u);
  3165.       eshup1 (u);
  3166.       eaddm (u, y);
  3167.       eiremain (t, y, ldp);
  3168.       digit = equot[NI - 1];
  3169.       expon -= 1;
  3170.     }
  3171.   s = string;
  3172.   if (sign)
  3173.     *s++ = '-';
  3174.   else
  3175.     *s++ = ' ';
  3176. /* Examine number of digits requested by caller. */
  3177.   if (outformat == 3)
  3178.     ndigs += expon;
  3179. /*
  3180. else if( ndigs < 0 )
  3181.         ndigs = 0;
  3182. */
  3183.   if (ndigs > NDEC)
  3184.     ndigs = NDEC;
  3185.   if (digit == 10)
  3186.     {
  3187.       *s++ = '1';
  3188.       *s++ = '.';
  3189.       if (ndigs > 0)
  3190.         {
  3191.           *s++ = '0';
  3192.           ndigs -= 1;
  3193.         }
  3194.       expon += 1;
  3195.       if (ndigs < 0)
  3196.         {
  3197.           ss = s;
  3198.           goto doexp;
  3199.         }
  3200.     }
  3201.   else
  3202.     {
  3203.       *s++ = (char) digit + '0';
  3204.       *s++ = '.';
  3205.     }
  3206. /* Generate digits after the decimal point. */
  3207.   for (k = 0; k <= ndigs; k++)
  3208.     {
  3209. /* multiply current number by 10, without normalizing */
  3210.       eshup1 (y);
  3211.       emovz (y, u);
  3212.       eshup1 (u);
  3213.       eshup1 (u);
  3214.       eaddm (u, y);
  3215.       eiremain (t, y, ldp);
  3216.       *s++ = (char) equot[NI - 1] + '0';
  3217.     }
  3218.   digit = equot[NI - 1];
  3219.   --s;
  3220.   ss = s;
  3221. /* round off the ASCII string */
  3222.   if (digit > 4)
  3223.     {
  3224. /* Test for critical rounding case in ASCII output. */
  3225.       if (digit == 5)
  3226.         {
  3227.           emovo (y, t, ldp);
  3228.           if (ecmp (t, ezero) != 0)
  3229.             goto roun;          /* round to nearest */
  3230.           if (ndigs < 0 || (*(s - 1 - (*(s - 1) == '.')) & 1) == 0)
  3231.             goto doexp;         /* round to even */
  3232.         }
  3233. /* Round up and propagate carry-outs */
  3234.     roun:
  3235.       --s;
  3236.       k = *s & 0x7f;
  3237. /* Carry out to most significant digit? */
  3238.       if (ndigs < 0)
  3239.         {
  3240.           /* This will print like "1E-6". */
  3241.           *s = '1';
  3242.           expon += 1;
  3243.           goto doexp;
  3244.         }
  3245.       else if (k == '.')
  3246.         {
  3247.           --s;
  3248.           k = *s;
  3249.           k += 1;
  3250.           *s = (char) k;
  3251. /* Most significant digit carries to 10? */
  3252.           if (k > '9')
  3253.             {
  3254.               expon += 1;
  3255.               *s = '1';
  3256.             }
  3257.           goto doexp;
  3258.         }
  3259. /* Round up and carry out from less significant digits */
  3260.       k += 1;
  3261.       *s = (char) k;
  3262.       if (k > '9')
  3263.         {
  3264.           *s = '0';
  3265.           goto roun;
  3266.         }
  3267.     }
  3268. doexp:
  3269. #ifdef __GO32__
  3270.   if (expon >= 0)
  3271.     sprintf (ss, "e+%02d", expon);
  3272.   else
  3273.     sprintf (ss, "e-%02d", -expon);
  3274. #else
  3275.   sprintf (ss, "E%d", expon);
  3276. #endif
  3277. bxit:
  3278.   ldp->rndprc = rndsav;
  3279.   ldp->outexpon = expon;
  3280. }
  3281.  
  3282.  
  3283. #if 0 /* Broken, unusable implementation of strtold */
  3284.  
  3285. /*
  3286. ;                                                               ASCTOQ
  3287. ;               ASCTOQ.MAC              LATEST REV: 11 JAN 84
  3288. ;                                       SLM, 3 JAN 78
  3289. ;
  3290. ;       Convert ASCII string to quadruple precision floating point
  3291. ;
  3292. ;               Numeric input is free field decimal number
  3293. ;               with max of 15 digits with or without
  3294. ;               decimal point entered as ASCII from teletype.
  3295. ;       Entering E after the number followed by a second
  3296. ;       number causes the second number to be interpreted
  3297. ;       as a power of 10 to be multiplied by the first number
  3298. ;       (i.e., "scientific" notation).
  3299. ;
  3300. ;       Usage:
  3301. ;               asctoq( string, q );
  3302. */
  3303.  
  3304. long double
  3305. _strtold (char *s, char **se)
  3306. {
  3307.   union uconv x;
  3308.   LDPARMS rnd;
  3309.   LDPARMS *ldp = &rnd;
  3310.   int lenldstr;
  3311.  
  3312.   rnd.rlast = -1;
  3313.   rnd.rndprc = NBITS;
  3314.  
  3315.   lenldstr = asctoeg (s, &x.pe, LDBL_MANT_DIG, ldp);
  3316.   if (se)
  3317.     *se = s + lenldstr;
  3318.   return x.d;
  3319. }
  3320.  
  3321. #define REASONABLE_LEN 200
  3322.  
  3323. static int
  3324. asctoeg (char *ss, short unsigned int *y, int oprec, LDPARMS * ldp)
  3325. {
  3326.   unsigned short yy[NI], xt[NI], tt[NI];
  3327.   int esign, decflg, sgnflg, nexp, exp, prec, lost;
  3328.   int k, trail, c, rndsav;
  3329.   long lexp;
  3330.   unsigned short nsign;
  3331.   _CONST unsigned short *p;
  3332.   char *sp, *s, *lstr;
  3333.   int lenldstr;
  3334.   int mflag = 0;
  3335.   char tmpstr[REASONABLE_LEN];
  3336.  
  3337. /* Copy the input string. */
  3338.   c = strlen (ss) + 2;
  3339.   if (c <= REASONABLE_LEN)
  3340.     lstr = tmpstr;
  3341.   else
  3342.     {
  3343.       lstr = (char *) calloc (c, 1);
  3344.       mflag = 1;
  3345.     }
  3346.   s = ss;
  3347.   lenldstr = 0;
  3348.   while (*s == ' ')             /* skip leading spaces */
  3349.     {
  3350.       ++s;
  3351.       ++lenldstr;
  3352.     }
  3353.   sp = lstr;
  3354.   for (k = 0; k < c; k++)
  3355.     {
  3356.       if ((*sp++ = *s++) == '\0')
  3357.         break;
  3358.     }
  3359.   *sp = '\0';
  3360.   s = lstr;
  3361.  
  3362.   rndsav = ldp->rndprc;
  3363.   ldp->rndprc = NBITS;          /* Set to full precision */
  3364.   lost = 0;
  3365.   nsign = 0;
  3366.   decflg = 0;
  3367.   sgnflg = 0;
  3368.   nexp = 0;
  3369.   exp = 0;
  3370.   prec = 0;
  3371.   ecleaz (yy);
  3372.   trail = 0;
  3373.  
  3374. nxtcom:
  3375.   k = *s - '0';
  3376.   if ((k >= 0) && (k <= 9))
  3377.     {
  3378. /* Ignore leading zeros */
  3379.       if ((prec == 0) && (decflg == 0) && (k == 0))
  3380.         goto donchr;
  3381. /* Identify and strip trailing zeros after the decimal point. */
  3382.       if ((trail == 0) && (decflg != 0))
  3383.         {
  3384.           sp = s;
  3385.           while ((*sp >= '0') && (*sp <= '9'))
  3386.             ++sp;
  3387. /* Check for syntax error */
  3388.           c = *sp & 0x7f;
  3389.           if ((c != 'e') && (c != 'E') && (c != '\0')
  3390.               && (c != '\n') && (c != '\r') && (c != ' ') && (c != ','))
  3391.             goto error;
  3392.           --sp;
  3393.           while (*sp == '0')
  3394.             *sp-- = 'z';
  3395.           trail = 1;
  3396.           if (*s == 'z')
  3397.             goto donchr;
  3398.         }
  3399. /* If enough digits were given to more than fill up the yy register,
  3400.  * continuing until overflow into the high guard word yy[2]
  3401.  * guarantees that there will be a roundoff bit at the top
  3402.  * of the low guard word after normalization.
  3403.  */
  3404.       if (yy[2] == 0)
  3405.         {
  3406.           if (decflg)
  3407.             nexp += 1;          /* count digits after decimal point */
  3408.           eshup1 (yy);          /* multiply current number by 10 */
  3409.           emovz (yy, xt);
  3410.           eshup1 (xt);
  3411.           eshup1 (xt);
  3412.           eaddm (xt, yy);
  3413.           ecleaz (xt);
  3414.           xt[NI - 2] = (unsigned short) k;
  3415.           eaddm (xt, yy);
  3416.         }
  3417.       else
  3418.         {
  3419.           /* Mark any lost non-zero digit.  */
  3420.           lost |= k;
  3421.           /* Count lost digits before the decimal point.  */
  3422.           if (decflg == 0)
  3423.             nexp -= 1;
  3424.         }
  3425.       prec += 1;
  3426.       goto donchr;
  3427.     }
  3428.  
  3429.   switch (*s)
  3430.     {
  3431.     case 'z':
  3432.       break;
  3433.     case 'E':
  3434.     case 'e':
  3435.       goto expnt;
  3436.     case '.':                   /* decimal point */
  3437.       if (decflg)
  3438.         goto error;
  3439.       ++decflg;
  3440.       break;
  3441.     case '-':
  3442.       nsign = 0xffff;
  3443.       if (sgnflg)
  3444.         goto error;
  3445.       ++sgnflg;
  3446.       break;
  3447.     case '+':
  3448.       if (sgnflg)
  3449.         goto error;
  3450.       ++sgnflg;
  3451.       break;
  3452.     case ',':
  3453.     case ' ':
  3454.     case '\0':
  3455.     case '\n':
  3456.     case '\r':
  3457.       goto daldone;
  3458.     case 'i':
  3459.     case 'I':
  3460.       goto infinite;
  3461.     default:
  3462.     error:
  3463. #ifdef NANS
  3464.       enan (yy, NI * 16);
  3465. #else
  3466.       mtherr ("asctoe", DOMAIN);
  3467.       ecleaz (yy);
  3468. #endif
  3469.       goto aexit;
  3470.     }
  3471. donchr:
  3472.   ++s;
  3473.   goto nxtcom;
  3474.  
  3475. /* Exponent interpretation */
  3476. expnt:
  3477.  
  3478.   esign = 1;
  3479.   exp = 0;
  3480.   ++s;
  3481. /* check for + or - */
  3482.   if (*s == '-')
  3483.     {
  3484.       esign = -1;
  3485.       ++s;
  3486.     }
  3487.   if (*s == '+')
  3488.     ++s;
  3489.   while ((*s >= '0') && (*s <= '9'))
  3490.     {
  3491.       exp *= 10;
  3492.       exp += *s++ - '0';
  3493.       if (exp > 4977)
  3494.         {
  3495.           if (esign < 0)
  3496.             goto zero;
  3497.           else
  3498.             goto infinite;
  3499.         }
  3500.     }
  3501.   if (esign < 0)
  3502.     exp = -exp;
  3503.   if (exp > 4932)
  3504.     {
  3505.     infinite:
  3506.       ecleaz (yy);
  3507.       yy[E] = 0x7fff;           /* infinity */
  3508.       goto aexit;
  3509.     }
  3510.   if (exp < -4977)
  3511.     {
  3512.     zero:
  3513.       ecleaz (yy);
  3514.       goto aexit;
  3515.     }
  3516.  
  3517. daldone:
  3518.   nexp = exp - nexp;
  3519. /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
  3520.   while ((nexp > 0) && (yy[2] == 0))
  3521.     {
  3522.       emovz (yy, xt);
  3523.       eshup1 (xt);
  3524.       eshup1 (xt);
  3525.       eaddm (yy, xt);
  3526.       eshup1 (xt);
  3527.       if (xt[2] != 0)
  3528.         break;
  3529.       nexp -= 1;
  3530.       emovz (xt, yy);
  3531.     }
  3532.   if ((k = enormlz (yy)) > NBITS)
  3533.     {
  3534.       ecleaz (yy);
  3535.       goto aexit;
  3536.     }
  3537.   lexp = (EXONE - 1 + NBITS) - k;
  3538.   emdnorm (yy, lost, 0, lexp, 64, ldp);
  3539. /* convert to external format */
  3540.  
  3541.  
  3542. /* Multiply by 10**nexp.  If precision is 64 bits,
  3543.  * the maximum relative error incurred in forming 10**n
  3544.  * for 0 <= n <= 324 is 8.2e-20, at 10**180.
  3545.  * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
  3546.  * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
  3547.  */
  3548.   lexp = yy[E];
  3549.   if (nexp == 0)
  3550.     {
  3551.       k = 0;
  3552.       goto expdon;
  3553.     }
  3554.   esign = 1;
  3555.   if (nexp < 0)
  3556.     {
  3557.       nexp = -nexp;
  3558.       esign = -1;
  3559.       if (nexp > 4096)
  3560.         {                       /* Punt.  Can't handle this without 2 divides. */
  3561.           emovi (etens[0], tt);
  3562.           lexp -= tt[E];
  3563.           k = edivm (tt, yy, ldp);
  3564.           lexp += EXONE;
  3565.           nexp -= 4096;
  3566.         }
  3567.     }
  3568.   p = &etens[NTEN][0];
  3569.   emov (eone, xt);
  3570.   exp = 1;
  3571.   do
  3572.     {
  3573.       if (exp & nexp)
  3574.         emul (p, xt, xt, ldp);
  3575.       p -= NE;
  3576.       exp = exp + exp;
  3577.     }
  3578.   while (exp <= MAXP);
  3579.  
  3580.   emovi (xt, tt);
  3581.   if (esign < 0)
  3582.     {
  3583.       lexp -= tt[E];
  3584.       k = edivm (tt, yy, ldp);
  3585.       lexp += EXONE;
  3586.     }
  3587.   else
  3588.     {
  3589.       lexp += tt[E];
  3590.       k = emulm (tt, yy, ldp);
  3591.       lexp -= EXONE - 1;
  3592.     }
  3593.  
  3594. expdon:
  3595.  
  3596. /* Round and convert directly to the destination type */
  3597.   if (oprec == 53)
  3598.     lexp -= EXONE - 0x3ff;
  3599.   else if (oprec == 24)
  3600.     lexp -= EXONE - 0177;
  3601. #ifdef DEC
  3602.   else if (oprec == 56)
  3603.     lexp -= EXONE - 0201;
  3604. #endif
  3605.   ldp->rndprc = oprec;
  3606.   emdnorm (yy, k, 0, lexp, 64, ldp);
  3607.  
  3608. aexit:
  3609.  
  3610.   ldp->rndprc = rndsav;
  3611.   yy[0] = nsign;
  3612.   switch (oprec)
  3613.     {
  3614. #ifdef DEC
  3615.     case 56:
  3616.       todec (yy, y);            /* see etodec.c */
  3617.       break;
  3618. #endif
  3619. #if LDBL_MANT_DIG == 53
  3620.     case 53:
  3621.       toe53 (yy, y);
  3622.       break;
  3623. #elif LDBL_MANT_DIG == 24
  3624.     case 24:
  3625.       toe24 (yy, y);
  3626.       break;
  3627. #elif LDBL_MANT_DIG == 64
  3628.     case 64:
  3629.       toe64 (yy, y);
  3630.       break;
  3631. #elif LDBL_MANT_DIG == 113
  3632.     case 113:
  3633.       toe113 (yy, y);
  3634.       break;
  3635. #else
  3636.     case NBITS:
  3637.       emovo (yy, y, ldp);
  3638.       break;
  3639. #endif
  3640.     }
  3641.   lenldstr += s - lstr;
  3642.   if (mflag)
  3643.     free (lstr);
  3644.   return lenldstr;
  3645. }
  3646.  
  3647. #endif
  3648.  
  3649. /* y = largest integer not greater than x
  3650.  * (truncated toward minus infinity)
  3651.  *
  3652.  * unsigned short x[NE], y[NE]
  3653.  * LDPARMS *ldp
  3654.  *
  3655.  * efloor( x, y, ldp );
  3656.  */
  3657. static _CONST unsigned short bmask[] = {
  3658.   0xffff,
  3659.   0xfffe,
  3660.   0xfffc,
  3661.   0xfff8,
  3662.   0xfff0,
  3663.   0xffe0,
  3664.   0xffc0,
  3665.   0xff80,
  3666.   0xff00,
  3667.   0xfe00,
  3668.   0xfc00,
  3669.   0xf800,
  3670.   0xf000,
  3671.   0xe000,
  3672.   0xc000,
  3673.   0x8000,
  3674.   0x0000,
  3675. };
  3676.  
  3677. static void
  3678. efloor (short unsigned int *x, short unsigned int *y, LDPARMS * ldp)
  3679. {
  3680.   register unsigned short *p;
  3681.   int e, expon, i;
  3682.   unsigned short f[NE];
  3683.  
  3684.   emov (x, f);                  /* leave in external format */
  3685.   expon = (int) f[NE - 1];
  3686.   e = (expon & 0x7fff) - (EXONE - 1);
  3687.   if (e <= 0)
  3688.     {
  3689.       eclear (y);
  3690.       goto isitneg;
  3691.     }
  3692. /* number of bits to clear out */
  3693.   e = NBITS - e;
  3694.   emov (f, y);
  3695.   if (e <= 0)
  3696.     return;
  3697.  
  3698.   p = &y[0];
  3699.   while (e >= 16)
  3700.     {
  3701.       *p++ = 0;
  3702.       e -= 16;
  3703.     }
  3704. /* clear the remaining bits */
  3705.   *p &= bmask[e];
  3706. /* truncate negatives toward minus infinity */
  3707. isitneg:
  3708.  
  3709.   if ((unsigned short) expon & (unsigned short) 0x8000)
  3710.     {
  3711.       for (i = 0; i < NE - 1; i++)
  3712.         {
  3713.           if (f[i] != y[i])
  3714.             {
  3715.               esub (eone, y, y, ldp);
  3716.               break;
  3717.             }
  3718.         }
  3719.     }
  3720. }
  3721.  
  3722.  
  3723.  
  3724. static void
  3725. eiremain (short unsigned int *den, short unsigned int *num, LDPARMS * ldp)
  3726. {
  3727.   long ld, ln;
  3728.   unsigned short j;
  3729.   unsigned short *equot = ldp->equot;
  3730.  
  3731.   ld = den[E];
  3732.   ld -= enormlz (den);
  3733.   ln = num[E];
  3734.   ln -= enormlz (num);
  3735.   ecleaz (equot);
  3736.   while (ln >= ld)
  3737.     {
  3738.       if (ecmpm (den, num) <= 0)
  3739.         {
  3740.           esubm (den, num);
  3741.           j = 1;
  3742.         }
  3743.       else
  3744.         {
  3745.           j = 0;
  3746.         }
  3747.       eshup1 (equot);
  3748.       equot[NI - 1] |= j;
  3749.       eshup1 (num);
  3750.       ln -= 1;
  3751.     }
  3752.   emdnorm (num, 0, 0, ln, 0, ldp);
  3753. }
  3754.  
  3755. /* NaN bit patterns
  3756.  */
  3757. #ifdef MIEEE
  3758. #if !defined(__mips)
  3759. static _CONST unsigned short nan113[8] = {
  3760.   0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff
  3761. };
  3762.  
  3763. static _CONST unsigned short nan64[6] = {
  3764.   0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff
  3765. };
  3766. static _CONST unsigned short nan53[4] = { 0x7fff, 0xffff, 0xffff, 0xffff };
  3767. static _CONST unsigned short nan24[2] = { 0x7fff, 0xffff };
  3768. #elif defined(__mips_nan2008)   /* __mips */
  3769. static _CONST unsigned short nan113[8] = { 0x7fff, 0x8000, 0, 0, 0, 0, 0, 0 };
  3770. static _CONST unsigned short nan64[6] = { 0x7fff, 0xc000, 0, 0, 0, 0 };
  3771. static _CONST unsigned short nan53[4] = { 0x7ff8, 0, 0, 0 };
  3772. static _CONST unsigned short nan24[2] = { 0x7fc0, 0 };
  3773. #else /* __mips && !__mips_nan2008 */
  3774. static _CONST unsigned short nan113[8] = {
  3775.   0x7fff, 0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff
  3776. };
  3777.  
  3778. static _CONST unsigned short nan64[6] = {
  3779.   0x7fff, 0xbfff, 0xffff, 0xffff, 0xffff, 0xffff
  3780. };
  3781. static _CONST unsigned short nan53[4] = { 0x7ff7, 0xffff, 0xffff, 0xffff };
  3782. static _CONST unsigned short nan24[2] = { 0x7fbf, 0xffff };
  3783. #endif /* __mips && !__mips_nan2008 */
  3784. #else /* !MIEEE */
  3785. #if !defined(__mips) || defined(__mips_nan2008)
  3786. static _CONST unsigned short nan113[8] = { 0, 0, 0, 0, 0, 0, 0x8000, 0x7fff };
  3787. static _CONST unsigned short nan64[6] = { 0, 0, 0, 0, 0xc000, 0x7fff };
  3788. static _CONST unsigned short nan53[4] = { 0, 0, 0, 0x7ff8 };
  3789. static _CONST unsigned short nan24[2] = { 0, 0x7fc0 };
  3790. #else /* __mips && !__mips_nan2008 */
  3791. static _CONST unsigned short nan113[8] = {
  3792.   0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0x7fff, 0x7fff
  3793. };
  3794.  
  3795. static _CONST unsigned short nan64[6] = {
  3796.   0xffff, 0xffff, 0xffff, 0xffff, 0xbfff, 0x7fff
  3797. };
  3798. static _CONST unsigned short nan53[4] = { 0xffff, 0xffff, 0xffff, 0x7ff7 };
  3799. static _CONST unsigned short nan24[2] = { 0xffff, 0x7fbf };
  3800. #endif /* __mips && !__mips_nan2008 */
  3801. #endif /* !MIEEE */
  3802.  
  3803.  
  3804. static void
  3805. enan (short unsigned int *nan, int size)
  3806. {
  3807.   int i, n;
  3808.   _CONST unsigned short *p;
  3809.  
  3810.   switch (size)
  3811.     {
  3812. #ifndef DEC
  3813.     case 113:
  3814.       n = 8;
  3815.       p = nan113;
  3816.       break;
  3817.  
  3818.     case 64:
  3819.       n = 6;
  3820.       p = nan64;
  3821.       break;
  3822.  
  3823.     case 53:
  3824.       n = 4;
  3825.       p = nan53;
  3826.       break;
  3827.  
  3828.     case 24:
  3829.       n = 2;
  3830.       p = nan24;
  3831.       break;
  3832.  
  3833.     case NBITS:
  3834. #if !defined(__mips) || defined(__mips_nan2008)
  3835.       for (i = 0; i < NE - 2; i++)
  3836.         *nan++ = 0;
  3837.       *nan++ = 0xc000;
  3838. #else /* __mips && !__mips_nan2008 */
  3839.       for (i = 0; i < NE - 2; i++)
  3840.         *nan++ = 0xffff;
  3841.       *nan++ = 0xbfff;
  3842. #endif /* __mips && !__mips_nan2008 */
  3843.       *nan++ = 0x7fff;
  3844.       return;
  3845.  
  3846.     case NI * 16:
  3847.       *nan++ = 0;
  3848.       *nan++ = 0x7fff;
  3849.       *nan++ = 0;
  3850. #if !defined(__mips) || defined(__mips_nan2008)
  3851.       *nan++ = 0xc000;
  3852.       for (i = 4; i < NI - 1; i++)
  3853.         *nan++ = 0;
  3854. #else /* __mips && !__mips_nan2008 */
  3855.       *nan++ = 0xbfff;
  3856.       for (i = 4; i < NI - 1; i++)
  3857.         *nan++ = 0xffff;
  3858. #endif /* __mips && !__mips_nan2008 */
  3859.       *nan++ = 0;
  3860.       return;
  3861. #endif
  3862.     default:
  3863.       mtherr ("enan", DOMAIN);
  3864.       return;
  3865.     }
  3866.   for (i = 0; i < n; i++)
  3867.     *nan++ = *p++;
  3868. }
  3869.