(store_expr): Call size_binop instead of doing the same thing with

[gcc.git] / gcc / real.c

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

   Copyright (C) 1993 Free Software Foundation, Inc.

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 "config.h"

#include "tree.h"

/* To enable support of XFmode extended real floating point, define
LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).

To support cross compilation between IEEE and VAX floating
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  */

/* Type of computer arithmetic.
 * Only one of DEC, MIEEE, IBMPC, or UNK should get defined.
 * The following modification converts gcc macros into the ones
 * used by ieee.c.
 *
 * Note: long double formats differ between IBMPC and MIEEE
 * by more than just endian-ness.
 */

/* 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 == IEEE_FLOAT_FORMAT
#if WORDS_BIG_ENDIAN
/* Motorola IEEE, high order words come first (Sun workstation): */
#define MIEEE 1
#else /* not big-endian */
/* Intel IEEE, low order words come first:
 */
#define IBMPC 1
#endif /*  big-endian */
#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 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 == IEEE_FLOAT_FORMAT
#ifdef HOST_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 VAX */

#endif /* REAL_ARITHMETIC not defined */

/* Define for support of infinity.  */
#ifndef DEC
#define INFINITY
#endif


/* ehead.h
 *
 * Include file for extended precision arithmetic programs.
 */

/* Number of 16 bit words in external e type format */
#define NE 6

/* Number of 16 bit words in internal format */
#define NI (NE+3)

/* Array offset to exponent */
#define E 1

/* Array offset to high guard word */
#define M 2

/* Number of bits of precision */
#define NBITS ((NI-4)*16)

/* Maximum number of decimal digits in ASCII conversion
 * = NBITS*log10(2)
 */
#define NDEC (NBITS*8/27)

/* The exponent of 1.0 */
#define EXONE (0x3fff)

/* 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
#define GET_REAL(r,e) bcopy (r, e, 2*NE)
#define PUT_REAL(e,r) bcopy (e, r, 2*NE)
#else /* no XFmode */

#ifdef REAL_ARITHMETIC
/* Emulator uses target format internally
   but host stores it in host endian-ness. */

#if defined (HOST_WORDS_BIG_ENDIAN) == WORDS_BIG_ENDIAN
#define GET_REAL(r,e) e53toe ((r), (e))
#define PUT_REAL(e,r) etoe53 ((e), (r))

#else /* endian-ness differs */
/* emulator uses target endian-ness internally */
#define GET_REAL(r,e)		\
do { EMUSHORT w[4];		\
 w[3] = ((EMUSHORT *) r)[0];	\
 w[2] = ((EMUSHORT *) r)[1];	\
 w[1] = ((EMUSHORT *) r)[2];	\
 w[0] = ((EMUSHORT *) r)[3];	\
 e53toe (w, (e)); } while (0)

#define PUT_REAL(e,r)		\
do { EMUSHORT w[4];		\
 etoe53 ((e), w);		\
 *((EMUSHORT *) r) = w[3];	\
 *((EMUSHORT *) r + 1) = w[2];	\
 *((EMUSHORT *) r + 2) = w[1];	\
 *((EMUSHORT *) r + 3) = w[0]; } while (0)

#endif /* endian-ness differs */

#else /* not REAL_ARITHMETIC */

/* emulator uses host format */
#define GET_REAL(r,e) e53toe ((r), (e))
#define PUT_REAL(e,r) etoe53 ((e), (r))

#endif /* not REAL_ARITHMETIC */
#endif /* no XFmode */

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

/* Pack output array with 32-bit numbers obtained from
   array containing 16-bit numbers, swapping ends if required. */
void 
endian (e, x, mode)
     unsigned EMUSHORT e[];
     long x[];
     enum machine_mode mode;
{
  unsigned long th, t;

#if WORDS_BIG_ENDIAN
  switch (mode)
    {

    case XFmode:

      /* Swap halfwords in the third long. */
      th = (unsigned long) e[4] & 0xffff;
      t = (unsigned long) e[5] & 0xffff;
      t |= th << 16;
      x[2] = (long) t;
      /* fall into the double case */

    case DFmode:

      /* swap halfwords in the second word */
      th = (unsigned long) e[2] & 0xffff;
      t = (unsigned long) e[3] & 0xffff;
      t |= th << 16;
      x[1] = (long) t;
      /* fall into the float case */

    case SFmode:

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

    default:
      abort ();
    }

#else

  /* Pack the output array without swapping. */

  switch (mode)
    {

    case XFmode:

      /* Pack the third long.
	 Each element of the input REAL_VALUE_TYPE array has 16 bit useful bits
	 in it.  */
      th = (unsigned long) e[5] & 0xffff;
      t = (unsigned long) e[4] & 0xffff;
      t |= th << 16;
      x[2] = (long) t;
      /* fall into the double case */

    case DFmode:

      /* pack the second long */
      th = (unsigned long) e[3] & 0xffff;
      t = (unsigned long) e[2] & 0xffff;
      t |= th << 16;
      x[1] = (long) t;
      /* fall into the float case */

    case SFmode:

      /* 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);
  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)
	abort ();
#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_FIX_TRUNCATE (x) (etrunci (x))
 */
REAL_VALUE_TYPE 
etrunci (x)
     REAL_VALUE_TYPE x;
{
  unsigned EMUSHORT f[NE], g[NE];
  REAL_VALUE_TYPE r;
  long l;

  GET_REAL (&x, g);
  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_FIX_TRUNCATE (x) (etruncui (x))
 */
REAL_VALUE_TYPE 
etruncui (x)
     REAL_VALUE_TYPE x;
{
  unsigned EMUSHORT f[NE], g[NE];
  REAL_VALUE_TYPE r;
  unsigned long l;

  GET_REAL (&x, g);
  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;
    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 to int
 * implements REAL_VALUE_FIX (x) (eroundi (x))
 * The type of rounding is left unspecified by real.h.
 * It is implemented here as round to nearest (add .5 and chop).
 */
int 
eroundi (x)
     REAL_VALUE_TYPE x;
{
  unsigned EMUSHORT f[NE], g[NE];
  EMULONG l;

  GET_REAL (&x, f);
  eround (f, g);
  eifrac (g, &l, f);
  return ((int) l);
}

/* Round real to nearest unsigned int
 * implements  REAL_VALUE_UNSIGNED_FIX (x) ((unsigned int) eroundi (x))
 * Negative input returns zero.
 * The type of rounding is left unspecified by real.h.
 * It is implemented here as round to nearest (add .5 and chop).
 */
unsigned int 
eroundui (x)
     REAL_VALUE_TYPE x;
{
  unsigned EMUSHORT f[NE], g[NE];
  unsigned EMULONG l;

  GET_REAL (&x, f);
  eround (f, g);
  euifrac (g, &l, f);
  return ((unsigned int)l);
}


/* REAL_VALUE_FROM_INT macro.
 */
void 
ereal_from_int (d, i, j)
     REAL_VALUE_TYPE *d;
     long i, j;
{
  unsigned EMUSHORT df[NE], dg[NE];
  long 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_LONG, df);
  ultoe (&high, dg);
  emul (dg, df, dg);
  ultoe (&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 long i, j;
{
  unsigned EMUSHORT df[NE], dg[NE];
  unsigned long low, high;

  low = i;
  high = j;
  eldexp (eone, HOST_BITS_PER_LONG, 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)
     long *low, *high;
     REAL_VALUE_TYPE rr;
{
  unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
  int s;

  GET_REAL (&rr, d);
  /* convert positive value */
  s = 0;
  if (eisneg (d))
    {
      eneg (d);
      s = 1;
    }
  eldexp (eone, HOST_BITS_PER_LONG, df);
  ediv (df, d, dg);		/* dg = d / 2^32 is the high word */
  euifrac (dg, high, dh);
  emul (df, dh, dg);		/* fractional part is the low word */
  euifrac (dg, 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);
  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 an IEEE double precision number is a NaN. */

int
target_isnan (x)
     REAL_VALUE_TYPE x;
{
  return (0);
}


/* Check for a negative IEEE double precision number.
 * this means strictly less than zero, not -0.
 */

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

  GET_REAL (&x, e);
  if (ecmp (e, ezero) < 0)
    return (1);
  return (0);
}

/* 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);
  eclear (t);
  switch (mode)
    {
    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 (e);
      return (r);

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

#endif /* REAL_ARITHMETIC defined */

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

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

/*
 * Revision history:
 *
 *  5 Jan 84	PDP-11 assembly language version
 *  2 Mar 86	fixed bug in asctoq
 *  6 Dec 86	C language version
 * 30 Aug 88	100 digit version, improved rounding
 * 15 May 92    80-bit long double support
 *
 * Author:  S. L. Moshier.
 */


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

#define DOMAIN		1	/* argument domain error */
#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 EDOM		33
#define ERANGE		34

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

/*include "ehead.h"*/
/* 0.0 */
unsigned EMUSHORT ezero[NE] =
{
  0, 0000000, 0000000, 0000000, 0000000, 0000000,};
extern unsigned EMUSHORT ezero[];

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

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

/* 2.0E0 */
unsigned EMUSHORT etwo[NE] =
{
  0, 0000000, 0000000, 0000000, 0100000, 0040000,};
extern unsigned EMUSHORT etwo[];

/* 3.2E1 */
unsigned EMUSHORT e32[NE] =
{
  0, 0000000, 0000000, 0000000, 0100000, 0040004,};
extern unsigned EMUSHORT e32[];

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

/* 1.41421356237309504880168872420969807856967187537695E0 */
unsigned EMUSHORT esqrt2[NE] =
{
  0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
extern unsigned EMUSHORT esqrt2[];

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

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

/* 5.7721566490153286060651209008240243104215933593992E-1 */
unsigned EMUSHORT eeul[NE] =
{
  0xd1be, 0xc7a4, 0076660, 0063743, 0111704, 0x3ffe,};
extern unsigned EMUSHORT eeul[];

/*
include "ehead.h"
include "mconf.h"
*/



/* Control register for rounding precision.
 * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
 */
int rndprc = NBITS;
extern int rndprc;

void eaddm (), esubm (), emdnorm (), asctoeg ();
static void toe24 (), toe53 (), toe64 ();
void eremain (), einit (), eiremain ();
int ecmpm (), edivm (), emulm ();
void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
void etodec (), todec (), dectoe ();




void 
einit ()
{
}

/*
; Clear out entire external format number.
;
; unsigned EMUSHORT x[];
; eclear (x);
*/

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.
 *
 * emov (a, b);
 */

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
;
;	EMUSHORT x[NE];
;	eabs (x);
*/

void 
eabs (x)
     unsigned EMUSHORT x[];	/* x is the memory address of a short */
{

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




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

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

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



/* Return 1 if external format number is negative,
 * else return zero.
 */
int 
eisneg (x)
     unsigned EMUSHORT x[];
{

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


/* Return 1 if external format number has maximum possible exponent,
 * else return zero.
 */
int 
eisinf (x)
     unsigned EMUSHORT x[];
{

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


/*
; Fill entire number, including exponent and significand, with
; largest possible number.  These programs implement a saturation
; value that is an ordinary, legal number.  A special value
; "infinity" may also be implemented; this would require tests
; for that value and implementation of special rules for arithmetic
; operations involving inifinity.
*/

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 == 64)
	{
	  *(x - 5) = 0;
	}
      if (rndprc == 53)
	{
	  *(x - 4) = 0xf800;
	}
      else
	{
	  *(x - 4) = 0;
	  *(x - 3) = 0;
	  *(x - 2) = 0xff00;
	}
    }
#endif
}



/* Move in external format number,
 * converting it to internal format.
 */
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)
    {
      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.
 */
void 
emovo (a, b)
     unsigned EMUSHORT *a, *b;
{
  register unsigned EMUSHORT *p, *q;
  unsigned EMUSHORT i;

  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)
    {
      einfin (b);
      return;
    }
#endif
  /* skip over guard word */
  ++p;
  /* move the significand */
  for (i = 0; i < NE - 1; i++)
    *q-- = *p++;
}




/* Clear out internal format number.
 */

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

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


/*
;	Compare significands of numbers in internal format.
;	Guard words are included in the comparison.
;
;	unsigned EMUSHORT a[NI], b[NI];
;	cmpm (a, b);
;
;	for the significands:
;	returns	+1 if a > b
;		 0 if a == b
;		-1 if a < b
*/
int
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
*/

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

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

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

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

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

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

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

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


/* Divide significands */

static unsigned EMUSHORT equot[NI];

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



/*
 * Normalize and round off.
 *
 * The internal format number to be rounded is "s".
 * Input "lost" indicates whether or not the number is exact.
 * This is the so-called sticky bit.
 *
 * Input "subflg" indicates whether the number was obtained
 * by a subtraction operation.  In that case if lost is nonzero
 * then the number is slightly smaller than indicated.
 *
 * Input "exp" is the biased exponent, which may be negative.
 * the exponent field of "s" is ignored but is replaced by
 * "exp" as adjusted by normalization and rounding.
 *
 * Input "rcntrl" is the rounding control.
 */

static int rlast = -1;
static int rw = 0;
static unsigned EMUSHORT rmsk = 0;
static unsigned EMUSHORT rmbit = 0;
static unsigned EMUSHORT rebit = 0;
static int re = 0;
static unsigned EMUSHORT rbit[NI];

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;
	  rbit[rw - 1] = 1;
	  re = NI - 2;
	  rebit = 1;
	  break;
	case 64:
	  rw = 7;
	  rmsk = 0xffff;
	  rmbit = 0x8000;
	  rbit[rw - 1] = 1;
	  re = rw - 1;
	  rebit = 1;
	  break;
	  /* For DEC arithmetic */
	case 56:
	  rw = 6;
	  rmsk = 0xff;
	  rmbit = 0x80;
	  rbit[rw] = 0x100;
	  re = rw;
	  rebit = 0x100;
	  break;
	case 53:
	  rw = 6;
	  rmsk = 0x7ff;
	  rmbit = 0x0400;
	  rbit[rw] = 0x800;
	  re = rw;
	  rebit = 0x800;
	  break;
	case 24:
	  rw = 4;
	  rmsk = 0xff;
	  rmbit = 0x80;
	  rbit[rw] = 0x100;
	  re = rw;
	  rebit = 0x100;
	  break;
	}
      rlast = rndprc;
    }

  if (rndprc >= 64)
    {
      r = s[rw] & rmsk;
      if (rndprc == 64)
	{
	  i = rw + 1;
	  while (i < NI)
	    {
	      if (s[i])
		r |= 1;
	      s[i] = 0;
	      ++i;
	    }
	}
    }
  else
    {
      if (exp <= 0)
	eshdn1 (s);
      r = s[rw] & rmsk;
      /* These tests assume NI = 8 */
      i = rw + 1;
      while (i < NI)
	{
	  if (s[i])
	    r |= 1;
	  s[i] = 0;
	  ++i;
	}
      /*
	 if (rndprc == 24)
	 {
	 if (s[5] || s[6])
	 r |= 1;
	 s[5] = 0;
	 s[6] = 0;
	 }
	 */
      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 ((rndprc < 64) && (exp <= 0))
    {
      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;
      pedwarn ("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)
	{
	  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.
;
;	unsigned EMUSHORT a[NE], b[NE], c[NE];
;	esub (a, b, c);	 c = b - a
*/

static int subflg = 0;

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

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


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

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

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.
;
;	unsigned EMUSHORT a[NE], b[NE], c[NE];
;	ediv (a, b, c);	c = b / a
*/
void 
ediv (a, b, c)
     unsigned EMUSHORT *a, *b, *c;
{
  unsigned EMUSHORT ai[NI], bi[NI];
  int i;
  EMULONG lt, lta, ltb;

#ifdef INFINITY
  if (eisinf (b))
    {
      if (eisneg (a) ^ eisneg (b))
	*(c + (NE - 1)) = 0x8000;
      else
	*(c + (NE - 1)) = 0;
      einfin (c);
      return;
    }
  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;
      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.
;
;	unsigned EMUSHORT a[NE], b[NE], c[NE];
;	emul (a, b, c);	c = b * a
*/
void 
emul (a, b, c)
     unsigned EMUSHORT *a, *b, *c;
{
  unsigned EMUSHORT ai[NI], bi[NI];
  int i, j;
  EMULONG lt, lta, ltb;

#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
;	double d;
;	unsigned EMUSHORT x[N+2];
;	e53toe (&d, x);
*/
void 
e53toe (e, y)
     unsigned EMUSHORT *e, *y;
{
#ifdef DEC

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

#else

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

  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)
    {
      einfin (y);
      if (r & 0x8000)
	eneg (y);
      return;
    }
#endif
  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
  (void) 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 DEC */
}

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

  p = yy;
  for (i = 0; i < NE - 5; i++)
    *p++ = 0;
#ifdef IBMPC
  for (i = 0; i < 5; i++)
    *p++ = *e++;
#endif
#ifdef DEC
  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)
    {
      einfin (y);
      if (*p & 0x8000)
	eneg (y);
      return;
    }
#endif
  for (i = 0; i < NE; i++)
    *q++ = *p++;
}


/*
; Convert IEEE single precision to e type
;	float d;
;	unsigned EMUSHORT x[N+2];
;	dtox (&d, x);
*/
void 
e24toe (e, y)
     unsigned EMUSHORT *e, *y;
{
  register unsigned EMUSHORT r;
  register unsigned EMUSHORT *p;
  unsigned EMUSHORT yy[NI];
  int denorm, k;

  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)
    {
      einfin (y);
      if (r & 0x8000)
	eneg (y);
      return;
    }
#endif
  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
  (void) 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);
}


void 
etoe64 (x, e)
     unsigned EMUSHORT *x, *e;
{
  unsigned EMUSHORT xi[NI];
  EMULONG exp;
  int rndsav;

  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;

  p = a;
#ifdef MIEEE
  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++;
#ifdef MIEEE
  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 */
#ifdef MIEEE
  for (i = 0; i < 4; i++)
    *q++ = *p++;
#else
  for (i = 0; i < 4; i++)
    *q-- = *p++;
#endif
}


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

#ifdef DEC

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

void 
etoe53 (x, e)
     unsigned EMUSHORT *x, *e;
{
  unsigned EMUSHORT xi[NI];
  EMULONG exp;
  int rndsav;

  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;


  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)
    {
      (void) eshift (x, 4);
    }
  else
    {
      i <<= 4;
      (void) 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 DEC */



/*
; e type to IEEE single precision
;	float d;
;	unsigned EMUSHORT x[N+2];
;	xtod (x, &d);
*/
void 
etoe24 (x, e)
     unsigned EMUSHORT *x, *e;
{
  EMULONG exp;
  unsigned EMUSHORT xi[NI];
  int rndsav;

  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;

  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++;
  if (i >= 255)
    {				/* Saturate at largest number less than infinity. */
#ifdef INFINITY
      *y |= (unsigned EMUSHORT) 0x7f80;
#ifdef IBMPC
      *(--y) = 0;
#endif
#ifdef DEC
      *(--y) = 0;
#endif
#ifdef MIEEE
      ++y;
      *y = 0;
#endif
#else
      *y |= (unsigned EMUSHORT) 0x7f7f;
#ifdef IBMPC
      *(--y) = 0xffff;
#endif
#ifdef DEC
      *(--y) = 0xffff;
#endif
#ifdef MIEEE
      ++y;
      *y = 0xffff;
#endif
#endif
      return;
    }
  if (i == 0)
    {
      (void) eshift (x, 7);
    }
  else
    {
      i <<= 7;
      (void) 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
}


/* Compare two e type numbers.
 *
 * unsigned EMUSHORT a[NE], b[NE];
 * ecmp (a, b);
 *
 *  returns +1 if a > b
 *           0 if a == b
 *          -1 if a < b
 */
int 
ecmp (a, b)
     unsigned EMUSHORT *a, *b;
{
  unsigned EMUSHORT ai[NI], bi[NI];
  register unsigned EMUSHORT *p, *q;
  register int i;
  int msign;

  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)
 *
 * unsigned EMUSHORT x[NE], y[NE]
 * eround (x, y);
 */
void 
eround (x, y)
     unsigned EMUSHORT *x, *y;
{
  eadd (ehalf, x, y);
  efloor (y, y);
}




/*
; convert long integer to e type
;
;	long l;
;	unsigned EMUSHORT x[NE];
;	ltoe (&l, x);
; note &l is the memory address of l
*/
void 
ltoe (lp, y)
     long *lp;			/* lp is the memory address of a long integer */
     unsigned EMUSHORT *y;		/* y is the address of a short */
{
  unsigned EMUSHORT yi[NI];
  unsigned long ll;
  int k;

  ecleaz (yi);
  if (*lp < 0)
    {
      /* make it positive */
      ll = (unsigned long) (-(*lp));
      yi[0] = 0xffff;		/* put correct sign in the e type number */
    }
  else
    {
      ll = (unsigned long) (*lp);
    }
  /* move the long integer to yi significand area */
  yi[M] = (unsigned EMUSHORT) (ll >> 16);
  yi[M + 1] = (unsigned EMUSHORT) ll;

  yi[E] = EXONE + 15;		/* 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, y);		/* output the answer */
}

/*
; convert unsigned long integer to e type
;
;	unsigned long l;
;	unsigned EMUSHORT x[NE];
;	ltox (&l, x);
; note &l is the memory address of l
*/
void 
ultoe (lp, y)
     unsigned long *lp;		/* lp is the memory address of a long integer */
     unsigned EMUSHORT *y;		/* y is the address of a short */
{
  unsigned EMUSHORT yi[NI];
  unsigned long ll;
  int k;

  ecleaz (yi);
  ll = *lp;

  /* move the long integer to ayi significand area */
  yi[M] = (unsigned EMUSHORT) (ll >> 16);
  yi[M + 1] = (unsigned EMUSHORT) ll;

  yi[E] = EXONE + 15;		/* 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, y);		/* output the answer */
}


/*
;	Find long integer and fractional parts

;	long i;
;	unsigned EMUSHORT x[NE], frac[NE];
;	xifrac (x, &i, frac);

  The integer output has the sign of the input.  The fraction is
the positive fractional part of abs (x).
*/
void 
eifrac (x, i, frac)
     unsigned EMUSHORT *x;
     long *i;
     unsigned EMUSHORT *frac;
{
  unsigned EMUSHORT xi[NI];
  int k;

  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_LONG - 1))
    {
      /*
	 ;	long integer overflow: output large integer
	 ;	and correct fraction
	 */
      if (xi[0])
	*i = ((unsigned long) 1) << (HOST_BITS_PER_LONG - 1);
      else
	*i = (((unsigned long) 1) << (HOST_BITS_PER_LONG - 1)) - 1;
      (void) eshift (xi, k);
      pedwarn ("Overflow on truncation to integer");
      goto lab11;
    }

  if (k > 16)
    {
      /*
	 ; shift more than 16 bits: shift up k-16, output the integer,
	 ; then complete the shift to get the fraction.
	 */
      k -= 16;
      (void) eshift (xi, k);

      *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
      eshup6 (xi);
      goto lab10;
    }

  /* shift not more than 16 bits */
  (void) eshift (xi, k);
  *i = (long) xi[M] & 0xffff;

 lab10:

  if (xi[0])
    *i = -(*i);
 lab11:

  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 long integer and fractional parts

;	unsigned long i;
;	unsigned EMUSHORT x[NE], frac[NE];
;	xifrac (x, &i, frac);

  A negative e type input yields integer output = 0
  but correct fraction.
*/
void 
euifrac (x, i, frac)
     unsigned EMUSHORT *x;
     long *i;
     unsigned EMUSHORT *frac;
{
  unsigned EMUSHORT xi[NI];
  int 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 > 32)
    {
      /*
	 ;	long integer overflow: output large integer
	 ;	and correct fraction
	 */
      *i = ~(0L);
      (void) eshift (xi, k);
      pedwarn ("Overflow on truncation to unsigned integer");
      goto lab10;
    }

  if (k > 16)
    {
      /*
	 ; shift more than 16 bits: shift up k-16, output the integer,
	 ; then complete the shift to get the fraction.
	 */
      k -= 16;
      (void) eshift (xi, k);

      *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
      eshup6 (xi);
      goto lab10;
    }

  /* shift not more than 16 bits */
  (void) eshift (xi, k);
  *i = (long) xi[M] & 0xffff;

 lab10:

  if (xi[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
;
;	Shifts significand area up or down by the number of bits
;	given by the variable sc.
*/
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);
}



/*
;	normalize
;
; Shift normalizes the significand area pointed to by argument
; shift count (up = positive) is returned.
*/
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

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 */
};

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

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


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

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


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

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


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

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;

  ss = string;
  s = wstring;
  while ((*s++ = *ss++) != '\0')
    ;
  rndsav = rndprc;
  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.  Don't bother with illegal infinities.
   */
  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 < 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];
      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++ = ' ';
  *s++ = (char) digit + '0';
  *s++ = '.';
  /* Examine number of digits requested by caller. */
  if (ndigs < 0)
    ndigs = 0;
  if (ndigs > NDEC)
    ndigs = NDEC;
  /* 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')
    ;
}




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

/* ASCII to single */
void 
asctoe24 (s, y)
     char *s;
     unsigned EMUSHORT *y;
{
  asctoeg (s, y, 24);
}


/* ASCII to double */
void 
asctoe53 (s, y)
     char *s;
     unsigned EMUSHORT *y;
{
#ifdef DEC
  asctoeg (s, y, 56);
#else
  asctoeg (s, y, 53);
#endif
}


/* ASCII to long double */
void 
asctoe64 (s, y)
     char *s;
     unsigned EMUSHORT *y;
{
  asctoeg (s, y, 64);
}

/* ASCII to super double */
void 
asctoe (s, y)
     char *s;
     unsigned EMUSHORT *y;
{
  asctoeg (s, y, NBITS);
}

/* Space to make a copy of the input string: */
static char lstr[82];

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;

  /* Copy the input string. */
  s = ss;
  while (*s == ' ')		/* skip leading spaces */
    ++s;
  sp = lstr;
  for (k = 0; k < 79; k++)
    {
      if ((*sp++ = *s++) == '\0')
	break;
    }
  *sp = '\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
	{
	  lost |= k;
	}
      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':
      ecleaz (yy);
      yy[E] = 0x7fff;		/* infinity */
      goto aexit;
    default:
    error:
      mtherr ("asctoe", DOMAIN);
      eclear (y);
      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 (esign < 0)
    exp = -exp;

 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;
  else if (oprec == 24)
    lexp -= EXONE - 0177;
#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
    case 53:
      toe53 (yy, y);
      break;
    case 24:
      toe24 (yy, y);
      break;
    case 64:
      toe64 (yy, y);
      break;
    case NBITS:
      emovo (yy, y);
      break;
    }
}



/* y = largest integer not greater than x
 * (truncated toward minus infinity)
 *
 * unsigned EMUSHORT x[NE], y[NE]
 *
 * efloor (x, y);
 */
static unsigned EMUSHORT bmask[] =
{
  0xffff,
  0xfffe,
  0xfffc,
  0xfff8,
  0xfff0,
  0xffe0,
  0xffc0,
  0xff80,
  0xff00,
  0xfe00,
  0xfc00,
  0xf800,
  0xf000,
  0xe000,
  0xc000,
  0x8000,
  0x0000,
};

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


/* unsigned EMUSHORT x[], s[];
 * int *exp;
 *
 * efrexp (x, exp, s);
 *
 * Returns s and exp such that  s * 2**exp = x and .5 <= s < 1.
 * For example, 1.1 = 0.55 * 2**1
 * Handles denormalized numbers properly using long integer exp.
 */
void 
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);
}



/* unsigned EMUSHORT x[], y[];
 * long pwr2;
 *
 * eldexp (x, pwr2, y);
 *
 * Returns y = x * 2**pwr2.
 */
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[].
 */
void 
eremain (a, b, c)
     unsigned EMUSHORT a[], b[], c[];
{
  unsigned EMUSHORT den[NI], num[NI];

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

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

/*							mtherr.c
 *
 *	Library common error handling routine
 *
 *
 *
 * SYNOPSIS:
 *
 * char *fctnam;
 * int code;
 * void mtherr ();
 *
 * mtherr (fctnam, code);
 *
 *
 *
 * DESCRIPTION:
 *
 * This routine may be called to report one of the following
 * error conditions (in the include file mconf.h).
 *
 *   Mnemonic        Value          Significance
 *
 *    DOMAIN            1       argument domain error
 *    SING              2       function singularity
 *    OVERFLOW          3       overflow range error
 *    UNDERFLOW         4       underflow range error
 *    TLOSS             5       total loss of precision
 *    PLOSS             6       partial loss of precision
 *    EDOM             33       Unix domain error code
 *    ERANGE           34       Unix range error code
 *
 * The default version of the file prints the function name,
 * passed to it by the pointer fctnam, followed by the
 * error condition.  The display is directed to the standard
 * output device.  The routine then returns to the calling
 * program.  Users may wish to modify the program to abort by
 * calling exit under severe error conditions such as domain
 * errors.
 *
 * Since all error conditions pass control to this function,
 * the display may be easily changed, eliminated, or directed
 * to an error logging device.
 *
 * SEE ALSO:
 *
 * mconf.h
 *
 */
\f
/*
Cephes Math Library Release 2.0:  April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

/* include "mconf.h" */

/* Notice: the order of appearance of the following
 * messages is bound to the error codes defined
 * in mconf.h.
 */
static char *ermsg[7] =
{
  "unknown",			/* error code 0 */
  "domain",			/* error code 1 */
  "singularity",		/* et seq.      */
  "overflow",
  "underflow",
  "total loss of precision",
  "partial loss of precision"
};

int merror = 0;
extern int merror;

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 >= 6))
    code = 0;
  sprintf (errstr, "\n%s %s error\n", name, ermsg[code]);
  pedwarn (errstr);
  /* Set global error message word */
  merror = code + 1;

  /* Return to calling
   * program
   */
}

/* Here is etodec.c .
 *
 */

/*
;	convert DEC double precision to e type
;	double d;
;	EMUSHORT e[NE];
;	dectoe (&d, e);
*/
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);
*/
#if 0
static unsigned EMUSHORT decbit[NI] = {0, 0, 0, 0, 0, 0, 0200, 0};

void 
etodec (x, d)
     unsigned EMUSHORT *x, *d;
{
  unsigned EMUSHORT xi[NI];
  register unsigned EMUSHORT r;
  int i, j;

  emovi (x, xi);
  *d = 0;
  if (xi[0] != 0)
    *d = 0100000;
  r = xi[E];
  if (r < (EXONE - 128))
    goto zout;
  i = xi[M + 4];
  if ((i & 0200) != 0)
    {
      if ((i & 0377) == 0200)
	{
	  if ((i & 0400) != 0)
	    {
	      /* check all less significant bits */
	      for (j = M + 5; j < NI; j++)
		{
		  if (xi[j] != 0)
		    goto yesrnd;
		}
	    }
	  goto nornd;
	}
    yesrnd:
      eaddm (decbit, xi);
      r -= enormlz (xi);
    }

 nornd:

  r -= EXONE;
  r += 0201;
  if (r < 0)
    {
    zout:
      *d++ = 0;
      *d++ = 0;
      *d++ = 0;
      *d++ = 0;
      return;
    }
  if (r >= 0377)
    {
      *d++ = 077777;
      *d++ = -1;
      *d++ = -1;
      *d++ = -1;
      return;
    }
  r &= 0377;
  r <<= 7;
  eshup8 (xi);
  xi[M] &= 0177;
  r |= xi[M];
  *d++ |= r;
  *d++ = xi[M + 1];
  *d++ = xi[M + 2];
  *d++ = xi[M + 3];
}

#else

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

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;
      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 /* not 0 */

#endif /* EMU_NON_COMPILE not defined */