(all insn and peephole patterns): Rewrite without using arm_output_asm_insn.

[official-gcc.git] / gcc / real.c

/* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
   and support for XFmode IEEE extended real floating point arithmetic.
   Copyright (C) 1993, 1994 Free Software Foundation, Inc.
   Contributed by Stephen L. Moshier (moshier@world.std.com).

This file is part of GNU CC.

GNU CC is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2, or (at your option)
any later version.

GNU CC is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
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
the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.  */

#include <stdio.h>
#include <errno.h>
#include "config.h"
#include "tree.h"

#ifndef errno
extern int errno;
#endif

/* 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, 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
that tries to use host floating point arithmetic to convert
REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
etc.  In cross-compile situations a REAL_VALUE_TYPE may not
be intelligible to the host computer's native arithmetic.

The emulator defaults to the host's floating point format so that
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  */
\f
/* Type of computer arithmetic.
   Only one of DEC, IBM, MIEEE, IBMPC, or UNK should get defined.

   `MIEEE' 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 MIEEE
   also invokes the particular XFmode (`long double' type) data
   structure used by the Motorola 680x0 series processors.

   `IBMPC' refers generally to little-endian IEEE machines. In this
   case, if LONG_DOUBLE_TYPE_SIZE has been defined to be 96, then
   IBMPC 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.  */
#ifdef REAL_ARITHMETIC

#if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
/* 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 FLOAT_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 */
#else /* it's not IEEE either */
/* UNKnown arithmetic.  We don't support this and can't go on. */
unknown arithmetic type
#define UNK 1
#endif /* not IEEE */
#endif /* not IBM */
#endif /* not VAX */

#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 == 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_WORDS_BIG_ENDIAN
#define MIEEE 1
#else /* not big-endian */
#define IBMPC 1
#endif /*  big-endian */
#else /* it's not IEEE either */
unknown arithmetic type
#define UNK 1
#endif /* not IEEE */
#endif /* not IBM */
#endif /* not VAX */

#endif /* REAL_ARITHMETIC not defined */

/* 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 NANS
#endif

/* 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,
   and another type at least twice whatever that size is. */

#if HOST_BITS_PER_CHAR >= 16
#define EMUSHORT char
#define EMUSHORT_SIZE HOST_BITS_PER_CHAR
#define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
#else
#if HOST_BITS_PER_SHORT >= 16
#define EMUSHORT short
#define EMUSHORT_SIZE HOST_BITS_PER_SHORT
#define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
#else
#if HOST_BITS_PER_INT >= 16
#define EMUSHORT int
#define EMUSHORT_SIZE HOST_BITS_PER_INT
#define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
#else
#if HOST_BITS_PER_LONG >= 16
#define EMUSHORT long
#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. */
#define EMU_NON_COMPILE
#endif
#endif
#endif
#endif

#if HOST_BITS_PER_SHORT >= EMULONG_SIZE
#define EMULONG short
#else
#if HOST_BITS_PER_INT >= EMULONG_SIZE
#define EMULONG int
#else
#if HOST_BITS_PER_LONG >= EMULONG_SIZE
#define EMULONG long
#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. */
#define EMU_NON_COMPILE
#endif
#endif
#endif
#endif


/* The host interface doesn't work if no 16-bit size exists. */
#if EMUSHORT_SIZE != 16
#define EMU_NON_COMPILE
#endif

/* OK to continue compilation. */
#ifndef EMU_NON_COMPILE

/* Construct macros to translate between REAL_VALUE_TYPE and e type.
   In GET_REAL and PUT_REAL, r and e are pointers.
   A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
   in memory, with no holes.  */

#if LONG_DOUBLE_TYPE_SIZE == 96
/* 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 (r, e, 2*NE)
#define PUT_REAL(e,r) bcopy (e, r, 2*NE)
#else /* no XFmode */
#if LONG_DOUBLE_TYPE_SIZE == 128
#define NE 10
#define MAXDECEXP 4932
#define MINDECEXP -4977
#define GET_REAL(r,e) bcopy (r, e, 2*NE)
#define PUT_REAL(e,r) bcopy (e, r, 2*NE)
#else
#define NE 6
#define MAXDECEXP 4932
#define MINDECEXP -4956
#ifdef REAL_ARITHMETIC
/* Emulator uses target format internally
   but host stores it in host endian-ness. */

#if HOST_FLOAT_WORDS_BIG_ENDIAN == FLOAT_WORDS_BIG_ENDIAN
#define GET_REAL(r,e) e53toe ((unsigned EMUSHORT*) (r), (e))
#define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))

#else /* endian-ness differs */
/* emulator uses target endian-ness internally */
#define GET_REAL(r,e)           \
do { 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 { 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)

#endif /* endian-ness differs */

#else /* not REAL_ARITHMETIC */

/* emulator uses host format */
#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 TFmode */
#endif /* no XFmode */


/* 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 unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
extern unsigned EMUSHORT elog2[], esqrt2[];

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[];
     enum machine_mode mode;
{
  unsigned long th, t;

#if FLOAT_WORDS_BIG_ENDIAN
  switch (mode)
    {

    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 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;

    default:
      abort ();
    }

#else

  /* Pack the output array without swapping. */

  switch (mode)
    {

    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 SFmode:

      /* pack the first long */
      th = (unsigned long) e[1] & 0xffff;
      t = (unsigned long) e[0] & 0xffff;
      t |= th << 16;
      x[0] = t;
      break;

    default:
      abort ();
    }

#endif
}


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

void 
earith (value, icode, r1, r2)
     REAL_VALUE_TYPE *value;
     int icode;
     REAL_VALUE_TYPE *r1;
     REAL_VALUE_TYPE *r2;
{
  unsigned EMUSHORT d1[NE], d2[NE], v[NE];
  enum tree_code code;

  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)
    {
    case PLUS_EXPR:
      eadd (d2, d1, v);
      break;

    case MINUS_EXPR:
      esub (d2, d1, v);         /* d1 - d2 */
      break;

    case MULT_EXPR:
      emul (d2, d1, v);
      break;

    case RDIV_EXPR:
#ifndef REAL_INFINITY
      if (ecmp (d2, ezero) == 0)
        {
#ifdef NANS
        enan (v, eisneg (d1) ^ eisneg (d2));
        break;
#else
        abort ();
#endif
        }
#endif
      ediv (d2, d1, v); /* d1/d2 */
      break;

    case MIN_EXPR:              /* min (d1,d2) */
      if (ecmp (d1, d2) < 0)
        emov (d1, v);
      else
        emov (d2, v);
      break;

    case MAX_EXPR:              /* max (d1,d2) */
      if (ecmp (d1, d2) > 0)
        emov (d1, v);
      else
        emov (d2, v);
      break;
    default:
      emov (ezero, v);
      break;
    }
PUT_REAL (v, value);
}


/* 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;
  HOST_WIDE_INT l;

  GET_REAL (&x, g);
#ifdef NANS
  if (eisnan (g))
    return (x);
#endif
  eifrac (g, &l, f);
  ltoe (&l, g);
  PUT_REAL (g, &r);
  return (r);
}


/* 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;
  unsigned HOST_WIDE_INT l;

  GET_REAL (&x, g);
#ifdef NANS
  if (eisnan (g))
    return (x);
#endif
  euifrac (g, &l, f);
  ultoe (&l, g);
  PUT_REAL (g, &r);
  return (r);
}


/* 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;
     enum machine_mode t;
{
  unsigned EMUSHORT tem[NE], e[NE];
  REAL_VALUE_TYPE r;

  switch (t)
    {
    case SFmode:
      asctoe24 (s, tem);
      e24toe (tem, e);
      break;
    case DFmode:
      asctoe53 (s, tem);
      e53toe (tem, e);
      break;
    case XFmode:
      asctoe64 (s, tem);
      e64toe (tem, e);
      break;
    case TFmode:
      asctoe113 (s, tem);
      e113toe (tem, e);
      break;
    default:
      asctoe (s, e);
    }
  PUT_REAL (e, &r);
  return (r);
}


/* Expansion of REAL_NEGATE.  */

REAL_VALUE_TYPE 
ereal_negate (x)
     REAL_VALUE_TYPE x;
{
  unsigned EMUSHORT e[NE];
  REAL_VALUE_TYPE r;

  GET_REAL (&x, e);
  eneg (e);
  PUT_REAL (e, &r);
  return (r);
}


/* 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];
  HOST_WIDE_INT l;

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

/* 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];
  unsigned HOST_WIDE_INT l;

  GET_REAL (&x, f);
#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.  */

void 
ereal_from_int (d, i, j)
     REAL_VALUE_TYPE *d;
     HOST_WIDE_INT i, j;
{
  unsigned EMUSHORT df[NE], dg[NE];
  HOST_WIDE_INT low, high;
  int sign;

  sign = 0;
  low = i;
  if ((high = j) < 0)
    {
      sign = 1;
      /* complement and add 1 */
      high = ~high;
      if (low)
        low = -low;
      else
        high += 1;
    }
  eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
  emul (dg, df, dg);
  ultoe ((unsigned HOST_WIDE_INT *) &low, df);
  eadd (df, dg, dg);
  if (sign)
    eneg (dg);
  PUT_REAL (dg, d);
}


/* REAL_VALUE_FROM_UNSIGNED_INT macro.   */

void 
ereal_from_uint (d, i, j)
     REAL_VALUE_TYPE *d;
     unsigned HOST_WIDE_INT i, j;
{
  unsigned EMUSHORT df[NE], dg[NE];
  unsigned HOST_WIDE_INT low, high;

  low = i;
  high = j;
  eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  ultoe (&high, dg);
  emul (dg, df, dg);
  ultoe (&low, df);
  eadd (df, dg, dg);
  PUT_REAL (dg, d);
}


/* REAL_VALUE_TO_INT macro.  */

void 
ereal_to_int (low, high, rr)
     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);
#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))
    {
      eneg (d);
      s = 1;
    }
  eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
  ediv (df, d, dg);             /* dg = d / 2^32 is the high word */
  euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
  emul (df, dh, dg);            /* fractional part is the low word */
  euifrac (dg, (unsigned HOST_WIDE_INT *)low, dh);
  if (s)
    {
      /* complement and add 1 */
      *high = ~(*high);
      if (*low)
        *low = -(*low);
      else
        *high += 1;
    }
}


/* REAL_VALUE_LDEXP macro.  */

REAL_VALUE_TYPE
ereal_ldexp (x, n)
     REAL_VALUE_TYPE x;
     int n;
{
  unsigned EMUSHORT e[NE], y[NE];
  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
   of the same names may be defined in fold-const.c.  */

#ifdef REAL_ARITHMETIC

/* Check for infinity in a REAL_VALUE_TYPE. */

int
target_isinf (x)
     REAL_VALUE_TYPE x;
{
  unsigned EMUSHORT e[NE];

#ifdef INFINITY
  GET_REAL (&x, e);
  return (eisinf (e));
#else
  return 0;
#endif
}


/* Check whether a REAL_VALUE_TYPE item is a NaN. */

int
target_isnan (x)
     REAL_VALUE_TYPE x;
{
  unsigned EMUSHORT e[NE];

#ifdef NANS
  GET_REAL (&x, e);
  return (eisnan (e));
#else
  return (0);
#endif
}


/* 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;
{
  return ereal_isneg (x);
}

/* Expansion of REAL_VALUE_TRUNCATE.
   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 arg;
{
  unsigned EMUSHORT e[NE], t[NE];
  REAL_VALUE_TYPE r;

  GET_REAL (&arg, e);
#ifdef NANS
  if (eisnan (e))
    return (arg);
#endif
  eclear (t);
  switch (mode)
    {
    case TFmode:
      etoe113 (e, t);
      e113toe (t, t);
      break;

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

    case DFmode:
      etoe53 (e, t);
      e53toe (t, t);
      break;

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

    case SImode:
      r = etrunci (arg);
      return (r);

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

    default:
      return (arg);
    }
  PUT_REAL (t, &r);
  return (r);
}

#endif /* REAL_ARITHMETIC defined */

/* 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
/* Target values are arrays of host longs. A long is guaranteed
   to be at least 32 bits wide. */

/* 128-bit long double */

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);
}

/* 80-bit long double */

void 
etarldouble (r, l)
     REAL_VALUE_TYPE r;
     long l[];
{
  unsigned EMUSHORT e[NE];

  GET_REAL (&r, e);
  etoe64 (e, e);
  endian (e, l, XFmode);
}

void 
etardouble (r, l)
     REAL_VALUE_TYPE r;
     long l[];
{
  unsigned EMUSHORT e[NE];

  GET_REAL (&r, e);
  etoe53 (e, e);
  endian (e, l, DFmode);
}

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

  GET_REAL (&r, e);
  etoe24 (e, e);
  endian (e, &l, SFmode);
  return ((long) l);
}

void
ereal_to_decimal (x, s)
     REAL_VALUE_TYPE x;
     char *s;
{
  unsigned EMUSHORT e[NE];

  GET_REAL (&x, e);
  etoasc (e, s, 20);
}

int
ereal_cmp (x, y)
     REAL_VALUE_TYPE x, y;
{
  unsigned EMUSHORT ex[NE], ey[NE];

  GET_REAL (&x, ex);
  GET_REAL (&y, ey);
  return (ecmp (ex, ey));
}

int
ereal_isneg (x)
     REAL_VALUE_TYPE x;
{
  unsigned EMUSHORT ex[NE];

  GET_REAL (&x, ex);
  return (eisneg (ex));
}

/* End of REAL_ARITHMETIC interface */
\f
/*
  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
        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 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 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. */

/* Constant definitions for math error conditions.  */

#define DOMAIN          1       /* argument domain error */
#define SING            2       /* argument singularity */
#define OVERFLOW        3       /* overflow range error */
#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 */

#if LONG_DOUBLE_TYPE_SIZE == 128
/* 0.0 */
unsigned EMUSHORT ezero[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
extern unsigned EMUSHORT ezero[];

/* 5.0E-1 */
unsigned EMUSHORT ehalf[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
extern unsigned EMUSHORT ehalf[];

/* 1.0E0 */
unsigned EMUSHORT eone[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
extern unsigned EMUSHORT eone[];

/* 2.0E0 */
unsigned EMUSHORT etwo[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
extern unsigned EMUSHORT etwo[];

/* 3.2E1 */
unsigned EMUSHORT e32[NE] =
 {0x0000, 0x0000, 0x0000, 0x0000,
  0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
extern unsigned EMUSHORT e32[];

/* 6.93147180559945309417232121458176568075500134360255E-1 */
unsigned EMUSHORT elog2[NE] =
 {0x40f3, 0xf6af, 0x03f2, 0xb398,
  0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
extern unsigned EMUSHORT elog2[];

/* 1.41421356237309504880168872420969807856967187537695E0 */
unsigned EMUSHORT esqrt2[NE] =
 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
  0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
extern unsigned EMUSHORT esqrt2[];

/* 3.14159265358979323846264338327950288419716939937511E0 */
unsigned EMUSHORT epi[NE] =
 {0x2902, 0x1cd1, 0x80dc, 0x628b,
  0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
extern unsigned EMUSHORT epi[];

#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.
   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;

/*  Clear out entire external format number.  */

static void 
eclear (x)
     register unsigned EMUSHORT *x;
{
  register int i;

  for (i = 0; i < NE; i++)
    *x++ = 0;
}



/* Move external format number from a to b.  */

static void 
emov (a, b)
     register unsigned EMUSHORT *a, *b;
{
  register int i;

  for (i = 0; i < NE; i++)
    *b++ = *a++;
}


/* Absolute value of external format number.  */

static void 
eabs (x)
     unsigned EMUSHORT x[];
{
  /* sign is top bit of last word of external format */
  x[NE - 1] &= 0x7fff;          
}

/* Negate external format number.  */

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

  x[NE - 1] ^= 0x8000;          /* Toggle the sign bit */
}



/* Return 1 if sign bit of external format number is nonzero, else zero.  */

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

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


/* Return 1 if external format number is infinity, else return zero.  */

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

#ifdef NANS
  if (eisnan (x))
    return (0);
#endif
  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.  */

static int 
eisnan (x)
     unsigned EMUSHORT x[];
{
#ifdef NANS
  int i;

  /* 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 external format number with infinity pattern (IEEE)
    or largest possible number (non-IEEE). */

static void 
einfin (x)
     register unsigned EMUSHORT *x;
{
  register int i;

#ifdef INFINITY
  for (i = 0; i < NE - 1; i++)
    *x++ = 0;
  *x |= 32767;
#else
  for (i = 0; i < NE - 1; i++)
    *x++ = 0xffff;
  *x |= 32766;
  if (rndprc < NBITS)
    {
      if (rndprc == 113)
        {
          *(x - 9) = 0;
          *(x - 8) = 0;
        }
      if (rndprc == 64)
        {
          *(x - 5) = 0;
        }
      if (rndprc == 53)
        {
          *(x - 4) = 0xf800;
        }
      else
        {
          *(x - 4) = 0;
          *(x - 3) = 0;
          *(x - 2) = 0xff00;
        }
    }
#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 external format number, converting it to internal format.  */

static void 
emovi (a, b)
     unsigned EMUSHORT *a, *b;
{
  register unsigned EMUSHORT *p, *q;
  int i;

  q = b;
  p = a + (NE - 1);             /* point to last word of external number */
  /* get the sign bit */
  if (*p & 0x8000)
    *q++ = 0xffff;
  else
    *q++ = 0;
  /* get the exponent */
  *q = *p--;
  *q++ &= 0x7fff;               /* delete the sign bit */
#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

  /* clear high guard word */
  *q++ = 0;
  /* move in the significand */
  for (i = 0; i < NE - 1; i++)
    *q++ = *p--;
  /* clear low guard word */
  *q = 0;
}


/* Move internal format number out, converting it to external format.  */

static void 
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 */
  /* combine sign and exponent */
  i = *p++;
  if (i)
    *q-- = *p++ | 0x8000;
  else
    *q-- = *p++;
#ifdef INFINITY
  if (*(p - 1) == 0x7fff)
    {
#ifdef NANS
      if (eiisnan (a))
        {
          enan (b, eiisneg (a));
          return;
        }
#endif
      einfin (b);
        return;
    }
#endif
  /* skip over guard word */
  ++p;
  /* move the significand */
  for (j = 0; j < NE - 1; j++)
    *q-- = *p++;
}

/* Clear out internal format number.  */

static void 
ecleaz (xi)
     register unsigned EMUSHORT *xi;
{
  register int i;

  for (i = 0; i < NI; i++)
    *xi++ = 0;
}


/* Same, but don't touch the sign. */

static void 
ecleazs (xi)
     register unsigned EMUSHORT *xi;
{
  register int i;

  ++xi;
  for (i = 0; i < NI - 1; i++)
    *xi++ = 0;
}



/* Move internal format number from a to b.  */

static void 
emovz (a, b)
     register unsigned EMUSHORT *a, *b;
{
  register int i;

  for (i = 0; i < NI - 1; i++)
    *b++ = *a++;
  /* clear low guard word */
  *b = 0;
}

/* Generate internal format NaN.
   The explicit pattern for this is maximum exponent and
   top two significant bits set.  */

static void
einan (x)
     unsigned EMUSHORT x[];
{

  ecleaz (x);
  x[E] = 0x7fff;
  x[M + 1] = 0xc000;
}

/* Return nonzero if internal format number 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 internal format number is nonzero.  */

static int 
eiisneg (x)
     unsigned EMUSHORT x[];
{

  return x[0] != 0;
}

/* Fill internal format number 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 internal format number 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 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;
{
  int i;

  a += M;                       /* skip up to significand area */
  b += M;
  for (i = M; i < NI; i++)
    {
      if (*a++ != *b++)
        goto difrnt;
    }
  return (0);

 difrnt:
  if (*(--a) > *(--b))
    return (1);
  else
    return (-1);
}


/* Shift significand down by 1 bit.  */

static void 
eshdn1 (x)
     register unsigned EMUSHORT *x;
{
  register unsigned EMUSHORT bits;
  int i;

  x += M;                       /* point to significand area */

  bits = 0;
  for (i = M; i < NI; i++)
    {
      if (*x & 1)
        bits |= 1;
      *x >>= 1;
      if (bits & 2)
        *x |= 0x8000;
      bits <<= 1;
      ++x;
    }
}



/* Shift significand up by 1 bit.  */

static void 
eshup1 (x)
     register unsigned EMUSHORT *x;
{
  register unsigned EMUSHORT bits;
  int i;

  x += NI - 1;
  bits = 0;

  for (i = M; i < NI; i++)
    {
      if (*x & 0x8000)
        bits |= 1;
      *x <<= 1;
      if (bits & 2)
        *x |= 1;
      bits <<= 1;
      --x;
    }
}


/* Shift significand down by 8 bits.  */

static void 
eshdn8 (x)
     register unsigned EMUSHORT *x;
{
  register unsigned EMUSHORT newbyt, oldbyt;
  int i;

  x += M;
  oldbyt = 0;
  for (i = M; i < NI; i++)
    {
      newbyt = *x << 8;
      *x >>= 8;
      *x |= oldbyt;
      oldbyt = newbyt;
      ++x;
    }
}

/* Shift significand up by 8 bits.  */

static void 
eshup8 (x)
     register unsigned EMUSHORT *x;
{
  int i;
  register unsigned EMUSHORT newbyt, oldbyt;

  x += NI - 1;
  oldbyt = 0;

  for (i = M; i < NI; i++)
    {
      newbyt = *x >> 8;
      *x <<= 8;
      *x |= oldbyt;
      oldbyt = newbyt;
      --x;
    }
}

/* Shift significand up by 16 bits.  */

static void 
eshup6 (x)
     register unsigned EMUSHORT *x;
{
  int i;
  register unsigned EMUSHORT *p;

  p = x + M;
  x += M + 1;

  for (i = M; i < NI - 1; i++)
    *p++ = *x++;

  *p = 0;
}

/* Shift significand down by 16 bits.  */

static void 
eshdn6 (x)
     register unsigned EMUSHORT *x;
{
  int i;
  register unsigned EMUSHORT *p;

  x += NI - 1;
  p = x + 1;

  for (i = M; i < NI - 1; i++)
    *(--p) = *(--x);

  *(--p) = 0;
}
\f
/* Add significands.  x + y replaces y.  */

static void 
eaddm (x, y)
     unsigned EMUSHORT *x, *y;
{
  register unsigned EMULONG a;
  int i;
  unsigned int carry;

  x += NI - 1;
  y += NI - 1;
  carry = 0;
  for (i = M; i < NI; i++)
    {
      a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
      if (a & 0x10000)
        carry = 1;
      else
        carry = 0;
      *y = (unsigned EMUSHORT) a;
      --x;
      --y;
    }
}

/* Subtract significands.  y - x replaces y.  */

static void 
esubm (x, y)
     unsigned EMUSHORT *x, *y;
{
  unsigned EMULONG a;
  int i;
  unsigned int carry;

  x += NI - 1;
  y += NI - 1;
  carry = 0;
  for (i = M; i < NI; i++)
    {
      a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
      if (a & 0x10000)
        carry = 1;
      else
        carry = 0;
      *y = (unsigned EMUSHORT) a;
      --x;
      --y;
    }
}


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 i;
  register unsigned EMUSHORT *p, *q;
  unsigned EMUSHORT j;

  p = &equot[0];
  *p++ = num[0];
  *p++ = num[1];

  for (i = M; i < NI; i++)
    {
      *p++ = 0;
    }

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

  p = &den[M + 2];
  if (*p++ == 0)
    {
      for (i = M + 3; i < NI; i++)
        {
          if (*p++ != 0)
            goto fulldiv;
        }
      if ((den[M + 1] & 1) != 0)
        goto fulldiv;
      eshdn1 (num);
      eshdn1 (den);

      p = &den[M + 1];
      q = &num[M + 1];

      for (i = 0; i < NBITS + 2; i++)
        {
          if (*p <= *q)
            {
              *q -= *p;
              j = 1;
            }
          else
            {
              j = 0;
            }
          eshup1 (equot);
          equot[NI - 2] |= j;
          eshup1 (num);
        }
      goto divdon;
    }

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

 fulldiv:

  p = &equot[NI - 2];
  for (i = 0; i < NBITS + 2; i++)
    {
      if (ecmpm (den, num) <= 0)
        {
          esubm (den, num);
          j = 1;                /* quotient bit = 1 */
        }
      else
        j = 0;
      eshup1 (equot);
      *p |= j;
      eshup1 (num);
    }

 divdon:

  eshdn1 (equot);
  eshdn1 (equot);

  /* 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 */
int 
emulm (a, b)
     unsigned EMUSHORT a[], b[];
{
  unsigned EMUSHORT *p, *q;
  int i, j, k;

  equot[0] = b[0];
  equot[1] = b[1];
  for (i = M; i < NI; i++)
    equot[i] = 0;

  p = &a[NI - 2];
  k = NBITS;
  while (*p == 0)               /* significand is not supposed to be zero */
    {
      eshdn6 (a);
      k -= 16;
    }
  if ((*p & 0xff) == 0)
    {
      eshdn8 (a);
      k -= 8;
    }

  q = &equot[NI - 1];
  j = 0;
  for (i = 0; i < k; i++)
    {
      if (*p & 1)
        eaddm (b, equot);
      /* remember if there were any nonzero bits shifted out */
      if (*q & 1)
        j |= 1;
      eshdn1 (a);
      eshdn1 (equot);
    }

  for (i = 0; i < NI; i++)
    b[i] = equot[i];

  /* return flag for lost nonzero bits */
  return (j);
}

#else

/* Radix 65536 versions of multiply and divide  */


/* Multiply significand of e-type number b
   by 16-bit quantity a, e-type result to c. */

static void
m16m (a, b, c)
     unsigned int a;
     unsigned short b[], c[];
{
  register unsigned short *pp;
  register unsigned long carry;
  unsigned short *ps;
  unsigned short p[NI];
  unsigned long aa, m;
  int i;

  aa = a;
  pp = &p[NI-2];
  *pp++ = 0;
  *pp = 0;
  ps = &b[NI-1];

  for (i=M+1; i<NI; i++)
    {
      if (*ps == 0)
        {
          --ps;
          --pp;
          *(pp-1) = 0;
        }
      else
        {
          m = (unsigned long) aa * *ps--;
          carry = (m & 0xffff) + *pp;
          *pp-- = (unsigned short)carry;
          carry = (carry >> 16) + (m >> 16) + *pp;
          *pp = (unsigned short)carry;
          *(pp-1) = carry >> 16;
        }
    }
  for (i=M; i<NI; i++)
    c[i] = p[i];
}


/* Divide significands. Neither the numerator nor the denominator
   is permitted to have its high guard word nonzero.  */

static int
edivm (den, num)
     unsigned short den[], num[];
{
  int i;
  register unsigned short *p;
  unsigned long tnum;
  unsigned short j, tdenm, tquot;
  unsigned short 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 long) 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 */
static int
emulm (a, b)
     unsigned short a[], b[];
{
  unsigned short *p, *q;
  unsigned short pprod[NI];
  unsigned short 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" 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.

  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);
      return;
    }
#endif
  if (exp < 0L)
    {
      if (exp > (EMULONG) (-NBITS - 1))
        {
          j = (int) exp;
          i = eshift (s, j);
          if (i)
            lost = 1;
        }
      else
        {
          ecleazs (s);
          return;
        }
    }
  /* Round off, unless told not to by rcntrl. */
  if (rcntrl == 0)
    goto mdfin;
  /* Set up rounding parameters if the control register changed. */
  if (rndprc != rlast)
    {
      ecleaz (rbit);
      switch (rndprc)
        {
        default:
        case NBITS:
          rw = NI - 1;          /* low guard word */
          rmsk = 0xffff;
          rmbit = 0x8000;
          re = rw - 1;
          rebit = 1;
          break;
        case 113:
          rw = 10;
          rmsk = 0x7fff;
          rmbit = 0x4000;
          rebit = 0x8000;
          re = rw;
          break;
        case 64:
          rw = 7;
          rmsk = 0xffff;
          rmbit = 0x8000;
          re = rw - 1;
          rebit = 1;
          break;
          /* For DEC or IBM arithmetic */
        case 56:
          rw = 6;
          rmsk = 0xff;
          rmbit = 0x80;
          rebit = 0x100;
          re = rw;
          break;
        case 53:
          rw = 6;
          rmsk = 0x7ff;
          rmbit = 0x0400;
          rebit = 0x800;
          re = rw;
          break;
        case 24:
          rw = 4;
          rmsk = 0xff;
          rmbit = 0x80;
          rebit = 0x100;
          re = rw;
          break;
        }
      rbit[re] = rebit;
      rlast = rndprc;
    }

  /* Shift down 1 temporarily if the data structure has an implied
     most significant bit and the number is denormal.  */
  if ((exp <= 0) && (rndprc != 64) && (rndprc != NBITS))
    {
      lost |= s[NI - 1] & 1;
      eshdn1 (s);
    }
  /* Clear out all bits below the rounding bit,
     remembering in r if any were nonzero.  */
  r = s[rw] & rmsk;
  if (rndprc < NBITS)
    {
      i = rw + 1;
      while (i < NI)
        {
          if (s[i])
            r |= 1;
          s[i] = 0;
          ++i;
        }
    }
  s[rw] &= ~rmsk;
  if ((r & rmbit) != 0)
    {
      if (r == rmbit)
        {
          if (lost == 0)
            {                   /* round to even */
              if ((s[re] & rebit) == 0)
                goto mddone;
            }
          else
            {
              if (subflg != 0)
                goto mddone;
            }
        }
      eaddm (rbit, s);
    }
 mddone:
  if ((exp <= 0) && (rndprc != 64) && (rndprc != NBITS))
    {
      eshup1 (s);
    }
  if (s[2] != 0)
    {                           /* overflow on roundoff */
      eshdn1 (s);
      exp += 1;
    }
 mdfin:
  s[NI - 1] = 0;
  if (exp >= 32767L)
    {
#ifndef INFINITY
    overf:
#endif
#ifdef INFINITY
      s[1] = 32767;
      for (i = 2; i < NI - 1; i++)
        s[i] = 0;
      if (extra_warnings)
        warning ("floating point overflow");
#else
      s[1] = 32766;
      s[2] = 0;
      for (i = M + 1; i < NI - 1; i++)
        s[i] = 0xffff;
      s[NI - 1] = 0;
      if ((rndprc < 64) || (rndprc == 113))
        {
          s[rw] &= ~rmsk;
          if (rndprc == 24)
            {
              s[5] = 0;
              s[6] = 0;
            }
        }
#endif
      return;
    }
  if (exp < 0)
    s[1] = 0;
  else
    s[1] = (unsigned EMUSHORT) exp;
}



/*  Subtract external format numbers.  */

static int subflg = 0;

static void 
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);
}


/* Add.  */

static void 
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);
}

static void 
eadd1 (a, b, c)
     unsigned EMUSHORT *a, *b, *c;
{
  unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
  int i, lost, j, k;
  EMULONG lt, lta, ltb;

#ifdef INFINITY
  if (eisinf (a))
    {
      emov (a, c);
      if (subflg)
        eneg (c);
      return;
    }
  if (eisinf (b))
    {
      emov (b, c);
      return;
    }
#endif
  emovi (a, ai);
  emovi (b, bi);
  if (subflg)
    ai[0] = ~ai[0];

  /* compare exponents */
  lta = ai[E];
  ltb = bi[E];
  lt = lta - ltb;
  if (lt > 0L)
    {                           /* put the larger number in bi */
      emovz (bi, ci);
      emovz (ai, bi);
      emovz (ci, ai);
      ltb = bi[E];
      lt = -lt;
    }
  lost = 0;
  if (lt != 0L)
    {
      if (lt < (EMULONG) (-NBITS - 1))
        goto done;              /* answer same as larger addend */
      k = (int) lt;
      lost = eshift (ai, k);    /* shift the smaller number down */
    }
  else
    {
      /* exponents were the same, so must compare significands */
      i = ecmpm (ai, bi);
      if (i == 0)
        {                       /* the numbers are identical in magnitude */
          /* if different signs, result is zero */
          if (ai[0] != bi[0])
            {
              eclear (c);
              return;
            }
          /* if same sign, result is double */
          /* double denomalized tiny number */
          if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
            {
              eshup1 (bi);
              goto done;
            }
          /* add 1 to exponent unless both are zero! */
          for (j = 1; j < NI - 1; j++)
            {
              if (bi[j] != 0)
                {
                  /* This could overflow, but let emovo take care of that. */
                  ltb += 1;
                  break;
                }
            }
          bi[E] = (unsigned EMUSHORT) ltb;
          goto done;
        }
      if (i > 0)
        {                       /* put the larger number in bi */
          emovz (bi, ci);
          emovz (ai, bi);
          emovz (ci, ai);
        }
    }
  if (ai[0] == bi[0])
    {
      eaddm (ai, bi);
      subflg = 0;
    }
  else
    {
      esubm (ai, bi);
      subflg = 1;
    }
  emdnorm (bi, lost, subflg, ltb, 64);

 done:
  emovo (bi, c);
}



/* Divide.  */

static void 
ediv (a, b, c)
     unsigned EMUSHORT *a, *b, *c;
{
  unsigned EMUSHORT ai[NI], bi[NI];
  int i;
  EMULONG lt, lta, ltb;

#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, eisneg (a) ^ eisneg (b));
    return;
    }
#endif
/* Infinity over anything else is infinity. */
#ifdef INFINITY
  if (eisinf (b))
    {
      if (eisneg (a) ^ eisneg (b))
        *(c + (NE - 1)) = 0x8000;
      else
        *(c + (NE - 1)) = 0;
      einfin (c);
      return;
    }
/* Anything else over infinity is zero. */
  if (eisinf (a))
    {
      eclear (c);
      return;
    }
#endif
  emovi (a, ai);
  emovi (b, bi);
  lta = ai[E];
  ltb = bi[E];
  if (bi[E] == 0)
    {                           /* See if numerator is zero. */
      for (i = 1; i < NI - 1; i++)
        {
          if (bi[i] != 0)
            {
              ltb -= enormlz (bi);
              goto dnzro1;
            }
        }
      eclear (c);
      return;
    }
 dnzro1:

  if (ai[E] == 0)
    {                           /* possible divide by zero */
      for (i = 1; i < NI - 1; i++)
        {
          if (ai[i] != 0)
            {
              lta -= enormlz (ai);
              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);
      return;
    }
 dnzro2:

  i = edivm (ai, bi);
  /* 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);
}



/* Multiply.  */

static void 
emul (a, b, c)
     unsigned EMUSHORT *a, *b, *c;
{
  unsigned EMUSHORT ai[NI], bi[NI];
  int i, j;
  EMULONG lt, lta, ltb;

#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, eisneg (a) ^ eisneg (b));
    return;
    }
#endif
/* Infinity times anything else is infinity. */
#ifdef INFINITY
  if (eisinf (a) || eisinf (b))
    {
      if (eisneg (a) ^ eisneg (b))
        *(c + (NE - 1)) = 0x8000;
      else
        *(c + (NE - 1)) = 0;
      einfin (c);
      return;
    }
#endif
  emovi (a, ai);
  emovi (b, bi);
  lta = ai[E];
  ltb = bi[E];
  if (ai[E] == 0)
    {
      for (i = 1; i < NI - 1; i++)
        {
          if (ai[i] != 0)
            {
              lta -= enormlz (ai);
              goto mnzer1;
            }
        }
      eclear (c);
      return;
    }
 mnzer1:

  if (bi[E] == 0)
    {
      for (i = 1; i < NI - 1; i++)
        {
          if (bi[i] != 0)
            {
              ltb -= enormlz (bi);
              goto mnzer2;
            }
        }
      eclear (c);
      return;
    }
 mnzer2:

  /* Multiply significands */
  j = emulm (ai, bi);
  /* 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);
}




/* Convert IEEE double precision to e type.  */

static void
e53toe (pe, y)
     unsigned EMUSHORT *pe, *y;
{
#ifdef DEC

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

#else
#ifdef IBM

  ibmtoe (pe, y, DFmode);

#else
  register unsigned EMUSHORT r;
  register unsigned EMUSHORT *e, *p;
  unsigned EMUSHORT yy[NI];
  int denorm, k;

  e = pe;
  denorm = 0;                   /* flag if denormalized number */
  ecleaz (yy);
#ifdef IBMPC
  e += 3;
#endif
  r = *e;
  yy[0] = 0;
  if (r & 0x8000)
    yy[0] = 0xffff;
  yy[M] = (r & 0x0f) | 0x10;
  r &= ~0x800f;                 /* strip sign and 4 significand bits */
#ifdef INFINITY
  if (r == 0x7ff0)
    {
#ifdef NANS
#ifdef IBMPC
      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
#endif  /* NANS */
      eclear (y);
      einfin (y);
      if (yy[0])
        eneg (y);
      return;
    }
#endif  /* INFINITY */
  r >>= 4;
  /* If zero exponent, then the significand is denormalized.
     So take back the understood high significand bit. */

  if (r == 0)
    {
      denorm = 1;
      yy[M] &= ~0x10;
    }
  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++;
#endif
  eshift (yy, -5);
  if (denorm)
    {                           /* if zero exponent, then normalize the significand */
      if ((k = enormlz (yy)) > NBITS)
        ecleazs (yy);
      else
        yy[E] -= (unsigned EMUSHORT) (k - 1);
    }
  emovo (yy, y);
#endif /* not IBM */
#endif /* not DEC */
}

static void 
e64toe (pe, y)
     unsigned EMUSHORT *pe, *y;
{
  unsigned EMUSHORT yy[NI];
  unsigned EMUSHORT *e, *p, *q;
  int i;

  e = pe;
  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 IBM
  p = &yy[0] + (NE - 1);
  *p-- = *e++;
  ++e;
  for (i = 0; i < 5; i++)
    *p-- = *e++;
#endif
#ifdef MIEEE
  p = &yy[0] + (NE - 1);
  *p-- = *e++;
  ++e;
  for (i = 0; i < 4; i++)
    *p-- = *e++;
#endif
  p = yy;
  q = y;
#ifdef INFINITY
  if (*p == 0x7fff)
    {
#ifdef NANS
#ifdef IBMPC
      for (i = 0; i < 4; i++)
        {
          if (pe[i] != 0)
            {
              enan (y, (*p & 0x8000) != 0);
              return;
            }
        }
#else
      for (i = 1; i <= 4; i++)
        {
          if (pe[i] != 0)
            {
              enan (y, (*p & 0x8000) != 0);
              return;
            }
        }
#endif
#endif /* NANS */
      eclear (y);
      einfin (y);
      if (*p & 0x8000)
        eneg (y);
      return;
    }
#endif  /* INFINITY */
  for (i = 0; i < NE; i++)
    *q++ = *p++;
}


static void 
e113toe (pe, y)
     unsigned EMUSHORT *pe, *y;
{
  register unsigned EMUSHORT r;
  unsigned EMUSHORT *e, *p;
  unsigned EMUSHORT yy[NI];
  int denorm, i;

  e = pe;
  denorm = 0;
  ecleaz (yy);
#ifdef IBMPC
  e += 7;
#endif
  r = *e;
  yy[0] = 0;
  if (r & 0x8000)
    yy[0] = 0xffff;
  r &= 0x7fff;
#ifdef INFINITY
  if (r == 0x7fff)
    {
#ifdef NANS
#ifdef IBMPC
      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
#endif /* NANS */
      eclear (y);
      einfin (y);
      if (yy[0])
        eneg (y);
      return;
    }
#endif  /* INFINITY */
  yy[E] = r;
  p = &yy[M + 1];
#ifdef IBMPC
  for (i = 0; i < 7; i++)
    *p++ = *(--e);
#endif
#ifdef MIEEE
  ++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 IEEE single precision to e type.  */

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;

  e = pe;
  denorm = 0;                   /* flag if denormalized number */
  ecleaz (yy);
#ifdef IBMPC
  e += 1;
#endif
#ifdef DEC
  e += 1;
#endif
  r = *e;
  yy[0] = 0;
  if (r & 0x8000)
    yy[0] = 0xffff;
  yy[M] = (r & 0x7f) | 0200;
  r &= ~0x807f;                 /* strip sign and 7 significand bits */
#ifdef INFINITY
  if (r == 0x7f80)
    {
#ifdef NANS
#ifdef MIEEE
      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
#endif  /* NANS */
      eclear (y);
      einfin (y);
      if (yy[0])
        eneg (y);
      return;
    }
#endif  /* INFINITY */
  r >>= 7;
  /* If zero exponent, then the significand is denormalized.
     So take back the understood high significand bit. */
  if (r == 0)
    {
      denorm = 1;
      yy[M] &= ~0200;
    }
  r += EXONE - 0177;
  yy[E] = r;
  p = &yy[M + 1];
#ifdef IBMPC
  *p++ = *(--e);
#endif
#ifdef DEC
  *p++ = *(--e);
#endif
#ifdef MIEEE
  ++e;
  *p++ = *e++;
#endif
  eshift (yy, -8);
  if (denorm)
    {                           /* if zero exponent, then normalize the significand */
      if ((k = enormlz (yy)) > NBITS)
        ecleazs (yy);
      else
        yy[E] -= (unsigned EMUSHORT) (k - 1);
    }
  emovo (yy, y);
#endif /* not IBM */
}


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);
}

/* Move out internal format to ieee long double */

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;
#ifdef MIEEE
  q = b;
#else
  q = b + 7;                    /* point to output exponent */
#endif

  /* If not denormal, delete the implied bit. */
  if (a[E] != 0)
    {
      eshup1 (a);
    }
  /* combine sign and exponent */
  i = *p++;
#ifdef MIEEE
  if (i)
    *q++ = *p++ | 0x8000;
  else
    *q++ = *p++;
#else
  if (i)
    *q-- = *p++ | 0x8000;
  else
    *q-- = *p++;
#endif
  /* skip over guard word */
  ++p;
  /* move the significand */
#ifdef MIEEE
  for (i = 0; i < 7; i++)
    *q++ = *p++;
#else
  for (i = 0; i < 7; i++)
    *q-- = *p++;
#endif
}

static void 
etoe64 (x, e)
     unsigned EMUSHORT *x, *e;
{
  unsigned EMUSHORT xi[NI];
  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];
#ifdef INFINITY
  if (eisinf (x))
    goto nonorm;
#endif
  /* round off to nearest or even */
  rndsav = rndprc;
  rndprc = 64;
  emdnorm (xi, 0, 0, exp, 64);
  rndprc = rndsav;
 nonorm:
  toe64 (xi, e);
}


/* Move out internal format to ieee long double. */

static void 
toe64 (a, b)
     unsigned EMUSHORT *a, *b;
{
  register unsigned EMUSHORT *p, *q;
  unsigned EMUSHORT i;

#ifdef NANS
  if (eiisnan (a))
    {
      make_nan (b, eiisneg (a), XFmode);
      return;
    }
#endif
  p = a;
#if defined(MIEEE) || defined(IBM)
  q = b;
#else
  q = b + 4;                    /* point to output exponent */
#if LONG_DOUBLE_TYPE_SIZE == 96
  /* Clear the last two bytes of 12-byte Intel format */
  *(q+1) = 0;
#endif
#endif

  /* combine sign and exponent */
  i = *p++;
#if defined(MIEEE) || defined(IBM)
  if (i)
    *q++ = *p++ | 0x8000;
  else
    *q++ = *p++;
  *q++ = 0;
#else
  if (i)
    *q-- = *p++ | 0x8000;
  else
    *q-- = *p++;
#endif
  /* skip over guard word */
  ++p;
  /* move the significand */
#if defined(MIEEE) || defined(IBM)
  for (i = 0; i < 4; i++)
    *q++ = *p++;
#else
  for (i = 0; i < 4; i++)
    *q-- = *p++;
#endif
}


/* e type to IEEE double precision.  */

#ifdef DEC

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

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

#else
#ifdef IBM

static void 
etoe53 (x, e)
     unsigned EMUSHORT *x, *e;
{
  etoibm (x, e, DFmode);
}

static void 
toe53 (x, y)
     unsigned EMUSHORT *x, *y;
{
  toibm (x, y, DFmode);
}

#else  /* it's neither DEC nor IBM */

static void 
etoe53 (x, e)
     unsigned EMUSHORT *x, *e;
{
  unsigned EMUSHORT xi[NI];
  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);
#ifdef INFINITY
  if (eisinf (x))
    goto nonorm;
#endif
  /* round off to nearest or even */
  rndsav = rndprc;
  rndprc = 53;
  emdnorm (xi, 0, 0, exp, 64);
  rndprc = rndsav;
 nonorm:
  toe53 (xi, e);
}


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

#ifdef NANS
  if (eiisnan (x))
    {
      make_nan (y, eiisneg (x), DFmode);
      return;
    }
#endif
  p = &x[0];
#ifdef IBMPC
  y += 3;
#endif
  *y = 0;                       /* output high order */
  if (*p++)
    *y = 0x8000;                /* output sign bit */

  i = *p++;
  if (i >= (unsigned int) 2047)
    {                           /* Saturate at largest number less than infinity. */
#ifdef INFINITY
      *y |= 0x7ff0;
#ifdef IBMPC
      *(--y) = 0;
      *(--y) = 0;
      *(--y) = 0;
#endif
#ifdef MIEEE
      ++y;
      *y++ = 0;
      *y++ = 0;
      *y++ = 0;
#endif
#else
      *y |= (unsigned EMUSHORT) 0x7fef;
#ifdef IBMPC
      *(--y) = 0xffff;
      *(--y) = 0xffff;
      *(--y) = 0xffff;
#endif
#ifdef MIEEE
      ++y;
      *y++ = 0xffff;
      *y++ = 0xffff;
      *y++ = 0xffff;
#endif
#endif
      return;
    }
  if (i == 0)
    {
      eshift (x, 4);
    }
  else
    {
      i <<= 4;
      eshift (x, 5);
    }
  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
}

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



/* e type to IEEE single precision.  */

#ifdef IBM

static void 
etoe24 (x, e)
     unsigned EMUSHORT *x, *e;
{
  etoibm (x, e, SFmode);
}

static void 
toe24 (x, y)
     unsigned EMUSHORT *x, *y;
{
  toibm (x, y, SFmode);
}

#else

static void 
etoe24 (x, e)
     unsigned EMUSHORT *x, *e;
{
  EMULONG exp;
  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);
#ifdef INFINITY
  if (eisinf (x))
    goto nonorm;
#endif
  /* round off to nearest or even */
  rndsav = rndprc;
  rndprc = 24;
  emdnorm (xi, 0, 0, exp, 64);
  rndprc = rndsav;
 nonorm:
  toe24 (xi, e);
}

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

#ifdef NANS
  if (eiisnan (x))
    {
      make_nan (y, eiisneg (x), SFmode);
      return;
    }
#endif
  p = &x[0];
#ifdef IBMPC
  y += 1;
#endif
#ifdef DEC
  y += 1;
#endif
  *y = 0;                       /* output high order */
  if (*p++)
    *y = 0x8000;                /* output sign bit */

  i = *p++;
/* Handle overflow cases. */
  if (i >= 255)
    {
#ifdef INFINITY
      *y |= (unsigned EMUSHORT) 0x7f80;
#ifdef IBMPC
      *(--y) = 0;
#endif
#ifdef DEC
      *(--y) = 0;
#endif
#ifdef MIEEE
      ++y;
      *y = 0;
#endif
#else  /* no INFINITY */
      *y |= (unsigned EMUSHORT) 0x7f7f;
#ifdef IBMPC
      *(--y) = 0xffff;
#endif
#ifdef DEC
      *(--y) = 0xffff;
#endif
#ifdef MIEEE
      ++y;
      *y = 0xffff;
#endif
#ifdef ERANGE
      errno = ERANGE;
#endif
#endif  /* no INFINITY */
      return;
    }
  if (i == 0)
    {
      eshift (x, 7);
    }
  else
    {
      i <<= 7;
      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
#ifdef DEC
  *(--y) = *p;
#endif
#ifdef MIEEE
  ++y;
  *y = *p;
#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.  */

static int 
ecmp (a, b)
     unsigned EMUSHORT *a, *b;
{
  unsigned EMUSHORT ai[NI], bi[NI];
  register unsigned EMUSHORT *p, *q;
  register int i;
  int msign;

#ifdef NANS
  if (eisnan (a)  || eisnan (b))
      return (-2);
#endif
  emovi (a, ai);
  p = ai;
  emovi (b, bi);
  q = bi;

  if (*p != *q)
    {                           /* the signs are different */
      /* -0 equals + 0 */
      for (i = 1; i < NI - 1; i++)
        {
          if (ai[i] != 0)
            goto nzro;
          if (bi[i] != 0)
            goto nzro;
        }
      return (0);
    nzro:
      if (*p == 0)
        return (1);
      else
        return (-1);
    }
  /* both are the same sign */
  if (*p == 0)
    msign = 1;
  else
    msign = -1;
  i = NI - 1;
  do
    {
      if (*p++ != *q++)
        {
          goto diff;
        }
    }
  while (--i > 0);

  return (0);                   /* equality */



 diff:

  if (*(--p) > *(--q))
    return (msign);             /* p is bigger */
  else
    return (-msign);            /* p is littler */
}




/* Find nearest integer to x = floor (x + 0.5).  */

static void 
eround (x, y)
     unsigned EMUSHORT *x, *y;
{
  eadd (ehalf, x, y);
  efloor (y, y);
}




/* Convert HOST_WIDE_INT to e type.  */

static void 
ltoe (lp, y)
     HOST_WIDE_INT *lp;
     unsigned EMUSHORT *y;
{
  unsigned EMUSHORT yi[NI];
  unsigned HOST_WIDE_INT ll;
  int k;

  ecleaz (yi);
  if (*lp < 0)
    {
      /* make it positive */
      ll = (unsigned HOST_WIDE_INT) (-(*lp));
      yi[0] = 0xffff;           /* put correct sign in the e type number */
    }
  else
    {
      ll = (unsigned HOST_WIDE_INT) (*lp);
    }
  /* 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[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
    yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
  emovo (yi, y);                /* output the answer */
}

/* Convert unsigned HOST_WIDE_INT to e type.  */

static void 
ultoe (lp, y)
     unsigned HOST_WIDE_INT *lp;
     unsigned EMUSHORT *y;
{
  unsigned EMUSHORT yi[NI];
  unsigned HOST_WIDE_INT ll;
  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[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
    yi[E] -= (unsigned EMUSHORT) k;  /* subtract shift count from exponent */
  emovo (yi, y);                /* output the answer */
}


/* Find signed HOST_WIDE_INT integer and floating point fractional
   parts 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).  */

static void 
eifrac (x, i, frac)
     unsigned EMUSHORT *x;
     HOST_WIDE_INT *i;
     unsigned EMUSHORT *frac;
{
  unsigned EMUSHORT xi[NI];
  int j, k;
  unsigned HOST_WIDE_INT ll;

  emovi (x, xi);
  k = (int) xi[E] - (EXONE - 1);
  if (k <= 0)
    {
      /* if exponent <= 0, integer = 0 and real output is fraction */
      *i = 0L;
      emovo (xi, frac);
      return;
    }
  if (k > (HOST_BITS_PER_WIDE_INT - 1))
    {
      /* long integer overflow: output large integer
         and correct fraction  */
      if (xi[0])
        *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
      else
        {
#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");
    }
  else 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);
      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);
    }
  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;
  if ((k = enormlz (xi)) > NBITS)
    ecleaz (xi);
  else
    xi[E] -= (unsigned EMUSHORT) k;

  emovo (xi, frac);
}


/* Find unsigned HOST_WIDE_INT integer and floating point fractional parts.
   A negative e type input yields integer output = 0
   but correct fraction.  */

static void 
euifrac (x, i, frac)
     unsigned EMUSHORT *x;
     unsigned HOST_WIDE_INT *i;
     unsigned EMUSHORT *frac;
{
  unsigned HOST_WIDE_INT ll;
  unsigned EMUSHORT xi[NI];
  int j, k;

  emovi (x, xi);
  k = (int) xi[E] - (EXONE - 1);
  if (k <= 0)
    {
      /* if exponent <= 0, integer = 0 and argument is fraction */
      *i = 0L;
      emovo (xi, frac);
      return;
    }
  if (k > HOST_BITS_PER_WIDE_INT)
    {
      /* 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");
    }
  else 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);
      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);
      *i = (HOST_WIDE_INT) xi[M] & 0xffff;
    }

  if (xi[0])  /* A negative value yields unsigned integer 0. */
    *i = 0L;

  xi[0] = 0;
  xi[E] = EXONE - 1;
  xi[M] = 0;
  if ((k = enormlz (xi)) > NBITS)
    ecleaz (xi);
  else
    xi[E] -= (unsigned EMUSHORT) k;

  emovo (xi, frac);
}



/* Shift significand area up or down by the number of bits given by SC.  */

static int 
eshift (x, sc)
     unsigned EMUSHORT *x;
     int sc;
{
  unsigned EMUSHORT lost;
  unsigned EMUSHORT *p;

  if (sc == 0)
    return (0);

  lost = 0;
  p = x + NI - 1;

  if (sc < 0)
    {
      sc = -sc;
      while (sc >= 16)
        {
          lost |= *p;           /* remember lost bits */
          eshdn6 (x);
          sc -= 16;
        }

      while (sc >= 8)
        {
          lost |= *p & 0xff;
          eshdn8 (x);
          sc -= 8;
        }

      while (sc > 0)
        {
          lost |= *p & 1;
          eshdn1 (x);
          sc -= 1;
        }
    }
  else
    {
      while (sc >= 16)
        {
          eshup6 (x);
          sc -= 16;
        }

      while (sc >= 8)
        {
          eshup8 (x);
          sc -= 8;
        }

      while (sc > 0)
        {
          eshup1 (x);
          sc -= 1;
        }
    }
  if (lost)
    lost = 1;
  return ((int) lost);
}



/* Shift normalize the significand area pointed to by argument.
   Shift count (up = positive) is returned.  */

static int 
enormlz (x)
     unsigned EMUSHORT x[];
{
  register unsigned EMUSHORT *p;
  int sc;

  sc = 0;
  p = &x[M];
  if (*p != 0)
    goto normdn;
  ++p;
  if (*p & 0x8000)
    return (0);                 /* already normalized */
  while (*p == 0)
    {
      eshup6 (x);
      sc += 16;

      /* With guard word, there are NBITS+16 bits available.
       Return true if all are zero.  */
      if (sc > NBITS)
        return (sc);
    }
  /* see if high byte is zero */
  while ((*p & 0xff00) == 0)
    {
      eshup8 (x);
      sc += 8;
    }
  /* now shift 1 bit at a time */
  while ((*p & 0x8000) == 0)
    {
      eshup1 (x);
      sc += 1;
      if (sc > NBITS)
        {
          mtherr ("enormlz", UNDERFLOW);
          return (sc);
        }
    }
  return (sc);

  /* Normalize by shifting down out of the high guard word
     of the significand */
 normdn:

  if (*p & 0xff00)
    {
      eshdn8 (x);
      sc -= 8;
    }
  while (*p != 0)
    {
      eshdn1 (x);
      sc -= 1;

      if (sc < -NBITS)
        {
          mtherr ("enormlz", OVERFLOW);
          return (sc);
        }
    }
  return (sc);
}




/* Convert e type number to decimal format ASCII string.
   The constants are for 64 bit precision.  */

#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 */
  {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,},    /* 10**2048 */
  {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
  {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
  {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
  {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
  {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
  {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
  {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
  {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
  {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
  {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
  {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,},    /* 10**1 */
};

static unsigned EMUSHORT emtens[NTEN + 1][NE] =
{
  {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,},    /* 10**-4096 */
  {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,},    /* 10**-2048 */
  {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
  {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
  {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
  {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
  {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
  {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
  {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
  {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
  {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
  {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
  {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,},    /* 10**-1 */
};
#endif

static void 
e24toasc (x, string, ndigs)
     unsigned EMUSHORT x[];
     char *string;
     int ndigs;
{
  unsigned EMUSHORT w[NI];

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


static void 
e53toasc (x, string, ndigs)
     unsigned EMUSHORT x[];
     char *string;
     int ndigs;
{
  unsigned EMUSHORT w[NI];

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


static void 
e64toasc (x, string, ndigs)
     unsigned EMUSHORT x[];
     char *string;
     int ndigs;
{
  unsigned EMUSHORT w[NI];

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

static void 
e113toasc (x, string, ndigs)
     unsigned EMUSHORT x[];
     char *string;
     int ndigs;
{
  unsigned EMUSHORT w[NI];

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


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

static void 
etoasc (x, string, ndigs)
     unsigned EMUSHORT x[];
     char *string;
     int ndigs;
{
  EMUSHORT digit;
  unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
  unsigned EMUSHORT *p, *r, *ten;
  unsigned EMUSHORT sign;
  int i, j, k, expon, rndsav;
  char *s, *ss;
  unsigned EMUSHORT m;


  rndsav = rndprc;
  ss = string;
  s = wstring;
  *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)
    {
      sign = 0xffff;
      y[NE - 1] &= 0x7fff;
    }
  else
    {
      sign = 0;
    }
  expon = 0;
  ten = &etens[NTEN][0];
  emov (eone, t);
  /* Test for zero exponent */
  if (y[NE - 1] == 0)
    {
      for (k = 0; k < NE - 1; k++)
        {
          if (y[k] != 0)
            goto tnzro;         /* denormalized number */
        }
      goto isone;               /* legal all zeros */
    }
 tnzro:

  /* Test for infinity. */
  if (y[NE - 1] == 0x7fff)
    {
      if (sign)
        sprintf (wstring, " -Infinity ");
      else
        sprintf (wstring, " Infinity ");
      goto bxit;
    }

  /* Test for exponent nonzero but significand denormalized.
   * This is an error condition.
   */
  if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
    {
      mtherr ("etoasc", DOMAIN);
      sprintf (wstring, "NaN");
      goto bxit;
    }

  /* Compare to 1.0 */
  i = ecmp (eone, y);
  if (i == 0)
    goto isone;

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

  if (i < 0)
    {                           /* Number is greater than 1 */
      /* Convert significand to an integer and strip trailing decimal zeros. */
      emov (y, u);
      u[NE - 1] = EXONE + NBITS - 1;

      p = &etens[NTEN - 4][0];
      m = 16;
      do
        {
          ediv (p, u, t);
          efloor (t, w);
          for (j = 0; j < NE - 1; j++)
            {
              if (t[j] != w[j])
                goto noint;
            }
          emov (t, u);
          expon += (int) m;
        noint:
          p += NE;
          m >>= 1;
        }
      while (m != 0);

      /* Rescale from integer significand */
      u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
      emov (u, y);
      /* Find power of 10 */
      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)
            {
              ediv (p, u, u);
              emul (p, t, t);
              expon += (int) m;
            }
          m >>= 1;
          if (m == 0)
            break;
          p += NE;
        }
    }
  else
    {                           /* Number is less than 1.0 */
      /* Pad significand with trailing decimal zeros. */
      if (y[NE - 1] == 0)
        {
          while ((y[NE - 2] & 0x8000) == 0)
            {
              emul (ten, y, y);
              expon -= 1;
            }
        }
      else
        {
          emovi (y, w);
          for (i = 0; i < NDEC + 1; i++)
            {
              if ((w[NI - 1] & 0x7) != 0)
                break;
              /* multiply by 10 */
              emovz (w, u);
              eshdn1 (u);
              eshdn1 (u);
              eaddm (w, u);
              u[1] += 3;
              while (u[2] != 0)
                {
                  eshdn1 (u);
                  u[1] += 1;
                }
              if (u[NI - 1] != 0)
                break;
              if (eone[NE - 1] <= u[1])
                break;
              emovz (u, w);
              expon -= 1;
            }
          emovo (w, y);
        }
      k = -MAXP;
      p = &emtens[0][0];
      r = &etens[0][0];
      emov (y, w);
      emov (eone, t);
      while (ecmp (eone, w) > 0)
        {
          if (ecmp (p, w) >= 0)
            {
              emul (r, w, w);
              emul (r, t, t);
              expon += k;
            }
          k /= 2;
          if (k == 0)
            break;
          p += NE;
          r += NE;
        }
      ediv (t, eone, t);
    }
 isone:
  /* Find the first (leading) digit. */
  emovi (t, w);
  emovz (w, t);
  emovi (y, w);
  emovz (w, y);
  eiremain (t, y);
  digit = equot[NI - 1];
  while ((digit == 0) && (ecmp (y, ezero) != 0))
    {
      eshup1 (y);
      emovz (y, u);
      eshup1 (u);
      eshup1 (u);
      eaddm (u, y);
      eiremain (t, y);
      digit = equot[NI - 1];
      expon -= 1;
    }
  s = wstring;
  if (sign)
    *s++ = '-';
  else
    *s++ = ' ';
  /* Examine number of digits requested by caller. */
  if (ndigs < 0)
    ndigs = 0;
  if (ndigs > NDEC)
    ndigs = NDEC;
  if (digit == 10)
    {
      *s++ = '1';
      *s++ = '.';
      if (ndigs > 0)
        {
          *s++ = '0';
          ndigs -= 1;
        }
      expon += 1;
    }
  else
    {
      *s++ = (char)digit + '0';
      *s++ = '.';
    }
  /* Generate digits after the decimal point. */
  for (k = 0; k <= ndigs; k++)
    {
      /* multiply current number by 10, without normalizing */
      eshup1 (y);
      emovz (y, u);
      eshup1 (u);
      eshup1 (u);
      eaddm (u, y);
      eiremain (t, y);
      *s++ = (char) equot[NI - 1] + '0';
    }
  digit = equot[NI - 1];
  --s;
  ss = s;
  /* round off the ASCII string */
  if (digit > 4)
    {
      /* Test for critical rounding case in ASCII output. */
      if (digit == 5)
        {
          emovo (y, t);
          if (ecmp (t, ezero) != 0)
            goto roun;          /* round to nearest */
          if ((*(s - 1) & 1) == 0)
            goto doexp;         /* round to even */
        }
      /* Round up and propagate carry-outs */
    roun:
      --s;
      k = *s & 0x7f;
      /* Carry out to most significant digit? */
      if (k == '.')
        {
          --s;
          k = *s;
          k += 1;
          *s = (char) k;
          /* Most significant digit carries to 10? */
          if (k > '9')
            {
              expon += 1;
              *s = '1';
            }
          goto doexp;
        }
      /* Round up and carry out from less significant digits */
      k += 1;
      *s = (char) k;
      if (k > '9')
        {
          *s = '0';
          goto roun;
        }
    }
 doexp:
  /*
     if (expon >= 0)
     sprintf (ss, "e+%d", expon);
     else
     sprintf (ss, "e%d", expon);
     */
  sprintf (ss, "e%d", expon);
 bxit:
  rndprc = rndsav;
  /* copy out the working string */
  s = string;
  ss = wstring;
  while (*ss == ' ')            /* strip possible leading space */
    ++ss;
  while ((*s++ = *ss++) != '\0')
    ;
}


/* 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).  */

/* ASCII to single */

static void 
asctoe24 (s, y)
     char *s;
     unsigned EMUSHORT *y;
{
  asctoeg (s, y, 24);
}


/* ASCII to double */

static void 
asctoe53 (s, y)
     char *s;
     unsigned EMUSHORT *y;
{
#if defined(DEC) || defined(IBM)
  asctoeg (s, y, 56);
#else
  asctoeg (s, y, 53);
#endif
}


/* ASCII to long double */

static void 
asctoe64 (s, y)
     char *s;
     unsigned EMUSHORT *y;
{
  asctoeg (s, y, 64);
}

/* ASCII to 128-bit long double */

static void 
asctoe113 (s, y)
     char *s;
     unsigned EMUSHORT *y;
{
  asctoeg (s, y, 113);
}

/* ASCII to super double */

static void 
asctoe (s, y)
     char *s;
     unsigned EMUSHORT *y;
{
  asctoeg (s, y, NBITS);
}


/* ASCII to e type, with specified rounding precision = oprec. */

static void 
asctoeg (ss, y, oprec)
     char *ss;
     unsigned EMUSHORT *y;
     int oprec;
{
  unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
  int esign, decflg, sgnflg, nexp, exp, prec, lost;
  int k, trail, c, rndsav;
  EMULONG lexp;
  unsigned EMUSHORT nsign, *p;
  char *sp, *s, *lstr;

  /* Copy the input string. */
  lstr = (char *) alloca (strlen (ss) + 1);
  s = ss;
  while (*s == ' ')             /* skip leading spaces */
    ++s;
  sp = lstr;
  while ((*sp++ = *s++) != '\0')
    ;
  s = lstr;

  rndsav = rndprc;
  rndprc = NBITS;               /* Set to full precision */
  lost = 0;
  nsign = 0;
  decflg = 0;
  sgnflg = 0;
  nexp = 0;
  exp = 0;
  prec = 0;
  ecleaz (yy);
  trail = 0;

 nxtcom:
  k = *s - '0';
  if ((k >= 0) && (k <= 9))
    {
      /* Ignore leading zeros */
      if ((prec == 0) && (decflg == 0) && (k == 0))
        goto donchr;
      /* Identify and strip trailing zeros after the decimal point. */
      if ((trail == 0) && (decflg != 0))
        {
          sp = s;
          while ((*sp >= '0') && (*sp <= '9'))
            ++sp;
          /* Check for syntax error */
          c = *sp & 0x7f;
          if ((c != 'e') && (c != 'E') && (c != '\0')
              && (c != '\n') && (c != '\r') && (c != ' ')
              && (c != ','))
            goto error;
          --sp;
          while (*sp == '0')
            *sp-- = 'z';
          trail = 1;
          if (*s == 'z')
            goto donchr;
        }

      /* 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.  */

      if (yy[2] == 0)
        {
          if (decflg)
            nexp += 1;          /* count digits after decimal point */
          eshup1 (yy);          /* multiply current number by 10 */
          emovz (yy, xt);
          eshup1 (xt);
          eshup1 (xt);
          eaddm (xt, yy);
          ecleaz (xt);
          xt[NI - 2] = (unsigned EMUSHORT) k;
          eaddm (xt, yy);
        }
      else
        {
          /* Mark any lost non-zero digit.  */
          lost |= k;
          /* Count lost digits before the decimal point.  */
          if (decflg == 0)
            nexp -= 1;
        }
      prec += 1;
      goto donchr;
    }

  switch (*s)
    {
    case 'z':
      break;
    case 'E':
    case 'e':
      goto expnt;
    case '.':                   /* decimal point */
      if (decflg)
        goto error;
      ++decflg;
      break;
    case '-':
      nsign = 0xffff;
      if (sgnflg)
        goto error;
      ++sgnflg;
      break;
    case '+':
      if (sgnflg)
        goto error;
      ++sgnflg;
      break;
    case ',':
    case ' ':
    case '\0':
    case '\n':
    case '\r':
      goto daldone;
    case 'i':
    case 'I':
      goto infinite;
    default:
    error:
#ifdef NANS
      einan (yy);
#else
      mtherr ("asctoe", DOMAIN);
      eclear (yy);
#endif
      goto aexit;
    }
 donchr:
  ++s;
  goto nxtcom;

  /* Exponent interpretation */
 expnt:

  esign = 1;
  exp = 0;
  ++s;
  /* check for + or - */
  if (*s == '-')
    {
      esign = -1;
      ++s;
    }
  if (*s == '+')
    ++s;
  while ((*s >= '0') && (*s <= '9'))
    {
      exp *= 10;
      exp += *s++ - '0';
      if (exp > -(MINDECEXP))
        {
          if (esign < 0)
            goto zero;
          else
            goto infinite;
        }
    }
  if (esign < 0)
    exp = -exp;
  if (exp > MAXDECEXP)
    {
 infinite:
      ecleaz (yy);
      yy[E] = 0x7fff;           /* infinity */
      goto aexit;
    }
  if (exp < MINDECEXP)
    {
 zero:
      ecleaz (yy);
      goto aexit;
    }

 daldone:
  nexp = exp - nexp;
  /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
  while ((nexp > 0) && (yy[2] == 0))
    {
      emovz (yy, xt);
      eshup1 (xt);
      eshup1 (xt);
      eaddm (yy, xt);
      eshup1 (xt);
      if (xt[2] != 0)
        break;
      nexp -= 1;
      emovz (xt, yy);
    }
  if ((k = enormlz (yy)) > NBITS)
    {
      ecleaz (yy);
      goto aexit;
    }
  lexp = (EXONE - 1 + NBITS) - k;
  emdnorm (yy, lost, 0, lexp, 64);

  /* 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.  */

  lexp = yy[E];
  if (nexp == 0)
    {
      k = 0;
      goto expdon;
    }
  esign = 1;
  if (nexp < 0)
    {
      nexp = -nexp;
      esign = -1;
      if (nexp > 4096)
        {
          /* Punt.  Can't handle this without 2 divides. */
          emovi (etens[0], tt);
          lexp -= tt[E];
          k = edivm (tt, yy);
          lexp += EXONE;
          nexp -= 4096;
        }
    }
  p = &etens[NTEN][0];
  emov (eone, xt);
  exp = 1;
  do
    {
      if (exp & nexp)
        emul (p, xt, xt);
      p -= NE;
      exp = exp + exp;
    }
  while (exp <= MAXP);

  emovi (xt, tt);
  if (esign < 0)
    {
      lexp -= tt[E];
      k = edivm (tt, yy);
      lexp += EXONE;
    }
  else
    {
      lexp += tt[E];
      k = emulm (tt, yy);
      lexp -= EXONE - 1;
    }

 expdon:

  /* 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;
#endif
#ifdef DEC
  else if (oprec == 56)
    lexp -= EXONE - 0201;
#endif
  rndprc = oprec;
  emdnorm (yy, k, 0, lexp, 64);

 aexit:

  rndprc = rndsav;
  yy[0] = nsign;
  switch (oprec)
    {
#ifdef DEC
    case 56:
      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 24:
      toe24 (yy, y);
      break;
    case 64:
      toe64 (yy, y);
      break;
    case 113:
      toe113 (yy, y);
      break;
    case NBITS:
      emovo (yy, y);
      break;
    }
}



/* y = largest integer not greater than x (truncated toward minus infinity)  */

static unsigned EMUSHORT bmask[] =
{
  0xffff,
  0xfffe,
  0xfffc,
  0xfff8,
  0xfff0,
  0xffe0,
  0xffc0,
  0xff80,
  0xff00,
  0xfe00,
  0xfc00,
  0xf800,
  0xf000,
  0xe000,
  0xc000,
  0x8000,
  0x0000,
};

static void 
efloor (x, y)
     unsigned EMUSHORT x[], y[];
{
  register unsigned EMUSHORT *p;
  int e, expon, i;
  unsigned EMUSHORT f[NE];

  emov (x, f);                  /* leave in external format */
  expon = (int) f[NE - 1];
  e = (expon & 0x7fff) - (EXONE - 1);
  if (e <= 0)
    {
      eclear (y);
      goto isitneg;
    }
  /* number of bits to clear out */
  e = NBITS - e;
  emov (f, y);
  if (e <= 0)
    return;

  p = &y[0];
  while (e >= 16)
    {
      *p++ = 0;
      e -= 16;
    }
  /* clear the remaining bits */
  *p &= bmask[e];
  /* truncate negatives toward minus infinity */
 isitneg:

  if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
    {
      for (i = 0; i < NE - 1; i++)
        {
          if (f[i] != y[i])
            {
              esub (eone, y, y);
              break;
            }
        }
    }
}


/* 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.  */

static void 
efrexp (x, exp, s)
     unsigned EMUSHORT x[];
     int *exp;
     unsigned EMUSHORT s[];
{
  unsigned EMUSHORT xi[NI];
  EMULONG li;

  emovi (x, xi);
  li = (EMULONG) ((EMUSHORT) xi[1]);

  if (li == 0)
    {
      li -= enormlz (xi);
    }
  xi[1] = 0x3ffe;
  emovo (xi, s);
  *exp = (int) (li - 0x3ffe);
}



/* Return y = x * 2**pwr2.  */

static void 
eldexp (x, pwr2, y)
     unsigned EMUSHORT x[];
     int pwr2;
     unsigned EMUSHORT y[];
{
  unsigned EMUSHORT xi[NI];
  EMULONG li;
  int i;

  emovi (x, xi);
  li = xi[1];
  li += pwr2;
  i = 0;
  emdnorm (xi, i, i, li, 64);
  emovo (xi, y);
}


/* c = remainder after dividing b by a
   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];

#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);
      eclear (c);
      return;
    }
  emovi (a, den);
  emovi (b, num);
  eiremain (den, num);
  /* Sign of remainder = sign of quotient */
  if (a[0] == b[0])
    num[0] = 0;
  else
    num[0] = 0xffff;
  emovo (num, c);
}

static void 
eiremain (den, num)
     unsigned EMUSHORT den[], num[];
{
  EMULONG ld, ln;
  unsigned EMUSHORT j;

  ld = den[E];
  ld -= enormlz (den);
  ln = num[E];
  ln -= enormlz (num);
  ecleaz (equot);
  while (ln >= ld)
    {
      if (ecmpm (den, num) <= 0)
        {
          esubm (den, num);
          j = 1;
        }
      else
        {
          j = 0;
        }
      eshup1 (equot);
      equot[NI - 1] |= j;
      eshup1 (num);
      ln -= 1;
    }
  emdnorm (num, 0, 0, ln, 0);
}

/* 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
     INVALID           7       NaN - producing operation
     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. */

/* Note: 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 */
  "singularity",                /* et seq.      */
  "overflow",
  "underflow",
  "total loss of precision",
  "partial loss of precision",
  "invalid operation"
};

int merror = 0;
extern int merror;

static void 
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.

     Display error message defined by the code argument.  */

  if ((code <= 0) || (code >= NMSGS))
    code = 0;
  sprintf (errstr, " %s %s error", name, ermsg[code]);
  if (extra_warnings)
    warning (errstr);
  /* Set global error message word */
  merror = code + 1;
}

#ifdef DEC
/* Convert DEC double precision to e type.  */

static void 
dectoe (d, e)
     unsigned EMUSHORT *d;
     unsigned EMUSHORT *e;
{
  unsigned EMUSHORT y[NI];
  register unsigned EMUSHORT r, *p;

  ecleaz (y);                   /* start with a zero */
  p = y;                        /* point to our number */
  r = *d;                       /* get DEC exponent word */
  if (*d & (unsigned int) 0x8000)
    *p = 0xffff;                /* fill in our sign */
  ++p;                          /* bump pointer to our exponent word */
  r &= 0x7fff;                  /* strip the sign bit */
  if (r == 0)                   /* answer = 0 if high order DEC word = 0 */
    goto done;


  r >>= 7;                      /* shift exponent word down 7 bits */
  r += EXONE - 0201;            /* subtract DEC exponent offset */
  /* add our e type exponent offset */
  *p++ = r;                     /* to form our exponent */

  r = *d++;                     /* now do the high order mantissa */
  r &= 0177;                    /* strip off the DEC exponent and sign bits */
  r |= 0200;                    /* the DEC understood high order mantissa bit */
  *p++ = r;                     /* put result in our high guard word */

  *p++ = *d++;                  /* fill in the rest of our mantissa */
  *p++ = *d++;
  *p = *d;

  eshdn8 (y);                   /* shift our mantissa down 8 bits */
 done:
  emovo (y, e);
}



/*
;       convert e type to DEC double precision
;       double d;
;       EMUSHORT e[NE];
;       etodec (e, &d);
*/

static void 
etodec (x, d)
     unsigned EMUSHORT *x, *d;
{
  unsigned EMUSHORT xi[NI];
  EMULONG exp;
  int rndsav;

  emovi (x, xi);
  exp = (EMULONG) xi[E] - (EXONE - 0201);       /* adjust exponent for offsets */
/* round off to nearest or even */
  rndsav = rndprc;
  rndprc = 56;
  emdnorm (xi, 0, 0, exp, 64);
  rndprc = rndsav;
  todec (xi, d);
}

static void 
todec (x, y)
     unsigned EMUSHORT *x, *y;
{
  unsigned EMUSHORT i;
  unsigned EMUSHORT *p;

  p = x;
  *y = 0;
  if (*p++)
    *y = 0100000;
  i = *p++;
  if (i == 0)
    {
      *y++ = 0;
      *y++ = 0;
      *y++ = 0;
      *y++ = 0;
      return;
    }
  if (i > 0377)
    {
      *y++ |= 077777;
      *y++ = 0xffff;
      *y++ = 0xffff;
      *y++ = 0xffff;
#ifdef ERANGE
      errno = ERANGE;
#endif
      return;
    }
  i &= 0377;
  i <<= 7;
  eshup8 (x);
  x[M] &= 0177;
  i |= x[M];
  *y++ |= i;
  *y++ = x[M + 1];
  *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;

  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 MIEEE
unsigned EMUSHORT TFnan[8] =
 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
#endif
#ifdef IBMPC
unsigned EMUSHORT TFnan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
#endif
#endif

#ifdef XFMODE_NAN
XFMODE_NAN;
#else
#ifdef MIEEE
unsigned EMUSHORT XFnan[6] = {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
#endif
#ifdef IBMPC
unsigned EMUSHORT XFnan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
#endif
#endif

#ifdef DFMODE_NAN
DFMODE_NAN;
#else
#ifdef MIEEE
unsigned EMUSHORT DFnan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
#endif
#ifdef IBMPC
unsigned EMUSHORT DFnan[4] = {0, 0, 0, 0xfff8};
#endif
#endif

#ifdef SFMODE_NAN
SFMODE_NAN;
#else
#ifdef MIEEE
unsigned EMUSHORT SFnan[2] = {0x7fff, 0xffff};
#endif
#ifdef IBMPC
unsigned EMUSHORT SFnan[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;
      p = TFnan;
      break;
    case XFmode:
      n = 6;
      p = XFnan;
      break;
    case DFmode:
      n = 4;
      p = DFnan;
      break;
    case SFmode:
      n = 2;
      p = SFnan;
      break;
#endif
    default:
      abort ();
    }
#ifdef MIEEE
  *nan++ = (sign << 15) | *p++;
#endif
  while (--n != 0)
    *nan++ = *p++;
#ifndef MIEEE
  *nan = (sign << 15) | *p;
#endif
}

/* 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 FLOAT_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);
#endif
  /* 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 FLOAT_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
#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 FLOAT_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++;
#endif
  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++;
#endif
  /* 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];
#endif
      k -= j;
      do
        {
          eshup6 (xi);
#if WORDS_BIG_ENDIAN
          *i++ = xi[M];
#else
          *i-- = xi[M];
#endif
        }
      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;
#endif
    }
}


/* 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];
#endif
      k -= j;
      do
        {
          eshup6 (xi);
#if WORDS_BIG_ENDIAN
          *i++ = xi[M];
#else
          *i-- = xi[M];
#endif
        }
      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;
#endif
    }
  /* Negate if negative */
  if (xi[0])
    {
      carry = 0;
#if WORDS_BIG_ENDIAN
      isave += 3;
#endif
      for (k = 0; k < 4; k++)
        {
          acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
#if WORDS_BIG_ENDIAN
          *isave-- = acc;
#else
          *isave++ = acc;
#endif
          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 */
\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;
{

switch (mode)
  {
  case SFmode:
    return 24;

  case DFmode:
#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 XFmode:
    return 64;
  case TFmode:
    return 113;

  default:
    abort ();
  }
}