X-Git-Url: http://git.sourceforge.jp/view?a=blobdiff_plain;f=gcc%2Freal.c;h=f7e22eae3bbee36c817fc9da86ffed201c76c1fc;hb=be44de9a82ac199b984024851a57a9dab31fb743;hp=518aaa8e6b000ce92ce914e916009b9639660c27;hpb=b6941bbc817c0d1e09c99a17b0d0a03c65969be6;p=pf3gnuchains%2Fgcc-fork.git diff --git a/gcc/real.c b/gcc/real.c index 518aaa8e6b0..f7e22eae3bb 100644 --- a/gcc/real.c +++ b/gcc/real.c @@ -1,8 +1,7 @@ /* 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. @@ -18,7 +17,8 @@ GNU General Public License for more details. You should have received a copy of the GNU General Public License along with GNU CC; see the file COPYING. If not, write to -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 #include @@ -32,7 +32,7 @@ extern int errno; /* To enable support of XFmode extended real floating point, define LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h). -To 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 @@ -45,27 +45,65 @@ The emulator defaults to the host's floating point format so that its decimal conversion functions can be used if desired (see real.h). -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. */ + /* 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 , + 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. */ @@ -75,81 +113,65 @@ research.att.com: netlib/cephes/ldouble.shar.Z */ /* 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 + /* 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 @@ -171,7 +193,7 @@ unknown arithmetic type #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 @@ -190,7 +212,7 @@ unknown arithmetic type #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 @@ -198,12 +220,12 @@ unknown arithmetic type #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. @@ -212,66 +234,199 @@ unknown arithmetic 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 *)); +static void esubm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void m16m PROTO((unsigned int, unsigned short *, + unsigned short *)); +static int edivm PROTO((unsigned short *, unsigned short *)); +static int emulm PROTO((unsigned short *, unsigned short *)); +static void emdnorm PROTO((unsigned EMUSHORT *, int, int, EMULONG, int)); +static void esub PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *, + unsigned EMUSHORT *)); +static void eadd PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *, + unsigned EMUSHORT *)); +static void eadd1 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *, + unsigned EMUSHORT *)); +static void ediv PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *, + unsigned EMUSHORT *)); +static void emul PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *, + unsigned EMUSHORT *)); +static void e53toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void e64toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void e113toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void e24toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void etoe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void toe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void etoe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void toe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void etoe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void toe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void etoe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void toe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static int ecmp PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void eround PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void ltoe PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *)); +static void ultoe PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *)); +static void eifrac PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *, + unsigned EMUSHORT *)); +static void euifrac PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *, + unsigned EMUSHORT *)); +static int eshift PROTO((unsigned EMUSHORT *, int)); +static int enormlz PROTO((unsigned EMUSHORT *)); +static void e24toasc PROTO((unsigned EMUSHORT *, char *, int)); +static void e53toasc PROTO((unsigned EMUSHORT *, char *, int)); +static void e64toasc PROTO((unsigned EMUSHORT *, char *, int)); +static void e113toasc PROTO((unsigned EMUSHORT *, char *, int)); +static void etoasc PROTO((unsigned EMUSHORT *, char *, int)); +static void asctoe24 PROTO((char *, unsigned EMUSHORT *)); +static void asctoe53 PROTO((char *, unsigned EMUSHORT *)); +static void asctoe64 PROTO((char *, unsigned EMUSHORT *)); +static void asctoe113 PROTO((char *, unsigned EMUSHORT *)); +static void asctoe PROTO((char *, unsigned EMUSHORT *)); +static void asctoeg PROTO((char *, unsigned EMUSHORT *, int)); +static void efloor PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void efrexp PROTO((unsigned EMUSHORT *, int *, + unsigned EMUSHORT *)); +static void eldexp PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *)); +static void eremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *, + unsigned EMUSHORT *)); +static void eiremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void mtherr PROTO((char *, int)); +static void dectoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void etodec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void todec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void ibmtoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *, + enum machine_mode)); +static void etoibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *, + enum machine_mode)); +static void toibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *, + enum machine_mode)); +static void make_nan PROTO((unsigned EMUSHORT *, int, enum machine_mode)); +static void uditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void ditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void etoudi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void etodi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); +static void esqrt PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *)); + +/* Copy 32-bit numbers obtained from array containing 16-bit numbers, + swapping ends if required, into output array of longs. The + result is normally passed to fprintf by the ASM_OUTPUT_ macros. */ + +static void endian (e, x, mode) unsigned EMUSHORT e[]; long x[]; @@ -279,87 +434,104 @@ endian (e, x, mode) { 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; @@ -372,6 +544,19 @@ earith (value, icode, r1, r2) 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) { @@ -390,8 +575,15 @@ earith (value, icode, r1, r2) 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; @@ -416,18 +608,22 @@ PUT_REAL (v, value); } -/* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT - * implements REAL_VALUE_FIX_TRUNCATE (x) (etrunci (x)) - */ +/* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT. + implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */ + REAL_VALUE_TYPE etrunci (x) REAL_VALUE_TYPE x; { unsigned EMUSHORT f[NE], g[NE]; REAL_VALUE_TYPE r; - 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); @@ -435,18 +631,22 @@ etrunci (x) } -/* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT - * implements REAL_VALUE_UNSIGNED_FIX_TRUNCATE (x) (etruncui (x)) - */ +/* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT; + implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */ + REAL_VALUE_TYPE etruncui (x) REAL_VALUE_TYPE x; { unsigned EMUSHORT f[NE], g[NE]; REAL_VALUE_TYPE r; - 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); @@ -454,11 +654,10 @@ etruncui (x) } -/* This is the REAL_VALUE_ATOF function. - * It converts a decimal string to binary, rounding off - * as indicated by the machine_mode argument. Then it - * promotes the rounded value to REAL_VALUE_TYPE. - */ +/* This is the REAL_VALUE_ATOF function. It converts a decimal string to + binary, rounding off as indicated by the machine_mode argument. Then it + promotes the rounded value to REAL_VALUE_TYPE. */ + REAL_VALUE_TYPE ereal_atof (s, t) char *s; @@ -469,6 +668,7 @@ ereal_atof (s, t) switch (t) { + case HFmode: case SFmode: asctoe24 (s, tem); e24toe (tem, e); @@ -481,6 +681,10 @@ ereal_atof (s, t) asctoe64 (s, tem); e64toe (tem, e); break; + case TFmode: + asctoe113 (s, tem); + e113toe (tem, e); + break; default: asctoe (s, e); } @@ -489,8 +693,8 @@ ereal_atof (s, t) } -/* Expansion of REAL_NEGATE. - */ +/* Expansion of REAL_NEGATE. */ + REAL_VALUE_TYPE ereal_negate (x) REAL_VALUE_TYPE x; @@ -505,55 +709,66 @@ ereal_negate (x) } -/* Round real to int - * implements REAL_VALUE_FIX (x) (eroundi (x)) - * The type of rounding is left unspecified by real.h. - * It is implemented here as round to nearest (add .5 and chop). - */ -int -eroundi (x) +/* Round real toward zero to HOST_WIDE_INT; + implements REAL_VALUE_FIX (x). */ + +HOST_WIDE_INT +efixi (x) REAL_VALUE_TYPE x; { unsigned EMUSHORT f[NE], g[NE]; - 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) @@ -566,49 +781,121 @@ ereal_from_int (d, i, j) 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)) @@ -616,11 +903,11 @@ ereal_to_int (low, high, rr) 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 */ @@ -633,8 +920,8 @@ ereal_to_int (low, high, rr) } -/* REAL_VALUE_LDEXP macro. - */ +/* REAL_VALUE_LDEXP macro. */ + REAL_VALUE_TYPE ereal_ldexp (x, n) REAL_VALUE_TYPE x; @@ -644,16 +931,22 @@ ereal_ldexp (x, n) 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; @@ -668,36 +961,36 @@ target_isinf (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; @@ -707,9 +1000,18 @@ real_value_truncate (mode, arg) 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); @@ -720,26 +1022,135 @@ real_value_truncate (mode, arg) 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); +} + + +/* 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; @@ -752,6 +1163,9 @@ etarldouble (r, l) 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; @@ -764,12 +1178,15 @@ etardouble (r, l) 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); @@ -777,6 +1194,11 @@ etarsingle (r) 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; @@ -788,6 +1210,9 @@ ereal_to_decimal (x, s) 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; @@ -799,6 +1224,8 @@ ereal_cmp (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; @@ -810,178 +1237,163 @@ ereal_isneg (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). - * - */ - + /* - * 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. - */ - /* Constant definitions for math error conditions. */ #define DOMAIN 1 /* argument domain error */ @@ -990,107 +1402,88 @@ ereal_isneg (x) #define UNDERFLOW 4 /* underflow range error */ #define TLOSS 5 /* total loss of precision */ #define PLOSS 6 /* partial loss of precision */ +#define 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; { @@ -1100,14 +1493,9 @@ eclear (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; { @@ -1118,32 +1506,19 @@ emov (a, b) } -/* -; Absolute value of external format number -; -; EMUSHORT x[NE]; -; eabs (x); -*/ +/* Absolute value of e-type X. */ -void +static void eabs (x) - 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[]; { @@ -1151,12 +1526,9 @@ eneg (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[]; { @@ -1167,32 +1539,51 @@ eisneg (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; { @@ -1208,6 +1599,11 @@ einfin (x) *x |= 32766; if (rndprc < NBITS) { + if (rndprc == 113) + { + *(x - 9) = 0; + *(x - 8) = 0; + } if (rndprc == 64) { *(x - 5) = 0; @@ -1226,12 +1622,26 @@ einfin (x) #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; { @@ -1251,11 +1661,22 @@ emovi (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 */ @@ -1265,16 +1686,15 @@ emovi (a, b) *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 */ @@ -1287,24 +1707,27 @@ emovo (a, b) #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; { @@ -1314,10 +1737,9 @@ ecleaz (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; { @@ -1328,11 +1750,9 @@ ecleazs (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; { @@ -1344,20 +1764,86 @@ emovz (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; { @@ -1379,12 +1865,9 @@ ecmpm (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; { @@ -1406,13 +1889,9 @@ eshdn1 (x) } } +/* Shift significand of exploded e-type X up by 1 bit. */ - -/* -; Shift significand up by 1 bit -*/ - -void +static void eshup1 (x) register unsigned EMUSHORT *x; { @@ -1435,12 +1914,9 @@ eshup1 (x) } +/* Shift significand of exploded e-type X down by 8 bits. */ -/* -; Shift significand down by 8 bits -*/ - -void +static void eshdn8 (x) register unsigned EMUSHORT *x; { @@ -1459,11 +1935,9 @@ eshdn8 (x) } } -/* -; Shift significand up by 8 bits -*/ +/* Shift significand of exploded e-type X up by 8 bits. */ -void +static void eshup8 (x) register unsigned EMUSHORT *x; { @@ -1483,11 +1957,9 @@ eshup8 (x) } } -/* -; Shift significand up by 16 bits -*/ +/* Shift significand of exploded e-type X up by 16 bits. */ -void +static void eshup6 (x) register unsigned EMUSHORT *x; { @@ -1503,11 +1975,9 @@ eshup6 (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; { @@ -1522,13 +1992,10 @@ eshdn6 (x) *(--p) = 0; } - -/* -; 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; { @@ -1552,12 +2019,9 @@ eaddm (x, y) } } -/* -; Subtract significands -; y - x replaces y -*/ +/* Subtract significands of exploded e-type X and Y. Y - X replaces Y. */ -void +static void esubm (x, y) unsigned EMUSHORT *x, *y; { @@ -1582,10 +2046,15 @@ esubm (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[]; @@ -1603,9 +2072,9 @@ edivm (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) { @@ -1640,9 +2109,9 @@ edivm (den, num) 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]; @@ -1683,6 +2152,7 @@ edivm (den, num) /* Multiply significands */ + int emulm (a, b) unsigned EMUSHORT a[], b[]; @@ -1697,7 +2167,7 @@ emulm (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; @@ -1728,61 +2198,225 @@ emulm (a, b) 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> 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 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 NBITS) + { + ecleazs (s); + return; + } +#endif + exp -= j; +#ifndef INFINITY + if (exp >= 32767L) + goto overf; +#else if ((j > NBITS) && (exp < 32767)) { ecleazs (s); @@ -1804,10 +2438,10 @@ emdnorm (s, lost, subflg, exp, rcntrl) 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); @@ -1818,68 +2452,64 @@ emdnorm (s, lost, subflg, exp, rcntrl) 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) { @@ -1888,17 +2518,8 @@ emdnorm (s, lost, subflg, exp, rcntrl) 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) @@ -1917,7 +2538,9 @@ emdnorm (s, lost, subflg, exp, rcntrl) 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); } @@ -1945,7 +2568,7 @@ emdnorm (s, lost, subflg, exp, rcntrl) for (i = M + 1; i < NI - 1; i++) s[i] = 0xffff; s[NI - 1] = 0; - if (rndprc < 64) + if ((rndprc < 64) || (rndprc == 113)) { s[rw] &= ~rmsk; if (rndprc == 24) @@ -1963,43 +2586,76 @@ emdnorm (s, lost, subflg, exp, rcntrl) 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; { @@ -2059,7 +2715,7 @@ eadd1 (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); @@ -2070,8 +2726,15 @@ eadd1 (a, b, c) { 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; } } @@ -2101,36 +2764,53 @@ eadd1 (a, b, c) emovo (bi, c); } +/* Divide: C = B/A, all e type. */ - -/* -; Divide. -; -; unsigned EMUSHORT a[NE], b[NE], c[NE]; -; ediv (a, b, c); c = b / a -*/ -void +static void ediv (a, b, c) unsigned EMUSHORT *a, *b, *c; { unsigned EMUSHORT ai[NI], bi[NI]; - 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); @@ -2138,7 +2818,7 @@ ediv (a, b, c) 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) @@ -2148,7 +2828,7 @@ ediv (a, b, c) } } eclear (c); - return; + goto divsign; } dnzro1: @@ -2162,13 +2842,11 @@ ediv (a, b, c) 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: @@ -2176,39 +2854,61 @@ ediv (a, b, c) /* 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); @@ -2226,7 +2926,7 @@ emul (a, b, c) } } eclear (c); - return; + goto mulsign; } mnzer1: @@ -2241,7 +2941,7 @@ emul (a, b, c) } } eclear (c); - return; + goto mulsign; } mnzer2: @@ -2250,43 +2950,46 @@ emul (a, b, c) /* 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) @@ -2296,15 +2999,37 @@ e53toe (e, y) #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; @@ -2313,16 +3038,20 @@ e53toe (e, y) 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) @@ -2333,70 +3062,239 @@ e53toe (e, y) 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; @@ -2410,15 +3308,34 @@ e24toe (e, y) #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; @@ -2427,15 +3344,17 @@ e24toe (e, y) 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) @@ -2446,10 +3365,103 @@ e24toe (e, y) yy[E] -= (unsigned EMUSHORT) (k - 1); } emovo (yy, y); +#endif /* not IBM */ } +/* Convert e-type X to IEEE 128-bit long double format E. */ -void +static void +etoe113 (x, e) + unsigned EMUSHORT *x, *e; +{ + unsigned EMUSHORT xi[NI]; + EMULONG exp; + int rndsav; + +#ifdef NANS + if (eisnan (x)) + { + make_nan (e, eisneg (x), TFmode); + return; + } +#endif + emovi (x, xi); + exp = (EMULONG) xi[E]; +#ifdef INFINITY + if (eisinf (x)) + goto nonorm; +#endif + /* round off to nearest or even */ + rndsav = rndprc; + rndprc = 113; + emdnorm (xi, 0, 0, exp, 64); + rndprc = rndsav; + nonorm: + toe113 (xi, e); +} + +/* Convert exploded e-type X, that has already been rounded to + 113-bit precision, to IEEE 128-bit long double format Y. */ + +static void +toe113 (a, b) + unsigned EMUSHORT *a, *b; +{ + register unsigned EMUSHORT *p, *q; + unsigned EMUSHORT i; + +#ifdef NANS + if (eiisnan (a)) + { + make_nan (b, eiisneg (a), TFmode); + return; + } +#endif + p = a; + if (REAL_WORDS_BIG_ENDIAN) + q = b; + else + q = b + 7; /* point to output exponent */ + + /* If not denormal, delete the implied bit. */ + if (a[E] != 0) + { + eshup1 (a); + } + /* combine sign and exponent */ + i = *p++; + if (REAL_WORDS_BIG_ENDIAN) + { + if (i) + *q++ = *p++ | 0x8000; + else + *q++ = *p++; + } + else + { + if (i) + *q-- = *p++ | 0x8000; + else + *q-- = *p++; + } + /* skip over guard word */ + ++p; + /* move the significand */ + if (REAL_WORDS_BIG_ENDIAN) + { + for (i = 0; i < 7; i++) + *q++ = *p++; + } + else + { + for (i = 0; i < 7; i++) + *q-- = *p++; + } +} + +/* Convert e-type X to IEEE double extended format E. */ + +static void etoe64 (x, e) unsigned EMUSHORT *x, *e; { @@ -2457,6 +3469,13 @@ etoe64 (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]; @@ -2473,7 +3492,9 @@ etoe64 (x, 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; @@ -2481,60 +3502,125 @@ toe64 (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; @@ -2543,8 +3629,31 @@ toe53 (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; { @@ -2552,6 +3661,13 @@ etoe53 (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); @@ -2568,6 +3684,8 @@ etoe53 (x, e) toe53 (xi, e); } +/* Convert exploded e-type X, that has already been rounded to + 53-bit precision, to IEEE double Y. */ static void toe53 (x, y) @@ -2576,10 +3694,17 @@ 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++) @@ -2587,33 +3712,38 @@ toe53 (x, y) 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; } @@ -2628,30 +3758,52 @@ toe53 (x, y) } i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */ *y |= (unsigned EMUSHORT) i; /* high order output already has sign bit set */ -#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; { @@ -2659,6 +3811,13 @@ etoe24 (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); @@ -2675,6 +3834,9 @@ etoe24 (x, e) toe24 (xi, e); } +/* Convert exploded e-type X, that has already been rounded to + float precision, to IEEE float Y. */ + static void toe24 (x, y) unsigned EMUSHORT *x, *y; @@ -2682,9 +3844,17 @@ toe24 (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; @@ -2694,32 +3864,36 @@ toe24 (x, y) *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; @@ -2737,30 +3911,30 @@ toe24 (x, y) 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; { @@ -2769,6 +3943,10 @@ ecmp (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); @@ -2808,8 +3986,6 @@ ecmp (a, b) return (0); /* equality */ - - diff: if (*(--p) > *(--q)) @@ -2818,15 +3994,9 @@ ecmp (a, b) 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; { @@ -2834,42 +4004,41 @@ eround (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 @@ -2877,31 +4046,33 @@ ltoe (lp, y) 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 @@ -2910,24 +4081,22 @@ ultoe (lp, y) } -/* -; Find long integer and fractional parts - -; long i; -; unsigned EMUSHORT x[NE], frac[NE]; -; xifrac (x, &i, frac); +/* Find signed HOST_WIDE_INT integer I and floating point fractional + part FRAC of e-type (packed internal format) floating point input X. + The integer output I has the sign of the input, except that + positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC. + The output e-type fraction FRAC is the positive fractional + part of abs (X). */ - The integer output has the sign of the input. The fraction is -the positive fractional part of abs (x). -*/ -void +static void eifrac (x, i, frac) unsigned EMUSHORT *x; - 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); @@ -2938,46 +4107,54 @@ eifrac (x, i, frac) 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; @@ -2990,24 +4167,19 @@ eifrac (x, i, frac) } -/* -; Find unsigned long integer and fractional parts - -; unsigned long i; -; unsigned EMUSHORT x[NE], frac[NE]; -; xifrac (x, &i, frac); +/* Find unsigned HOST_WIDE_INT integer I and floating point fractional part + FRAC of e-type X. A negative input yields integer output = 0 but + correct fraction. */ - A negative e type input yields integer output = 0 - but correct fraction. -*/ -void +static void euifrac (x, i, frac) unsigned EMUSHORT *x; - 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); @@ -3018,40 +4190,41 @@ euifrac (x, i, frac) 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; @@ -3065,15 +4238,9 @@ euifrac (x, i, frac) emovo (xi, frac); } +/* Shift the significand of exploded e-type X up or down by SC bits. */ - -/* -; Shift significand -; -; Shifts significand area up or down by the number of bits -; given by the variable sc. -*/ -int +static int eshift (x, sc) unsigned EMUSHORT *x; int sc; @@ -3136,15 +4303,10 @@ eshift (x, 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[]; { @@ -3162,9 +4324,9 @@ enormlz (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); } @@ -3210,16 +4372,73 @@ enormlz (x) 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 */ @@ -3253,8 +4472,12 @@ static unsigned EMUSHORT emtens[NTEN + 1][NE] = {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,}, {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */ }; +#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; @@ -3262,30 +4485,14 @@ e24toasc (x, string, ndigs) { unsigned EMUSHORT w[NI]; -#ifdef INFINITY -#ifdef IBMPC - if ((x[1] & 0x7f80) == 0x7f80) -#else - if ((x[0] & 0x7f80) == 0x7f80) -#endif - { -#ifdef IBMPC - if (x[1] & 0x8000) -#else - if (x[0] & 0x8000) -#endif - sprintf (string, " -Infinity "); - else - sprintf (string, " Infinity "); - return; - } -#endif e24toe (x, w); etoasc (w, string, ndigs); } +/* 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; @@ -3293,30 +4500,14 @@ e53toasc (x, string, ndigs) { unsigned EMUSHORT w[NI]; -#ifdef INFINITY -#ifdef IBMPC - if ((x[3] & 0x7ff0) == 0x7ff0) -#else - if ((x[0] & 0x7ff0) == 0x7ff0) -#endif - { -#ifdef IBMPC - if (x[3] & 0x8000) -#else - if (x[0] & 0x8000) -#endif - sprintf (string, " -Infinity "); - else - sprintf (string, " Infinity "); - return; - } -#endif e53toe (x, w); etoasc (w, string, ndigs); } +/* 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; @@ -3324,32 +4515,31 @@ e64toasc (x, string, ndigs) { unsigned EMUSHORT w[NI]; -#ifdef INFINITY -#ifdef IBMPC - if ((x[4] & 0x7fff) == 0x7fff) -#else - if ((x[0] & 0x7fff) == 0x7fff) -#endif - { -#ifdef IBMPC - if (x[4] & 0x8000) -#else - if (x[0] & 0x8000) -#endif - sprintf (string, " -Infinity "); - else - sprintf (string, " Infinity "); - return; - } -#endif e64toe (x, w); etoasc (w, string, ndigs); } +/* 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; @@ -3363,11 +4553,19 @@ etoasc (x, string, ndigs) 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) @@ -3390,12 +4588,11 @@ etoasc (x, string, ndigs) 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) @@ -3420,9 +4617,12 @@ etoasc (x, string, ndigs) 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; @@ -3452,6 +4652,7 @@ etoasc (x, string, ndigs) 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) @@ -3468,7 +4669,7 @@ etoasc (x, string, ndigs) } 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) @@ -3526,7 +4727,7 @@ etoasc (x, string, ndigs) ediv (t, eone, t); } isone: - /* Find the first (leading) digit. */ + /* Find the first (leading) digit. */ emovi (t, w); emovz (w, t); emovi (y, w); @@ -3549,7 +4750,7 @@ etoasc (x, string, ndigs) *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) @@ -3567,10 +4768,10 @@ etoasc (x, string, ndigs) } 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 */ @@ -3588,7 +4789,7 @@ etoasc (x, string, ndigs) /* 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); @@ -3645,29 +4846,16 @@ etoasc (x, string, ndigs) } +/* Convert ASCII string to floating point. + Numeric input is a free format decimal number of any length, with + or without decimal point. Entering E after the number followed by an + integer number causes the second number to be interpreted as a power of + 10 to be multiplied by the first number (i.e., "scientific" notation). */ -/* -; ASCTOQ -; ASCTOQ.MAC LATEST REV: 11 JAN 84 -; SLM, 3 JAN 78 -; -; Convert ASCII string to quadruple precision floating point -; -; Numeric input is free field decimal number -; with max of 15 digits with or without -; decimal point entered as ASCII from teletype. -; Entering E after the number followed by a second -; number causes the second number to be interpreted -; as a power of 10 to be multiplied by the first number -; (i.e., "scientific" notation). -; -; Usage: -; asctoq (string, q); -*/ - -/* ASCII to single */ -void +/* Convert ASCII string S to single precision float value Y. */ + +static void asctoe24 (s, y) char *s; unsigned EMUSHORT *y; @@ -3676,13 +4864,14 @@ asctoe24 (s, y) } -/* ASCII to double */ -void +/* Convert ASCII string S to double precision value Y. */ + +static void asctoe53 (s, y) char *s; unsigned EMUSHORT *y; { -#ifdef DEC +#if defined(DEC) || defined(IBM) asctoeg (s, y, 56); #else asctoeg (s, y, 53); @@ -3690,8 +4879,9 @@ asctoe53 (s, y) } -/* ASCII to long double */ -void +/* Convert ASCII string S to double extended value Y. */ + +static void asctoe64 (s, y) char *s; unsigned EMUSHORT *y; @@ -3699,8 +4889,19 @@ asctoe64 (s, 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; @@ -3708,10 +4909,10 @@ asctoe (s, 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; @@ -3722,19 +4923,16 @@ asctoeg (ss, y, oprec) 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; @@ -3756,7 +4954,7 @@ asctoeg (ss, y, oprec) /* 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; @@ -3775,11 +4973,12 @@ asctoeg (ss, y, oprec) 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) @@ -3795,7 +4994,11 @@ asctoeg (ss, y, oprec) } 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; @@ -3835,8 +5038,12 @@ asctoeg (ss, y, oprec) goto infinite; default: error: +#ifdef NANS + einan (yy); +#else mtherr ("asctoe", DOMAIN); - eclear (y); + eclear (yy); +#endif goto aexit; } donchr: @@ -3845,7 +5052,15 @@ asctoeg (ss, y, oprec) /* 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; @@ -3861,7 +5076,7 @@ asctoeg (ss, y, oprec) { exp *= 10; exp += *s++ - '0'; - if (exp > 4956) + if (exp > -(MINDECEXP)) { if (esign < 0) goto zero; @@ -3871,14 +5086,14 @@ asctoeg (ss, y, oprec) } 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); @@ -3887,7 +5102,7 @@ asctoeg (ss, y, oprec) 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); @@ -3907,15 +5122,15 @@ asctoeg (ss, y, oprec) } 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) { @@ -3928,7 +5143,8 @@ asctoeg (ss, y, oprec) 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); @@ -3967,8 +5183,13 @@ asctoeg (ss, y, oprec) /* 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; @@ -3987,6 +5208,11 @@ asctoeg (ss, y, oprec) todec (yy, y); /* see etodec.c */ break; #endif +#ifdef IBM + case 56: + toibm (yy, y, DFmode); + break; +#endif case 53: toe53 (yy, y); break; @@ -3996,6 +5222,9 @@ asctoeg (ss, y, oprec) case 64: toe64 (yy, y); break; + case 113: + toe113 (yy, y); + break; case NBITS: emovo (yy, y); break; @@ -4004,13 +5233,9 @@ asctoeg (ss, y, oprec) -/* y = largest integer not greater than x - * (truncated toward minus infinity) - * - * unsigned EMUSHORT x[NE], y[NE] - * - * efloor (x, y); - */ +/* Return Y = largest integer not greater than X (truncated toward minus + infinity). */ + static unsigned EMUSHORT bmask[] = { 0xffff, @@ -4032,7 +5257,7 @@ static unsigned EMUSHORT bmask[] = 0x0000, }; -void +static void efloor (x, y) unsigned EMUSHORT x[], y[]; { @@ -4079,16 +5304,10 @@ efloor (x, y) } -/* unsigned EMUSHORT x[], s[]; - * int *exp; - * - * efrexp (x, exp, s); - * - * Returns s and exp such that s * 2**exp = x and .5 <= s < 1. - * For example, 1.1 = 0.55 * 2**1 - * Handles denormalized numbers properly using long integer exp. - */ -void +/* Return S and EXP such that S * 2^EXP = X and .5 <= S < 1. + For example, 1.1 = 0.55 * 2^1. */ + +static void efrexp (x, exp, s) unsigned EMUSHORT x[]; int *exp; @@ -4098,6 +5317,7 @@ efrexp (x, exp, s) EMULONG li; emovi (x, xi); + /* Handle denormalized numbers properly using long integer exponent. */ li = (EMULONG) ((EMUSHORT) xi[1]); if (li == 0) @@ -4109,16 +5329,9 @@ efrexp (x, exp, s) *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; @@ -4137,15 +5350,25 @@ eldexp (x, pwr2, y) } -/* c = remainder after dividing b by a - * Least significant integer quotient bits left in equot[]. - */ -void +/* C = remainder after dividing B by A, all e type values. + Least significant integer quotient bits left in EQUOT. */ + +static void eremain (a, b, c) unsigned EMUSHORT a[], b[], c[]; { unsigned EMUSHORT den[NI], num[NI]; +#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); @@ -4163,7 +5386,10 @@ eremain (a, b, c) emovo (num, c); } -void +/* Return quotient of exploded e-types NUM / DEN in EQUOT, + remainder in NUM. */ + +static void eiremain (den, num) unsigned EMUSHORT den[], num[]; { @@ -4183,9 +5409,7 @@ eiremain (den, num) j = 1; } else - { j = 0; - } eshup1 (equot); equot[NI - 1] |= j; eshup1 (num); @@ -4194,69 +5418,26 @@ eiremain (den, 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 - * - */ - -/* -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 */ @@ -4264,51 +5445,37 @@ static char *ermsg[7] = "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; @@ -4346,88 +5513,9 @@ dectoe (d, 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; { @@ -4436,8 +5524,9 @@ etodec (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); @@ -4445,7 +5534,10 @@ etodec (x, d) todec (xi, d); } -void +/* Convert exploded e-type X, that has already been rounded to + 56-bit precision, to DEC format double Y. */ + +static void todec (x, y) unsigned EMUSHORT *x, *y; { @@ -4486,7 +5578,723 @@ todec (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 */ + +/* 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 (); + } +}