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