/* 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 Free Software Foundation, Inc.
+ Contributed by Stephen L. Moshier (moshier@world.std.com).
This file is part of GNU CC.
or Simon & Schuster Int'l, 1989. A library of XFmode elementary
transcendental functions can be obtained by ftp from
research.att.com: netlib/cephes/ldouble.shar.Z */
-
+\f
/* Type of computer arithmetic.
- * Only one of DEC, IBM, MIEEE, IBMPC, or UNK should get defined.
- */
+ Only one of DEC, IBM, MIEEE, IBMPC, or UNK should get defined.
-/* `MIEEE' refers generically to big-endian IEEE floating-point data
+ `MIEEE' refers generically to big-endian IEEE floating-point data
structure. This definition should work in SFmode `float' type and
DFmode `double' type on virtually all big-endian IEEE machines.
If LONG_DOUBLE_TYPE_SIZE has been defined to be 96, then MIEEE
The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
and may deactivate XFmode since `long double' is used to refer
- to both modes. */
+ to both modes.
+
+ The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
+ contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
+ separate the floating point unit's endian-ness from that of
+ the integer addressing. This permits one to define a big-endian
+ FPU on a little-endian machine (e.g., ARM). An extension to
+ BYTES_BIG_ENDIAN may be required for some machines in the future.
+ These optional macros may be defined in tm.h. In real.h, they
+ default to WORDS_BIG_ENDIAN, etc., so there is no need to define
+ them for any normal host or target machine on which the floats
+ and the integers have the same endian-ness. */
+
/* The following converts gcc macros into the ones used by this file. */
#define IBM 1
#else /* it's also not an IBM */
#if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
-#if WORDS_BIG_ENDIAN
+#if FLOAT_WORDS_BIG_ENDIAN
/* Motorola IEEE, high order words come first (Sun workstation): */
#define MIEEE 1
#else /* not big-endian */
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 */
#define IBM 1
#else /* it's also not an IBM */
#if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
-#ifdef HOST_WORDS_BIG_ENDIAN
+#if HOST_FLOAT_WORDS_BIG_ENDIAN
#define MIEEE 1
#else /* not big-endian */
#define IBMPC 1
#define INFINITY
#endif
#endif
-
+\f
/* Find a host integer type that is at least 16 bits wide,
and another type at least twice whatever that size is. */
/* 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))
+#if HOST_FLOAT_WORDS_BIG_ENDIAN == FLOAT_WORDS_BIG_ENDIAN
+#define GET_REAL(r,e) e53toe ((unsigned EMUSHORT*) (r), (e))
+#define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))
#else /* endian-ness differs */
/* emulator uses target endian-ness internally */
#define GET_REAL(r,e) \
-do { EMUSHORT w[4]; \
+do { unsigned EMUSHORT w[4]; \
w[3] = ((EMUSHORT *) r)[0]; \
w[2] = ((EMUSHORT *) r)[1]; \
w[1] = ((EMUSHORT *) r)[2]; \
e53toe (w, (e)); } while (0)
#define PUT_REAL(e,r) \
-do { EMUSHORT w[4]; \
+do { unsigned EMUSHORT w[4]; \
etoe53 ((e), w); \
*((EMUSHORT *) r) = w[3]; \
*((EMUSHORT *) r + 1) = w[2]; \
#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 */
/* The exponent of 1.0 */
#define EXONE (0x3fff)
-void warning ();
extern int extra_warnings;
-int ecmp (), enormlz (), eshift ();
-int eisneg (), eisinf (), eisnan (), eiisinf (), eiisnan ();
-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 ereal_to_decimal (), eiinfin (), einan ();
-void esqrt (), elog (), eexp (), etanh (), epow ();
-void asctoe (), asctoe24 (), asctoe53 (), asctoe64 (), asctoe113 ();
-void etoasc (), e24toasc (), e53toasc (), e64toasc (), e113toasc ();
-void etoe64 (), etoe53 (), etoe24 (), e64toe (), e53toe (), e24toe ();
-void etoe113 (), e113toe ();
-void mtherr (), make_nan ();
-void enan ();
extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
extern unsigned EMUSHORT elog2[], esqrt2[];
-/* Pack output array with 32-bit numbers obtained from
- array containing 16-bit numbers, swapping ends if required. */
-void
+static void endian PROTO((unsigned EMUSHORT *, long *,
+ enum machine_mode));
+static void eclear PROTO((unsigned EMUSHORT *));
+static void emov PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void eabs PROTO((unsigned EMUSHORT *));
+static void eneg PROTO((unsigned EMUSHORT *));
+static int eisneg PROTO((unsigned EMUSHORT *));
+static int eisinf PROTO((unsigned EMUSHORT *));
+static int eisnan PROTO((unsigned EMUSHORT *));
+static void einfin PROTO((unsigned EMUSHORT *));
+static void enan PROTO((unsigned EMUSHORT *, int));
+static void emovi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void emovo PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ecleaz PROTO((unsigned EMUSHORT *));
+static void ecleazs PROTO((unsigned EMUSHORT *));
+static void emovz PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void einan PROTO((unsigned EMUSHORT *));
+static int eiisnan PROTO((unsigned EMUSHORT *));
+static int eiisneg PROTO((unsigned EMUSHORT *));
+static void eiinfin PROTO((unsigned EMUSHORT *));
+static int eiisinf PROTO((unsigned EMUSHORT *));
+static int ecmpm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void eshdn1 PROTO((unsigned EMUSHORT *));
+static void eshup1 PROTO((unsigned EMUSHORT *));
+static void eshdn8 PROTO((unsigned EMUSHORT *));
+static void eshup8 PROTO((unsigned EMUSHORT *));
+static void eshup6 PROTO((unsigned EMUSHORT *));
+static void eshdn6 PROTO((unsigned EMUSHORT *));
+static void eaddm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));\f
+static void esubm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void m16m PROTO((unsigned int, unsigned short *,
+ unsigned short *));
+static int edivm PROTO((unsigned short *, unsigned short *));
+static int emulm PROTO((unsigned short *, unsigned short *));
+static void emdnorm PROTO((unsigned EMUSHORT *, int, int, EMULONG, int));
+static void esub PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void eadd PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void eadd1 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void ediv PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void emul PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void e53toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void e64toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void e113toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void e24toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void toe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static int ecmp PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void eround PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ltoe PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *));
+static void ultoe PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *));
+static void eifrac PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *,
+ unsigned EMUSHORT *));
+static void euifrac PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *,
+ unsigned EMUSHORT *));
+static int eshift PROTO((unsigned EMUSHORT *, int));
+static int enormlz PROTO((unsigned EMUSHORT *));
+static void e24toasc PROTO((unsigned EMUSHORT *, char *, int));
+static void e53toasc PROTO((unsigned EMUSHORT *, char *, int));
+static void e64toasc PROTO((unsigned EMUSHORT *, char *, int));
+static void e113toasc PROTO((unsigned EMUSHORT *, char *, int));
+static void etoasc PROTO((unsigned EMUSHORT *, char *, int));
+static void asctoe24 PROTO((char *, unsigned EMUSHORT *));
+static void asctoe53 PROTO((char *, unsigned EMUSHORT *));
+static void asctoe64 PROTO((char *, unsigned EMUSHORT *));
+static void asctoe113 PROTO((char *, unsigned EMUSHORT *));
+static void asctoe PROTO((char *, unsigned EMUSHORT *));
+static void asctoeg PROTO((char *, unsigned EMUSHORT *, int));
+static void efloor PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void efrexp PROTO((unsigned EMUSHORT *, int *,
+ unsigned EMUSHORT *));
+static void eldexp PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *));
+static void eremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ unsigned EMUSHORT *));
+static void eiremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void mtherr PROTO((char *, int));
+static void dectoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etodec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void todec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ibmtoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ enum machine_mode));
+static void etoibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ enum machine_mode));
+static void toibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
+ enum machine_mode));
+static void make_nan PROTO((unsigned EMUSHORT *, int, enum machine_mode));
+static void uditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void ditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etoudi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void etodi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+static void esqrt PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
+\f
+/* Copy 32-bit numbers obtained from array containing 16-bit numbers,
+ swapping ends if required, into output array of longs. The
+ result is normally passed to fprintf by the ASM_OUTPUT_ macros. */
+
+static void
endian (e, x, mode)
unsigned EMUSHORT e[];
long x[];
{
unsigned long th, t;
-#if WORDS_BIG_ENDIAN
+#if FLOAT_WORDS_BIG_ENDIAN
switch (mode)
{
}
-/* 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;
if (ecmp (d2, ezero) == 0)
{
#ifdef NANS
- enan (v);
+ enan (v, eisneg (d1) ^ eisneg (d2));
break;
#else
abort ();
}
-/* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT
- * implements REAL_VALUE_RNDZINT (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
}
-/* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT
- * implements REAL_VALUE_UNSIGNED_RNDZINT (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
}
-/* 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;
}
-/* Expansion of REAL_NEGATE.
- */
+/* Expansion of REAL_NEGATE. */
+
REAL_VALUE_TYPE
ereal_negate (x)
REAL_VALUE_TYPE x;
REAL_VALUE_TYPE r;
GET_REAL (&x, e);
-#ifdef NANS
- if (eisnan (e))
- return (x);
-#endif
eneg (e);
PUT_REAL (e, &r);
return (r);
}
-/* Round real toward zero to HOST_WIDE_INT
- * implements REAL_VALUE_FIX (x).
- */
-long
+/* 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];
- long l;
+ HOST_WIDE_INT l;
GET_REAL (&x, f);
#ifdef NANS
}
/* Round real toward zero to unsigned HOST_WIDE_INT
- * implements REAL_VALUE_UNSIGNED_FIX (x).
- * Negative input returns zero.
- */
-unsigned long
+ 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 long l;
+ unsigned HOST_WIDE_INT l;
GET_REAL (&x, f);
#ifdef NANS
}
-/* REAL_VALUE_FROM_INT macro.
- */
+/* REAL_VALUE_FROM_INT macro. */
+
void
ereal_from_int (d, i, j)
REAL_VALUE_TYPE *d;
- long i, j;
+ HOST_WIDE_INT i, j;
{
unsigned EMUSHORT df[NE], dg[NE];
- long low, high;
+ HOST_WIDE_INT low, high;
int sign;
sign = 0;
else
high += 1;
}
- 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);
}
-/* REAL_VALUE_FROM_UNSIGNED_INT macro.
- */
+/* REAL_VALUE_FROM_UNSIGNED_INT macro. */
+
void
ereal_from_uint (d, i, j)
REAL_VALUE_TYPE *d;
- unsigned long i, j;
+ unsigned HOST_WIDE_INT i, j;
{
unsigned EMUSHORT df[NE], dg[NE];
- unsigned long low, high;
+ unsigned HOST_WIDE_INT low, high;
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);
}
-/* 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];
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);
emul (df, dh, dg); /* fractional part is the low word */
}
-/* REAL_VALUE_LDEXP macro.
- */
+/* REAL_VALUE_LDEXP macro. */
+
REAL_VALUE_TYPE
ereal_ldexp (x, n)
REAL_VALUE_TYPE x;
}
/* 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. */
+
int
target_isinf (x)
REAL_VALUE_TYPE x;
/* Check for a negative REAL_VALUE_TYPE number.
- * this means strictly less than zero, not -0.
- */
+ 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) == -1)
- 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;
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);
#endif /* REAL_ARITHMETIC defined */
+/* Used for debugging--print the value of R in human-readable format
+ on stderr. */
+
+void
+debug_real (r)
+ REAL_VALUE_TYPE r;
+{
+ char dstr[30];
+
+ REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
+ fprintf (stderr, "%s", dstr);
+}
+
+\f
/* Target values are arrays of host longs. A long is guaranteed
to be at least 32 bits wide. */
/* 128-bit long double */
+
void
etartdouble (r, l)
REAL_VALUE_TYPE r;
}
/* 80-bit long double */
+
void
etarldouble (r, l)
REAL_VALUE_TYPE r;
}
/* 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
- * 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 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) 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) 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 (non-IEEE)
- * or is infinite (IEEE)
- * eisnan (e) 1 if e is a NaN
- *
- *
- * Routines for internal format numbers
- *
- * eaddm (ai, bi) add significands, bi = bi + ai
- * ecleaz (ei) ei = 0
- * ecleazs (ei) set ei = 0 but leave its sign alone
- * ecmpm (ai, bi) compare significands, return 1, 0, or -1
- * edivm (ai, bi) divide significands, bi = bi / ai
- * emdnorm (ai,l,s,exp) normalize and round off
- * emovi (a, ai) convert external a to internal ai
- * emovo (ai, a) convert internal ai to external a
- * emovz (ai, bi) bi = ai, low guard word of bi = 0
- * emulm (ai, bi) multiply significands, bi = bi * ai
- * enormlz (ei) left-justify the significand
- * eshdn1 (ai) shift significand and guards down 1 bit
- * eshdn8 (ai) shift down 8 bits
- * eshdn6 (ai) shift down 16 bits
- * eshift (ai, n) shift ai n bits up (or down if n < 0)
- * eshup1 (ai) shift significand and guards up 1 bit
- * eshup8 (ai) shift up 8 bits
- * eshup6 (ai) shift up 16 bits
- * esubm (ai, bi) subtract significands, bi = bi - ai
- * eiisinf (ai) 1 if infinite
- * eiisnan (ai) 1 if a NaN
- * 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).
- *
- */
-
-
-/* mconf.h
- *
- * Common include file for math routines
- *
- *
- *
- * SYNOPSIS:
- *
- * #include "mconf.h"
- *
- *
- *
- * DESCRIPTION:
- *
- * This file contains definitions for error codes that are
- * passed to the common error handling routine mtherr
- * (which see).
- *
- * The file also includes a conditional assembly definition
- * for the type of computer arithmetic (Intel IEEE, DEC, Motorola
- * IEEE, or UNKnown).
- *
- * For Digital Equipment PDP-11 and VAX computers, certain
- * IBM systems, and others that use numbers with a 56-bit
- * significand, the symbol DEC should be defined. In this
- * mode, most floating point constants are given as arrays
- * of octal integers to eliminate decimal to binary conversion
- * errors that might be introduced by the compiler.
- *
- * For computers, such as IBM PC, that follow the IEEE
- * Standard for Binary Floating Point Arithmetic (ANSI/IEEE
- * Std 754-1985), the symbol IBMPC or MIEEE should be defined.
- * These numbers have 53-bit significands. In this mode, constants
- * are provided as arrays of hexadecimal 16 bit integers.
- *
- * To accommodate other types of computer arithmetic, all
- * constants are also provided in a normal decimal radix
- * which one can hope are correctly converted to a suitable
- * format by the available C language compiler. To invoke
- * this mode, the symbol UNK is defined.
- *
- * An important difference among these modes is a predefined
- * set of machine arithmetic constants for each. The numbers
- * MACHEP (the machine roundoff error), MAXNUM (largest number
- * represented), and several other parameters are preset by
- * the configuration symbol. Check the file const.c to
- * ensure that these values are correct for your computer.
- *
- * For ANSI C compatibility, define ANSIC equal to 1. Currently
- * this affects only the atan2 function and others that use it.
- */
\f
+/*
+ Extended precision IEEE binary floating point arithmetic routines
+
+ Numbers are stored in C language as arrays of 16-bit unsigned
+ short integers. The arguments of the routines are pointers to
+ the arrays.
+
+ External e type data structure, simulates Intel 8087 chip
+ temporary real format but possibly with a larger significand:
+
+ NE-1 significand words (least significant word first,
+ most significant bit is normally set)
+ exponent (value = EXONE for 1.0,
+ top bit is the sign)
+
+
+ Internal data structure of a number (a "word" is 16 bits):
+
+ ei[0] sign word (0 for positive, 0xffff for negative)
+ ei[1] biased exponent (value = EXONE for the number 1.0)
+ ei[2] high guard word (always zero after normalization)
+ ei[3]
+ to ei[NI-2] significand (NI-4 significand words,
+ most significant word first,
+ most significant bit is set)
+ ei[NI-1] low guard word (0x8000 bit is rounding place)
+
+
+
+ Routines for external format numbers
+
+ asctoe (string, e) ASCII string to extended double e type
+ asctoe64 (string, &d) ASCII string to long double
+ asctoe53 (string, &d) ASCII string to double
+ asctoe24 (string, &f) ASCII string to single
+ asctoeg (string, e, prec) ASCII string to specified precision
+ e24toe (&f, e) IEEE single precision to e type
+ e53toe (&d, e) IEEE double precision to e type
+ e64toe (&d, e) IEEE long double precision to e type
+ e113toe (&d, e) 128-bit long double precision to e type
+ eabs (e) absolute value
+ eadd (a, b, c) c = b + a
+ eclear (e) e = 0
+ ecmp (a, b) Returns 1 if a > b, 0 if a == b,
+ -1 if a < b, -2 if either a or b is a NaN.
+ ediv (a, b, c) c = b / a
+ efloor (a, b) truncate to integer, toward -infinity
+ efrexp (a, exp, s) extract exponent and significand
+ eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
+ euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
+ einfin (e) set e to infinity, leaving its sign alone
+ eldexp (a, n, b) multiply by 2**n
+ emov (a, b) b = a
+ emul (a, b, c) c = b * a
+ eneg (e) e = -e
+ eround (a, b) b = nearest integer value to a
+ esub (a, b, c) c = b - a
+ e24toasc (&f, str, n) single to ASCII string, n digits after decimal
+ e53toasc (&d, str, n) double to ASCII string, n digits after decimal
+ e64toasc (&d, str, n) 80-bit long double to ASCII string
+ e113toasc (&d, str, n) 128-bit long double to ASCII string
+ etoasc (e, str, n) e to ASCII string, n digits after decimal
+ etoe24 (e, &f) convert e type to IEEE single precision
+ etoe53 (e, &d) convert e type to IEEE double precision
+ etoe64 (e, &d) convert e type to IEEE long double precision
+ ltoe (&l, e) HOST_WIDE_INT to e type
+ ultoe (&l, e) unsigned HOST_WIDE_INT to e type
+ eisneg (e) 1 if sign bit of e != 0, else 0
+ eisinf (e) 1 if e has maximum exponent (non-IEEE)
+ or is infinite (IEEE)
+ eisnan (e) 1 if e is a NaN
+
+
+ Routines for internal format numbers
+
+ eaddm (ai, bi) add significands, bi = bi + ai
+ ecleaz (ei) ei = 0
+ ecleazs (ei) set ei = 0 but leave its sign alone
+ ecmpm (ai, bi) compare significands, return 1, 0, or -1
+ edivm (ai, bi) divide significands, bi = bi / ai
+ emdnorm (ai,l,s,exp) normalize and round off
+ emovi (a, ai) convert external a to internal ai
+ emovo (ai, a) convert internal ai to external a
+ emovz (ai, bi) bi = ai, low guard word of bi = 0
+ emulm (ai, bi) multiply significands, bi = bi * ai
+ enormlz (ei) left-justify the significand
+ eshdn1 (ai) shift significand and guards down 1 bit
+ eshdn8 (ai) shift down 8 bits
+ eshdn6 (ai) shift down 16 bits
+ eshift (ai, n) shift ai n bits up (or down if n < 0)
+ eshup1 (ai) shift significand and guards up 1 bit
+ eshup8 (ai) shift up 8 bits
+ eshup6 (ai) shift up 16 bits
+ esubm (ai, bi) subtract significands, bi = bi - ai
+ eiisinf (ai) 1 if infinite
+ eiisnan (ai) 1 if a NaN
+ eiisneg (ai) 1 if sign bit of ai != 0, else 0
+ einan (ai) set ai = NaN
+ eiinfin (ai) set ai = infinity
+
+ The result is always normalized and rounded to NI-4 word precision
+ after each arithmetic operation.
+
+ Exception flags are NOT fully supported.
+
+ Signaling NaN's are NOT supported; they are treated the same
+ as quiet NaN's.
+
+ Define INFINITY for support of infinity; otherwise a
+ saturation arithmetic is implemented.
+
+ Define NANS for support of Not-a-Number items; otherwise the
+ arithmetic will never produce a NaN output, and might be confused
+ by a NaN input.
+ If NaN's are supported, the output of `ecmp (a,b)' is -2 if
+ either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
+ may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
+ if in doubt.
+
+ Denormals are always supported here where appropriate (e.g., not
+ for conversion to DEC numbers). */
+
+/* Definitions for error codes that are passed to the common error handling
+ routine mtherr.
+
+ For Digital Equipment PDP-11 and VAX computers, certain
+ IBM systems, and others that use numbers with a 56-bit
+ significand, the symbol DEC should be defined. In this
+ mode, most floating point constants are given as arrays
+ of octal integers to eliminate decimal to binary conversion
+ errors that might be introduced by the compiler.
+
+ For computers, such as IBM PC, that follow the IEEE
+ Standard for Binary Floating Point Arithmetic (ANSI/IEEE
+ Std 754-1985), the symbol IBMPC or MIEEE should be defined.
+ These numbers have 53-bit significands. In this mode, constants
+ are provided as arrays of hexadecimal 16 bit integers.
+
+ To accommodate other types of computer arithmetic, all
+ constants are also provided in a normal decimal radix
+ which one can hope are correctly converted to a suitable
+ format by the available C language compiler. To invoke
+ this mode, the symbol UNK is defined.
+
+ An important difference among these modes is a predefined
+ set of machine arithmetic constants for each. The numbers
+ MACHEP (the machine roundoff error), MAXNUM (largest number
+ represented), and several other parameters are preset by
+ the configuration symbol. Check the file const.c to
+ ensure that these values are correct for your computer.
+
+ For ANSI C compatibility, define ANSIC equal to 1. Currently
+ this affects only the atan2 function and others that use it. */
+
/* Constant definitions for math error conditions. */
#define DOMAIN 1 /* argument domain error */
/* Control register for rounding precision.
- * This can be set to 113 (if NE=10), 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 (), toe113 ();
-void eremain (), einit (), eiremain ();
-int ecmpm (), edivm (), emulm ();
-void emovi (), emovo (), emovz (), ecleaz (), ecleazs (), eadd1 ();
-#ifdef DEC
-void etodec (), todec (), dectoe ();
-#endif
-#ifdef IBM
-void etoibm (), toibm (), ibmtoe ();
-#endif
+/* Clear out entire external format number. */
-
-void
-einit ()
-{
-}
-
-/*
-; Clear out entire external format number.
-;
-; unsigned EMUSHORT x[];
-; eclear (x);
-*/
-
-void
+static void
eclear (x)
register unsigned EMUSHORT *x;
{
-/* Move external format number from a to b.
- *
- * emov (a, b);
- */
+/* Move external format number from a to b. */
-void
+static void
emov (a, b)
register unsigned EMUSHORT *a, *b;
{
}
-/*
-; Absolute value of external format number
-;
-; EMUSHORT x[NE];
-; eabs (x);
-*/
+/* Absolute value of external format number. */
-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 external format number. */
-
-
-/*
-; Negate external format number
-;
-; unsigned EMUSHORT x[NE];
-; eneg (x);
-*/
-
-void
+static void
eneg (x)
unsigned EMUSHORT x[];
{
-#ifdef NANS
- if (eisnan (x))
- return;
-#endif
x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
}
-/* Return 1 if external format number is negative,
- * else return zero, including when it is a NaN.
- */
-int
+/* Return 1 if sign bit of external format number is nonzero, else zero. */
+
+static int
eisneg (x)
unsigned EMUSHORT x[];
{
-#ifdef NANS
- if (eisnan (x))
- return (0);
-#endif
if (x[NE - 1] & 0x8000)
return (1);
else
}
-/* Return 1 if external format number is infinity.
- * else return zero.
- */
-int
+/* Return 1 if external format number is infinity, else return zero. */
+
+static int
eisinf (x)
unsigned EMUSHORT x[];
{
}
-/* 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. */
+/* 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. */
-int
+static int
eisnan (x)
unsigned EMUSHORT x[];
{
-
#ifdef NANS
int i;
-/* NaN has maximum exponent */
+
+ /* NaN has maximum exponent */
if ((x[NE - 1] & 0x7fff) != 0x7fff)
return (0);
-/* ... and non-zero significand field. */
+ /* ... and non-zero significand field. */
for (i = 0; i < NE - 1; i++)
{
if (*x++ != 0)
return (1);
}
#endif
+
return (0);
}
/* Fill external format number with infinity pattern (IEEE)
or largest possible number (non-IEEE). */
-void
+static void
einfin (x)
register unsigned EMUSHORT *x;
{
This generates Intel's quiet NaN pattern for extended real.
The exponent is 7fff, the leading mantissa word is c000. */
-void
-enan (x)
+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 = 0x7fff;
+ *x = (sign << 15) | 0x7fff;
}
-/* Move in external format number,
- * converting it to internal format.
- */
-void
+/* Move in external format number, converting it to internal format. */
+
+static void
emovi (a, b)
unsigned EMUSHORT *a, *b;
{
return;
}
#endif
+
for (i = 2; i < NI; i++)
*q++ = 0;
return;
}
#endif
+
/* clear high guard word */
*q++ = 0;
/* move in the significand */
}
-/* Move internal format number out,
- * converting it to external format.
- */
-void
+/* Move internal format number out, converting it to external format. */
+
+static void
emovo (a, b)
unsigned EMUSHORT *a, *b;
{
register unsigned EMUSHORT *p, *q;
unsigned EMUSHORT i;
+ int j;
p = a;
q = b + (NE - 1); /* point to output exponent */
#ifdef NANS
if (eiisnan (a))
{
- enan (b);
+ enan (b, eiisneg (a));
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 internal format number. */
-
-
-/* Clear out internal format number.
- */
-
-void
+static void
ecleaz (xi)
register unsigned EMUSHORT *xi;
{
}
-/* same, but don't touch the sign. */
+/* Same, but don't touch the sign. */
-void
+static void
ecleazs (xi)
register unsigned EMUSHORT *xi;
{
-/* Move internal format number from a to b.
- */
-void
+/* Move internal format number from a to b. */
+
+static void
emovz (a, b)
register unsigned EMUSHORT *a, *b;
{
/* Generate internal format NaN.
The explicit pattern for this is maximum exponent and
- top two significand bits set. */
+ top two significant bits set. */
-void
+static void
einan (x)
unsigned EMUSHORT x[];
{
/* Return nonzero if internal format number is a NaN. */
-int
+static int
eiisnan (x)
unsigned EMUSHORT x[];
{
return (0);
}
+/* Return nonzero if sign of internal format number is nonzero. */
+
+static int
+eiisneg (x)
+ unsigned EMUSHORT x[];
+{
+
+ return x[0] != 0;
+}
+
/* Fill internal format number with infinity pattern.
This has maximum exponent and significand all zeros. */
-void
+static void
eiinfin (x)
unsigned EMUSHORT x[];
{
/* Return nonzero if internal format number is infinite. */
-int
+static int
eiisinf (x)
unsigned EMUSHORT x[];
{
}
-/*
-; 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
+/* Compare significands of numbers in internal format.
+ Guard words are included in the comparison.
+
+ Returns +1 if a > b
+ 0 if a == b
+ -1 if a < b */
+
+static int
ecmpm (a, b)
register unsigned EMUSHORT *a, *b;
{
}
-/*
-; Shift significand down by 1 bit
-*/
+/* Shift significand down by 1 bit. */
-void
+static void
eshdn1 (x)
register unsigned EMUSHORT *x;
{
-/*
-; Shift significand up by 1 bit
-*/
+/* Shift significand up by 1 bit. */
-void
+static void
eshup1 (x)
register unsigned EMUSHORT *x;
{
}
+/* Shift significand down by 8 bits. */
-/*
-; Shift significand down by 8 bits
-*/
-
-void
+static void
eshdn8 (x)
register unsigned EMUSHORT *x;
{
}
}
-/*
-; Shift significand up by 8 bits
-*/
+/* Shift significand up by 8 bits. */
-void
+static void
eshup8 (x)
register unsigned EMUSHORT *x;
{
}
}
-/*
-; Shift significand up by 16 bits
-*/
+/* Shift significand up by 16 bits. */
-void
+static void
eshup6 (x)
register unsigned EMUSHORT *x;
{
*p = 0;
}
-/*
-; Shift significand down by 16 bits
-*/
+/* Shift significand down by 16 bits. */
-void
+static void
eshdn6 (x)
register unsigned EMUSHORT *x;
{
*(--p) = 0;
}
\f
-/*
-; Add significands
-; x + y replaces y
-*/
+/* Add significands. x + y replaces y. */
-void
+static void
eaddm (x, y)
unsigned EMUSHORT *x, *y;
{
}
}
-/*
-; Subtract significands
-; y - x replaces y
-*/
+/* Subtract significands. y - x replaces y. */
-void
+static void
esubm (x, y)
unsigned EMUSHORT *x, *y;
{
*p++ = 0;
}
- /* Use faster compare and subtraction if denominator
- * has only 15 bits of significance.
- */
+ /* Use faster compare and subtraction if denominator has only 15 bits of
+ significance. */
+
p = &den[M + 2];
if (*p++ == 0)
{
goto divdon;
}
- /* The number of quotient bits to calculate is
- * NBITS + 1 scaling guard bit + 1 roundoff bit.
- */
+ /* The number of quotient bits to calculate is NBITS + 1 scaling guard
+ bit + 1 roundoff bit. */
+
fulldiv:
p = &equot[NI - 2];
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;
/* Multiply significand of e-type number b
-by 16-bit quantity a, e-type result to c. */
+ by 16-bit quantity a, e-type result to c. */
-void m16m( a, b, c )
-unsigned short a;
-unsigned short b[], c[];
+static void
+m16m (a, b, c)
+ unsigned int a;
+ unsigned short b[], c[];
{
-register unsigned short *pp;
-register unsigned long carry;
-unsigned short *ps;
-unsigned short p[NI];
-unsigned long aa, m;
-int i;
-
-aa = a;
-pp = &p[NI-2];
-*pp++ = 0;
-*pp = 0;
-ps = &b[NI-1];
-
-for( i=M+1; i<NI; i++ )
+ register unsigned short *pp;
+ register unsigned long carry;
+ unsigned short *ps;
+ unsigned short p[NI];
+ unsigned long aa, m;
+ int i;
+
+ aa = a;
+ pp = &p[NI-2];
+ *pp++ = 0;
+ *pp = 0;
+ ps = &b[NI-1];
+
+ for (i=M+1; i<NI; i++)
+ {
+ if (*ps == 0)
{
- if( *ps == 0 )
- {
- --ps;
- --pp;
- *(pp-1) = 0;
- }
- else
- {
- m = (unsigned long) aa * *ps--;
- carry = (m & 0xffff) + *pp;
- *pp-- = (unsigned short )carry;
- carry = (carry >> 16) + (m >> 16) + *pp;
- *pp = (unsigned short )carry;
- *(pp-1) = carry >> 16;
- }
+ --ps;
+ --pp;
+ *(pp-1) = 0;
}
-for( i=M; i<NI; i++ )
- c[i] = p[i];
+ else
+ {
+ m = (unsigned long) aa * *ps--;
+ carry = (m & 0xffff) + *pp;
+ *pp-- = (unsigned short)carry;
+ carry = (carry >> 16) + (m >> 16) + *pp;
+ *pp = (unsigned short)carry;
+ *(pp-1) = carry >> 16;
+ }
+ }
+ for (i=M; i<NI; i++)
+ c[i] = p[i];
}
/* Divide significands. Neither the numerator nor the denominator
-is permitted to have its high guard word nonzero. */
-
+ is permitted to have its high guard word nonzero. */
-int edivm( den, num )
-unsigned short den[], num[];
+static int
+edivm (den, num)
+ unsigned short den[], num[];
{
-int i;
-register unsigned short *p;
-unsigned long tnum;
-unsigned short j, tdenm, tquot;
-unsigned short tprod[NI+1];
+ int i;
+ register unsigned short *p;
+ unsigned long tnum;
+ unsigned short j, tdenm, tquot;
+ unsigned short tprod[NI+1];
-p = &equot[0];
-*p++ = num[0];
-*p++ = num[1];
+ p = &equot[0];
+ *p++ = num[0];
+ *p++ = num[1];
-for( i=M; i<NI; i++ )
- {
- *p++ = 0;
- }
-eshdn1( num );
-tdenm = den[M+1];
-for( i=M; i<NI; i++ )
- {
- /* Find trial quotient digit (the radix is 65536). */
- tnum = (((unsigned long) num[M]) << 16) + num[M+1];
-
- /* Do not execute the divide instruction if it will overflow. */
- if( (tdenm * 0xffffL) < tnum )
- tquot = 0xffff;
- else
- tquot = tnum / tdenm;
- /* Multiply denominator by trial quotient digit. */
- m16m( tquot, den, tprod );
- /* The quotient digit may have been overestimated. */
- if( ecmpm( tprod, num ) > 0 )
- {
- tquot -= 1;
- esubm( den, tprod );
- if( ecmpm( tprod, num ) > 0 )
- {
- tquot -= 1;
- esubm( den, tprod );
- }
- }
- esubm( tprod, num );
- equot[i] = tquot;
- eshup6(num);
- }
-/* test for nonzero remainder after roundoff bit */
-p = &num[M];
-j = 0;
-for( i=M; i<NI; i++ )
+ for (i=M; i<NI; i++)
+ {
+ *p++ = 0;
+ }
+ eshdn1 (num);
+ tdenm = den[M+1];
+ for (i=M; i<NI; i++)
+ {
+ /* Find trial quotient digit (the radix is 65536). */
+ tnum = (((unsigned long) num[M]) << 16) + num[M+1];
+
+ /* Do not execute the divide instruction if it will overflow. */
+ if ((tdenm * 0xffffL) < tnum)
+ tquot = 0xffff;
+ else
+ tquot = tnum / tdenm;
+ /* Multiply denominator by trial quotient digit. */
+ m16m ((unsigned int)tquot, den, tprod);
+ /* The quotient digit may have been overestimated. */
+ if (ecmpm (tprod, num) > 0)
{
- j |= *p++;
+ tquot -= 1;
+ esubm (den, tprod);
+ if (ecmpm (tprod, num) > 0)
+ {
+ tquot -= 1;
+ esubm (den, tprod);
+ }
}
-if( j )
- j = 1;
+ esubm (tprod, num);
+ equot[i] = tquot;
+ eshup6(num);
+ }
+ /* test for nonzero remainder after roundoff bit */
+ p = &num[M];
+ j = 0;
+ for (i=M; i<NI; i++)
+ {
+ j |= *p++;
+ }
+ if (j)
+ j = 1;
-for( i=0; i<NI; i++ )
- num[i] = equot[i];
+ for (i=0; i<NI; i++)
+ num[i] = equot[i];
-return( (int )j );
+ return ((int)j);
}
/* Multiply significands */
-int emulm( a, b )
-unsigned short a[], b[];
+static int
+emulm (a, b)
+ unsigned short a[], b[];
{
-unsigned short *p, *q;
-unsigned short pprod[NI];
-unsigned short j;
-int i;
-
-equot[0] = b[0];
-equot[1] = b[1];
-for( i=M; i<NI; i++ )
- equot[i] = 0;
-
-j = 0;
-p = &a[NI-1];
-q = &equot[NI-1];
-for( i=M+1; i<NI; i++ )
+ unsigned short *p, *q;
+ unsigned short pprod[NI];
+ unsigned short j;
+ int i;
+
+ equot[0] = b[0];
+ equot[1] = b[1];
+ for (i=M; i<NI; i++)
+ equot[i] = 0;
+
+ j = 0;
+ p = &a[NI-1];
+ q = &equot[NI-1];
+ for (i=M+1; i<NI; i++)
+ {
+ if (*p == 0)
{
- if( *p == 0 )
- {
- --p;
- }
- else
- {
- m16m( *p--, b, pprod );
- eaddm(pprod, equot);
- }
- j |= *q;
- eshdn6(equot);
+ --p;
+ }
+ else
+ {
+ m16m ((unsigned int) *p--, b, pprod);
+ eaddm(pprod, equot);
}
+ j |= *q;
+ eshdn6(equot);
+ }
-for( i=0; i<NI; i++ )
- b[i] = equot[i];
+ for (i=0; i<NI; i++)
+ b[i] = equot[i];
-/* return flag for lost nonzero bits */
-return( (int)j );
+ /* return flag for lost nonzero bits */
+ return ((int)j);
}
#endif
-/*
- * Normalize and round off.
- *
- * The internal format number to be rounded is "s".
- * Input "lost" indicates whether or not the number is exact.
- * This is the so-called sticky bit.
- *
- * Input "subflg" indicates whether the number was obtained
- * by a subtraction operation. In that case if lost is nonzero
- * then the number is slightly smaller than indicated.
- *
- * Input "exp" is the biased exponent, which may be negative.
- * the exponent field of "s" is ignored but is replaced by
- * "exp" as adjusted by normalization and rounding.
- *
- * Input "rcntrl" is the rounding control.
- */
+/* Normalize and round off.
-/* For future reference: In order for emdnorm to round off denormal
+ The internal format number to be rounded is "s".
+ Input "lost" indicates whether or not the number is exact.
+ This is the so-called sticky bit.
+
+ Input "subflg" indicates whether the number was obtained
+ by a subtraction operation. In that case if lost is nonzero
+ then the number is slightly smaller than indicated.
+
+ Input "exp" is the biased exponent, which may be negative.
+ the exponent field of "s" is ignored but is replaced by
+ "exp" as adjusted by normalization and rounding.
+
+ Input "rcntrl" is the rounding control.
+
+ For future reference: In order for emdnorm to round off denormal
significands at the right point, the input exponent must be
adjusted to be the actual value it would have after conversion to
the final floating point type. This adjustment has been
static int re = 0;
static unsigned EMUSHORT rbit[NI];
-void
+static void
emdnorm (s, lost, subflg, exp, rcntrl)
unsigned EMUSHORT s[];
int lost;
-/*
-; Subtract external format numbers.
-;
-; unsigned EMUSHORT a[NE], b[NE], c[NE];
-; esub (a, b, c); c = b - a
-*/
+/* Subtract external format numbers. */
static int subflg = 0;
-void
+static void
esub (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
&& ((eisneg (a) ^ eisneg (b)) == 0))
{
mtherr ("esub", INVALID);
- enan (c);
+ enan (c, 0);
return;
}
#endif
}
-/*
-; Add.
-;
-; unsigned EMUSHORT a[NE], b[NE], c[NE];
-; eadd (a, b, c); c = b + a
-*/
-void
+/* Add. */
+
+static void
eadd (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
&& ((eisneg (a) ^ eisneg (b)) != 0))
{
mtherr ("esub", INVALID);
- enan (c);
+ enan (c, 0);
return;
}
#endif
eadd1 (a, b, c);
}
-void
+static void
eadd1 (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
-/*
-; Divide.
-;
-; unsigned EMUSHORT a[NE], b[NE], c[NE];
-; ediv (a, b, c); c = b / a
-*/
-void
+/* Divide. */
+
+static void
ediv (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
|| (eisinf (a) && eisinf (b)))
{
mtherr ("ediv", INVALID);
- enan (c);
+ enan (c, eisneg (a) ^ eisneg (b));
return;
}
#endif
-/*
-; Multiply.
-;
-; unsigned EMUSHORT a[NE], b[NE], c[NE];
-; emul (a, b, c); c = b * a
-*/
-void
+/* Multiply. */
+
+static void
emul (a, b, c)
unsigned EMUSHORT *a, *b, *c;
{
|| (eisinf (b) && (ecmp (a, ezero) == 0)))
{
mtherr ("emul", INVALID);
- enan (c);
+ enan (c, eisneg (a) ^ eisneg (b));
return;
}
#endif
-/*
-; Convert IEEE double precision to e type
-; double d;
-; unsigned EMUSHORT x[N+2];
-; e53toe (&d, x);
-*/
-void
+/* Convert IEEE double precision to e type. */
+
+static void
e53toe (pe, y)
unsigned EMUSHORT *pe, *y;
{
if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
|| (pe[1] != 0) || (pe[0] != 0))
{
- enan (y);
+ enan (y, yy[0] != 0);
return;
}
#else
if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
|| (pe[2] != 0) || (pe[3] != 0))
{
- enan (y);
+ enan (y, yy[0] != 0);
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;
#endif /* not DEC */
}
-void
+static void
e64toe (pe, y)
unsigned EMUSHORT *pe, *y;
{
{
if (pe[i] != 0)
{
- enan (y);
+ enan (y, (*p & 0x8000) != 0);
return;
}
}
{
if (pe[i] != 0)
{
- enan (y);
+ enan (y, (*p & 0x8000) != 0);
return;
}
}
}
-void
+static void
e113toe (pe, y)
unsigned EMUSHORT *pe, *y;
{
{
if (pe[i] != 0)
{
- enan (y);
+ enan (y, yy[0] != 0);
return;
}
}
{
if (pe[i] != 0)
{
- enan (y);
+ enan (y, yy[0] != 0);
return;
}
}
}
-/*
-; Convert IEEE single precision to e type
-; float d;
-; unsigned EMUSHORT x[N+2];
-; dtox (&d, x);
-*/
-void
+/* Convert IEEE single precision to e type. */
+
+static void
e24toe (pe, y)
unsigned EMUSHORT *pe, *y;
{
#ifdef MIEEE
if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
{
- enan (y);
+ enan (y, yy[0] != 0);
return;
}
#else
if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
{
- enan (y);
+ enan (y, yy[0] != 0);
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;
}
-void
+static void
etoe113 (x, e)
unsigned EMUSHORT *x, *e;
{
#ifdef NANS
if (eisnan (x))
{
- make_nan (e, TFmode);
+ make_nan (e, eisneg (x), TFmode);
return;
}
#endif
toe113 (xi, e);
}
-/* move out internal format to ieee long double */
+/* Move out internal format to ieee long double */
+
static void
toe113 (a, b)
unsigned EMUSHORT *a, *b;
#ifdef NANS
if (eiisnan (a))
{
- make_nan (b, TFmode);
+ make_nan (b, eiisneg (a), TFmode);
return;
}
#endif
#endif
}
-void
+static void
etoe64 (x, e)
unsigned EMUSHORT *x, *e;
{
#ifdef NANS
if (eisnan (x))
{
- make_nan (e, XFmode);
+ make_nan (e, eisneg (x), XFmode);
return;
}
#endif
toe64 (xi, e);
}
-/* move out internal format to ieee long double */
+
+/* Move out internal format to ieee long double. */
+
static void
toe64 (a, b)
unsigned EMUSHORT *a, *b;
#ifdef NANS
if (eiisnan (a))
{
- make_nan (b, XFmode);
+ make_nan (b, eiisneg (a), XFmode);
return;
}
#endif
}
-/*
-; e type to IEEE double precision
-; double d;
-; unsigned EMUSHORT x[NE];
-; etoe53 (x, &d);
-*/
+/* e type to IEEE double precision. */
#ifdef DEC
-void
+static void
etoe53 (x, e)
unsigned EMUSHORT *x, *e;
{
#else
#ifdef IBM
-void
+static void
etoe53 (x, e)
unsigned EMUSHORT *x, *e;
{
#else /* it's neither DEC nor IBM */
-void
+static void
etoe53 (x, e)
unsigned EMUSHORT *x, *e;
{
#ifdef NANS
if (eisnan (x))
{
- make_nan (e, DFmode);
+ make_nan (e, eisneg (x), DFmode);
return;
}
#endif
#ifdef NANS
if (eiisnan (x))
{
- make_nan (y, DFmode);
+ make_nan (y, eiisneg (x), DFmode);
return;
}
#endif
-/*
-; e type to IEEE single precision
-; float d;
-; unsigned EMUSHORT x[N+2];
-; xtod (x, &d);
-*/
+/* e type to IEEE single precision. */
+
#ifdef IBM
-void
+static void
etoe24 (x, e)
unsigned EMUSHORT *x, *e;
{
#else
-void
+static void
etoe24 (x, e)
unsigned EMUSHORT *x, *e;
{
#ifdef NANS
if (eisnan (x))
{
- make_nan (e, SFmode);
+ make_nan (e, eisneg (x), SFmode);
return;
}
#endif
#ifdef NANS
if (eiisnan (x))
{
- make_nan (y, SFmode);
+ make_nan (y, eiisneg (x), SFmode);
return;
}
#endif
}
#endif /* not IBM */
-/* 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
- * -2 if either a or b is a NaN.
- */
-int
+/* Compare two e type numbers.
+ Return +1 if a > b
+ 0 if a == b
+ -1 if a < b
+ -2 if either a or b is a NaN. */
+
+static int
ecmp (a, b)
unsigned EMUSHORT *a, *b;
{
-/* Find nearest integer to x = floor (x + 0.5)
- *
- * unsigned EMUSHORT x[NE], y[NE]
- * eround (x, y);
- */
-void
+/* Find nearest integer to x = floor (x + 0.5). */
+
+static void
eround (x, y)
unsigned EMUSHORT *x, *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
+/* Convert HOST_WIDE_INT to e type. */
+
+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_LONG == 64
+#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);
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 to e type. */
+
+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_LONG == 64
+#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);
}
-/*
-; Find long integer and fractional parts
-
-; long i;
-; unsigned EMUSHORT x[NE], frac[NE];
-; xifrac (x, &i, frac);
+/* Find signed HOST_WIDE_INT integer and floating point fractional
+ parts of e-type (packed internal format) floating point input X.
+ The integer output I has the sign of the input, except that
+ positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
+ The output e-type fraction FRAC is the positive fractional
+ part of abs (X). */
- 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 j, k;
- unsigned long ll;
+ unsigned HOST_WIDE_INT ll;
emovi (x, xi);
k = (int) xi[E] - (EXONE - 1);
emovo (xi, frac);
return;
}
- if (k > (HOST_BITS_PER_LONG - 1))
+ if (k > (HOST_BITS_PER_WIDE_INT - 1))
{
/* long integer overflow: output large integer
and correct fraction */
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");
{
/* shift not more than 16 bits */
eshift (xi, k);
- *i = (long) xi[M] & 0xffff;
+ *i = (HOST_WIDE_INT) xi[M] & 0xffff;
if (xi[0])
*i = -(*i);
}
}
-/* Find unsigned long integer and fractional parts.
+/* Find unsigned HOST_WIDE_INT integer and floating point fractional parts.
A negative e type input yields integer output = 0
but correct fraction. */
-void
+static void
euifrac (x, i, frac)
unsigned EMUSHORT *x;
- unsigned long *i;
+ unsigned HOST_WIDE_INT *i;
unsigned EMUSHORT *frac;
{
- unsigned long ll;
+ unsigned HOST_WIDE_INT ll;
unsigned EMUSHORT xi[NI];
int j, k;
emovo (xi, frac);
return;
}
- if (k > HOST_BITS_PER_LONG)
+ if (k > HOST_BITS_PER_WIDE_INT)
{
/* Long integer overflow: output large integer
and correct fraction.
{
/* shift not more than 16 bits */
eshift (xi, k);
- *i = (long) xi[M] & 0xffff;
+ *i = (HOST_WIDE_INT) xi[M] & 0xffff;
}
if (xi[0]) /* A negative value yields unsigned integer 0. */
-/*
-; Shift significand
-;
-; Shifts significand area up or down by the number of bits
-; given by the variable sc.
-*/
-int
+/* Shift significand area up or down by the number of bits given by SC. */
+
+static int
eshift (x, sc)
unsigned EMUSHORT *x;
int sc;
-/*
-; normalize
-;
-; Shift normalizes the significand area pointed to by argument
-; shift count (up = positive) is returned.
-*/
-int
+/* Shift normalize the significand area pointed to by argument.
+ Shift count (up = positive) is returned. */
+
+static int
enormlz (x)
unsigned EMUSHORT x[];
{
{
eshup6 (x);
sc += 16;
+
/* With guard word, there are NBITS+16 bits available.
- * return true if all are zero.
- */
+ Return true if all are zero. */
if (sc > NBITS)
return (sc);
}
/* Convert e type number to decimal format ASCII string.
- * The constants are for 64 bit precision.
- */
+ The constants are for 64 bit precision. */
#define NTEN 12
#define MAXP 4096
};
#endif
-void
+static void
e24toasc (x, string, ndigs)
unsigned EMUSHORT x[];
char *string;
}
-void
+static void
e53toasc (x, string, ndigs)
unsigned EMUSHORT x[];
char *string;
}
-void
+static void
e64toasc (x, string, ndigs)
unsigned EMUSHORT x[];
char *string;
etoasc (w, string, ndigs);
}
-void
+static void
e113toasc (x, string, ndigs)
unsigned EMUSHORT x[];
char *string;
static char wstring[80]; /* working storage for ASCII output */
-void
+static void
etoasc (x, string, ndigs)
unsigned EMUSHORT x[];
char *string;
}
else
{
- *s++ = (char )digit + '0';
+ *s++ = (char)digit + '0';
*s++ = '.';
}
/* Generate digits after the decimal point. */
}
+/* Convert ASCII string to quadruple precision floating point
-
-/*
-; 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);
-*/
+ Numeric input is free field decimal number with max of 15 digits with or
+ without decimal point entered as ASCII from teletype. Entering E after
+ the number followed by a second number causes the second number to be
+ interpreted as a power of 10 to be multiplied by the first number
+ (i.e., "scientific" notation). */
/* ASCII to single */
-void
+
+static void
asctoe24 (s, y)
char *s;
unsigned EMUSHORT *y;
/* ASCII to double */
-void
+
+static void
asctoe53 (s, y)
char *s;
unsigned EMUSHORT *y;
/* ASCII to long double */
-void
+
+static void
asctoe64 (s, y)
char *s;
unsigned EMUSHORT *y;
}
/* ASCII to 128-bit long double */
-void
+
+static void
asctoe113 (s, y)
char *s;
unsigned EMUSHORT *y;
}
/* ASCII to super double */
-void
+
+static void
asctoe (s, y)
char *s;
unsigned EMUSHORT *y;
asctoeg (s, y, NBITS);
}
-/* Space to make a copy of the input string: */
-static char lstr[82];
-void
+/* ASCII to e type, with specified rounding precision = oprec. */
+
+static void
asctoeg (ss, y, oprec)
char *ss;
unsigned EMUSHORT *y;
int k, trail, c, rndsav;
EMULONG lexp;
unsigned EMUSHORT nsign, *p;
- char *sp, *s;
+ char *sp, *s, *lstr;
/* Copy the input string. */
+ 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;
if (*s == 'z')
goto donchr;
}
+
/* If enough digits were given to more than fill up the yy register,
- * continuing until overflow into the high guard word yy[2]
- * guarantees that there will be a roundoff bit at the top
- * of the low guard word after normalization.
- */
+ continuing until overflow into the high guard word yy[2]
+ guarantees that there will be a roundoff bit at the top
+ of the low guard word after normalization. */
+
if (yy[2] == 0)
{
if (decflg)
}
else
{
+ /* Mark any lost non-zero digit. */
lost |= k;
+ /* Count lost digits before the decimal point. */
+ if (decflg == 0)
+ nexp -= 1;
}
prec += 1;
goto donchr;
}
lexp = (EXONE - 1 + NBITS) - k;
emdnorm (yy, lost, 0, lexp, 64);
- /* convert to external format */
+ /* Convert to external format:
+
+ Multiply by 10**nexp. If precision is 64 bits,
+ the maximum relative error incurred in forming 10**n
+ for 0 <= n <= 324 is 8.2e-20, at 10**180.
+ For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
+ For 0 >= n >= -999, it is -1.55e-19 at 10**-435. */
- /* Multiply by 10**nexp. If precision is 64 bits,
- * the maximum relative error incurred in forming 10**n
- * for 0 <= n <= 324 is 8.2e-20, at 10**180.
- * For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
- * For 0 >= n >= -999, it is -1.55e-19 at 10**-435.
- */
lexp = yy[E];
if (nexp == 0)
{
nexp = -nexp;
esign = -1;
if (nexp > 4096)
- { /* Punt. Can't handle this without 2 divides. */
+ {
+ /* Punt. Can't handle this without 2 divides. */
emovi (etens[0], tt);
lexp -= tt[E];
k = edivm (tt, yy);
-/* y = largest integer not greater than x
- * (truncated toward minus infinity)
- *
- * unsigned EMUSHORT x[NE], y[NE]
- *
- * efloor (x, y);
- */
+/* y = largest integer not greater than x (truncated toward minus infinity) */
+
static unsigned EMUSHORT bmask[] =
{
0xffff,
0x0000,
};
-void
+static void
efloor (x, y)
unsigned EMUSHORT x[], y[];
{
}
-/* unsigned EMUSHORT x[], s[];
- * int *exp;
- *
- * efrexp (x, exp, s);
- *
- * Returns s and exp such that s * 2**exp = x and .5 <= s < 1.
- * For example, 1.1 = 0.55 * 2**1
- * Handles denormalized numbers properly using long integer exp.
- */
-void
+/* Returns s and exp such that s * 2**exp = x and .5 <= s < 1.
+ For example, 1.1 = 0.55 * 2**1
+ Handles denormalized numbers properly using long integer exp. */
+
+static void
efrexp (x, exp, s)
unsigned EMUSHORT x[];
int *exp;
-/* unsigned EMUSHORT x[], y[];
- * long pwr2;
- *
- * eldexp (x, pwr2, y);
- *
- * Returns y = x * 2**pwr2.
- */
-void
+/* Return y = x * 2**pwr2. */
+
+static void
eldexp (x, pwr2, y)
unsigned EMUSHORT x[];
int pwr2;
/* c = remainder after dividing b by a
- * Least significant integer quotient bits left in equot[].
- */
-void
+ 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))
+ if (eisinf (b)
+ || (ecmp (a, ezero) == 0)
+ || eisnan (a)
+ || eisnan (b))
{
- enan (c);
+ enan (c, 0);
return;
}
#endif
emovo (num, c);
}
-void
+static void
eiremain (den, num)
unsigned EMUSHORT 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
- * INVALID 7 NaN - producing operation
- * EDOM 33 Unix domain error code
- * ERANGE 34 Unix range error code
- *
- * The default version of the file prints the function name,
- * passed to it by the pointer fctnam, followed by the
- * error condition. The display is directed to the standard
- * output device. The routine then returns to the calling
- * program. Users may wish to modify the program to abort by
- * calling exit under severe error conditions such as domain
- * errors.
- *
- * Since all error conditions pass control to this function,
- * the display may be easily changed, eliminated, or directed
- * to an error logging device.
- *
- * SEE ALSO:
- *
- * mconf.h
- *
- */
-\f
-/*
-Cephes Math Library Release 2.0: April, 1987
-Copyright 1984, 1987 by Stephen L. Moshier
-Direct inquiries to 30 Frost Street, Cambridge, MA 02140
-*/
-
-/* include "mconf.h" */
+/* This routine may be called to report one of the following
+ error conditions (in the include file mconf.h).
+
+ Mnemonic Value Significance
+
+ DOMAIN 1 argument domain error
+ SING 2 function singularity
+ OVERFLOW 3 overflow range error
+ UNDERFLOW 4 underflow range error
+ TLOSS 5 total loss of precision
+ PLOSS 6 partial loss of precision
+ INVALID 7 NaN - producing operation
+ EDOM 33 Unix domain error code
+ ERANGE 34 Unix range error code
+
+ The default version of the file prints the function name,
+ passed to it by the pointer fctnam, followed by the
+ error condition. The display is directed to the standard
+ output device. The routine then returns to the calling
+ program. Users may wish to modify the program to abort by
+ calling exit under severe error conditions such as domain
+ errors.
+
+ Since all error conditions pass control to this function,
+ the display may be easily changed, eliminated, or directed
+ to an error logging device. */
+
+/* Note: the order of appearance of the following messages is bound to the
+ error codes defined above. */
-/* Notice: the order of appearance of the following
- * messages is bound to the error codes defined
- * in mconf.h.
- */
#define NMSGS 8
static char *ermsg[NMSGS] =
{
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.
- */
+ /* Display string passed by calling program, which is supposed to be the
+ name of the function in which the error occurred.
+
+ Display error message defined by the code argument. */
- /* Display error message defined
- * by the code argument.
- */
if ((code <= 0) || (code >= NMSGS))
code = 0;
sprintf (errstr, " %s %s error", name, ermsg[code]);
warning (errstr);
/* Set global error message word */
merror = code + 1;
-
- /* Return to calling
- * program
- */
}
#ifdef DEC
-/* Here is etodec.c .
- *
- */
+/* Convert DEC double precision to e type. */
-/*
-; 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;
; etodec (e, &d);
*/
-void
+static void
etodec (x, d)
unsigned EMUSHORT *x, *d;
{
todec (xi, d);
}
-void
+static void
todec (x, y)
unsigned EMUSHORT *x, *y;
{
#endif /* DEC */
#ifdef IBM
-/* Here is etoibm
- *
- */
+/* Convert IBM single/double precision to e type. */
-/*
-; convert IBM single/double precision to e type
-; single/double d;
-; EMUSHORT e[NE];
-; enum machine_mode mode; SFmode/DFmode
-; ibmtoe (&d, e, mode);
-*/
-void
+static void
ibmtoe (d, e, mode)
unsigned EMUSHORT *d;
unsigned EMUSHORT *e;
-/*
-; convert e type to IBM single/double precision
-; single/double d;
-; EMUSHORT e[NE];
-; enum machine_mode mode; SFmode/DFmode
-; etoibm (e, &d, mode);
-*/
+/* Convert e type to IBM single/double precision. */
-void
+static void
etoibm (x, d, mode)
unsigned EMUSHORT *x, *d;
enum machine_mode mode;
toibm (xi, d, mode);
}
-void
+static void
toibm (x, y, mode)
unsigned EMUSHORT *x, *y;
enum machine_mode mode;
#endif
-void
-make_nan (nan, mode)
-unsigned EMUSHORT *nan;
-enum machine_mode mode;
+static void
+make_nan (nan, sign, mode)
+ unsigned EMUSHORT *nan;
+ int sign;
+ enum machine_mode mode;
{
- int i, n;
+ int n;
unsigned EMUSHORT *p;
- n = 0;
switch (mode)
{
/* Possibly the `reserved operand' patterns on a VAX can be
default:
abort ();
}
- for (i=0; i < n; i++)
+#ifdef MIEEE
+ *nan++ = (sign << 15) | *p++;
+#endif
+ while (--n != 0)
*nan++ = *p++;
+#ifndef MIEEE
+ *nan = (sign << 15) | *p;
+#endif
}
/* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
REAL_VALUE_TYPE
ereal_from_float (f)
- unsigned long f;
+ HOST_WIDE_INT f;
{
REAL_VALUE_TYPE r;
unsigned EMUSHORT s[2];
/* 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 WORDS_BIG_ENDIAN
+#if FLOAT_WORDS_BIG_ENDIAN
s[0] = (unsigned EMUSHORT) (f >> 16);
s[1] = (unsigned EMUSHORT) f;
#else
This is the inverse of the function `etardouble' invoked by
REAL_VALUE_TO_TARGET_DOUBLE.
- The DFmode is stored as an array of longs (i.e., HOST_WIDE_INTs)
- with 32 bits of the value per each long. The first element
+ 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)
- unsigned long d[];
+ HOST_WIDE_INT d[];
{
REAL_VALUE_TYPE r;
unsigned EMUSHORT s[4];
unsigned EMUSHORT e[NE];
- /* Convert array of 32 bit pieces to equivalent array of 16 bit pieces.
- This is the inverse of `endian'. */
-#if WORDS_BIG_ENDIAN
+ /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
+#if FLOAT_WORDS_BIG_ENDIAN
s[0] = (unsigned EMUSHORT) (d[0] >> 16);
s[1] = (unsigned EMUSHORT) d[0];
- s[2] = (unsigned EMUSHORT) (d[1] >> 16);
- s[3] = (unsigned EMUSHORT) d[1];
+ if (HOST_BITS_PER_WIDE_INT >= 64)
+ {
+ /* 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);
+ }
+ else
+ {
+ s[2] = (unsigned EMUSHORT) (d[1] >> 16);
+ s[3] = (unsigned EMUSHORT) d[1];
+ }
#else
+/* Target float words are little-endian. */
s[0] = (unsigned EMUSHORT) d[0];
s[1] = (unsigned EMUSHORT) (d[0] >> 16);
- s[2] = (unsigned EMUSHORT) d[1];
- s[3] = (unsigned EMUSHORT) (d[1] >> 16);
+ if (HOST_BITS_PER_WIDE_INT >= 64)
+ {
+ s[2] = (unsigned EMUSHORT) (d[0] >> 32);
+ s[3] = (unsigned EMUSHORT) (d[0] >> 48);
+ }
+ else
+ {
+ s[2] = (unsigned EMUSHORT) d[1];
+ s[3] = (unsigned EMUSHORT) (d[1] >> 16);
+ }
#endif
/* Convert target double to E-type. */
e53toe (s, e);
}
-/* Convert target computer unsigned 64-bit integer to e-type. */
+/* Convert target computer unsigned 64-bit integer to e-type.
+ The endian-ness of DImode follows the convention for integers,
+ so we use WORDS_BIG_ENDIAN here, not FLOAT_WORDS_BIG_ENDIAN. */
-void
+static void
uditoe (di, e)
unsigned EMUSHORT *di; /* Address of the 64-bit int. */
unsigned EMUSHORT *e;
/* Convert target computer signed 64-bit integer to e-type. */
-void
+static void
ditoe (di, e)
unsigned EMUSHORT *di; /* Address of the 64-bit int. */
unsigned EMUSHORT *e;
/* Convert e-type to unsigned 64-bit int. */
-void
+static void
etoudi (x, i)
unsigned EMUSHORT *x;
unsigned EMUSHORT *i;
/* Convert e-type to signed 64-bit int. */
-void
+static void
etodi (x, i)
unsigned EMUSHORT *x;
unsigned EMUSHORT *i;
static int esqinited = 0;
static unsigned short sqrndbit[NI];
-void
+static void
esqrt (x, y)
unsigned EMUSHORT *x, *y;
{
i = ecmp (x, ezero);
if (i <= 0)
{
-#ifdef NANS
- if (i == -2)
+ if (i == -1)
{
- enan (y);
- return;
+ mtherr ("esqrt", DOMAIN);
+ eclear (y);
}
-#endif
- eclear (y);
- if (i < 0)
- mtherr ("esqrt", DOMAIN);
+ else
+ emov (x, y);
return;
}
projects / platform / upstream / gcc.git / blobdiff