(output_function_{pro,epi}logue): Use lea instead of add.w when

[pf3gnuchains/gcc-fork.git] / gcc / real.c

diff --git a/gcc/real.c b/gcc/real.c
index 518aaa8..f7e22ea 100644 (file)
--- a/gcc/real.c
+++ b/gcc/real.c
@@ -1,8 +1,7 @@
 /* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
 /* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,

-and support for XFmode IEEE extended real floating point arithmetic.
-Contributed by Stephen L. Moshier (moshier@world.std.com).
-
-   Copyright (C) 1993 Free Software Foundation, Inc.
+   and support for XFmode IEEE extended real floating point arithmetic.
+   Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
+   Contributed by Stephen L. Moshier (moshier@world.std.com).

 
 This file is part of GNU CC.
 
 
 This file is part of GNU CC.
 


@@ -18,7 +17,8 @@ GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with GNU CC; see the file COPYING.  If not, write to
 
 You should have received a copy of the GNU General Public License
 along with GNU CC; see the file COPYING.  If not, write to

-the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */
+the Free Software Foundation, 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA.  */

 
 #include <stdio.h>
 #include <errno.h>
 
 #include <stdio.h>
 #include <errno.h>


@@ -32,7 +32,7 @@ extern int errno;
 /* To enable support of XFmode extended real floating point, define
 LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
 
 /* To enable support of XFmode extended real floating point, define
 LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
 

-To support cross compilation between IEEE and VAX floating
+To support cross compilation between IEEE, VAX and IBM floating

 point formats, define REAL_ARITHMETIC in the tm.h file.
 
 In either case the machine files (tm.h) must not contain any code
 point formats, define REAL_ARITHMETIC in the tm.h file.
 
 In either case the machine files (tm.h) must not contain any code


@@ -45,27 +45,65 @@ The emulator defaults to the host's floating point format so that
 its decimal conversion functions can be used if desired (see
 real.h).
 
 its decimal conversion functions can be used if desired (see
 real.h).
 

-The first part of this file interfaces gcc to ieee.c, which is a
-floating point arithmetic suite that was not written with gcc in
-mind.  The interface is followed by ieee.c itself and related
-items. Avoid changing ieee.c unless you have suitable test
-programs available.  A special version of the PARANOIA floating
-point arithmetic tester, modified for this purpose, can be found
-on usc.edu : /pub/C-numanal/ieeetest.zoo.  Some tutorial
-information on ieee.c is given in my book: S. L. Moshier,
-_Methods and Programs for Mathematical Functions_, Prentice-Hall
-or Simon & Schuster Int'l, 1989.  A library of XFmode elementary
-transcendental functions can be obtained by ftp from
-research.att.com: netlib/cephes/ldouble.shar.Z  */
-
+The first part of this file interfaces gcc to a floating point
+arithmetic suite that was not written with gcc in mind.  Avoid
+changing the low-level arithmetic routines unless you have suitable
+test programs available.  A special version of the PARANOIA floating
+point arithmetic tester, modified for this purpose, can be found on
+usc.edu: /pub/C-numanal/ieeetest.zoo.  Other tests, and libraries of
+XFmode and TFmode transcendental functions, can be obtained by ftp from
+netlib.att.com: netlib/cephes.   */
+\f

 /* Type of computer arithmetic.
 /* Type of computer arithmetic.

- * Only one of DEC, MIEEE, IBMPC, or UNK should get defined.
- * The following modification converts gcc macros into the ones
- * used by ieee.c.
- *
- * Note: long double formats differ between IBMPC and MIEEE
- * by more than just endian-ness.
- */
+   Only one of DEC, IBM, IEEE, or UNK should get defined.
+
+   `IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
+   to big-endian IEEE floating-point data structure.  This definition
+   should work in SFmode `float' type and DFmode `double' type on
+   virtually all big-endian IEEE machines.  If LONG_DOUBLE_TYPE_SIZE
+   has been defined to be 96, then IEEE also invokes the particular
+   XFmode (`long double' type) data structure used by the Motorola
+   680x0 series processors.
+
+   `IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
+   little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
+   has been defined to be 96, then IEEE also invokes the particular
+   XFmode `long double' data structure used by the Intel 80x86 series
+   processors.
+
+   `DEC' refers specifically to the Digital Equipment Corp PDP-11
+   and VAX floating point data structure.  This model currently
+   supports no type wider than DFmode.
+
+   `IBM' refers specifically to the IBM System/370 and compatible
+   floating point data structure.  This model currently supports
+   no type wider than DFmode.  The IBM conversions were contributed by
+   frank@atom.ansto.gov.au (Frank Crawford).
+
+   If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
+   then `long double' and `double' are both implemented, but they
+   both mean DFmode.  In this case, the software floating-point
+   support available here is activated by writing
+      #define REAL_ARITHMETIC
+   in tm.h. 
+
+   The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
+   and may deactivate XFmode since `long double' is used to refer
+   to both modes.
+
+   The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
+   contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
+   separate the floating point unit's endian-ness from that of
+   the integer addressing.  This permits one to define a big-endian
+   FPU on a little-endian machine (e.g., ARM).  An extension to
+   BYTES_BIG_ENDIAN may be required for some machines in the future.
+   These optional macros may be defined in tm.h.  In real.h, they
+   default to WORDS_BIG_ENDIAN, etc., so there is no need to define
+   them for any normal host or target machine on which the floats
+   and the integers have the same endian-ness.   */
+
+
+/* The following converts gcc macros into the ones used by this file.  */

 
 /* REAL_ARITHMETIC defined means that macros in real.h are
    defined to call emulator functions.  */
 
 /* REAL_ARITHMETIC defined means that macros in real.h are
    defined to call emulator functions.  */


@@ -75,81 +113,65 @@ research.att.com: netlib/cephes/ldouble.shar.Z  */
 /* PDP-11, Pro350, VAX: */
 #define DEC 1
 #else /* it's not VAX */
 /* PDP-11, Pro350, VAX: */
 #define DEC 1
 #else /* it's not VAX */

+#if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
+/* IBM System/370 style */
+#define IBM 1
+#else /* it's also not an IBM */

 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT

-#if WORDS_BIG_ENDIAN
-/* Motorola IEEE, high order words come first (Sun workstation): */
-#define MIEEE 1
-#else /* not big-endian */
-/* Intel IEEE, low order words come first:
- */
-#define IBMPC 1
-#endif /*  big-endian */
+#define IEEE

 #else /* it's not IEEE either */
 #else /* it's not IEEE either */

-/* UNKnown arithmetic.  We don't support this and can't go on. */
+/* UNKnown arithmetic.  We don't support this and can't go on.  */

 unknown arithmetic type
 #define UNK 1
 #endif /* not IEEE */
 unknown arithmetic type
 #define UNK 1
 #endif /* not IEEE */

+#endif /* not IBM */

 #endif /* not VAX */
 
 #endif /* not VAX */
 

+#define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
+

 #else
 /* REAL_ARITHMETIC not defined means that the *host's* data
    structure will be used.  It may differ by endian-ness from the
    target machine's structure and will get its ends swapped
    accordingly (but not here).  Probably only the decimal <-> binary
    functions in this file will actually be used in this case.  */
 #else
 /* REAL_ARITHMETIC not defined means that the *host's* data
    structure will be used.  It may differ by endian-ness from the
    target machine's structure and will get its ends swapped
    accordingly (but not here).  Probably only the decimal <-> binary
    functions in this file will actually be used in this case.  */

+

 #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
 #define DEC 1
 #else /* it's not VAX */
 #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
 #define DEC 1
 #else /* it's not VAX */

+#if HOST_FLOAT_FORMAT == IBM_FLOAT_FORMAT
+/* IBM System/370 style */
+#define IBM 1
+#else /* it's also not an IBM */

 #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
 #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT

-#ifdef HOST_WORDS_BIG_ENDIAN
-#define MIEEE 1
-#else /* not big-endian */
-#define IBMPC 1
-#endif /*  big-endian */
+#define IEEE

 #else /* it's not IEEE either */
 unknown arithmetic type
 #define UNK 1
 #endif /* not IEEE */
 #else /* it's not IEEE either */
 unknown arithmetic type
 #define UNK 1
 #endif /* not IEEE */

+#endif /* not IBM */

 #endif /* not VAX */
 
 #endif /* not VAX */
 

+#define REAL_WORDS_BIG_ENDIAN HOST_FLOAT_WORDS_BIG_ENDIAN
+

 #endif /* REAL_ARITHMETIC not defined */
 
 #endif /* REAL_ARITHMETIC not defined */
 

-/* Define for support of infinity.  */
-#ifndef DEC
+/* Define INFINITY for support of infinity.
+   Define NANS for support of Not-a-Number's (NaN's).  */
+#if !defined(DEC) && !defined(IBM)

 #define INFINITY
 #define INFINITY

+#define NANS

 #endif
 
 #endif
 

-
-/* ehead.h
- *
- * Include file for extended precision arithmetic programs.
- */
-
-/* Number of 16 bit words in external e type format */
-#define NE 6
-
-/* Number of 16 bit words in internal format */
-#define NI (NE+3)
-
-/* Array offset to exponent */
-#define E 1
-
-/* Array offset to high guard word */
-#define M 2
-
-/* Number of bits of precision */
-#define NBITS ((NI-4)*16)
-
-/* Maximum number of decimal digits in ASCII conversion
- * = NBITS*log10(2)
- */
-#define NDEC (NBITS*8/27)
-
-/* The exponent of 1.0 */
-#define EXONE (0x3fff)
-
+/* Support of NaNs requires support of infinity.  */
+#ifdef NANS
+#ifndef INFINITY
+#define INFINITY
+#endif
+#endif
+\f

 /* Find a host integer type that is at least 16 bits wide,
 /* Find a host integer type that is at least 16 bits wide,

-   and another type at least twice whatever that size is. */
+   and another type at least twice whatever that size is.  */

 
 #if HOST_BITS_PER_CHAR >= 16
 #define EMUSHORT char
 
 #if HOST_BITS_PER_CHAR >= 16
 #define EMUSHORT char


@@ -171,7 +193,7 @@ unknown arithmetic type
 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
 #else
 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
 #else

-/*  You will have to modify this program to have a smaller unit size. */
+/*  You will have to modify this program to have a smaller unit size.  */

 #define EMU_NON_COMPILE
 #endif
 #endif
 #define EMU_NON_COMPILE
 #endif
 #endif


@@ -190,7 +212,7 @@ unknown arithmetic type
 #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
 #define EMULONG long long int
 #else
 #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
 #define EMULONG long long int
 #else

-/*  You will have to modify this program to have a smaller unit size. */
+/*  You will have to modify this program to have a smaller unit size.  */

 #define EMU_NON_COMPILE
 #endif
 #endif
 #define EMU_NON_COMPILE
 #endif
 #endif


@@ -198,12 +220,12 @@ unknown arithmetic type
 #endif
 
 
 #endif
 
 

-/* The host interface doesn't work if no 16-bit size exists. */
+/* The host interface doesn't work if no 16-bit size exists.  */

 #if EMUSHORT_SIZE != 16
 #define EMU_NON_COMPILE
 #endif
 
 #if EMUSHORT_SIZE != 16
 #define EMU_NON_COMPILE
 #endif
 

-/* OK to continue compilation. */
+/* OK to continue compilation.  */

 #ifndef EMU_NON_COMPILE
 
 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
 #ifndef EMU_NON_COMPILE
 
 /* Construct macros to translate between REAL_VALUE_TYPE and e type.


@@ -212,66 +234,199 @@ unknown arithmetic type
    in memory, with no holes.  */
 
 #if LONG_DOUBLE_TYPE_SIZE == 96
    in memory, with no holes.  */
 
 #if LONG_DOUBLE_TYPE_SIZE == 96

-#define GET_REAL(r,e) bcopy (r, e, 2*NE)
-#define PUT_REAL(e,r) bcopy (e, r, 2*NE)
+/* Number of 16 bit words in external e type format */
+#define NE 6
+#define MAXDECEXP 4932
+#define MINDECEXP -4956
+#define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
+#define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)

 #else /* no XFmode */
 #else /* no XFmode */

-
+#if LONG_DOUBLE_TYPE_SIZE == 128
+#define NE 10
+#define MAXDECEXP 4932
+#define MINDECEXP -4977
+#define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
+#define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
+#else
+#define NE 6
+#define MAXDECEXP 4932
+#define MINDECEXP -4956

 #ifdef REAL_ARITHMETIC
 /* Emulator uses target format internally
 #ifdef REAL_ARITHMETIC
 /* Emulator uses target format internally

-   but host stores it in host endian-ness. */
-
-#if defined (HOST_WORDS_BIG_ENDIAN) == WORDS_BIG_ENDIAN
-#define GET_REAL(r,e) e53toe ((r), (e))
-#define PUT_REAL(e,r) etoe53 ((e), (r))
-
-#else /* endian-ness differs */
-/* emulator uses target endian-ness internally */
-#define GET_REAL(r,e)          \
-do { EMUSHORT w[4];            \
- w[3] = ((EMUSHORT *) r)[0];   \
- w[2] = ((EMUSHORT *) r)[1];   \
- w[1] = ((EMUSHORT *) r)[2];   \
- w[0] = ((EMUSHORT *) r)[3];   \
- e53toe (w, (e)); } while (0)
-
-#define PUT_REAL(e,r)          \
-do { EMUSHORT w[4];            \
- etoe53 ((e), w);              \
- *((EMUSHORT *) r) = w[3];     \
- *((EMUSHORT *) r + 1) = w[2]; \
- *((EMUSHORT *) r + 2) = w[1]; \
- *((EMUSHORT *) r + 3) = w[0]; } while (0)
-
-#endif /* endian-ness differs */
+   but host stores it in host endian-ness.  */
+
+#define GET_REAL(r,e)                                          \
+do {                                                           \
+     if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
+       e53toe ((unsigned EMUSHORT *) (r), (e));                        \
+     else                                                      \
+       {                                                       \
+        unsigned EMUSHORT w[4];                                \
+        w[3] = ((EMUSHORT *) r)[0];                            \
+        w[2] = ((EMUSHORT *) r)[1];                            \
+        w[1] = ((EMUSHORT *) r)[2];                            \
+        w[0] = ((EMUSHORT *) r)[3];                            \
+        e53toe (w, (e));                                       \
+       }                                                       \
+   } while (0)
+
+#define PUT_REAL(e,r)                                          \
+do {                                                           \
+     if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
+       etoe53 ((e), (unsigned EMUSHORT *) (r));                        \
+     else                                                      \
+       {                                                       \
+        unsigned EMUSHORT w[4];                                \
+        etoe53 ((e), w);                                       \
+        *((EMUSHORT *) r) = w[3];                              \
+        *((EMUSHORT *) r + 1) = w[2];                          \
+        *((EMUSHORT *) r + 2) = w[1];                          \
+        *((EMUSHORT *) r + 3) = w[0];                          \
+       }                                                       \
+   } while (0)

 
 #else /* not REAL_ARITHMETIC */
 
 /* emulator uses host format */
 
 #else /* not REAL_ARITHMETIC */
 
 /* emulator uses host format */

-#define GET_REAL(r,e) e53toe ((r), (e))
-#define PUT_REAL(e,r) etoe53 ((e), (r))
+#define GET_REAL(r,e) e53toe ((unsigned EMUSHORT *) (r), (e))
+#define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))

 
 #endif /* not REAL_ARITHMETIC */
 
 #endif /* not REAL_ARITHMETIC */

+#endif /* not TFmode */

 #endif /* no XFmode */
 
 #endif /* no XFmode */
 

-void warning ();
+
+/* Number of 16 bit words in internal format */
+#define NI (NE+3)
+
+/* Array offset to exponent */
+#define E 1
+
+/* Array offset to high guard word */
+#define M 2
+
+/* Number of bits of precision */
+#define NBITS ((NI-4)*16)
+
+/* Maximum number of decimal digits in ASCII conversion
+ * = NBITS*log10(2)
+ */
+#define NDEC (NBITS*8/27)
+
+/* The exponent of 1.0 */
+#define EXONE (0x3fff)
+

 extern int extra_warnings;
 extern int extra_warnings;

-int ecmp (), enormlz (), eshift (), eisneg (), eisinf ();
-void eadd (), esub (), emul (), ediv ();
-void eshup1 (), eshup8 (), eshup6 (), eshdn1 (), eshdn8 (), eshdn6 ();
-void eabs (), eneg (), emov (), eclear (), einfin (), efloor ();
-void eldexp (), efrexp (), eifrac (), euifrac (), ltoe (), ultoe ();
-void eround (), ereal_to_decimal ();
-void esqrt (), elog (), eexp (), etanh (), epow ();
-void asctoe (), asctoe24 (), asctoe53 (), asctoe64 ();
-void etoasc (), e24toasc (), e53toasc (), e64toasc ();
-void etoe64 (), etoe53 (), etoe24 (), e64toe (), e53toe (), e24toe ();
-void mtherr ();

 extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
 extern unsigned EMUSHORT elog2[], esqrt2[];
 
 extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
 extern unsigned EMUSHORT elog2[], esqrt2[];
 

-/* Pack output array with 32-bit numbers obtained from
-   array containing 16-bit numbers, swapping ends if required. */
-void 
+static void endian     PROTO((unsigned EMUSHORT *, long *,
+                              enum machine_mode));
+static void eclear     PROTO((unsigned EMUSHORT *));
+static void emov       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void eabs       PROTO((unsigned EMUSHORT *));
+static void eneg       PROTO((unsigned EMUSHORT *));
+static int eisneg      PROTO((unsigned EMUSHORT *));
+static int eisinf      PROTO((unsigned EMUSHORT *));
+static int eisnan      PROTO((unsigned EMUSHORT *));
+static void einfin     PROTO((unsigned EMUSHORT *));
+static void enan       PROTO((unsigned EMUSHORT *, int));
+static void emovi      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void emovo      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ecleaz     PROTO((unsigned EMUSHORT *));
+static void ecleazs    PROTO((unsigned EMUSHORT *));
+static void emovz      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void einan      PROTO((unsigned EMUSHORT *));
+static int eiisnan     PROTO((unsigned EMUSHORT *));
+static int eiisneg     PROTO((unsigned EMUSHORT *));
+static void eiinfin    PROTO((unsigned EMUSHORT *));
+static int eiisinf     PROTO((unsigned EMUSHORT *));
+static int ecmpm       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void eshdn1     PROTO((unsigned EMUSHORT *));
+static void eshup1     PROTO((unsigned EMUSHORT *));
+static void eshdn8     PROTO((unsigned EMUSHORT *));
+static void eshup8     PROTO((unsigned EMUSHORT *));
+static void eshup6     PROTO((unsigned EMUSHORT *));
+static void eshdn6     PROTO((unsigned EMUSHORT *));
+static void eaddm      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));\f
+static void esubm      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void m16m       PROTO((unsigned int, unsigned short *,
+                              unsigned short *));
+static int edivm       PROTO((unsigned short *, unsigned short *));
+static int emulm       PROTO((unsigned short *, unsigned short *));
+static void emdnorm    PROTO((unsigned EMUSHORT *, int, int, EMULONG, int));
+static void esub       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+                              unsigned EMUSHORT *));
+static void eadd       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+                              unsigned EMUSHORT *));
+static void eadd1      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+                              unsigned EMUSHORT *));
+static void ediv       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+                              unsigned EMUSHORT *));
+static void emul       PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+                              unsigned EMUSHORT *));
+static void e53toe     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void e64toe     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void e113toe    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void e24toe     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe113    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe113     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe64     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe64      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe53     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe53      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe24     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe24      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static int ecmp                PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void eround     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ltoe       PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *));
+static void ultoe      PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *));
+static void eifrac     PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *,
+                              unsigned EMUSHORT *));
+static void euifrac    PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *,
+                              unsigned EMUSHORT *));
+static int eshift      PROTO((unsigned EMUSHORT *, int));
+static int enormlz     PROTO((unsigned EMUSHORT *));
+static void e24toasc   PROTO((unsigned EMUSHORT *, char *, int));
+static void e53toasc   PROTO((unsigned EMUSHORT *, char *, int));
+static void e64toasc   PROTO((unsigned EMUSHORT *, char *, int));
+static void e113toasc  PROTO((unsigned EMUSHORT *, char *, int));
+static void etoasc     PROTO((unsigned EMUSHORT *, char *, int));
+static void asctoe24   PROTO((char *, unsigned EMUSHORT *));
+static void asctoe53   PROTO((char *, unsigned EMUSHORT *));
+static void asctoe64   PROTO((char *, unsigned EMUSHORT *));
+static void asctoe113  PROTO((char *, unsigned EMUSHORT *));
+static void asctoe     PROTO((char *, unsigned EMUSHORT *));
+static void asctoeg    PROTO((char *, unsigned EMUSHORT *, int));
+static void efloor     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void efrexp     PROTO((unsigned EMUSHORT *, int *,
+                              unsigned EMUSHORT *));
+static void eldexp     PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *));
+static void eremain    PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+                              unsigned EMUSHORT *));
+static void eiremain   PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void mtherr     PROTO((char *, int));
+static void dectoe     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etodec     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void todec      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ibmtoe     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+                              enum machine_mode));
+static void etoibm     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+                              enum machine_mode));
+static void toibm      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+                              enum machine_mode));
+static void make_nan   PROTO((unsigned EMUSHORT *, int, enum machine_mode));
+static void uditoe     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ditoe      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoudi     PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etodi      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void esqrt      PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+\f
+/* Copy 32-bit numbers obtained from array containing 16-bit numbers,
+   swapping ends if required, into output array of longs.  The
+   result is normally passed to fprintf by the ASM_OUTPUT_ macros.   */
+
+static void 

 endian (e, x, mode)
      unsigned EMUSHORT e[];
      long x[];
 endian (e, x, mode)
      unsigned EMUSHORT e[];
      long x[];


@@ -279,87 +434,104 @@ endian (e, x, mode)
 {
   unsigned long th, t;
 
 {
   unsigned long th, t;
 

-#if WORDS_BIG_ENDIAN
-  switch (mode)
+  if (REAL_WORDS_BIG_ENDIAN)

     {
     {

+      switch (mode)
+       {

 
 

-    case XFmode:
-
-      /* Swap halfwords in the third long. */
-      th = (unsigned long) e[4] & 0xffff;
-      t = (unsigned long) e[5] & 0xffff;
-      t |= th << 16;
-      x[2] = (long) t;
-      /* fall into the double case */
-
-    case DFmode:
-
-      /* swap halfwords in the second word */
-      th = (unsigned long) e[2] & 0xffff;
-      t = (unsigned long) e[3] & 0xffff;
-      t |= th << 16;
-      x[1] = (long) t;
-      /* fall into the float case */
-
-    case SFmode:
-
-      /* swap halfwords in the first word */
-      th = (unsigned long) e[0] & 0xffff;
-      t = (unsigned long) e[1] & 0xffff;
-      t |= th << 16;
-      x[0] = t;
-      break;
+       case TFmode:
+         /* Swap halfwords in the fourth long.  */
+         th = (unsigned long) e[6] & 0xffff;
+         t = (unsigned long) e[7] & 0xffff;
+         t |= th << 16;
+         x[3] = (long) t;
+
+       case XFmode:
+
+         /* Swap halfwords in the third long.  */
+         th = (unsigned long) e[4] & 0xffff;
+         t = (unsigned long) e[5] & 0xffff;
+         t |= th << 16;
+         x[2] = (long) t;
+         /* fall into the double case */
+
+       case DFmode:
+
+         /* swap halfwords in the second word */
+         th = (unsigned long) e[2] & 0xffff;
+         t = (unsigned long) e[3] & 0xffff;
+         t |= th << 16;
+         x[1] = (long) t;
+         /* fall into the float case */
+
+       case HFmode:
+       case SFmode:
+
+         /* swap halfwords in the first word */
+         th = (unsigned long) e[0] & 0xffff;
+         t = (unsigned long) e[1] & 0xffff;
+         t |= th << 16;
+         x[0] = (long) t;
+         break;

 
 

-    default:
-      abort ();
+       default:
+         abort ();
+       }

     }
     }

-
-#else
-
-  /* Pack the output array without swapping. */
-
-  switch (mode)
+  else

     {
     {

+      /* Pack the output array without swapping.  */

 
 

-    case XFmode:
-
-      /* Pack the third long.
-        Each element of the input REAL_VALUE_TYPE array has 16 bit useful bits
-        in it.  */
-      th = (unsigned long) e[5] & 0xffff;
-      t = (unsigned long) e[4] & 0xffff;
-      t |= th << 16;
-      x[2] = (long) t;
-      /* fall into the double case */
-
-    case DFmode:
-
-      /* pack the second long */
-      th = (unsigned long) e[3] & 0xffff;
-      t = (unsigned long) e[2] & 0xffff;
-      t |= th << 16;
-      x[1] = (long) t;
-      /* fall into the float case */
-
-    case SFmode:
+      switch (mode)
+       {

 
 

-      /* pack the first long */
-      th = (unsigned long) e[1] & 0xffff;
-      t = (unsigned long) e[0] & 0xffff;
-      t |= th << 16;
-      x[0] = t;
-      break;
+       case TFmode:
+
+         /* Pack the fourth long.  */
+         th = (unsigned long) e[7] & 0xffff;
+         t = (unsigned long) e[6] & 0xffff;
+         t |= th << 16;
+         x[3] = (long) t;
+
+       case XFmode:
+
+         /* Pack the third long.
+            Each element of the input REAL_VALUE_TYPE array has 16 useful bits
+            in it.  */
+         th = (unsigned long) e[5] & 0xffff;
+         t = (unsigned long) e[4] & 0xffff;
+         t |= th << 16;
+         x[2] = (long) t;
+         /* fall into the double case */
+
+       case DFmode:
+
+         /* pack the second long */
+         th = (unsigned long) e[3] & 0xffff;
+         t = (unsigned long) e[2] & 0xffff;
+         t |= th << 16;
+         x[1] = (long) t;
+         /* fall into the float case */
+
+       case HFmode:
+       case SFmode:
+
+         /* pack the first long */
+         th = (unsigned long) e[1] & 0xffff;
+         t = (unsigned long) e[0] & 0xffff;
+         t |= th << 16;
+         x[0] = (long) t;
+         break;

 
 

-    default:
-      abort ();
+       default:
+         abort ();
+       }

     }
     }

-
-#endif

 }
 
 
 }
 
 

-/* This is the implementation of the REAL_ARITHMETIC macro.
- */
+/* This is the implementation of the REAL_ARITHMETIC macro.  */
+

 void 
 earith (value, icode, r1, r2)
      REAL_VALUE_TYPE *value;
 void 
 earith (value, icode, r1, r2)
      REAL_VALUE_TYPE *value;


@@ -372,6 +544,19 @@ earith (value, icode, r1, r2)
 
   GET_REAL (r1, d1);
   GET_REAL (r2, d2);
 
   GET_REAL (r1, d1);
   GET_REAL (r2, d2);

+#ifdef NANS
+/*  Return NaN input back to the caller.  */
+  if (eisnan (d1))
+    {
+      PUT_REAL (d1, value);
+      return;
+    }
+  if (eisnan (d2))
+    {
+      PUT_REAL (d2, value);
+      return;
+    }
+#endif

   code = (enum tree_code) icode;
   switch (code)
     {
   code = (enum tree_code) icode;
   switch (code)
     {


@@ -390,8 +575,15 @@ earith (value, icode, r1, r2)
     case RDIV_EXPR:
 #ifndef REAL_INFINITY
       if (ecmp (d2, ezero) == 0)
     case RDIV_EXPR:
 #ifndef REAL_INFINITY
       if (ecmp (d2, ezero) == 0)

+       {
+#ifdef NANS
+       enan (v, eisneg (d1) ^ eisneg (d2));
+       break;
+#else

        abort ();
 #endif
        abort ();
 #endif

+       }
+#endif

       ediv (d2, d1, v);        /* d1/d2 */
       break;
 
       ediv (d2, d1, v);        /* d1/d2 */
       break;
 


@@ -416,18 +608,22 @@ PUT_REAL (v, value);
 }
 
 
 }
 
 

-/* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT
- * implements REAL_VALUE_FIX_TRUNCATE (x) (etrunci (x))
- */
+/* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
+   implements REAL_VALUE_RNDZINT (x) (etrunci (x)).  */
+

 REAL_VALUE_TYPE 
 etrunci (x)
      REAL_VALUE_TYPE x;
 {
   unsigned EMUSHORT f[NE], g[NE];
   REAL_VALUE_TYPE r;
 REAL_VALUE_TYPE 
 etrunci (x)
      REAL_VALUE_TYPE x;
 {
   unsigned EMUSHORT f[NE], g[NE];
   REAL_VALUE_TYPE r;

-  long l;
+  HOST_WIDE_INT l;

 
   GET_REAL (&x, g);
 
   GET_REAL (&x, g);

+#ifdef NANS
+  if (eisnan (g))
+    return (x);
+#endif

   eifrac (g, &l, f);
   ltoe (&l, g);
   PUT_REAL (g, &r);
   eifrac (g, &l, f);
   ltoe (&l, g);
   PUT_REAL (g, &r);


@@ -435,18 +631,22 @@ etrunci (x)
 }
 
 
 }
 
 

-/* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT
- * implements REAL_VALUE_UNSIGNED_FIX_TRUNCATE (x) (etruncui (x))
- */
+/* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
+   implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)).  */
+

 REAL_VALUE_TYPE 
 etruncui (x)
      REAL_VALUE_TYPE x;
 {
   unsigned EMUSHORT f[NE], g[NE];
   REAL_VALUE_TYPE r;
 REAL_VALUE_TYPE 
 etruncui (x)
      REAL_VALUE_TYPE x;
 {
   unsigned EMUSHORT f[NE], g[NE];
   REAL_VALUE_TYPE r;

-  unsigned long l;
+  unsigned HOST_WIDE_INT l;

 
   GET_REAL (&x, g);
 
   GET_REAL (&x, g);

+#ifdef NANS
+  if (eisnan (g))
+    return (x);
+#endif

   euifrac (g, &l, f);
   ultoe (&l, g);
   PUT_REAL (g, &r);
   euifrac (g, &l, f);
   ultoe (&l, g);
   PUT_REAL (g, &r);


@@ -454,11 +654,10 @@ etruncui (x)
 }
 
 
 }
 
 

-/* This is the REAL_VALUE_ATOF function.
- * It converts a decimal string to binary, rounding off
- * as indicated by the machine_mode argument.  Then it
- * promotes the rounded value to REAL_VALUE_TYPE.
- */
+/* This is the REAL_VALUE_ATOF function.  It converts a decimal string to
+   binary, rounding off as indicated by the machine_mode argument.  Then it
+   promotes the rounded value to REAL_VALUE_TYPE.  */
+

 REAL_VALUE_TYPE 
 ereal_atof (s, t)
      char *s;
 REAL_VALUE_TYPE 
 ereal_atof (s, t)
      char *s;


@@ -469,6 +668,7 @@ ereal_atof (s, t)
 
   switch (t)
     {
 
   switch (t)
     {

+    case HFmode:

     case SFmode:
       asctoe24 (s, tem);
       e24toe (tem, e);
     case SFmode:
       asctoe24 (s, tem);
       e24toe (tem, e);


@@ -481,6 +681,10 @@ ereal_atof (s, t)
       asctoe64 (s, tem);
       e64toe (tem, e);
       break;
       asctoe64 (s, tem);
       e64toe (tem, e);
       break;

+    case TFmode:
+      asctoe113 (s, tem);
+      e113toe (tem, e);
+      break;

     default:
       asctoe (s, e);
     }
     default:
       asctoe (s, e);
     }


@@ -489,8 +693,8 @@ ereal_atof (s, t)
 }
 
 
 }
 
 

-/* Expansion of REAL_NEGATE.
- */
+/* Expansion of REAL_NEGATE.  */
+

 REAL_VALUE_TYPE 
 ereal_negate (x)
      REAL_VALUE_TYPE x;
 REAL_VALUE_TYPE 
 ereal_negate (x)
      REAL_VALUE_TYPE x;


@@ -505,55 +709,66 @@ ereal_negate (x)
 }
 
 
 }
 
 

-/* Round real to int
- * implements REAL_VALUE_FIX (x) (eroundi (x))
- * The type of rounding is left unspecified by real.h.
- * It is implemented here as round to nearest (add .5 and chop).
- */
-int 
-eroundi (x)
+/* Round real toward zero to HOST_WIDE_INT;
+   implements REAL_VALUE_FIX (x).  */
+
+HOST_WIDE_INT
+efixi (x)

      REAL_VALUE_TYPE x;
 {
   unsigned EMUSHORT f[NE], g[NE];
      REAL_VALUE_TYPE x;
 {
   unsigned EMUSHORT f[NE], g[NE];

-  EMULONG l;
+  HOST_WIDE_INT l;

 
   GET_REAL (&x, f);
 
   GET_REAL (&x, f);

-  eround (f, g);
-  eifrac (g, &l, f);
-  return ((int) l);
+#ifdef NANS
+  if (eisnan (f))
+    {
+      warning ("conversion from NaN to int");
+      return (-1);
+    }
+#endif
+  eifrac (f, &l, g);
+  return l;

 }
 
 }
 

-/* Round real to nearest unsigned int
- * implements  REAL_VALUE_UNSIGNED_FIX (x) ((unsigned int) eroundi (x))
- * Negative input returns zero.
- * The type of rounding is left unspecified by real.h.
- * It is implemented here as round to nearest (add .5 and chop).
- */
-unsigned int 
-eroundui (x)
+/* Round real toward zero to unsigned HOST_WIDE_INT
+   implements  REAL_VALUE_UNSIGNED_FIX (x).
+   Negative input returns zero.  */
+
+unsigned HOST_WIDE_INT
+efixui (x)

      REAL_VALUE_TYPE x;
 {
   unsigned EMUSHORT f[NE], g[NE];
      REAL_VALUE_TYPE x;
 {
   unsigned EMUSHORT f[NE], g[NE];

-  unsigned EMULONG l;
+  unsigned HOST_WIDE_INT l;

 
   GET_REAL (&x, f);
 
   GET_REAL (&x, f);

-  eround (f, g);
-  euifrac (g, &l, f);
-  return ((unsigned int)l);
+#ifdef NANS
+  if (eisnan (f))
+    {
+      warning ("conversion from NaN to unsigned int");
+      return (-1);
+    }
+#endif
+  euifrac (f, &l, g);
+  return l;

 }
 
 
 }
 
 

-/* REAL_VALUE_FROM_INT macro.
- */
+/* REAL_VALUE_FROM_INT macro.  */
+

 void 
 void 

-ereal_from_int (d, i, j)
+ereal_from_int (d, i, j, mode)

      REAL_VALUE_TYPE *d;
      REAL_VALUE_TYPE *d;

-     long i, j;
+     HOST_WIDE_INT i, j;
+     enum machine_mode mode;

 {
   unsigned EMUSHORT df[NE], dg[NE];
 {
   unsigned EMUSHORT df[NE], dg[NE];

-  long low, high;
+  HOST_WIDE_INT low, high;

   int sign;
 
   int sign;
 

+  if (GET_MODE_CLASS (mode) != MODE_FLOAT)
+    abort ();

   sign = 0;
   low = i;
   if ((high = j) < 0)
   sign = 0;
   low = i;
   if ((high = j) < 0)


@@ -566,49 +781,121 @@ ereal_from_int (d, i, j)
       else
        high += 1;
     }
       else
        high += 1;
     }

-  eldexp (eone, HOST_BITS_PER_LONG, df);
-  ultoe (&high, dg);
+  eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
+  ultoe ((unsigned HOST_WIDE_INT *) &high, dg);

   emul (dg, df, dg);
   emul (dg, df, dg);

-  ultoe (&low, df);
+  ultoe ((unsigned HOST_WIDE_INT *) &low, df);

   eadd (df, dg, dg);
   if (sign)
     eneg (dg);
   eadd (df, dg, dg);
   if (sign)
     eneg (dg);

+
+  /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
+     Avoid double-rounding errors later by rounding off now from the
+     extra-wide internal format to the requested precision.  */
+  switch (GET_MODE_BITSIZE (mode))
+    {
+    case 32:
+      etoe24 (dg, df);
+      e24toe (df, dg);
+      break;
+
+    case 64:
+      etoe53 (dg, df);
+      e53toe (df, dg);
+      break;
+
+    case 96:
+      etoe64 (dg, df);
+      e64toe (df, dg);
+      break;
+
+    case 128:
+      etoe113 (dg, df);
+      e113toe (df, dg);
+      break;
+
+    default:
+      abort ();
+  }
+

   PUT_REAL (dg, d);
 }
 
 
   PUT_REAL (dg, d);
 }
 
 

-/* REAL_VALUE_FROM_UNSIGNED_INT macro.
- */
+/* REAL_VALUE_FROM_UNSIGNED_INT macro.   */
+

 void 
 void 

-ereal_from_uint (d, i, j)
+ereal_from_uint (d, i, j, mode)

      REAL_VALUE_TYPE *d;
      REAL_VALUE_TYPE *d;

-     unsigned long i, j;
+     unsigned HOST_WIDE_INT i, j;
+     enum machine_mode mode;

 {
   unsigned EMUSHORT df[NE], dg[NE];
 {
   unsigned EMUSHORT df[NE], dg[NE];

-  unsigned long low, high;
+  unsigned HOST_WIDE_INT low, high;

 
 

+  if (GET_MODE_CLASS (mode) != MODE_FLOAT)
+    abort ();

   low = i;
   high = j;
   low = i;
   high = j;

-  eldexp (eone, HOST_BITS_PER_LONG, df);
+  eldexp (eone, HOST_BITS_PER_WIDE_INT, df);

   ultoe (&high, dg);
   emul (dg, df, dg);
   ultoe (&low, df);
   eadd (df, dg, dg);
   ultoe (&high, dg);
   emul (dg, df, dg);
   ultoe (&low, df);
   eadd (df, dg, dg);

+
+  /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
+     Avoid double-rounding errors later by rounding off now from the
+     extra-wide internal format to the requested precision.  */
+  switch (GET_MODE_BITSIZE (mode))
+    {
+    case 32:
+      etoe24 (dg, df);
+      e24toe (df, dg);
+      break;
+
+    case 64:
+      etoe53 (dg, df);
+      e53toe (df, dg);
+      break;
+
+    case 96:
+      etoe64 (dg, df);
+      e64toe (df, dg);
+      break;
+
+    case 128:
+      etoe113 (dg, df);
+      e113toe (df, dg);
+      break;
+
+    default:
+      abort ();
+  }
+

   PUT_REAL (dg, d);
 }
 
 
   PUT_REAL (dg, d);
 }
 
 

-/* REAL_VALUE_TO_INT macro
- */
+/* REAL_VALUE_TO_INT macro.  */
+

 void 
 ereal_to_int (low, high, rr)
 void 
 ereal_to_int (low, high, rr)

-     long *low, *high;
+     HOST_WIDE_INT *low, *high;

      REAL_VALUE_TYPE rr;
 {
   unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
   int s;
 
   GET_REAL (&rr, d);
      REAL_VALUE_TYPE rr;
 {
   unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
   int s;
 
   GET_REAL (&rr, d);

+#ifdef NANS
+  if (eisnan (d))
+    {
+      warning ("conversion from NaN to int");
+      *low = -1;
+      *high = -1;
+      return;
+    }
+#endif

   /* convert positive value */
   s = 0;
   if (eisneg (d))
   /* convert positive value */
   s = 0;
   if (eisneg (d))


@@ -616,11 +903,11 @@ ereal_to_int (low, high, rr)
       eneg (d);
       s = 1;
     }
       eneg (d);
       s = 1;
     }

-  eldexp (eone, HOST_BITS_PER_LONG, df);
+  eldexp (eone, HOST_BITS_PER_WIDE_INT, df);

   ediv (df, d, dg);            /* dg = d / 2^32 is the high word */
   ediv (df, d, dg);            /* dg = d / 2^32 is the high word */

-  euifrac (dg, high, dh);
+  euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);

   emul (df, dh, dg);           /* fractional part is the low word */
   emul (df, dh, dg);           /* fractional part is the low word */

-  euifrac (dg, low, dh);
+  euifrac (dg, (unsigned HOST_WIDE_INT *)low, dh);

   if (s)
     {
       /* complement and add 1 */
   if (s)
     {
       /* complement and add 1 */


@@ -633,8 +920,8 @@ ereal_to_int (low, high, rr)
 }
 
 
 }
 
 

-/* REAL_VALUE_LDEXP macro.
- */
+/* REAL_VALUE_LDEXP macro.  */
+

 REAL_VALUE_TYPE
 ereal_ldexp (x, n)
      REAL_VALUE_TYPE x;
 REAL_VALUE_TYPE
 ereal_ldexp (x, n)
      REAL_VALUE_TYPE x;


@@ -644,16 +931,22 @@ ereal_ldexp (x, n)
   REAL_VALUE_TYPE r;
 
   GET_REAL (&x, e);
   REAL_VALUE_TYPE r;
 
   GET_REAL (&x, e);

+#ifdef NANS
+  if (eisnan (e))
+    return (x);
+#endif

   eldexp (e, n, y);
   PUT_REAL (y, &r);
   return (r);
 }
 
 /* These routines are conditionally compiled because functions
   eldexp (e, n, y);
   PUT_REAL (y, &r);
   return (r);
 }
 
 /* These routines are conditionally compiled because functions

- * of the same names may be defined in fold-const.c.  */
+   of the same names may be defined in fold-const.c.  */
+

 #ifdef REAL_ARITHMETIC
 
 #ifdef REAL_ARITHMETIC
 

-/* Check for infinity in a REAL_VALUE_TYPE. */
+/* Check for infinity in a REAL_VALUE_TYPE.  */
+

 int
 target_isinf (x)
      REAL_VALUE_TYPE x;
 int
 target_isinf (x)
      REAL_VALUE_TYPE x;


@@ -668,36 +961,36 @@ target_isinf (x)
 #endif
 }
 
 #endif
 }
 

-
-/* Check whether an IEEE double precision number is a NaN. */
+/* Check whether a REAL_VALUE_TYPE item is a NaN.  */

 
 int
 target_isnan (x)
      REAL_VALUE_TYPE x;
 {
 
 int
 target_isnan (x)
      REAL_VALUE_TYPE x;
 {

+  unsigned EMUSHORT e[NE];
+
+#ifdef NANS
+  GET_REAL (&x, e);
+  return (eisnan (e));
+#else

   return (0);
   return (0);

+#endif

 }
 
 
 }
 
 

-/* Check for a negative IEEE double precision number.
- * this means strictly less than zero, not -0.
- */
+/* Check for a negative REAL_VALUE_TYPE number.
+   This just checks the sign bit, so that -0 counts as negative.  */

 
 int
 target_negative (x)
      REAL_VALUE_TYPE x;
 {
 
 int
 target_negative (x)
      REAL_VALUE_TYPE x;
 {

-  unsigned EMUSHORT e[NE];
-
-  GET_REAL (&x, e);
-  if (ecmp (e, ezero) < 0)
-    return (1);
-  return (0);
+  return ereal_isneg (x);

 }
 
 /* Expansion of REAL_VALUE_TRUNCATE.
 }
 
 /* Expansion of REAL_VALUE_TRUNCATE.

- * The result is in floating point, rounded to nearest or even.
- */
+   The result is in floating point, rounded to nearest or even.  */
+

 REAL_VALUE_TYPE
 real_value_truncate (mode, arg)
      enum machine_mode mode;
 REAL_VALUE_TYPE
 real_value_truncate (mode, arg)
      enum machine_mode mode;


@@ -707,9 +1000,18 @@ real_value_truncate (mode, arg)
   REAL_VALUE_TYPE r;
 
   GET_REAL (&arg, e);
   REAL_VALUE_TYPE r;
 
   GET_REAL (&arg, e);

+#ifdef NANS
+  if (eisnan (e))
+    return (arg);
+#endif

   eclear (t);
   switch (mode)
     {
   eclear (t);
   switch (mode)
     {

+    case TFmode:
+      etoe113 (e, t);
+      e113toe (t, t);
+      break;
+

     case XFmode:
       etoe64 (e, t);
       e64toe (t, t);
     case XFmode:
       etoe64 (e, t);
       e64toe (t, t);


@@ -720,26 +1022,135 @@ real_value_truncate (mode, arg)
       e53toe (t, t);
       break;
 
       e53toe (t, t);
       break;
 

+    case HFmode:

     case SFmode:
       etoe24 (e, t);
       e24toe (t, t);
       break;
 
     case SImode:
     case SFmode:
       etoe24 (e, t);
       e24toe (t, t);
       break;
 
     case SImode:

-      r = etrunci (e);
+      r = etrunci (arg);

       return (r);
 
       return (r);
 

+    /* If an unsupported type was requested, presume that
+       the machine files know something useful to do with
+       the unmodified value.  */
+

     default:
     default:

-      abort ();
+      return (arg);

     }
   PUT_REAL (t, &r);
   return (r);
 }
 
     }
   PUT_REAL (t, &r);
   return (r);
 }
 

+/* Try to change R into its exact multiplicative inverse in machine mode
+   MODE.  Return nonzero function value if successful.  */
+
+int
+exact_real_inverse (mode, r)
+     enum machine_mode mode;
+     REAL_VALUE_TYPE *r;
+{
+  unsigned EMUSHORT e[NE], einv[NE];
+  REAL_VALUE_TYPE rinv;
+  int i;
+
+  GET_REAL (r, e);
+
+  /* Test for input in range.  Don't transform IEEE special values.  */
+  if (eisinf (e) || eisnan (e) || (ecmp (e, ezero) == 0))
+    return 0;
+
+  /* Test for a power of 2: all significand bits zero except the MSB.
+     We are assuming the target has binary (or hex) arithmetic.  */
+  if (e[NE - 2] != 0x8000)
+    return 0;
+
+  for (i = 0; i < NE - 2; i++)
+    {
+      if (e[i] != 0)
+       return 0;
+    }
+
+  /* Compute the inverse and truncate it to the required mode.  */
+  ediv (e, eone, einv);
+  PUT_REAL (einv, &rinv);
+  rinv = real_value_truncate (mode, rinv);
+
+#ifdef CHECK_FLOAT_VALUE
+  /* This check is not redundant.  It may, for example, flush
+     a supposedly IEEE denormal value to zero.  */
+  i = 0;
+  if (CHECK_FLOAT_VALUE (mode, rinv, i))
+    return 0;
+#endif
+  GET_REAL (&rinv, einv);
+
+  /* Check the bits again, because the truncation might have
+     generated an arbitrary saturation value on overflow.  */
+  if (einv[NE - 2] != 0x8000)
+    return 0;
+
+  for (i = 0; i < NE - 2; i++)
+    {
+      if (einv[i] != 0)
+       return 0;
+    }
+
+  /* Fail if the computed inverse is out of range.  */
+  if (eisinf (einv) || eisnan (einv) || (ecmp (einv, ezero) == 0))
+    return 0;
+
+  /* Output the reciprocal and return success flag.  */
+  PUT_REAL (einv, r);
+  return 1;
+}

 #endif /* REAL_ARITHMETIC defined */
 
 #endif /* REAL_ARITHMETIC defined */
 

-/* Target values are arrays of host longs. A long is guaranteed
-   to be at least 32 bits wide. */
+/* Used for debugging--print the value of R in human-readable format
+   on stderr.  */
+
+void
+debug_real (r)
+     REAL_VALUE_TYPE r;
+{
+  char dstr[30];
+
+  REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
+  fprintf (stderr, "%s", dstr);
+}  
+
+\f
+/* The following routines convert REAL_VALUE_TYPE to the various floating
+   point formats that are meaningful to supported computers.
+
+   The results are returned in 32-bit pieces, each piece stored in a `long'.  
+   This is so they can be printed by statements like
+ 
+      fprintf (file, "%lx, %lx", L[0],  L[1]);
+
+   that will work on both narrow- and wide-word host computers.  */
+
+/* Convert R to a 128-bit long double precision value.  The output array L
+   contains four 32-bit pieces of the result, in the order they would appear
+   in memory.  */
+
+void 
+etartdouble (r, l)
+     REAL_VALUE_TYPE r;
+     long l[];
+{
+  unsigned EMUSHORT e[NE];
+
+  GET_REAL (&r, e);
+  etoe113 (e, e);
+  endian (e, l, TFmode);
+}
+
+/* Convert R to a double extended precision value.  The output array L
+   contains three 32-bit pieces of the result, in the order they would
+   appear in memory.  */
+

 void 
 etarldouble (r, l)
      REAL_VALUE_TYPE r;
 void 
 etarldouble (r, l)
      REAL_VALUE_TYPE r;


@@ -752,6 +1163,9 @@ etarldouble (r, l)
   endian (e, l, XFmode);
 }
 
   endian (e, l, XFmode);
 }
 

+/* Convert R to a double precision value.  The output array L contains two
+   32-bit pieces of the result, in the order they would appear in memory.  */
+

 void 
 etardouble (r, l)
      REAL_VALUE_TYPE r;
 void 
 etardouble (r, l)
      REAL_VALUE_TYPE r;


@@ -764,12 +1178,15 @@ etardouble (r, l)
   endian (e, l, DFmode);
 }
 
   endian (e, l, DFmode);
 }
 

+/* Convert R to a single precision float value stored in the least-significant
+   bits of a `long'.  */
+

 long
 etarsingle (r)
      REAL_VALUE_TYPE r;
 {
   unsigned EMUSHORT e[NE];
 long
 etarsingle (r)
      REAL_VALUE_TYPE r;
 {
   unsigned EMUSHORT e[NE];

-  unsigned long l;
+  long l;

 
   GET_REAL (&r, e);
   etoe24 (e, e);
 
   GET_REAL (&r, e);
   etoe24 (e, e);


@@ -777,6 +1194,11 @@ etarsingle (r)
   return ((long) l);
 }
 
   return ((long) l);
 }
 

+/* Convert X to a decimal ASCII string S for output to an assembly
+   language file.  Note, there is no standard way to spell infinity or
+   a NaN, so these values may require special treatment in the tm.h
+   macros.  */
+

 void
 ereal_to_decimal (x, s)
      REAL_VALUE_TYPE x;
 void
 ereal_to_decimal (x, s)
      REAL_VALUE_TYPE x;


@@ -788,6 +1210,9 @@ ereal_to_decimal (x, s)
   etoasc (e, s, 20);
 }
 
   etoasc (e, s, 20);
 }
 

+/* Compare X and Y.  Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
+   or -2 if either is a NaN.   */
+

 int
 ereal_cmp (x, y)
      REAL_VALUE_TYPE x, y;
 int
 ereal_cmp (x, y)
      REAL_VALUE_TYPE x, y;


@@ -799,6 +1224,8 @@ ereal_cmp (x, y)
   return (ecmp (ex, ey));
 }
 
   return (ecmp (ex, ey));
 }
 

+/*  Return 1 if the sign bit of X is set, else return 0.  */
+

 int
 ereal_isneg (x)
      REAL_VALUE_TYPE x;
 int
 ereal_isneg (x)
      REAL_VALUE_TYPE x;


@@ -810,178 +1237,163 @@ ereal_isneg (x)
 }
 
 /* End of REAL_ARITHMETIC interface */
 }
 
 /* End of REAL_ARITHMETIC interface */

-
-/*                                                     ieee.c
- *
- *    Extended precision IEEE binary floating point arithmetic routines
- *
- * Numbers are stored in C language as arrays of 16-bit unsigned
- * short integers.  The arguments of the routines are pointers to
- * the arrays.
- *
- *
- * External e type data structure, simulates Intel 8087 chip
- * temporary real format but possibly with a larger significand:
- *
- *     NE-1 significand words  (least significant word first,
- *                              most significant bit is normally set)
- *     exponent                (value = EXONE for 1.0,
- *                             top bit is the sign)
- *
- *
- * Internal data structure of a number (a "word" is 16 bits):
- *
- * ei[0]       sign word       (0 for positive, 0xffff for negative)
- * ei[1]       biased exponent (value = EXONE for the number 1.0)
- * ei[2]       high guard word (always zero after normalization)
- * ei[3]
- * to ei[NI-2] significand     (NI-4 significand words,
- *                              most significant word first,
- *                              most significant bit is set)
- * ei[NI-1]    low guard word  (0x8000 bit is rounding place)
- *
- *
- *
- *             Routines for external format numbers
- *
- *     asctoe (string, e)      ASCII string to extended double e type
- *     asctoe64 (string, &d)   ASCII string to long double
- *     asctoe53 (string, &d)   ASCII string to double
- *     asctoe24 (string, &f)   ASCII string to single
- *     asctoeg (string, e, prec) ASCII string to specified precision
- *     e24toe (&f, e)          IEEE single precision to e type
- *     e53toe (&d, e)          IEEE double precision to e type
- *     e64toe (&d, e)          IEEE long double precision to e type
- *     eabs (e)                        absolute value
- *     eadd (a, b, c)          c = b + a
- *     eclear (e)              e = 0
- *     ecmp (a, b)             compare a to b, return 1, 0, or -1
- *     ediv (a, b, c)          c = b / a
- *     efloor (a, b)           truncate to integer, toward -infinity
- *     efrexp (a, exp, s)      extract exponent and significand
- *     eifrac (e, &l, frac)    e to long integer and e type fraction
- *     euifrac (e, &l, frac)   e to unsigned long integer and e type fraction
- *     einfin (e)              set e to infinity, leaving its sign alone
- *     eldexp (a, n, b)        multiply by 2**n
- *     emov (a, b)             b = a
- *     emul (a, b, c)          c = b * a
- *     eneg (e)                        e = -e
- *     eround (a, b)           b = nearest integer value to a
- *     esub (a, b, c)          c = b - a
- *     e24toasc (&f, str, n)   single to ASCII string, n digits after decimal
- *     e53toasc (&d, str, n)   double to ASCII string, n digits after decimal
- *     e64toasc (&d, str, n)   long double to ASCII string
- *     etoasc (e, str, n)      e to ASCII string, n digits after decimal
- *     etoe24 (e, &f)          convert e type to IEEE single precision
- *     etoe53 (e, &d)          convert e type to IEEE double precision
- *     etoe64 (e, &d)          convert e type to IEEE long double precision
- *     ltoe (&l, e)            long (32 bit) integer to e type
- *     ultoe (&l, e)           unsigned long (32 bit) integer to e type
- *      eisneg (e)              1 if sign bit of e != 0, else 0
- *      eisinf (e)              1 if e has maximum exponent
- *
- *
- *             Routines for internal format numbers
- *
- *     eaddm (ai, bi)          add significands, bi = bi + ai
- *     ecleaz (ei)             ei = 0
- *     ecleazs (ei)            set ei = 0 but leave its sign alone
- *     ecmpm (ai, bi)          compare significands, return 1, 0, or -1
- *     edivm (ai, bi)          divide  significands, bi = bi / ai
- *     emdnorm (ai,l,s,exp)    normalize and round off
- *     emovi (a, ai)           convert external a to internal ai
- *     emovo (ai, a)           convert internal ai to external a
- *     emovz (ai, bi)          bi = ai, low guard word of bi = 0
- *     emulm (ai, bi)          multiply significands, bi = bi * ai
- *     enormlz (ei)            left-justify the significand
- *     eshdn1 (ai)             shift significand and guards down 1 bit
- *     eshdn8 (ai)             shift down 8 bits
- *     eshdn6 (ai)             shift down 16 bits
- *     eshift (ai, n)          shift ai n bits up (or down if n < 0)
- *     eshup1 (ai)             shift significand and guards up 1 bit
- *     eshup8 (ai)             shift up 8 bits
- *     eshup6 (ai)             shift up 16 bits
- *     esubm (ai, bi)          subtract significands, bi = bi - ai
- *
- *
- * The result is always normalized and rounded to NI-4 word precision
- * after each arithmetic operation.
- *
- * Exception flags and NaNs are NOT fully supported.
- * This arithmetic should never produce a NaN output, but it might
- * be confused by a NaN input.
- * Define INFINITY in mconf.h for support of infinity; otherwise a
- * saturation arithmetic is implemented.
- * Denormals are always supported here where appropriate (e.g., not
- * for conversion to DEC numbers).
- *
- */
-
+\f

 /*
 /*

- * Revision history:
- *
- *  5 Jan 84   PDP-11 assembly language version
- *  2 Mar 86   fixed bug in asctoq
- *  6 Dec 86   C language version
- * 30 Aug 88   100 digit version, improved rounding
- * 15 May 92    80-bit long double support
- *
- * Author:  S. L. Moshier.
- */
+  Extended precision IEEE binary floating point arithmetic routines
+
+  Numbers are stored in C language as arrays of 16-bit unsigned
+  short integers.  The arguments of the routines are pointers to
+  the arrays.
+
+  External e type data structure, similar to Intel 8087 chip
+  temporary real format but possibly with a larger significand:
+
+       NE-1 significand words  (least significant word first,
+                                most significant bit is normally set)
+       exponent                (value = EXONE for 1.0,
+                               top bit is the sign)
+
+
+  Internal exploded e-type data structure of a number (a "word" is 16 bits):
+
+  ei[0]        sign word       (0 for positive, 0xffff for negative)
+  ei[1]        biased exponent (value = EXONE for the number 1.0)
+  ei[2]        high guard word (always zero after normalization)
+  ei[3]
+  to ei[NI-2]  significand     (NI-4 significand words,
+                                most significant word first,
+                                most significant bit is set)
+  ei[NI-1]     low guard word  (0x8000 bit is rounding place)
+ 
+ 
+ 
+               Routines for external format e-type numbers
+ 
+       asctoe (string, e)      ASCII string to extended double e type
+       asctoe64 (string, &d)   ASCII string to long double
+       asctoe53 (string, &d)   ASCII string to double
+       asctoe24 (string, &f)   ASCII string to single
+       asctoeg (string, e, prec) ASCII string to specified precision
+       e24toe (&f, e)          IEEE single precision to e type
+       e53toe (&d, e)          IEEE double precision to e type
+       e64toe (&d, e)          IEEE long double precision to e type
+       e113toe (&d, e)         128-bit long double precision to e type
+       eabs (e)                        absolute value
+       eadd (a, b, c)          c = b + a
+       eclear (e)              e = 0
+       ecmp (a, b)             Returns 1 if a > b, 0 if a == b,
+                               -1 if a < b, -2 if either a or b is a NaN.
+       ediv (a, b, c)          c = b / a
+       efloor (a, b)           truncate to integer, toward -infinity
+       efrexp (a, exp, s)      extract exponent and significand
+       eifrac (e, &l, frac)    e to HOST_WIDE_INT and e type fraction
+       euifrac (e, &l, frac)   e to unsigned HOST_WIDE_INT and e type fraction
+       einfin (e)              set e to infinity, leaving its sign alone
+       eldexp (a, n, b)        multiply by 2**n
+       emov (a, b)             b = a
+       emul (a, b, c)          c = b * a
+       eneg (e)                        e = -e
+       eround (a, b)           b = nearest integer value to a
+       esub (a, b, c)          c = b - a
+       e24toasc (&f, str, n)   single to ASCII string, n digits after decimal
+       e53toasc (&d, str, n)   double to ASCII string, n digits after decimal
+       e64toasc (&d, str, n)   80-bit long double to ASCII string
+       e113toasc (&d, str, n)  128-bit long double to ASCII string
+       etoasc (e, str, n)      e to ASCII string, n digits after decimal
+       etoe24 (e, &f)          convert e type to IEEE single precision
+       etoe53 (e, &d)          convert e type to IEEE double precision
+       etoe64 (e, &d)          convert e type to IEEE long double precision
+       ltoe (&l, e)            HOST_WIDE_INT to e type
+       ultoe (&l, e)           unsigned HOST_WIDE_INT to e type
+       eisneg (e)              1 if sign bit of e != 0, else 0
+       eisinf (e)              1 if e has maximum exponent (non-IEEE)
+                               or is infinite (IEEE)
+        eisnan (e)              1 if e is a NaN
+ 
+
+               Routines for internal format exploded e-type numbers
+ 
+       eaddm (ai, bi)          add significands, bi = bi + ai
+       ecleaz (ei)             ei = 0
+       ecleazs (ei)            set ei = 0 but leave its sign alone
+       ecmpm (ai, bi)          compare significands, return 1, 0, or -1
+       edivm (ai, bi)          divide  significands, bi = bi / ai
+       emdnorm (ai,l,s,exp)    normalize and round off
+       emovi (a, ai)           convert external a to internal ai
+       emovo (ai, a)           convert internal ai to external a
+       emovz (ai, bi)          bi = ai, low guard word of bi = 0
+       emulm (ai, bi)          multiply significands, bi = bi * ai
+       enormlz (ei)            left-justify the significand
+       eshdn1 (ai)             shift significand and guards down 1 bit
+       eshdn8 (ai)             shift down 8 bits
+       eshdn6 (ai)             shift down 16 bits
+       eshift (ai, n)          shift ai n bits up (or down if n < 0)
+       eshup1 (ai)             shift significand and guards up 1 bit
+       eshup8 (ai)             shift up 8 bits
+       eshup6 (ai)             shift up 16 bits
+       esubm (ai, bi)          subtract significands, bi = bi - ai
+        eiisinf (ai)            1 if infinite
+        eiisnan (ai)            1 if a NaN
+       eiisneg (ai)            1 if sign bit of ai != 0, else 0
+        einan (ai)              set ai = NaN
+        eiinfin (ai)            set ai = infinity
+
+  The result is always normalized and rounded to NI-4 word precision
+  after each arithmetic operation.
+
+  Exception flags are NOT fully supported.
+ 
+  Signaling NaN's are NOT supported; they are treated the same
+  as quiet NaN's.
+ 
+  Define INFINITY for support of infinity; otherwise a
+  saturation arithmetic is implemented.
+ 
+  Define NANS for support of Not-a-Number items; otherwise the
+  arithmetic will never produce a NaN output, and might be confused
+  by a NaN input.
+  If NaN's are supported, the output of `ecmp (a,b)' is -2 if
+  either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
+  may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
+  if in doubt.
+ 
+  Denormals are always supported here where appropriate (e.g., not
+  for conversion to DEC numbers).  */
+
+/* Definitions for error codes that are passed to the common error handling
+   routine mtherr.
+
+   For Digital Equipment PDP-11 and VAX computers, certain
+  IBM systems, and others that use numbers with a 56-bit
+  significand, the symbol DEC should be defined.  In this
+  mode, most floating point constants are given as arrays
+  of octal integers to eliminate decimal to binary conversion
+  errors that might be introduced by the compiler.
+ 
+  For computers, such as IBM PC, that follow the IEEE
+  Standard for Binary Floating Point Arithmetic (ANSI/IEEE
+  Std 754-1985), the symbol IEEE should be defined.
+  These numbers have 53-bit significands.  In this mode, constants
+  are provided as arrays of hexadecimal 16 bit integers.
+  The endian-ness of generated values is controlled by
+  REAL_WORDS_BIG_ENDIAN.
+ 
+  To accommodate other types of computer arithmetic, all
+  constants are also provided in a normal decimal radix
+  which one can hope are correctly converted to a suitable
+  format by the available C language compiler.  To invoke
+  this mode, the symbol UNK is defined.
+ 
+  An important difference among these modes is a predefined
+  set of machine arithmetic constants for each.  The numbers
+  MACHEP (the machine roundoff error), MAXNUM (largest number
+  represented), and several other parameters are preset by
+  the configuration symbol.  Check the file const.c to
+  ensure that these values are correct for your computer.
+ 
+  For ANSI C compatibility, define ANSIC equal to 1.  Currently
+  this affects only the atan2 function and others that use it.  */

 
 

-
-/*                                                     mconf.h
- *
- *     Common include file for math routines
- *
- *
- *
- * SYNOPSIS:
- *
- * #include "mconf.h"
- *
- *
- *
- * DESCRIPTION:
- *
- * This file contains definitions for error codes that are
- * passed to the common error handling routine mtherr
- * (which see).
- *
- * The file also includes a conditional assembly definition
- * for the type of computer arithmetic (Intel IEEE, DEC, Motorola
- * IEEE, or UNKnown).
- *
- * For Digital Equipment PDP-11 and VAX computers, certain
- * IBM systems, and others that use numbers with a 56-bit
- * significand, the symbol DEC should be defined.  In this
- * mode, most floating point constants are given as arrays
- * of octal integers to eliminate decimal to binary conversion
- * errors that might be introduced by the compiler.
- *
- * For computers, such as IBM PC, that follow the IEEE
- * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
- * Std 754-1985), the symbol IBMPC or MIEEE should be defined.
- * These numbers have 53-bit significands.  In this mode, constants
- * are provided as arrays of hexadecimal 16 bit integers.
- *
- * To accommodate other types of computer arithmetic, all
- * constants are also provided in a normal decimal radix
- * which one can hope are correctly converted to a suitable
- * format by the available C language compiler.  To invoke
- * this mode, the symbol UNK is defined.
- *
- * An important difference among these modes is a predefined
- * set of machine arithmetic constants for each.  The numbers
- * MACHEP (the machine roundoff error), MAXNUM (largest number
- * represented), and several other parameters are preset by
- * the configuration symbol.  Check the file const.c to
- * ensure that these values are correct for your computer.
- *
- * For ANSI C compatibility, define ANSIC equal to 1.  Currently
- * this affects only the atan2 function and others that use it.
- */
-\f

 /* Constant definitions for math error conditions.  */
 
 #define DOMAIN         1       /* argument domain error */
 /* Constant definitions for math error conditions.  */
 
 #define DOMAIN         1       /* argument domain error */


@@ -990,107 +1402,88 @@ ereal_isneg (x)
 #define UNDERFLOW      4       /* underflow range error */
 #define TLOSS          5       /* total loss of precision */
 #define PLOSS          6       /* partial loss of precision */
 #define UNDERFLOW      4       /* underflow range error */
 #define TLOSS          5       /* total loss of precision */
 #define PLOSS          6       /* partial loss of precision */

+#define INVALID                7       /* NaN-producing operation */

 
 /*  e type constants used by high precision check routines */
 
 
 /*  e type constants used by high precision check routines */
 

-/*include "ehead.h"*/
+#if LONG_DOUBLE_TYPE_SIZE == 128

 /* 0.0 */
 unsigned EMUSHORT ezero[NE] =
 /* 0.0 */
 unsigned EMUSHORT ezero[NE] =

-{
-  0, 0000000, 0000000, 0000000, 0000000, 0000000,};
+ {0x0000, 0x0000, 0x0000, 0x0000,
+  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};

 extern unsigned EMUSHORT ezero[];
 
 /* 5.0E-1 */
 unsigned EMUSHORT ehalf[NE] =
 extern unsigned EMUSHORT ezero[];
 
 /* 5.0E-1 */
 unsigned EMUSHORT ehalf[NE] =

-{
-  0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
+ {0x0000, 0x0000, 0x0000, 0x0000,
+  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};

 extern unsigned EMUSHORT ehalf[];
 
 /* 1.0E0 */
 unsigned EMUSHORT eone[NE] =
 extern unsigned EMUSHORT ehalf[];
 
 /* 1.0E0 */
 unsigned EMUSHORT eone[NE] =

-{
-  0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
+ {0x0000, 0x0000, 0x0000, 0x0000,
+  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};

 extern unsigned EMUSHORT eone[];
 
 /* 2.0E0 */
 unsigned EMUSHORT etwo[NE] =
 extern unsigned EMUSHORT eone[];
 
 /* 2.0E0 */
 unsigned EMUSHORT etwo[NE] =

-{
-  0, 0000000, 0000000, 0000000, 0100000, 0040000,};
+ {0x0000, 0x0000, 0x0000, 0x0000,
+  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};

 extern unsigned EMUSHORT etwo[];
 
 /* 3.2E1 */
 unsigned EMUSHORT e32[NE] =
 extern unsigned EMUSHORT etwo[];
 
 /* 3.2E1 */
 unsigned EMUSHORT e32[NE] =

-{
-  0, 0000000, 0000000, 0000000, 0100000, 0040004,};
+ {0x0000, 0x0000, 0x0000, 0x0000,
+  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};

 extern unsigned EMUSHORT e32[];
 
 /* 6.93147180559945309417232121458176568075500134360255E-1 */
 unsigned EMUSHORT elog2[NE] =
 extern unsigned EMUSHORT e32[];
 
 /* 6.93147180559945309417232121458176568075500134360255E-1 */
 unsigned EMUSHORT elog2[NE] =

-{
-  0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
+ {0x40f3, 0xf6af, 0x03f2, 0xb398,
+  0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};

 extern unsigned EMUSHORT elog2[];
 
 /* 1.41421356237309504880168872420969807856967187537695E0 */
 unsigned EMUSHORT esqrt2[NE] =
 extern unsigned EMUSHORT elog2[];
 
 /* 1.41421356237309504880168872420969807856967187537695E0 */
 unsigned EMUSHORT esqrt2[NE] =

-{
-  0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
+ {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
+  0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};

 extern unsigned EMUSHORT esqrt2[];
 
 extern unsigned EMUSHORT esqrt2[];
 

-/* 2/sqrt (PI) =
- * 1.12837916709551257389615890312154517168810125865800E0 */
-unsigned EMUSHORT eoneopi[NE] =
-{
-  0x71d5, 0x688d, 0012333, 0135202, 0110156, 0x3fff,};
-extern unsigned EMUSHORT eoneopi[];
-

 /* 3.14159265358979323846264338327950288419716939937511E0 */
 unsigned EMUSHORT epi[NE] =
 /* 3.14159265358979323846264338327950288419716939937511E0 */
 unsigned EMUSHORT epi[NE] =

-{
+ {0x2902, 0x1cd1, 0x80dc, 0x628b,

   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
 extern unsigned EMUSHORT epi[];
 
   0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
 extern unsigned EMUSHORT epi[];
 

-/* 5.7721566490153286060651209008240243104215933593992E-1 */
-unsigned EMUSHORT eeul[NE] =
-{
-  0xd1be, 0xc7a4, 0076660, 0063743, 0111704, 0x3ffe,};
-extern unsigned EMUSHORT eeul[];
-
-/*
-include "ehead.h"
-include "mconf.h"
-*/
-
-
+#else
+/* LONG_DOUBLE_TYPE_SIZE is other than 128 */
+unsigned EMUSHORT ezero[NE] =
+ {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
+unsigned EMUSHORT ehalf[NE] =
+ {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
+unsigned EMUSHORT eone[NE] =
+ {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
+unsigned EMUSHORT etwo[NE] =
+ {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
+unsigned EMUSHORT e32[NE] =
+ {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
+unsigned EMUSHORT elog2[NE] =
+ {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
+unsigned EMUSHORT esqrt2[NE] =
+ {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
+unsigned EMUSHORT epi[NE] =
+ {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
+#endif

 
 /* Control register for rounding precision.
 
 /* Control register for rounding precision.

- * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
- */
+   This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits.  */
+

 int rndprc = NBITS;
 extern int rndprc;
 
 int rndprc = NBITS;
 extern int rndprc;
 

-void eaddm (), esubm (), emdnorm (), asctoeg ();
-static void toe24 (), toe53 (), toe64 ();
-void eremain (), einit (), eiremain ();
-int ecmpm (), edivm (), emulm ();
-void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
-void etodec (), todec (), dectoe ();
-
-
-
-
-void 
-einit ()
-{
-}
-
-/*
-; Clear out entire external format number.
-;
-; unsigned EMUSHORT x[];
-; eclear (x);
-*/
+/*  Clear out entire e-type number X.  */

 
 

-void 
+static void 

 eclear (x)
      register unsigned EMUSHORT *x;
 {
 eclear (x)
      register unsigned EMUSHORT *x;
 {


@@ -1100,14 +1493,9 @@ eclear (x)
     *x++ = 0;
 }
 
     *x++ = 0;
 }
 

+/* Move e-type number from A to B.  */

 
 

-
-/* Move external format number from a to b.
- *
- * emov (a, b);
- */
-
-void 
+static void 

 emov (a, b)
      register unsigned EMUSHORT *a, *b;
 {
 emov (a, b)
      register unsigned EMUSHORT *a, *b;
 {


@@ -1118,32 +1506,19 @@ emov (a, b)
 }
 
 
 }
 
 

-/*
-;      Absolute value of external format number
-;
-;      EMUSHORT x[NE];
-;      eabs (x);
-*/
+/* Absolute value of e-type X.  */

 
 

-void 
+static void 

 eabs (x)
 eabs (x)

-     unsigned EMUSHORT x[];    /* x is the memory address of a short */
+     unsigned EMUSHORT x[];

 {
 {

-
-  x[NE - 1] &= 0x7fff;         /* sign is top bit of last word of external format */
+  /* sign is top bit of last word of external format */
+  x[NE - 1] &= 0x7fff;         

 }
 
 }
 

+/* Negate the e-type number X.  */

 
 

-
-
-/*
-;      Negate external format number
-;
-;      unsigned EMUSHORT x[NE];
-;      eneg (x);
-*/
-
-void 
+static void 

 eneg (x)
      unsigned EMUSHORT x[];
 {
 eneg (x)
      unsigned EMUSHORT x[];
 {


@@ -1151,12 +1526,9 @@ eneg (x)
   x[NE - 1] ^= 0x8000;         /* Toggle the sign bit */
 }
 
   x[NE - 1] ^= 0x8000;         /* Toggle the sign bit */
 }
 

+/* Return 1 if sign bit of e-type number X is nonzero, else zero.  */

 
 

-
-/* Return 1 if external format number is negative,
- * else return zero.
- */
-int 
+static int 

 eisneg (x)
      unsigned EMUSHORT x[];
 {
 eisneg (x)
      unsigned EMUSHORT x[];
 {


@@ -1167,32 +1539,51 @@ eisneg (x)
     return (0);
 }
 
     return (0);
 }
 

+/* Return 1 if e-type number X is infinity, else return zero.  */

 
 

-/* Return 1 if external format number has maximum possible exponent,
- * else return zero.
- */
-int 
+static int 

 eisinf (x)
      unsigned EMUSHORT x[];
 {
 
 eisinf (x)
      unsigned EMUSHORT x[];
 {
 

+#ifdef NANS
+  if (eisnan (x))
+    return (0);
+#endif

   if ((x[NE - 1] & 0x7fff) == 0x7fff)
     return (1);
   else
     return (0);
 }
 
   if ((x[NE - 1] & 0x7fff) == 0x7fff)
     return (1);
   else
     return (0);
 }
 

+/* Check if e-type number is not a number.  The bit pattern is one that we
+   defined, so we know for sure how to detect it.  */

 
 

-/*
-; Fill entire number, including exponent and significand, with
-; largest possible number.  These programs implement a saturation
-; value that is an ordinary, legal number.  A special value
-; "infinity" may also be implemented; this would require tests
-; for that value and implementation of special rules for arithmetic
-; operations involving inifinity.
-*/
+static int 
+eisnan (x)
+     unsigned EMUSHORT x[];
+{
+#ifdef NANS
+  int i;

 
 

-void 
+  /* NaN has maximum exponent */
+  if ((x[NE - 1] & 0x7fff) != 0x7fff)
+    return (0);
+  /* ... and non-zero significand field.  */
+  for (i = 0; i < NE - 1; i++)
+    {
+      if (*x++ != 0)
+        return (1);
+    }
+#endif
+
+  return (0);
+}
+
+/*  Fill e-type number X with infinity pattern (IEEE)
+    or largest possible number (non-IEEE).  */
+
+static void 

 einfin (x)
      register unsigned EMUSHORT *x;
 {
 einfin (x)
      register unsigned EMUSHORT *x;
 {


@@ -1208,6 +1599,11 @@ einfin (x)
   *x |= 32766;
   if (rndprc < NBITS)
     {
   *x |= 32766;
   if (rndprc < NBITS)
     {

+      if (rndprc == 113)
+       {
+         *(x - 9) = 0;
+         *(x - 8) = 0;
+       }

       if (rndprc == 64)
        {
          *(x - 5) = 0;
       if (rndprc == 64)
        {
          *(x - 5) = 0;


@@ -1226,12 +1622,26 @@ einfin (x)
 #endif
 }
 
 #endif
 }
 

+/* Output an e-type NaN.
+   This generates Intel's quiet NaN pattern for extended real.
+   The exponent is 7fff, the leading mantissa word is c000.  */
+
+static void 
+enan (x, sign)
+     register unsigned EMUSHORT *x;
+     int sign;
+{
+  register int i;
+
+  for (i = 0; i < NE - 2; i++)
+    *x++ = 0;
+  *x++ = 0xc000;
+  *x = (sign << 15) | 0x7fff;
+}

 
 

+/* Move in an e-type number A, converting it to exploded e-type B.  */

 
 

-/* Move in external format number,
- * converting it to internal format.
- */
-void 
+static void 

 emovi (a, b)
      unsigned EMUSHORT *a, *b;
 {
 emovi (a, b)
      unsigned EMUSHORT *a, *b;
 {


@@ -1251,11 +1661,22 @@ emovi (a, b)
 #ifdef INFINITY
   if ((*(q - 1) & 0x7fff) == 0x7fff)
     {
 #ifdef INFINITY
   if ((*(q - 1) & 0x7fff) == 0x7fff)
     {

+#ifdef NANS
+      if (eisnan (a))
+       {
+         *q++ = 0;
+         for (i = 3; i < NI; i++)
+           *q++ = *p--;
+         return;
+       }
+#endif
+

       for (i = 2; i < NI; i++)
        *q++ = 0;
       return;
     }
 #endif
       for (i = 2; i < NI; i++)
        *q++ = 0;
       return;
     }
 #endif

+

   /* clear high guard word */
   *q++ = 0;
   /* move in the significand */
   /* clear high guard word */
   *q++ = 0;
   /* move in the significand */


@@ -1265,16 +1686,15 @@ emovi (a, b)
   *q = 0;
 }
 
   *q = 0;
 }
 

+/* Move out exploded e-type number A, converting it to e type B.  */

 
 

-/* Move internal format number out,
- * converting it to external format.
- */
-void 
+static void 

 emovo (a, b)
      unsigned EMUSHORT *a, *b;
 {
   register unsigned EMUSHORT *p, *q;
   unsigned EMUSHORT i;
 emovo (a, b)
      unsigned EMUSHORT *a, *b;
 {
   register unsigned EMUSHORT *p, *q;
   unsigned EMUSHORT i;

+  int j;

 
   p = a;
   q = b + (NE - 1);            /* point to output exponent */
 
   p = a;
   q = b + (NE - 1);            /* point to output exponent */


@@ -1287,24 +1707,27 @@ emovo (a, b)
 #ifdef INFINITY
   if (*(p - 1) == 0x7fff)
     {
 #ifdef INFINITY
   if (*(p - 1) == 0x7fff)
     {

+#ifdef NANS
+      if (eiisnan (a))
+       {
+         enan (b, eiisneg (a));
+         return;
+       }
+#endif

       einfin (b);
       einfin (b);

-      return;
+       return;

     }
 #endif
   /* skip over guard word */
   ++p;
   /* move the significand */
     }
 #endif
   /* skip over guard word */
   ++p;
   /* move the significand */

-  for (i = 0; i < NE - 1; i++)
+  for (j = 0; j < NE - 1; j++)

     *q-- = *p++;
 }
 
     *q-- = *p++;
 }
 

+/* Clear out exploded e-type number XI.  */

 
 

-
-
-/* Clear out internal format number.
- */
-
-void 
+static void 

 ecleaz (xi)
      register unsigned EMUSHORT *xi;
 {
 ecleaz (xi)
      register unsigned EMUSHORT *xi;
 {


@@ -1314,10 +1737,9 @@ ecleaz (xi)
     *xi++ = 0;
 }
 
     *xi++ = 0;
 }
 

+/* Clear out exploded e-type XI, but don't touch the sign.  */

 
 

-/* same, but don't touch the sign. */
-
-void 
+static void 

 ecleazs (xi)
      register unsigned EMUSHORT *xi;
 {
 ecleazs (xi)
      register unsigned EMUSHORT *xi;
 {


@@ -1328,11 +1750,9 @@ ecleazs (xi)
     *xi++ = 0;
 }
 
     *xi++ = 0;
 }
 

+/* Move exploded e-type number from A to B.  */

 
 

-
-/* Move internal format number from a to b.
- */
-void 
+static void 

 emovz (a, b)
      register unsigned EMUSHORT *a, *b;
 {
 emovz (a, b)
      register unsigned EMUSHORT *a, *b;
 {


@@ -1344,20 +1764,86 @@ emovz (a, b)
   *b = 0;
 }
 
   *b = 0;
 }
 

+/* Generate exploded e-type NaN.
+   The explicit pattern for this is maximum exponent and
+   top two significant bits set.  */

 
 

-/*
-;      Compare significands of numbers in internal format.
-;      Guard words are included in the comparison.
-;
-;      unsigned EMUSHORT a[NI], b[NI];
-;      cmpm (a, b);
-;
-;      for the significands:
-;      returns +1 if a > b
-;               0 if a == b
-;              -1 if a < b
-*/
-int
+static void
+einan (x)
+     unsigned EMUSHORT x[];
+{
+
+  ecleaz (x);
+  x[E] = 0x7fff;
+  x[M + 1] = 0xc000;
+}
+
+/* Return nonzero if exploded e-type X is a NaN.  */
+
+static int 
+eiisnan (x)
+     unsigned EMUSHORT x[];
+{
+  int i;
+
+  if ((x[E] & 0x7fff) == 0x7fff)
+    {
+      for (i = M + 1; i < NI; i++)
+       {
+         if (x[i] != 0)
+           return (1);
+       }
+    }
+  return (0);
+}
+
+/* Return nonzero if sign of exploded e-type X is nonzero.  */
+
+static int 
+eiisneg (x)
+     unsigned EMUSHORT x[];
+{
+
+  return x[0] != 0;
+}
+
+/* Fill exploded e-type X with infinity pattern.
+   This has maximum exponent and significand all zeros.  */
+
+static void
+eiinfin (x)
+     unsigned EMUSHORT x[];
+{
+
+  ecleaz (x);
+  x[E] = 0x7fff;
+}
+
+/* Return nonzero if exploded e-type X is infinite.  */
+
+static int 
+eiisinf (x)
+     unsigned EMUSHORT x[];
+{
+
+#ifdef NANS
+  if (eiisnan (x))
+    return (0);
+#endif
+  if ((x[E] & 0x7fff) == 0x7fff)
+    return (1);
+  return (0);
+}
+
+
+/* Compare significands of numbers in internal exploded e-type format.
+   Guard words are included in the comparison.
+
+   Returns     +1 if a > b
+                0 if a == b
+               -1 if a < b   */
+
+static int

 ecmpm (a, b)
      register unsigned EMUSHORT *a, *b;
 {
 ecmpm (a, b)
      register unsigned EMUSHORT *a, *b;
 {


@@ -1379,12 +1865,9 @@ ecmpm (a, b)
     return (-1);
 }
 
     return (-1);
 }
 

+/* Shift significand of exploded e-type X down by 1 bit.  */

 
 

-/*
-;      Shift significand down by 1 bit
-*/
-
-void 
+static void 

 eshdn1 (x)
      register unsigned EMUSHORT *x;
 {
 eshdn1 (x)
      register unsigned EMUSHORT *x;
 {


@@ -1406,13 +1889,9 @@ eshdn1 (x)
     }
 }
 
     }
 }
 

+/* Shift significand of exploded e-type X up by 1 bit.  */

 
 

-
-/*
-;      Shift significand up by 1 bit
-*/
-
-void 
+static void 

 eshup1 (x)
      register unsigned EMUSHORT *x;
 {
 eshup1 (x)
      register unsigned EMUSHORT *x;
 {


@@ -1435,12 +1914,9 @@ eshup1 (x)
 }
 
 
 }
 
 

+/* Shift significand of exploded e-type X down by 8 bits.  */

 
 

-/*
-;      Shift significand down by 8 bits
-*/
-
-void 
+static void 

 eshdn8 (x)
      register unsigned EMUSHORT *x;
 {
 eshdn8 (x)
      register unsigned EMUSHORT *x;
 {


@@ -1459,11 +1935,9 @@ eshdn8 (x)
     }
 }
 
     }
 }
 

-/*
-;      Shift significand up by 8 bits
-*/
+/* Shift significand of exploded e-type X up by 8 bits.  */

 
 

-void 
+static void 

 eshup8 (x)
      register unsigned EMUSHORT *x;
 {
 eshup8 (x)
      register unsigned EMUSHORT *x;
 {


@@ -1483,11 +1957,9 @@ eshup8 (x)
     }
 }
 
     }
 }
 

-/*
-;      Shift significand up by 16 bits
-*/
+/* Shift significand of exploded e-type X up by 16 bits.  */

 
 

-void 
+static void 

 eshup6 (x)
      register unsigned EMUSHORT *x;
 {
 eshup6 (x)
      register unsigned EMUSHORT *x;
 {


@@ -1503,11 +1975,9 @@ eshup6 (x)
   *p = 0;
 }
 
   *p = 0;
 }
 

-/*
-;      Shift significand down by 16 bits
-*/
+/* Shift significand of exploded e-type X down by 16 bits.  */

 
 

-void 
+static void 

 eshdn6 (x)
      register unsigned EMUSHORT *x;
 {
 eshdn6 (x)
      register unsigned EMUSHORT *x;
 {


@@ -1522,13 +1992,10 @@ eshdn6 (x)
 
   *(--p) = 0;
 }
 
   *(--p) = 0;
 }

-\f
-/*
-;      Add significands
-;      x + y replaces y
-*/

 
 

-void 
+/* Add significands of exploded e-type X and Y.  X + Y replaces Y.  */
+
+static void 

 eaddm (x, y)
      unsigned EMUSHORT *x, *y;
 {
 eaddm (x, y)
      unsigned EMUSHORT *x, *y;
 {


@@ -1552,12 +2019,9 @@ eaddm (x, y)
     }
 }
 
     }
 }
 

-/*
-;      Subtract significands
-;      y - x replaces y
-*/
+/* Subtract significands of exploded e-type X and Y.  Y - X replaces Y.  */

 
 

-void 
+static void 

 esubm (x, y)
      unsigned EMUSHORT *x, *y;
 {
 esubm (x, y)
      unsigned EMUSHORT *x, *y;
 {


@@ -1582,10 +2046,15 @@ esubm (x, y)
 }
 
 
 }
 
 

-/* Divide significands */
-

 static unsigned EMUSHORT equot[NI];
 
 static unsigned EMUSHORT equot[NI];
 

+
+#if 0
+/* Radix 2 shift-and-add versions of multiply and divide  */
+
+
+/* Divide significands */
+

 int 
 edivm (den, num)
      unsigned EMUSHORT den[], num[];
 int 
 edivm (den, num)
      unsigned EMUSHORT den[], num[];


@@ -1603,9 +2072,9 @@ edivm (den, num)
       *p++ = 0;
     }
 
       *p++ = 0;
     }
 

-  /* Use faster compare and subtraction if denominator
-   * has only 15 bits of significance.
-   */
+  /* Use faster compare and subtraction if denominator has only 15 bits of
+     significance.  */
+

   p = &den[M + 2];
   if (*p++ == 0)
     {
   p = &den[M + 2];
   if (*p++ == 0)
     {


@@ -1640,9 +2109,9 @@ edivm (den, num)
       goto divdon;
     }
 
       goto divdon;
     }
 

-  /* The number of quotient bits to calculate is
-   * NBITS + 1 scaling guard bit + 1 roundoff bit.
-   */
+  /* The number of quotient bits to calculate is NBITS + 1 scaling guard
+     bit + 1 roundoff bit.  */
+

  fulldiv:
 
   p = &equot[NI - 2];
  fulldiv:
 
   p = &equot[NI - 2];


@@ -1683,6 +2152,7 @@ edivm (den, num)
 
 
 /* Multiply significands */
 
 
 /* Multiply significands */

+

 int 
 emulm (a, b)
      unsigned EMUSHORT a[], b[];
 int 
 emulm (a, b)
      unsigned EMUSHORT a[], b[];


@@ -1697,7 +2167,7 @@ emulm (a, b)
 
   p = &a[NI - 2];
   k = NBITS;
 
   p = &a[NI - 2];
   k = NBITS;

-  while (*p == 0)              /* significand is not supposed to be all zero */
+  while (*p == 0)              /* significand is not supposed to be zero */

     {
       eshdn6 (a);
       k -= 16;
     {
       eshdn6 (a);
       k -= 16;


@@ -1728,61 +2198,225 @@ emulm (a, b)
   return (j);
 }
 
   return (j);
 }
 

+#else

 
 

+/* Radix 65536 versions of multiply and divide.  */

 
 

-/*
- * Normalize and round off.
- *
- * The internal format number to be rounded is "s".
- * Input "lost" indicates whether or not the number is exact.
- * This is the so-called sticky bit.
- *
- * Input "subflg" indicates whether the number was obtained
- * by a subtraction operation.  In that case if lost is nonzero
- * then the number is slightly smaller than indicated.
- *
- * Input "exp" is the biased exponent, which may be negative.
- * the exponent field of "s" is ignored but is replaced by
- * "exp" as adjusted by normalization and rounding.
- *
- * Input "rcntrl" is the rounding control.
- */
-
-static int rlast = -1;
-static int rw = 0;
-static unsigned EMUSHORT rmsk = 0;
-static unsigned EMUSHORT rmbit = 0;
-static unsigned EMUSHORT rebit = 0;
-static int re = 0;
-static unsigned EMUSHORT rbit[NI];
+/* Multiply significand of e-type number B
+   by 16-bit quantity A, return e-type result to C.  */

 
 

-void 
-emdnorm (s, lost, subflg, exp, rcntrl)
-     unsigned EMUSHORT s[];
-     int lost;
-     int subflg;
-     EMULONG exp;
-     int rcntrl;
+static void
+m16m (a, b, c)
+     unsigned int a;
+     unsigned EMUSHORT b[], c[];

 {
 {

-  int i, j;
-  unsigned EMUSHORT r;
+  register unsigned EMUSHORT *pp;
+  register unsigned EMULONG carry;
+  unsigned EMUSHORT *ps;
+  unsigned EMUSHORT p[NI];
+  unsigned EMULONG aa, m;
+  int i;

 
 

-  /* Normalize */
-  j = enormlz (s);
+  aa = a;
+  pp = &p[NI-2];
+  *pp++ = 0;
+  *pp = 0;
+  ps = &b[NI-1];

 
 

-  /* a blank significand could mean either zero or infinity. */
-#ifndef INFINITY
-  if (j > NBITS)
+  for (i=M+1; i<NI; i++)

     {
     {

-      ecleazs (s);
-      return;
+      if (*ps == 0)
+       {
+         --ps;
+         --pp;
+         *(pp-1) = 0;
+       }
+      else
+       {
+         m = (unsigned EMULONG) aa * *ps--;
+         carry = (m & 0xffff) + *pp;
+         *pp-- = (unsigned EMUSHORT)carry;
+         carry = (carry >> 16) + (m >> 16) + *pp;
+         *pp = (unsigned EMUSHORT)carry;
+         *(pp-1) = carry >> 16;
+       }

     }
     }

-#endif
-  exp -= j;
-#ifndef INFINITY
-  if (exp >= 32767L)
-    goto overf;
-#else
+  for (i=M; i<NI; i++)
+    c[i] = p[i];
+}
+
+/* Divide significands of exploded e-types NUM / DEN.  Neither the
+   numerator NUM nor the denominator DEN is permitted to have its high guard
+   word nonzero.  */
+
+static int
+edivm (den, num)
+     unsigned EMUSHORT den[], num[];
+{
+  int i;
+  register unsigned EMUSHORT *p;
+  unsigned EMULONG tnum;
+  unsigned EMUSHORT j, tdenm, tquot;
+  unsigned EMUSHORT tprod[NI+1];
+
+  p = &equot[0];
+  *p++ = num[0];
+  *p++ = num[1];
+
+  for (i=M; i<NI; i++)
+    {
+      *p++ = 0;
+    }
+  eshdn1 (num);
+  tdenm = den[M+1];
+  for (i=M; i<NI; i++)
+    {
+      /* Find trial quotient digit (the radix is 65536).  */
+      tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
+
+      /* Do not execute the divide instruction if it will overflow.  */
+      if ((tdenm * 0xffffL) < tnum)
+       tquot = 0xffff;
+      else
+       tquot = tnum / tdenm;
+      /* Multiply denominator by trial quotient digit.  */
+      m16m ((unsigned int)tquot, den, tprod);
+      /* The quotient digit may have been overestimated.  */
+      if (ecmpm (tprod, num) > 0)
+       {
+         tquot -= 1;
+         esubm (den, tprod);
+         if (ecmpm (tprod, num) > 0)
+           {
+             tquot -= 1;
+             esubm (den, tprod);
+           }
+       }
+      esubm (tprod, num);
+      equot[i] = tquot;
+      eshup6(num);
+    }
+  /* test for nonzero remainder after roundoff bit */
+  p = &num[M];
+  j = 0;
+  for (i=M; i<NI; i++)
+    {
+      j |= *p++;
+    }
+  if (j)
+    j = 1;
+
+  for (i=0; i<NI; i++)
+    num[i] = equot[i];
+
+  return ((int)j);
+}
+
+/* Multiply significands of exploded e-type A and B, result in B.  */
+
+static int
+emulm (a, b)
+     unsigned EMUSHORT a[], b[];
+{
+  unsigned EMUSHORT *p, *q;
+  unsigned EMUSHORT pprod[NI];
+  unsigned EMUSHORT j;
+  int i;
+
+  equot[0] = b[0];
+  equot[1] = b[1];
+  for (i=M; i<NI; i++)
+    equot[i] = 0;
+
+  j = 0;
+  p = &a[NI-1];
+  q = &equot[NI-1];
+  for (i=M+1; i<NI; i++)
+    {
+      if (*p == 0)
+       {
+         --p;
+       }
+      else
+       {
+         m16m ((unsigned int) *p--, b, pprod);
+         eaddm(pprod, equot);
+       }
+      j |= *q;
+      eshdn6(equot);
+    }
+
+  for (i=0; i<NI; i++)
+    b[i] = equot[i];
+
+  /* return flag for lost nonzero bits */
+  return ((int)j);
+}
+#endif
+
+
+/* Normalize and round off.
+
+  The internal format number to be rounded is S.
+  Input LOST is 0 if the value is exact.  This is the so-called sticky bit.
+ 
+  Input SUBFLG indicates whether the number was obtained
+  by a subtraction operation.  In that case if LOST is nonzero
+  then the number is slightly smaller than indicated.
+ 
+  Input EXP is the biased exponent, which may be negative.
+  the exponent field of S is ignored but is replaced by
+  EXP as adjusted by normalization and rounding.
+ 
+  Input RCNTRL is the rounding control.  If it is nonzero, the
+  returned value will be rounded to RNDPRC bits.
+
+  For future reference:  In order for emdnorm to round off denormal
+   significands at the right point, the input exponent must be
+   adjusted to be the actual value it would have after conversion to
+   the final floating point type.  This adjustment has been
+   implemented for all type conversions (etoe53, etc.) and decimal
+   conversions, but not for the arithmetic functions (eadd, etc.). 
+   Data types having standard 15-bit exponents are not affected by
+   this, but SFmode and DFmode are affected. For example, ediv with
+   rndprc = 24 will not round correctly to 24-bit precision if the
+   result is denormal.   */
+
+static int rlast = -1;
+static int rw = 0;
+static unsigned EMUSHORT rmsk = 0;
+static unsigned EMUSHORT rmbit = 0;
+static unsigned EMUSHORT rebit = 0;
+static int re = 0;
+static unsigned EMUSHORT rbit[NI];
+
+static void 
+emdnorm (s, lost, subflg, exp, rcntrl)
+     unsigned EMUSHORT s[];
+     int lost;
+     int subflg;
+     EMULONG exp;
+     int rcntrl;
+{
+  int i, j;
+  unsigned EMUSHORT r;
+
+  /* Normalize */
+  j = enormlz (s);
+
+  /* a blank significand could mean either zero or infinity.  */
+#ifndef INFINITY
+  if (j > NBITS)
+    {
+      ecleazs (s);
+      return;
+    }
+#endif
+  exp -= j;
+#ifndef INFINITY
+  if (exp >= 32767L)
+    goto overf;
+#else

   if ((j > NBITS) && (exp < 32767))
     {
       ecleazs (s);
   if ((j > NBITS) && (exp < 32767))
     {
       ecleazs (s);


@@ -1804,10 +2438,10 @@ emdnorm (s, lost, subflg, exp, rcntrl)
          return;
        }
     }
          return;
        }
     }

-  /* Round off, unless told not to by rcntrl. */
+  /* Round off, unless told not to by rcntrl.  */

   if (rcntrl == 0)
     goto mdfin;
   if (rcntrl == 0)
     goto mdfin;

-  /* Set up rounding parameters if the control register changed. */
+  /* Set up rounding parameters if the control register changed.  */

   if (rndprc != rlast)
     {
       ecleaz (rbit);
   if (rndprc != rlast)
     {
       ecleaz (rbit);


@@ -1818,68 +2452,64 @@ emdnorm (s, lost, subflg, exp, rcntrl)
          rw = NI - 1;          /* low guard word */
          rmsk = 0xffff;
          rmbit = 0x8000;
          rw = NI - 1;          /* low guard word */
          rmsk = 0xffff;
          rmbit = 0x8000;

-         rbit[rw - 1] = 1;
-         re = NI - 2;
+         re = rw - 1;

          rebit = 1;
          break;
          rebit = 1;
          break;

+       case 113:
+         rw = 10;
+         rmsk = 0x7fff;
+         rmbit = 0x4000;
+         rebit = 0x8000;
+         re = rw;
+         break;

        case 64:
          rw = 7;
          rmsk = 0xffff;
          rmbit = 0x8000;
        case 64:
          rw = 7;
          rmsk = 0xffff;
          rmbit = 0x8000;

-         rbit[rw - 1] = 1;

          re = rw - 1;
          rebit = 1;
          break;
          re = rw - 1;
          rebit = 1;
          break;

-         /* For DEC arithmetic */
+         /* For DEC or IBM arithmetic */

        case 56:
          rw = 6;
          rmsk = 0xff;
          rmbit = 0x80;
        case 56:
          rw = 6;
          rmsk = 0xff;
          rmbit = 0x80;

-         rbit[rw] = 0x100;
-         re = rw;

          rebit = 0x100;
          rebit = 0x100;

+         re = rw;

          break;
        case 53:
          rw = 6;
          rmsk = 0x7ff;
          rmbit = 0x0400;
          break;
        case 53:
          rw = 6;
          rmsk = 0x7ff;
          rmbit = 0x0400;

-         rbit[rw] = 0x800;
-         re = rw;

          rebit = 0x800;
          rebit = 0x800;

+         re = rw;

          break;
        case 24:
          rw = 4;
          rmsk = 0xff;
          rmbit = 0x80;
          break;
        case 24:
          rw = 4;
          rmsk = 0xff;
          rmbit = 0x80;

-         rbit[rw] = 0x100;
-         re = rw;

          rebit = 0x100;
          rebit = 0x100;

+         re = rw;

          break;
        }
          break;
        }

+      rbit[re] = rebit;

       rlast = rndprc;
     }
 
       rlast = rndprc;
     }
 

-  if (rndprc >= 64)
+  /* Shift down 1 temporarily if the data structure has an implied
+     most significant bit and the number is denormal.
+     Intel long double denormals also lose one bit of precision.  */
+  if ((exp <= 0) && (rndprc != NBITS)
+      && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))

     {
     {

-      r = s[rw] & rmsk;
-      if (rndprc == 64)
-       {
-         i = rw + 1;
-         while (i < NI)
-           {
-             if (s[i])
-               r |= 1;
-             s[i] = 0;
-             ++i;
-           }
-       }
+      lost |= s[NI - 1] & 1;
+      eshdn1 (s);

     }
     }

-  else
+  /* Clear out all bits below the rounding bit,
+     remembering in r if any were nonzero.  */
+  r = s[rw] & rmsk;
+  if (rndprc < NBITS)

     {
     {

-      if (exp <= 0)
-       eshdn1 (s);
-      r = s[rw] & rmsk;
-      /* These tests assume NI = 8 */

       i = rw + 1;
       while (i < NI)
        {
       i = rw + 1;
       while (i < NI)
        {


@@ -1888,17 +2518,8 @@ emdnorm (s, lost, subflg, exp, rcntrl)
          s[i] = 0;
          ++i;
        }
          s[i] = 0;
          ++i;
        }

-      /*
-        if (rndprc == 24)
-        {
-        if (s[5] || s[6])
-        r |= 1;
-        s[5] = 0;
-        s[6] = 0;
-        }
-        */
-      s[rw] &= ~rmsk;

     }
     }

+  s[rw] &= ~rmsk;

   if ((r & rmbit) != 0)
     {
       if (r == rmbit)
   if ((r & rmbit) != 0)
     {
       if (r == rmbit)


@@ -1917,7 +2538,9 @@ emdnorm (s, lost, subflg, exp, rcntrl)
       eaddm (rbit, s);
     }
  mddone:
       eaddm (rbit, s);
     }
  mddone:

-  if ((rndprc < 64) && (exp <= 0))
+/* Undo the temporary shift for denormal values.  */
+  if ((exp <= 0) && (rndprc != NBITS)
+      && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))

     {
       eshup1 (s);
     }
     {
       eshup1 (s);
     }


@@ -1945,7 +2568,7 @@ emdnorm (s, lost, subflg, exp, rcntrl)
       for (i = M + 1; i < NI - 1; i++)
        s[i] = 0xffff;
       s[NI - 1] = 0;
       for (i = M + 1; i < NI - 1; i++)
        s[i] = 0xffff;
       s[NI - 1] = 0;

-      if (rndprc < 64)
+      if ((rndprc < 64) || (rndprc == 113))

        {
          s[rw] &= ~rmsk;
          if (rndprc == 24)
        {
          s[rw] &= ~rmsk;
          if (rndprc == 24)


@@ -1963,43 +2586,76 @@ emdnorm (s, lost, subflg, exp, rcntrl)
     s[1] = (unsigned EMUSHORT) exp;
 }
 
     s[1] = (unsigned EMUSHORT) exp;
 }
 

-
-
-/*
-;      Subtract external format numbers.
-;
-;      unsigned EMUSHORT a[NE], b[NE], c[NE];
-;      esub (a, b, c);  c = b - a
-*/
+/*  Subtract.  C = B - A, all e type numbers.  */

 
 static int subflg = 0;
 
 
 static int subflg = 0;
 

-void 
+static void 

 esub (a, b, c)
      unsigned EMUSHORT *a, *b, *c;
 {
 
 esub (a, b, c)
      unsigned EMUSHORT *a, *b, *c;
 {
 

+#ifdef NANS
+  if (eisnan (a))
+    {
+      emov (a, c);
+      return;
+    }
+  if (eisnan (b))
+    {
+      emov (b, c);
+      return;
+    }
+/* Infinity minus infinity is a NaN.
+   Test for subtracting infinities of the same sign.  */
+  if (eisinf (a) && eisinf (b)
+      && ((eisneg (a) ^ eisneg (b)) == 0))
+    {
+      mtherr ("esub", INVALID);
+      enan (c, 0);
+      return;
+    }
+#endif

   subflg = 1;
   eadd1 (a, b, c);
 }
 
   subflg = 1;
   eadd1 (a, b, c);
 }
 

+/* Add.  C = A + B, all e type.  */

 
 

-/*
-;      Add.
-;
-;      unsigned EMUSHORT a[NE], b[NE], c[NE];
-;      eadd (a, b, c);  c = b + a
-*/
-void 
+static void 

 eadd (a, b, c)
      unsigned EMUSHORT *a, *b, *c;
 {
 
 eadd (a, b, c)
      unsigned EMUSHORT *a, *b, *c;
 {
 

+#ifdef NANS
+/* NaN plus anything is a NaN.  */
+  if (eisnan (a))
+    {
+      emov (a, c);
+      return;
+    }
+  if (eisnan (b))
+    {
+      emov (b, c);
+      return;
+    }
+/* Infinity minus infinity is a NaN.
+   Test for adding infinities of opposite signs.  */
+  if (eisinf (a) && eisinf (b)
+      && ((eisneg (a) ^ eisneg (b)) != 0))
+    {
+      mtherr ("esub", INVALID);
+      enan (c, 0);
+      return;
+    }
+#endif

   subflg = 0;
   eadd1 (a, b, c);
 }
 
   subflg = 0;
   eadd1 (a, b, c);
 }
 

-void 
+/* Arithmetic common to both addition and subtraction.  */
+
+static void 

 eadd1 (a, b, c)
      unsigned EMUSHORT *a, *b, *c;
 {
 eadd1 (a, b, c)
      unsigned EMUSHORT *a, *b, *c;
 {


@@ -2059,7 +2715,7 @@ eadd1 (a, b, c)
              return;
            }
          /* if same sign, result is double */
              return;
            }
          /* if same sign, result is double */

-         /* double denomalized tiny number */
+         /* double denormalized tiny number */

          if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
            {
              eshup1 (bi);
          if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
            {
              eshup1 (bi);


@@ -2070,8 +2726,15 @@ eadd1 (a, b, c)
            {
              if (bi[j] != 0)
                {
            {
              if (bi[j] != 0)
                {

-                 /* This could overflow, but let emovo take care of that. */

                  ltb += 1;
                  ltb += 1;

+                 if (ltb >= 0x7fff)
+                   {
+                     eclear (c);
+                     if (ai[0] != 0)
+                       eneg (c);
+                     einfin (c);
+                     return;
+                   }

                  break;
                }
            }
                  break;
                }
            }


@@ -2101,36 +2764,53 @@ eadd1 (a, b, c)
   emovo (bi, c);
 }
 
   emovo (bi, c);
 }
 

+/* Divide: C = B/A, all e type.  */

 
 

-
-/*
-;      Divide.
-;
-;      unsigned EMUSHORT a[NE], b[NE], c[NE];
-;      ediv (a, b, c); c = b / a
-*/
-void 
+static void 

 ediv (a, b, c)
      unsigned EMUSHORT *a, *b, *c;
 {
   unsigned EMUSHORT ai[NI], bi[NI];
 ediv (a, b, c)
      unsigned EMUSHORT *a, *b, *c;
 {
   unsigned EMUSHORT ai[NI], bi[NI];

-  int i;
+  int i, sign;

   EMULONG lt, lta, ltb;
 
   EMULONG lt, lta, ltb;
 

+/* IEEE says if result is not a NaN, the sign is "-" if and only if
+   operands have opposite signs -- but flush -0 to 0 later if not IEEE.  */
+  sign = eisneg(a) ^ eisneg(b);
+
+#ifdef NANS
+/* Return any NaN input.  */
+  if (eisnan (a))
+    {
+    emov (a, c);
+    return;
+    }
+  if (eisnan (b))
+    {
+    emov (b, c);
+    return;
+    }
+/* Zero over zero, or infinity over infinity, is a NaN.  */
+  if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
+      || (eisinf (a) && eisinf (b)))
+    {
+    mtherr ("ediv", INVALID);
+    enan (c, sign);
+    return;
+    }
+#endif
+/* Infinity over anything else is infinity.  */

 #ifdef INFINITY
   if (eisinf (b))
     {
 #ifdef INFINITY
   if (eisinf (b))
     {

-      if (eisneg (a) ^ eisneg (b))
-       *(c + (NE - 1)) = 0x8000;
-      else
-       *(c + (NE - 1)) = 0;

       einfin (c);
       einfin (c);

-      return;
+      goto divsign;

     }
     }

+/* Anything else over infinity is zero.  */

   if (eisinf (a))
     {
       eclear (c);
   if (eisinf (a))
     {
       eclear (c);

-      return;
+      goto divsign;

     }
 #endif
   emovi (a, ai);
     }
 #endif
   emovi (a, ai);


@@ -2138,7 +2818,7 @@ ediv (a, b, c)
   lta = ai[E];
   ltb = bi[E];
   if (bi[E] == 0)
   lta = ai[E];
   ltb = bi[E];
   if (bi[E] == 0)

-    {                          /* See if numerator is zero. */
+    {                          /* See if numerator is zero.  */

       for (i = 1; i < NI - 1; i++)
        {
          if (bi[i] != 0)
       for (i = 1; i < NI - 1; i++)
        {
          if (bi[i] != 0)


@@ -2148,7 +2828,7 @@ ediv (a, b, c)
            }
        }
       eclear (c);
            }
        }
       eclear (c);

-      return;
+      goto divsign;

     }
  dnzro1:
 
     }
  dnzro1:
 


@@ -2162,13 +2842,11 @@ ediv (a, b, c)
              goto dnzro2;
            }
        }
              goto dnzro2;
            }
        }

-      if (ai[0] == bi[0])
-       *(c + (NE - 1)) = 0;
-      else
-       *(c + (NE - 1)) = 0x8000;
+/* Divide by zero is not an invalid operation.
+   It is a divide-by-zero operation!   */

       einfin (c);
       mtherr ("ediv", SING);
       einfin (c);
       mtherr ("ediv", SING);

-      return;
+      goto divsign;

     }
  dnzro2:
 
     }
  dnzro2:
 


@@ -2176,39 +2854,61 @@ ediv (a, b, c)
   /* calculate exponent */
   lt = ltb - lta + EXONE;
   emdnorm (bi, i, 0, lt, 64);
   /* calculate exponent */
   lt = ltb - lta + EXONE;
   emdnorm (bi, i, 0, lt, 64);

-  /* set the sign */
-  if (ai[0] == bi[0])
-    bi[0] = 0;
-  else
-    bi[0] = 0Xffff;

   emovo (bi, c);
   emovo (bi, c);

-}

 
 

+ divsign:
+
+  if (sign
+#ifndef IEEE
+      && (ecmp (c, ezero) != 0)
+#endif
+      )
+     *(c+(NE-1)) |= 0x8000;
+  else
+     *(c+(NE-1)) &= ~0x8000;
+}

 
 

+/* Multiply e-types A and B, return e-type product C.   */

 
 

-/*
-;      Multiply.
-;
-;      unsigned EMUSHORT a[NE], b[NE], c[NE];
-;      emul (a, b, c); c = b * a
-*/
-void 
+static void 

 emul (a, b, c)
      unsigned EMUSHORT *a, *b, *c;
 {
   unsigned EMUSHORT ai[NI], bi[NI];
 emul (a, b, c)
      unsigned EMUSHORT *a, *b, *c;
 {
   unsigned EMUSHORT ai[NI], bi[NI];

-  int i, j;
+  int i, j, sign;

   EMULONG lt, lta, ltb;
 
   EMULONG lt, lta, ltb;
 

+/* IEEE says if result is not a NaN, the sign is "-" if and only if
+   operands have opposite signs -- but flush -0 to 0 later if not IEEE.  */
+  sign = eisneg(a) ^ eisneg(b);
+
+#ifdef NANS
+/* NaN times anything is the same NaN.  */
+  if (eisnan (a))
+    {
+    emov (a, c);
+    return;
+    }
+  if (eisnan (b))
+    {
+    emov (b, c);
+    return;
+    }
+/* Zero times infinity is a NaN.  */
+  if ((eisinf (a) && (ecmp (b, ezero) == 0))
+      || (eisinf (b) && (ecmp (a, ezero) == 0)))
+    {
+    mtherr ("emul", INVALID);
+    enan (c, sign);
+    return;
+    }
+#endif
+/* Infinity times anything else is infinity.  */

 #ifdef INFINITY
   if (eisinf (a) || eisinf (b))
     {
 #ifdef INFINITY
   if (eisinf (a) || eisinf (b))
     {

-      if (eisneg (a) ^ eisneg (b))
-       *(c + (NE - 1)) = 0x8000;
-      else
-       *(c + (NE - 1)) = 0;

       einfin (c);
       einfin (c);

-      return;
+      goto mulsign;

     }
 #endif
   emovi (a, ai);
     }
 #endif
   emovi (a, ai);


@@ -2226,7 +2926,7 @@ emul (a, b, c)
            }
        }
       eclear (c);
            }
        }
       eclear (c);

-      return;
+      goto mulsign;

     }
  mnzer1:
 
     }
  mnzer1:
 


@@ -2241,7 +2941,7 @@ emul (a, b, c)
            }
        }
       eclear (c);
            }
        }
       eclear (c);

-      return;
+      goto mulsign;

     }
  mnzer2:
 
     }
  mnzer2:
 


@@ -2250,43 +2950,46 @@ emul (a, b, c)
   /* calculate exponent */
   lt = lta + ltb - (EXONE - 1);
   emdnorm (bi, j, 0, lt, 64);
   /* calculate exponent */
   lt = lta + ltb - (EXONE - 1);
   emdnorm (bi, j, 0, lt, 64);

-  /* calculate sign of product */
-  if (ai[0] == bi[0])
-    bi[0] = 0;
-  else
-    bi[0] = 0xffff;

   emovo (bi, c);
   emovo (bi, c);

-}

 
 

+ mulsign:

 
 

+  if (sign
+#ifndef IEEE
+      && (ecmp (c, ezero) != 0)
+#endif
+      )
+     *(c+(NE-1)) |= 0x8000;
+  else
+     *(c+(NE-1)) &= ~0x8000;
+}

 
 

+/* Convert double precision PE to e-type Y.  */

 
 

-/*
-; Convert IEEE double precision to e type
-;      double d;
-;      unsigned EMUSHORT x[N+2];
-;      e53toe (&d, x);
-*/
-void 
-e53toe (e, y)
-     unsigned EMUSHORT *e, *y;
+static void
+e53toe (pe, y)
+     unsigned EMUSHORT *pe, *y;

 {
 #ifdef DEC
 
 {
 #ifdef DEC
 

-  dectoe (e, y);               /* see etodec.c */
+  dectoe (pe, y);

 
 #else
 
 #else

+#ifdef IBM
+
+  ibmtoe (pe, y, DFmode);

 
 

+#else

   register unsigned EMUSHORT r;
   register unsigned EMUSHORT r;

-  register unsigned EMUSHORT *p;
+  register unsigned EMUSHORT *e, *p;

   unsigned EMUSHORT yy[NI];
   int denorm, k;
 
   unsigned EMUSHORT yy[NI];
   int denorm, k;
 

+  e = pe;

   denorm = 0;                  /* flag if denormalized number */
   ecleaz (yy);
   denorm = 0;                  /* flag if denormalized number */
   ecleaz (yy);

-#ifdef IBMPC
-  e += 3;
-#endif
+  if (! REAL_WORDS_BIG_ENDIAN)
+    e += 3;

   r = *e;
   yy[0] = 0;
   if (r & 0x8000)
   r = *e;
   yy[0] = 0;
   if (r & 0x8000)


@@ -2296,15 +2999,37 @@ e53toe (e, y)
 #ifdef INFINITY
   if (r == 0x7ff0)
     {
 #ifdef INFINITY
   if (r == 0x7ff0)
     {

+#ifdef NANS
+      if (! REAL_WORDS_BIG_ENDIAN)
+       {
+         if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
+             || (pe[1] != 0) || (pe[0] != 0))
+           {
+             enan (y, yy[0] != 0);
+             return;
+           }
+       }
+      else
+       {
+         if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
+             || (pe[2] != 0) || (pe[3] != 0))
+           {
+             enan (y, yy[0] != 0);
+             return;
+           }
+       }
+#endif  /* NANS */
+      eclear (y);

       einfin (y);
       einfin (y);

-      if (r & 0x8000)
+      if (yy[0])

        eneg (y);
       return;
     }
        eneg (y);
       return;
     }

-#endif
+#endif  /* INFINITY */

   r >>= 4;
   /* If zero exponent, then the significand is denormalized.
   r >>= 4;
   /* If zero exponent, then the significand is denormalized.

-   * So, take back the understood high significand bit. */
+     So take back the understood high significand bit.  */
+

   if (r == 0)
     {
       denorm = 1;
   if (r == 0)
     {
       denorm = 1;


@@ -2313,16 +3038,20 @@ e53toe (e, y)
   r += EXONE - 01777;
   yy[E] = r;
   p = &yy[M + 1];
   r += EXONE - 01777;
   yy[E] = r;
   p = &yy[M + 1];

-#ifdef IBMPC
-  *p++ = *(--e);
-  *p++ = *(--e);
-  *p++ = *(--e);
-#endif
-#ifdef MIEEE
-  ++e;
-  *p++ = *e++;
-  *p++ = *e++;
-  *p++ = *e++;
+#ifdef IEEE
+  if (! REAL_WORDS_BIG_ENDIAN)
+    {
+      *p++ = *(--e);
+      *p++ = *(--e);
+      *p++ = *(--e);
+    }
+  else
+    {
+      ++e;
+      *p++ = *e++;
+      *p++ = *e++;
+      *p++ = *e++;
+    }

 #endif
   eshift (yy, -5);
   if (denorm)
 #endif
   eshift (yy, -5);
   if (denorm)


@@ -2333,70 +3062,239 @@ e53toe (e, y)
        yy[E] -= (unsigned EMUSHORT) (k - 1);
     }
   emovo (yy, y);
        yy[E] -= (unsigned EMUSHORT) (k - 1);
     }
   emovo (yy, y);

+#endif /* not IBM */

 #endif /* not DEC */
 }
 
 #endif /* not DEC */
 }
 

-void 
-e64toe (e, y)
-     unsigned EMUSHORT *e, *y;
+/* Convert double extended precision float PE to e type Y.  */
+
+static void 
+e64toe (pe, y)
+     unsigned EMUSHORT *pe, *y;

 {
   unsigned EMUSHORT yy[NI];
 {
   unsigned EMUSHORT yy[NI];

-  unsigned EMUSHORT *p, *q;
+  unsigned EMUSHORT *e, *p, *q;

   int i;
 
   int i;
 

+  e = pe;

   p = yy;
   for (i = 0; i < NE - 5; i++)
     *p++ = 0;
   p = yy;
   for (i = 0; i < NE - 5; i++)
     *p++ = 0;

-#ifdef IBMPC
-  for (i = 0; i < 5; i++)
-    *p++ = *e++;
-#endif
+/* This precision is not ordinarily supported on DEC or IBM.  */

 #ifdef DEC
   for (i = 0; i < 5; i++)
     *p++ = *e++;
 #endif
 #ifdef DEC
   for (i = 0; i < 5; i++)
     *p++ = *e++;
 #endif

-#ifdef MIEEE
+#ifdef IBM

   p = &yy[0] + (NE - 1);
   *p-- = *e++;
   ++e;
   p = &yy[0] + (NE - 1);
   *p-- = *e++;
   ++e;

-  for (i = 0; i < 4; i++)
+  for (i = 0; i < 5; i++)

     *p-- = *e++;
 #endif
     *p-- = *e++;
 #endif

-  p = yy;
-  q = y;
-#ifdef INFINITY
-  if (*p == 0x7fff)
+#ifdef IEEE
+  if (! REAL_WORDS_BIG_ENDIAN)

     {
     {

-      einfin (y);
-      if (*p & 0x8000)
-       eneg (y);
-      return;
-    }
-#endif
-  for (i = 0; i < NE; i++)
+      for (i = 0; i < 5; i++)
+       *p++ = *e++;
+
+      /* For denormal long double Intel format, shift significand up one
+        -- but only if the top significand bit is zero.  A top bit of 1
+        is "pseudodenormal" when the exponent is zero.  */
+      if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
+       {
+         unsigned EMUSHORT temp[NI];
+
+         emovi(yy, temp);
+         eshup1(temp);
+         emovo(temp,y);
+         return;
+       }
+    }
+  else
+    {
+      p = &yy[0] + (NE - 1);
+#ifdef ARM_EXTENDED_IEEE_FORMAT
+      /* For ARMs, the exponent is in the lowest 15 bits of the word.  */
+      *p-- = (e[0] & 0x8000) | (e[1] & 0x7ffff);
+      e += 2;
+#else
+      *p-- = *e++;
+      ++e;
+#endif
+      for (i = 0; i < 4; i++)
+       *p-- = *e++;
+    }
+#endif
+#ifdef INFINITY
+  /* Point to the exponent field and check max exponent cases.  */
+  p = &yy[NE - 1];
+  if ((*p & 0x7fff) == 0x7fff)
+    {
+#ifdef NANS
+      if (! REAL_WORDS_BIG_ENDIAN)
+       {
+         for (i = 0; i < 4; i++)
+           {
+             if ((i != 3 && pe[i] != 0)
+                 /* Anything but 0x8000 here, including 0, is a NaN.  */
+                 || (i == 3 && pe[i] != 0x8000))
+               {
+                 enan (y, (*p & 0x8000) != 0);
+                 return;
+               }
+           }
+       }
+      else
+       {
+#ifdef ARM_EXTENDED_IEEE_FORMAT
+         for (i = 2; i <= 5; i++)
+           {
+             if (pe[i] != 0)
+               {
+                 enan (y, (*p & 0x8000) != 0);
+                 return;
+               }
+           }
+#else /* not ARM */
+         /* In Motorola extended precision format, the most significant
+            bit of an infinity mantissa could be either 1 or 0.  It is
+            the lower order bits that tell whether the value is a NaN.  */
+         if ((pe[2] & 0x7fff) != 0)
+           goto bigend_nan;
+
+         for (i = 3; i <= 5; i++)
+           {
+             if (pe[i] != 0)
+               {
+bigend_nan:
+                 enan (y, (*p & 0x8000) != 0);
+                 return;
+               }
+           }
+#endif /* not ARM */
+       }
+#endif /* NANS */
+      eclear (y);
+      einfin (y);
+      if (*p & 0x8000)
+       eneg (y);
+      return;
+    }
+#endif  /* INFINITY */
+  p = yy;
+  q = y;
+  for (i = 0; i < NE; i++)

     *q++ = *p++;
 }
 
     *q++ = *p++;
 }
 

+/* Convert 128-bit long double precision float PE to e type Y.  */

 
 

-/*
-; Convert IEEE single precision to e type
-;      float d;
-;      unsigned EMUSHORT x[N+2];
-;      dtox (&d, x);
-*/
-void 
-e24toe (e, y)
-     unsigned EMUSHORT *e, *y;
+static void 
+e113toe (pe, y)
+     unsigned EMUSHORT *pe, *y;

 {
   register unsigned EMUSHORT r;
 {
   register unsigned EMUSHORT r;

-  register unsigned EMUSHORT *p;
+  unsigned EMUSHORT *e, *p;
+  unsigned EMUSHORT yy[NI];
+  int denorm, i;
+
+  e = pe;
+  denorm = 0;
+  ecleaz (yy);
+#ifdef IEEE
+  if (! REAL_WORDS_BIG_ENDIAN)
+    e += 7;
+#endif
+  r = *e;
+  yy[0] = 0;
+  if (r & 0x8000)
+    yy[0] = 0xffff;
+  r &= 0x7fff;
+#ifdef INFINITY
+  if (r == 0x7fff)
+    {
+#ifdef NANS
+      if (! REAL_WORDS_BIG_ENDIAN)
+       {
+         for (i = 0; i < 7; i++)
+           {
+             if (pe[i] != 0)
+               {
+                 enan (y, yy[0] != 0);
+                 return;
+               }
+           }
+       }
+      else
+       {
+         for (i = 1; i < 8; i++)
+           {
+             if (pe[i] != 0)
+               {
+                 enan (y, yy[0] != 0);
+                 return;
+               }
+           }
+       }
+#endif /* NANS */
+      eclear (y);
+      einfin (y);
+      if (yy[0])
+       eneg (y);
+      return;
+    }
+#endif  /* INFINITY */
+  yy[E] = r;
+  p = &yy[M + 1];
+#ifdef IEEE
+  if (! REAL_WORDS_BIG_ENDIAN)
+    {
+      for (i = 0; i < 7; i++)
+       *p++ = *(--e);
+    }
+  else
+    {
+      ++e;
+      for (i = 0; i < 7; i++)
+       *p++ = *e++;
+    }
+#endif
+/* If denormal, remove the implied bit; else shift down 1.  */
+  if (r == 0)
+    {
+      yy[M] = 0;
+    }
+  else
+    {
+      yy[M] = 1;
+      eshift (yy, -1);
+    }
+  emovo (yy, y);
+}
+
+/* Convert single precision float PE to e type Y.  */
+
+static void 
+e24toe (pe, y)
+     unsigned EMUSHORT *pe, *y;
+{
+#ifdef IBM
+
+  ibmtoe (pe, y, SFmode);
+
+#else
+  register unsigned EMUSHORT r;
+  register unsigned EMUSHORT *e, *p;

   unsigned EMUSHORT yy[NI];
   int denorm, k;
 
   unsigned EMUSHORT yy[NI];
   int denorm, k;
 

+  e = pe;

   denorm = 0;                  /* flag if denormalized number */
   ecleaz (yy);
   denorm = 0;                  /* flag if denormalized number */
   ecleaz (yy);

-#ifdef IBMPC
-  e += 1;
+#ifdef IEEE
+  if (! REAL_WORDS_BIG_ENDIAN)
+    e += 1;

 #endif
 #ifdef DEC
   e += 1;
 #endif
 #ifdef DEC
   e += 1;


@@ -2410,15 +3308,34 @@ e24toe (e, y)
 #ifdef INFINITY
   if (r == 0x7f80)
     {
 #ifdef INFINITY
   if (r == 0x7f80)
     {

+#ifdef NANS
+      if (REAL_WORDS_BIG_ENDIAN)
+       {
+         if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
+           {
+             enan (y, yy[0] != 0);
+             return;
+           }
+       }
+      else
+       {
+         if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
+           {
+             enan (y, yy[0] != 0);
+             return;
+           }
+       }
+#endif  /* NANS */
+      eclear (y);

       einfin (y);
       einfin (y);

-      if (r & 0x8000)
+      if (yy[0])

        eneg (y);
       return;
     }
        eneg (y);
       return;
     }

-#endif
+#endif  /* INFINITY */

   r >>= 7;
   /* If zero exponent, then the significand is denormalized.
   r >>= 7;
   /* If zero exponent, then the significand is denormalized.

-   * So, take back the understood high significand bit. */
+     So take back the understood high significand bit.  */

   if (r == 0)
     {
       denorm = 1;
   if (r == 0)
     {
       denorm = 1;


@@ -2427,15 +3344,17 @@ e24toe (e, y)
   r += EXONE - 0177;
   yy[E] = r;
   p = &yy[M + 1];
   r += EXONE - 0177;
   yy[E] = r;
   p = &yy[M + 1];

-#ifdef IBMPC
-  *p++ = *(--e);
-#endif

 #ifdef DEC
   *p++ = *(--e);
 #endif
 #ifdef DEC
   *p++ = *(--e);
 #endif

-#ifdef MIEEE
-  ++e;
-  *p++ = *e++;
+#ifdef IEEE
+  if (! REAL_WORDS_BIG_ENDIAN)
+    *p++ = *(--e);
+  else
+    {
+      ++e;
+      *p++ = *e++;
+    }

 #endif
   eshift (yy, -8);
   if (denorm)
 #endif
   eshift (yy, -8);
   if (denorm)


@@ -2446,10 +3365,103 @@ e24toe (e, y)
        yy[E] -= (unsigned EMUSHORT) (k - 1);
     }
   emovo (yy, y);
        yy[E] -= (unsigned EMUSHORT) (k - 1);
     }
   emovo (yy, y);

+#endif /* not IBM */

 }
 
 }
 

+/* Convert e-type X to IEEE 128-bit long double format E.  */

 
 

-void 
+static void 
+etoe113 (x, e)
+     unsigned EMUSHORT *x, *e;
+{
+  unsigned EMUSHORT xi[NI];
+  EMULONG exp;
+  int rndsav;
+
+#ifdef NANS
+  if (eisnan (x))
+    {
+      make_nan (e, eisneg (x), TFmode);
+      return;
+    }
+#endif
+  emovi (x, xi);
+  exp = (EMULONG) xi[E];
+#ifdef INFINITY
+  if (eisinf (x))
+    goto nonorm;
+#endif
+  /* round off to nearest or even */
+  rndsav = rndprc;
+  rndprc = 113;
+  emdnorm (xi, 0, 0, exp, 64);
+  rndprc = rndsav;
+ nonorm:
+  toe113 (xi, e);
+}
+
+/* Convert exploded e-type X, that has already been rounded to
+   113-bit precision, to IEEE 128-bit long double format Y.  */
+
+static void 
+toe113 (a, b)
+     unsigned EMUSHORT *a, *b;
+{
+  register unsigned EMUSHORT *p, *q;
+  unsigned EMUSHORT i;
+
+#ifdef NANS
+  if (eiisnan (a))
+    {
+      make_nan (b, eiisneg (a), TFmode);
+      return;
+    }
+#endif
+  p = a;
+  if (REAL_WORDS_BIG_ENDIAN)
+    q = b;
+  else
+    q = b + 7;                 /* point to output exponent */
+
+  /* If not denormal, delete the implied bit.  */
+  if (a[E] != 0)
+    {
+      eshup1 (a);
+    }
+  /* combine sign and exponent */
+  i = *p++;
+  if (REAL_WORDS_BIG_ENDIAN)
+    {
+      if (i)
+       *q++ = *p++ | 0x8000;
+      else
+       *q++ = *p++;
+    }
+  else
+    {
+      if (i)
+       *q-- = *p++ | 0x8000;
+      else
+       *q-- = *p++;
+    }
+  /* skip over guard word */
+  ++p;
+  /* move the significand */
+  if (REAL_WORDS_BIG_ENDIAN)
+    {
+      for (i = 0; i < 7; i++)
+       *q++ = *p++;
+    }
+  else
+    {
+      for (i = 0; i < 7; i++)
+       *q-- = *p++;
+    }
+}
+
+/* Convert e-type X to IEEE double extended format E.  */
+
+static void 

 etoe64 (x, e)
      unsigned EMUSHORT *x, *e;
 {
 etoe64 (x, e)
      unsigned EMUSHORT *x, *e;
 {


@@ -2457,6 +3469,13 @@ etoe64 (x, e)
   EMULONG exp;
   int rndsav;
 
   EMULONG exp;
   int rndsav;
 

+#ifdef NANS
+  if (eisnan (x))
+    {
+      make_nan (e, eisneg (x), XFmode);
+      return;
+    }
+#endif

   emovi (x, xi);
   /* adjust exponent for offset */
   exp = (EMULONG) xi[E];
   emovi (x, xi);
   /* adjust exponent for offset */
   exp = (EMULONG) xi[E];


@@ -2473,7 +3492,9 @@ etoe64 (x, e)
   toe64 (xi, e);
 }
 
   toe64 (xi, e);
 }
 

-/* move out internal format to ieee long double */
+/* Convert exploded e-type X, that has already been rounded to
+   64-bit precision, to IEEE double extended format Y.  */
+

 static void 
 toe64 (a, b)
      unsigned EMUSHORT *a, *b;
 static void 
 toe64 (a, b)
      unsigned EMUSHORT *a, *b;


@@ -2481,60 +3502,125 @@ toe64 (a, b)
   register unsigned EMUSHORT *p, *q;
   unsigned EMUSHORT i;
 
   register unsigned EMUSHORT *p, *q;
   unsigned EMUSHORT i;
 

+#ifdef NANS
+  if (eiisnan (a))
+    {
+      make_nan (b, eiisneg (a), XFmode);
+      return;
+    }
+#endif
+  /* Shift denormal long double Intel format significand down one bit.  */
+  if ((a[E] == 0) && ! REAL_WORDS_BIG_ENDIAN)
+    eshdn1 (a);

   p = a;
   p = a;

-#ifdef MIEEE
+#ifdef IBM

   q = b;
   q = b;

-#else
-  q = b + 4;                   /* point to output exponent */
+#endif
+#ifdef DEC
+  q = b + 4;
+#endif
+#ifdef IEEE
+  if (REAL_WORDS_BIG_ENDIAN)
+    q = b;
+  else
+    {
+      q = b + 4;                       /* point to output exponent */

 #if LONG_DOUBLE_TYPE_SIZE == 96
 #if LONG_DOUBLE_TYPE_SIZE == 96

-  /* Clear the last two bytes of 12-byte Intel format */
-  *(q+1) = 0;
+      /* Clear the last two bytes of 12-byte Intel format */
+      *(q+1) = 0;

 #endif
 #endif

+    }

 #endif
 
   /* combine sign and exponent */
   i = *p++;
 #endif
 
   /* combine sign and exponent */
   i = *p++;

-#ifdef MIEEE
+#ifdef IBM

   if (i)
     *q++ = *p++ | 0x8000;
   else
     *q++ = *p++;
   *q++ = 0;
   if (i)
     *q++ = *p++ | 0x8000;
   else
     *q++ = *p++;
   *q++ = 0;

-#else
+#endif
+#ifdef DEC

   if (i)
     *q-- = *p++ | 0x8000;
   else
     *q-- = *p++;
 #endif
   if (i)
     *q-- = *p++ | 0x8000;
   else
     *q-- = *p++;
 #endif

+#ifdef IEEE
+  if (REAL_WORDS_BIG_ENDIAN)
+    {
+#ifdef ARM_EXTENDED_IEEE_FORMAT
+      /* The exponent is in the lowest 15 bits of the first word.  */
+      *q++ = i ? 0x8000 : 0;
+      *q++ = *p++;
+#else
+      if (i)
+       *q++ = *p++ | 0x8000;
+      else
+       *q++ = *p++;
+      *q++ = 0;
+#endif
+    }
+  else
+    {
+      if (i)
+       *q-- = *p++ | 0x8000;
+      else
+       *q-- = *p++;
+    }
+#endif

   /* skip over guard word */
   ++p;
   /* move the significand */
   /* skip over guard word */
   ++p;
   /* move the significand */

-#ifdef MIEEE
+#ifdef IBM

   for (i = 0; i < 4; i++)
     *q++ = *p++;
   for (i = 0; i < 4; i++)
     *q++ = *p++;

-#else
+#endif
+#ifdef DEC

   for (i = 0; i < 4; i++)
     *q-- = *p++;
 #endif
   for (i = 0; i < 4; i++)
     *q-- = *p++;
 #endif

+#ifdef IEEE
+  if (REAL_WORDS_BIG_ENDIAN)
+    {
+      for (i = 0; i < 4; i++)
+       *q++ = *p++;
+    }
+  else
+    {
+#ifdef INFINITY
+      if (eiisinf (a))
+       {
+         /* Intel long double infinity significand.  */
+         *q-- = 0x8000;
+         *q-- = 0;
+         *q-- = 0;
+         *q = 0;
+         return;
+       }
+#endif
+      for (i = 0; i < 4; i++)
+       *q-- = *p++;
+    }
+#endif

 }
 
 }
 

-
-/*
-; e type to IEEE double precision
-;      double d;
-;      unsigned EMUSHORT x[NE];
-;      etoe53 (x, &d);
-*/
+/* e type to double precision.  */

 
 #ifdef DEC
 
 #ifdef DEC

+/* Convert e-type X to DEC-format double E.  */

 
 

-void 
+static void 

 etoe53 (x, e)
      unsigned EMUSHORT *x, *e;
 {
   etodec (x, e);               /* see etodec.c */
 }
 
 etoe53 (x, e)
      unsigned EMUSHORT *x, *e;
 {
   etodec (x, e);               /* see etodec.c */
 }
 

+/* Convert exploded e-type X, that has already been rounded to
+   56-bit double precision, to DEC double Y.  */
+

 static void 
 toe53 (x, y)
      unsigned EMUSHORT *x, *y;
 static void 
 toe53 (x, y)
      unsigned EMUSHORT *x, *y;


@@ -2543,8 +3629,31 @@ toe53 (x, y)
 }
 
 #else
 }
 
 #else

+#ifdef IBM
+/* Convert e-type X to IBM 370-format double E.  */

 
 

-void 
+static void 
+etoe53 (x, e)
+     unsigned EMUSHORT *x, *e;
+{
+  etoibm (x, e, DFmode);
+}
+
+/* Convert exploded e-type X, that has already been rounded to
+   56-bit precision, to IBM 370 double Y.  */
+
+static void 
+toe53 (x, y)
+     unsigned EMUSHORT *x, *y;
+{
+  toibm (x, y, DFmode);
+}
+
+#else  /* it's neither DEC nor IBM */
+
+/* Convert e-type X to IEEE double E.  */
+
+static void 

 etoe53 (x, e)
      unsigned EMUSHORT *x, *e;
 {
 etoe53 (x, e)
      unsigned EMUSHORT *x, *e;
 {


@@ -2552,6 +3661,13 @@ etoe53 (x, e)
   EMULONG exp;
   int rndsav;
 
   EMULONG exp;
   int rndsav;
 

+#ifdef NANS
+  if (eisnan (x))
+    {
+      make_nan (e, eisneg (x), DFmode);
+      return;
+    }
+#endif

   emovi (x, xi);
   /* adjust exponent for offsets */
   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
   emovi (x, xi);
   /* adjust exponent for offsets */
   exp = (EMULONG) xi[E] - (EXONE - 0x3ff);


@@ -2568,6 +3684,8 @@ etoe53 (x, e)
   toe53 (xi, e);
 }
 
   toe53 (xi, e);
 }
 

+/* Convert exploded e-type X, that has already been rounded to
+   53-bit precision, to IEEE double Y.  */

 
 static void 
 toe53 (x, y)
 
 static void 
 toe53 (x, y)


@@ -2576,10 +3694,17 @@ toe53 (x, y)
   unsigned EMUSHORT i;
   unsigned EMUSHORT *p;
 
   unsigned EMUSHORT i;
   unsigned EMUSHORT *p;
 

-
+#ifdef NANS
+  if (eiisnan (x))
+    {
+      make_nan (y, eiisneg (x), DFmode);
+      return;
+    }
+#endif

   p = &x[0];
   p = &x[0];

-#ifdef IBMPC
-  y += 3;
+#ifdef IEEE
+  if (! REAL_WORDS_BIG_ENDIAN)
+    y += 3;

 #endif
   *y = 0;                      /* output high order */
   if (*p++)
 #endif
   *y = 0;                      /* output high order */
   if (*p++)


@@ -2587,33 +3712,38 @@ toe53 (x, y)
 
   i = *p++;
   if (i >= (unsigned int) 2047)
 
   i = *p++;
   if (i >= (unsigned int) 2047)

-    {                          /* Saturate at largest number less than infinity. */
+    {
+      /* Saturate at largest number less than infinity.  */

 #ifdef INFINITY
       *y |= 0x7ff0;
 #ifdef INFINITY
       *y |= 0x7ff0;

-#ifdef IBMPC
-      *(--y) = 0;
-      *(--y) = 0;
-      *(--y) = 0;
-#endif
-#ifdef MIEEE
-      ++y;
-      *y++ = 0;
-      *y++ = 0;
-      *y++ = 0;
-#endif
+      if (! REAL_WORDS_BIG_ENDIAN)
+       {
+         *(--y) = 0;
+         *(--y) = 0;
+         *(--y) = 0;
+       }
+      else
+       {
+         ++y;
+         *y++ = 0;
+         *y++ = 0;
+         *y++ = 0;
+       }

 #else
       *y |= (unsigned EMUSHORT) 0x7fef;
 #else
       *y |= (unsigned EMUSHORT) 0x7fef;

-#ifdef IBMPC
-      *(--y) = 0xffff;
-      *(--y) = 0xffff;
-      *(--y) = 0xffff;
-#endif
-#ifdef MIEEE
-      ++y;
-      *y++ = 0xffff;
-      *y++ = 0xffff;
-      *y++ = 0xffff;
-#endif
+      if (! REAL_WORDS_BIG_ENDIAN)
+       {
+         *(--y) = 0xffff;
+         *(--y) = 0xffff;
+         *(--y) = 0xffff;
+       }
+      else
+       {
+         ++y;
+         *y++ = 0xffff;
+         *y++ = 0xffff;
+         *y++ = 0xffff;
+       }

 #endif
       return;
     }
 #endif
       return;
     }


@@ -2628,30 +3758,52 @@ toe53 (x, y)
     }
   i |= *p++ & (unsigned EMUSHORT) 0x0f;        /* *p = xi[M] */
   *y |= (unsigned EMUSHORT) i; /* high order output already has sign bit set */
     }
   i |= *p++ & (unsigned EMUSHORT) 0x0f;        /* *p = xi[M] */
   *y |= (unsigned EMUSHORT) i; /* high order output already has sign bit set */

-#ifdef IBMPC
-  *(--y) = *p++;
-  *(--y) = *p++;
-  *(--y) = *p;
-#endif
-#ifdef MIEEE
-  ++y;
-  *y++ = *p++;
-  *y++ = *p++;
-  *y++ = *p++;
-#endif
+  if (! REAL_WORDS_BIG_ENDIAN)
+    {
+      *(--y) = *p++;
+      *(--y) = *p++;
+      *(--y) = *p;
+    }
+  else
+    {
+      ++y;
+      *y++ = *p++;
+      *y++ = *p++;
+      *y++ = *p++;
+    }

 }
 
 }
 

+#endif /* not IBM */

 #endif /* not DEC */
 
 
 
 #endif /* not DEC */
 
 
 

-/*
-; e type to IEEE single precision
-;      float d;
-;      unsigned EMUSHORT x[N+2];
-;      xtod (x, &d);
-*/
-void 
+/* e type to single precision.  */
+
+#ifdef IBM
+/* Convert e-type X to IBM 370 float E.  */
+
+static void 
+etoe24 (x, e)
+     unsigned EMUSHORT *x, *e;
+{
+  etoibm (x, e, SFmode);
+}
+
+/* Convert exploded e-type X, that has already been rounded to
+   float precision, to IBM 370 float Y.  */
+
+static void 
+toe24 (x, y)
+     unsigned EMUSHORT *x, *y;
+{
+  toibm (x, y, SFmode);
+}
+
+#else
+/* Convert e-type X to IEEE float E.  DEC float is the same as IEEE float.  */
+
+static void 

 etoe24 (x, e)
      unsigned EMUSHORT *x, *e;
 {
 etoe24 (x, e)
      unsigned EMUSHORT *x, *e;
 {


@@ -2659,6 +3811,13 @@ etoe24 (x, e)
   unsigned EMUSHORT xi[NI];
   int rndsav;
 
   unsigned EMUSHORT xi[NI];
   int rndsav;
 

+#ifdef NANS
+  if (eisnan (x))
+    {
+      make_nan (e, eisneg (x), SFmode);
+      return;
+    }
+#endif

   emovi (x, xi);
   /* adjust exponent for offsets */
   exp = (EMULONG) xi[E] - (EXONE - 0177);
   emovi (x, xi);
   /* adjust exponent for offsets */
   exp = (EMULONG) xi[E] - (EXONE - 0177);


@@ -2675,6 +3834,9 @@ etoe24 (x, e)
   toe24 (xi, e);
 }
 
   toe24 (xi, e);
 }
 

+/* Convert exploded e-type X, that has already been rounded to
+   float precision, to IEEE float Y.  */
+

 static void 
 toe24 (x, y)
      unsigned EMUSHORT *x, *y;
 static void 
 toe24 (x, y)
      unsigned EMUSHORT *x, *y;


@@ -2682,9 +3844,17 @@ toe24 (x, y)
   unsigned EMUSHORT i;
   unsigned EMUSHORT *p;
 
   unsigned EMUSHORT i;
   unsigned EMUSHORT *p;
 

+#ifdef NANS
+  if (eiisnan (x))
+    {
+      make_nan (y, eiisneg (x), SFmode);
+      return;
+    }
+#endif

   p = &x[0];
   p = &x[0];

-#ifdef IBMPC
-  y += 1;
+#ifdef IEEE
+  if (! REAL_WORDS_BIG_ENDIAN)
+    y += 1;

 #endif
 #ifdef DEC
   y += 1;
 #endif
 #ifdef DEC
   y += 1;


@@ -2694,32 +3864,36 @@ toe24 (x, y)
     *y = 0x8000;               /* output sign bit */
 
   i = *p++;
     *y = 0x8000;               /* output sign bit */
 
   i = *p++;

-/* Handle overflow cases. */
+/* Handle overflow cases.  */

   if (i >= 255)
     {
 #ifdef INFINITY
       *y |= (unsigned EMUSHORT) 0x7f80;
   if (i >= 255)
     {
 #ifdef INFINITY
       *y |= (unsigned EMUSHORT) 0x7f80;

-#ifdef IBMPC
-      *(--y) = 0;
-#endif

 #ifdef DEC
       *(--y) = 0;
 #endif
 #ifdef DEC
       *(--y) = 0;
 #endif

-#ifdef MIEEE
-      ++y;
-      *y = 0;
+#ifdef IEEE
+      if (! REAL_WORDS_BIG_ENDIAN)
+       *(--y) = 0;
+      else
+       {
+         ++y;
+         *y = 0;
+       }

 #endif
 #else  /* no INFINITY */
       *y |= (unsigned EMUSHORT) 0x7f7f;
 #endif
 #else  /* no INFINITY */
       *y |= (unsigned EMUSHORT) 0x7f7f;

-#ifdef IBMPC
-      *(--y) = 0xffff;
-#endif

 #ifdef DEC
       *(--y) = 0xffff;
 #endif
 #ifdef DEC
       *(--y) = 0xffff;
 #endif

-#ifdef MIEEE
-      ++y;
-      *y = 0xffff;
+#ifdef IEEE
+      if (! REAL_WORDS_BIG_ENDIAN)
+       *(--y) = 0xffff;
+      else
+       {
+         ++y;
+         *y = 0xffff;
+       }

 #endif
 #ifdef ERANGE
       errno = ERANGE;
 #endif
 #ifdef ERANGE
       errno = ERANGE;


@@ -2737,30 +3911,30 @@ toe24 (x, y)
       eshift (x, 8);
     }
   i |= *p++ & (unsigned EMUSHORT) 0x7f;        /* *p = xi[M] */
       eshift (x, 8);
     }
   i |= *p++ & (unsigned EMUSHORT) 0x7f;        /* *p = xi[M] */

-  *y |= i;                     /* high order output already has sign bit set */
-#ifdef IBMPC
-  *(--y) = *p;
-#endif
+  /* High order output already has sign bit set.  */
+  *y |= i;

 #ifdef DEC
   *(--y) = *p;
 #endif
 #ifdef DEC
   *(--y) = *p;
 #endif

-#ifdef MIEEE
-  ++y;
-  *y = *p;
+#ifdef IEEE
+  if (! REAL_WORDS_BIG_ENDIAN)
+    *(--y) = *p;
+  else
+    {
+      ++y;
+      *y = *p;
+    }

 #endif
 }
 #endif
 }

+#endif  /* not IBM */

 
 

+/* Compare two e type numbers. 
+   Return +1 if a > b
+           0 if a == b
+          -1 if a < b
+          -2 if either a or b is a NaN.  */

 
 

-/* Compare two e type numbers.
- *
- * unsigned EMUSHORT a[NE], b[NE];
- * ecmp (a, b);
- *
- *  returns +1 if a > b
- *           0 if a == b
- *          -1 if a < b
- */
-int 
+static int 

 ecmp (a, b)
      unsigned EMUSHORT *a, *b;
 {
 ecmp (a, b)
      unsigned EMUSHORT *a, *b;
 {


@@ -2769,6 +3943,10 @@ ecmp (a, b)
   register int i;
   int msign;
 
   register int i;
   int msign;
 

+#ifdef NANS
+  if (eisnan (a)  || eisnan (b))
+      return (-2);
+#endif

   emovi (a, ai);
   p = ai;
   emovi (b, bi);
   emovi (a, ai);
   p = ai;
   emovi (b, bi);


@@ -2808,8 +3986,6 @@ ecmp (a, b)
 
   return (0);                  /* equality */
 
 
   return (0);                  /* equality */
 

-
-

  diff:
 
   if (*(--p) > *(--q))
  diff:
 
   if (*(--p) > *(--q))


@@ -2818,15 +3994,9 @@ ecmp (a, b)
     return (-msign);           /* p is littler */
 }
 
     return (-msign);           /* p is littler */
 }
 

+/* Find e-type nearest integer to X, as floor (X + 0.5).  */

 
 

-
-
-/* Find nearest integer to x = floor (x + 0.5)
- *
- * unsigned EMUSHORT x[NE], y[NE]
- * eround (x, y);
- */
-void 
+static void 

 eround (x, y)
      unsigned EMUSHORT *x, *y;
 {
 eround (x, y)
      unsigned EMUSHORT *x, *y;
 {


@@ -2834,42 +4004,41 @@ eround (x, y)
   efloor (y, y);
 }
 
   efloor (y, y);
 }
 

+/* Convert HOST_WIDE_INT LP to e type Y.  */

 
 

-
-
-/*
-; convert long integer to e type
-;
-;      long l;
-;      unsigned EMUSHORT x[NE];
-;      ltoe (&l, x);
-; note &l is the memory address of l
-*/
-void 
+static void 

 ltoe (lp, y)
 ltoe (lp, y)

-     long *lp;                 /* lp is the memory address of a long integer */
-     unsigned EMUSHORT *y;             /* y is the address of a short */
+     HOST_WIDE_INT *lp;
+     unsigned EMUSHORT *y;

 {
   unsigned EMUSHORT yi[NI];
 {
   unsigned EMUSHORT yi[NI];

-  unsigned long ll;
+  unsigned HOST_WIDE_INT ll;

   int k;
 
   ecleaz (yi);
   if (*lp < 0)
     {
       /* make it positive */
   int k;
 
   ecleaz (yi);
   if (*lp < 0)
     {
       /* make it positive */

-      ll = (unsigned long) (-(*lp));
+      ll = (unsigned HOST_WIDE_INT) (-(*lp));

       yi[0] = 0xffff;          /* put correct sign in the e type number */
     }
   else
     {
       yi[0] = 0xffff;          /* put correct sign in the e type number */
     }
   else
     {

-      ll = (unsigned long) (*lp);
+      ll = (unsigned HOST_WIDE_INT) (*lp);

     }
   /* move the long integer to yi significand area */
     }
   /* move the long integer to yi significand area */

+#if HOST_BITS_PER_WIDE_INT == 64
+  yi[M] = (unsigned EMUSHORT) (ll >> 48);
+  yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
+  yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
+  yi[M + 3] = (unsigned EMUSHORT) ll;
+  yi[E] = EXONE + 47;          /* exponent if normalize shift count were 0 */
+#else

   yi[M] = (unsigned EMUSHORT) (ll >> 16);
   yi[M + 1] = (unsigned EMUSHORT) ll;
   yi[M] = (unsigned EMUSHORT) (ll >> 16);
   yi[M + 1] = (unsigned EMUSHORT) ll;

-

   yi[E] = EXONE + 15;          /* exponent if normalize shift count were 0 */
   yi[E] = EXONE + 15;          /* exponent if normalize shift count were 0 */

+#endif
+

   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
     ecleaz (yi);               /* it was zero */
   else
   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
     ecleaz (yi);               /* it was zero */
   else


@@ -2877,31 +4046,33 @@ ltoe (lp, y)
   emovo (yi, y);               /* output the answer */
 }
 
   emovo (yi, y);               /* output the answer */
 }
 

-/*
-; convert unsigned long integer to e type
-;
-;      unsigned long l;
-;      unsigned EMUSHORT x[NE];
-;      ltox (&l, x);
-; note &l is the memory address of l
-*/
-void 
+/* Convert unsigned HOST_WIDE_INT LP to e type Y.  */
+
+static void 

 ultoe (lp, y)
 ultoe (lp, y)

-     unsigned long *lp;                /* lp is the memory address of a long integer */
-     unsigned EMUSHORT *y;             /* y is the address of a short */
+     unsigned HOST_WIDE_INT *lp;
+     unsigned EMUSHORT *y;

 {
   unsigned EMUSHORT yi[NI];
 {
   unsigned EMUSHORT yi[NI];

-  unsigned long ll;
+  unsigned HOST_WIDE_INT ll;

   int k;
 
   ecleaz (yi);
   ll = *lp;
 
   /* move the long integer to ayi significand area */
   int k;
 
   ecleaz (yi);
   ll = *lp;
 
   /* move the long integer to ayi significand area */

+#if HOST_BITS_PER_WIDE_INT == 64
+  yi[M] = (unsigned EMUSHORT) (ll >> 48);
+  yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
+  yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
+  yi[M + 3] = (unsigned EMUSHORT) ll;
+  yi[E] = EXONE + 47;          /* exponent if normalize shift count were 0 */
+#else

   yi[M] = (unsigned EMUSHORT) (ll >> 16);
   yi[M + 1] = (unsigned EMUSHORT) ll;
   yi[M] = (unsigned EMUSHORT) (ll >> 16);
   yi[M + 1] = (unsigned EMUSHORT) ll;

-

   yi[E] = EXONE + 15;          /* exponent if normalize shift count were 0 */
   yi[E] = EXONE + 15;          /* exponent if normalize shift count were 0 */

+#endif
+

   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
     ecleaz (yi);               /* it was zero */
   else
   if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
     ecleaz (yi);               /* it was zero */
   else


@@ -2910,24 +4081,22 @@ ultoe (lp, y)
 }
 
 
 }
 
 

-/*
-;      Find long integer and fractional parts
-
-;      long i;
-;      unsigned EMUSHORT x[NE], frac[NE];
-;      xifrac (x, &i, frac);
+/* Find signed HOST_WIDE_INT integer I and floating point fractional
+   part FRAC of e-type (packed internal format) floating point input X.
+   The integer output I has the sign of the input, except that
+   positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
+   The output e-type fraction FRAC is the positive fractional
+   part of abs (X).  */

 
 

-  The integer output has the sign of the input.  The fraction is
-the positive fractional part of abs (x).
-*/
-void 
+static void 

 eifrac (x, i, frac)
      unsigned EMUSHORT *x;
 eifrac (x, i, frac)
      unsigned EMUSHORT *x;

-     long *i;
+     HOST_WIDE_INT *i;

      unsigned EMUSHORT *frac;
 {
   unsigned EMUSHORT xi[NI];
      unsigned EMUSHORT *frac;
 {
   unsigned EMUSHORT xi[NI];

-  int k;
+  int j, k;
+  unsigned HOST_WIDE_INT ll;

 
   emovi (x, xi);
   k = (int) xi[E] - (EXONE - 1);
 
   emovi (x, xi);
   k = (int) xi[E] - (EXONE - 1);


@@ -2938,46 +4107,54 @@ eifrac (x, i, frac)
       emovo (xi, frac);
       return;
     }
       emovo (xi, frac);
       return;
     }

-  if (k > (HOST_BITS_PER_LONG - 1))
+  if (k > (HOST_BITS_PER_WIDE_INT - 1))

     {
     {

-      /*
-        ;      long integer overflow: output large integer
-        ;      and correct fraction
-        */
+      /* long integer overflow: output large integer
+        and correct fraction  */

       if (xi[0])
       if (xi[0])

-       *i = ((unsigned long) 1) << (HOST_BITS_PER_LONG - 1);
+       *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);

       else
       else

-       *i = (((unsigned long) 1) << (HOST_BITS_PER_LONG - 1)) - 1;
+       {
+#ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
+         /* In this case, let it overflow and convert as if unsigned.  */
+         euifrac (x, &ll, frac);
+         *i = (HOST_WIDE_INT) ll;
+         return;
+#else
+         /* In other cases, return the largest positive integer.  */
+         *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
+#endif
+       }

       eshift (xi, k);
       if (extra_warnings)
        warning ("overflow on truncation to integer");
       eshift (xi, k);
       if (extra_warnings)
        warning ("overflow on truncation to integer");

-      goto lab11;

     }
     }

-
-  if (k > 16)
+  else if (k > 16)

     {
     {

-      /*
-        ; shift more than 16 bits: shift up k-16, output the integer,
-        ; then complete the shift to get the fraction.
-        */
-      k -= 16;
-      eshift (xi, k);
-
-      *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
-      eshup6 (xi);
-      goto lab10;
+      /* Shift more than 16 bits: first shift up k-16 mod 16,
+        then shift up by 16's.  */
+      j = k - ((k >> 4) << 4);
+      eshift (xi, j);
+      ll = xi[M];
+      k -= j;
+      do
+       {
+         eshup6 (xi);
+         ll = (ll << 16) | xi[M];
+       }
+      while ((k -= 16) > 0);
+      *i = ll;
+      if (xi[0])
+       *i = -(*i);

     }
     }

-
-  /* shift not more than 16 bits */
-  eshift (xi, k);
-  *i = (long) xi[M] & 0xffff;
-
- lab10:
-
-  if (xi[0])
-    *i = -(*i);
- lab11:
-
+  else
+      {
+        /* shift not more than 16 bits */
+          eshift (xi, k);
+        *i = (HOST_WIDE_INT) xi[M] & 0xffff;
+        if (xi[0])
+         *i = -(*i);
+      }

   xi[0] = 0;
   xi[E] = EXONE - 1;
   xi[M] = 0;
   xi[0] = 0;
   xi[E] = EXONE - 1;
   xi[M] = 0;


@@ -2990,24 +4167,19 @@ eifrac (x, i, frac)
 }
 
 
 }
 
 

-/*
-;      Find unsigned long integer and fractional parts
-
-;      unsigned long i;
-;      unsigned EMUSHORT x[NE], frac[NE];
-;      xifrac (x, &i, frac);
+/* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
+   FRAC of e-type X.  A negative input yields integer output = 0 but
+   correct fraction.  */

 
 

-  A negative e type input yields integer output = 0
-  but correct fraction.
-*/
-void 
+static void 

 euifrac (x, i, frac)
      unsigned EMUSHORT *x;
 euifrac (x, i, frac)
      unsigned EMUSHORT *x;

-     long *i;
+     unsigned HOST_WIDE_INT *i;

      unsigned EMUSHORT *frac;
 {
      unsigned EMUSHORT *frac;
 {

+  unsigned HOST_WIDE_INT ll;

   unsigned EMUSHORT xi[NI];
   unsigned EMUSHORT xi[NI];

-  int k;
+  int j, k;

 
   emovi (x, xi);
   k = (int) xi[E] - (EXONE - 1);
 
   emovi (x, xi);
   k = (int) xi[E] - (EXONE - 1);


@@ -3018,40 +4190,41 @@ euifrac (x, i, frac)
       emovo (xi, frac);
       return;
     }
       emovo (xi, frac);
       return;
     }

-  if (k > 32)
+  if (k > HOST_BITS_PER_WIDE_INT)

     {
     {

-      /*
-        ;      long integer overflow: output large integer
-        ;      and correct fraction
-        */
+      /* Long integer overflow: output large integer
+        and correct fraction.
+        Note, the BSD microvax compiler says that ~(0UL)
+        is a syntax error.  */

       *i = ~(0L);
       eshift (xi, k);
       if (extra_warnings)
        warning ("overflow on truncation to unsigned integer");
       *i = ~(0L);
       eshift (xi, k);
       if (extra_warnings)
        warning ("overflow on truncation to unsigned integer");

-      goto lab10;

     }
     }

-
-  if (k > 16)
+  else if (k > 16)

     {
     {

-      /*
-        ; shift more than 16 bits: shift up k-16, output the integer,
-        ; then complete the shift to get the fraction.
-        */
-      k -= 16;
+      /* Shift more than 16 bits: first shift up k-16 mod 16,
+        then shift up by 16's.  */
+      j = k - ((k >> 4) << 4);
+      eshift (xi, j);
+      ll = xi[M];
+      k -= j;
+      do
+       {
+         eshup6 (xi);
+         ll = (ll << 16) | xi[M];
+       }
+      while ((k -= 16) > 0);
+      *i = ll;
+    }
+  else
+    {
+      /* shift not more than 16 bits */

       eshift (xi, k);
       eshift (xi, k);

-
-      *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
-      eshup6 (xi);
-      goto lab10;
+      *i = (HOST_WIDE_INT) xi[M] & 0xffff;

     }
 
     }
 

-  /* shift not more than 16 bits */
-  eshift (xi, k);
-  *i = (long) xi[M] & 0xffff;
-
- lab10:
-
-  if (xi[0])
+  if (xi[0])  /* A negative value yields unsigned integer 0.  */

     *i = 0L;
 
   xi[0] = 0;
     *i = 0L;
 
   xi[0] = 0;


@@ -3065,15 +4238,9 @@ euifrac (x, i, frac)
   emovo (xi, frac);
 }
 
   emovo (xi, frac);
 }
 

+/* Shift the significand of exploded e-type X up or down by SC bits.  */

 
 

-
-/*
-;      Shift significand
-;
-;      Shifts significand area up or down by the number of bits
-;      given by the variable sc.
-*/
-int 
+static int 

 eshift (x, sc)
      unsigned EMUSHORT *x;
      int sc;
 eshift (x, sc)
      unsigned EMUSHORT *x;
      int sc;


@@ -3136,15 +4303,10 @@ eshift (x, sc)
   return ((int) lost);
 }
 
   return ((int) lost);
 }
 

+/* Shift normalize the significand area of exploded e-type X.
+   Return the shift count (up = positive).  */

 
 

-
-/*
-;      normalize
-;
-; Shift normalizes the significand area pointed to by argument
-; shift count (up = positive) is returned.
-*/
-int 
+static int 

 enormlz (x)
      unsigned EMUSHORT x[];
 {
 enormlz (x)
      unsigned EMUSHORT x[];
 {


@@ -3162,9 +4324,9 @@ enormlz (x)
     {
       eshup6 (x);
       sc += 16;
     {
       eshup6 (x);
       sc += 16;

+

       /* With guard word, there are NBITS+16 bits available.
       /* With guard word, there are NBITS+16 bits available.

-       * return true if all are zero.
-       */
+       Return true if all are zero.  */

       if (sc > NBITS)
        return (sc);
     }
       if (sc > NBITS)
        return (sc);
     }


@@ -3210,16 +4372,73 @@ enormlz (x)
   return (sc);
 }
 
   return (sc);
 }
 

-
-
-
-/* Convert e type number to decimal format ASCII string.
- * The constants are for 64 bit precision.
- */
+/* Powers of ten used in decimal <-> binary conversions.  */

 
 #define NTEN 12
 #define MAXP 4096
 
 
 #define NTEN 12
 #define MAXP 4096
 

+#if LONG_DOUBLE_TYPE_SIZE == 128
+static unsigned EMUSHORT etens[NTEN + 1][NE] =
+{
+  {0x6576, 0x4a92, 0x804a, 0x153f,
+   0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},   /* 10**4096 */
+  {0x6a32, 0xce52, 0x329a, 0x28ce,
+   0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},   /* 10**2048 */
+  {0x526c, 0x50ce, 0xf18b, 0x3d28,
+   0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
+  {0x9c66, 0x58f8, 0xbc50, 0x5c54,
+   0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
+  {0x851e, 0xeab7, 0x98fe, 0x901b,
+   0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
+  {0x0235, 0x0137, 0x36b1, 0x336c,
+   0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
+  {0x50f8, 0x25fb, 0xc76b, 0x6b71,
+   0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
+  {0x0000, 0x0000, 0x0000, 0x0000,
+   0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
+  {0x0000, 0x0000, 0x0000, 0x0000,
+   0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
+  {0x0000, 0x0000, 0x0000, 0x0000,
+   0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
+  {0x0000, 0x0000, 0x0000, 0x0000,
+   0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
+  {0x0000, 0x0000, 0x0000, 0x0000,
+   0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
+  {0x0000, 0x0000, 0x0000, 0x0000,
+   0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},   /* 10**1 */
+};
+
+static unsigned EMUSHORT emtens[NTEN + 1][NE] =
+{
+  {0x2030, 0xcffc, 0xa1c3, 0x8123,
+   0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},   /* 10**-4096 */
+  {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
+   0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},   /* 10**-2048 */
+  {0xf53f, 0xf698, 0x6bd3, 0x0158,
+   0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
+  {0xe731, 0x04d4, 0xe3f2, 0xd332,
+   0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
+  {0xa23e, 0x5308, 0xfefb, 0x1155,
+   0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
+  {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
+   0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
+  {0x2a20, 0x6224, 0x47b3, 0x98d7,
+   0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
+  {0x0b5b, 0x4af2, 0xa581, 0x18ed,
+   0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
+  {0xbf71, 0xa9b3, 0x7989, 0xbe68,
+   0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
+  {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
+   0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
+  {0xc155, 0xa4a8, 0x404e, 0x6113,
+   0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
+  {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
+   0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
+  {0xcccd, 0xcccc, 0xcccc, 0xcccc,
+   0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},   /* 10**-1 */
+};
+#else
+/* LONG_DOUBLE_TYPE_SIZE is other than 128 */

 static unsigned EMUSHORT etens[NTEN + 1][NE] =
 {
   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},   /* 10**4096 */
 static unsigned EMUSHORT etens[NTEN + 1][NE] =
 {
   {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,},   /* 10**4096 */


@@ -3253,8 +4472,12 @@ static unsigned EMUSHORT emtens[NTEN + 1][NE] =
   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},   /* 10**-1 */
 };
   {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
   {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},   /* 10**-1 */
 };

+#endif

 
 

-void 
+/* Convert float value X to ASCII string STRING with NDIG digits after
+   the decimal point.  */
+
+static void 

 e24toasc (x, string, ndigs)
      unsigned EMUSHORT x[];
      char *string;
 e24toasc (x, string, ndigs)
      unsigned EMUSHORT x[];
      char *string;


@@ -3262,30 +4485,14 @@ e24toasc (x, string, ndigs)
 {
   unsigned EMUSHORT w[NI];
 
 {
   unsigned EMUSHORT w[NI];
 

-#ifdef INFINITY
-#ifdef IBMPC
-  if ((x[1] & 0x7f80) == 0x7f80)
-#else
-  if ((x[0] & 0x7f80) == 0x7f80)
-#endif
-    {
-#ifdef IBMPC
-      if (x[1] & 0x8000)
-#else
-      if (x[0] & 0x8000)
-#endif
-        sprintf (string, " -Infinity ");
-      else
-        sprintf (string, " Infinity ");
-      return;
-    }
-#endif

   e24toe (x, w);
   etoasc (w, string, ndigs);
 }
 
   e24toe (x, w);
   etoasc (w, string, ndigs);
 }
 

+/* Convert double value X to ASCII string STRING with NDIG digits after
+   the decimal point.  */

 
 

-void 
+static void 

 e53toasc (x, string, ndigs)
      unsigned EMUSHORT x[];
      char *string;
 e53toasc (x, string, ndigs)
      unsigned EMUSHORT x[];
      char *string;


@@ -3293,30 +4500,14 @@ e53toasc (x, string, ndigs)
 {
   unsigned EMUSHORT w[NI];
 
 {
   unsigned EMUSHORT w[NI];
 

-#ifdef INFINITY
-#ifdef IBMPC
-  if ((x[3] & 0x7ff0) == 0x7ff0)
-#else
-  if ((x[0] & 0x7ff0) == 0x7ff0)
-#endif
-    {
-#ifdef IBMPC
-      if (x[3] & 0x8000)
-#else
-      if (x[0] & 0x8000)
-#endif
-        sprintf (string, " -Infinity ");
-      else
-        sprintf (string, " Infinity ");
-      return;
-    }
-#endif

   e53toe (x, w);
   etoasc (w, string, ndigs);
 }
 
   e53toe (x, w);
   etoasc (w, string, ndigs);
 }
 

+/* Convert double extended value X to ASCII string STRING with NDIG digits
+   after the decimal point.  */

 
 

-void 
+static void 

 e64toasc (x, string, ndigs)
      unsigned EMUSHORT x[];
      char *string;
 e64toasc (x, string, ndigs)
      unsigned EMUSHORT x[];
      char *string;


@@ -3324,32 +4515,31 @@ e64toasc (x, string, ndigs)
 {
   unsigned EMUSHORT w[NI];
 
 {
   unsigned EMUSHORT w[NI];
 

-#ifdef INFINITY
-#ifdef IBMPC
-  if ((x[4] & 0x7fff) == 0x7fff)
-#else
-  if ((x[0] & 0x7fff) == 0x7fff)
-#endif
-    {
-#ifdef IBMPC
-      if (x[4] & 0x8000)
-#else
-      if (x[0] & 0x8000)
-#endif
-        sprintf (string, " -Infinity ");
-      else
-        sprintf (string, " Infinity ");
-      return;
-    }
-#endif

   e64toe (x, w);
   etoasc (w, string, ndigs);
 }
 
   e64toe (x, w);
   etoasc (w, string, ndigs);
 }
 

+/* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
+   after the decimal point.  */
+
+static void 
+e113toasc (x, string, ndigs)
+     unsigned EMUSHORT x[];
+     char *string;
+     int ndigs;
+{
+  unsigned EMUSHORT w[NI];
+
+  e113toe (x, w);
+  etoasc (w, string, ndigs);
+}
+
+/* Convert e-type X to ASCII string STRING with NDIGS digits after
+   the decimal point.  */

 
 static char wstring[80];       /* working storage for ASCII output */
 
 
 static char wstring[80];       /* working storage for ASCII output */
 

-void 
+static void 

 etoasc (x, string, ndigs)
      unsigned EMUSHORT x[];
      char *string;
 etoasc (x, string, ndigs)
      unsigned EMUSHORT x[];
      char *string;


@@ -3363,11 +4553,19 @@ etoasc (x, string, ndigs)
   char *s, *ss;
   unsigned EMUSHORT m;
 
   char *s, *ss;
   unsigned EMUSHORT m;
 

+
+  rndsav = rndprc;

   ss = string;
   s = wstring;
   ss = string;
   s = wstring;

-  while ((*s++ = *ss++) != '\0')
-    ;
-  rndsav = rndprc;
+  *ss = '\0';
+  *s = '\0';
+#ifdef NANS
+  if (eisnan (x))
+    {
+      sprintf (wstring, " NaN ");
+      goto bxit;
+    }
+#endif

   rndprc = NBITS;              /* set to full precision */
   emov (x, y);                 /* retain external format */
   if (y[NE - 1] & 0x8000)
   rndprc = NBITS;              /* set to full precision */
   emov (x, y);                 /* retain external format */
   if (y[NE - 1] & 0x8000)


@@ -3390,12 +4588,11 @@ etoasc (x, string, ndigs)
          if (y[k] != 0)
            goto tnzro;         /* denormalized number */
        }
          if (y[k] != 0)
            goto tnzro;         /* denormalized number */
        }

-      goto isone;              /* legal all zeros */
+      goto isone;              /* valid all zeros */

     }
  tnzro:
 
     }
  tnzro:
 

-  /* Test for infinity.  Don't bother with illegal infinities.
-   */
+  /* Test for infinity.  */

   if (y[NE - 1] == 0x7fff)
     {
       if (sign)
   if (y[NE - 1] == 0x7fff)
     {
       if (sign)


@@ -3420,9 +4617,12 @@ etoasc (x, string, ndigs)
   if (i == 0)
     goto isone;
 
   if (i == 0)
     goto isone;
 

+  if (i == -2)
+    abort ();
+

   if (i < 0)
     {                          /* Number is greater than 1 */
   if (i < 0)
     {                          /* Number is greater than 1 */

-      /* Convert significand to an integer and strip trailing decimal zeros. */
+      /* Convert significand to an integer and strip trailing decimal zeros.  */

       emov (y, u);
       u[NE - 1] = EXONE + NBITS - 1;
 
       emov (y, u);
       u[NE - 1] = EXONE + NBITS - 1;
 


@@ -3452,6 +4652,7 @@ etoasc (x, string, ndigs)
       emov (eone, t);
       m = MAXP;
       p = &etens[0][0];
       emov (eone, t);
       m = MAXP;
       p = &etens[0][0];

+      /* An unordered compare result shouldn't happen here.  */

       while (ecmp (ten, u) <= 0)
        {
          if (ecmp (p, u) <= 0)
       while (ecmp (ten, u) <= 0)
        {
          if (ecmp (p, u) <= 0)


@@ -3468,7 +4669,7 @@ etoasc (x, string, ndigs)
     }
   else
     {                          /* Number is less than 1.0 */
     }
   else
     {                          /* Number is less than 1.0 */

-      /* Pad significand with trailing decimal zeros. */
+      /* Pad significand with trailing decimal zeros.  */

       if (y[NE - 1] == 0)
        {
          while ((y[NE - 2] & 0x8000) == 0)
       if (y[NE - 1] == 0)
        {
          while ((y[NE - 2] & 0x8000) == 0)


@@ -3526,7 +4727,7 @@ etoasc (x, string, ndigs)
       ediv (t, eone, t);
     }
  isone:
       ediv (t, eone, t);
     }
  isone:

-  /* Find the first (leading) digit. */
+  /* Find the first (leading) digit.  */

   emovi (t, w);
   emovz (w, t);
   emovi (y, w);
   emovi (t, w);
   emovz (w, t);
   emovi (y, w);


@@ -3549,7 +4750,7 @@ etoasc (x, string, ndigs)
     *s++ = '-';
   else
     *s++ = ' ';
     *s++ = '-';
   else
     *s++ = ' ';

-  /* Examine number of digits requested by caller. */
+  /* Examine number of digits requested by caller.  */

   if (ndigs < 0)
     ndigs = 0;
   if (ndigs > NDEC)
   if (ndigs < 0)
     ndigs = 0;
   if (ndigs > NDEC)


@@ -3567,10 +4768,10 @@ etoasc (x, string, ndigs)
     }
   else
     {
     }
   else
     {

-      *s++ = (char )digit + '0';
+      *s++ = (char)digit + '0';

       *s++ = '.';
     }
       *s++ = '.';
     }

-  /* Generate digits after the decimal point. */
+  /* Generate digits after the decimal point.  */

   for (k = 0; k <= ndigs; k++)
     {
       /* multiply current number by 10, without normalizing */
   for (k = 0; k <= ndigs; k++)
     {
       /* multiply current number by 10, without normalizing */


@@ -3588,7 +4789,7 @@ etoasc (x, string, ndigs)
   /* round off the ASCII string */
   if (digit > 4)
     {
   /* round off the ASCII string */
   if (digit > 4)
     {

-      /* Test for critical rounding case in ASCII output. */
+      /* Test for critical rounding case in ASCII output.  */

       if (digit == 5)
        {
          emovo (y, t);
       if (digit == 5)
        {
          emovo (y, t);


@@ -3645,29 +4846,16 @@ etoasc (x, string, ndigs)
 }
 
 
 }
 
 

+/* Convert ASCII string to floating point.

 
 

+   Numeric input is a free format decimal number of any length, with
+   or without decimal point.  Entering E after the number followed by an
+   integer number causes the second number to be interpreted as a power of
+   10 to be multiplied by the first number (i.e., "scientific" notation).  */

 
 

-/*
-;                                                              ASCTOQ
-;              ASCTOQ.MAC              LATEST REV: 11 JAN 84
-;                                      SLM, 3 JAN 78
-;
-;      Convert ASCII string to quadruple precision floating point
-;
-;              Numeric input is free field decimal number
-;              with max of 15 digits with or without
-;              decimal point entered as ASCII from teletype.
-;      Entering E after the number followed by a second
-;      number causes the second number to be interpreted
-;      as a power of 10 to be multiplied by the first number
-;      (i.e., "scientific" notation).
-;
-;      Usage:
-;              asctoq (string, q);
-*/
-
-/* ASCII to single */
-void 
+/* Convert ASCII string S to single precision float value Y.  */
+
+static void 

 asctoe24 (s, y)
      char *s;
      unsigned EMUSHORT *y;
 asctoe24 (s, y)
      char *s;
      unsigned EMUSHORT *y;


@@ -3676,13 +4864,14 @@ asctoe24 (s, y)
 }
 
 
 }
 
 

-/* ASCII to double */
-void 
+/* Convert ASCII string S to double precision value Y.  */
+
+static void 

 asctoe53 (s, y)
      char *s;
      unsigned EMUSHORT *y;
 {
 asctoe53 (s, y)
      char *s;
      unsigned EMUSHORT *y;
 {

-#ifdef DEC
+#if defined(DEC) || defined(IBM)

   asctoeg (s, y, 56);
 #else
   asctoeg (s, y, 53);
   asctoeg (s, y, 56);
 #else
   asctoeg (s, y, 53);


@@ -3690,8 +4879,9 @@ asctoe53 (s, y)
 }
 
 
 }
 
 

-/* ASCII to long double */
-void 
+/* Convert ASCII string S to double extended value Y.  */
+
+static void 

 asctoe64 (s, y)
      char *s;
      unsigned EMUSHORT *y;
 asctoe64 (s, y)
      char *s;
      unsigned EMUSHORT *y;


@@ -3699,8 +4889,19 @@ asctoe64 (s, y)
   asctoeg (s, y, 64);
 }
 
   asctoeg (s, y, 64);
 }
 

-/* ASCII to super double */
-void 
+/* Convert ASCII string S to 128-bit long double Y.  */
+
+static void 
+asctoe113 (s, y)
+     char *s;
+     unsigned EMUSHORT *y;
+{
+  asctoeg (s, y, 113);
+}
+
+/* Convert ASCII string S to e type Y.  */
+
+static void 

 asctoe (s, y)
      char *s;
      unsigned EMUSHORT *y;
 asctoe (s, y)
      char *s;
      unsigned EMUSHORT *y;


@@ -3708,10 +4909,10 @@ asctoe (s, y)
   asctoeg (s, y, NBITS);
 }
 
   asctoeg (s, y, NBITS);
 }
 

-/* Space to make a copy of the input string: */
-static char lstr[82];
+/* Convert ASCII string SS to e type Y, with a specified rounding precision
+   of OPREC bits.  */

 
 

-void 
+static void 

 asctoeg (ss, y, oprec)
      char *ss;
      unsigned EMUSHORT *y;
 asctoeg (ss, y, oprec)
      char *ss;
      unsigned EMUSHORT *y;


@@ -3722,19 +4923,16 @@ asctoeg (ss, y, oprec)
   int k, trail, c, rndsav;
   EMULONG lexp;
   unsigned EMUSHORT nsign, *p;
   int k, trail, c, rndsav;
   EMULONG lexp;
   unsigned EMUSHORT nsign, *p;

-  char *sp, *s;
+  char *sp, *s, *lstr;

 
 

-  /* Copy the input string. */
+  /* Copy the input string.  */
+  lstr = (char *) alloca (strlen (ss) + 1);

   s = ss;
   while (*s == ' ')            /* skip leading spaces */
     ++s;
   sp = lstr;
   s = ss;
   while (*s == ' ')            /* skip leading spaces */
     ++s;
   sp = lstr;

-  for (k = 0; k < 79; k++)
-    {
-      if ((*sp++ = *s++) == '\0')
-       break;
-    }
-  *sp = '\0';
+  while ((*sp++ = *s++) != '\0')
+    ;

   s = lstr;
 
   rndsav = rndprc;
   s = lstr;
 
   rndsav = rndprc;


@@ -3756,7 +4954,7 @@ asctoeg (ss, y, oprec)
       /* Ignore leading zeros */
       if ((prec == 0) && (decflg == 0) && (k == 0))
        goto donchr;
       /* Ignore leading zeros */
       if ((prec == 0) && (decflg == 0) && (k == 0))
        goto donchr;

-      /* Identify and strip trailing zeros after the decimal point. */
+      /* Identify and strip trailing zeros after the decimal point.  */

       if ((trail == 0) && (decflg != 0))
        {
          sp = s;
       if ((trail == 0) && (decflg != 0))
        {
          sp = s;


@@ -3775,11 +4973,12 @@ asctoeg (ss, y, oprec)
          if (*s == 'z')
            goto donchr;
        }
          if (*s == 'z')
            goto donchr;
        }

+

       /* If enough digits were given to more than fill up the yy register,
       /* If enough digits were given to more than fill up the yy register,

-       * continuing until overflow into the high guard word yy[2]
-       * guarantees that there will be a roundoff bit at the top
-       * of the low guard word after normalization.
-       */
+        continuing until overflow into the high guard word yy[2]
+        guarantees that there will be a roundoff bit at the top
+        of the low guard word after normalization.  */
+

       if (yy[2] == 0)
        {
          if (decflg)
       if (yy[2] == 0)
        {
          if (decflg)


@@ -3795,7 +4994,11 @@ asctoeg (ss, y, oprec)
        }
       else
        {
        }
       else
        {

+         /* Mark any lost non-zero digit.  */

          lost |= k;
          lost |= k;

+         /* Count lost digits before the decimal point.  */
+         if (decflg == 0)
+           nexp -= 1;

        }
       prec += 1;
       goto donchr;
        }
       prec += 1;
       goto donchr;


@@ -3835,8 +5038,12 @@ asctoeg (ss, y, oprec)
       goto infinite;
     default:
     error:
       goto infinite;
     default:
     error:

+#ifdef NANS
+      einan (yy);
+#else

       mtherr ("asctoe", DOMAIN);
       mtherr ("asctoe", DOMAIN);

-      eclear (y);
+      eclear (yy);
+#endif

       goto aexit;
     }
  donchr:
       goto aexit;
     }
  donchr:


@@ -3845,7 +5052,15 @@ asctoeg (ss, y, oprec)
 
   /* Exponent interpretation */
  expnt:
 
   /* Exponent interpretation */
  expnt:

+  /* 0.0eXXX is zero, regardless of XXX.  Check for the 0.0. */
+  for (k = 0; k < NI; k++)
+    {
+      if (yy[k] != 0)
+       goto read_expnt;
+    }
+  goto aexit;

 
 

+read_expnt:

   esign = 1;
   exp = 0;
   ++s;
   esign = 1;
   exp = 0;
   ++s;


@@ -3861,7 +5076,7 @@ asctoeg (ss, y, oprec)
     {
       exp *= 10;
       exp += *s++ - '0';
     {
       exp *= 10;
       exp += *s++ - '0';

-      if (exp > 4956)
+      if (exp > -(MINDECEXP))

        {
          if (esign < 0)
            goto zero;
        {
          if (esign < 0)
            goto zero;


@@ -3871,14 +5086,14 @@ asctoeg (ss, y, oprec)
     }
   if (esign < 0)
     exp = -exp;
     }
   if (esign < 0)
     exp = -exp;

-  if (exp > 4932)
+  if (exp > MAXDECEXP)

     {
  infinite:
       ecleaz (yy);
       yy[E] = 0x7fff;          /* infinity */
       goto aexit;
     }
     {
  infinite:
       ecleaz (yy);
       yy[E] = 0x7fff;          /* infinity */
       goto aexit;
     }

-  if (exp < -4956)
+  if (exp < MINDECEXP)

     {
  zero:
       ecleaz (yy);
     {
  zero:
       ecleaz (yy);


@@ -3887,7 +5102,7 @@ asctoeg (ss, y, oprec)
 
  daldone:
   nexp = exp - nexp;
 
  daldone:
   nexp = exp - nexp;

-  /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
+  /* Pad trailing zeros to minimize power of 10, per IEEE spec.  */

   while ((nexp > 0) && (yy[2] == 0))
     {
       emovz (yy, xt);
   while ((nexp > 0) && (yy[2] == 0))
     {
       emovz (yy, xt);


@@ -3907,15 +5122,15 @@ asctoeg (ss, y, oprec)
     }
   lexp = (EXONE - 1 + NBITS) - k;
   emdnorm (yy, lost, 0, lexp, 64);
     }
   lexp = (EXONE - 1 + NBITS) - k;
   emdnorm (yy, lost, 0, lexp, 64);

-  /* convert to external format */

 
 

+  /* Convert to external format:
+
+     Multiply by 10**nexp.  If precision is 64 bits,
+     the maximum relative error incurred in forming 10**n
+     for 0 <= n <= 324 is 8.2e-20, at 10**180.
+     For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
+     For 0 >= n >= -999, it is -1.55e-19 at 10**-435.  */

 
 

-  /* Multiply by 10**nexp.  If precision is 64 bits,
-   * the maximum relative error incurred in forming 10**n
-   * for 0 <= n <= 324 is 8.2e-20, at 10**180.
-   * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
-   * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
-   */

   lexp = yy[E];
   if (nexp == 0)
     {
   lexp = yy[E];
   if (nexp == 0)
     {


@@ -3928,7 +5143,8 @@ asctoeg (ss, y, oprec)
       nexp = -nexp;
       esign = -1;
       if (nexp > 4096)
       nexp = -nexp;
       esign = -1;
       if (nexp > 4096)

-       {                       /* Punt.  Can't handle this without 2 divides. */
+       {
+         /* Punt.  Can't handle this without 2 divides.  */

          emovi (etens[0], tt);
          lexp -= tt[E];
          k = edivm (tt, yy);
          emovi (etens[0], tt);
          lexp -= tt[E];
          k = edivm (tt, yy);


@@ -3967,8 +5183,13 @@ asctoeg (ss, y, oprec)
   /* Round and convert directly to the destination type */
   if (oprec == 53)
     lexp -= EXONE - 0x3ff;
   /* Round and convert directly to the destination type */
   if (oprec == 53)
     lexp -= EXONE - 0x3ff;

+#ifdef IBM
+  else if (oprec == 24 || oprec == 56)
+    lexp -= EXONE - (0x41 << 2);
+#else

   else if (oprec == 24)
     lexp -= EXONE - 0177;
   else if (oprec == 24)
     lexp -= EXONE - 0177;

+#endif

 #ifdef DEC
   else if (oprec == 56)
     lexp -= EXONE - 0201;
 #ifdef DEC
   else if (oprec == 56)
     lexp -= EXONE - 0201;


@@ -3987,6 +5208,11 @@ asctoeg (ss, y, oprec)
       todec (yy, y);           /* see etodec.c */
       break;
 #endif
       todec (yy, y);           /* see etodec.c */
       break;
 #endif

+#ifdef IBM
+    case 56:
+      toibm (yy, y, DFmode);
+      break;
+#endif

     case 53:
       toe53 (yy, y);
       break;
     case 53:
       toe53 (yy, y);
       break;


@@ -3996,6 +5222,9 @@ asctoeg (ss, y, oprec)
     case 64:
       toe64 (yy, y);
       break;
     case 64:
       toe64 (yy, y);
       break;

+    case 113:
+      toe113 (yy, y);
+      break;

     case NBITS:
       emovo (yy, y);
       break;
     case NBITS:
       emovo (yy, y);
       break;


@@ -4004,13 +5233,9 @@ asctoeg (ss, y, oprec)
 
 
 
 
 
 

-/* y = largest integer not greater than x
- * (truncated toward minus infinity)
- *
- * unsigned EMUSHORT x[NE], y[NE]
- *
- * efloor (x, y);
- */
+/* Return Y = largest integer not greater than X (truncated toward minus
+   infinity).  */
+

 static unsigned EMUSHORT bmask[] =
 {
   0xffff,
 static unsigned EMUSHORT bmask[] =
 {
   0xffff,


@@ -4032,7 +5257,7 @@ static unsigned EMUSHORT bmask[] =
   0x0000,
 };
 
   0x0000,
 };
 

-void 
+static void 

 efloor (x, y)
      unsigned EMUSHORT x[], y[];
 {
 efloor (x, y)
      unsigned EMUSHORT x[], y[];
 {


@@ -4079,16 +5304,10 @@ efloor (x, y)
 }
 
 
 }
 
 

-/* unsigned EMUSHORT x[], s[];
- * int *exp;
- *
- * efrexp (x, exp, s);
- *
- * Returns s and exp such that  s * 2**exp = x and .5 <= s < 1.
- * For example, 1.1 = 0.55 * 2**1
- * Handles denormalized numbers properly using long integer exp.
- */
-void 
+/* Return S and EXP such that  S * 2^EXP = X and .5 <= S < 1.
+   For example, 1.1 = 0.55 * 2^1.  */
+
+static void 

 efrexp (x, exp, s)
      unsigned EMUSHORT x[];
      int *exp;
 efrexp (x, exp, s)
      unsigned EMUSHORT x[];
      int *exp;


@@ -4098,6 +5317,7 @@ efrexp (x, exp, s)
   EMULONG li;
 
   emovi (x, xi);
   EMULONG li;
 
   emovi (x, xi);

+  /*  Handle denormalized numbers properly using long integer exponent.  */

   li = (EMULONG) ((EMUSHORT) xi[1]);
 
   if (li == 0)
   li = (EMULONG) ((EMUSHORT) xi[1]);
 
   if (li == 0)


@@ -4109,16 +5329,9 @@ efrexp (x, exp, s)
   *exp = (int) (li - 0x3ffe);
 }
 
   *exp = (int) (li - 0x3ffe);
 }
 

+/* Return e type Y = X * 2^PWR2.  */

 
 

-
-/* unsigned EMUSHORT x[], y[];
- * long pwr2;
- *
- * eldexp (x, pwr2, y);
- *
- * Returns y = x * 2**pwr2.
- */
-void 
+static void 

 eldexp (x, pwr2, y)
      unsigned EMUSHORT x[];
      int pwr2;
 eldexp (x, pwr2, y)
      unsigned EMUSHORT x[];
      int pwr2;


@@ -4137,15 +5350,25 @@ eldexp (x, pwr2, y)
 }
 
 
 }
 
 

-/* c = remainder after dividing b by a
- * Least significant integer quotient bits left in equot[].
- */
-void 
+/* C = remainder after dividing B by A, all e type values.
+   Least significant integer quotient bits left in EQUOT.  */
+
+static void 

 eremain (a, b, c)
      unsigned EMUSHORT a[], b[], c[];
 {
   unsigned EMUSHORT den[NI], num[NI];
 
 eremain (a, b, c)
      unsigned EMUSHORT a[], b[], c[];
 {
   unsigned EMUSHORT den[NI], num[NI];
 

+#ifdef NANS
+  if (eisinf (b)
+      || (ecmp (a, ezero) == 0)
+      || eisnan (a)
+      || eisnan (b))
+    {
+      enan (c, 0);
+      return;
+    }
+#endif

   if (ecmp (a, ezero) == 0)
     {
       mtherr ("eremain", SING);
   if (ecmp (a, ezero) == 0)
     {
       mtherr ("eremain", SING);


@@ -4163,7 +5386,10 @@ eremain (a, b, c)
   emovo (num, c);
 }
 
   emovo (num, c);
 }
 

-void 
+/*  Return quotient of exploded e-types NUM / DEN in EQUOT,
+    remainder in NUM.  */
+
+static void 

 eiremain (den, num)
      unsigned EMUSHORT den[], num[];
 {
 eiremain (den, num)
      unsigned EMUSHORT den[], num[];
 {


@@ -4183,9 +5409,7 @@ eiremain (den, num)
          j = 1;
        }
       else
          j = 1;
        }
       else

-       {

          j = 0;
          j = 0;

-       }

       eshup1 (equot);
       equot[NI - 1] |= j;
       eshup1 (num);
       eshup1 (equot);
       equot[NI - 1] |= j;
       eshup1 (num);


@@ -4194,69 +5418,26 @@ eiremain (den, num)
   emdnorm (num, 0, 0, ln, 0);
 }
 
   emdnorm (num, 0, 0, ln, 0);
 }
 

-/*                                                     mtherr.c
- *
- *     Library common error handling routine
- *
- *
- *
- * SYNOPSIS:
- *
- * char *fctnam;
- * int code;
- * void mtherr ();
- *
- * mtherr (fctnam, code);
- *
- *
- *
- * DESCRIPTION:
- *
- * This routine may be called to report one of the following
- * error conditions (in the include file mconf.h).
- *
- *   Mnemonic        Value          Significance
- *
- *    DOMAIN            1       argument domain error
- *    SING              2       function singularity
- *    OVERFLOW          3       overflow range error
- *    UNDERFLOW         4       underflow range error
- *    TLOSS             5       total loss of precision
- *    PLOSS             6       partial loss of precision
- *    EDOM             33       Unix domain error code
- *    ERANGE           34       Unix range error code
- *
- * The default version of the file prints the function name,
- * passed to it by the pointer fctnam, followed by the
- * error condition.  The display is directed to the standard
- * output device.  The routine then returns to the calling
- * program.  Users may wish to modify the program to abort by
- * calling exit under severe error conditions such as domain
- * errors.
- *
- * Since all error conditions pass control to this function,
- * the display may be easily changed, eliminated, or directed
- * to an error logging device.
- *
- * SEE ALSO:
- *
- * mconf.h
- *
- */
-\f
-/*
-Cephes Math Library Release 2.0:  April, 1987
-Copyright 1984, 1987 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-/* include "mconf.h" */
-
-/* Notice: the order of appearance of the following
- * messages is bound to the error codes defined
- * in mconf.h.
- */
-static char *ermsg[7] =
+/* Report an error condition CODE encountered in function NAME.
+   CODE is one of the following:
+
+    Mnemonic        Value          Significance
+ 
+     DOMAIN            1       argument domain error
+     SING              2       function singularity
+     OVERFLOW          3       overflow range error
+     UNDERFLOW         4       underflow range error
+     TLOSS             5       total loss of precision
+     PLOSS             6       partial loss of precision
+     INVALID           7       NaN - producing operation
+     EDOM             33       Unix domain error code
+     ERANGE           34       Unix range error code
+ 
+   The order of appearance of the following messages is bound to the
+   error codes defined above.  */
+
+#define NMSGS 8
+static char *ermsg[NMSGS] =

 {
   "unknown",                   /* error code 0 */
   "domain",                    /* error code 1 */
 {
   "unknown",                   /* error code 0 */
   "domain",                    /* error code 1 */


@@ -4264,51 +5445,37 @@ static char *ermsg[7] =
   "overflow",
   "underflow",
   "total loss of precision",
   "overflow",
   "underflow",
   "total loss of precision",

-  "partial loss of precision"
+  "partial loss of precision",
+  "invalid operation"

 };
 
 int merror = 0;
 extern int merror;
 
 };
 
 int merror = 0;
 extern int merror;
 

-void 
+static void 

 mtherr (name, code)
      char *name;
      int code;
 {
   char errstr[80];
 
 mtherr (name, code)
      char *name;
      int code;
 {
   char errstr[80];
 

-  /* Display string passed by calling program,
-   * which is supposed to be the name of the
-   * function in which the error occurred.
-   */
+  /* The string passed by the calling program is supposed to be the
+     name of the function in which the error occurred.
+     The code argument selects which error message string will be printed.  */

 
 

-  /* Display error message defined
-   * by the code argument.
-   */
-  if ((code <= 0) || (code >= 6))
+  if ((code <= 0) || (code >= NMSGS))

     code = 0;
     code = 0;

-  sprintf (errstr, "\n%s %s error\n", name, ermsg[code]);
+  sprintf (errstr, " %s %s error", name, ermsg[code]);

   if (extra_warnings)
     warning (errstr);
   /* Set global error message word */
   merror = code + 1;
   if (extra_warnings)
     warning (errstr);
   /* Set global error message word */
   merror = code + 1;

-
-  /* Return to calling
-   * program
-   */

 }
 
 }
 

-/* Here is etodec.c .
- *
- */
+#ifdef DEC
+/* Convert DEC double precision D to e type E.  */

 
 

-/*
-;      convert DEC double precision to e type
-;      double d;
-;      EMUSHORT e[NE];
-;      dectoe (&d, e);
-*/
-void 
+static void 

 dectoe (d, e)
      unsigned EMUSHORT *d;
      unsigned EMUSHORT *e;
 dectoe (d, e)
      unsigned EMUSHORT *d;
      unsigned EMUSHORT *e;


@@ -4346,88 +5513,9 @@ dectoe (d, e)
   emovo (y, e);
 }
 
   emovo (y, e);
 }
 

+/* Convert e type X to DEC double precision D.  */

 
 

-
-/*
-;      convert e type to DEC double precision
-;      double d;
-;      EMUSHORT e[NE];
-;      etodec (e, &d);
-*/
-#if 0
-static unsigned EMUSHORT decbit[NI] = {0, 0, 0, 0, 0, 0, 0200, 0};
-
-void 
-etodec (x, d)
-     unsigned EMUSHORT *x, *d;
-{
-  unsigned EMUSHORT xi[NI];
-  register unsigned EMUSHORT r;
-  int i, j;
-
-  emovi (x, xi);
-  *d = 0;
-  if (xi[0] != 0)
-    *d = 0100000;
-  r = xi[E];
-  if (r < (EXONE - 128))
-    goto zout;
-  i = xi[M + 4];
-  if ((i & 0200) != 0)
-    {
-      if ((i & 0377) == 0200)
-       {
-         if ((i & 0400) != 0)
-           {
-             /* check all less significant bits */
-             for (j = M + 5; j < NI; j++)
-               {
-                 if (xi[j] != 0)
-                   goto yesrnd;
-               }
-           }
-         goto nornd;
-       }
-    yesrnd:
-      eaddm (decbit, xi);
-      r -= enormlz (xi);
-    }
-
- nornd:
-
-  r -= EXONE;
-  r += 0201;
-  if (r < 0)
-    {
-    zout:
-      *d++ = 0;
-      *d++ = 0;
-      *d++ = 0;
-      *d++ = 0;
-      return;
-    }
-  if (r >= 0377)
-    {
-      *d++ = 077777;
-      *d++ = -1;
-      *d++ = -1;
-      *d++ = -1;
-      return;
-    }
-  r &= 0377;
-  r <<= 7;
-  eshup8 (xi);
-  xi[M] &= 0177;
-  r |= xi[M];
-  *d++ |= r;
-  *d++ = xi[M + 1];
-  *d++ = xi[M + 2];
-  *d++ = xi[M + 3];
-}
-
-#else
-
-void 
+static void 

 etodec (x, d)
      unsigned EMUSHORT *x, *d;
 {
 etodec (x, d)
      unsigned EMUSHORT *x, *d;
 {


@@ -4436,8 +5524,9 @@ etodec (x, d)
   int rndsav;
 
   emovi (x, xi);
   int rndsav;
 
   emovi (x, xi);

-  exp = (EMULONG) xi[E] - (EXONE - 0201);      /* adjust exponent for offsets */
-/* round off to nearest or even */
+  /* Adjust exponent for offsets.  */
+  exp = (EMULONG) xi[E] - (EXONE - 0201);
+  /* Round off to nearest or even.  */

   rndsav = rndprc;
   rndprc = 56;
   emdnorm (xi, 0, 0, exp, 64);
   rndsav = rndprc;
   rndprc = 56;
   emdnorm (xi, 0, 0, exp, 64);


@@ -4445,7 +5534,10 @@ etodec (x, d)
   todec (xi, d);
 }
 
   todec (xi, d);
 }
 

-void 
+/* Convert exploded e-type X, that has already been rounded to
+   56-bit precision, to DEC format double Y.  */
+
+static void 

 todec (x, y)
      unsigned EMUSHORT *x, *y;
 {
 todec (x, y)
      unsigned EMUSHORT *x, *y;
 {


@@ -4486,7 +5578,723 @@ todec (x, y)
   *y++ = x[M + 2];
   *y++ = x[M + 3];
 }
   *y++ = x[M + 2];
   *y++ = x[M + 3];
 }

+#endif /* DEC */
+
+#ifdef IBM
+/* Convert IBM single/double precision to e type.  */
+
+static void 
+ibmtoe (d, e, mode)
+     unsigned EMUSHORT *d;
+     unsigned EMUSHORT *e;
+     enum machine_mode mode;
+{
+  unsigned EMUSHORT y[NI];
+  register unsigned EMUSHORT r, *p;
+  int rndsav;

 
 

-#endif /* not 0 */
+  ecleaz (y);                  /* start with a zero */
+  p = y;                       /* point to our number */
+  r = *d;                      /* get IBM exponent word */
+  if (*d & (unsigned int) 0x8000)
+    *p = 0xffff;               /* fill in our sign */
+  ++p;                         /* bump pointer to our exponent word */
+  r &= 0x7f00;                 /* strip the sign bit */
+  r >>= 6;                     /* shift exponent word down 6 bits */
+                               /* in fact shift by 8 right and 2 left */
+  r += EXONE - (0x41 << 2);    /* subtract IBM exponent offset */
+                               /* add our e type exponent offset */
+  *p++ = r;                    /* to form our exponent */
+
+  *p++ = *d++ & 0xff;          /* now do the high order mantissa */
+                               /* strip off the IBM exponent and sign bits */
+  if (mode != SFmode)          /* there are only 2 words in SFmode */
+    {
+      *p++ = *d++;             /* fill in the rest of our mantissa */
+      *p++ = *d++;
+    }
+  *p = *d;
+
+  if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
+    y[0] = y[E] = 0;
+  else
+    y[E] -= 5 + enormlz (y);   /* now normalise the mantissa */
+                             /* handle change in RADIX */
+  emovo (y, e);
+}
+
+
+
+/* Convert e type to IBM single/double precision.  */
+
+static void 
+etoibm (x, d, mode)
+     unsigned EMUSHORT *x, *d;
+     enum machine_mode mode;
+{
+  unsigned EMUSHORT xi[NI];
+  EMULONG exp;
+  int rndsav;
+
+  emovi (x, xi);
+  exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2));       /* adjust exponent for offsets */
+                                                       /* round off to nearest or even */
+  rndsav = rndprc;
+  rndprc = 56;
+  emdnorm (xi, 0, 0, exp, 64);
+  rndprc = rndsav;
+  toibm (xi, d, mode);
+}
+
+static void 
+toibm (x, y, mode)
+     unsigned EMUSHORT *x, *y;
+     enum machine_mode mode;
+{
+  unsigned EMUSHORT i;
+  unsigned EMUSHORT *p;
+  int r;
+
+  p = x;
+  *y = 0;
+  if (*p++)
+    *y = 0x8000;
+  i = *p++;
+  if (i == 0)
+    {
+      *y++ = 0;
+      *y++ = 0;
+      if (mode != SFmode)
+       {
+         *y++ = 0;
+         *y++ = 0;
+       }
+      return;
+    }
+  r = i & 0x3;
+  i >>= 2;
+  if (i > 0x7f)
+    {
+      *y++ |= 0x7fff;
+      *y++ = 0xffff;
+      if (mode != SFmode)
+       {
+         *y++ = 0xffff;
+         *y++ = 0xffff;
+       }
+#ifdef ERANGE
+      errno = ERANGE;
+#endif
+      return;
+    }
+  i &= 0x7f;
+  *y |= (i << 8);
+  eshift (x, r + 5);
+  *y++ |= x[M];
+  *y++ = x[M + 1];
+  if (mode != SFmode)
+    {
+      *y++ = x[M + 2];
+      *y++ = x[M + 3];
+    }
+}
+#endif /* IBM */
+
+/* Output a binary NaN bit pattern in the target machine's format.  */
+
+/* If special NaN bit patterns are required, define them in tm.h
+   as arrays of unsigned 16-bit shorts.  Otherwise, use the default
+   patterns here.  */
+#ifdef TFMODE_NAN
+TFMODE_NAN;
+#else
+#ifdef IEEE
+unsigned EMUSHORT TFbignan[8] =
+ {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
+unsigned EMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
+#endif
+#endif
+
+#ifdef XFMODE_NAN
+XFMODE_NAN;
+#else
+#ifdef IEEE
+unsigned EMUSHORT XFbignan[6] =
+ {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
+unsigned EMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
+#endif
+#endif
+
+#ifdef DFMODE_NAN
+DFMODE_NAN;
+#else
+#ifdef IEEE
+unsigned EMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
+unsigned EMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
+#endif
+#endif
+
+#ifdef SFMODE_NAN
+SFMODE_NAN;
+#else
+#ifdef IEEE
+unsigned EMUSHORT SFbignan[2] = {0x7fff, 0xffff};
+unsigned EMUSHORT SFlittlenan[2] = {0, 0xffc0};
+#endif
+#endif
+
+
+static void
+make_nan (nan, sign, mode)
+     unsigned EMUSHORT *nan;
+     int sign;
+     enum machine_mode mode;
+{
+  int n;
+  unsigned EMUSHORT *p;
+
+  switch (mode)
+    {
+/* Possibly the `reserved operand' patterns on a VAX can be
+   used like NaN's, but probably not in the same way as IEEE.  */
+#if !defined(DEC) && !defined(IBM)
+    case TFmode:
+      n = 8;
+      if (REAL_WORDS_BIG_ENDIAN)
+       p = TFbignan;
+      else
+       p = TFlittlenan;
+      break;
+    case XFmode:
+      n = 6;
+      if (REAL_WORDS_BIG_ENDIAN)
+       p = XFbignan;
+      else
+       p = XFlittlenan;
+      break;
+    case DFmode:
+      n = 4;
+      if (REAL_WORDS_BIG_ENDIAN)
+       p = DFbignan;
+      else
+       p = DFlittlenan;
+      break;
+    case HFmode:
+    case SFmode:
+      n = 2;
+      if (REAL_WORDS_BIG_ENDIAN)
+       p = SFbignan;
+      else
+       p = SFlittlenan;
+      break;
+#endif
+    default:
+      abort ();
+    }
+  if (REAL_WORDS_BIG_ENDIAN)
+    *nan++ = (sign << 15) | *p++;
+  while (--n != 0)
+    *nan++ = *p++;
+  if (! REAL_WORDS_BIG_ENDIAN)
+    *nan = (sign << 15) | *p;
+}
+
+/* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
+   This is the inverse of the function `etarsingle' invoked by
+   REAL_VALUE_TO_TARGET_SINGLE.  */
+
+REAL_VALUE_TYPE
+ereal_from_float (f)
+     HOST_WIDE_INT f;
+{
+  REAL_VALUE_TYPE r;
+  unsigned EMUSHORT s[2];
+  unsigned EMUSHORT e[NE];

 
 

+  /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
+   This is the inverse operation to what the function `endian' does.  */
+  if (REAL_WORDS_BIG_ENDIAN)
+    {
+      s[0] = (unsigned EMUSHORT) (f >> 16);
+      s[1] = (unsigned EMUSHORT) f;
+    }
+  else
+    {
+      s[0] = (unsigned EMUSHORT) f;
+      s[1] = (unsigned EMUSHORT) (f >> 16);
+    }
+  /* Convert and promote the target float to E-type.  */
+  e24toe (s, e);
+  /* Output E-type to REAL_VALUE_TYPE.  */
+  PUT_REAL (e, &r);
+  return r;
+}
+
+
+/* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
+   This is the inverse of the function `etardouble' invoked by
+   REAL_VALUE_TO_TARGET_DOUBLE.
+
+   The DFmode is stored as an array of HOST_WIDE_INT in the target's
+   data format, with no holes in the bit packing.  The first element
+   of the input array holds the bits that would come first in the
+   target computer's memory.  */
+
+REAL_VALUE_TYPE
+ereal_from_double (d)
+     HOST_WIDE_INT d[];
+{
+  REAL_VALUE_TYPE r;
+  unsigned EMUSHORT s[4];
+  unsigned EMUSHORT e[NE];
+
+  /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces.  */
+  if (REAL_WORDS_BIG_ENDIAN)
+    {
+      s[0] = (unsigned EMUSHORT) (d[0] >> 16);
+      s[1] = (unsigned EMUSHORT) d[0];
+#if HOST_BITS_PER_WIDE_INT == 32
+      s[2] = (unsigned EMUSHORT) (d[1] >> 16);
+      s[3] = (unsigned EMUSHORT) d[1];
+#else
+      /* In this case the entire target double is contained in the
+        first array element.  The second element of the input is
+        ignored.  */
+      s[2] = (unsigned EMUSHORT) (d[0] >> 48);
+      s[3] = (unsigned EMUSHORT) (d[0] >> 32);
+#endif
+    }
+  else
+    {
+      /* Target float words are little-endian.  */
+      s[0] = (unsigned EMUSHORT) d[0];
+      s[1] = (unsigned EMUSHORT) (d[0] >> 16);
+#if HOST_BITS_PER_WIDE_INT == 32
+      s[2] = (unsigned EMUSHORT) d[1];
+      s[3] = (unsigned EMUSHORT) (d[1] >> 16);
+#else
+      s[2] = (unsigned EMUSHORT) (d[0] >> 32);
+      s[3] = (unsigned EMUSHORT) (d[0] >> 48);
+#endif
+    }
+  /* Convert target double to E-type.  */
+  e53toe (s, e);
+  /* Output E-type to REAL_VALUE_TYPE.  */
+  PUT_REAL (e, &r);
+  return r;
+}
+
+
+/* Convert target computer unsigned 64-bit integer to e-type.
+   The endian-ness of DImode follows the convention for integers,
+   so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN.  */
+
+static void
+uditoe (di, e)
+     unsigned EMUSHORT *di;  /* Address of the 64-bit int.  */
+     unsigned EMUSHORT *e;
+{
+  unsigned EMUSHORT yi[NI];
+  int k;
+
+  ecleaz (yi);
+  if (WORDS_BIG_ENDIAN)
+    {
+      for (k = M; k < M + 4; k++)
+       yi[k] = *di++;
+    }
+  else
+    {
+      for (k = M + 3; k >= M; k--)
+       yi[k] = *di++;
+    }
+  yi[E] = EXONE + 47;  /* exponent if normalize shift count were 0 */
+  if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
+    ecleaz (yi);               /* it was zero */
+  else
+    yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
+  emovo (yi, e);
+}
+
+/* Convert target computer signed 64-bit integer to e-type.  */
+
+static void
+ditoe (di, e)
+     unsigned EMUSHORT *di;  /* Address of the 64-bit int.  */
+     unsigned EMUSHORT *e;
+{
+  unsigned EMULONG acc;
+  unsigned EMUSHORT yi[NI];
+  unsigned EMUSHORT carry;
+  int k, sign;
+
+  ecleaz (yi);
+  if (WORDS_BIG_ENDIAN)
+    {
+      for (k = M; k < M + 4; k++)
+       yi[k] = *di++;
+    }
+  else
+    {
+      for (k = M + 3; k >= M; k--)
+       yi[k] = *di++;
+    }
+  /* Take absolute value */
+  sign = 0;
+  if (yi[M] & 0x8000)
+    {
+      sign = 1;
+      carry = 0;
+      for (k = M + 3; k >= M; k--)
+       {
+         acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
+         yi[k] = acc;
+         carry = 0;
+         if (acc & 0x10000)
+           carry = 1;
+       }
+    }
+  yi[E] = EXONE + 47;  /* exponent if normalize shift count were 0 */
+  if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
+    ecleaz (yi);               /* it was zero */
+  else
+    yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
+  emovo (yi, e);
+  if (sign)
+       eneg (e);
+}
+
+
+/* Convert e-type to unsigned 64-bit int.  */
+
+static void 
+etoudi (x, i)
+     unsigned EMUSHORT *x;
+     unsigned EMUSHORT *i;
+{
+  unsigned EMUSHORT xi[NI];
+  int j, k;
+
+  emovi (x, xi);
+  if (xi[0])
+    {
+      xi[M] = 0;
+      goto noshift;
+    }
+  k = (int) xi[E] - (EXONE - 1);
+  if (k <= 0)
+    {
+      for (j = 0; j < 4; j++)
+       *i++ = 0;
+      return;
+    }
+  if (k > 64)
+    {
+      for (j = 0; j < 4; j++)
+       *i++ = 0xffff;
+      if (extra_warnings)
+       warning ("overflow on truncation to integer");
+      return;
+    }
+  if (k > 16)
+    {
+      /* Shift more than 16 bits: first shift up k-16 mod 16,
+        then shift up by 16's.  */
+      j = k - ((k >> 4) << 4);
+      if (j == 0)
+       j = 16;
+      eshift (xi, j);
+      if (WORDS_BIG_ENDIAN)
+       *i++ = xi[M];
+      else
+       {
+         i += 3;
+         *i-- = xi[M];
+       }
+      k -= j;
+      do
+       {
+         eshup6 (xi);
+         if (WORDS_BIG_ENDIAN)
+           *i++ = xi[M];
+         else
+           *i-- = xi[M];
+       }
+      while ((k -= 16) > 0);
+    }
+  else
+    {
+        /* shift not more than 16 bits */
+      eshift (xi, k);
+
+noshift:
+
+      if (WORDS_BIG_ENDIAN)
+       {
+         i += 3;
+         *i-- = xi[M];
+         *i-- = 0;
+         *i-- = 0;
+         *i = 0;
+       }
+      else
+       {
+         *i++ = xi[M];
+         *i++ = 0;
+         *i++ = 0;
+         *i = 0;
+       }
+    }
+}
+
+
+/* Convert e-type to signed 64-bit int.  */
+
+static void 
+etodi (x, i)
+     unsigned EMUSHORT *x;
+     unsigned EMUSHORT *i;
+{
+  unsigned EMULONG acc;
+  unsigned EMUSHORT xi[NI];
+  unsigned EMUSHORT carry;
+  unsigned EMUSHORT *isave;
+  int j, k;
+
+  emovi (x, xi);
+  k = (int) xi[E] - (EXONE - 1);
+  if (k <= 0)
+    {
+      for (j = 0; j < 4; j++)
+       *i++ = 0;
+      return;
+    }
+  if (k > 64)
+    {
+      for (j = 0; j < 4; j++)
+       *i++ = 0xffff;
+      if (extra_warnings)
+       warning ("overflow on truncation to integer");
+      return;
+    }
+  isave = i;
+  if (k > 16)
+    {
+      /* Shift more than 16 bits: first shift up k-16 mod 16,
+        then shift up by 16's.  */
+      j = k - ((k >> 4) << 4);
+      if (j == 0)
+       j = 16;
+      eshift (xi, j);
+      if (WORDS_BIG_ENDIAN)
+       *i++ = xi[M];
+      else
+       {
+         i += 3;
+         *i-- = xi[M];
+       }
+      k -= j;
+      do
+       {
+         eshup6 (xi);
+         if (WORDS_BIG_ENDIAN)
+           *i++ = xi[M];
+         else
+           *i-- = xi[M];
+       }
+      while ((k -= 16) > 0);
+    }
+  else
+    {
+        /* shift not more than 16 bits */
+      eshift (xi, k);
+
+      if (WORDS_BIG_ENDIAN)
+       {
+         i += 3;
+         *i = xi[M];
+         *i-- = 0;
+         *i-- = 0;
+         *i = 0;
+       }
+      else
+       {
+         *i++ = xi[M];
+         *i++ = 0;
+         *i++ = 0;
+         *i = 0;
+       }
+    }
+  /* Negate if negative */
+  if (xi[0])
+    {
+      carry = 0;
+      if (WORDS_BIG_ENDIAN)
+       isave += 3;
+      for (k = 0; k < 4; k++)
+       {
+         acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
+         if (WORDS_BIG_ENDIAN)
+           *isave-- = acc;
+         else
+           *isave++ = acc;
+         carry = 0;
+         if (acc & 0x10000)
+           carry = 1;
+       }
+    }
+}
+
+
+/* Longhand square root routine.  */
+
+
+static int esqinited = 0;
+static unsigned short sqrndbit[NI];
+
+static void 
+esqrt (x, y)
+     unsigned EMUSHORT *x, *y;
+{
+  unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
+  EMULONG m, exp;
+  int i, j, k, n, nlups;
+
+  if (esqinited == 0)
+    {
+      ecleaz (sqrndbit);
+      sqrndbit[NI - 2] = 1;
+      esqinited = 1;
+    }
+  /* Check for arg <= 0 */
+  i = ecmp (x, ezero);
+  if (i <= 0)
+    {
+      if (i == -1)
+       {
+         mtherr ("esqrt", DOMAIN);
+         eclear (y);
+       }
+      else
+       emov (x, y);
+      return;
+    }
+
+#ifdef INFINITY
+  if (eisinf (x))
+    {
+      eclear (y);
+      einfin (y);
+      return;
+    }
+#endif
+  /* Bring in the arg and renormalize if it is denormal.  */
+  emovi (x, xx);
+  m = (EMULONG) xx[1];         /* local long word exponent */
+  if (m == 0)
+    m -= enormlz (xx);
+
+  /* Divide exponent by 2 */
+  m -= 0x3ffe;
+  exp = (unsigned short) ((m / 2) + 0x3ffe);
+
+  /* Adjust if exponent odd */
+  if ((m & 1) != 0)
+    {
+      if (m > 0)
+       exp += 1;
+      eshdn1 (xx);
+    }
+
+  ecleaz (sq);
+  ecleaz (num);
+  n = 8;                       /* get 8 bits of result per inner loop */
+  nlups = rndprc;
+  j = 0;
+
+  while (nlups > 0)
+    {
+      /* bring in next word of arg */
+      if (j < NE)
+       num[NI - 1] = xx[j + 3];
+      /* Do additional bit on last outer loop, for roundoff.  */
+      if (nlups <= 8)
+       n = nlups + 1;
+      for (i = 0; i < n; i++)
+       {
+         /* Next 2 bits of arg */
+         eshup1 (num);
+         eshup1 (num);
+         /* Shift up answer */
+         eshup1 (sq);
+         /* Make trial divisor */
+         for (k = 0; k < NI; k++)
+           temp[k] = sq[k];
+         eshup1 (temp);
+         eaddm (sqrndbit, temp);
+         /* Subtract and insert answer bit if it goes in */
+         if (ecmpm (temp, num) <= 0)
+           {
+             esubm (temp, num);
+             sq[NI - 2] |= 1;
+           }
+       }
+      nlups -= n;
+      j += 1;
+    }
+
+  /* Adjust for extra, roundoff loop done.  */
+  exp += (NBITS - 1) - rndprc;
+
+  /* Sticky bit = 1 if the remainder is nonzero.  */
+  k = 0;
+  for (i = 3; i < NI; i++)
+    k |= (int) num[i];
+
+  /* Renormalize and round off.  */
+  emdnorm (sq, k, 0, exp, 64);
+  emovo (sq, y);
+}

 #endif /* EMU_NON_COMPILE not defined */
 #endif /* EMU_NON_COMPILE not defined */

+\f
+/* Return the binary precision of the significand for a given
+   floating point mode.  The mode can hold an integer value
+   that many bits wide, without losing any bits.  */
+
+int
+significand_size (mode)
+     enum machine_mode mode;
+{
+
+/* Don't test the modes, but their sizes, lest this
+   code won't work for BITS_PER_UNIT != 8 .  */
+
+switch (GET_MODE_BITSIZE (mode))
+  {
+  case 32:
+    return 24;
+
+  case 64:
+#if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
+    return 53;
+#else
+#if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
+    return 56;
+#else
+#if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
+    return 56;
+#else
+    abort ();
+#endif
+#endif
+#endif
+
+  case 96:
+    return 64;
+  case 128:
+    return 113;
+
+  default:
+    abort ();
+  }
+}