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