Google Git
Sign in
gnu / gcc / 737cee66f221eb8de21f0e41affc4b13cb7c0feb / . / gcc / real.c
blob: f7e22eae3bbee36c817fc9da86ffed201c76c1fc [file] [log] [blame]
/* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
and support for XFmode IEEE extended real floating point arithmetic.
Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
Contributed by Stephen L. Moshier (moshier@world.std.com).
This file is part of GNU CC.
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, 59 Temple Place - Suite 330,
Boston, MA 02111-1307, 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 a floating point
arithmetic suite that was not written with gcc in mind. Avoid
changing the low-level arithmetic routines unless you have suitable
test programs available. A special version of the PARANOIA floating
point arithmetic tester, modified for this purpose, can be found on
usc.edu: /pub/C-numanal/ieeetest.zoo. Other tests, and libraries of
XFmode and TFmode transcendental functions, can be obtained by ftp from
netlib.att.com: netlib/cephes. */
/* Type of computer arithmetic.
Only one of DEC, IBM, IEEE, or UNK should get defined.
`IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
to big-endian IEEE floating-point data structure. This definition
should work in SFmode `float' type and DFmode `double' type on
virtually all big-endian IEEE machines. If LONG_DOUBLE_TYPE_SIZE
has been defined to be 96, then IEEE also invokes the particular
XFmode (`long double' type) data structure used by the Motorola
680x0 series processors.
`IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
has been defined to be 96, then IEEE also invokes the particular
XFmode `long double' data structure used by the Intel 80x86 series
processors.
`DEC' refers specifically to the Digital Equipment Corp PDP-11
and VAX floating point data structure. This model currently
supports no type wider than DFmode.
`IBM' refers specifically to the IBM System/370 and compatible
floating point data structure. This model currently supports
no type wider than DFmode. The IBM conversions were contributed by
frank@atom.ansto.gov.au (Frank Crawford).
If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
then `long double' and `double' are both implemented, but they
both mean DFmode. In this case, the software floating-point
support available here is activated by writing
#define REAL_ARITHMETIC
in tm.h.
The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
and may deactivate XFmode since `long double' is used to refer
to both modes.
The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
separate the floating point unit's endian-ness from that of
the integer addressing. This permits one to define a big-endian
FPU on a little-endian machine (e.g., ARM). An extension to
BYTES_BIG_ENDIAN may be required for some machines in the future.
These optional macros may be defined in tm.h. In real.h, they
default to WORDS_BIG_ENDIAN, etc., so there is no need to define
them for any normal host or target machine on which the floats
and the integers have the same endian-ness. */
/* The following converts gcc macros into the ones used by this file. */
/* REAL_ARITHMETIC defined means that macros in real.h are
defined to call emulator functions. */
#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
#define IEEE
#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 */
#define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
#else
/* REAL_ARITHMETIC not defined means that the *host's* data
structure will be used. It may differ by endian-ness from the
target machine's structure and will get its ends swapped
accordingly (but not here). Probably only the decimal <-> binary
functions in this file will actually be used in this case. */
#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
#define IEEE
#else /* it's not IEEE either */
unknown arithmetic type
#define UNK 1
#endif /* not IEEE */
#endif /* not IBM */
#endif /* not VAX */
#define REAL_WORDS_BIG_ENDIAN HOST_FLOAT_WORDS_BIG_ENDIAN
#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
/* 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 ((char *) r, (char *) e, 2*NE)
#define PUT_REAL(e,r) bcopy ((char *) e, (char *) 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 ((char *) r, (char *) e, 2*NE)
#define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
#else
#define NE 6
#define MAXDECEXP 4932
#define MINDECEXP -4956
#ifdef REAL_ARITHMETIC
/* Emulator uses target format internally
but host stores it in host endian-ness. */
#define GET_REAL(r,e) \
do { \
if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
e53toe ((unsigned EMUSHORT *) (r), (e)); \
else \
{ \
unsigned EMUSHORT w[4]; \
w[3] = ((EMUSHORT *) r)[0]; \
w[2] = ((EMUSHORT *) r)[1]; \
w[1] = ((EMUSHORT *) r)[2]; \
w[0] = ((EMUSHORT *) r)[3]; \
e53toe (w, (e)); \
} \
} while (0)
#define PUT_REAL(e,r) \
do { \
if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
etoe53 ((e), (unsigned EMUSHORT *) (r)); \
else \
{ \
unsigned EMUSHORT w[4]; \
etoe53 ((e), w); \
*((EMUSHORT *) r) = w[3]; \
*((EMUSHORT *) r + 1) = w[2]; \
*((EMUSHORT *) r + 2) = w[1]; \
*((EMUSHORT *) r + 3) = w[0]; \
} \
} while (0)
#else /* not REAL_ARITHMETIC */
/* emulator uses host format */
#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 *));
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 *));
/* 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 (REAL_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 HFmode:
case SFmode:
/* swap halfwords in the first word */
th = (unsigned long) e[0] & 0xffff;
t = (unsigned long) e[1] & 0xffff;
t |= th << 16;
x[0] = (long) t;
break;
default:
abort ();
}
}
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 HFmode:
case SFmode:
/* pack the first long */
th = (unsigned long) e[1] & 0xffff;
t = (unsigned long) e[0] & 0xffff;
t |= th << 16;
x[0] = (long) t;
break;
default:
abort ();
}
}
}
/* 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 HFmode:
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, mode)
REAL_VALUE_TYPE *d;
HOST_WIDE_INT i, j;
enum machine_mode mode;
{
unsigned EMUSHORT df[NE], dg[NE];
HOST_WIDE_INT low, high;
int sign;
if (GET_MODE_CLASS (mode) != MODE_FLOAT)
abort ();
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);
/* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
Avoid double-rounding errors later by rounding off now from the
extra-wide internal format to the requested precision. */
switch (GET_MODE_BITSIZE (mode))
{
case 32:
etoe24 (dg, df);
e24toe (df, dg);
break;
case 64:
etoe53 (dg, df);
e53toe (df, dg);
break;
case 96:
etoe64 (dg, df);
e64toe (df, dg);
break;
case 128:
etoe113 (dg, df);
e113toe (df, dg);
break;
default:
abort ();
}
PUT_REAL (dg, d);
}
/* REAL_VALUE_FROM_UNSIGNED_INT macro. */
void
ereal_from_uint (d, i, j, mode)
REAL_VALUE_TYPE *d;
unsigned HOST_WIDE_INT i, j;
enum machine_mode mode;
{
unsigned EMUSHORT df[NE], dg[NE];
unsigned HOST_WIDE_INT low, high;
if (GET_MODE_CLASS (mode) != MODE_FLOAT)
abort ();
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);
/* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
Avoid double-rounding errors later by rounding off now from the
extra-wide internal format to the requested precision. */
switch (GET_MODE_BITSIZE (mode))
{
case 32:
etoe24 (dg, df);
e24toe (df, dg);
break;
case 64:
etoe53 (dg, df);
e53toe (df, dg);
break;
case 96:
etoe64 (dg, df);
e64toe (df, dg);
break;
case 128:
etoe113 (dg, df);
e113toe (df, dg);
break;
default:
abort ();
}
PUT_REAL (dg, d);
}
/* 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 HFmode:
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);
}
/* Try to change R into its exact multiplicative inverse in machine mode
MODE. Return nonzero function value if successful. */
int
exact_real_inverse (mode, r)
enum machine_mode mode;
REAL_VALUE_TYPE *r;
{
unsigned EMUSHORT e[NE], einv[NE];
REAL_VALUE_TYPE rinv;
int i;
GET_REAL (r, e);
/* Test for input in range. Don't transform IEEE special values. */
if (eisinf (e) || eisnan (e) || (ecmp (e, ezero) == 0))
return 0;
/* Test for a power of 2: all significand bits zero except the MSB.
We are assuming the target has binary (or hex) arithmetic. */
if (e[NE - 2] != 0x8000)
return 0;
for (i = 0; i < NE - 2; i++)
{
if (e[i] != 0)
return 0;
}
/* Compute the inverse and truncate it to the required mode. */
ediv (e, eone, einv);
PUT_REAL (einv, &rinv);
rinv = real_value_truncate (mode, rinv);
#ifdef CHECK_FLOAT_VALUE
/* This check is not redundant. It may, for example, flush
a supposedly IEEE denormal value to zero. */
i = 0;
if (CHECK_FLOAT_VALUE (mode, rinv, i))
return 0;
#endif
GET_REAL (&rinv, einv);
/* Check the bits again, because the truncation might have
generated an arbitrary saturation value on overflow. */
if (einv[NE - 2] != 0x8000)
return 0;
for (i = 0; i < NE - 2; i++)
{
if (einv[i] != 0)
return 0;
}
/* Fail if the computed inverse is out of range. */
if (eisinf (einv) || eisnan (einv) || (ecmp (einv, ezero) == 0))
return 0;
/* Output the reciprocal and return success flag. */
PUT_REAL (einv, r);
return 1;
}
#endif /* REAL_ARITHMETIC defined */
/* 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);
}
/* The following routines convert REAL_VALUE_TYPE to the various floating
point formats that are meaningful to supported computers.
The results are returned in 32-bit pieces, each piece stored in a `long'.
This is so they can be printed by statements like
fprintf (file, "%lx, %lx", L[0], L[1]);
that will work on both narrow- and wide-word host computers. */
/* Convert R to a 128-bit long double precision value. The output array L
contains four 32-bit pieces of the result, in the order they would appear
in memory. */
void
etartdouble (r, l)
REAL_VALUE_TYPE r;
long l[];
{
unsigned EMUSHORT e[NE];
GET_REAL (&r, e);
etoe113 (e, e);
endian (e, l, TFmode);
}
/* Convert R to a double extended precision value. The output array L
contains three 32-bit pieces of the result, in the order they would
appear in memory. */
void
etarldouble (r, l)
REAL_VALUE_TYPE r;
long l[];
{
unsigned EMUSHORT e[NE];
GET_REAL (&r, e);
etoe64 (e, e);
endian (e, l, XFmode);
}
/* Convert R to a double precision value. The output array L contains two
32-bit pieces of the result, in the order they would appear in memory. */
void
etardouble (r, l)
REAL_VALUE_TYPE r;
long l[];
{
unsigned EMUSHORT e[NE];
GET_REAL (&r, e);
etoe53 (e, e);
endian (e, l, DFmode);
}
/* Convert R to a single precision float value stored in the least-significant
bits of a `long'. */
long
etarsingle (r)
REAL_VALUE_TYPE r;
{
unsigned EMUSHORT e[NE];
long l;
GET_REAL (&r, e);
etoe24 (e, e);
endian (e, &l, SFmode);
return ((long) l);
}
/* Convert X to a decimal ASCII string S for output to an assembly
language file. Note, there is no standard way to spell infinity or
a NaN, so these values may require special treatment in the tm.h
macros. */
void
ereal_to_decimal (x, s)
REAL_VALUE_TYPE x;
char *s;
{
unsigned EMUSHORT e[NE];
GET_REAL (&x, e);
etoasc (e, s, 20);
}
/* Compare X and Y. Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
or -2 if either is a NaN. */
int
ereal_cmp (x, y)
REAL_VALUE_TYPE x, y;
{
unsigned EMUSHORT ex[NE], ey[NE];
GET_REAL (&x, ex);
GET_REAL (&y, ey);
return (ecmp (ex, ey));
}
/* Return 1 if the sign bit of X is set, else return 0. */
int
ereal_isneg (x)
REAL_VALUE_TYPE x;
{
unsigned EMUSHORT ex[NE];
GET_REAL (&x, ex);
return (eisneg (ex));
}
/* End of REAL_ARITHMETIC interface */
/*
Extended precision IEEE binary floating point arithmetic routines
Numbers are stored in C language as arrays of 16-bit unsigned
short integers. The arguments of the routines are pointers to
the arrays.
External e type data structure, similar to Intel 8087 chip
temporary real format but possibly with a larger significand:
NE-1 significand words (least significant word first,
most significant bit is normally set)
exponent (value = EXONE for 1.0,
top bit is the sign)
Internal exploded e-type data structure of a number (a "word" is 16 bits):
ei[0] sign word (0 for positive, 0xffff for negative)
ei[1] biased exponent (value = EXONE for the number 1.0)
ei[2] high guard word (always zero after normalization)
ei[3]
to ei[NI-2] significand (NI-4 significand words,
most significant word first,
most significant bit is set)
ei[NI-1] low guard word (0x8000 bit is rounding place)
Routines for external format e-type numbers
asctoe (string, e) ASCII string to extended double e type
asctoe64 (string, &d) ASCII string to long double
asctoe53 (string, &d) ASCII string to double
asctoe24 (string, &f) ASCII string to single
asctoeg (string, e, prec) ASCII string to specified precision
e24toe (&f, e) IEEE single precision to e type
e53toe (&d, e) IEEE double precision to e type
e64toe (&d, e) IEEE long double precision to e type
e113toe (&d, e) 128-bit long double precision to e type
eabs (e) absolute value
eadd (a, b, c) c = b + a
eclear (e) e = 0
ecmp (a, b) Returns 1 if a > b, 0 if a == b,
-1 if a < b, -2 if either a or b is a NaN.
ediv (a, b, c) c = b / a
efloor (a, b) truncate to integer, toward -infinity
efrexp (a, exp, s) extract exponent and significand
eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
einfin (e) set e to infinity, leaving its sign alone
eldexp (a, n, b) multiply by 2**n
emov (a, b) b = a
emul (a, b, c) c = b * a
eneg (e) e = -e
eround (a, b) b = nearest integer value to a
esub (a, b, c) c = b - a
e24toasc (&f, str, n) single to ASCII string, n digits after decimal
e53toasc (&d, str, n) double to ASCII string, n digits after decimal
e64toasc (&d, str, n) 80-bit long double to ASCII string
e113toasc (&d, str, n) 128-bit long double to ASCII string
etoasc (e, str, n) e to ASCII string, n digits after decimal
etoe24 (e, &f) convert e type to IEEE single precision
etoe53 (e, &d) convert e type to IEEE double precision
etoe64 (e, &d) convert e type to IEEE long double precision
ltoe (&l, e) HOST_WIDE_INT to e type
ultoe (&l, e) unsigned HOST_WIDE_INT to e type
eisneg (e) 1 if sign bit of e != 0, else 0
eisinf (e) 1 if e has maximum exponent (non-IEEE)
or is infinite (IEEE)
eisnan (e) 1 if e is a NaN
Routines for internal format exploded e-type numbers
eaddm (ai, bi) add significands, bi = bi + ai
ecleaz (ei) ei = 0
ecleazs (ei) set ei = 0 but leave its sign alone
ecmpm (ai, bi) compare significands, return 1, 0, or -1
edivm (ai, bi) divide significands, bi = bi / ai
emdnorm (ai,l,s,exp) normalize and round off
emovi (a, ai) convert external a to internal ai
emovo (ai, a) convert internal ai to external a
emovz (ai, bi) bi = ai, low guard word of bi = 0
emulm (ai, bi) multiply significands, bi = bi * ai
enormlz (ei) left-justify the significand
eshdn1 (ai) shift significand and guards down 1 bit
eshdn8 (ai) shift down 8 bits
eshdn6 (ai) shift down 16 bits
eshift (ai, n) shift ai n bits up (or down if n < 0)
eshup1 (ai) shift significand and guards up 1 bit
eshup8 (ai) shift up 8 bits
eshup6 (ai) shift up 16 bits
esubm (ai, bi) subtract significands, bi = bi - ai
eiisinf (ai) 1 if infinite
eiisnan (ai) 1 if a NaN
eiisneg (ai) 1 if sign bit of ai != 0, else 0
einan (ai) set ai = NaN
eiinfin (ai) set ai = infinity
The result is always normalized and rounded to NI-4 word precision
after each arithmetic operation.
Exception flags are NOT fully supported.
Signaling NaN's are NOT supported; they are treated the same
as quiet NaN's.
Define INFINITY for support of infinity; otherwise a
saturation arithmetic is implemented.
Define NANS for support of Not-a-Number items; otherwise the
arithmetic will never produce a NaN output, and might be confused
by a NaN input.
If NaN's are supported, the output of `ecmp (a,b)' is -2 if
either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
if in doubt.
Denormals are always supported here where appropriate (e.g., not
for conversion to DEC numbers). */
/* Definitions for error codes that are passed to the common error handling
routine mtherr.
For Digital Equipment PDP-11 and VAX computers, certain
IBM systems, and others that use numbers with a 56-bit
significand, the symbol DEC should be defined. In this
mode, most floating point constants are given as arrays
of octal integers to eliminate decimal to binary conversion
errors that might be introduced by the compiler.
For computers, such as IBM PC, that follow the IEEE
Standard for Binary Floating Point Arithmetic (ANSI/IEEE
Std 754-1985), the symbol IEEE should be defined.
These numbers have 53-bit significands. In this mode, constants
are provided as arrays of hexadecimal 16 bit integers.
The endian-ness of generated values is controlled by
REAL_WORDS_BIG_ENDIAN.
To accommodate other types of computer arithmetic, all
constants are also provided in a normal decimal radix
which one can hope are correctly converted to a suitable
format by the available C language compiler. To invoke
this mode, the symbol UNK is defined.
An important difference among these modes is a predefined
set of machine arithmetic constants for each. The numbers
MACHEP (the machine roundoff error), MAXNUM (largest number
represented), and several other parameters are preset by
the configuration symbol. Check the file const.c to
ensure that these values are correct for your computer.
For ANSI C compatibility, define ANSIC equal to 1. Currently
this affects only the atan2 function and others that use it. */
/* 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 e-type number X. */
static void
eclear (x)
register unsigned EMUSHORT *x;
{
register int i;
for (i = 0; i < NE; i++)
*x++ = 0;
}
/* Move e-type 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 e-type X. */
static void
eabs (x)
unsigned EMUSHORT x[];
{
/* sign is top bit of last word of external format */
x[NE - 1] &= 0x7fff;
}
/* Negate the e-type number X. */
static void
eneg (x)
unsigned EMUSHORT x[];
{
x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
}
/* Return 1 if sign bit of e-type number X is nonzero, else zero. */
static int
eisneg (x)
unsigned EMUSHORT x[];
{
if (x[NE - 1] & 0x8000)
return (1);
else
return (0);
}
/* Return 1 if e-type number X 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 e-type number X 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 an e-type number A, converting it to exploded e-type B. */
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 out exploded e-type number A, converting it to e type B. */
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 exploded e-type number XI. */
static void
ecleaz (xi)
register unsigned EMUSHORT *xi;
{
register int i;
for (i = 0; i < NI; i++)
*xi++ = 0;
}
/* Clear out exploded e-type XI, 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 exploded e-type 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 exploded e-type 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 exploded e-type X is a NaN. */
static int
eiisnan (x)
unsigned EMUSHORT x[];
{
int i;
if ((x[E] & 0x7fff) == 0x7fff)
{
for (i = M + 1; i < NI; i++)
{
if (x[i] != 0)
return (1);
}
}
return (0);
}
/* Return nonzero if sign of exploded e-type X is nonzero. */
static int
eiisneg (x)
unsigned EMUSHORT x[];
{
return x[0] != 0;
}
/* Fill exploded e-type X with infinity pattern.
This has maximum exponent and significand all zeros. */
static void
eiinfin (x)
unsigned EMUSHORT x[];
{
ecleaz (x);
x[E] = 0x7fff;
}
/* Return nonzero if exploded e-type X is infinite. */
static int
eiisinf (x)
unsigned EMUSHORT x[];
{
#ifdef NANS
if (eiisnan (x))
return (0);
#endif
if ((x[E] & 0x7fff) == 0x7fff)
return (1);
return (0);
}
/* Compare significands of numbers in internal exploded e-type format.
Guard words are included in the comparison.
Returns +1 if a > b
0 if a == b
-1 if a < b */
static int
ecmpm (a, b)
register unsigned EMUSHORT *a, *b;
{
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 of exploded e-type X 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 of exploded e-type X 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 of exploded e-type X 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 of exploded e-type X 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 of exploded e-type X 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 of exploded e-type X 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;
}
/* Add significands of exploded e-type X and Y. 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 of exploded e-type X and Y. 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, return e-type result to C. */
static void
m16m (a, b, c)
unsigned int a;
unsigned EMUSHORT b[], c[];
{
register unsigned EMUSHORT *pp;
register unsigned EMULONG carry;
unsigned EMUSHORT *ps;
unsigned EMUSHORT p[NI];
unsigned EMULONG 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 EMULONG) aa * *ps--;
carry = (m & 0xffff) + *pp;
*pp-- = (unsigned EMUSHORT)carry;
carry = (carry >> 16) + (m >> 16) + *pp;
*pp = (unsigned EMUSHORT)carry;
*(pp-1) = carry >> 16;
}
}
for (i=M; i<NI; i++)
c[i] = p[i];
}
/* Divide significands of exploded e-types NUM / DEN. Neither the
numerator NUM nor the denominator DEN is permitted to have its high guard
word nonzero. */
static int
edivm (den, num)
unsigned EMUSHORT den[], num[];
{
int i;
register unsigned EMUSHORT *p;
unsigned EMULONG tnum;
unsigned EMUSHORT j, tdenm, tquot;
unsigned EMUSHORT tprod[NI+1];
p = &equot[0];
*p++ = num[0];
*p++ = num[1];
for (i=M; i<NI; i++)
{
*p++ = 0;
}
eshdn1 (num);
tdenm = den[M+1];
for (i=M; i<NI; i++)
{
/* Find trial quotient digit (the radix is 65536). */
tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
/* Do not execute the divide instruction if it will overflow. */
if ((tdenm * 0xffffL) < tnum)
tquot = 0xffff;
else
tquot = tnum / tdenm;
/* Multiply denominator by trial quotient digit. */
m16m ((unsigned int)tquot, den, tprod);
/* The quotient digit may have been overestimated. */
if (ecmpm (tprod, num) > 0)
{
tquot -= 1;
esubm (den, tprod);
if (ecmpm (tprod, num) > 0)
{
tquot -= 1;
esubm (den, tprod);
}
}
esubm (tprod, num);
equot[i] = tquot;
eshup6(num);
}
/* test for nonzero remainder after roundoff bit */
p = &num[M];
j = 0;
for (i=M; i<NI; i++)
{
j |= *p++;
}
if (j)
j = 1;
for (i=0; i<NI; i++)
num[i] = equot[i];
return ((int)j);
}
/* Multiply significands of exploded e-type A and B, result in B. */
static int
emulm (a, b)
unsigned EMUSHORT a[], b[];
{
unsigned EMUSHORT *p, *q;
unsigned EMUSHORT pprod[NI];
unsigned EMUSHORT j;
int i;
equot[0] = b[0];
equot[1] = b[1];
for (i=M; i<NI; i++)
equot[i] = 0;
j = 0;
p = &a[NI-1];
q = &equot[NI-1];
for (i=M+1; i<NI; i++)
{
if (*p == 0)
{
--p;
}
else
{
m16m ((unsigned int) *p--, b, pprod);
eaddm(pprod, equot);
}
j |= *q;
eshdn6(equot);
}
for (i=0; i<NI; i++)
b[i] = equot[i];
/* return flag for lost nonzero bits */
return ((int)j);
}
#endif
/* Normalize and round off.
The internal format number to be rounded is S.
Input LOST is 0 if the value is exact. This is the so-called sticky bit.
Input SUBFLG indicates whether the number was obtained
by a subtraction operation. In that case if LOST is nonzero
then the number is slightly smaller than indicated.
Input EXP is the biased exponent, which may be negative.
the exponent field of S is ignored but is replaced by
EXP as adjusted by normalization and rounding.
Input RCNTRL is the rounding control. If it is nonzero, the
returned value will be rounded to RNDPRC bits.
For future reference: In order for emdnorm to round off denormal
significands at the right point, the input exponent must be
adjusted to be the actual value it would have after conversion to
the final floating point type. This adjustment has been
implemented for all type conversions (etoe53, etc.) and decimal
conversions, but not for the arithmetic functions (eadd, etc.).
Data types having standard 15-bit exponents are not affected by
this, but SFmode and DFmode are affected. For example, ediv with
rndprc = 24 will not round correctly to 24-bit precision if the
result is denormal. */
static int rlast = -1;
static int rw = 0;
static unsigned EMUSHORT rmsk = 0;
static unsigned EMUSHORT rmbit = 0;
static unsigned EMUSHORT rebit = 0;
static int re = 0;
static unsigned EMUSHORT rbit[NI];
static void
emdnorm (s, lost, subflg, exp, rcntrl)
unsigned EMUSHORT s[];
int lost;
int subflg;
EMULONG exp;
int rcntrl;
{
int i, j;
unsigned EMUSHORT r;
/* Normalize */
j = enormlz (s);
/* a blank significand could mean either zero or infinity. */
#ifndef INFINITY
if (j > NBITS)
{
ecleazs (s);
return;
}
#endif
exp -= j;
#ifndef INFINITY
if (exp >= 32767L)
goto overf;
#else
if ((j > NBITS) && (exp < 32767))
{
ecleazs (s);
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.
Intel long double denormals also lose one bit of precision. */
if ((exp <= 0) && (rndprc != NBITS)
&& ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
{
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:
/* Undo the temporary shift for denormal values. */
if ((exp <= 0) && (rndprc != NBITS)
&& ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
{
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. C = B - A, all e type 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. C = A + B, all e type. */
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);
}
/* Arithmetic common to both addition and subtraction. */
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 denormalized 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)
{
ltb += 1;
if (ltb >= 0x7fff)
{
eclear (c);
if (ai[0] != 0)
eneg (c);
einfin (c);
return;
}
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: C = B/A, all e type. */
static void
ediv (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
unsigned EMUSHORT ai[NI], bi[NI];
int i, sign;
EMULONG lt, lta, ltb;
/* IEEE says if result is not a NaN, the sign is "-" if and only if
operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
sign = eisneg(a) ^ eisneg(b);
#ifdef NANS
/* Return any NaN input. */
if (eisnan (a))
{
emov (a, c);
return;
}
if (eisnan (b))
{
emov (b, c);
return;
}
/* Zero over zero, or infinity over infinity, is a NaN. */
if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
|| (eisinf (a) && eisinf (b)))
{
mtherr ("ediv", INVALID);
enan (c, sign);
return;
}
#endif
/* Infinity over anything else is infinity. */
#ifdef INFINITY
if (eisinf (b))
{
einfin (c);
goto divsign;
}
/* Anything else over infinity is zero. */
if (eisinf (a))
{
eclear (c);
goto divsign;
}
#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);
goto divsign;
}
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;
}
}
/* Divide by zero is not an invalid operation.
It is a divide-by-zero operation! */
einfin (c);
mtherr ("ediv", SING);
goto divsign;
}
dnzro2:
i = edivm (ai, bi);
/* calculate exponent */
lt = ltb - lta + EXONE;
emdnorm (bi, i, 0, lt, 64);
emovo (bi, c);
divsign:
if (sign
#ifndef IEEE
&& (ecmp (c, ezero) != 0)
#endif
)
*(c+(NE-1)) |= 0x8000;
else
*(c+(NE-1)) &= ~0x8000;
}
/* Multiply e-types A and B, return e-type product C. */
static void
emul (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
unsigned EMUSHORT ai[NI], bi[NI];
int i, j, sign;
EMULONG lt, lta, ltb;
/* IEEE says if result is not a NaN, the sign is "-" if and only if
operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
sign = eisneg(a) ^ eisneg(b);
#ifdef NANS
/* NaN times anything is the same NaN. */
if (eisnan (a))
{
emov (a, c);
return;
}
if (eisnan (b))
{
emov (b, c);
return;
}
/* Zero times infinity is a NaN. */
if ((eisinf (a) && (ecmp (b, ezero) == 0))
|| (eisinf (b) && (ecmp (a, ezero) == 0)))
{
mtherr ("emul", INVALID);
enan (c, sign);
return;
}
#endif
/* Infinity times anything else is infinity. */
#ifdef INFINITY
if (eisinf (a) || eisinf (b))
{
einfin (c);
goto mulsign;
}
#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);
goto mulsign;
}
mnzer1:
if (bi[E] == 0)
{
for (i = 1; i < NI - 1; i++)
{
if (bi[i] != 0)
{
ltb -= enormlz (bi);
goto mnzer2;
}
}
eclear (c);
goto mulsign;
}
mnzer2:
/* Multiply significands */
j = emulm (ai, bi);
/* calculate exponent */
lt = lta + ltb - (EXONE - 1);
emdnorm (bi, j, 0, lt, 64);
emovo (bi, c);
mulsign:
if (sign
#ifndef IEEE
&& (ecmp (c, ezero) != 0)
#endif
)
*(c+(NE-1)) |= 0x8000;
else
*(c+(NE-1)) &= ~0x8000;
}
/* Convert double precision PE to e-type Y. */
static void
e53toe (pe, y)
unsigned EMUSHORT *pe, *y;
{
#ifdef DEC
dectoe (pe, y);
#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);
if (! REAL_WORDS_BIG_ENDIAN)
e += 3;
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
if (! REAL_WORDS_BIG_ENDIAN)
{
if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
|| (pe[1] != 0) || (pe[0] != 0))
{
enan (y, yy[0] != 0);
return;
}
}
else
{
if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
|| (pe[2] != 0) || (pe[3] != 0))
{
enan (y, yy[0] != 0);
return;
}
}
#endif /* NANS */
eclear (y);
einfin (y);
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 IEEE
if (! REAL_WORDS_BIG_ENDIAN)
{
*p++ = *(--e);
*p++ = *(--e);
*p++ = *(--e);
}
else
{
++e;
*p++ = *e++;
*p++ = *e++;
*p++ = *e++;
}
#endif
eshift (yy, -5);
if (denorm)
{ /* 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 */
}
/* Convert double extended precision float PE to e type Y. */
static void
e64toe (pe, y)
unsigned EMUSHORT *pe, *y;
{
unsigned EMUSHORT yy[NI];
unsigned EMUSHORT *e, *p, *q;
int i;
e = pe;
p = yy;
for (i = 0; i < NE - 5; i++)
*p++ = 0;
/* 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 IEEE
if (! REAL_WORDS_BIG_ENDIAN)
{
for (i = 0; i < 5; i++)
*p++ = *e++;
/* For denormal long double Intel format, shift significand up one
-- but only if the top significand bit is zero. A top bit of 1
is "pseudodenormal" when the exponent is zero. */
if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
{
unsigned EMUSHORT temp[NI];
emovi(yy, temp);
eshup1(temp);
emovo(temp,y);
return;
}
}
else
{
p = &yy[0] + (NE - 1);
#ifdef ARM_EXTENDED_IEEE_FORMAT
/* For ARMs, the exponent is in the lowest 15 bits of the word. */
*p-- = (e[0] & 0x8000) | (e[1] & 0x7ffff);
e += 2;
#else
*p-- = *e++;
++e;
#endif
for (i = 0; i < 4; i++)
*p-- = *e++;
}
#endif
#ifdef INFINITY
/* Point to the exponent field and check max exponent cases. */
p = &yy[NE - 1];
if ((*p & 0x7fff) == 0x7fff)
{
#ifdef NANS
if (! REAL_WORDS_BIG_ENDIAN)
{
for (i = 0; i < 4; i++)
{
if ((i != 3 && pe[i] != 0)
/* Anything but 0x8000 here, including 0, is a NaN. */
|| (i == 3 && pe[i] != 0x8000))
{
enan (y, (*p & 0x8000) != 0);
return;
}
}
}
else
{
#ifdef ARM_EXTENDED_IEEE_FORMAT
for (i = 2; i <= 5; i++)
{
if (pe[i] != 0)
{
enan (y, (*p & 0x8000) != 0);
return;
}
}
#else /* not ARM */
/* In Motorola extended precision format, the most significant
bit of an infinity mantissa could be either 1 or 0. It is
the lower order bits that tell whether the value is a NaN. */
if ((pe[2] & 0x7fff) != 0)
goto bigend_nan;
for (i = 3; i <= 5; i++)
{
if (pe[i] != 0)
{
bigend_nan:
enan (y, (*p & 0x8000) != 0);
return;
}
}
#endif /* not ARM */
}
#endif /* NANS */
eclear (y);
einfin (y);
if (*p & 0x8000)
eneg (y);
return;
}
#endif /* INFINITY */
p = yy;
q = y;
for (i = 0; i < NE; i++)
*q++ = *p++;
}
/* Convert 128-bit long double precision float PE to e type Y. */
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 IEEE
if (! REAL_WORDS_BIG_ENDIAN)
e += 7;
#endif
r = *e;
yy[0] = 0;
if (r & 0x8000)
yy[0] = 0xffff;
r &= 0x7fff;
#ifdef INFINITY
if (r == 0x7fff)
{
#ifdef NANS
if (! REAL_WORDS_BIG_ENDIAN)
{
for (i = 0; i < 7; i++)
{
if (pe[i] != 0)
{
enan (y, yy[0] != 0);
return;
}
}
}
else
{
for (i = 1; i < 8; i++)
{
if (pe[i] != 0)
{
enan (y, yy[0] != 0);
return;
}
}
}
#endif /* NANS */
eclear (y);
einfin (y);
if (yy[0])
eneg (y);
return;
}
#endif /* INFINITY */
yy[E] = r;
p = &yy[M + 1];
#ifdef IEEE
if (! REAL_WORDS_BIG_ENDIAN)
{
for (i = 0; i < 7; i++)
*p++ = *(--e);<