/* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
-and support for XFmode IEEE extended real floating point arithmetic.
-Contributed by Stephen L. Moshier (moshier@world.std.com).
-
- Copyright (C) 1993 Free Software Foundation, Inc.
+ and support for XFmode IEEE extended real floating point arithmetic.
+ Copyright (C) 1993, 1994, 1995, 1996 Free Software Foundation, Inc.
+ Contributed by Stephen L. Moshier (moshier@world.std.com).
This file is part of GNU CC.
You should have received a copy of the GNU General Public License
along with GNU CC; see the file COPYING. If not, write to
-the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
+the Free Software Foundation, 59 Temple Place - Suite 330,
+Boston, MA 02111-1307, USA. */
#include <stdio.h>
#include <errno.h>
/* To enable support of XFmode extended real floating point, define
LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
-To support cross compilation between IEEE and VAX floating
+To support cross compilation between IEEE, VAX and IBM floating
point formats, define REAL_ARITHMETIC in the tm.h file.
In either case the machine files (tm.h) must not contain any code
its decimal conversion functions can be used if desired (see
real.h).
-The first part of this file interfaces gcc to ieee.c, which is a
-floating point arithmetic suite that was not written with gcc in
-mind. The interface is followed by ieee.c itself and related
-items. Avoid changing ieee.c unless you have suitable test
-programs available. A special version of the PARANOIA floating
-point arithmetic tester, modified for this purpose, can be found
-on usc.edu : /pub/C-numanal/ieeetest.zoo. Some tutorial
-information on ieee.c is given in my book: S. L. Moshier,
-_Methods and Programs for Mathematical Functions_, Prentice-Hall
-or Simon & Schuster Int'l, 1989. A library of XFmode elementary
-transcendental functions can be obtained by ftp from
-research.att.com: netlib/cephes/ldouble.shar.Z */
-
+The first part of this file interfaces gcc to a floating point
+arithmetic suite that was not written with gcc in mind. Avoid
+changing the low-level arithmetic routines unless you have suitable
+test programs available. A special version of the PARANOIA floating
+point arithmetic tester, modified for this purpose, can be found on
+usc.edu: /pub/C-numanal/ieeetest.zoo. Other tests, and libraries of
+XFmode and TFmode transcendental functions, can be obtained by ftp from
+netlib.att.com: netlib/cephes. */
+\f
/* Type of computer arithmetic.
- * Only one of DEC, MIEEE, IBMPC, or UNK should get defined.
- * The following modification converts gcc macros into the ones
- * used by ieee.c.
- *
- * Note: long double formats differ between IBMPC and MIEEE
- * by more than just endian-ness.
- */
+ Only one of DEC, IBM, IEEE, or UNK should get defined.
+
+ `IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
+ to big-endian IEEE floating-point data structure. This definition
+ should work in SFmode `float' type and DFmode `double' type on
+ virtually all big-endian IEEE machines. If LONG_DOUBLE_TYPE_SIZE
+ has been defined to be 96, then IEEE also invokes the particular
+ XFmode (`long double' type) data structure used by the Motorola
+ 680x0 series processors.
+
+ `IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
+ little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
+ has been defined to be 96, then IEEE also invokes the particular
+ XFmode `long double' data structure used by the Intel 80x86 series
+ processors.
+
+ `DEC' refers specifically to the Digital Equipment Corp PDP-11
+ and VAX floating point data structure. This model currently
+ supports no type wider than DFmode.
+
+ `IBM' refers specifically to the IBM System/370 and compatible
+ floating point data structure. This model currently supports
+ no type wider than DFmode. The IBM conversions were contributed by
+ frank@atom.ansto.gov.au (Frank Crawford).
+
+ If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
+ then `long double' and `double' are both implemented, but they
+ both mean DFmode. In this case, the software floating-point
+ support available here is activated by writing
+ #define REAL_ARITHMETIC
+ in tm.h.
+
+ The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
+ and may deactivate XFmode since `long double' is used to refer
+ to both modes.
+
+ The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
+ contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
+ separate the floating point unit's endian-ness from that of
+ the integer addressing. This permits one to define a big-endian
+ FPU on a little-endian machine (e.g., ARM). An extension to
+ BYTES_BIG_ENDIAN may be required for some machines in the future.
+ These optional macros may be defined in tm.h. In real.h, they
+ default to WORDS_BIG_ENDIAN, etc., so there is no need to define
+ them for any normal host or target machine on which the floats
+ and the integers have the same endian-ness. */
+
+
+/* The following converts gcc macros into the ones used by this file. */
/* REAL_ARITHMETIC defined means that macros in real.h are
defined to call emulator functions. */
/* PDP-11, Pro350, VAX: */
#define DEC 1
#else /* it's not VAX */
+#if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
+/* IBM System/370 style */
+#define IBM 1
+#else /* it's also not an IBM */
#if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
-#if WORDS_BIG_ENDIAN
-/* Motorola IEEE, high order words come first (Sun workstation): */
-#define MIEEE 1
-#else /* not big-endian */
-/* Intel IEEE, low order words come first:
- */
-#define IBMPC 1
-#endif /* big-endian */
+#define IEEE
#else /* it's not IEEE either */
-/* UNKnown arithmetic. We don't support this and can't go on. */
+/* UNKnown arithmetic. We don't support this and can't go on. */
unknown arithmetic type
#define UNK 1
#endif /* not IEEE */
+#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
-#ifdef HOST_WORDS_BIG_ENDIAN
-#define MIEEE 1
-#else /* not big-endian */
-#define IBMPC 1
-#endif /* big-endian */
+#define IEEE
#else /* it's not IEEE either */
unknown arithmetic type
#define UNK 1
#endif /* not IEEE */
+#endif /* not IBM */
#endif /* not VAX */
+#define REAL_WORDS_BIG_ENDIAN HOST_FLOAT_WORDS_BIG_ENDIAN
+
#endif /* REAL_ARITHMETIC not defined */
-/* Define for support of infinity. */
-#ifndef DEC
+/* Define INFINITY for support of infinity.
+ Define NANS for support of Not-a-Number's (NaN's). */
+#if !defined(DEC) && !defined(IBM)
#define INFINITY
+#define NANS
#endif
-
-/* ehead.h
- *
- * Include file for extended precision arithmetic programs.
- */
-
-/* Number of 16 bit words in external e type format */
-#define NE 6
-
-/* Number of 16 bit words in internal format */
-#define NI (NE+3)
-
-/* Array offset to exponent */
-#define E 1
-
-/* Array offset to high guard word */
-#define M 2
-
-/* Number of bits of precision */
-#define NBITS ((NI-4)*16)
-
-/* Maximum number of decimal digits in ASCII conversion
- * = NBITS*log10(2)
- */
-#define NDEC (NBITS*8/27)
-
-/* The exponent of 1.0 */
-#define EXONE (0x3fff)
-
+/* Support of NaNs requires support of infinity. */
+#ifdef NANS
+#ifndef INFINITY
+#define INFINITY
+#endif
+#endif
+\f
/* Find a host integer type that is at least 16 bits wide,
- and another type at least twice whatever that size is. */
+ and another type at least twice whatever that size is. */
#if HOST_BITS_PER_CHAR >= 16
#define EMUSHORT char
#define EMUSHORT_SIZE HOST_BITS_PER_LONG
#define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
#else
-/* You will have to modify this program to have a smaller unit size. */
+/* You will have to modify this program to have a smaller unit size. */
#define EMU_NON_COMPILE
#endif
#endif
#if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
#define EMULONG long long int
#else
-/* You will have to modify this program to have a smaller unit size. */
+/* You will have to modify this program to have a smaller unit size. */
#define EMU_NON_COMPILE
#endif
#endif
#endif
-/* The host interface doesn't work if no 16-bit size exists. */
+/* The host interface doesn't work if no 16-bit size exists. */
#if EMUSHORT_SIZE != 16
#define EMU_NON_COMPILE
#endif
-/* OK to continue compilation. */
+/* OK to continue compilation. */
#ifndef EMU_NON_COMPILE
/* Construct macros to translate between REAL_VALUE_TYPE and e type.
in memory, with no holes. */
#if LONG_DOUBLE_TYPE_SIZE == 96
-#define GET_REAL(r,e) bcopy (r, e, 2*NE)
-#define PUT_REAL(e,r) bcopy (e, r, 2*NE)
+/* Number of 16 bit words in external e type format */
+#define NE 6
+#define MAXDECEXP 4932
+#define MINDECEXP -4956
+#define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
+#define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
#else /* no XFmode */
-
+#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. */
-
-#if defined (HOST_WORDS_BIG_ENDIAN) == WORDS_BIG_ENDIAN
-#define GET_REAL(r,e) e53toe ((r), (e))
-#define PUT_REAL(e,r) etoe53 ((e), (r))
-
-#else /* endian-ness differs */
-/* emulator uses target endian-ness internally */
-#define GET_REAL(r,e) \
-do { EMUSHORT w[4]; \
- w[3] = ((EMUSHORT *) r)[0]; \
- w[2] = ((EMUSHORT *) r)[1]; \
- w[1] = ((EMUSHORT *) r)[2]; \
- w[0] = ((EMUSHORT *) r)[3]; \
- e53toe (w, (e)); } while (0)
-
-#define PUT_REAL(e,r) \
-do { EMUSHORT w[4]; \
- etoe53 ((e), w); \
- *((EMUSHORT *) r) = w[3]; \
- *((EMUSHORT *) r + 1) = w[2]; \
- *((EMUSHORT *) r + 2) = w[1]; \
- *((EMUSHORT *) r + 3) = w[0]; } while (0)
-
-#endif /* endian-ness differs */
+ but host stores it in host endian-ness. */
+
+#define GET_REAL(r,e) \
+do { \
+ if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
+ e53toe ((unsigned EMUSHORT *) (r), (e)); \
+ else \
+ { \
+ unsigned EMUSHORT w[4]; \
+ w[3] = ((EMUSHORT *) r)[0]; \
+ w[2] = ((EMUSHORT *) r)[1]; \
+ w[1] = ((EMUSHORT *) r)[2]; \
+ w[0] = ((EMUSHORT *) r)[3]; \
+ e53toe (w, (e)); \
+ } \
+ } while (0)
+
+#define PUT_REAL(e,r) \
+do { \
+ if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
+ etoe53 ((e), (unsigned EMUSHORT *) (r)); \
+ else \
+ { \
+ unsigned EMUSHORT w[4]; \
+ etoe53 ((e), w); \
+ *((EMUSHORT *) r) = w[3]; \
+ *((EMUSHORT *) r + 1) = w[2]; \
+ *((EMUSHORT *) r + 2) = w[1]; \
+ *((EMUSHORT *) r + 3) = w[0]; \
+ } \
+ } while (0)
#else /* not REAL_ARITHMETIC */
/* emulator uses host format */
-#define GET_REAL(r,e) e53toe ((r), (e))
-#define PUT_REAL(e,r) etoe53 ((e), (r))
+#define GET_REAL(r,e) e53toe ((unsigned EMUSHORT *) (r), (e))
+#define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))
#endif /* not REAL_ARITHMETIC */
+#endif /* not TFmode */
#endif /* no XFmode */
-void warning ();
+
+/* Number of 16 bit words in internal format */
+#define NI (NE+3)
+
+/* Array offset to exponent */
+#define E 1
+
+/* Array offset to high guard word */
+#define M 2
+
+/* Number of bits of precision */
+#define NBITS ((NI-4)*16)
+
+/* Maximum number of decimal digits in ASCII conversion
+ * = NBITS*log10(2)
+ */
+#define NDEC (NBITS*8/27)
+
+/* The exponent of 1.0 */
+#define EXONE (0x3fff)
+
extern int extra_warnings;
-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
+static void endian PROTO((unsigned EMUSHORT *, long *,
+ enum machine_mode));
+static void eclear PROTO((unsigned EMUSHORT *));
+static void emov PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void eabs PROTO((unsigned EMUSHORT *));
+static void eneg PROTO((unsigned EMUSHORT *));
+static int eisneg PROTO((unsigned EMUSHORT *));
+static int eisinf PROTO((unsigned EMUSHORT *));
+static int eisnan PROTO((unsigned EMUSHORT *));
+static void einfin PROTO((unsigned EMUSHORT *));
+static void enan PROTO((unsigned EMUSHORT *, int));
+static void emovi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void emovo PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ecleaz PROTO((unsigned EMUSHORT *));
+static void ecleazs PROTO((unsigned EMUSHORT *));
+static void emovz PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void einan PROTO((unsigned EMUSHORT *));
+static int eiisnan PROTO((unsigned EMUSHORT *));
+static int eiisneg PROTO((unsigned EMUSHORT *));
+static void eiinfin PROTO((unsigned EMUSHORT *));
+static int eiisinf PROTO((unsigned EMUSHORT *));
+static int ecmpm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void eshdn1 PROTO((unsigned EMUSHORT *));
+static void eshup1 PROTO((unsigned EMUSHORT *));
+static void eshdn8 PROTO((unsigned EMUSHORT *));
+static void eshup8 PROTO((unsigned EMUSHORT *));
+static void eshup6 PROTO((unsigned EMUSHORT *));
+static void eshdn6 PROTO((unsigned EMUSHORT *));
+static void eaddm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));\f
+static void esubm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void m16m PROTO((unsigned int, unsigned short *,
+ unsigned short *));
+static int edivm PROTO((unsigned short *, unsigned short *));
+static int emulm PROTO((unsigned short *, unsigned short *));
+static void emdnorm PROTO((unsigned EMUSHORT *, int, int, EMULONG, int));
+static void esub PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void eadd PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void eadd1 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void ediv PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void emul PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void e53toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void e64toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void e113toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void e24toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static int ecmp PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void eround PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ltoe PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *));
+static void ultoe PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *));
+static void eifrac PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *,
+ unsigned EMUSHORT *));
+static void euifrac PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *,
+ unsigned EMUSHORT *));
+static int eshift PROTO((unsigned EMUSHORT *, int));
+static int enormlz PROTO((unsigned EMUSHORT *));
+static void e24toasc PROTO((unsigned EMUSHORT *, char *, int));
+static void e53toasc PROTO((unsigned EMUSHORT *, char *, int));
+static void e64toasc PROTO((unsigned EMUSHORT *, char *, int));
+static void e113toasc PROTO((unsigned EMUSHORT *, char *, int));
+static void etoasc PROTO((unsigned EMUSHORT *, char *, int));
+static void asctoe24 PROTO((char *, unsigned EMUSHORT *));
+static void asctoe53 PROTO((char *, unsigned EMUSHORT *));
+static void asctoe64 PROTO((char *, unsigned EMUSHORT *));
+static void asctoe113 PROTO((char *, unsigned EMUSHORT *));
+static void asctoe PROTO((char *, unsigned EMUSHORT *));
+static void asctoeg PROTO((char *, unsigned EMUSHORT *, int));
+static void efloor PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void efrexp PROTO((unsigned EMUSHORT *, int *,
+ unsigned EMUSHORT *));
+static void eldexp PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *));
+static void eremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void eiremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void mtherr PROTO((char *, int));
+static void dectoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etodec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void todec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ibmtoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ enum machine_mode));
+static void etoibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ enum machine_mode));
+static void toibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ enum machine_mode));
+static void make_nan PROTO((unsigned EMUSHORT *, int, enum machine_mode));
+static void uditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoudi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etodi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void esqrt PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+\f
+/* Copy 32-bit numbers obtained from array containing 16-bit numbers,
+ swapping ends if required, into output array of longs. The
+ result is normally passed to fprintf by the ASM_OUTPUT_ macros. */
+
+static void
endian (e, x, mode)
unsigned EMUSHORT e[];
long x[];
{
unsigned long th, t;
-#if WORDS_BIG_ENDIAN
- switch (mode)
+ if (REAL_WORDS_BIG_ENDIAN)
{
+ switch (mode)
+ {
- case XFmode:
-
- /* Swap halfwords in the third long. */
- th = (unsigned long) e[4] & 0xffff;
- t = (unsigned long) e[5] & 0xffff;
- t |= th << 16;
- x[2] = (long) t;
- /* fall into the double case */
-
- case DFmode:
-
- /* swap halfwords in the second word */
- th = (unsigned long) e[2] & 0xffff;
- t = (unsigned long) e[3] & 0xffff;
- t |= th << 16;
- x[1] = (long) t;
- /* fall into the float case */
-
- case SFmode:
-
- /* swap halfwords in the first word */
- th = (unsigned long) e[0] & 0xffff;
- t = (unsigned long) e[1] & 0xffff;
- t |= th << 16;
- x[0] = t;
- break;
+ case TFmode:
+ /* Swap halfwords in the fourth long. */
+ th = (unsigned long) e[6] & 0xffff;
+ t = (unsigned long) e[7] & 0xffff;
+ t |= th << 16;
+ x[3] = (long) t;
+
+ case XFmode:
+
+ /* Swap halfwords in the third long. */
+ th = (unsigned long) e[4] & 0xffff;
+ t = (unsigned long) e[5] & 0xffff;
+ t |= th << 16;
+ x[2] = (long) t;
+ /* fall into the double case */
+
+ case DFmode:
+
+ /* swap halfwords in the second word */
+ th = (unsigned long) e[2] & 0xffff;
+ t = (unsigned long) e[3] & 0xffff;
+ t |= th << 16;
+ x[1] = (long) t;
+ /* fall into the float case */
+
+ case HFmode:
+ case SFmode:
+
+ /* swap halfwords in the first word */
+ th = (unsigned long) e[0] & 0xffff;
+ t = (unsigned long) e[1] & 0xffff;
+ t |= th << 16;
+ x[0] = (long) t;
+ break;
- default:
- abort ();
+ default:
+ abort ();
+ }
}
-
-#else
-
- /* Pack the output array without swapping. */
-
- switch (mode)
+ else
{
+ /* Pack the output array without swapping. */
- case XFmode:
-
- /* Pack the third long.
- Each element of the input REAL_VALUE_TYPE array has 16 bit useful bits
- in it. */
- th = (unsigned long) e[5] & 0xffff;
- t = (unsigned long) e[4] & 0xffff;
- t |= th << 16;
- x[2] = (long) t;
- /* fall into the double case */
-
- case DFmode:
-
- /* pack the second long */
- th = (unsigned long) e[3] & 0xffff;
- t = (unsigned long) e[2] & 0xffff;
- t |= th << 16;
- x[1] = (long) t;
- /* fall into the float case */
-
- case SFmode:
+ switch (mode)
+ {
- /* pack the first long */
- th = (unsigned long) e[1] & 0xffff;
- t = (unsigned long) e[0] & 0xffff;
- t |= th << 16;
- x[0] = t;
- break;
+ case TFmode:
+
+ /* Pack the fourth long. */
+ th = (unsigned long) e[7] & 0xffff;
+ t = (unsigned long) e[6] & 0xffff;
+ t |= th << 16;
+ x[3] = (long) t;
+
+ case XFmode:
+
+ /* Pack the third long.
+ Each element of the input REAL_VALUE_TYPE array has 16 useful bits
+ in it. */
+ th = (unsigned long) e[5] & 0xffff;
+ t = (unsigned long) e[4] & 0xffff;
+ t |= th << 16;
+ x[2] = (long) t;
+ /* fall into the double case */
+
+ case DFmode:
+
+ /* pack the second long */
+ th = (unsigned long) e[3] & 0xffff;
+ t = (unsigned long) e[2] & 0xffff;
+ t |= th << 16;
+ x[1] = (long) t;
+ /* fall into the float case */
+
+ case HFmode:
+ case SFmode:
+
+ /* pack the first long */
+ th = (unsigned long) e[1] & 0xffff;
+ t = (unsigned long) e[0] & 0xffff;
+ t |= th << 16;
+ x[0] = (long) t;
+ break;
- default:
- abort ();
+ default:
+ abort ();
+ }
}
-
-#endif
}
-/* This is the implementation of the REAL_ARITHMETIC macro.
- */
+/* This is the implementation of the REAL_ARITHMETIC macro. */
+
void
earith (value, icode, r1, r2)
REAL_VALUE_TYPE *value;
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 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;
}
-/* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT
- * implements REAL_VALUE_FIX_TRUNCATE (x) (etrunci (x))
- */
+/* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
+ implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */
+
REAL_VALUE_TYPE
etrunci (x)
REAL_VALUE_TYPE x;
{
unsigned EMUSHORT f[NE], g[NE];
REAL_VALUE_TYPE r;
- long l;
+ 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);
}
-/* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT
- * implements REAL_VALUE_UNSIGNED_FIX_TRUNCATE (x) (etruncui (x))
- */
+/* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
+ implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */
+
REAL_VALUE_TYPE
etruncui (x)
REAL_VALUE_TYPE x;
{
unsigned EMUSHORT f[NE], g[NE];
REAL_VALUE_TYPE r;
- unsigned long l;
+ 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);
}
-/* This is the REAL_VALUE_ATOF function.
- * It converts a decimal string to binary, rounding off
- * as indicated by the machine_mode argument. Then it
- * promotes the rounded value to REAL_VALUE_TYPE.
- */
+/* This is the REAL_VALUE_ATOF function. It converts a decimal string to
+ binary, rounding off as indicated by the machine_mode argument. Then it
+ promotes the rounded value to REAL_VALUE_TYPE. */
+
REAL_VALUE_TYPE
ereal_atof (s, t)
char *s;
switch (t)
{
+ case HFmode:
case SFmode:
asctoe24 (s, tem);
e24toe (tem, e);
asctoe64 (s, tem);
e64toe (tem, e);
break;
+ case TFmode:
+ asctoe113 (s, tem);
+ e113toe (tem, e);
+ break;
default:
asctoe (s, e);
}
}
-/* Expansion of REAL_NEGATE.
- */
+/* Expansion of REAL_NEGATE. */
+
REAL_VALUE_TYPE
ereal_negate (x)
REAL_VALUE_TYPE x;
}
-/* Round real to int
- * implements REAL_VALUE_FIX (x) (eroundi (x))
- * The type of rounding is left unspecified by real.h.
- * It is implemented here as round to nearest (add .5 and chop).
- */
-int
-eroundi (x)
+/* Round real toward zero to HOST_WIDE_INT;
+ implements REAL_VALUE_FIX (x). */
+
+HOST_WIDE_INT
+efixi (x)
REAL_VALUE_TYPE x;
{
unsigned EMUSHORT f[NE], g[NE];
- EMULONG l;
+ HOST_WIDE_INT l;
GET_REAL (&x, f);
- eround (f, g);
- eifrac (g, &l, f);
- return ((int) l);
+#ifdef NANS
+ if (eisnan (f))
+ {
+ warning ("conversion from NaN to int");
+ return (-1);
+ }
+#endif
+ eifrac (f, &l, g);
+ return l;
}
-/* Round real to nearest unsigned int
- * implements REAL_VALUE_UNSIGNED_FIX (x) ((unsigned int) eroundi (x))
- * Negative input returns zero.
- * The type of rounding is left unspecified by real.h.
- * It is implemented here as round to nearest (add .5 and chop).
- */
-unsigned int
-eroundui (x)
+/* Round real toward zero to unsigned HOST_WIDE_INT
+ implements REAL_VALUE_UNSIGNED_FIX (x).
+ Negative input returns zero. */
+
+unsigned HOST_WIDE_INT
+efixui (x)
REAL_VALUE_TYPE x;
{
unsigned EMUSHORT f[NE], g[NE];
- unsigned EMULONG l;
+ unsigned HOST_WIDE_INT l;
GET_REAL (&x, f);
- eround (f, g);
- euifrac (g, &l, f);
- return ((unsigned int)l);
+#ifdef NANS
+ if (eisnan (f))
+ {
+ warning ("conversion from NaN to unsigned int");
+ return (-1);
+ }
+#endif
+ euifrac (f, &l, g);
+ return l;
}
-/* REAL_VALUE_FROM_INT macro.
- */
+/* REAL_VALUE_FROM_INT macro. */
+
void
-ereal_from_int (d, i, j)
+ereal_from_int (d, i, j, mode)
REAL_VALUE_TYPE *d;
- long i, j;
+ HOST_WIDE_INT i, j;
+ enum machine_mode mode;
{
unsigned EMUSHORT df[NE], dg[NE];
- long low, high;
+ HOST_WIDE_INT low, high;
int sign;
+ if (GET_MODE_CLASS (mode) != MODE_FLOAT)
+ abort ();
sign = 0;
low = i;
if ((high = j) < 0)
else
high += 1;
}
- eldexp (eone, HOST_BITS_PER_LONG, df);
- ultoe (&high, dg);
+ eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
+ ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
emul (dg, df, dg);
- ultoe (&low, df);
+ 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.
- */
+/* REAL_VALUE_FROM_UNSIGNED_INT macro. */
+
void
-ereal_from_uint (d, i, j)
+ereal_from_uint (d, i, j, mode)
REAL_VALUE_TYPE *d;
- unsigned long i, j;
+ unsigned HOST_WIDE_INT i, j;
+ enum machine_mode mode;
{
unsigned EMUSHORT df[NE], dg[NE];
- unsigned long low, high;
+ unsigned HOST_WIDE_INT low, high;
+ if (GET_MODE_CLASS (mode) != MODE_FLOAT)
+ abort ();
low = i;
high = j;
- eldexp (eone, HOST_BITS_PER_LONG, df);
+ eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
ultoe (&high, dg);
emul (dg, df, dg);
ultoe (&low, df);
eadd (df, dg, dg);
+
+ /* 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
- */
+/* REAL_VALUE_TO_INT macro. */
+
void
ereal_to_int (low, high, rr)
- long *low, *high;
+ HOST_WIDE_INT *low, *high;
REAL_VALUE_TYPE rr;
{
unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
int s;
GET_REAL (&rr, d);
+#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_LONG, df);
+ eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
ediv (df, d, dg); /* dg = d / 2^32 is the high word */
- euifrac (dg, high, dh);
+ euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
emul (df, dh, dg); /* fractional part is the low word */
- euifrac (dg, low, dh);
+ euifrac (dg, (unsigned HOST_WIDE_INT *)low, dh);
if (s)
{
/* complement and add 1 */
}
-/* REAL_VALUE_LDEXP macro.
- */
+/* REAL_VALUE_LDEXP macro. */
+
REAL_VALUE_TYPE
ereal_ldexp (x, n)
REAL_VALUE_TYPE x;
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. */
+ of the same names may be defined in fold-const.c. */
+
#ifdef REAL_ARITHMETIC
-/* Check for infinity in a REAL_VALUE_TYPE. */
+/* Check for infinity in a REAL_VALUE_TYPE. */
+
int
target_isinf (x)
REAL_VALUE_TYPE x;
#endif
}
-
-/* Check whether an IEEE double precision number is a NaN. */
+/* Check whether a REAL_VALUE_TYPE item is a NaN. */
int
target_isnan (x)
REAL_VALUE_TYPE x;
{
+ unsigned EMUSHORT e[NE];
+
+#ifdef NANS
+ GET_REAL (&x, e);
+ return (eisnan (e));
+#else
return (0);
+#endif
}
-/* Check for a negative IEEE double precision number.
- * this means strictly less than zero, not -0.
- */
+/* Check for a negative REAL_VALUE_TYPE number.
+ This just checks the sign bit, so that -0 counts as negative. */
int
target_negative (x)
REAL_VALUE_TYPE x;
{
- unsigned EMUSHORT e[NE];
-
- GET_REAL (&x, e);
- if (ecmp (e, ezero) < 0)
- return (1);
- return (0);
+ return ereal_isneg (x);
}
/* Expansion of REAL_VALUE_TRUNCATE.
- * The result is in floating point, rounded to nearest or even.
- */
+ The result is in floating point, rounded to nearest or even. */
+
REAL_VALUE_TYPE
real_value_truncate (mode, arg)
enum machine_mode mode;
REAL_VALUE_TYPE 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);
e53toe (t, t);
break;
+ case HFmode:
case SFmode:
etoe24 (e, t);
e24toe (t, t);
break;
case SImode:
- r = etrunci (e);
+ 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:
- abort ();
+ 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 */
-/* Target values are arrays of host longs. A long is guaranteed
- to be at least 32 bits wide. */
+/* Used for debugging--print the value of R in human-readable format
+ on stderr. */
+
+void
+debug_real (r)
+ REAL_VALUE_TYPE r;
+{
+ char dstr[30];
+
+ REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
+ fprintf (stderr, "%s", dstr);
+}
+
+\f
+/* The following routines convert REAL_VALUE_TYPE to the various floating
+ point formats that are meaningful to supported computers.
+
+ The results are returned in 32-bit pieces, each piece stored in a `long'.
+ This is so they can be printed by statements like
+
+ fprintf (file, "%lx, %lx", L[0], L[1]);
+
+ that will work on both narrow- and wide-word host computers. */
+
+/* Convert R to a 128-bit long double precision value. The output array L
+ contains four 32-bit pieces of the result, in the order they would appear
+ in memory. */
+
+void
+etartdouble (r, l)
+ REAL_VALUE_TYPE r;
+ long l[];
+{
+ unsigned EMUSHORT e[NE];
+
+ GET_REAL (&r, e);
+ etoe113 (e, e);
+ endian (e, l, TFmode);
+}
+
+/* Convert R to a double extended precision value. The output array L
+ contains three 32-bit pieces of the result, in the order they would
+ appear in memory. */
+
void
etarldouble (r, l)
REAL_VALUE_TYPE r;
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;
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];
- unsigned long l;
+ long l;
GET_REAL (&r, e);
etoe24 (e, e);
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;
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;
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;
}
/* End of REAL_ARITHMETIC interface */
-
-/* ieee.c
- *
- * Extended precision IEEE binary floating point arithmetic routines
- *
- * Numbers are stored in C language as arrays of 16-bit unsigned
- * short integers. The arguments of the routines are pointers to
- * the arrays.
- *
- *
- * External e type data structure, simulates Intel 8087 chip
- * temporary real format but possibly with a larger significand:
- *
- * NE-1 significand words (least significant word first,
- * most significant bit is normally set)
- * exponent (value = EXONE for 1.0,
- * top bit is the sign)
- *
- *
- * Internal data structure of a number (a "word" is 16 bits):
- *
- * ei[0] sign word (0 for positive, 0xffff for negative)
- * ei[1] biased exponent (value = EXONE for the number 1.0)
- * ei[2] high guard word (always zero after normalization)
- * ei[3]
- * to ei[NI-2] significand (NI-4 significand words,
- * most significant word first,
- * most significant bit is set)
- * ei[NI-1] low guard word (0x8000 bit is rounding place)
- *
- *
- *
- * Routines for external format numbers
- *
- * asctoe (string, e) ASCII string to extended double e type
- * asctoe64 (string, &d) ASCII string to long double
- * asctoe53 (string, &d) ASCII string to double
- * asctoe24 (string, &f) ASCII string to single
- * asctoeg (string, e, prec) ASCII string to specified precision
- * e24toe (&f, e) IEEE single precision to e type
- * e53toe (&d, e) IEEE double precision to e type
- * e64toe (&d, e) IEEE long double precision to e type
- * eabs (e) absolute value
- * eadd (a, b, c) c = b + a
- * eclear (e) e = 0
- * ecmp (a, b) compare a to b, return 1, 0, or -1
- * ediv (a, b, c) c = b / a
- * efloor (a, b) truncate to integer, toward -infinity
- * efrexp (a, exp, s) extract exponent and significand
- * eifrac (e, &l, frac) e to long integer and e type fraction
- * euifrac (e, &l, frac) e to unsigned long integer and e type fraction
- * einfin (e) set e to infinity, leaving its sign alone
- * eldexp (a, n, b) multiply by 2**n
- * emov (a, b) b = a
- * emul (a, b, c) c = b * a
- * eneg (e) e = -e
- * eround (a, b) b = nearest integer value to a
- * esub (a, b, c) c = b - a
- * e24toasc (&f, str, n) single to ASCII string, n digits after decimal
- * e53toasc (&d, str, n) double to ASCII string, n digits after decimal
- * e64toasc (&d, str, n) long double to ASCII string
- * etoasc (e, str, n) e to ASCII string, n digits after decimal
- * etoe24 (e, &f) convert e type to IEEE single precision
- * etoe53 (e, &d) convert e type to IEEE double precision
- * etoe64 (e, &d) convert e type to IEEE long double precision
- * ltoe (&l, e) long (32 bit) integer to e type
- * ultoe (&l, e) unsigned long (32 bit) integer to e type
- * eisneg (e) 1 if sign bit of e != 0, else 0
- * eisinf (e) 1 if e has maximum exponent
- *
- *
- * Routines for internal format numbers
- *
- * eaddm (ai, bi) add significands, bi = bi + ai
- * ecleaz (ei) ei = 0
- * ecleazs (ei) set ei = 0 but leave its sign alone
- * ecmpm (ai, bi) compare significands, return 1, 0, or -1
- * edivm (ai, bi) divide significands, bi = bi / ai
- * emdnorm (ai,l,s,exp) normalize and round off
- * emovi (a, ai) convert external a to internal ai
- * emovo (ai, a) convert internal ai to external a
- * emovz (ai, bi) bi = ai, low guard word of bi = 0
- * emulm (ai, bi) multiply significands, bi = bi * ai
- * enormlz (ei) left-justify the significand
- * eshdn1 (ai) shift significand and guards down 1 bit
- * eshdn8 (ai) shift down 8 bits
- * eshdn6 (ai) shift down 16 bits
- * eshift (ai, n) shift ai n bits up (or down if n < 0)
- * eshup1 (ai) shift significand and guards up 1 bit
- * eshup8 (ai) shift up 8 bits
- * eshup6 (ai) shift up 16 bits
- * esubm (ai, bi) subtract significands, bi = bi - ai
- *
- *
- * The result is always normalized and rounded to NI-4 word precision
- * after each arithmetic operation.
- *
- * Exception flags and NaNs are NOT fully supported.
- * This arithmetic should never produce a NaN output, but it might
- * be confused by a NaN input.
- * Define INFINITY in mconf.h for support of infinity; otherwise a
- * saturation arithmetic is implemented.
- * Denormals are always supported here where appropriate (e.g., not
- * for conversion to DEC numbers).
- *
- */
-
+\f
/*
- * Revision history:
- *
- * 5 Jan 84 PDP-11 assembly language version
- * 2 Mar 86 fixed bug in asctoq
- * 6 Dec 86 C language version
- * 30 Aug 88 100 digit version, improved rounding
- * 15 May 92 80-bit long double support
- *
- * Author: S. L. Moshier.
- */
+ Extended precision IEEE binary floating point arithmetic routines
+
+ Numbers are stored in C language as arrays of 16-bit unsigned
+ short integers. The arguments of the routines are pointers to
+ the arrays.
+
+ External e type data structure, similar to Intel 8087 chip
+ temporary real format but possibly with a larger significand:
+
+ NE-1 significand words (least significant word first,
+ most significant bit is normally set)
+ exponent (value = EXONE for 1.0,
+ top bit is the sign)
+
+
+ Internal exploded e-type data structure of a number (a "word" is 16 bits):
+
+ ei[0] sign word (0 for positive, 0xffff for negative)
+ ei[1] biased exponent (value = EXONE for the number 1.0)
+ ei[2] high guard word (always zero after normalization)
+ ei[3]
+ to ei[NI-2] significand (NI-4 significand words,
+ most significant word first,
+ most significant bit is set)
+ ei[NI-1] low guard word (0x8000 bit is rounding place)
+
+
+
+ Routines for external format e-type numbers
+
+ asctoe (string, e) ASCII string to extended double e type
+ asctoe64 (string, &d) ASCII string to long double
+ asctoe53 (string, &d) ASCII string to double
+ asctoe24 (string, &f) ASCII string to single
+ asctoeg (string, e, prec) ASCII string to specified precision
+ e24toe (&f, e) IEEE single precision to e type
+ e53toe (&d, e) IEEE double precision to e type
+ e64toe (&d, e) IEEE long double precision to e type
+ e113toe (&d, e) 128-bit long double precision to e type
+ eabs (e) absolute value
+ eadd (a, b, c) c = b + a
+ eclear (e) e = 0
+ ecmp (a, b) Returns 1 if a > b, 0 if a == b,
+ -1 if a < b, -2 if either a or b is a NaN.
+ ediv (a, b, c) c = b / a
+ efloor (a, b) truncate to integer, toward -infinity
+ efrexp (a, exp, s) extract exponent and significand
+ eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
+ euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
+ einfin (e) set e to infinity, leaving its sign alone
+ eldexp (a, n, b) multiply by 2**n
+ emov (a, b) b = a
+ emul (a, b, c) c = b * a
+ eneg (e) e = -e
+ eround (a, b) b = nearest integer value to a
+ esub (a, b, c) c = b - a
+ e24toasc (&f, str, n) single to ASCII string, n digits after decimal
+ e53toasc (&d, str, n) double to ASCII string, n digits after decimal
+ e64toasc (&d, str, n) 80-bit long double to ASCII string
+ e113toasc (&d, str, n) 128-bit long double to ASCII string
+ etoasc (e, str, n) e to ASCII string, n digits after decimal
+ etoe24 (e, &f) convert e type to IEEE single precision
+ etoe53 (e, &d) convert e type to IEEE double precision
+ etoe64 (e, &d) convert e type to IEEE long double precision
+ ltoe (&l, e) HOST_WIDE_INT to e type
+ ultoe (&l, e) unsigned HOST_WIDE_INT to e type
+ eisneg (e) 1 if sign bit of e != 0, else 0
+ eisinf (e) 1 if e has maximum exponent (non-IEEE)
+ or is infinite (IEEE)
+ eisnan (e) 1 if e is a NaN
+
+
+ Routines for internal format exploded e-type numbers
+
+ eaddm (ai, bi) add significands, bi = bi + ai
+ ecleaz (ei) ei = 0
+ ecleazs (ei) set ei = 0 but leave its sign alone
+ ecmpm (ai, bi) compare significands, return 1, 0, or -1
+ edivm (ai, bi) divide significands, bi = bi / ai
+ emdnorm (ai,l,s,exp) normalize and round off
+ emovi (a, ai) convert external a to internal ai
+ emovo (ai, a) convert internal ai to external a
+ emovz (ai, bi) bi = ai, low guard word of bi = 0
+ emulm (ai, bi) multiply significands, bi = bi * ai
+ enormlz (ei) left-justify the significand
+ eshdn1 (ai) shift significand and guards down 1 bit
+ eshdn8 (ai) shift down 8 bits
+ eshdn6 (ai) shift down 16 bits
+ eshift (ai, n) shift ai n bits up (or down if n < 0)
+ eshup1 (ai) shift significand and guards up 1 bit
+ eshup8 (ai) shift up 8 bits
+ eshup6 (ai) shift up 16 bits
+ esubm (ai, bi) subtract significands, bi = bi - ai
+ eiisinf (ai) 1 if infinite
+ eiisnan (ai) 1 if a NaN
+ eiisneg (ai) 1 if sign bit of ai != 0, else 0
+ einan (ai) set ai = NaN
+ eiinfin (ai) set ai = infinity
+
+ The result is always normalized and rounded to NI-4 word precision
+ after each arithmetic operation.
+
+ Exception flags are NOT fully supported.
+
+ Signaling NaN's are NOT supported; they are treated the same
+ as quiet NaN's.
+
+ Define INFINITY for support of infinity; otherwise a
+ saturation arithmetic is implemented.
+
+ Define NANS for support of Not-a-Number items; otherwise the
+ arithmetic will never produce a NaN output, and might be confused
+ by a NaN input.
+ If NaN's are supported, the output of `ecmp (a,b)' is -2 if
+ either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
+ may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
+ if in doubt.
+
+ Denormals are always supported here where appropriate (e.g., not
+ for conversion to DEC numbers). */
+
+/* Definitions for error codes that are passed to the common error handling
+ routine mtherr.
+
+ For Digital Equipment PDP-11 and VAX computers, certain
+ IBM systems, and others that use numbers with a 56-bit
+ significand, the symbol DEC should be defined. In this
+ mode, most floating point constants are given as arrays
+ of octal integers to eliminate decimal to binary conversion
+ errors that might be introduced by the compiler.
+
+ For computers, such as IBM PC, that follow the IEEE
+ Standard for Binary Floating Point Arithmetic (ANSI/IEEE
+ Std 754-1985), the symbol IEEE should be defined.
+ These numbers have 53-bit significands. In this mode, constants
+ are provided as arrays of hexadecimal 16 bit integers.
+ The endian-ness of generated values is controlled by
+ REAL_WORDS_BIG_ENDIAN.
+
+ To accommodate other types of computer arithmetic, all
+ constants are also provided in a normal decimal radix
+ which one can hope are correctly converted to a suitable
+ format by the available C language compiler. To invoke
+ this mode, the symbol UNK is defined.
+
+ An important difference among these modes is a predefined
+ set of machine arithmetic constants for each. The numbers
+ MACHEP (the machine roundoff error), MAXNUM (largest number
+ represented), and several other parameters are preset by
+ the configuration symbol. Check the file const.c to
+ ensure that these values are correct for your computer.
+
+ For ANSI C compatibility, define ANSIC equal to 1. Currently
+ this affects only the atan2 function and others that use it. */
-
-/* mconf.h
- *
- * Common include file for math routines
- *
- *
- *
- * SYNOPSIS:
- *
- * #include "mconf.h"
- *
- *
- *
- * DESCRIPTION:
- *
- * This file contains definitions for error codes that are
- * passed to the common error handling routine mtherr
- * (which see).
- *
- * The file also includes a conditional assembly definition
- * for the type of computer arithmetic (Intel IEEE, DEC, Motorola
- * IEEE, or UNKnown).
- *
- * For Digital Equipment PDP-11 and VAX computers, certain
- * IBM systems, and others that use numbers with a 56-bit
- * significand, the symbol DEC should be defined. In this
- * mode, most floating point constants are given as arrays
- * of octal integers to eliminate decimal to binary conversion
- * errors that might be introduced by the compiler.
- *
- * For computers, such as IBM PC, that follow the IEEE
- * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
- * Std 754-1985), the symbol IBMPC or MIEEE should be defined.
- * These numbers have 53-bit significands. In this mode, constants
- * are provided as arrays of hexadecimal 16 bit integers.
- *
- * To accommodate other types of computer arithmetic, all
- * constants are also provided in a normal decimal radix
- * which one can hope are correctly converted to a suitable
- * format by the available C language compiler. To invoke
- * this mode, the symbol UNK is defined.
- *
- * An important difference among these modes is a predefined
- * set of machine arithmetic constants for each. The numbers
- * MACHEP (the machine roundoff error), MAXNUM (largest number
- * represented), and several other parameters are preset by
- * the configuration symbol. Check the file const.c to
- * ensure that these values are correct for your computer.
- *
- * For ANSI C compatibility, define ANSIC equal to 1. Currently
- * this affects only the atan2 function and others that use it.
- */
-\f
/* Constant definitions for math error conditions. */
#define DOMAIN 1 /* argument domain error */
#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 */
-/*include "ehead.h"*/
+#if LONG_DOUBLE_TYPE_SIZE == 128
/* 0.0 */
unsigned EMUSHORT ezero[NE] =
-{
- 0, 0000000, 0000000, 0000000, 0000000, 0000000,};
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
extern unsigned EMUSHORT ezero[];
/* 5.0E-1 */
unsigned EMUSHORT ehalf[NE] =
-{
- 0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
extern unsigned EMUSHORT ehalf[];
/* 1.0E0 */
unsigned EMUSHORT eone[NE] =
-{
- 0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
extern unsigned EMUSHORT eone[];
/* 2.0E0 */
unsigned EMUSHORT etwo[NE] =
-{
- 0, 0000000, 0000000, 0000000, 0100000, 0040000,};
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
extern unsigned EMUSHORT etwo[];
/* 3.2E1 */
unsigned EMUSHORT e32[NE] =
-{
- 0, 0000000, 0000000, 0000000, 0100000, 0040004,};
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
extern unsigned EMUSHORT e32[];
/* 6.93147180559945309417232121458176568075500134360255E-1 */
unsigned EMUSHORT elog2[NE] =
-{
- 0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
+ {0x40f3, 0xf6af, 0x03f2, 0xb398,
+ 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
extern unsigned EMUSHORT elog2[];
/* 1.41421356237309504880168872420969807856967187537695E0 */
unsigned EMUSHORT esqrt2[NE] =
-{
- 0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
+ {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
+ 0x597d, 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] =
-{
+ {0x2902, 0x1cd1, 0x80dc, 0x628b,
0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
extern unsigned EMUSHORT epi[];
-/* 5.7721566490153286060651209008240243104215933593992E-1 */
-unsigned EMUSHORT eeul[NE] =
-{
- 0xd1be, 0xc7a4, 0076660, 0063743, 0111704, 0x3ffe,};
-extern unsigned EMUSHORT eeul[];
-
-/*
-include "ehead.h"
-include "mconf.h"
-*/
-
-
+#else
+/* LONG_DOUBLE_TYPE_SIZE is other than 128 */
+unsigned EMUSHORT ezero[NE] =
+ {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
+unsigned EMUSHORT ehalf[NE] =
+ {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
+unsigned EMUSHORT eone[NE] =
+ {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
+unsigned EMUSHORT etwo[NE] =
+ {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
+unsigned EMUSHORT e32[NE] =
+ {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
+unsigned EMUSHORT elog2[NE] =
+ {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
+unsigned EMUSHORT esqrt2[NE] =
+ {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
+unsigned EMUSHORT epi[NE] =
+ {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
+#endif
/* Control register for rounding precision.
- * This can be set to 80 (if NE=6), 64, 56, 53, or 24 bits.
- */
+ This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits. */
+
int rndprc = NBITS;
extern int rndprc;
-void eaddm (), esubm (), emdnorm (), asctoeg ();
-static void toe24 (), toe53 (), toe64 ();
-void eremain (), einit (), eiremain ();
-int ecmpm (), edivm (), emulm ();
-void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
-void etodec (), todec (), dectoe ();
-
-
-
-
-void
-einit ()
-{
-}
-
-/*
-; Clear out entire external format number.
-;
-; unsigned EMUSHORT x[];
-; eclear (x);
-*/
+/* Clear out entire e-type number X. */
-void
+static void
eclear (x)
register unsigned EMUSHORT *x;
{
*x++ = 0;
}
+/* Move e-type number from A to B. */
-
-/* Move external format number from a to b.
- *
- * emov (a, b);
- */
-
-void
+static void
emov (a, b)
register unsigned EMUSHORT *a, *b;
{
}
-/*
-; Absolute value of external format number
-;
-; EMUSHORT x[NE];
-; eabs (x);
-*/
+/* Absolute value of e-type X. */
-void
+static void
eabs (x)
- unsigned EMUSHORT x[]; /* x is the memory address of a short */
+ unsigned EMUSHORT x[];
{
-
- x[NE - 1] &= 0x7fff; /* sign is top bit of last word of external format */
+ /* sign is top bit of last word of external format */
+ x[NE - 1] &= 0x7fff;
}
+/* Negate the e-type number X. */
-
-
-/*
-; Negate external format number
-;
-; unsigned EMUSHORT x[NE];
-; eneg (x);
-*/
-
-void
+static void
eneg (x)
unsigned EMUSHORT x[];
{
x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
}
+/* Return 1 if sign bit of e-type number X is nonzero, else zero. */
-
-/* Return 1 if external format number is negative,
- * else return zero.
- */
-int
+static int
eisneg (x)
unsigned EMUSHORT x[];
{
return (0);
}
+/* Return 1 if e-type number X is infinity, else return zero. */
-/* Return 1 if external format number has maximum possible exponent,
- * else return zero.
- */
-int
+static int
eisinf (x)
unsigned EMUSHORT x[];
{
+#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. */
-/*
-; Fill entire number, including exponent and significand, with
-; largest possible number. These programs implement a saturation
-; value that is an ordinary, legal number. A special value
-; "infinity" may also be implemented; this would require tests
-; for that value and implementation of special rules for arithmetic
-; operations involving inifinity.
-*/
+static int
+eisnan (x)
+ unsigned EMUSHORT x[];
+{
+#ifdef NANS
+ int i;
-void
+ /* NaN has maximum exponent */
+ if ((x[NE - 1] & 0x7fff) != 0x7fff)
+ return (0);
+ /* ... and non-zero significand field. */
+ for (i = 0; i < NE - 1; i++)
+ {
+ if (*x++ != 0)
+ return (1);
+ }
+#endif
+
+ return (0);
+}
+
+/* Fill e-type number X with infinity pattern (IEEE)
+ or largest possible number (non-IEEE). */
+
+static void
einfin (x)
register unsigned EMUSHORT *x;
{
*x |= 32766;
if (rndprc < NBITS)
{
+ if (rndprc == 113)
+ {
+ *(x - 9) = 0;
+ *(x - 8) = 0;
+ }
if (rndprc == 64)
{
*(x - 5) = 0;
#endif
}
+/* Output an e-type NaN.
+ This generates Intel's quiet NaN pattern for extended real.
+ The exponent is 7fff, the leading mantissa word is c000. */
+
+static void
+enan (x, sign)
+ register unsigned EMUSHORT *x;
+ int sign;
+{
+ register int i;
+
+ for (i = 0; i < NE - 2; i++)
+ *x++ = 0;
+ *x++ = 0xc000;
+ *x = (sign << 15) | 0x7fff;
+}
+/* Move in an e-type number A, converting it to exploded e-type B. */
-/* Move in external format number,
- * converting it to internal format.
- */
-void
+static void
emovi (a, b)
unsigned EMUSHORT *a, *b;
{
#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 */
*q = 0;
}
+/* Move out exploded e-type number A, converting it to e type B. */
-/* Move internal format number out,
- * converting it to external format.
- */
-void
+static void
emovo (a, b)
unsigned EMUSHORT *a, *b;
{
register unsigned EMUSHORT *p, *q;
unsigned EMUSHORT i;
+ int j;
p = a;
q = b + (NE - 1); /* point to output exponent */
#ifdef INFINITY
if (*(p - 1) == 0x7fff)
{
+#ifdef NANS
+ if (eiisnan (a))
+ {
+ enan (b, eiisneg (a));
+ return;
+ }
+#endif
einfin (b);
- return;
+ return;
}
#endif
/* skip over guard word */
++p;
/* move the significand */
- for (i = 0; i < NE - 1; i++)
+ for (j = 0; j < NE - 1; j++)
*q-- = *p++;
}
+/* Clear out exploded e-type number XI. */
-
-
-/* Clear out internal format number.
- */
-
-void
+static void
ecleaz (xi)
register unsigned EMUSHORT *xi;
{
*xi++ = 0;
}
+/* Clear out exploded e-type XI, but don't touch the sign. */
-/* same, but don't touch the sign. */
-
-void
+static void
ecleazs (xi)
register unsigned EMUSHORT *xi;
{
*xi++ = 0;
}
+/* Move exploded e-type number from A to B. */
-
-/* Move internal format number from a to b.
- */
-void
+static void
emovz (a, b)
register unsigned EMUSHORT *a, *b;
{
*b = 0;
}
+/* Generate exploded e-type NaN.
+ The explicit pattern for this is maximum exponent and
+ top two significant bits set. */
-/*
-; Compare significands of numbers in internal format.
-; Guard words are included in the comparison.
-;
-; unsigned EMUSHORT a[NI], b[NI];
-; cmpm (a, b);
-;
-; for the significands:
-; returns +1 if a > b
-; 0 if a == b
-; -1 if a < b
-*/
-int
+static void
+einan (x)
+ unsigned EMUSHORT x[];
+{
+
+ ecleaz (x);
+ x[E] = 0x7fff;
+ x[M + 1] = 0xc000;
+}
+
+/* Return nonzero if exploded e-type X is a NaN. */
+
+static int
+eiisnan (x)
+ unsigned EMUSHORT x[];
+{
+ int i;
+
+ if ((x[E] & 0x7fff) == 0x7fff)
+ {
+ for (i = M + 1; i < NI; i++)
+ {
+ if (x[i] != 0)
+ return (1);
+ }
+ }
+ return (0);
+}
+
+/* Return nonzero if sign of exploded e-type X is nonzero. */
+
+static int
+eiisneg (x)
+ unsigned EMUSHORT x[];
+{
+
+ return x[0] != 0;
+}
+
+/* Fill exploded e-type X with infinity pattern.
+ This has maximum exponent and significand all zeros. */
+
+static void
+eiinfin (x)
+ unsigned EMUSHORT x[];
+{
+
+ ecleaz (x);
+ x[E] = 0x7fff;
+}
+
+/* Return nonzero if exploded e-type X is infinite. */
+
+static int
+eiisinf (x)
+ unsigned EMUSHORT x[];
+{
+
+#ifdef NANS
+ if (eiisnan (x))
+ return (0);
+#endif
+ if ((x[E] & 0x7fff) == 0x7fff)
+ return (1);
+ return (0);
+}
+
+
+/* Compare significands of numbers in internal exploded e-type format.
+ Guard words are included in the comparison.
+
+ Returns +1 if a > b
+ 0 if a == b
+ -1 if a < b */
+
+static int
ecmpm (a, b)
register unsigned EMUSHORT *a, *b;
{
return (-1);
}
+/* Shift significand of exploded e-type X down by 1 bit. */
-/*
-; Shift significand down by 1 bit
-*/
-
-void
+static void
eshdn1 (x)
register unsigned EMUSHORT *x;
{
}
}
+/* Shift significand of exploded e-type X up by 1 bit. */
-
-/*
-; Shift significand up by 1 bit
-*/
-
-void
+static void
eshup1 (x)
register unsigned EMUSHORT *x;
{
}
+/* Shift significand of exploded e-type X down by 8 bits. */
-/*
-; Shift significand down by 8 bits
-*/
-
-void
+static void
eshdn8 (x)
register unsigned EMUSHORT *x;
{
}
}
-/*
-; Shift significand up by 8 bits
-*/
+/* Shift significand of exploded e-type X up by 8 bits. */
-void
+static void
eshup8 (x)
register unsigned EMUSHORT *x;
{
}
}
-/*
-; Shift significand up by 16 bits
-*/
+/* Shift significand of exploded e-type X up by 16 bits. */
-void
+static void
eshup6 (x)
register unsigned EMUSHORT *x;
{
*p = 0;
}
-/*
-; Shift significand down by 16 bits
-*/
+/* Shift significand of exploded e-type X down by 16 bits. */
-void
+static void
eshdn6 (x)
register unsigned EMUSHORT *x;
{
*(--p) = 0;
}
-\f
-/*
-; Add significands
-; x + y replaces y
-*/
-void
+/* Add significands of exploded e-type X and Y. X + Y replaces Y. */
+
+static void
eaddm (x, y)
unsigned EMUSHORT *x, *y;
{
}
}
-/*
-; Subtract significands
-; y - x replaces y
-*/
+/* Subtract significands of exploded e-type X and Y. Y - X replaces Y. */
-void
+static void
esubm (x, y)
unsigned EMUSHORT *x, *y;
{
}
-/* Divide significands */
-
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[];
*p++ = 0;
}
- /* Use faster compare and subtraction if denominator
- * has only 15 bits of significance.
- */
+ /* Use faster compare and subtraction if denominator has only 15 bits of
+ significance. */
+
p = &den[M + 2];
if (*p++ == 0)
{
goto divdon;
}
- /* The number of quotient bits to calculate is
- * NBITS + 1 scaling guard bit + 1 roundoff bit.
- */
+ /* The number of quotient bits to calculate is NBITS + 1 scaling guard
+ bit + 1 roundoff bit. */
+
fulldiv:
p = &equot[NI - 2];
/* Multiply significands */
+
int
emulm (a, b)
unsigned EMUSHORT a[], b[];
p = &a[NI - 2];
k = NBITS;
- while (*p == 0) /* significand is not supposed to be all zero */
+ while (*p == 0) /* significand is not supposed to be zero */
{
eshdn6 (a);
k -= 16;
return (j);
}
+#else
+/* Radix 65536 versions of multiply and divide. */
-/*
- * Normalize and round off.
- *
- * The internal format number to be rounded is "s".
- * Input "lost" indicates whether or not the number is exact.
- * This is the so-called sticky bit.
- *
- * Input "subflg" indicates whether the number was obtained
- * by a subtraction operation. In that case if lost is nonzero
- * then the number is slightly smaller than indicated.
- *
- * Input "exp" is the biased exponent, which may be negative.
- * the exponent field of "s" is ignored but is replaced by
- * "exp" as adjusted by normalization and rounding.
- *
- * Input "rcntrl" is the rounding control.
- */
-
-static int rlast = -1;
-static int rw = 0;
-static unsigned EMUSHORT rmsk = 0;
-static unsigned EMUSHORT rmbit = 0;
-static unsigned EMUSHORT rebit = 0;
-static int re = 0;
-static unsigned EMUSHORT rbit[NI];
+/* Multiply significand of e-type number B
+ by 16-bit quantity A, return e-type result to C. */
-void
-emdnorm (s, lost, subflg, exp, rcntrl)
- unsigned EMUSHORT s[];
- int lost;
- int subflg;
- EMULONG exp;
- int rcntrl;
+static void
+m16m (a, b, c)
+ unsigned int a;
+ unsigned EMUSHORT b[], c[];
{
- int i, j;
- unsigned EMUSHORT r;
+ register unsigned EMUSHORT *pp;
+ register unsigned EMULONG carry;
+ unsigned EMUSHORT *ps;
+ unsigned EMUSHORT p[NI];
+ unsigned EMULONG aa, m;
+ int i;
- /* Normalize */
- j = enormlz (s);
+ aa = a;
+ pp = &p[NI-2];
+ *pp++ = 0;
+ *pp = 0;
+ ps = &b[NI-1];
- /* a blank significand could mean either zero or infinity. */
-#ifndef INFINITY
- if (j > NBITS)
+ for (i=M+1; i<NI; i++)
{
- ecleazs (s);
- return;
+ if (*ps == 0)
+ {
+ --ps;
+ --pp;
+ *(pp-1) = 0;
+ }
+ else
+ {
+ m = (unsigned EMULONG) aa * *ps--;
+ carry = (m & 0xffff) + *pp;
+ *pp-- = (unsigned EMUSHORT)carry;
+ carry = (carry >> 16) + (m >> 16) + *pp;
+ *pp = (unsigned EMUSHORT)carry;
+ *(pp-1) = carry >> 16;
+ }
}
-#endif
- exp -= j;
-#ifndef INFINITY
- if (exp >= 32767L)
- goto overf;
-#else
+ for (i=M; i<NI; i++)
+ c[i] = p[i];
+}
+
+/* Divide significands of exploded e-types NUM / DEN. Neither the
+ numerator NUM nor the denominator DEN is permitted to have its high guard
+ word nonzero. */
+
+static int
+edivm (den, num)
+ unsigned EMUSHORT den[], num[];
+{
+ int i;
+ register unsigned EMUSHORT *p;
+ unsigned EMULONG tnum;
+ unsigned EMUSHORT j, tdenm, tquot;
+ unsigned EMUSHORT tprod[NI+1];
+
+ p = &equot[0];
+ *p++ = num[0];
+ *p++ = num[1];
+
+ for (i=M; i<NI; i++)
+ {
+ *p++ = 0;
+ }
+ eshdn1 (num);
+ tdenm = den[M+1];
+ for (i=M; i<NI; i++)
+ {
+ /* Find trial quotient digit (the radix is 65536). */
+ tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
+
+ /* Do not execute the divide instruction if it will overflow. */
+ if ((tdenm * 0xffffL) < tnum)
+ tquot = 0xffff;
+ else
+ tquot = tnum / tdenm;
+ /* Multiply denominator by trial quotient digit. */
+ m16m ((unsigned int)tquot, den, tprod);
+ /* The quotient digit may have been overestimated. */
+ if (ecmpm (tprod, num) > 0)
+ {
+ tquot -= 1;
+ esubm (den, tprod);
+ if (ecmpm (tprod, num) > 0)
+ {
+ tquot -= 1;
+ esubm (den, tprod);
+ }
+ }
+ esubm (tprod, num);
+ equot[i] = tquot;
+ eshup6(num);
+ }
+ /* test for nonzero remainder after roundoff bit */
+ p = &num[M];
+ j = 0;
+ for (i=M; i<NI; i++)
+ {
+ j |= *p++;
+ }
+ if (j)
+ j = 1;
+
+ for (i=0; i<NI; i++)
+ num[i] = equot[i];
+
+ return ((int)j);
+}
+
+/* Multiply significands of exploded e-type A and B, result in B. */
+
+static int
+emulm (a, b)
+ unsigned EMUSHORT a[], b[];
+{
+ unsigned EMUSHORT *p, *q;
+ unsigned EMUSHORT pprod[NI];
+ unsigned EMUSHORT j;
+ int i;
+
+ equot[0] = b[0];
+ equot[1] = b[1];
+ for (i=M; i<NI; i++)
+ equot[i] = 0;
+
+ j = 0;
+ p = &a[NI-1];
+ q = &equot[NI-1];
+ for (i=M+1; i<NI; i++)
+ {
+ if (*p == 0)
+ {
+ --p;
+ }
+ else
+ {
+ m16m ((unsigned int) *p--, b, pprod);
+ eaddm(pprod, equot);
+ }
+ j |= *q;
+ eshdn6(equot);
+ }
+
+ for (i=0; i<NI; i++)
+ b[i] = equot[i];
+
+ /* return flag for lost nonzero bits */
+ return ((int)j);
+}
+#endif
+
+
+/* Normalize and round off.
+
+ The internal format number to be rounded is S.
+ Input LOST is 0 if the value is exact. This is the so-called sticky bit.
+
+ Input SUBFLG indicates whether the number was obtained
+ by a subtraction operation. In that case if LOST is nonzero
+ then the number is slightly smaller than indicated.
+
+ Input EXP is the biased exponent, which may be negative.
+ the exponent field of S is ignored but is replaced by
+ EXP as adjusted by normalization and rounding.
+
+ Input RCNTRL is the rounding control. If it is nonzero, the
+ returned value will be rounded to RNDPRC bits.
+
+ For future reference: In order for emdnorm to round off denormal
+ significands at the right point, the input exponent must be
+ adjusted to be the actual value it would have after conversion to
+ the final floating point type. This adjustment has been
+ implemented for all type conversions (etoe53, etc.) and decimal
+ conversions, but not for the arithmetic functions (eadd, etc.).
+ Data types having standard 15-bit exponents are not affected by
+ this, but SFmode and DFmode are affected. For example, ediv with
+ rndprc = 24 will not round correctly to 24-bit precision if the
+ result is denormal. */
+
+static int rlast = -1;
+static int rw = 0;
+static unsigned EMUSHORT rmsk = 0;
+static unsigned EMUSHORT rmbit = 0;
+static unsigned EMUSHORT rebit = 0;
+static int re = 0;
+static unsigned EMUSHORT rbit[NI];
+
+static void
+emdnorm (s, lost, subflg, exp, rcntrl)
+ unsigned EMUSHORT s[];
+ int lost;
+ int subflg;
+ EMULONG exp;
+ int rcntrl;
+{
+ int i, j;
+ unsigned EMUSHORT r;
+
+ /* Normalize */
+ j = enormlz (s);
+
+ /* a blank significand could mean either zero or infinity. */
+#ifndef INFINITY
+ if (j > NBITS)
+ {
+ ecleazs (s);
+ return;
+ }
+#endif
+ exp -= j;
+#ifndef INFINITY
+ if (exp >= 32767L)
+ goto overf;
+#else
if ((j > NBITS) && (exp < 32767))
{
ecleazs (s);
return;
}
}
- /* Round off, unless told not to by rcntrl. */
+ /* Round off, unless told not to by rcntrl. */
if (rcntrl == 0)
goto mdfin;
- /* Set up rounding parameters if the control register changed. */
+ /* Set up rounding parameters if the control register changed. */
if (rndprc != rlast)
{
ecleaz (rbit);
rw = NI - 1; /* low guard word */
rmsk = 0xffff;
rmbit = 0x8000;
- rbit[rw - 1] = 1;
- re = NI - 2;
+ 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;
- rbit[rw - 1] = 1;
re = rw - 1;
rebit = 1;
break;
- /* For DEC arithmetic */
+ /* For DEC or IBM arithmetic */
case 56:
rw = 6;
rmsk = 0xff;
rmbit = 0x80;
- rbit[rw] = 0x100;
- re = rw;
rebit = 0x100;
+ re = rw;
break;
case 53:
rw = 6;
rmsk = 0x7ff;
rmbit = 0x0400;
- rbit[rw] = 0x800;
- re = rw;
rebit = 0x800;
+ re = rw;
break;
case 24:
rw = 4;
rmsk = 0xff;
rmbit = 0x80;
- rbit[rw] = 0x100;
- re = rw;
rebit = 0x100;
+ re = rw;
break;
}
+ rbit[re] = rebit;
rlast = rndprc;
}
- if (rndprc >= 64)
+ /* Shift down 1 temporarily if the data structure has an implied
+ most significant bit and the number is denormal.
+ Intel long double denormals also lose one bit of precision. */
+ if ((exp <= 0) && (rndprc != NBITS)
+ && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
{
- r = s[rw] & rmsk;
- if (rndprc == 64)
- {
- i = rw + 1;
- while (i < NI)
- {
- if (s[i])
- r |= 1;
- s[i] = 0;
- ++i;
- }
- }
+ lost |= s[NI - 1] & 1;
+ eshdn1 (s);
}
- else
+ /* Clear out all bits below the rounding bit,
+ remembering in r if any were nonzero. */
+ r = s[rw] & rmsk;
+ if (rndprc < NBITS)
{
- if (exp <= 0)
- eshdn1 (s);
- r = s[rw] & rmsk;
- /* These tests assume NI = 8 */
i = rw + 1;
while (i < NI)
{
s[i] = 0;
++i;
}
- /*
- if (rndprc == 24)
- {
- if (s[5] || s[6])
- r |= 1;
- s[5] = 0;
- s[6] = 0;
- }
- */
- s[rw] &= ~rmsk;
}
+ s[rw] &= ~rmsk;
if ((r & rmbit) != 0)
{
if (r == rmbit)
eaddm (rbit, s);
}
mddone:
- if ((rndprc < 64) && (exp <= 0))
+/* Undo the temporary shift for denormal values. */
+ if ((exp <= 0) && (rndprc != NBITS)
+ && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
{
eshup1 (s);
}
for (i = M + 1; i < NI - 1; i++)
s[i] = 0xffff;
s[NI - 1] = 0;
- if (rndprc < 64)
+ if ((rndprc < 64) || (rndprc == 113))
{
s[rw] &= ~rmsk;
if (rndprc == 24)
s[1] = (unsigned EMUSHORT) exp;
}
-
-
-/*
-; Subtract external format numbers.
-;
-; unsigned EMUSHORT a[NE], b[NE], c[NE];
-; esub (a, b, c); c = b - a
-*/
+/* Subtract. C = B - A, all e type numbers. */
static int subflg = 0;
-void
+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. */
-/*
-; Add.
-;
-; unsigned EMUSHORT a[NE], b[NE], c[NE];
-; eadd (a, b, c); c = b + a
-*/
-void
+static void
eadd (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
+#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);
}
-void
+/* Arithmetic common to both addition and subtraction. */
+
+static void
eadd1 (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
return;
}
/* if same sign, result is double */
- /* double denomalized tiny number */
+ /* double denormalized tiny number */
if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
{
eshup1 (bi);
{
if (bi[j] != 0)
{
- /* This could overflow, but let emovo take care of that. */
ltb += 1;
+ if (ltb >= 0x7fff)
+ {
+ eclear (c);
+ if (ai[0] != 0)
+ eneg (c);
+ einfin (c);
+ return;
+ }
break;
}
}
emovo (bi, c);
}
+/* Divide: C = B/A, all e type. */
-
-/*
-; Divide.
-;
-; unsigned EMUSHORT a[NE], b[NE], c[NE];
-; ediv (a, b, c); c = b / a
-*/
-void
+static void
ediv (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
unsigned EMUSHORT ai[NI], bi[NI];
- int i;
+ 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))
{
- if (eisneg (a) ^ eisneg (b))
- *(c + (NE - 1)) = 0x8000;
- else
- *(c + (NE - 1)) = 0;
einfin (c);
- return;
+ goto divsign;
}
+/* Anything else over infinity is zero. */
if (eisinf (a))
{
eclear (c);
- return;
+ goto divsign;
}
#endif
emovi (a, ai);
lta = ai[E];
ltb = bi[E];
if (bi[E] == 0)
- { /* See if numerator is zero. */
+ { /* See if numerator is zero. */
for (i = 1; i < NI - 1; i++)
{
if (bi[i] != 0)
}
}
eclear (c);
- return;
+ goto divsign;
}
dnzro1:
goto dnzro2;
}
}
- if (ai[0] == bi[0])
- *(c + (NE - 1)) = 0;
- else
- *(c + (NE - 1)) = 0x8000;
+/* Divide by zero is not an invalid operation.
+ It is a divide-by-zero operation! */
einfin (c);
mtherr ("ediv", SING);
- return;
+ goto divsign;
}
dnzro2:
/* 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);
-}
+ divsign:
+
+ if (sign
+#ifndef IEEE
+ && (ecmp (c, ezero) != 0)
+#endif
+ )
+ *(c+(NE-1)) |= 0x8000;
+ else
+ *(c+(NE-1)) &= ~0x8000;
+}
+/* Multiply e-types A and B, return e-type product C. */
-/*
-; Multiply.
-;
-; unsigned EMUSHORT a[NE], b[NE], c[NE];
-; emul (a, b, c); c = b * a
-*/
-void
+static void
emul (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
unsigned EMUSHORT ai[NI], bi[NI];
- int i, j;
+ 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))
{
- if (eisneg (a) ^ eisneg (b))
- *(c + (NE - 1)) = 0x8000;
- else
- *(c + (NE - 1)) = 0;
einfin (c);
- return;
+ goto mulsign;
}
#endif
emovi (a, ai);
}
}
eclear (c);
- return;
+ goto mulsign;
}
mnzer1:
}
}
eclear (c);
- return;
+ goto mulsign;
}
mnzer2:
/* 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);
-}
+ mulsign:
+ if (sign
+#ifndef IEEE
+ && (ecmp (c, ezero) != 0)
+#endif
+ )
+ *(c+(NE-1)) |= 0x8000;
+ else
+ *(c+(NE-1)) &= ~0x8000;
+}
+/* Convert double precision PE to e-type Y. */
-/*
-; Convert IEEE double precision to e type
-; double d;
-; unsigned EMUSHORT x[N+2];
-; e53toe (&d, x);
-*/
-void
-e53toe (e, y)
- unsigned EMUSHORT *e, *y;
+static void
+e53toe (pe, y)
+ unsigned EMUSHORT *pe, *y;
{
#ifdef DEC
- dectoe (e, y); /* see etodec.c */
+ dectoe (pe, y);
#else
+#ifdef IBM
+
+ ibmtoe (pe, y, DFmode);
+#else
register unsigned EMUSHORT r;
- register unsigned EMUSHORT *p;
+ register unsigned EMUSHORT *e, *p;
unsigned EMUSHORT yy[NI];
int denorm, k;
+ e = pe;
denorm = 0; /* flag if denormalized number */
ecleaz (yy);
-#ifdef IBMPC
- e += 3;
-#endif
+ if (! REAL_WORDS_BIG_ENDIAN)
+ e += 3;
r = *e;
yy[0] = 0;
if (r & 0x8000)
#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 (r & 0x8000)
+ if (yy[0])
eneg (y);
return;
}
-#endif
+#endif /* INFINITY */
r >>= 4;
/* If zero exponent, then the significand is denormalized.
- * So, take back the understood high significand bit. */
+ So take back the understood high significand bit. */
+
if (r == 0)
{
denorm = 1;
r += EXONE - 01777;
yy[E] = r;
p = &yy[M + 1];
-#ifdef IBMPC
- *p++ = *(--e);
- *p++ = *(--e);
- *p++ = *(--e);
-#endif
-#ifdef MIEEE
- ++e;
- *p++ = *e++;
- *p++ = *e++;
- *p++ = *e++;
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
+ {
+ *p++ = *(--e);
+ *p++ = *(--e);
+ *p++ = *(--e);
+ }
+ else
+ {
+ ++e;
+ *p++ = *e++;
+ *p++ = *e++;
+ *p++ = *e++;
+ }
#endif
eshift (yy, -5);
if (denorm)
yy[E] -= (unsigned EMUSHORT) (k - 1);
}
emovo (yy, y);
+#endif /* not IBM */
#endif /* not DEC */
}
-void
-e64toe (e, y)
- unsigned EMUSHORT *e, *y;
+/* Convert double extended precision float PE to e type Y. */
+
+static void
+e64toe (pe, y)
+ unsigned EMUSHORT *pe, *y;
{
unsigned EMUSHORT yy[NI];
- unsigned EMUSHORT *p, *q;
+ unsigned EMUSHORT *e, *p, *q;
int i;
+ e = pe;
p = yy;
for (i = 0; i < NE - 5; i++)
*p++ = 0;
-#ifdef IBMPC
- for (i = 0; i < 5; i++)
- *p++ = *e++;
-#endif
+/* This precision is not ordinarily supported on DEC or IBM. */
#ifdef DEC
for (i = 0; i < 5; i++)
*p++ = *e++;
#endif
-#ifdef MIEEE
+#ifdef IBM
p = &yy[0] + (NE - 1);
*p-- = *e++;
++e;
- for (i = 0; i < 4; i++)
+ for (i = 0; i < 5; i++)
*p-- = *e++;
#endif
- p = yy;
- q = y;
-#ifdef INFINITY
- if (*p == 0x7fff)
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
{
- einfin (y);
- if (*p & 0x8000)
- eneg (y);
- return;
- }
-#endif
- for (i = 0; i < NE; i++)
+ for (i = 0; i < 5; i++)
+ *p++ = *e++;
+
+ /* For denormal long double Intel format, shift significand up one
+ -- but only if the top significand bit is zero. A top bit of 1
+ is "pseudodenormal" when the exponent is zero. */
+ if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
+ {
+ unsigned EMUSHORT temp[NI];
+
+ emovi(yy, temp);
+ eshup1(temp);
+ emovo(temp,y);
+ return;
+ }
+ }
+ else
+ {
+ p = &yy[0] + (NE - 1);
+#ifdef ARM_EXTENDED_IEEE_FORMAT
+ /* For ARMs, the exponent is in the lowest 15 bits of the word. */
+ *p-- = (e[0] & 0x8000) | (e[1] & 0x7ffff);
+ e += 2;
+#else
+ *p-- = *e++;
+ ++e;
+#endif
+ for (i = 0; i < 4; i++)
+ *p-- = *e++;
+ }
+#endif
+#ifdef INFINITY
+ /* Point to the exponent field and check max exponent cases. */
+ p = &yy[NE - 1];
+ if ((*p & 0x7fff) == 0x7fff)
+ {
+#ifdef NANS
+ if (! REAL_WORDS_BIG_ENDIAN)
+ {
+ for (i = 0; i < 4; i++)
+ {
+ if ((i != 3 && pe[i] != 0)
+ /* Anything but 0x8000 here, including 0, is a NaN. */
+ || (i == 3 && pe[i] != 0x8000))
+ {
+ enan (y, (*p & 0x8000) != 0);
+ return;
+ }
+ }
+ }
+ else
+ {
+#ifdef ARM_EXTENDED_IEEE_FORMAT
+ for (i = 2; i <= 5; i++)
+ {
+ if (pe[i] != 0)
+ {
+ enan (y, (*p & 0x8000) != 0);
+ return;
+ }
+ }
+#else /* not ARM */
+ /* In Motorola extended precision format, the most significant
+ bit of an infinity mantissa could be either 1 or 0. It is
+ the lower order bits that tell whether the value is a NaN. */
+ if ((pe[2] & 0x7fff) != 0)
+ goto bigend_nan;
+
+ for (i = 3; i <= 5; i++)
+ {
+ if (pe[i] != 0)
+ {
+bigend_nan:
+ enan (y, (*p & 0x8000) != 0);
+ return;
+ }
+ }
+#endif /* not ARM */
+ }
+#endif /* NANS */
+ eclear (y);
+ einfin (y);
+ if (*p & 0x8000)
+ eneg (y);
+ return;
+ }
+#endif /* INFINITY */
+ p = yy;
+ q = y;
+ for (i = 0; i < NE; i++)
*q++ = *p++;
}
+/* Convert 128-bit long double precision float PE to e type Y. */
-/*
-; Convert IEEE single precision to e type
-; float d;
-; unsigned EMUSHORT x[N+2];
-; dtox (&d, x);
-*/
-void
-e24toe (e, y)
- unsigned EMUSHORT *e, *y;
+static void
+e113toe (pe, y)
+ unsigned EMUSHORT *pe, *y;
{
register unsigned EMUSHORT r;
- register unsigned EMUSHORT *p;
+ unsigned EMUSHORT *e, *p;
+ unsigned EMUSHORT yy[NI];
+ int denorm, i;
+
+ e = pe;
+ denorm = 0;
+ ecleaz (yy);
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
+ e += 7;
+#endif
+ r = *e;
+ yy[0] = 0;
+ if (r & 0x8000)
+ yy[0] = 0xffff;
+ r &= 0x7fff;
+#ifdef INFINITY
+ if (r == 0x7fff)
+ {
+#ifdef NANS
+ if (! REAL_WORDS_BIG_ENDIAN)
+ {
+ for (i = 0; i < 7; i++)
+ {
+ if (pe[i] != 0)
+ {
+ enan (y, yy[0] != 0);
+ return;
+ }
+ }
+ }
+ else
+ {
+ for (i = 1; i < 8; i++)
+ {
+ if (pe[i] != 0)
+ {
+ enan (y, yy[0] != 0);
+ return;
+ }
+ }
+ }
+#endif /* NANS */
+ eclear (y);
+ einfin (y);
+ if (yy[0])
+ eneg (y);
+ return;
+ }
+#endif /* INFINITY */
+ yy[E] = r;
+ p = &yy[M + 1];
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
+ {
+ for (i = 0; i < 7; i++)
+ *p++ = *(--e);
+ }
+ else
+ {
+ ++e;
+ for (i = 0; i < 7; i++)
+ *p++ = *e++;
+ }
+#endif
+/* If denormal, remove the implied bit; else shift down 1. */
+ if (r == 0)
+ {
+ yy[M] = 0;
+ }
+ else
+ {
+ yy[M] = 1;
+ eshift (yy, -1);
+ }
+ emovo (yy, y);
+}
+
+/* Convert single precision float PE to e type Y. */
+
+static void
+e24toe (pe, y)
+ unsigned EMUSHORT *pe, *y;
+{
+#ifdef IBM
+
+ ibmtoe (pe, y, SFmode);
+
+#else
+ register unsigned EMUSHORT r;
+ register unsigned EMUSHORT *e, *p;
unsigned EMUSHORT yy[NI];
int denorm, k;
+ e = pe;
denorm = 0; /* flag if denormalized number */
ecleaz (yy);
-#ifdef IBMPC
- e += 1;
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
+ e += 1;
#endif
#ifdef DEC
e += 1;
#ifdef INFINITY
if (r == 0x7f80)
{
+#ifdef NANS
+ if (REAL_WORDS_BIG_ENDIAN)
+ {
+ if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
+ {
+ enan (y, yy[0] != 0);
+ return;
+ }
+ }
+ else
+ {
+ if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
+ {
+ enan (y, yy[0] != 0);
+ return;
+ }
+ }
+#endif /* NANS */
+ eclear (y);
einfin (y);
- if (r & 0x8000)
+ if (yy[0])
eneg (y);
return;
}
-#endif
+#endif /* INFINITY */
r >>= 7;
/* If zero exponent, then the significand is denormalized.
- * So, take back the understood high significand bit. */
+ So take back the understood high significand bit. */
if (r == 0)
{
denorm = 1;
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++;
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
+ *p++ = *(--e);
+ else
+ {
+ ++e;
+ *p++ = *e++;
+ }
#endif
eshift (yy, -8);
if (denorm)
yy[E] -= (unsigned EMUSHORT) (k - 1);
}
emovo (yy, y);
+#endif /* not IBM */
}
+/* Convert e-type X to IEEE 128-bit long double format E. */
-void
+static void
+etoe113 (x, e)
+ unsigned EMUSHORT *x, *e;
+{
+ unsigned EMUSHORT xi[NI];
+ EMULONG exp;
+ int rndsav;
+
+#ifdef NANS
+ if (eisnan (x))
+ {
+ make_nan (e, eisneg (x), TFmode);
+ return;
+ }
+#endif
+ emovi (x, xi);
+ exp = (EMULONG) xi[E];
+#ifdef INFINITY
+ if (eisinf (x))
+ goto nonorm;
+#endif
+ /* round off to nearest or even */
+ rndsav = rndprc;
+ rndprc = 113;
+ emdnorm (xi, 0, 0, exp, 64);
+ rndprc = rndsav;
+ nonorm:
+ toe113 (xi, e);
+}
+
+/* Convert exploded e-type X, that has already been rounded to
+ 113-bit precision, to IEEE 128-bit long double format Y. */
+
+static void
+toe113 (a, b)
+ unsigned EMUSHORT *a, *b;
+{
+ register unsigned EMUSHORT *p, *q;
+ unsigned EMUSHORT i;
+
+#ifdef NANS
+ if (eiisnan (a))
+ {
+ make_nan (b, eiisneg (a), TFmode);
+ return;
+ }
+#endif
+ p = a;
+ if (REAL_WORDS_BIG_ENDIAN)
+ q = b;
+ else
+ q = b + 7; /* point to output exponent */
+
+ /* If not denormal, delete the implied bit. */
+ if (a[E] != 0)
+ {
+ eshup1 (a);
+ }
+ /* combine sign and exponent */
+ i = *p++;
+ if (REAL_WORDS_BIG_ENDIAN)
+ {
+ if (i)
+ *q++ = *p++ | 0x8000;
+ else
+ *q++ = *p++;
+ }
+ else
+ {
+ if (i)
+ *q-- = *p++ | 0x8000;
+ else
+ *q-- = *p++;
+ }
+ /* skip over guard word */
+ ++p;
+ /* move the significand */
+ if (REAL_WORDS_BIG_ENDIAN)
+ {
+ for (i = 0; i < 7; i++)
+ *q++ = *p++;
+ }
+ else
+ {
+ for (i = 0; i < 7; i++)
+ *q-- = *p++;
+ }
+}
+
+/* Convert e-type X to IEEE double extended format E. */
+
+static void
etoe64 (x, e)
unsigned EMUSHORT *x, *e;
{
EMULONG exp;
int rndsav;
+#ifdef NANS
+ if (eisnan (x))
+ {
+ make_nan (e, eisneg (x), XFmode);
+ return;
+ }
+#endif
emovi (x, xi);
/* adjust exponent for offset */
exp = (EMULONG) xi[E];
toe64 (xi, e);
}
-/* move out internal format to ieee long double */
+/* Convert exploded e-type X, that has already been rounded to
+ 64-bit precision, to IEEE double extended format Y. */
+
static void
toe64 (a, b)
unsigned EMUSHORT *a, *b;
register unsigned EMUSHORT *p, *q;
unsigned EMUSHORT i;
+#ifdef NANS
+ if (eiisnan (a))
+ {
+ make_nan (b, eiisneg (a), XFmode);
+ return;
+ }
+#endif
+ /* Shift denormal long double Intel format significand down one bit. */
+ if ((a[E] == 0) && ! REAL_WORDS_BIG_ENDIAN)
+ eshdn1 (a);
p = a;
-#ifdef MIEEE
+#ifdef IBM
q = b;
-#else
- q = b + 4; /* point to output exponent */
+#endif
+#ifdef DEC
+ q = b + 4;
+#endif
+#ifdef IEEE
+ if (REAL_WORDS_BIG_ENDIAN)
+ q = b;
+ else
+ {
+ q = b + 4; /* point to output exponent */
#if LONG_DOUBLE_TYPE_SIZE == 96
- /* Clear the last two bytes of 12-byte Intel format */
- *(q+1) = 0;
+ /* Clear the last two bytes of 12-byte Intel format */
+ *(q+1) = 0;
#endif
+ }
#endif
/* combine sign and exponent */
i = *p++;
-#ifdef MIEEE
+#ifdef IBM
if (i)
*q++ = *p++ | 0x8000;
else
*q++ = *p++;
*q++ = 0;
-#else
+#endif
+#ifdef DEC
if (i)
*q-- = *p++ | 0x8000;
else
*q-- = *p++;
#endif
+#ifdef IEEE
+ if (REAL_WORDS_BIG_ENDIAN)
+ {
+#ifdef ARM_EXTENDED_IEEE_FORMAT
+ /* The exponent is in the lowest 15 bits of the first word. */
+ *q++ = i ? 0x8000 : 0;
+ *q++ = *p++;
+#else
+ if (i)
+ *q++ = *p++ | 0x8000;
+ else
+ *q++ = *p++;
+ *q++ = 0;
+#endif
+ }
+ else
+ {
+ if (i)
+ *q-- = *p++ | 0x8000;
+ else
+ *q-- = *p++;
+ }
+#endif
/* skip over guard word */
++p;
/* move the significand */
-#ifdef MIEEE
+#ifdef IBM
for (i = 0; i < 4; i++)
*q++ = *p++;
-#else
+#endif
+#ifdef DEC
for (i = 0; i < 4; i++)
*q-- = *p++;
#endif
+#ifdef IEEE
+ if (REAL_WORDS_BIG_ENDIAN)
+ {
+ for (i = 0; i < 4; i++)
+ *q++ = *p++;
+ }
+ else
+ {
+#ifdef INFINITY
+ if (eiisinf (a))
+ {
+ /* Intel long double infinity significand. */
+ *q-- = 0x8000;
+ *q-- = 0;
+ *q-- = 0;
+ *q = 0;
+ return;
+ }
+#endif
+ for (i = 0; i < 4; i++)
+ *q-- = *p++;
+ }
+#endif
}
-
-/*
-; e type to IEEE double precision
-; double d;
-; unsigned EMUSHORT x[NE];
-; etoe53 (x, &d);
-*/
+/* e type to double precision. */
#ifdef DEC
+/* Convert e-type X to DEC-format double E. */
-void
+static void
etoe53 (x, e)
unsigned EMUSHORT *x, *e;
{
etodec (x, e); /* see etodec.c */
}
+/* Convert exploded e-type X, that has already been rounded to
+ 56-bit double precision, to DEC double Y. */
+
static void
toe53 (x, y)
unsigned EMUSHORT *x, *y;
}
#else
+#ifdef IBM
+/* Convert e-type X to IBM 370-format double E. */
-void
+static void
+etoe53 (x, e)
+ unsigned EMUSHORT *x, *e;
+{
+ etoibm (x, e, DFmode);
+}
+
+/* Convert exploded e-type X, that has already been rounded to
+ 56-bit precision, to IBM 370 double Y. */
+
+static void
+toe53 (x, y)
+ unsigned EMUSHORT *x, *y;
+{
+ toibm (x, y, DFmode);
+}
+
+#else /* it's neither DEC nor IBM */
+
+/* Convert e-type X to IEEE double E. */
+
+static void
etoe53 (x, e)
unsigned EMUSHORT *x, *e;
{
EMULONG exp;
int rndsav;
+#ifdef NANS
+ if (eisnan (x))
+ {
+ make_nan (e, eisneg (x), DFmode);
+ return;
+ }
+#endif
emovi (x, xi);
/* adjust exponent for offsets */
exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
toe53 (xi, e);
}
+/* Convert exploded e-type X, that has already been rounded to
+ 53-bit precision, to IEEE double Y. */
static void
toe53 (x, y)
unsigned EMUSHORT i;
unsigned EMUSHORT *p;
-
+#ifdef NANS
+ if (eiisnan (x))
+ {
+ make_nan (y, eiisneg (x), DFmode);
+ return;
+ }
+#endif
p = &x[0];
-#ifdef IBMPC
- y += 3;
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
+ y += 3;
#endif
*y = 0; /* output high order */
if (*p++)
i = *p++;
if (i >= (unsigned int) 2047)
- { /* Saturate at largest number less than infinity. */
+ {
+ /* Saturate at largest number less than infinity. */
#ifdef INFINITY
*y |= 0x7ff0;
-#ifdef IBMPC
- *(--y) = 0;
- *(--y) = 0;
- *(--y) = 0;
-#endif
-#ifdef MIEEE
- ++y;
- *y++ = 0;
- *y++ = 0;
- *y++ = 0;
-#endif
+ if (! REAL_WORDS_BIG_ENDIAN)
+ {
+ *(--y) = 0;
+ *(--y) = 0;
+ *(--y) = 0;
+ }
+ else
+ {
+ ++y;
+ *y++ = 0;
+ *y++ = 0;
+ *y++ = 0;
+ }
#else
*y |= (unsigned EMUSHORT) 0x7fef;
-#ifdef IBMPC
- *(--y) = 0xffff;
- *(--y) = 0xffff;
- *(--y) = 0xffff;
-#endif
-#ifdef MIEEE
- ++y;
- *y++ = 0xffff;
- *y++ = 0xffff;
- *y++ = 0xffff;
-#endif
+ if (! REAL_WORDS_BIG_ENDIAN)
+ {
+ *(--y) = 0xffff;
+ *(--y) = 0xffff;
+ *(--y) = 0xffff;
+ }
+ else
+ {
+ ++y;
+ *y++ = 0xffff;
+ *y++ = 0xffff;
+ *y++ = 0xffff;
+ }
#endif
return;
}
}
i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
*y |= (unsigned EMUSHORT) i; /* high order output already has sign bit set */
-#ifdef IBMPC
- *(--y) = *p++;
- *(--y) = *p++;
- *(--y) = *p;
-#endif
-#ifdef MIEEE
- ++y;
- *y++ = *p++;
- *y++ = *p++;
- *y++ = *p++;
-#endif
+ if (! REAL_WORDS_BIG_ENDIAN)
+ {
+ *(--y) = *p++;
+ *(--y) = *p++;
+ *(--y) = *p;
+ }
+ else
+ {
+ ++y;
+ *y++ = *p++;
+ *y++ = *p++;
+ *y++ = *p++;
+ }
}
+#endif /* not IBM */
#endif /* not DEC */
-/*
-; e type to IEEE single precision
-; float d;
-; unsigned EMUSHORT x[N+2];
-; xtod (x, &d);
-*/
-void
+/* e type to single precision. */
+
+#ifdef IBM
+/* Convert e-type X to IBM 370 float E. */
+
+static void
+etoe24 (x, e)
+ unsigned EMUSHORT *x, *e;
+{
+ etoibm (x, e, SFmode);
+}
+
+/* Convert exploded e-type X, that has already been rounded to
+ float precision, to IBM 370 float Y. */
+
+static void
+toe24 (x, y)
+ unsigned EMUSHORT *x, *y;
+{
+ toibm (x, y, SFmode);
+}
+
+#else
+/* Convert e-type X to IEEE float E. DEC float is the same as IEEE float. */
+
+static void
etoe24 (x, e)
unsigned EMUSHORT *x, *e;
{
unsigned EMUSHORT xi[NI];
int rndsav;
+#ifdef NANS
+ if (eisnan (x))
+ {
+ make_nan (e, eisneg (x), SFmode);
+ return;
+ }
+#endif
emovi (x, xi);
/* adjust exponent for offsets */
exp = (EMULONG) xi[E] - (EXONE - 0177);
toe24 (xi, e);
}
+/* Convert exploded e-type X, that has already been rounded to
+ float precision, to IEEE float Y. */
+
static void
toe24 (x, y)
unsigned EMUSHORT *x, *y;
unsigned EMUSHORT i;
unsigned EMUSHORT *p;
+#ifdef NANS
+ if (eiisnan (x))
+ {
+ make_nan (y, eiisneg (x), SFmode);
+ return;
+ }
+#endif
p = &x[0];
-#ifdef IBMPC
- y += 1;
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
+ y += 1;
#endif
#ifdef DEC
y += 1;
*y = 0x8000; /* output sign bit */
i = *p++;
-/* Handle overflow cases. */
+/* Handle overflow cases. */
if (i >= 255)
{
#ifdef INFINITY
*y |= (unsigned EMUSHORT) 0x7f80;
-#ifdef IBMPC
- *(--y) = 0;
-#endif
#ifdef DEC
*(--y) = 0;
#endif
-#ifdef MIEEE
- ++y;
- *y = 0;
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
+ *(--y) = 0;
+ else
+ {
+ ++y;
+ *y = 0;
+ }
#endif
#else /* no INFINITY */
*y |= (unsigned EMUSHORT) 0x7f7f;
-#ifdef IBMPC
- *(--y) = 0xffff;
-#endif
#ifdef DEC
*(--y) = 0xffff;
#endif
-#ifdef MIEEE
- ++y;
- *y = 0xffff;
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
+ *(--y) = 0xffff;
+ else
+ {
+ ++y;
+ *y = 0xffff;
+ }
#endif
#ifdef ERANGE
errno = ERANGE;
eshift (x, 8);
}
i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
- *y |= i; /* high order output already has sign bit set */
-#ifdef IBMPC
- *(--y) = *p;
-#endif
+ /* High order output already has sign bit set. */
+ *y |= i;
#ifdef DEC
*(--y) = *p;
#endif
-#ifdef MIEEE
- ++y;
- *y = *p;
+#ifdef IEEE
+ if (! REAL_WORDS_BIG_ENDIAN)
+ *(--y) = *p;
+ else
+ {
+ ++y;
+ *y = *p;
+ }
#endif
}
+#endif /* not IBM */
+/* Compare two e type numbers.
+ Return +1 if a > b
+ 0 if a == b
+ -1 if a < b
+ -2 if either a or b is a NaN. */
-/* Compare two e type numbers.
- *
- * unsigned EMUSHORT a[NE], b[NE];
- * ecmp (a, b);
- *
- * returns +1 if a > b
- * 0 if a == b
- * -1 if a < b
- */
-int
+static int
ecmp (a, b)
unsigned EMUSHORT *a, *b;
{
register int i;
int msign;
+#ifdef NANS
+ if (eisnan (a) || eisnan (b))
+ return (-2);
+#endif
emovi (a, ai);
p = ai;
emovi (b, bi);
return (0); /* equality */
-
-
diff:
if (*(--p) > *(--q))
return (-msign); /* p is littler */
}
+/* Find e-type nearest integer to X, as floor (X + 0.5). */
-
-
-/* Find nearest integer to x = floor (x + 0.5)
- *
- * unsigned EMUSHORT x[NE], y[NE]
- * eround (x, y);
- */
-void
+static void
eround (x, y)
unsigned EMUSHORT *x, *y;
{
efloor (y, y);
}
+/* Convert HOST_WIDE_INT LP to e type Y. */
-
-
-/*
-; convert long integer to e type
-;
-; long l;
-; unsigned EMUSHORT x[NE];
-; ltoe (&l, x);
-; note &l is the memory address of l
-*/
-void
+static void
ltoe (lp, y)
- long *lp; /* lp is the memory address of a long integer */
- unsigned EMUSHORT *y; /* y is the address of a short */
+ HOST_WIDE_INT *lp;
+ unsigned EMUSHORT *y;
{
unsigned EMUSHORT yi[NI];
- unsigned long ll;
+ unsigned HOST_WIDE_INT ll;
int k;
ecleaz (yi);
if (*lp < 0)
{
/* make it positive */
- ll = (unsigned long) (-(*lp));
+ ll = (unsigned HOST_WIDE_INT) (-(*lp));
yi[0] = 0xffff; /* put correct sign in the e type number */
}
else
{
- ll = (unsigned long) (*lp);
+ ll = (unsigned HOST_WIDE_INT) (*lp);
}
/* move the long integer to yi significand area */
+#if HOST_BITS_PER_WIDE_INT == 64
+ yi[M] = (unsigned EMUSHORT) (ll >> 48);
+ yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
+ yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
+ yi[M + 3] = (unsigned EMUSHORT) ll;
+ yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
+#else
yi[M] = (unsigned EMUSHORT) (ll >> 16);
yi[M + 1] = (unsigned EMUSHORT) ll;
-
yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
+#endif
+
if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
ecleaz (yi); /* it was zero */
else
emovo (yi, y); /* output the answer */
}
-/*
-; convert unsigned long integer to e type
-;
-; unsigned long l;
-; unsigned EMUSHORT x[NE];
-; ltox (&l, x);
-; note &l is the memory address of l
-*/
-void
+/* Convert unsigned HOST_WIDE_INT LP to e type Y. */
+
+static void
ultoe (lp, y)
- unsigned long *lp; /* lp is the memory address of a long integer */
- unsigned EMUSHORT *y; /* y is the address of a short */
+ unsigned HOST_WIDE_INT *lp;
+ unsigned EMUSHORT *y;
{
unsigned EMUSHORT yi[NI];
- unsigned long ll;
+ unsigned HOST_WIDE_INT ll;
int k;
ecleaz (yi);
ll = *lp;
/* move the long integer to ayi significand area */
+#if HOST_BITS_PER_WIDE_INT == 64
+ yi[M] = (unsigned EMUSHORT) (ll >> 48);
+ yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
+ yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
+ yi[M + 3] = (unsigned EMUSHORT) ll;
+ yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
+#else
yi[M] = (unsigned EMUSHORT) (ll >> 16);
yi[M + 1] = (unsigned EMUSHORT) ll;
-
yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
+#endif
+
if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
ecleaz (yi); /* it was zero */
else
}
-/*
-; Find long integer and fractional parts
-
-; long i;
-; unsigned EMUSHORT x[NE], frac[NE];
-; xifrac (x, &i, frac);
+/* Find signed HOST_WIDE_INT integer I and floating point fractional
+ part FRAC of e-type (packed internal format) floating point input X.
+ The integer output I has the sign of the input, except that
+ positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
+ The output e-type fraction FRAC is the positive fractional
+ part of abs (X). */
- The integer output has the sign of the input. The fraction is
-the positive fractional part of abs (x).
-*/
-void
+static void
eifrac (x, i, frac)
unsigned EMUSHORT *x;
- long *i;
+ HOST_WIDE_INT *i;
unsigned EMUSHORT *frac;
{
unsigned EMUSHORT xi[NI];
- int k;
+ int j, k;
+ unsigned HOST_WIDE_INT ll;
emovi (x, xi);
k = (int) xi[E] - (EXONE - 1);
emovo (xi, frac);
return;
}
- if (k > (HOST_BITS_PER_LONG - 1))
+ if (k > (HOST_BITS_PER_WIDE_INT - 1))
{
- /*
- ; long integer overflow: output large integer
- ; and correct fraction
- */
+ /* long integer overflow: output large integer
+ and correct fraction */
if (xi[0])
- *i = ((unsigned long) 1) << (HOST_BITS_PER_LONG - 1);
+ *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
else
- *i = (((unsigned long) 1) << (HOST_BITS_PER_LONG - 1)) - 1;
+ {
+#ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
+ /* In this case, let it overflow and convert as if unsigned. */
+ euifrac (x, &ll, frac);
+ *i = (HOST_WIDE_INT) ll;
+ return;
+#else
+ /* In other cases, return the largest positive integer. */
+ *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
+#endif
+ }
eshift (xi, k);
if (extra_warnings)
warning ("overflow on truncation to integer");
- goto lab11;
}
-
- if (k > 16)
+ else if (k > 16)
{
- /*
- ; shift more than 16 bits: shift up k-16, output the integer,
- ; then complete the shift to get the fraction.
- */
- k -= 16;
- eshift (xi, k);
-
- *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
- eshup6 (xi);
- goto lab10;
+ /* Shift more than 16 bits: first shift up k-16 mod 16,
+ then shift up by 16's. */
+ j = k - ((k >> 4) << 4);
+ eshift (xi, j);
+ ll = xi[M];
+ k -= j;
+ do
+ {
+ eshup6 (xi);
+ ll = (ll << 16) | xi[M];
+ }
+ while ((k -= 16) > 0);
+ *i = ll;
+ if (xi[0])
+ *i = -(*i);
}
-
- /* shift not more than 16 bits */
- eshift (xi, k);
- *i = (long) xi[M] & 0xffff;
-
- lab10:
-
- if (xi[0])
- *i = -(*i);
- lab11:
-
+ else
+ {
+ /* shift not more than 16 bits */
+ eshift (xi, k);
+ *i = (HOST_WIDE_INT) xi[M] & 0xffff;
+ if (xi[0])
+ *i = -(*i);
+ }
xi[0] = 0;
xi[E] = EXONE - 1;
xi[M] = 0;
}
-/*
-; Find unsigned long integer and fractional parts
-
-; unsigned long i;
-; unsigned EMUSHORT x[NE], frac[NE];
-; xifrac (x, &i, frac);
+/* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
+ FRAC of e-type X. A negative input yields integer output = 0 but
+ correct fraction. */
- A negative e type input yields integer output = 0
- but correct fraction.
-*/
-void
+static void
euifrac (x, i, frac)
unsigned EMUSHORT *x;
- long *i;
+ unsigned HOST_WIDE_INT *i;
unsigned EMUSHORT *frac;
{
+ unsigned HOST_WIDE_INT ll;
unsigned EMUSHORT xi[NI];
- int k;
+ int j, k;
emovi (x, xi);
k = (int) xi[E] - (EXONE - 1);
emovo (xi, frac);
return;
}
- if (k > 32)
+ if (k > HOST_BITS_PER_WIDE_INT)
{
- /*
- ; long integer overflow: output large integer
- ; and correct fraction
- */
+ /* Long integer overflow: output large integer
+ and correct fraction.
+ Note, the BSD microvax compiler says that ~(0UL)
+ is a syntax error. */
*i = ~(0L);
eshift (xi, k);
if (extra_warnings)
warning ("overflow on truncation to unsigned integer");
- goto lab10;
}
-
- if (k > 16)
+ else if (k > 16)
{
- /*
- ; shift more than 16 bits: shift up k-16, output the integer,
- ; then complete the shift to get the fraction.
- */
- k -= 16;
+ /* Shift more than 16 bits: first shift up k-16 mod 16,
+ then shift up by 16's. */
+ j = k - ((k >> 4) << 4);
+ eshift (xi, j);
+ ll = xi[M];
+ k -= j;
+ do
+ {
+ eshup6 (xi);
+ ll = (ll << 16) | xi[M];
+ }
+ while ((k -= 16) > 0);
+ *i = ll;
+ }
+ else
+ {
+ /* shift not more than 16 bits */
eshift (xi, k);
-
- *i = (long) (((unsigned long) xi[M] << 16) | xi[M + 1]);
- eshup6 (xi);
- goto lab10;
+ *i = (HOST_WIDE_INT) xi[M] & 0xffff;
}
- /* shift not more than 16 bits */
- eshift (xi, k);
- *i = (long) xi[M] & 0xffff;
-
- lab10:
-
- if (xi[0])
+ if (xi[0]) /* A negative value yields unsigned integer 0. */
*i = 0L;
xi[0] = 0;
emovo (xi, frac);
}
+/* Shift the significand of exploded e-type X up or down by SC bits. */
-
-/*
-; Shift significand
-;
-; Shifts significand area up or down by the number of bits
-; given by the variable sc.
-*/
-int
+static int
eshift (x, sc)
unsigned EMUSHORT *x;
int sc;
return ((int) lost);
}
+/* Shift normalize the significand area of exploded e-type X.
+ Return the shift count (up = positive). */
-
-/*
-; normalize
-;
-; Shift normalizes the significand area pointed to by argument
-; shift count (up = positive) is returned.
-*/
-int
+static int
enormlz (x)
unsigned EMUSHORT x[];
{
{
eshup6 (x);
sc += 16;
+
/* With guard word, there are NBITS+16 bits available.
- * return true if all are zero.
- */
+ Return true if all are zero. */
if (sc > NBITS)
return (sc);
}
return (sc);
}
-
-
-
-/* Convert e type number to decimal format ASCII string.
- * The constants are for 64 bit precision.
- */
+/* Powers of ten used in decimal <-> binary conversions. */
#define NTEN 12
#define MAXP 4096
+#if LONG_DOUBLE_TYPE_SIZE == 128
+static unsigned EMUSHORT etens[NTEN + 1][NE] =
+{
+ {0x6576, 0x4a92, 0x804a, 0x153f,
+ 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
+ {0x6a32, 0xce52, 0x329a, 0x28ce,
+ 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
+ {0x526c, 0x50ce, 0xf18b, 0x3d28,
+ 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
+ {0x9c66, 0x58f8, 0xbc50, 0x5c54,
+ 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
+ {0x851e, 0xeab7, 0x98fe, 0x901b,
+ 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
+ {0x0235, 0x0137, 0x36b1, 0x336c,
+ 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
+ {0x50f8, 0x25fb, 0xc76b, 0x6b71,
+ 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
+ {0x0000, 0x0000, 0x0000, 0x0000,
+ 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
+};
+
+static unsigned EMUSHORT emtens[NTEN + 1][NE] =
+{
+ {0x2030, 0xcffc, 0xa1c3, 0x8123,
+ 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
+ {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
+ 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
+ {0xf53f, 0xf698, 0x6bd3, 0x0158,
+ 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
+ {0xe731, 0x04d4, 0xe3f2, 0xd332,
+ 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
+ {0xa23e, 0x5308, 0xfefb, 0x1155,
+ 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
+ {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
+ 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
+ {0x2a20, 0x6224, 0x47b3, 0x98d7,
+ 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
+ {0x0b5b, 0x4af2, 0xa581, 0x18ed,
+ 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
+ {0xbf71, 0xa9b3, 0x7989, 0xbe68,
+ 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
+ {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
+ 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
+ {0xc155, 0xa4a8, 0x404e, 0x6113,
+ 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
+ {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
+ 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
+ {0xcccd, 0xcccc, 0xcccc, 0xcccc,
+ 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
+};
+#else
+/* LONG_DOUBLE_TYPE_SIZE is other than 128 */
static unsigned EMUSHORT etens[NTEN + 1][NE] =
{
{0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
{0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
{0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
};
+#endif
-void
+/* Convert float value X to ASCII string STRING with NDIG digits after
+ the decimal point. */
+
+static void
e24toasc (x, string, ndigs)
unsigned EMUSHORT x[];
char *string;
{
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);
}
+/* Convert double value X to ASCII string STRING with NDIG digits after
+ the decimal point. */
-void
+static void
e53toasc (x, string, ndigs)
unsigned EMUSHORT x[];
char *string;
{
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);
}
+/* Convert double extended value X to ASCII string STRING with NDIG digits
+ after the decimal point. */
-void
+static void
e64toasc (x, string, ndigs)
unsigned EMUSHORT x[];
char *string;
{
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);
}
+/* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
+ after the decimal point. */
+
+static void
+e113toasc (x, string, ndigs)
+ unsigned EMUSHORT x[];
+ char *string;
+ int ndigs;
+{
+ unsigned EMUSHORT w[NI];
+
+ e113toe (x, w);
+ etoasc (w, string, ndigs);
+}
+
+/* Convert e-type X to ASCII string STRING with NDIGS digits after
+ the decimal point. */
static char wstring[80]; /* working storage for ASCII output */
-void
+static void
etoasc (x, string, ndigs)
unsigned EMUSHORT x[];
char *string;
char *s, *ss;
unsigned EMUSHORT m;
+
+ rndsav = rndprc;
ss = string;
s = wstring;
- while ((*s++ = *ss++) != '\0')
- ;
- rndsav = rndprc;
+ *ss = '\0';
+ *s = '\0';
+#ifdef NANS
+ if (eisnan (x))
+ {
+ sprintf (wstring, " NaN ");
+ goto bxit;
+ }
+#endif
rndprc = NBITS; /* set to full precision */
emov (x, y); /* retain external format */
if (y[NE - 1] & 0x8000)
if (y[k] != 0)
goto tnzro; /* denormalized number */
}
- goto isone; /* legal all zeros */
+ goto isone; /* valid all zeros */
}
tnzro:
- /* Test for infinity. Don't bother with illegal infinities.
- */
+ /* Test for infinity. */
if (y[NE - 1] == 0x7fff)
{
if (sign)
if (i == 0)
goto isone;
+ if (i == -2)
+ abort ();
+
if (i < 0)
{ /* Number is greater than 1 */
- /* Convert significand to an integer and strip trailing decimal zeros. */
+ /* Convert significand to an integer and strip trailing decimal zeros. */
emov (y, u);
u[NE - 1] = EXONE + NBITS - 1;
emov (eone, t);
m = MAXP;
p = &etens[0][0];
+ /* An unordered compare result shouldn't happen here. */
while (ecmp (ten, u) <= 0)
{
if (ecmp (p, u) <= 0)
}
else
{ /* Number is less than 1.0 */
- /* Pad significand with trailing decimal zeros. */
+ /* Pad significand with trailing decimal zeros. */
if (y[NE - 1] == 0)
{
while ((y[NE - 2] & 0x8000) == 0)
ediv (t, eone, t);
}
isone:
- /* Find the first (leading) digit. */
+ /* Find the first (leading) digit. */
emovi (t, w);
emovz (w, t);
emovi (y, w);
*s++ = '-';
else
*s++ = ' ';
- /* Examine number of digits requested by caller. */
+ /* Examine number of digits requested by caller. */
if (ndigs < 0)
ndigs = 0;
if (ndigs > NDEC)
}
else
{
- *s++ = (char )digit + '0';
+ *s++ = (char)digit + '0';
*s++ = '.';
}
- /* Generate digits after the decimal point. */
+ /* Generate digits after the decimal point. */
for (k = 0; k <= ndigs; k++)
{
/* multiply current number by 10, without normalizing */
/* round off the ASCII string */
if (digit > 4)
{
- /* Test for critical rounding case in ASCII output. */
+ /* Test for critical rounding case in ASCII output. */
if (digit == 5)
{
emovo (y, t);
}
+/* Convert ASCII string to floating point.
+ Numeric input is a free format decimal number of any length, with
+ or without decimal point. Entering E after the number followed by an
+ integer number causes the second number to be interpreted as a power of
+ 10 to be multiplied by the first number (i.e., "scientific" notation). */
-/*
-; ASCTOQ
-; ASCTOQ.MAC LATEST REV: 11 JAN 84
-; SLM, 3 JAN 78
-;
-; Convert ASCII string to quadruple precision floating point
-;
-; Numeric input is free field decimal number
-; with max of 15 digits with or without
-; decimal point entered as ASCII from teletype.
-; Entering E after the number followed by a second
-; number causes the second number to be interpreted
-; as a power of 10 to be multiplied by the first number
-; (i.e., "scientific" notation).
-;
-; Usage:
-; asctoq (string, q);
-*/
-
-/* ASCII to single */
-void
+/* Convert ASCII string S to single precision float value Y. */
+
+static void
asctoe24 (s, y)
char *s;
unsigned EMUSHORT *y;
}
-/* ASCII to double */
-void
+/* Convert ASCII string S to double precision value Y. */
+
+static void
asctoe53 (s, y)
char *s;
unsigned EMUSHORT *y;
{
-#ifdef DEC
+#if defined(DEC) || defined(IBM)
asctoeg (s, y, 56);
#else
asctoeg (s, y, 53);
}
-/* ASCII to long double */
-void
+/* Convert ASCII string S to double extended value Y. */
+
+static void
asctoe64 (s, y)
char *s;
unsigned EMUSHORT *y;
asctoeg (s, y, 64);
}
-/* ASCII to super double */
-void
+/* Convert ASCII string S to 128-bit long double Y. */
+
+static void
+asctoe113 (s, y)
+ char *s;
+ unsigned EMUSHORT *y;
+{
+ asctoeg (s, y, 113);
+}
+
+/* Convert ASCII string S to e type Y. */
+
+static void
asctoe (s, y)
char *s;
unsigned EMUSHORT *y;
asctoeg (s, y, NBITS);
}
-/* Space to make a copy of the input string: */
-static char lstr[82];
+/* Convert ASCII string SS to e type Y, with a specified rounding precision
+ of OPREC bits. */
-void
+static void
asctoeg (ss, y, oprec)
char *ss;
unsigned EMUSHORT *y;
int k, trail, c, rndsav;
EMULONG lexp;
unsigned EMUSHORT nsign, *p;
- char *sp, *s;
+ char *sp, *s, *lstr;
- /* Copy the input string. */
+ /* Copy the input string. */
+ lstr = (char *) alloca (strlen (ss) + 1);
s = ss;
while (*s == ' ') /* skip leading spaces */
++s;
sp = lstr;
- for (k = 0; k < 79; k++)
- {
- if ((*sp++ = *s++) == '\0')
- break;
- }
- *sp = '\0';
+ while ((*sp++ = *s++) != '\0')
+ ;
s = lstr;
rndsav = rndprc;
/* Ignore leading zeros */
if ((prec == 0) && (decflg == 0) && (k == 0))
goto donchr;
- /* Identify and strip trailing zeros after the decimal point. */
+ /* Identify and strip trailing zeros after the decimal point. */
if ((trail == 0) && (decflg != 0))
{
sp = s;
if (*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.
- */
+ 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)
}
else
{
+ /* Mark any lost non-zero digit. */
lost |= k;
+ /* Count lost digits before the decimal point. */
+ if (decflg == 0)
+ nexp -= 1;
}
prec += 1;
goto donchr;
goto infinite;
default:
error:
+#ifdef NANS
+ einan (yy);
+#else
mtherr ("asctoe", DOMAIN);
- eclear (y);
+ eclear (yy);
+#endif
goto aexit;
}
donchr:
/* Exponent interpretation */
expnt:
+ /* 0.0eXXX is zero, regardless of XXX. Check for the 0.0. */
+ for (k = 0; k < NI; k++)
+ {
+ if (yy[k] != 0)
+ goto read_expnt;
+ }
+ goto aexit;
+read_expnt:
esign = 1;
exp = 0;
++s;
{
exp *= 10;
exp += *s++ - '0';
- if (exp > 4956)
+ if (exp > -(MINDECEXP))
{
if (esign < 0)
goto zero;
}
if (esign < 0)
exp = -exp;
- if (exp > 4932)
+ if (exp > MAXDECEXP)
{
infinite:
ecleaz (yy);
yy[E] = 0x7fff; /* infinity */
goto aexit;
}
- if (exp < -4956)
+ if (exp < MINDECEXP)
{
zero:
ecleaz (yy);
daldone:
nexp = exp - nexp;
- /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
+ /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
while ((nexp > 0) && (yy[2] == 0))
{
emovz (yy, xt);
}
lexp = (EXONE - 1 + NBITS) - k;
emdnorm (yy, lost, 0, lexp, 64);
- /* convert to external format */
+ /* Convert to external format:
+
+ Multiply by 10**nexp. If precision is 64 bits,
+ the maximum relative error incurred in forming 10**n
+ for 0 <= n <= 324 is 8.2e-20, at 10**180.
+ For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
+ For 0 >= n >= -999, it is -1.55e-19 at 10**-435. */
- /* Multiply by 10**nexp. If precision is 64 bits,
- * the maximum relative error incurred in forming 10**n
- * for 0 <= n <= 324 is 8.2e-20, at 10**180.
- * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
- * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
- */
lexp = yy[E];
if (nexp == 0)
{
nexp = -nexp;
esign = -1;
if (nexp > 4096)
- { /* Punt. Can't handle this without 2 divides. */
+ {
+ /* Punt. Can't handle this without 2 divides. */
emovi (etens[0], tt);
lexp -= tt[E];
k = edivm (tt, yy);
/* Round and convert directly to the destination type */
if (oprec == 53)
lexp -= EXONE - 0x3ff;
+#ifdef IBM
+ else if (oprec == 24 || oprec == 56)
+ lexp -= EXONE - (0x41 << 2);
+#else
else if (oprec == 24)
lexp -= EXONE - 0177;
+#endif
#ifdef DEC
else if (oprec == 56)
lexp -= EXONE - 0201;
todec (yy, y); /* see etodec.c */
break;
#endif
+#ifdef IBM
+ case 56:
+ toibm (yy, y, DFmode);
+ break;
+#endif
case 53:
toe53 (yy, y);
break;
case 64:
toe64 (yy, y);
break;
+ case 113:
+ toe113 (yy, y);
+ break;
case NBITS:
emovo (yy, y);
break;
-/* y = largest integer not greater than x
- * (truncated toward minus infinity)
- *
- * unsigned EMUSHORT x[NE], y[NE]
- *
- * efloor (x, y);
- */
+/* Return Y = largest integer not greater than X (truncated toward minus
+ infinity). */
+
static unsigned EMUSHORT bmask[] =
{
0xffff,
0x0000,
};
-void
+static void
efloor (x, y)
unsigned EMUSHORT x[], y[];
{
}
-/* unsigned EMUSHORT x[], s[];
- * int *exp;
- *
- * efrexp (x, exp, s);
- *
- * Returns s and exp such that s * 2**exp = x and .5 <= s < 1.
- * For example, 1.1 = 0.55 * 2**1
- * Handles denormalized numbers properly using long integer exp.
- */
-void
+/* Return S and EXP such that S * 2^EXP = X and .5 <= S < 1.
+ For example, 1.1 = 0.55 * 2^1. */
+
+static void
efrexp (x, exp, s)
unsigned EMUSHORT x[];
int *exp;
EMULONG li;
emovi (x, xi);
+ /* Handle denormalized numbers properly using long integer exponent. */
li = (EMULONG) ((EMUSHORT) xi[1]);
if (li == 0)
*exp = (int) (li - 0x3ffe);
}
+/* Return e type Y = X * 2^PWR2. */
-
-/* unsigned EMUSHORT x[], y[];
- * long pwr2;
- *
- * eldexp (x, pwr2, y);
- *
- * Returns y = x * 2**pwr2.
- */
-void
+static void
eldexp (x, pwr2, y)
unsigned EMUSHORT x[];
int pwr2;
}
-/* c = remainder after dividing b by a
- * Least significant integer quotient bits left in equot[].
- */
-void
+/* C = remainder after dividing B by A, all e type values.
+ Least significant integer quotient bits left in EQUOT. */
+
+static void
eremain (a, b, c)
unsigned EMUSHORT a[], b[], c[];
{
unsigned EMUSHORT den[NI], num[NI];
+#ifdef NANS
+ if (eisinf (b)
+ || (ecmp (a, ezero) == 0)
+ || eisnan (a)
+ || eisnan (b))
+ {
+ enan (c, 0);
+ return;
+ }
+#endif
if (ecmp (a, ezero) == 0)
{
mtherr ("eremain", SING);
emovo (num, c);
}
-void
+/* Return quotient of exploded e-types NUM / DEN in EQUOT,
+ remainder in NUM. */
+
+static void
eiremain (den, num)
unsigned EMUSHORT den[], num[];
{
j = 1;
}
else
- {
j = 0;
- }
eshup1 (equot);
equot[NI - 1] |= j;
eshup1 (num);
emdnorm (num, 0, 0, ln, 0);
}
-/* mtherr.c
- *
- * Library common error handling routine
- *
- *
- *
- * SYNOPSIS:
- *
- * char *fctnam;
- * int code;
- * void mtherr ();
- *
- * mtherr (fctnam, code);
- *
- *
- *
- * DESCRIPTION:
- *
- * This routine may be called to report one of the following
- * error conditions (in the include file mconf.h).
- *
- * Mnemonic Value Significance
- *
- * DOMAIN 1 argument domain error
- * SING 2 function singularity
- * OVERFLOW 3 overflow range error
- * UNDERFLOW 4 underflow range error
- * TLOSS 5 total loss of precision
- * PLOSS 6 partial loss of precision
- * EDOM 33 Unix domain error code
- * ERANGE 34 Unix range error code
- *
- * The default version of the file prints the function name,
- * passed to it by the pointer fctnam, followed by the
- * error condition. The display is directed to the standard
- * output device. The routine then returns to the calling
- * program. Users may wish to modify the program to abort by
- * calling exit under severe error conditions such as domain
- * errors.
- *
- * Since all error conditions pass control to this function,
- * the display may be easily changed, eliminated, or directed
- * to an error logging device.
- *
- * SEE ALSO:
- *
- * mconf.h
- *
- */
-\f
-/*
-Cephes Math Library Release 2.0: April, 1987
-Copyright 1984, 1987 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-/* include "mconf.h" */
-
-/* Notice: the order of appearance of the following
- * messages is bound to the error codes defined
- * in mconf.h.
- */
-static char *ermsg[7] =
+/* Report an error condition CODE encountered in function NAME.
+ CODE is one of the following:
+
+ Mnemonic Value Significance
+
+ DOMAIN 1 argument domain error
+ SING 2 function singularity
+ OVERFLOW 3 overflow range error
+ UNDERFLOW 4 underflow range error
+ TLOSS 5 total loss of precision
+ PLOSS 6 partial loss of precision
+ INVALID 7 NaN - producing operation
+ EDOM 33 Unix domain error code
+ ERANGE 34 Unix range error code
+
+ The order of appearance of the following messages is bound to the
+ error codes defined above. */
+
+#define NMSGS 8
+static char *ermsg[NMSGS] =
{
"unknown", /* error code 0 */
"domain", /* error code 1 */
"overflow",
"underflow",
"total loss of precision",
- "partial loss of precision"
+ "partial loss of precision",
+ "invalid operation"
};
int merror = 0;
extern int merror;
-void
+static void
mtherr (name, code)
char *name;
int code;
{
char errstr[80];
- /* Display string passed by calling program,
- * which is supposed to be the name of the
- * function in which the error occurred.
- */
+ /* The string passed by the calling program is supposed to be the
+ name of the function in which the error occurred.
+ The code argument selects which error message string will be printed. */
- /* Display error message defined
- * by the code argument.
- */
- if ((code <= 0) || (code >= 6))
+ if ((code <= 0) || (code >= NMSGS))
code = 0;
- sprintf (errstr, "\n%s %s error\n", name, ermsg[code]);
+ sprintf (errstr, " %s %s error", name, ermsg[code]);
if (extra_warnings)
warning (errstr);
/* Set global error message word */
merror = code + 1;
-
- /* Return to calling
- * program
- */
}
-/* Here is etodec.c .
- *
- */
+#ifdef DEC
+/* Convert DEC double precision D to e type E. */
-/*
-; convert DEC double precision to e type
-; double d;
-; EMUSHORT e[NE];
-; dectoe (&d, e);
-*/
-void
+static void
dectoe (d, e)
unsigned EMUSHORT *d;
unsigned EMUSHORT *e;
emovo (y, e);
}
+/* Convert e type X to DEC double precision D. */
-
-/*
-; convert e type to DEC double precision
-; double d;
-; EMUSHORT e[NE];
-; etodec (e, &d);
-*/
-#if 0
-static unsigned EMUSHORT decbit[NI] = {0, 0, 0, 0, 0, 0, 0200, 0};
-
-void
-etodec (x, d)
- unsigned EMUSHORT *x, *d;
-{
- unsigned EMUSHORT xi[NI];
- register unsigned EMUSHORT r;
- int i, j;
-
- emovi (x, xi);
- *d = 0;
- if (xi[0] != 0)
- *d = 0100000;
- r = xi[E];
- if (r < (EXONE - 128))
- goto zout;
- i = xi[M + 4];
- if ((i & 0200) != 0)
- {
- if ((i & 0377) == 0200)
- {
- if ((i & 0400) != 0)
- {
- /* check all less significant bits */
- for (j = M + 5; j < NI; j++)
- {
- if (xi[j] != 0)
- goto yesrnd;
- }
- }
- goto nornd;
- }
- yesrnd:
- eaddm (decbit, xi);
- r -= enormlz (xi);
- }
-
- nornd:
-
- r -= EXONE;
- r += 0201;
- if (r < 0)
- {
- zout:
- *d++ = 0;
- *d++ = 0;
- *d++ = 0;
- *d++ = 0;
- return;
- }
- if (r >= 0377)
- {
- *d++ = 077777;
- *d++ = -1;
- *d++ = -1;
- *d++ = -1;
- return;
- }
- r &= 0377;
- r <<= 7;
- eshup8 (xi);
- xi[M] &= 0177;
- r |= xi[M];
- *d++ |= r;
- *d++ = xi[M + 1];
- *d++ = xi[M + 2];
- *d++ = xi[M + 3];
-}
-
-#else
-
-void
+static void
etodec (x, d)
unsigned EMUSHORT *x, *d;
{
int rndsav;
emovi (x, xi);
- exp = (EMULONG) xi[E] - (EXONE - 0201); /* adjust exponent for offsets */
-/* round off to nearest or even */
+ /* Adjust exponent for offsets. */
+ exp = (EMULONG) xi[E] - (EXONE - 0201);
+ /* Round off to nearest or even. */
rndsav = rndprc;
rndprc = 56;
emdnorm (xi, 0, 0, exp, 64);
todec (xi, d);
}
-void
+/* Convert exploded e-type X, that has already been rounded to
+ 56-bit precision, to DEC format double Y. */
+
+static void
todec (x, y)
unsigned EMUSHORT *x, *y;
{
*y++ = x[M + 2];
*y++ = x[M + 3];
}
+#endif /* DEC */
+
+#ifdef IBM
+/* Convert IBM single/double precision to e type. */
+
+static void
+ibmtoe (d, e, mode)
+ unsigned EMUSHORT *d;
+ unsigned EMUSHORT *e;
+ enum machine_mode mode;
+{
+ unsigned EMUSHORT y[NI];
+ register unsigned EMUSHORT r, *p;
+ int rndsav;
-#endif /* not 0 */
+ ecleaz (y); /* start with a zero */
+ p = y; /* point to our number */
+ r = *d; /* get IBM exponent word */
+ if (*d & (unsigned int) 0x8000)
+ *p = 0xffff; /* fill in our sign */
+ ++p; /* bump pointer to our exponent word */
+ r &= 0x7f00; /* strip the sign bit */
+ r >>= 6; /* shift exponent word down 6 bits */
+ /* in fact shift by 8 right and 2 left */
+ r += EXONE - (0x41 << 2); /* subtract IBM exponent offset */
+ /* add our e type exponent offset */
+ *p++ = r; /* to form our exponent */
+
+ *p++ = *d++ & 0xff; /* now do the high order mantissa */
+ /* strip off the IBM exponent and sign bits */
+ if (mode != SFmode) /* there are only 2 words in SFmode */
+ {
+ *p++ = *d++; /* fill in the rest of our mantissa */
+ *p++ = *d++;
+ }
+ *p = *d;
+
+ if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
+ y[0] = y[E] = 0;
+ else
+ y[E] -= 5 + enormlz (y); /* now normalise the mantissa */
+ /* handle change in RADIX */
+ emovo (y, e);
+}
+
+
+
+/* Convert e type to IBM single/double precision. */
+
+static void
+etoibm (x, d, mode)
+ unsigned EMUSHORT *x, *d;
+ enum machine_mode mode;
+{
+ unsigned EMUSHORT xi[NI];
+ EMULONG exp;
+ int rndsav;
+
+ emovi (x, xi);
+ exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2)); /* adjust exponent for offsets */
+ /* round off to nearest or even */
+ rndsav = rndprc;
+ rndprc = 56;
+ emdnorm (xi, 0, 0, exp, 64);
+ rndprc = rndsav;
+ toibm (xi, d, mode);
+}
+
+static void
+toibm (x, y, mode)
+ unsigned EMUSHORT *x, *y;
+ enum machine_mode mode;
+{
+ unsigned EMUSHORT i;
+ unsigned EMUSHORT *p;
+ int r;
+
+ p = x;
+ *y = 0;
+ if (*p++)
+ *y = 0x8000;
+ i = *p++;
+ if (i == 0)
+ {
+ *y++ = 0;
+ *y++ = 0;
+ if (mode != SFmode)
+ {
+ *y++ = 0;
+ *y++ = 0;
+ }
+ return;
+ }
+ r = i & 0x3;
+ i >>= 2;
+ if (i > 0x7f)
+ {
+ *y++ |= 0x7fff;
+ *y++ = 0xffff;
+ if (mode != SFmode)
+ {
+ *y++ = 0xffff;
+ *y++ = 0xffff;
+ }
+#ifdef ERANGE
+ errno = ERANGE;
+#endif
+ return;
+ }
+ i &= 0x7f;
+ *y |= (i << 8);
+ eshift (x, r + 5);
+ *y++ |= x[M];
+ *y++ = x[M + 1];
+ if (mode != SFmode)
+ {
+ *y++ = x[M + 2];
+ *y++ = x[M + 3];
+ }
+}
+#endif /* IBM */
+
+/* Output a binary NaN bit pattern in the target machine's format. */
+
+/* If special NaN bit patterns are required, define them in tm.h
+ as arrays of unsigned 16-bit shorts. Otherwise, use the default
+ patterns here. */
+#ifdef TFMODE_NAN
+TFMODE_NAN;
+#else
+#ifdef IEEE
+unsigned EMUSHORT TFbignan[8] =
+ {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
+unsigned EMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
+#endif
+#endif
+
+#ifdef XFMODE_NAN
+XFMODE_NAN;
+#else
+#ifdef IEEE
+unsigned EMUSHORT XFbignan[6] =
+ {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
+unsigned EMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
+#endif
+#endif
+
+#ifdef DFMODE_NAN
+DFMODE_NAN;
+#else
+#ifdef IEEE
+unsigned EMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
+unsigned EMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
+#endif
+#endif
+
+#ifdef SFMODE_NAN
+SFMODE_NAN;
+#else
+#ifdef IEEE
+unsigned EMUSHORT SFbignan[2] = {0x7fff, 0xffff};
+unsigned EMUSHORT SFlittlenan[2] = {0, 0xffc0};
+#endif
+#endif
+
+
+static void
+make_nan (nan, sign, mode)
+ unsigned EMUSHORT *nan;
+ int sign;
+ enum machine_mode mode;
+{
+ int n;
+ unsigned EMUSHORT *p;
+
+ switch (mode)
+ {
+/* Possibly the `reserved operand' patterns on a VAX can be
+ used like NaN's, but probably not in the same way as IEEE. */
+#if !defined(DEC) && !defined(IBM)
+ case TFmode:
+ n = 8;
+ if (REAL_WORDS_BIG_ENDIAN)
+ p = TFbignan;
+ else
+ p = TFlittlenan;
+ break;
+ case XFmode:
+ n = 6;
+ if (REAL_WORDS_BIG_ENDIAN)
+ p = XFbignan;
+ else
+ p = XFlittlenan;
+ break;
+ case DFmode:
+ n = 4;
+ if (REAL_WORDS_BIG_ENDIAN)
+ p = DFbignan;
+ else
+ p = DFlittlenan;
+ break;
+ case HFmode:
+ case SFmode:
+ n = 2;
+ if (REAL_WORDS_BIG_ENDIAN)
+ p = SFbignan;
+ else
+ p = SFlittlenan;
+ break;
+#endif
+ default:
+ abort ();
+ }
+ if (REAL_WORDS_BIG_ENDIAN)
+ *nan++ = (sign << 15) | *p++;
+ while (--n != 0)
+ *nan++ = *p++;
+ if (! REAL_WORDS_BIG_ENDIAN)
+ *nan = (sign << 15) | *p;
+}
+
+/* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
+ This is the inverse of the function `etarsingle' invoked by
+ REAL_VALUE_TO_TARGET_SINGLE. */
+
+REAL_VALUE_TYPE
+ereal_from_float (f)
+ HOST_WIDE_INT f;
+{
+ REAL_VALUE_TYPE r;
+ unsigned EMUSHORT s[2];
+ unsigned EMUSHORT e[NE];
+ /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
+ This is the inverse operation to what the function `endian' does. */
+ if (REAL_WORDS_BIG_ENDIAN)
+ {
+ s[0] = (unsigned EMUSHORT) (f >> 16);
+ s[1] = (unsigned EMUSHORT) f;
+ }
+ else
+ {
+ s[0] = (unsigned EMUSHORT) f;
+ s[1] = (unsigned EMUSHORT) (f >> 16);
+ }
+ /* Convert and promote the target float to E-type. */
+ e24toe (s, e);
+ /* Output E-type to REAL_VALUE_TYPE. */
+ PUT_REAL (e, &r);
+ return r;
+}
+
+
+/* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
+ This is the inverse of the function `etardouble' invoked by
+ REAL_VALUE_TO_TARGET_DOUBLE.
+
+ The DFmode is stored as an array of HOST_WIDE_INT in the target's
+ data format, with no holes in the bit packing. The first element
+ of the input array holds the bits that would come first in the
+ target computer's memory. */
+
+REAL_VALUE_TYPE
+ereal_from_double (d)
+ HOST_WIDE_INT d[];
+{
+ REAL_VALUE_TYPE r;
+ unsigned EMUSHORT s[4];
+ unsigned EMUSHORT e[NE];
+
+ /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
+ if (REAL_WORDS_BIG_ENDIAN)
+ {
+ s[0] = (unsigned EMUSHORT) (d[0] >> 16);
+ s[1] = (unsigned EMUSHORT) d[0];
+#if HOST_BITS_PER_WIDE_INT == 32
+ s[2] = (unsigned EMUSHORT) (d[1] >> 16);
+ s[3] = (unsigned EMUSHORT) d[1];
+#else
+ /* In this case the entire target double is contained in the
+ first array element. The second element of the input is
+ ignored. */
+ s[2] = (unsigned EMUSHORT) (d[0] >> 48);
+ s[3] = (unsigned EMUSHORT) (d[0] >> 32);
+#endif
+ }
+ else
+ {
+ /* Target float words are little-endian. */
+ s[0] = (unsigned EMUSHORT) d[0];
+ s[1] = (unsigned EMUSHORT) (d[0] >> 16);
+#if HOST_BITS_PER_WIDE_INT == 32
+ s[2] = (unsigned EMUSHORT) d[1];
+ s[3] = (unsigned EMUSHORT) (d[1] >> 16);
+#else
+ s[2] = (unsigned EMUSHORT) (d[0] >> 32);
+ s[3] = (unsigned EMUSHORT) (d[0] >> 48);
+#endif
+ }
+ /* Convert target double to E-type. */
+ e53toe (s, e);
+ /* Output E-type to REAL_VALUE_TYPE. */
+ PUT_REAL (e, &r);
+ return r;
+}
+
+
+/* Convert target computer unsigned 64-bit integer to e-type.
+ The endian-ness of DImode follows the convention for integers,
+ so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN. */
+
+static void
+uditoe (di, e)
+ unsigned EMUSHORT *di; /* Address of the 64-bit int. */
+ unsigned EMUSHORT *e;
+{
+ unsigned EMUSHORT yi[NI];
+ int k;
+
+ ecleaz (yi);
+ if (WORDS_BIG_ENDIAN)
+ {
+ for (k = M; k < M + 4; k++)
+ yi[k] = *di++;
+ }
+ else
+ {
+ for (k = M + 3; k >= M; k--)
+ yi[k] = *di++;
+ }
+ yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
+ if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
+ ecleaz (yi); /* it was zero */
+ else
+ yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
+ emovo (yi, e);
+}
+
+/* Convert target computer signed 64-bit integer to e-type. */
+
+static void
+ditoe (di, e)
+ unsigned EMUSHORT *di; /* Address of the 64-bit int. */
+ unsigned EMUSHORT *e;
+{
+ unsigned EMULONG acc;
+ unsigned EMUSHORT yi[NI];
+ unsigned EMUSHORT carry;
+ int k, sign;
+
+ ecleaz (yi);
+ if (WORDS_BIG_ENDIAN)
+ {
+ for (k = M; k < M + 4; k++)
+ yi[k] = *di++;
+ }
+ else
+ {
+ for (k = M + 3; k >= M; k--)
+ yi[k] = *di++;
+ }
+ /* Take absolute value */
+ sign = 0;
+ if (yi[M] & 0x8000)
+ {
+ sign = 1;
+ carry = 0;
+ for (k = M + 3; k >= M; k--)
+ {
+ acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
+ yi[k] = acc;
+ carry = 0;
+ if (acc & 0x10000)
+ carry = 1;
+ }
+ }
+ yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
+ if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
+ ecleaz (yi); /* it was zero */
+ else
+ yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
+ emovo (yi, e);
+ if (sign)
+ eneg (e);
+}
+
+
+/* Convert e-type to unsigned 64-bit int. */
+
+static void
+etoudi (x, i)
+ unsigned EMUSHORT *x;
+ unsigned EMUSHORT *i;
+{
+ unsigned EMUSHORT xi[NI];
+ int j, k;
+
+ emovi (x, xi);
+ if (xi[0])
+ {
+ xi[M] = 0;
+ goto noshift;
+ }
+ k = (int) xi[E] - (EXONE - 1);
+ if (k <= 0)
+ {
+ for (j = 0; j < 4; j++)
+ *i++ = 0;
+ return;
+ }
+ if (k > 64)
+ {
+ for (j = 0; j < 4; j++)
+ *i++ = 0xffff;
+ if (extra_warnings)
+ warning ("overflow on truncation to integer");
+ return;
+ }
+ if (k > 16)
+ {
+ /* Shift more than 16 bits: first shift up k-16 mod 16,
+ then shift up by 16's. */
+ j = k - ((k >> 4) << 4);
+ if (j == 0)
+ j = 16;
+ eshift (xi, j);
+ if (WORDS_BIG_ENDIAN)
+ *i++ = xi[M];
+ else
+ {
+ i += 3;
+ *i-- = xi[M];
+ }
+ k -= j;
+ do
+ {
+ eshup6 (xi);
+ if (WORDS_BIG_ENDIAN)
+ *i++ = xi[M];
+ else
+ *i-- = xi[M];
+ }
+ while ((k -= 16) > 0);
+ }
+ else
+ {
+ /* shift not more than 16 bits */
+ eshift (xi, k);
+
+noshift:
+
+ if (WORDS_BIG_ENDIAN)
+ {
+ i += 3;
+ *i-- = xi[M];
+ *i-- = 0;
+ *i-- = 0;
+ *i = 0;
+ }
+ else
+ {
+ *i++ = xi[M];
+ *i++ = 0;
+ *i++ = 0;
+ *i = 0;
+ }
+ }
+}
+
+
+/* Convert e-type to signed 64-bit int. */
+
+static void
+etodi (x, i)
+ unsigned EMUSHORT *x;
+ unsigned EMUSHORT *i;
+{
+ unsigned EMULONG acc;
+ unsigned EMUSHORT xi[NI];
+ unsigned EMUSHORT carry;
+ unsigned EMUSHORT *isave;
+ int j, k;
+
+ emovi (x, xi);
+ k = (int) xi[E] - (EXONE - 1);
+ if (k <= 0)
+ {
+ for (j = 0; j < 4; j++)
+ *i++ = 0;
+ return;
+ }
+ if (k > 64)
+ {
+ for (j = 0; j < 4; j++)
+ *i++ = 0xffff;
+ if (extra_warnings)
+ warning ("overflow on truncation to integer");
+ return;
+ }
+ isave = i;
+ if (k > 16)
+ {
+ /* Shift more than 16 bits: first shift up k-16 mod 16,
+ then shift up by 16's. */
+ j = k - ((k >> 4) << 4);
+ if (j == 0)
+ j = 16;
+ eshift (xi, j);
+ if (WORDS_BIG_ENDIAN)
+ *i++ = xi[M];
+ else
+ {
+ i += 3;
+ *i-- = xi[M];
+ }
+ k -= j;
+ do
+ {
+ eshup6 (xi);
+ if (WORDS_BIG_ENDIAN)
+ *i++ = xi[M];
+ else
+ *i-- = xi[M];
+ }
+ while ((k -= 16) > 0);
+ }
+ else
+ {
+ /* shift not more than 16 bits */
+ eshift (xi, k);
+
+ if (WORDS_BIG_ENDIAN)
+ {
+ i += 3;
+ *i = xi[M];
+ *i-- = 0;
+ *i-- = 0;
+ *i = 0;
+ }
+ else
+ {
+ *i++ = xi[M];
+ *i++ = 0;
+ *i++ = 0;
+ *i = 0;
+ }
+ }
+ /* Negate if negative */
+ if (xi[0])
+ {
+ carry = 0;
+ if (WORDS_BIG_ENDIAN)
+ isave += 3;
+ for (k = 0; k < 4; k++)
+ {
+ acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
+ if (WORDS_BIG_ENDIAN)
+ *isave-- = acc;
+ else
+ *isave++ = acc;
+ carry = 0;
+ if (acc & 0x10000)
+ carry = 1;
+ }
+ }
+}
+
+
+/* Longhand square root routine. */
+
+
+static int esqinited = 0;
+static unsigned short sqrndbit[NI];
+
+static void
+esqrt (x, y)
+ unsigned EMUSHORT *x, *y;
+{
+ unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
+ EMULONG m, exp;
+ int i, j, k, n, nlups;
+
+ if (esqinited == 0)
+ {
+ ecleaz (sqrndbit);
+ sqrndbit[NI - 2] = 1;
+ esqinited = 1;
+ }
+ /* Check for arg <= 0 */
+ i = ecmp (x, ezero);
+ if (i <= 0)
+ {
+ if (i == -1)
+ {
+ mtherr ("esqrt", DOMAIN);
+ eclear (y);
+ }
+ else
+ emov (x, y);
+ return;
+ }
+
+#ifdef INFINITY
+ if (eisinf (x))
+ {
+ eclear (y);
+ einfin (y);
+ return;
+ }
+#endif
+ /* Bring in the arg and renormalize if it is denormal. */
+ emovi (x, xx);
+ m = (EMULONG) xx[1]; /* local long word exponent */
+ if (m == 0)
+ m -= enormlz (xx);
+
+ /* Divide exponent by 2 */
+ m -= 0x3ffe;
+ exp = (unsigned short) ((m / 2) + 0x3ffe);
+
+ /* Adjust if exponent odd */
+ if ((m & 1) != 0)
+ {
+ if (m > 0)
+ exp += 1;
+ eshdn1 (xx);
+ }
+
+ ecleaz (sq);
+ ecleaz (num);
+ n = 8; /* get 8 bits of result per inner loop */
+ nlups = rndprc;
+ j = 0;
+
+ while (nlups > 0)
+ {
+ /* bring in next word of arg */
+ if (j < NE)
+ num[NI - 1] = xx[j + 3];
+ /* Do additional bit on last outer loop, for roundoff. */
+ if (nlups <= 8)
+ n = nlups + 1;
+ for (i = 0; i < n; i++)
+ {
+ /* Next 2 bits of arg */
+ eshup1 (num);
+ eshup1 (num);
+ /* Shift up answer */
+ eshup1 (sq);
+ /* Make trial divisor */
+ for (k = 0; k < NI; k++)
+ temp[k] = sq[k];
+ eshup1 (temp);
+ eaddm (sqrndbit, temp);
+ /* Subtract and insert answer bit if it goes in */
+ if (ecmpm (temp, num) <= 0)
+ {
+ esubm (temp, num);
+ sq[NI - 2] |= 1;
+ }
+ }
+ nlups -= n;
+ j += 1;
+ }
+
+ /* Adjust for extra, roundoff loop done. */
+ exp += (NBITS - 1) - rndprc;
+
+ /* Sticky bit = 1 if the remainder is nonzero. */
+ k = 0;
+ for (i = 3; i < NI; i++)
+ k |= (int) num[i];
+
+ /* Renormalize and round off. */
+ emdnorm (sq, k, 0, exp, 64);
+ emovo (sq, y);
+}
#endif /* EMU_NON_COMPILE not defined */
+\f
+/* Return the binary precision of the significand for a given
+ floating point mode. The mode can hold an integer value
+ that many bits wide, without losing any bits. */
+
+int
+significand_size (mode)
+ enum machine_mode mode;
+{
+
+/* Don't test the modes, but their sizes, lest this
+ code won't work for BITS_PER_UNIT != 8 . */
+
+switch (GET_MODE_BITSIZE (mode))
+ {
+ case 32:
+ return 24;
+
+ case 64:
+#if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
+ return 53;
+#else
+#if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
+ return 56;
+#else
+#if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
+ return 56;
+#else
+ abort ();
+#endif
+#endif
+#endif
+
+ case 96:
+ return 64;
+ case 128:
+ return 113;
+
+ default:
+ abort ();
+ }
+}
OSDN Git Service
