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