1 /* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
2 and support for XFmode IEEE extended real floating point arithmetic.
3 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998,
4 1999, 2000, 2002 Free Software Foundation, Inc.
5 Contributed by Stephen L. Moshier (moshier@world.std.com).
7 This file is part of GCC.
9 GCC is free software; you can redistribute it and/or modify it under
10 the terms of the GNU General Public License as published by the Free
11 Software Foundation; either version 2, or (at your option) any later
14 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
15 WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19 You should have received a copy of the GNU General Public License
20 along with GCC; see the file COPYING. If not, write to the Free
21 Software Foundation, 59 Temple Place - Suite 330, Boston, MA
30 /* To enable support of XFmode extended real floating point, define
31 LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
33 Machine files (tm.h etc) must not contain any code
34 that tries to use host floating point arithmetic to convert
35 REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
36 etc. In cross-compile situations a REAL_VALUE_TYPE may not
37 be intelligible to the host computer's native arithmetic.
39 The first part of this file interfaces gcc to a floating point
40 arithmetic suite that was not written with gcc in mind. Avoid
41 changing the low-level arithmetic routines unless you have suitable
42 test programs available. A special version of the PARANOIA floating
43 point arithmetic tester, modified for this purpose, can be found on
44 usc.edu: /pub/C-numanal/ieeetest.zoo. Other tests, and libraries of
45 XFmode and TFmode transcendental functions, can be obtained by ftp from
46 netlib.att.com: netlib/cephes. */
48 /* Type of computer arithmetic.
49 Only one of DEC, IBM, IEEE, C4X, or UNK should get defined.
51 `IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
52 to big-endian IEEE floating-point data structure. This definition
53 should work in SFmode `float' type and DFmode `double' type on
54 virtually all big-endian IEEE machines. If LONG_DOUBLE_TYPE_SIZE
55 has been defined to be 96, then IEEE also invokes the particular
56 XFmode (`long double' type) data structure used by the Motorola
57 680x0 series processors.
59 `IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
60 little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
61 has been defined to be 96, then IEEE also invokes the particular
62 XFmode `long double' data structure used by the Intel 80x86 series
65 `DEC' refers specifically to the Digital Equipment Corp PDP-11
66 and VAX floating point data structure. This model currently
67 supports no type wider than DFmode.
69 `IBM' refers specifically to the IBM System/370 and compatible
70 floating point data structure. This model currently supports
71 no type wider than DFmode. The IBM conversions were contributed by
72 frank@atom.ansto.gov.au (Frank Crawford).
74 `C4X' refers specifically to the floating point format used on
75 Texas Instruments TMS320C3x and TMS320C4x digital signal
76 processors. This supports QFmode (32-bit float, double) and HFmode
77 (40-bit long double) where BITS_PER_BYTE is 32. Unlike IEEE
78 floats, C4x floats are not rounded to be even. The C4x conversions
79 were contributed by m.hayes@elec.canterbury.ac.nz (Michael Hayes) and
80 Haj.Ten.Brugge@net.HCC.nl (Herman ten Brugge).
82 If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
83 then `long double' and `double' are both implemented, but they
86 The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
87 and may deactivate XFmode since `long double' is used to refer
88 to both modes. Defining INTEL_EXTENDED_IEEE_FORMAT to non-zero
89 at the same time enables 80387-style 80-bit floats in a 128-bit
90 padded image, as seen on IA-64.
92 The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
93 contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
94 separate the floating point unit's endian-ness from that of
95 the integer addressing. This permits one to define a big-endian
96 FPU on a little-endian machine (e.g., ARM). An extension to
97 BYTES_BIG_ENDIAN may be required for some machines in the future.
98 These optional macros may be defined in tm.h. In real.h, they
99 default to WORDS_BIG_ENDIAN, etc., so there is no need to define
100 them for any normal host or target machine on which the floats
101 and the integers have the same endian-ness. */
104 /* The following converts gcc macros into the ones used by this file. */
106 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
107 /* PDP-11, Pro350, VAX: */
109 #else /* it's not VAX */
110 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
111 /* IBM System/370 style */
113 #else /* it's also not an IBM */
114 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
115 /* TMS320C3x/C4x style */
117 #else /* it's also not a C4X */
118 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
120 #else /* it's not IEEE either */
121 /* UNKnown arithmetic. We don't support this and can't go on. */
122 unknown arithmetic type
124 #endif /* not IEEE */
129 #define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
131 /* Define INFINITY for support of infinity.
132 Define NANS for support of Not-a-Number's (NaN's). */
133 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
138 /* Support of NaNs requires support of infinity. */
145 /* Find a host integer type that is at least 16 bits wide,
146 and another type at least twice whatever that size is. */
148 #if HOST_BITS_PER_CHAR >= 16
149 #define EMUSHORT char
150 #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
151 #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
153 #if HOST_BITS_PER_SHORT >= 16
154 #define EMUSHORT short
155 #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
156 #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
158 #if HOST_BITS_PER_INT >= 16
160 #define EMUSHORT_SIZE HOST_BITS_PER_INT
161 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
163 #if HOST_BITS_PER_LONG >= 16
164 #define EMUSHORT long
165 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
166 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
168 #error "You will have to modify this program to have a smaller unit size."
174 /* If no 16-bit type has been found and the compiler is GCC, try HImode. */
175 #if defined(__GNUC__) && EMUSHORT_SIZE != 16
176 typedef int HItype __attribute__ ((mode (HI)));
177 typedef unsigned int UHItype __attribute__ ((mode (HI)));
181 #define EMUSHORT HItype
182 #define UEMUSHORT UHItype
183 #define EMUSHORT_SIZE 16
184 #define EMULONG_SIZE 32
186 #define UEMUSHORT unsigned EMUSHORT
189 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
190 #define EMULONG short
192 #if HOST_BITS_PER_INT >= EMULONG_SIZE
195 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
198 #if HOST_BITS_PER_LONGLONG >= EMULONG_SIZE
199 #define EMULONG long long int
201 #error "You will have to modify this program to have a smaller unit size."
207 #if EMUSHORT_SIZE != 16
208 #error "The host interface doesn't work if no 16-bit size exists."
211 /* Calculate the size of the generic "e" type. This always has
212 identical in-memory size and representation to REAL_VALUE_TYPE.
213 There are only two supported sizes: ten and six 16-bit words (160
216 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && !INTEL_EXTENDED_IEEE_FORMAT
219 # define MAXDECEXP 4932
220 # define MINDECEXP -4977
223 # define MAXDECEXP 4932
224 # define MINDECEXP -4956
227 /* Fail compilation if 2*NE is not the appropriate size. */
229 struct compile_test_dummy {
230 char twice_NE_must_equal_sizeof_REAL_VALUE_TYPE
231 [(sizeof (REAL_VALUE_TYPE) == 2*NE) ? 1 : -1];
234 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
235 In GET_REAL and PUT_REAL, r and e are pointers.
236 A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
237 in memory, with no holes. */
238 #define GET_REAL(r, e) memcpy ((e), (r), 2*NE)
239 #define PUT_REAL(e, r) memcpy ((r), (e), 2*NE)
241 /* Number of 16 bit words in internal format */
244 /* Array offset to exponent */
247 /* Array offset to high guard word */
250 /* Number of bits of precision */
251 #define NBITS ((NI-4)*16)
253 /* Maximum number of decimal digits in ASCII conversion
256 #define NDEC (NBITS*8/27)
258 /* The exponent of 1.0 */
259 #define EXONE (0x3fff)
261 #if defined(HOST_EBCDIC)
262 /* bit 8 is significant in EBCDIC */
263 #define CHARMASK 0xff
265 #define CHARMASK 0x7f
268 extern int extra_warnings;
269 extern const UEMUSHORT ezero[NE], ehalf[NE], eone[NE], etwo[NE];
270 extern const UEMUSHORT elog2[NE], esqrt2[NE];
272 static void endian PARAMS ((const UEMUSHORT *, long *,
274 static void eclear PARAMS ((UEMUSHORT *));
275 static void emov PARAMS ((const UEMUSHORT *, UEMUSHORT *));
277 static void eabs PARAMS ((UEMUSHORT *));
279 static void eneg PARAMS ((UEMUSHORT *));
280 static int eisneg PARAMS ((const UEMUSHORT *));
281 static int eisinf PARAMS ((const UEMUSHORT *));
282 static int eisnan PARAMS ((const UEMUSHORT *));
283 static void einfin PARAMS ((UEMUSHORT *));
285 static void enan PARAMS ((UEMUSHORT *, int));
286 static void einan PARAMS ((UEMUSHORT *));
287 static int eiisnan PARAMS ((const UEMUSHORT *));
288 static void make_nan PARAMS ((UEMUSHORT *, int, enum machine_mode));
290 static int eiisneg PARAMS ((const UEMUSHORT *));
291 static void saturate PARAMS ((UEMUSHORT *, int, int, int));
292 static void emovi PARAMS ((const UEMUSHORT *, UEMUSHORT *));
293 static void emovo PARAMS ((const UEMUSHORT *, UEMUSHORT *));
294 static void ecleaz PARAMS ((UEMUSHORT *));
295 static void ecleazs PARAMS ((UEMUSHORT *));
296 static void emovz PARAMS ((const UEMUSHORT *, UEMUSHORT *));
298 static void eiinfin PARAMS ((UEMUSHORT *));
301 static int eiisinf PARAMS ((const UEMUSHORT *));
303 static int ecmpm PARAMS ((const UEMUSHORT *, const UEMUSHORT *));
304 static void eshdn1 PARAMS ((UEMUSHORT *));
305 static void eshup1 PARAMS ((UEMUSHORT *));
306 static void eshdn8 PARAMS ((UEMUSHORT *));
307 static void eshup8 PARAMS ((UEMUSHORT *));
308 static void eshup6 PARAMS ((UEMUSHORT *));
309 static void eshdn6 PARAMS ((UEMUSHORT *));
310 static void eaddm PARAMS ((const UEMUSHORT *, UEMUSHORT *));
\f
311 static void esubm PARAMS ((const UEMUSHORT *, UEMUSHORT *));
312 static void m16m PARAMS ((unsigned int, const UEMUSHORT *, UEMUSHORT *));
313 static int edivm PARAMS ((const UEMUSHORT *, UEMUSHORT *));
314 static int emulm PARAMS ((const UEMUSHORT *, UEMUSHORT *));
315 static void emdnorm PARAMS ((UEMUSHORT *, int, int, EMULONG, int));
316 static void esub PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
318 static void eadd PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
320 static void eadd1 PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
322 static void ediv PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
324 static void emul PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
326 static void e53toe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
327 static void e64toe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
328 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
329 static void e113toe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
331 static void e24toe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
332 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
333 static void etoe113 PARAMS ((const UEMUSHORT *, UEMUSHORT *));
334 static void toe113 PARAMS ((UEMUSHORT *, UEMUSHORT *));
336 static void etoe64 PARAMS ((const UEMUSHORT *, UEMUSHORT *));
337 static void toe64 PARAMS ((UEMUSHORT *, UEMUSHORT *));
338 static void etoe53 PARAMS ((const UEMUSHORT *, UEMUSHORT *));
339 static void toe53 PARAMS ((UEMUSHORT *, UEMUSHORT *));
340 static void etoe24 PARAMS ((const UEMUSHORT *, UEMUSHORT *));
341 static void toe24 PARAMS ((UEMUSHORT *, UEMUSHORT *));
342 static int ecmp PARAMS ((const UEMUSHORT *, const UEMUSHORT *));
344 static void eround PARAMS ((const UEMUSHORT *, UEMUSHORT *));
346 static void ltoe PARAMS ((const HOST_WIDE_INT *, UEMUSHORT *));
347 static void ultoe PARAMS ((const unsigned HOST_WIDE_INT *, UEMUSHORT *));
348 static void eifrac PARAMS ((const UEMUSHORT *, HOST_WIDE_INT *,
350 static void euifrac PARAMS ((const UEMUSHORT *, unsigned HOST_WIDE_INT *,
352 static int eshift PARAMS ((UEMUSHORT *, int));
353 static int enormlz PARAMS ((UEMUSHORT *));
355 static void e24toasc PARAMS ((const UEMUSHORT *, char *, int));
356 static void e53toasc PARAMS ((const UEMUSHORT *, char *, int));
357 static void e64toasc PARAMS ((const UEMUSHORT *, char *, int));
358 static void e113toasc PARAMS ((const UEMUSHORT *, char *, int));
360 static void etoasc PARAMS ((const UEMUSHORT *, char *, int));
361 static void asctoe24 PARAMS ((const char *, UEMUSHORT *));
362 static void asctoe53 PARAMS ((const char *, UEMUSHORT *));
363 static void asctoe64 PARAMS ((const char *, UEMUSHORT *));
364 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
365 static void asctoe113 PARAMS ((const char *, UEMUSHORT *));
367 static void asctoe PARAMS ((const char *, UEMUSHORT *));
368 static void asctoeg PARAMS ((const char *, UEMUSHORT *, int));
369 static void efloor PARAMS ((const UEMUSHORT *, UEMUSHORT *));
371 static void efrexp PARAMS ((const UEMUSHORT *, int *,
374 static void eldexp PARAMS ((const UEMUSHORT *, int, UEMUSHORT *));
376 static void eremain PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
379 static void eiremain PARAMS ((UEMUSHORT *, UEMUSHORT *));
380 static void mtherr PARAMS ((const char *, int));
382 static void dectoe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
383 static void etodec PARAMS ((const UEMUSHORT *, UEMUSHORT *));
384 static void todec PARAMS ((UEMUSHORT *, UEMUSHORT *));
387 static void ibmtoe PARAMS ((const UEMUSHORT *, UEMUSHORT *,
389 static void etoibm PARAMS ((const UEMUSHORT *, UEMUSHORT *,
391 static void toibm PARAMS ((UEMUSHORT *, UEMUSHORT *,
395 static void c4xtoe PARAMS ((const UEMUSHORT *, UEMUSHORT *,
397 static void etoc4x PARAMS ((const UEMUSHORT *, UEMUSHORT *,
399 static void toc4x PARAMS ((UEMUSHORT *, UEMUSHORT *,
403 static void uditoe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
404 static void ditoe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
405 static void etoudi PARAMS ((const UEMUSHORT *, UEMUSHORT *));
406 static void etodi PARAMS ((const UEMUSHORT *, UEMUSHORT *));
407 static void esqrt PARAMS ((const UEMUSHORT *, UEMUSHORT *));
410 /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
411 swapping ends if required, into output array of longs. The
412 result is normally passed to fprintf by the ASM_OUTPUT_ macros. */
418 enum machine_mode mode;
422 if (REAL_WORDS_BIG_ENDIAN)
427 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
428 /* Swap halfwords in the fourth long. */
429 th = (unsigned long) e[6] & 0xffff;
430 t = (unsigned long) e[7] & 0xffff;
439 /* Swap halfwords in the third long. */
440 th = (unsigned long) e[4] & 0xffff;
441 t = (unsigned long) e[5] & 0xffff;
447 /* Swap halfwords in the second word. */
448 th = (unsigned long) e[2] & 0xffff;
449 t = (unsigned long) e[3] & 0xffff;
456 /* Swap halfwords in the first word. */
457 th = (unsigned long) e[0] & 0xffff;
458 t = (unsigned long) e[1] & 0xffff;
469 /* Pack the output array without swapping. */
474 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
475 /* Pack the fourth long. */
476 th = (unsigned long) e[7] & 0xffff;
477 t = (unsigned long) e[6] & 0xffff;
486 /* Pack the third long.
487 Each element of the input REAL_VALUE_TYPE array has 16 useful bits
489 th = (unsigned long) e[5] & 0xffff;
490 t = (unsigned long) e[4] & 0xffff;
496 /* Pack the second long */
497 th = (unsigned long) e[3] & 0xffff;
498 t = (unsigned long) e[2] & 0xffff;
505 /* Pack the first long */
506 th = (unsigned long) e[1] & 0xffff;
507 t = (unsigned long) e[0] & 0xffff;
519 /* This is the implementation of the REAL_ARITHMETIC macro. */
522 earith (value, icode, r1, r2)
523 REAL_VALUE_TYPE *value;
528 UEMUSHORT d1[NE], d2[NE], v[NE];
534 /* Return NaN input back to the caller. */
537 PUT_REAL (d1, value);
542 PUT_REAL (d2, value);
546 code = (enum tree_code) icode;
554 esub (d2, d1, v); /* d1 - d2 */
563 if (ecmp (d2, ezero) == 0)
566 ediv (d2, d1, v); /* d1/d2 */
569 case MIN_EXPR: /* min (d1,d2) */
570 if (ecmp (d1, d2) < 0)
576 case MAX_EXPR: /* max (d1,d2) */
577 if (ecmp (d1, d2) > 0)
590 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
591 implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */
597 UEMUSHORT f[NE], g[NE];
613 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
614 implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */
620 UEMUSHORT f[NE], g[NE];
622 unsigned HOST_WIDE_INT l;
636 /* This is the REAL_VALUE_ATOF function. It converts a decimal or hexadecimal
637 string to binary, rounding off as indicated by the machine_mode argument.
638 Then it promotes the rounded value to REAL_VALUE_TYPE. */
645 UEMUSHORT tem[NE], e[NE];
671 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
691 /* Expansion of REAL_NEGATE. */
707 /* Round real toward zero to HOST_WIDE_INT;
708 implements REAL_VALUE_FIX (x). */
714 UEMUSHORT f[NE], g[NE];
721 warning ("conversion from NaN to int");
729 /* Round real toward zero to unsigned HOST_WIDE_INT
730 implements REAL_VALUE_UNSIGNED_FIX (x).
731 Negative input returns zero. */
733 unsigned HOST_WIDE_INT
737 UEMUSHORT f[NE], g[NE];
738 unsigned HOST_WIDE_INT l;
744 warning ("conversion from NaN to unsigned int");
753 /* REAL_VALUE_FROM_INT macro. */
756 ereal_from_int (d, i, j, mode)
759 enum machine_mode mode;
761 UEMUSHORT df[NE], dg[NE];
762 HOST_WIDE_INT low, high;
765 if (GET_MODE_CLASS (mode) != MODE_FLOAT)
772 /* complement and add 1 */
779 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
780 ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
782 ultoe ((unsigned HOST_WIDE_INT *) &low, df);
787 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
788 Avoid double-rounding errors later by rounding off now from the
789 extra-wide internal format to the requested precision. */
790 switch (GET_MODE_BITSIZE (mode))
808 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
825 /* REAL_VALUE_FROM_UNSIGNED_INT macro. */
828 ereal_from_uint (d, i, j, mode)
830 unsigned HOST_WIDE_INT i, j;
831 enum machine_mode mode;
833 UEMUSHORT df[NE], dg[NE];
834 unsigned HOST_WIDE_INT low, high;
836 if (GET_MODE_CLASS (mode) != MODE_FLOAT)
840 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
846 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
847 Avoid double-rounding errors later by rounding off now from the
848 extra-wide internal format to the requested precision. */
849 switch (GET_MODE_BITSIZE (mode))
867 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
884 /* REAL_VALUE_TO_INT macro. */
887 ereal_to_int (low, high, rr)
888 HOST_WIDE_INT *low, *high;
891 UEMUSHORT d[NE], df[NE], dg[NE], dh[NE];
898 warning ("conversion from NaN to int");
904 /* convert positive value */
911 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
912 ediv (df, d, dg); /* dg = d / 2^32 is the high word */
913 euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
914 emul (df, dh, dg); /* fractional part is the low word */
915 euifrac (dg, (unsigned HOST_WIDE_INT *) low, dh);
918 /* complement and add 1 */
928 /* REAL_VALUE_LDEXP macro. */
935 UEMUSHORT e[NE], y[NE];
948 /* Check for infinity in a REAL_VALUE_TYPE. */
952 REAL_VALUE_TYPE x ATTRIBUTE_UNUSED;
964 /* Check whether a REAL_VALUE_TYPE item is a NaN. */
968 REAL_VALUE_TYPE x ATTRIBUTE_UNUSED;
981 /* Check for a negative REAL_VALUE_TYPE number.
982 This just checks the sign bit, so that -0 counts as negative. */
988 return ereal_isneg (x);
991 /* Expansion of REAL_VALUE_TRUNCATE.
992 The result is in floating point, rounded to nearest or even. */
995 real_value_truncate (mode, arg)
996 enum machine_mode mode;
999 UEMUSHORT e[NE], t[NE];
1011 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
1048 /* If an unsupported type was requested, presume that
1049 the machine files know something useful to do with
1050 the unmodified value. */
1059 /* Try to change R into its exact multiplicative inverse in machine mode
1060 MODE. Return nonzero function value if successful. */
1063 exact_real_inverse (mode, r)
1064 enum machine_mode mode;
1067 UEMUSHORT e[NE], einv[NE];
1068 REAL_VALUE_TYPE rinv;
1073 /* Test for input in range. Don't transform IEEE special values. */
1074 if (eisinf (e) || eisnan (e) || (ecmp (e, ezero) == 0))
1077 /* Test for a power of 2: all significand bits zero except the MSB.
1078 We are assuming the target has binary (or hex) arithmetic. */
1079 if (e[NE - 2] != 0x8000)
1082 for (i = 0; i < NE - 2; i++)
1088 /* Compute the inverse and truncate it to the required mode. */
1089 ediv (e, eone, einv);
1090 PUT_REAL (einv, &rinv);
1091 rinv = real_value_truncate (mode, rinv);
1093 #ifdef CHECK_FLOAT_VALUE
1094 /* This check is not redundant. It may, for example, flush
1095 a supposedly IEEE denormal value to zero. */
1097 if (CHECK_FLOAT_VALUE (mode, rinv, i))
1100 GET_REAL (&rinv, einv);
1102 /* Check the bits again, because the truncation might have
1103 generated an arbitrary saturation value on overflow. */
1104 if (einv[NE - 2] != 0x8000)
1107 for (i = 0; i < NE - 2; i++)
1113 /* Fail if the computed inverse is out of range. */
1114 if (eisinf (einv) || eisnan (einv) || (ecmp (einv, ezero) == 0))
1117 /* Output the reciprocal and return success flag. */
1122 /* Used for debugging--print the value of R in human-readable format
1131 REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
1132 fprintf (stderr, "%s", dstr);
1136 /* The following routines convert REAL_VALUE_TYPE to the various floating
1137 point formats that are meaningful to supported computers.
1139 The results are returned in 32-bit pieces, each piece stored in a `long'.
1140 This is so they can be printed by statements like
1142 fprintf (file, "%lx, %lx", L[0], L[1]);
1144 that will work on both narrow- and wide-word host computers. */
1146 /* Convert R to a 128-bit long double precision value. The output array L
1147 contains four 32-bit pieces of the result, in the order they would appear
1158 #if INTEL_EXTENDED_IEEE_FORMAT == 0
1163 endian (e, l, TFmode);
1166 /* Convert R to a double extended precision value. The output array L
1167 contains three 32-bit pieces of the result, in the order they would
1168 appear in memory. */
1179 endian (e, l, XFmode);
1182 /* Convert R to a double precision value. The output array L contains two
1183 32-bit pieces of the result, in the order they would appear in memory. */
1194 endian (e, l, DFmode);
1197 /* Convert R to a single precision float value stored in the least-significant
1198 bits of a `long'. */
1209 endian (e, &l, SFmode);
1213 /* Convert X to a decimal ASCII string S for output to an assembly
1214 language file. Note, there is no standard way to spell infinity or
1215 a NaN, so these values may require special treatment in the tm.h
1219 ereal_to_decimal (x, s)
1229 /* Compare X and Y. Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
1230 or -2 if either is a NaN. */
1234 REAL_VALUE_TYPE x, y;
1236 UEMUSHORT ex[NE], ey[NE];
1240 return (ecmp (ex, ey));
1243 /* Return 1 if the sign bit of X is set, else return 0. */
1252 return (eisneg (ex));
1257 Extended precision IEEE binary floating point arithmetic routines
1259 Numbers are stored in C language as arrays of 16-bit unsigned
1260 short integers. The arguments of the routines are pointers to
1263 External e type data structure, similar to Intel 8087 chip
1264 temporary real format but possibly with a larger significand:
1266 NE-1 significand words (least significant word first,
1267 most significant bit is normally set)
1268 exponent (value = EXONE for 1.0,
1269 top bit is the sign)
1272 Internal exploded e-type data structure of a number (a "word" is 16 bits):
1274 ei[0] sign word (0 for positive, 0xffff for negative)
1275 ei[1] biased exponent (value = EXONE for the number 1.0)
1276 ei[2] high guard word (always zero after normalization)
1278 to ei[NI-2] significand (NI-4 significand words,
1279 most significant word first,
1280 most significant bit is set)
1281 ei[NI-1] low guard word (0x8000 bit is rounding place)
1285 Routines for external format e-type numbers
1287 asctoe (string, e) ASCII string to extended double e type
1288 asctoe64 (string, &d) ASCII string to long double
1289 asctoe53 (string, &d) ASCII string to double
1290 asctoe24 (string, &f) ASCII string to single
1291 asctoeg (string, e, prec) ASCII string to specified precision
1292 e24toe (&f, e) IEEE single precision to e type
1293 e53toe (&d, e) IEEE double precision to e type
1294 e64toe (&d, e) IEEE long double precision to e type
1295 e113toe (&d, e) 128-bit long double precision to e type
1297 eabs (e) absolute value
1299 eadd (a, b, c) c = b + a
1301 ecmp (a, b) Returns 1 if a > b, 0 if a == b,
1302 -1 if a < b, -2 if either a or b is a NaN.
1303 ediv (a, b, c) c = b / a
1304 efloor (a, b) truncate to integer, toward -infinity
1305 efrexp (a, exp, s) extract exponent and significand
1306 eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
1307 euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
1308 einfin (e) set e to infinity, leaving its sign alone
1309 eldexp (a, n, b) multiply by 2**n
1311 emul (a, b, c) c = b * a
1314 eround (a, b) b = nearest integer value to a
1316 esub (a, b, c) c = b - a
1318 e24toasc (&f, str, n) single to ASCII string, n digits after decimal
1319 e53toasc (&d, str, n) double to ASCII string, n digits after decimal
1320 e64toasc (&d, str, n) 80-bit long double to ASCII string
1321 e113toasc (&d, str, n) 128-bit long double to ASCII string
1323 etoasc (e, str, n) e to ASCII string, n digits after decimal
1324 etoe24 (e, &f) convert e type to IEEE single precision
1325 etoe53 (e, &d) convert e type to IEEE double precision
1326 etoe64 (e, &d) convert e type to IEEE long double precision
1327 ltoe (&l, e) HOST_WIDE_INT to e type
1328 ultoe (&l, e) unsigned HOST_WIDE_INT to e type
1329 eisneg (e) 1 if sign bit of e != 0, else 0
1330 eisinf (e) 1 if e has maximum exponent (non-IEEE)
1331 or is infinite (IEEE)
1332 eisnan (e) 1 if e is a NaN
1335 Routines for internal format exploded e-type numbers
1337 eaddm (ai, bi) add significands, bi = bi + ai
1339 ecleazs (ei) set ei = 0 but leave its sign alone
1340 ecmpm (ai, bi) compare significands, return 1, 0, or -1
1341 edivm (ai, bi) divide significands, bi = bi / ai
1342 emdnorm (ai,l,s,exp) normalize and round off
1343 emovi (a, ai) convert external a to internal ai
1344 emovo (ai, a) convert internal ai to external a
1345 emovz (ai, bi) bi = ai, low guard word of bi = 0
1346 emulm (ai, bi) multiply significands, bi = bi * ai
1347 enormlz (ei) left-justify the significand
1348 eshdn1 (ai) shift significand and guards down 1 bit
1349 eshdn8 (ai) shift down 8 bits
1350 eshdn6 (ai) shift down 16 bits
1351 eshift (ai, n) shift ai n bits up (or down if n < 0)
1352 eshup1 (ai) shift significand and guards up 1 bit
1353 eshup8 (ai) shift up 8 bits
1354 eshup6 (ai) shift up 16 bits
1355 esubm (ai, bi) subtract significands, bi = bi - ai
1356 eiisinf (ai) 1 if infinite
1357 eiisnan (ai) 1 if a NaN
1358 eiisneg (ai) 1 if sign bit of ai != 0, else 0
1359 einan (ai) set ai = NaN
1361 eiinfin (ai) set ai = infinity
1364 The result is always normalized and rounded to NI-4 word precision
1365 after each arithmetic operation.
1367 Exception flags are NOT fully supported.
1369 Signaling NaN's are NOT supported; they are treated the same
1372 Define INFINITY for support of infinity; otherwise a
1373 saturation arithmetic is implemented.
1375 Define NANS for support of Not-a-Number items; otherwise the
1376 arithmetic will never produce a NaN output, and might be confused
1378 If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1379 either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1380 may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1383 Denormals are always supported here where appropriate (e.g., not
1384 for conversion to DEC numbers). */
1386 /* Definitions for error codes that are passed to the common error handling
1389 For Digital Equipment PDP-11 and VAX computers, certain
1390 IBM systems, and others that use numbers with a 56-bit
1391 significand, the symbol DEC should be defined. In this
1392 mode, most floating point constants are given as arrays
1393 of octal integers to eliminate decimal to binary conversion
1394 errors that might be introduced by the compiler.
1396 For computers, such as IBM PC, that follow the IEEE
1397 Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1398 Std 754-1985), the symbol IEEE should be defined.
1399 These numbers have 53-bit significands. In this mode, constants
1400 are provided as arrays of hexadecimal 16 bit integers.
1401 The endian-ness of generated values is controlled by
1402 REAL_WORDS_BIG_ENDIAN.
1404 To accommodate other types of computer arithmetic, all
1405 constants are also provided in a normal decimal radix
1406 which one can hope are correctly converted to a suitable
1407 format by the available C language compiler. To invoke
1408 this mode, the symbol UNK is defined.
1410 An important difference among these modes is a predefined
1411 set of machine arithmetic constants for each. The numbers
1412 MACHEP (the machine roundoff error), MAXNUM (largest number
1413 represented), and several other parameters are preset by
1414 the configuration symbol. Check the file const.c to
1415 ensure that these values are correct for your computer.
1417 For ANSI C compatibility, define ANSIC equal to 1. Currently
1418 this affects only the atan2 function and others that use it. */
1420 /* Constant definitions for math error conditions. */
1422 #define DOMAIN 1 /* argument domain error */
1423 #define SING 2 /* argument singularity */
1424 #define OVERFLOW 3 /* overflow range error */
1425 #define UNDERFLOW 4 /* underflow range error */
1426 #define TLOSS 5 /* total loss of precision */
1427 #define PLOSS 6 /* partial loss of precision */
1428 #define INVALID 7 /* NaN-producing operation */
1430 /* e type constants used by high precision check routines */
1432 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
1434 const UEMUSHORT ezero[NE] =
1435 {0x0000, 0x0000, 0x0000, 0x0000,
1436 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
1439 const UEMUSHORT ehalf[NE] =
1440 {0x0000, 0x0000, 0x0000, 0x0000,
1441 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
1444 const UEMUSHORT eone[NE] =
1445 {0x0000, 0x0000, 0x0000, 0x0000,
1446 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
1449 const UEMUSHORT etwo[NE] =
1450 {0x0000, 0x0000, 0x0000, 0x0000,
1451 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
1454 const UEMUSHORT e32[NE] =
1455 {0x0000, 0x0000, 0x0000, 0x0000,
1456 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
1458 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1459 const UEMUSHORT elog2[NE] =
1460 {0x40f3, 0xf6af, 0x03f2, 0xb398,
1461 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1463 /* 1.41421356237309504880168872420969807856967187537695E0 */
1464 const UEMUSHORT esqrt2[NE] =
1465 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1466 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1468 /* 3.14159265358979323846264338327950288419716939937511E0 */
1469 const UEMUSHORT epi[NE] =
1470 {0x2902, 0x1cd1, 0x80dc, 0x628b,
1471 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1474 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1475 const UEMUSHORT ezero[NE] =
1476 {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1477 const UEMUSHORT ehalf[NE] =
1478 {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1479 const UEMUSHORT eone[NE] =
1480 {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1481 const UEMUSHORT etwo[NE] =
1482 {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1483 const UEMUSHORT e32[NE] =
1484 {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1485 const UEMUSHORT elog2[NE] =
1486 {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1487 const UEMUSHORT esqrt2[NE] =
1488 {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1489 const UEMUSHORT epi[NE] =
1490 {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1493 /* Control register for rounding precision.
1494 This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits. */
1499 /* Clear out entire e-type number X. */
1507 for (i = 0; i < NE; i++)
1511 /* Move e-type number from A to B. */
1520 for (i = 0; i < NE; i++)
1526 /* Absolute value of e-type X. */
1532 /* sign is top bit of last word of external format */
1533 x[NE - 1] &= 0x7fff;
1537 /* Negate the e-type number X. */
1544 x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
1547 /* Return 1 if sign bit of e-type number X is nonzero, else zero. */
1551 const UEMUSHORT x[];
1554 if (x[NE - 1] & 0x8000)
1560 /* Return 1 if e-type number X is infinity, else return zero. */
1564 const UEMUSHORT x[];
1571 if ((x[NE - 1] & 0x7fff) == 0x7fff)
1577 /* Check if e-type number is not a number. The bit pattern is one that we
1578 defined, so we know for sure how to detect it. */
1582 const UEMUSHORT x[] ATTRIBUTE_UNUSED;
1587 /* NaN has maximum exponent */
1588 if ((x[NE - 1] & 0x7fff) != 0x7fff)
1590 /* ... and non-zero significand field. */
1591 for (i = 0; i < NE - 1; i++)
1601 /* Fill e-type number X with infinity pattern (IEEE)
1602 or largest possible number (non-IEEE). */
1611 for (i = 0; i < NE - 1; i++)
1615 for (i = 0; i < NE - 1; i++)
1643 /* Output an e-type NaN.
1644 This generates Intel's quiet NaN pattern for extended real.
1645 The exponent is 7fff, the leading mantissa word is c000. */
1655 for (i = 0; i < NE - 2; i++)
1658 *x = (sign << 15) | 0x7fff;
1662 /* Move in an e-type number A, converting it to exploded e-type B. */
1674 p = a + (NE - 1); /* point to last word of external number */
1675 /* get the sign bit */
1680 /* get the exponent */
1682 *q++ &= 0x7fff; /* delete the sign bit */
1684 if ((*(q - 1) & 0x7fff) == 0x7fff)
1690 for (i = 3; i < NI; i++)
1696 for (i = 2; i < NI; i++)
1702 /* clear high guard word */
1704 /* move in the significand */
1705 for (i = 0; i < NE - 1; i++)
1707 /* clear low guard word */
1711 /* Move out exploded e-type number A, converting it to e type B. */
1724 q = b + (NE - 1); /* point to output exponent */
1725 /* combine sign and exponent */
1728 *q-- = *p++ | 0x8000;
1732 if (*(p - 1) == 0x7fff)
1737 enan (b, eiisneg (a));
1745 /* skip over guard word */
1747 /* move the significand */
1748 for (j = 0; j < NE - 1; j++)
1752 /* Clear out exploded e-type number XI. */
1760 for (i = 0; i < NI; i++)
1764 /* Clear out exploded e-type XI, but don't touch the sign. */
1773 for (i = 0; i < NI - 1; i++)
1777 /* Move exploded e-type number from A to B. */
1786 for (i = 0; i < NI - 1; i++)
1788 /* clear low guard word */
1792 /* Generate exploded e-type NaN.
1793 The explicit pattern for this is maximum exponent and
1794 top two significant bits set. */
1808 /* Return nonzero if exploded e-type X is a NaN. */
1813 const UEMUSHORT x[];
1817 if ((x[E] & 0x7fff) == 0x7fff)
1819 for (i = M + 1; i < NI; i++)
1829 /* Return nonzero if sign of exploded e-type X is nonzero. */
1833 const UEMUSHORT x[];
1840 /* Fill exploded e-type X with infinity pattern.
1841 This has maximum exponent and significand all zeros. */
1853 /* Return nonzero if exploded e-type X is infinite. */
1858 const UEMUSHORT x[];
1865 if ((x[E] & 0x7fff) == 0x7fff)
1869 #endif /* INFINITY */
1871 /* Compare significands of numbers in internal exploded e-type format.
1872 Guard words are included in the comparison.
1880 const UEMUSHORT *a, *b;
1884 a += M; /* skip up to significand area */
1886 for (i = M; i < NI; i++)
1894 if (*(--a) > *(--b))
1900 /* Shift significand of exploded e-type X down by 1 bit. */
1909 x += M; /* point to significand area */
1912 for (i = M; i < NI; i++)
1924 /* Shift significand of exploded e-type X up by 1 bit. */
1936 for (i = M; i < NI; i++)
1949 /* Shift significand of exploded e-type X down by 8 bits. */
1955 UEMUSHORT newbyt, oldbyt;
1960 for (i = M; i < NI; i++)
1970 /* Shift significand of exploded e-type X up by 8 bits. */
1977 UEMUSHORT newbyt, oldbyt;
1982 for (i = M; i < NI; i++)
1992 /* Shift significand of exploded e-type X up by 16 bits. */
2004 for (i = M; i < NI - 1; i++)
2010 /* Shift significand of exploded e-type X down by 16 bits. */
2022 for (i = M; i < NI - 1; i++)
2028 /* Add significands of exploded e-type X and Y. X + Y replaces Y. */
2042 for (i = M; i < NI; i++)
2044 a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
2055 /* Subtract significands of exploded e-type X and Y. Y - X replaces Y. */
2069 for (i = M; i < NI; i++)
2071 a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
2083 static UEMUSHORT equot[NI];
2087 /* Radix 2 shift-and-add versions of multiply and divide */
2090 /* Divide significands */
2094 UEMUSHORT den[], num[];
2104 for (i = M; i < NI; i++)
2109 /* Use faster compare and subtraction if denominator has only 15 bits of
2115 for (i = M + 3; i < NI; i++)
2120 if ((den[M + 1] & 1) != 0)
2128 for (i = 0; i < NBITS + 2; i++)
2146 /* The number of quotient bits to calculate is NBITS + 1 scaling guard
2147 bit + 1 roundoff bit. */
2152 for (i = 0; i < NBITS + 2; i++)
2154 if (ecmpm (den, num) <= 0)
2157 j = 1; /* quotient bit = 1 */
2171 /* test for nonzero remainder after roundoff bit */
2174 for (i = M; i < NI; i++)
2182 for (i = 0; i < NI; i++)
2188 /* Multiply significands */
2199 for (i = M; i < NI; i++)
2204 while (*p == 0) /* significand is not supposed to be zero */
2209 if ((*p & 0xff) == 0)
2217 for (i = 0; i < k; i++)
2221 /* remember if there were any nonzero bits shifted out */
2228 for (i = 0; i < NI; i++)
2231 /* return flag for lost nonzero bits */
2237 /* Radix 65536 versions of multiply and divide. */
2239 /* Multiply significand of e-type number B
2240 by 16-bit quantity A, return e-type result to C. */
2245 const UEMUSHORT b[];
2249 unsigned EMULONG carry;
2250 const UEMUSHORT *ps;
2252 unsigned EMULONG aa, m;
2261 for (i=M+1; i<NI; i++)
2271 m = (unsigned EMULONG) aa * *ps--;
2272 carry = (m & 0xffff) + *pp;
2273 *pp-- = (UEMUSHORT) carry;
2274 carry = (carry >> 16) + (m >> 16) + *pp;
2275 *pp = (UEMUSHORT) carry;
2276 *(pp-1) = carry >> 16;
2279 for (i=M; i<NI; i++)
2283 /* Divide significands of exploded e-types NUM / DEN. Neither the
2284 numerator NUM nor the denominator DEN is permitted to have its high guard
2289 const UEMUSHORT den[];
2294 unsigned EMULONG tnum;
2295 UEMUSHORT j, tdenm, tquot;
2296 UEMUSHORT tprod[NI+1];
2302 for (i=M; i<NI; i++)
2308 for (i=M; i<NI; i++)
2310 /* Find trial quotient digit (the radix is 65536). */
2311 tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
2313 /* Do not execute the divide instruction if it will overflow. */
2314 if ((tdenm * (unsigned long) 0xffff) < tnum)
2317 tquot = tnum / tdenm;
2318 /* Multiply denominator by trial quotient digit. */
2319 m16m ((unsigned int) tquot, den, tprod);
2320 /* The quotient digit may have been overestimated. */
2321 if (ecmpm (tprod, num) > 0)
2325 if (ecmpm (tprod, num) > 0)
2335 /* test for nonzero remainder after roundoff bit */
2338 for (i=M; i<NI; i++)
2345 for (i=0; i<NI; i++)
2351 /* Multiply significands of exploded e-type A and B, result in B. */
2355 const UEMUSHORT a[];
2360 UEMUSHORT pprod[NI];
2366 for (i=M; i<NI; i++)
2372 for (i=M+1; i<NI; i++)
2380 m16m ((unsigned int) *p--, b, pprod);
2381 eaddm (pprod, equot);
2387 for (i=0; i<NI; i++)
2390 /* return flag for lost nonzero bits */
2396 /* Normalize and round off.
2398 The internal format number to be rounded is S.
2399 Input LOST is 0 if the value is exact. This is the so-called sticky bit.
2401 Input SUBFLG indicates whether the number was obtained
2402 by a subtraction operation. In that case if LOST is nonzero
2403 then the number is slightly smaller than indicated.
2405 Input EXP is the biased exponent, which may be negative.
2406 the exponent field of S is ignored but is replaced by
2407 EXP as adjusted by normalization and rounding.
2409 Input RCNTRL is the rounding control. If it is nonzero, the
2410 returned value will be rounded to RNDPRC bits.
2412 For future reference: In order for emdnorm to round off denormal
2413 significands at the right point, the input exponent must be
2414 adjusted to be the actual value it would have after conversion to
2415 the final floating point type. This adjustment has been
2416 implemented for all type conversions (etoe53, etc.) and decimal
2417 conversions, but not for the arithmetic functions (eadd, etc.).
2418 Data types having standard 15-bit exponents are not affected by
2419 this, but SFmode and DFmode are affected. For example, ediv with
2420 rndprc = 24 will not round correctly to 24-bit precision if the
2421 result is denormal. */
2423 static int rlast = -1;
2425 static UEMUSHORT rmsk = 0;
2426 static UEMUSHORT rmbit = 0;
2427 static UEMUSHORT rebit = 0;
2429 static UEMUSHORT rbit[NI];
2432 emdnorm (s, lost, subflg, exp, rcntrl)
2445 /* a blank significand could mean either zero or infinity. */
2458 if ((j > NBITS) && (exp < 32767))
2466 if (exp > (EMULONG) (-NBITS - 1))
2479 /* Round off, unless told not to by rcntrl. */
2482 /* Set up rounding parameters if the control register changed. */
2483 if (rndprc != rlast)
2490 rw = NI - 1; /* low guard word */
2513 /* For DEC or IBM arithmetic */
2530 /* For C4x arithmetic */
2551 /* Shift down 1 temporarily if the data structure has an implied
2552 most significant bit and the number is denormal.
2553 Intel long double denormals also lose one bit of precision. */
2554 if ((exp <= 0) && (rndprc != NBITS)
2555 && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2557 lost |= s[NI - 1] & 1;
2560 /* Clear out all bits below the rounding bit,
2561 remembering in r if any were nonzero. */
2575 if ((r & rmbit) != 0)
2581 { /* round to even */
2582 if ((s[re] & rebit) == 0)
2595 /* Undo the temporary shift for denormal values. */
2596 if ((exp <= 0) && (rndprc != NBITS)
2597 && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2602 { /* overflow on roundoff */
2615 for (i = 2; i < NI - 1; i++)
2618 warning ("floating point overflow");
2622 for (i = M + 1; i < NI - 1; i++)
2625 if ((rndprc < 64) || (rndprc == 113))
2640 s[1] = (UEMUSHORT) exp;
2643 /* Subtract. C = B - A, all e type numbers. */
2645 static int subflg = 0;
2649 const UEMUSHORT *a, *b;
2664 /* Infinity minus infinity is a NaN.
2665 Test for subtracting infinities of the same sign. */
2666 if (eisinf (a) && eisinf (b)
2667 && ((eisneg (a) ^ eisneg (b)) == 0))
2669 mtherr ("esub", INVALID);
2678 /* Add. C = A + B, all e type. */
2682 const UEMUSHORT *a, *b;
2687 /* NaN plus anything is a NaN. */
2698 /* Infinity minus infinity is a NaN.
2699 Test for adding infinities of opposite signs. */
2700 if (eisinf (a) && eisinf (b)
2701 && ((eisneg (a) ^ eisneg (b)) != 0))
2703 mtherr ("esub", INVALID);
2712 /* Arithmetic common to both addition and subtraction. */
2716 const UEMUSHORT *a, *b;
2719 UEMUSHORT ai[NI], bi[NI], ci[NI];
2721 EMULONG lt, lta, ltb;
2742 /* compare exponents */
2747 { /* put the larger number in bi */
2757 if (lt < (EMULONG) (-NBITS - 1))
2758 goto done; /* answer same as larger addend */
2760 lost = eshift (ai, k); /* shift the smaller number down */
2764 /* exponents were the same, so must compare significands */
2767 { /* the numbers are identical in magnitude */
2768 /* if different signs, result is zero */
2774 /* if same sign, result is double */
2775 /* double denormalized tiny number */
2776 if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2781 /* add 1 to exponent unless both are zero! */
2782 for (j = 1; j < NI - 1; j++)
2798 bi[E] = (UEMUSHORT) ltb;
2802 { /* put the larger number in bi */
2818 emdnorm (bi, lost, subflg, ltb, !ROUND_TOWARDS_ZERO);
2824 /* Divide: C = B/A, all e type. */
2828 const UEMUSHORT *a, *b;
2831 UEMUSHORT ai[NI], bi[NI];
2833 EMULONG lt, lta, ltb;
2835 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2836 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
2837 sign = eisneg (a) ^ eisneg (b);
2840 /* Return any NaN input. */
2851 /* Zero over zero, or infinity over infinity, is a NaN. */
2852 if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
2853 || (eisinf (a) && eisinf (b)))
2855 mtherr ("ediv", INVALID);
2860 /* Infinity over anything else is infinity. */
2867 /* Anything else over infinity is zero. */
2879 { /* See if numerator is zero. */
2880 for (i = 1; i < NI - 1; i++)
2884 ltb -= enormlz (bi);
2894 { /* possible divide by zero */
2895 for (i = 1; i < NI - 1; i++)
2899 lta -= enormlz (ai);
2903 /* Divide by zero is not an invalid operation.
2904 It is a divide-by-zero operation! */
2906 mtherr ("ediv", SING);
2912 /* calculate exponent */
2913 lt = ltb - lta + EXONE;
2914 emdnorm (bi, i, 0, lt, !ROUND_TOWARDS_ZERO);
2921 && (ecmp (c, ezero) != 0)
2924 *(c+(NE-1)) |= 0x8000;
2926 *(c+(NE-1)) &= ~0x8000;
2929 /* Multiply e-types A and B, return e-type product C. */
2933 const UEMUSHORT *a, *b;
2936 UEMUSHORT ai[NI], bi[NI];
2938 EMULONG lt, lta, ltb;
2940 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2941 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
2942 sign = eisneg (a) ^ eisneg (b);
2945 /* NaN times anything is the same NaN. */
2956 /* Zero times infinity is a NaN. */
2957 if ((eisinf (a) && (ecmp (b, ezero) == 0))
2958 || (eisinf (b) && (ecmp (a, ezero) == 0)))
2960 mtherr ("emul", INVALID);
2965 /* Infinity times anything else is infinity. */
2967 if (eisinf (a) || eisinf (b))
2979 for (i = 1; i < NI - 1; i++)
2983 lta -= enormlz (ai);
2994 for (i = 1; i < NI - 1; i++)
2998 ltb -= enormlz (bi);
3007 /* Multiply significands */
3009 /* calculate exponent */
3010 lt = lta + ltb - (EXONE - 1);
3011 emdnorm (bi, j, 0, lt, !ROUND_TOWARDS_ZERO);
3018 && (ecmp (c, ezero) != 0)
3021 *(c+(NE-1)) |= 0x8000;
3023 *(c+(NE-1)) &= ~0x8000;
3026 /* Convert double precision PE to e-type Y. */
3030 const UEMUSHORT *pe;
3040 ibmtoe (pe, y, DFmode);
3045 c4xtoe (pe, y, HFmode);
3055 denorm = 0; /* flag if denormalized number */
3057 if (! REAL_WORDS_BIG_ENDIAN)
3063 yy[M] = (r & 0x0f) | 0x10;
3064 r &= ~0x800f; /* strip sign and 4 significand bits */
3069 if (! REAL_WORDS_BIG_ENDIAN)
3071 if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
3072 || (pe[1] != 0) || (pe[0] != 0))
3074 enan (y, yy[0] != 0);
3080 if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
3081 || (pe[2] != 0) || (pe[3] != 0))
3083 enan (y, yy[0] != 0);
3094 #endif /* INFINITY */
3096 /* If zero exponent, then the significand is denormalized.
3097 So take back the understood high significand bit. */
3108 if (! REAL_WORDS_BIG_ENDIAN)
3125 /* If zero exponent, then normalize the significand. */
3126 if ((k = enormlz (yy)) > NBITS)
3129 yy[E] -= (UEMUSHORT) (k - 1);
3132 #endif /* not C4X */
3133 #endif /* not IBM */
3134 #endif /* not DEC */
3137 /* Convert double extended precision float PE to e type Y. */
3141 const UEMUSHORT *pe;
3151 for (i = 0; i < NE - 5; i++)
3153 /* This precision is not ordinarily supported on DEC or IBM. */
3155 for (i = 0; i < 5; i++)
3159 p = &yy[0] + (NE - 1);
3162 for (i = 0; i < 5; i++)
3166 if (! REAL_WORDS_BIG_ENDIAN)
3168 for (i = 0; i < 5; i++)
3171 /* For denormal long double Intel format, shift significand up one
3172 -- but only if the top significand bit is zero. A top bit of 1
3173 is "pseudodenormal" when the exponent is zero. */
3174 if ((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
3186 p = &yy[0] + (NE - 1);
3187 #ifdef ARM_EXTENDED_IEEE_FORMAT
3188 /* For ARMs, the exponent is in the lowest 15 bits of the word. */
3189 *p-- = (e[0] & 0x8000) | (e[1] & 0x7ffff);
3195 for (i = 0; i < 4; i++)
3200 /* Point to the exponent field and check max exponent cases. */
3202 if ((*p & 0x7fff) == 0x7fff)
3205 if (! REAL_WORDS_BIG_ENDIAN)
3207 for (i = 0; i < 4; i++)
3209 if ((i != 3 && pe[i] != 0)
3210 /* Anything but 0x8000 here, including 0, is a NaN. */
3211 || (i == 3 && pe[i] != 0x8000))
3213 enan (y, (*p & 0x8000) != 0);
3220 #ifdef ARM_EXTENDED_IEEE_FORMAT
3221 for (i = 2; i <= 5; i++)
3225 enan (y, (*p & 0x8000) != 0);
3230 /* In Motorola extended precision format, the most significant
3231 bit of an infinity mantissa could be either 1 or 0. It is
3232 the lower order bits that tell whether the value is a NaN. */
3233 if ((pe[2] & 0x7fff) != 0)
3236 for (i = 3; i <= 5; i++)
3241 enan (y, (*p & 0x8000) != 0);
3245 #endif /* not ARM */
3254 #endif /* INFINITY */
3257 for (i = 0; i < NE; i++)
3261 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3262 /* Convert 128-bit long double precision float PE to e type Y. */
3266 const UEMUSHORT *pe;
3279 if (! REAL_WORDS_BIG_ENDIAN)
3291 if (! REAL_WORDS_BIG_ENDIAN)
3293 for (i = 0; i < 7; i++)
3297 enan (y, yy[0] != 0);
3304 for (i = 1; i < 8; i++)
3308 enan (y, yy[0] != 0);
3320 #endif /* INFINITY */
3324 if (! REAL_WORDS_BIG_ENDIAN)
3326 for (i = 0; i < 7; i++)
3332 for (i = 0; i < 7; i++)
3336 /* If denormal, remove the implied bit; else shift down 1. */
3350 /* Convert single precision float PE to e type Y. */
3354 const UEMUSHORT *pe;
3359 ibmtoe (pe, y, SFmode);
3365 c4xtoe (pe, y, QFmode);
3376 denorm = 0; /* flag if denormalized number */
3379 if (! REAL_WORDS_BIG_ENDIAN)
3389 yy[M] = (r & 0x7f) | 0200;
3390 r &= ~0x807f; /* strip sign and 7 significand bits */
3392 if (!LARGEST_EXPONENT_IS_NORMAL (32) && r == 0x7f80)
3395 if (REAL_WORDS_BIG_ENDIAN)
3397 if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
3399 enan (y, yy[0] != 0);
3405 if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
3407 enan (y, yy[0] != 0);
3418 #endif /* INFINITY */
3420 /* If zero exponent, then the significand is denormalized.
3421 So take back the understood high significand bit. */
3434 if (! REAL_WORDS_BIG_ENDIAN)
3444 { /* if zero exponent, then normalize the significand */
3445 if ((k = enormlz (yy)) > NBITS)
3448 yy[E] -= (UEMUSHORT) (k - 1);
3451 #endif /* not C4X */
3452 #endif /* not IBM */
3455 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3456 /* Convert e-type X to IEEE 128-bit long double format E. */
3470 make_nan (e, eisneg (x), TFmode);
3475 exp = (EMULONG) xi[E];
3480 /* round off to nearest or even */
3483 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
3491 /* Convert exploded e-type X, that has already been rounded to
3492 113-bit precision, to IEEE 128-bit long double format Y. */
3504 make_nan (b, eiisneg (a), TFmode);
3509 if (REAL_WORDS_BIG_ENDIAN)
3512 q = b + 7; /* point to output exponent */
3514 /* If not denormal, delete the implied bit. */
3519 /* combine sign and exponent */
3521 if (REAL_WORDS_BIG_ENDIAN)
3524 *q++ = *p++ | 0x8000;
3531 *q-- = *p++ | 0x8000;
3535 /* skip over guard word */
3537 /* move the significand */
3538 if (REAL_WORDS_BIG_ENDIAN)
3540 for (i = 0; i < 7; i++)
3545 for (i = 0; i < 7; i++)
3551 /* Convert e-type X to IEEE double extended format E. */
3565 make_nan (e, eisneg (x), XFmode);
3570 /* adjust exponent for offset */
3571 exp = (EMULONG) xi[E];
3576 /* round off to nearest or even */
3579 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
3587 /* Convert exploded e-type X, that has already been rounded to
3588 64-bit precision, to IEEE double extended format Y. */
3600 make_nan (b, eiisneg (a), XFmode);
3604 /* Shift denormal long double Intel format significand down one bit. */
3605 if ((a[E] == 0) && ! REAL_WORDS_BIG_ENDIAN)
3615 if (REAL_WORDS_BIG_ENDIAN)
3619 q = b + 4; /* point to output exponent */
3620 /* Clear the last two bytes of 12-byte Intel format. q is pointing
3621 into an array of size 6 (e.g. x[NE]), so the last two bytes are
3622 always there, and there are never more bytes, even when we are using
3623 INTEL_EXTENDED_IEEE_FORMAT. */
3628 /* combine sign and exponent */
3632 *q++ = *p++ | 0x8000;
3639 *q-- = *p++ | 0x8000;
3644 if (REAL_WORDS_BIG_ENDIAN)
3646 #ifdef ARM_EXTENDED_IEEE_FORMAT
3647 /* The exponent is in the lowest 15 bits of the first word. */
3648 *q++ = i ? 0x8000 : 0;
3652 *q++ = *p++ | 0x8000;
3661 *q-- = *p++ | 0x8000;
3666 /* skip over guard word */
3668 /* move the significand */
3670 for (i = 0; i < 4; i++)
3674 for (i = 0; i < 4; i++)
3678 if (REAL_WORDS_BIG_ENDIAN)
3680 for (i = 0; i < 4; i++)
3688 /* Intel long double infinity significand. */
3696 for (i = 0; i < 4; i++)
3702 /* e type to double precision. */
3705 /* Convert e-type X to DEC-format double E. */
3712 etodec (x, e); /* see etodec.c */
3715 /* Convert exploded e-type X, that has already been rounded to
3716 56-bit double precision, to DEC double Y. */
3727 /* Convert e-type X to IBM 370-format double E. */
3734 etoibm (x, e, DFmode);
3737 /* Convert exploded e-type X, that has already been rounded to
3738 56-bit precision, to IBM 370 double Y. */
3744 toibm (x, y, DFmode);
3747 #else /* it's neither DEC nor IBM */
3749 /* Convert e-type X to C4X-format long double E. */
3756 etoc4x (x, e, HFmode);
3759 /* Convert exploded e-type X, that has already been rounded to
3760 56-bit precision, to IBM 370 double Y. */
3766 toc4x (x, y, HFmode);
3769 #else /* it's neither DEC nor IBM nor C4X */
3771 /* Convert e-type X to IEEE double E. */
3785 make_nan (e, eisneg (x), DFmode);
3790 /* adjust exponent for offsets */
3791 exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
3796 /* round off to nearest or even */
3799 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
3807 /* Convert exploded e-type X, that has already been rounded to
3808 53-bit precision, to IEEE double Y. */
3820 make_nan (y, eiisneg (x), DFmode);
3824 if (LARGEST_EXPONENT_IS_NORMAL (64) && x[1] > 2047)
3826 saturate (y, eiisneg (x), 64, 1);
3831 if (! REAL_WORDS_BIG_ENDIAN)
3834 *y = 0; /* output high order */
3836 *y = 0x8000; /* output sign bit */
3839 if (i >= (unsigned int) 2047)
3841 /* Saturate at largest number less than infinity. */
3844 if (! REAL_WORDS_BIG_ENDIAN)
3858 *y |= (UEMUSHORT) 0x7fef;
3859 if (! REAL_WORDS_BIG_ENDIAN)
3884 i |= *p++ & (UEMUSHORT) 0x0f; /* *p = xi[M] */
3885 *y |= (UEMUSHORT) i; /* high order output already has sign bit set */
3886 if (! REAL_WORDS_BIG_ENDIAN)
3901 #endif /* not C4X */
3902 #endif /* not IBM */
3903 #endif /* not DEC */
3907 /* e type to single precision. */
3910 /* Convert e-type X to IBM 370 float E. */
3917 etoibm (x, e, SFmode);
3920 /* Convert exploded e-type X, that has already been rounded to
3921 float precision, to IBM 370 float Y. */
3927 toibm (x, y, SFmode);
3933 /* Convert e-type X to C4X float E. */
3940 etoc4x (x, e, QFmode);
3943 /* Convert exploded e-type X, that has already been rounded to
3944 float precision, to IBM 370 float Y. */
3950 toc4x (x, y, QFmode);
3955 /* Convert e-type X to IEEE float E. DEC float is the same as IEEE float. */
3969 make_nan (e, eisneg (x), SFmode);
3974 /* adjust exponent for offsets */
3975 exp = (EMULONG) xi[E] - (EXONE - 0177);
3980 /* round off to nearest or even */
3983 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
3991 /* Convert exploded e-type X, that has already been rounded to
3992 float precision, to IEEE float Y. */
4004 make_nan (y, eiisneg (x), SFmode);
4008 if (LARGEST_EXPONENT_IS_NORMAL (32) && x[1] > 255)
4010 saturate (y, eiisneg (x), 32, 1);
4015 if (! REAL_WORDS_BIG_ENDIAN)
4021 *y = 0; /* output high order */
4023 *y = 0x8000; /* output sign bit */
4026 /* Handle overflow cases. */
4027 if (!LARGEST_EXPONENT_IS_NORMAL (32) && i >= 255)
4030 *y |= (UEMUSHORT) 0x7f80;
4035 if (! REAL_WORDS_BIG_ENDIAN)
4043 #else /* no INFINITY */
4044 *y |= (UEMUSHORT) 0x7f7f;
4049 if (! REAL_WORDS_BIG_ENDIAN)
4060 #endif /* no INFINITY */
4072 i |= *p++ & (UEMUSHORT) 0x7f; /* *p = xi[M] */
4073 /* High order output already has sign bit set. */
4079 if (! REAL_WORDS_BIG_ENDIAN)
4088 #endif /* not C4X */
4089 #endif /* not IBM */
4091 /* Compare two e type numbers.
4095 -2 if either a or b is a NaN. */
4099 const UEMUSHORT *a, *b;
4101 UEMUSHORT ai[NI], bi[NI];
4107 if (eisnan (a) || eisnan (b))
4116 { /* the signs are different */
4118 for (i = 1; i < NI - 1; i++)
4132 /* both are the same sign */
4147 return (0); /* equality */
4151 if (*(--p) > *(--q))
4152 return (msign); /* p is bigger */
4154 return (-msign); /* p is littler */
4158 /* Find e-type nearest integer to X, as floor (X + 0.5). */
4170 /* Convert HOST_WIDE_INT LP to e type Y. */
4174 const HOST_WIDE_INT *lp;
4178 unsigned HOST_WIDE_INT ll;
4184 /* make it positive */
4185 ll = (unsigned HOST_WIDE_INT) (-(*lp));
4186 yi[0] = 0xffff; /* put correct sign in the e type number */
4190 ll = (unsigned HOST_WIDE_INT) (*lp);
4192 /* move the long integer to yi significand area */
4193 #if HOST_BITS_PER_WIDE_INT == 64
4194 yi[M] = (UEMUSHORT) (ll >> 48);
4195 yi[M + 1] = (UEMUSHORT) (ll >> 32);
4196 yi[M + 2] = (UEMUSHORT) (ll >> 16);
4197 yi[M + 3] = (UEMUSHORT) ll;
4198 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
4200 yi[M] = (UEMUSHORT) (ll >> 16);
4201 yi[M + 1] = (UEMUSHORT) ll;
4202 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
4205 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4206 ecleaz (yi); /* it was zero */
4208 yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
4209 emovo (yi, y); /* output the answer */
4212 /* Convert unsigned HOST_WIDE_INT LP to e type Y. */
4216 const unsigned HOST_WIDE_INT *lp;
4220 unsigned HOST_WIDE_INT ll;
4226 /* move the long integer to ayi significand area */
4227 #if HOST_BITS_PER_WIDE_INT == 64
4228 yi[M] = (UEMUSHORT) (ll >> 48);
4229 yi[M + 1] = (UEMUSHORT) (ll >> 32);
4230 yi[M + 2] = (UEMUSHORT) (ll >> 16);
4231 yi[M + 3] = (UEMUSHORT) ll;
4232 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
4234 yi[M] = (UEMUSHORT) (ll >> 16);
4235 yi[M + 1] = (UEMUSHORT) ll;
4236 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
4239 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4240 ecleaz (yi); /* it was zero */
4242 yi[E] -= (UEMUSHORT) k; /* subtract shift count from exponent */
4243 emovo (yi, y); /* output the answer */
4247 /* Find signed HOST_WIDE_INT integer I and floating point fractional
4248 part FRAC of e-type (packed internal format) floating point input X.
4249 The integer output I has the sign of the input, except that
4250 positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
4251 The output e-type fraction FRAC is the positive fractional
4262 unsigned HOST_WIDE_INT ll;
4265 k = (int) xi[E] - (EXONE - 1);
4268 /* if exponent <= 0, integer = 0 and real output is fraction */
4273 if (k > (HOST_BITS_PER_WIDE_INT - 1))
4275 /* long integer overflow: output large integer
4276 and correct fraction */
4278 *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
4281 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
4282 /* In this case, let it overflow and convert as if unsigned. */
4283 euifrac (x, &ll, frac);
4284 *i = (HOST_WIDE_INT) ll;
4287 /* In other cases, return the largest positive integer. */
4288 *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
4293 warning ("overflow on truncation to integer");
4297 /* Shift more than 16 bits: first shift up k-16 mod 16,
4298 then shift up by 16's. */
4299 j = k - ((k >> 4) << 4);
4306 ll = (ll << 16) | xi[M];
4308 while ((k -= 16) > 0);
4315 /* shift not more than 16 bits */
4317 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4324 if ((k = enormlz (xi)) > NBITS)
4327 xi[E] -= (UEMUSHORT) k;
4333 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
4334 FRAC of e-type X. A negative input yields integer output = 0 but
4335 correct fraction. */
4338 euifrac (x, i, frac)
4340 unsigned HOST_WIDE_INT *i;
4343 unsigned HOST_WIDE_INT ll;
4348 k = (int) xi[E] - (EXONE - 1);
4351 /* if exponent <= 0, integer = 0 and argument is fraction */
4356 if (k > HOST_BITS_PER_WIDE_INT)
4358 /* Long integer overflow: output large integer
4359 and correct fraction.
4360 Note, the BSD MicroVAX compiler says that ~(0UL)
4361 is a syntax error. */
4365 warning ("overflow on truncation to unsigned integer");
4369 /* Shift more than 16 bits: first shift up k-16 mod 16,
4370 then shift up by 16's. */
4371 j = k - ((k >> 4) << 4);
4378 ll = (ll << 16) | xi[M];
4380 while ((k -= 16) > 0);
4385 /* shift not more than 16 bits */
4387 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4390 if (xi[0]) /* A negative value yields unsigned integer 0. */
4396 if ((k = enormlz (xi)) > NBITS)
4399 xi[E] -= (UEMUSHORT) k;
4404 /* Shift the significand of exploded e-type X up or down by SC bits. */
4425 lost |= *p; /* remember lost bits */
4466 return ((int) lost);
4469 /* Shift normalize the significand area of exploded e-type X.
4470 Return the shift count (up = positive). */
4485 return (0); /* already normalized */
4491 /* With guard word, there are NBITS+16 bits available.
4492 Return true if all are zero. */
4496 /* see if high byte is zero */
4497 while ((*p & 0xff00) == 0)
4502 /* now shift 1 bit at a time */
4503 while ((*p & 0x8000) == 0)
4509 mtherr ("enormlz", UNDERFLOW);
4515 /* Normalize by shifting down out of the high guard word
4516 of the significand */
4531 mtherr ("enormlz", OVERFLOW);
4538 /* Powers of ten used in decimal <-> binary conversions. */
4543 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
4544 static const UEMUSHORT etens[NTEN + 1][NE] =
4546 {0x6576, 0x4a92, 0x804a, 0x153f,
4547 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4548 {0x6a32, 0xce52, 0x329a, 0x28ce,
4549 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4550 {0x526c, 0x50ce, 0xf18b, 0x3d28,
4551 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4552 {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4553 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4554 {0x851e, 0xeab7, 0x98fe, 0x901b,
4555 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4556 {0x0235, 0x0137, 0x36b1, 0x336c,
4557 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4558 {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4559 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4560 {0x0000, 0x0000, 0x0000, 0x0000,
4561 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4562 {0x0000, 0x0000, 0x0000, 0x0000,
4563 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4564 {0x0000, 0x0000, 0x0000, 0x0000,
4565 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4566 {0x0000, 0x0000, 0x0000, 0x0000,
4567 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4568 {0x0000, 0x0000, 0x0000, 0x0000,
4569 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4570 {0x0000, 0x0000, 0x0000, 0x0000,
4571 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4574 static const UEMUSHORT emtens[NTEN + 1][NE] =
4576 {0x2030, 0xcffc, 0xa1c3, 0x8123,
4577 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4578 {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4579 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4580 {0xf53f, 0xf698, 0x6bd3, 0x0158,
4581 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4582 {0xe731, 0x04d4, 0xe3f2, 0xd332,
4583 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4584 {0xa23e, 0x5308, 0xfefb, 0x1155,
4585 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4586 {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4587 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4588 {0x2a20, 0x6224, 0x47b3, 0x98d7,
4589 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4590 {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4591 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4592 {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4593 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4594 {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4595 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4596 {0xc155, 0xa4a8, 0x404e, 0x6113,
4597 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4598 {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4599 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4600 {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4601 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4604 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4605 static const UEMUSHORT etens[NTEN + 1][NE] =
4607 {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4608 {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4609 {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4610 {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4611 {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4612 {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4613 {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4614 {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4615 {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4616 {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4617 {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4618 {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4619 {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4622 static const UEMUSHORT emtens[NTEN + 1][NE] =
4624 {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4625 {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4626 {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4627 {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4628 {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4629 {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4630 {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4631 {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4632 {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4633 {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4634 {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4635 {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4636 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4641 /* Convert float value X to ASCII string STRING with NDIG digits after
4642 the decimal point. */
4645 e24toasc (x, string, ndigs)
4646 const UEMUSHORT x[];
4653 etoasc (w, string, ndigs);
4656 /* Convert double value X to ASCII string STRING with NDIG digits after
4657 the decimal point. */
4660 e53toasc (x, string, ndigs)
4661 const UEMUSHORT x[];
4668 etoasc (w, string, ndigs);
4671 /* Convert double extended value X to ASCII string STRING with NDIG digits
4672 after the decimal point. */
4675 e64toasc (x, string, ndigs)
4676 const UEMUSHORT x[];
4683 etoasc (w, string, ndigs);
4686 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4687 after the decimal point. */
4690 e113toasc (x, string, ndigs)
4691 const UEMUSHORT x[];
4698 etoasc (w, string, ndigs);
4702 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4703 the decimal point. */
4705 static char wstring[80]; /* working storage for ASCII output */
4708 etoasc (x, string, ndigs)
4709 const UEMUSHORT x[];
4714 UEMUSHORT y[NI], t[NI], u[NI], w[NI];
4715 const UEMUSHORT *p, *r, *ten;
4717 int i, j, k, expon, rndsav;
4730 sprintf (wstring, " NaN ");
4734 rndprc = NBITS; /* set to full precision */
4735 emov (x, y); /* retain external format */
4736 if (y[NE - 1] & 0x8000)
4739 y[NE - 1] &= 0x7fff;
4746 ten = &etens[NTEN][0];
4748 /* Test for zero exponent */
4751 for (k = 0; k < NE - 1; k++)
4754 goto tnzro; /* denormalized number */
4756 goto isone; /* valid all zeros */
4760 /* Test for infinity. */
4761 if (y[NE - 1] == 0x7fff)
4764 sprintf (wstring, " -Infinity ");
4766 sprintf (wstring, " Infinity ");
4770 /* Test for exponent nonzero but significand denormalized.
4771 * This is an error condition.
4773 if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4775 mtherr ("etoasc", DOMAIN);
4776 sprintf (wstring, "NaN");
4780 /* Compare to 1.0 */
4789 { /* Number is greater than 1 */
4790 /* Convert significand to an integer and strip trailing decimal zeros. */
4792 u[NE - 1] = EXONE + NBITS - 1;
4794 p = &etens[NTEN - 4][0];
4800 for (j = 0; j < NE - 1; j++)
4813 /* Rescale from integer significand */
4814 u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
4816 /* Find power of 10 */
4820 /* An unordered compare result shouldn't happen here. */
4821 while (ecmp (ten, u) <= 0)
4823 if (ecmp (p, u) <= 0)
4836 { /* Number is less than 1.0 */
4837 /* Pad significand with trailing decimal zeros. */
4840 while ((y[NE - 2] & 0x8000) == 0)
4849 for (i = 0; i < NDEC + 1; i++)
4851 if ((w[NI - 1] & 0x7) != 0)
4853 /* multiply by 10 */
4866 if (eone[NE - 1] <= u[1])
4878 while (ecmp (eone, w) > 0)
4880 if (ecmp (p, w) >= 0)
4895 /* Find the first (leading) digit. */
4901 digit = equot[NI - 1];
4902 while ((digit == 0) && (ecmp (y, ezero) != 0))
4910 digit = equot[NI - 1];
4918 /* Examine number of digits requested by caller. */
4936 *s++ = (char) digit + '0';
4939 /* Generate digits after the decimal point. */
4940 for (k = 0; k <= ndigs; k++)
4942 /* multiply current number by 10, without normalizing */
4949 *s++ = (char) equot[NI - 1] + '0';
4951 digit = equot[NI - 1];
4954 /* round off the ASCII string */
4957 /* Test for critical rounding case in ASCII output. */
4961 if (ecmp (t, ezero) != 0)
4962 goto roun; /* round to nearest */
4964 if ((*(s - 1) & 1) == 0)
4965 goto doexp; /* round to even */
4968 /* Round up and propagate carry-outs */
4972 /* Carry out to most significant digit? */
4979 /* Most significant digit carries to 10? */
4987 /* Round up and carry out from less significant digits */
4999 sprintf (ss, "e+%d", expon);
5001 sprintf (ss, "e%d", expon);
5003 sprintf (ss, "e%d", expon);
5006 /* copy out the working string */
5009 while (*ss == ' ') /* strip possible leading space */
5011 while ((*s++ = *ss++) != '\0')
5016 /* Convert ASCII string to floating point.
5018 Numeric input is a free format decimal number of any length, with
5019 or without decimal point. Entering E after the number followed by an
5020 integer number causes the second number to be interpreted as a power of
5021 10 to be multiplied by the first number (i.e., "scientific" notation). */
5023 /* Convert ASCII string S to single precision float value Y. */
5034 /* Convert ASCII string S to double precision value Y. */
5041 #if defined(DEC) || defined(IBM)
5053 /* Convert ASCII string S to double extended value Y. */
5063 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
5064 /* Convert ASCII string S to 128-bit long double Y. */
5071 asctoeg (s, y, 113);
5075 /* Convert ASCII string S to e type Y. */
5082 asctoeg (s, y, NBITS);
5085 /* Convert ASCII string SS to e type Y, with a specified rounding precision
5086 of OPREC bits. BASE is 16 for C99 hexadecimal floating constants. */
5089 asctoeg (ss, y, oprec)
5094 UEMUSHORT yy[NI], xt[NI], tt[NI];
5095 int esign, decflg, sgnflg, nexp, exp, prec, lost;
5096 int i, k, trail, c, rndsav;
5099 char *sp, *s, *lstr;
5102 /* Copy the input string. */
5103 lstr = (char *) alloca (strlen (ss) + 1);
5105 while (*ss == ' ') /* skip leading spaces */
5109 while ((*sp++ = *ss++) != '\0')
5113 if (s[0] == '0' && (s[1] == 'x' || s[1] == 'X'))
5120 rndprc = NBITS; /* Set to full precision */
5133 if ((k >= 0) && (k < base))
5135 /* Ignore leading zeros */
5136 if ((prec == 0) && (decflg == 0) && (k == 0))
5138 /* Identify and strip trailing zeros after the decimal point. */
5139 if ((trail == 0) && (decflg != 0))
5142 while (ISDIGIT (*sp) || (base == 16 && ISXDIGIT (*sp)))
5144 /* Check for syntax error */
5146 if ((base != 10 || ((c != 'e') && (c != 'E')))
5147 && (base != 16 || ((c != 'p') && (c != 'P')))
5149 && (c != '\n') && (c != '\r') && (c != ' ')
5151 goto unexpected_char_error;
5160 /* If enough digits were given to more than fill up the yy register,
5161 continuing until overflow into the high guard word yy[2]
5162 guarantees that there will be a roundoff bit at the top
5163 of the low guard word after normalization. */
5170 nexp += 4; /* count digits after decimal point */
5172 eshup1 (yy); /* multiply current number by 16 */
5180 nexp += 1; /* count digits after decimal point */
5182 eshup1 (yy); /* multiply current number by 10 */
5188 /* Insert the current digit. */
5190 xt[NI - 2] = (UEMUSHORT) k;
5195 /* Mark any lost non-zero digit. */
5197 /* Count lost digits before the decimal point. */
5219 case '.': /* decimal point */
5221 goto unexpected_char_error;
5227 goto unexpected_char_error;
5232 goto unexpected_char_error;
5245 unexpected_char_error:
5249 mtherr ("asctoe", DOMAIN);
5258 /* Exponent interpretation */
5260 /* 0.0eXXX is zero, regardless of XXX. Check for the 0.0. */
5261 for (k = 0; k < NI; k++)
5272 /* check for + or - */
5280 while (ISDIGIT (*s))
5289 if ((exp > MAXDECEXP) && (base == 10))
5293 yy[E] = 0x7fff; /* infinity */
5296 if ((exp < MINDECEXP) && (base == 10))
5306 /* Base 16 hexadecimal floating constant. */
5307 if ((k = enormlz (yy)) > NBITS)
5312 /* Adjust the exponent. NEXP is the number of hex digits,
5313 EXP is a power of 2. */
5314 lexp = (EXONE - 1 + NBITS) - k + yy[E] + exp - nexp;
5324 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
5325 while ((nexp > 0) && (yy[2] == 0))
5337 if ((k = enormlz (yy)) > NBITS)
5342 lexp = (EXONE - 1 + NBITS) - k;
5343 emdnorm (yy, lost, 0, lexp, 64);
5346 /* Convert to external format:
5348 Multiply by 10**nexp. If precision is 64 bits,
5349 the maximum relative error incurred in forming 10**n
5350 for 0 <= n <= 324 is 8.2e-20, at 10**180.
5351 For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5352 For 0 >= n >= -999, it is -1.55e-19 at 10**-435. */
5367 /* Punt. Can't handle this without 2 divides. */
5368 emovi (etens[0], tt);
5381 emul (etens[i], xt, xt);
5385 while (exp <= MAXP);
5404 /* Round and convert directly to the destination type */
5406 lexp -= EXONE - 0x3ff;
5408 else if (oprec == 24 || oprec == 32)
5409 lexp -= (EXONE - 0x7f);
5412 else if (oprec == 24 || oprec == 56)
5413 lexp -= EXONE - (0x41 << 2);
5415 else if (oprec == 24)
5416 lexp -= EXONE - 0177;
5420 else if (oprec == 56)
5421 lexp -= EXONE - 0201;
5424 emdnorm (yy, lost, 0, lexp, 64);
5434 todec (yy, y); /* see etodec.c */
5439 toibm (yy, y, DFmode);
5444 toc4x (yy, y, HFmode);
5457 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
5470 /* Return Y = largest integer not greater than X (truncated toward minus
5473 static const UEMUSHORT bmask[] =
5496 const UEMUSHORT x[];
5503 emov (x, f); /* leave in external format */
5504 expon = (int) f[NE - 1];
5505 e = (expon & 0x7fff) - (EXONE - 1);
5511 /* number of bits to clear out */
5523 /* clear the remaining bits */
5525 /* truncate negatives toward minus infinity */
5528 if ((UEMUSHORT) expon & (UEMUSHORT) 0x8000)
5530 for (i = 0; i < NE - 1; i++)
5543 /* Return S and EXP such that S * 2^EXP = X and .5 <= S < 1.
5544 For example, 1.1 = 0.55 * 2^1. */
5548 const UEMUSHORT x[];
5556 /* Handle denormalized numbers properly using long integer exponent. */
5557 li = (EMULONG) ((EMUSHORT) xi[1]);
5565 *exp = (int) (li - 0x3ffe);
5569 /* Return e type Y = X * 2^PWR2. */
5573 const UEMUSHORT x[];
5585 emdnorm (xi, i, i, li, !ROUND_TOWARDS_ZERO);
5591 /* C = remainder after dividing B by A, all e type values.
5592 Least significant integer quotient bits left in EQUOT. */
5596 const UEMUSHORT a[], b[];
5599 UEMUSHORT den[NI], num[NI];
5603 || (ecmp (a, ezero) == 0)
5611 if (ecmp (a, ezero) == 0)
5613 mtherr ("eremain", SING);
5619 eiremain (den, num);
5620 /* Sign of remainder = sign of quotient */
5629 /* Return quotient of exploded e-types NUM / DEN in EQUOT,
5630 remainder in NUM. */
5634 UEMUSHORT den[], num[];
5640 ld -= enormlz (den);
5642 ln -= enormlz (num);
5646 if (ecmpm (den, num) <= 0)
5658 emdnorm (num, 0, 0, ln, 0);
5661 /* Report an error condition CODE encountered in function NAME.
5663 Mnemonic Value Significance
5665 DOMAIN 1 argument domain error
5666 SING 2 function singularity
5667 OVERFLOW 3 overflow range error
5668 UNDERFLOW 4 underflow range error
5669 TLOSS 5 total loss of precision
5670 PLOSS 6 partial loss of precision
5671 INVALID 7 NaN - producing operation
5672 EDOM 33 Unix domain error code
5673 ERANGE 34 Unix range error code
5675 The order of appearance of the following messages is bound to the
5676 error codes defined above. */
5686 /* The string passed by the calling program is supposed to be the
5687 name of the function in which the error occurred.
5688 The code argument selects which error message string will be printed. */
5690 if (strcmp (name, "esub") == 0)
5691 name = "subtraction";
5692 else if (strcmp (name, "ediv") == 0)
5694 else if (strcmp (name, "emul") == 0)
5695 name = "multiplication";
5696 else if (strcmp (name, "enormlz") == 0)
5697 name = "normalization";
5698 else if (strcmp (name, "etoasc") == 0)
5699 name = "conversion to text";
5700 else if (strcmp (name, "asctoe") == 0)
5702 else if (strcmp (name, "eremain") == 0)
5704 else if (strcmp (name, "esqrt") == 0)
5705 name = "square root";
5710 case DOMAIN: warning ("%s: argument domain error" , name); break;
5711 case SING: warning ("%s: function singularity" , name); break;
5712 case OVERFLOW: warning ("%s: overflow range error" , name); break;
5713 case UNDERFLOW: warning ("%s: underflow range error" , name); break;
5714 case TLOSS: warning ("%s: total loss of precision" , name); break;
5715 case PLOSS: warning ("%s: partial loss of precision", name); break;
5716 case INVALID: warning ("%s: NaN - producing operation", name); break;
5721 /* Set global error message word */
5726 /* Convert DEC double precision D to e type E. */
5736 ecleaz (y); /* start with a zero */
5737 p = y; /* point to our number */
5738 r = *d; /* get DEC exponent word */
5739 if (*d & (unsigned int) 0x8000)
5740 *p = 0xffff; /* fill in our sign */
5741 ++p; /* bump pointer to our exponent word */
5742 r &= 0x7fff; /* strip the sign bit */
5743 if (r == 0) /* answer = 0 if high order DEC word = 0 */
5747 r >>= 7; /* shift exponent word down 7 bits */
5748 r += EXONE - 0201; /* subtract DEC exponent offset */
5749 /* add our e type exponent offset */
5750 *p++ = r; /* to form our exponent */
5752 r = *d++; /* now do the high order mantissa */
5753 r &= 0177; /* strip off the DEC exponent and sign bits */
5754 r |= 0200; /* the DEC understood high order mantissa bit */
5755 *p++ = r; /* put result in our high guard word */
5757 *p++ = *d++; /* fill in the rest of our mantissa */
5761 eshdn8 (y); /* shift our mantissa down 8 bits */
5766 /* Convert e type X to DEC double precision D. */
5778 /* Adjust exponent for offsets. */
5779 exp = (EMULONG) xi[E] - (EXONE - 0201);
5780 /* Round off to nearest or even. */
5783 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
5788 /* Convert exploded e-type X, that has already been rounded to
5789 56-bit precision, to DEC format double Y. */
5835 /* Convert IBM single/double precision to e type. */
5841 enum machine_mode mode;
5846 ecleaz (y); /* start with a zero */
5847 p = y; /* point to our number */
5848 r = *d; /* get IBM exponent word */
5849 if (*d & (unsigned int) 0x8000)
5850 *p = 0xffff; /* fill in our sign */
5851 ++p; /* bump pointer to our exponent word */
5852 r &= 0x7f00; /* strip the sign bit */
5853 r >>= 6; /* shift exponent word down 6 bits */
5854 /* in fact shift by 8 right and 2 left */
5855 r += EXONE - (0x41 << 2); /* subtract IBM exponent offset */
5856 /* add our e type exponent offset */
5857 *p++ = r; /* to form our exponent */
5859 *p++ = *d++ & 0xff; /* now do the high order mantissa */
5860 /* strip off the IBM exponent and sign bits */
5861 if (mode != SFmode) /* there are only 2 words in SFmode */
5863 *p++ = *d++; /* fill in the rest of our mantissa */
5868 if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5871 y[E] -= 5 + enormlz (y); /* now normalise the mantissa */
5872 /* handle change in RADIX */
5878 /* Convert e type to IBM single/double precision. */
5884 enum machine_mode mode;
5891 exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2)); /* adjust exponent for offsets */
5892 /* round off to nearest or even */
5895 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
5897 toibm (xi, d, mode);
5903 enum machine_mode mode;
5956 /* Convert C4X single/double precision to e type. */
5962 enum machine_mode mode;
5980 /* Short-circuit the zero case. */
5981 if ((dn[0] == 0x8000)
5982 && (dn[1] == 0x0000)
5983 && ((mode == QFmode) || ((dn[2] == 0x0000) && (dn[3] == 0x0000))))
5994 ecleaz (y); /* start with a zero */
5995 r = dn[0]; /* get sign/exponent part */
5996 if (r & (unsigned int) 0x0080)
5998 y[0] = 0xffff; /* fill in our sign */
6004 r >>= 8; /* Shift exponent word down 8 bits. */
6005 if (r & 0x80) /* Make the exponent negative if it is. */
6006 r = r | (~0 & ~0xff);
6010 /* Now do the high order mantissa. We don't "or" on the high bit
6011 because it is 2 (not 1) and is handled a little differently
6013 y[M] = dn[0] & 0x7f;
6016 if (mode != QFmode) /* There are only 2 words in QFmode. */
6018 y[M+2] = dn[2]; /* Fill in the rest of our mantissa. */
6026 /* Now do the two's complement on the data. */
6028 carry = 1; /* Initially add 1 for the two's complement. */
6029 for (i=size + M; i > M; i--)
6031 if (carry && (y[i] == 0x0000))
6032 /* We overflowed into the next word, carry is the same. */
6033 y[i] = carry ? 0x0000 : 0xffff;
6036 /* No overflow, just invert and add carry. */
6037 y[i] = ((~y[i]) + carry) & 0xffff;
6052 /* Add our e type exponent offset to form our exponent. */
6056 /* Now do the high order mantissa strip off the exponent and sign
6057 bits and add the high 1 bit. */
6058 y[M] = (dn[0] & 0x7f) | 0x80;
6061 if (mode != QFmode) /* There are only 2 words in QFmode. */
6063 y[M+2] = dn[2]; /* Fill in the rest of our mantissa. */
6073 /* Convert e type to C4X single/double precision. */
6079 enum machine_mode mode;
6087 /* Adjust exponent for offsets. */
6088 exp = (EMULONG) xi[E] - (EXONE - 0x7f);
6090 /* Round off to nearest or even. */
6092 rndprc = mode == QFmode ? 24 : 32;
6093 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
6095 toc4x (xi, d, mode);
6101 enum machine_mode mode;
6107 /* Short-circuit the zero case */
6108 if ((x[0] == 0) /* Zero exponent and sign */
6110 && (x[M] == 0) /* The rest is for zero mantissa */
6112 /* Only check for double if necessary */
6113 && ((mode == QFmode) || ((x[M+2] == 0) && (x[M+3] == 0))))
6115 /* We have a zero. Put it into the output and return. */
6128 /* Negative number require a two's complement conversion of the
6134 i = ((int) x[1]) - 0x7f;
6136 /* Now add 1 to the inverted data to do the two's complement. */
6145 x[v] = carry ? 0x0000 : 0xffff;
6148 x[v] = ((~x[v]) + carry) & 0xffff;
6154 /* The following is a special case. The C4X negative float requires
6155 a zero in the high bit (because the format is (2 - x) x 2^m), so
6156 if a one is in that bit, we have to shift left one to get rid
6157 of it. This only occurs if the number is -1 x 2^m. */
6158 if (x[M+1] & 0x8000)
6160 /* This is the case of -1 x 2^m, we have to rid ourselves of the
6161 high sign bit and shift the exponent. */
6167 i = ((int) x[1]) - 0x7f;
6169 if ((i < -128) || (i > 127))
6177 y[3] = (y[1] << 8) | ((y[2] >> 8) & 0xff);
6178 y[2] = (y[0] << 8) | ((y[1] >> 8) & 0xff);
6186 y[0] |= ((i & 0xff) << 8);
6190 y[0] |= x[M] & 0x7f;
6196 y[3] = (y[1] << 8) | ((y[2] >> 8) & 0xff);
6197 y[2] = (y[0] << 8) | ((y[1] >> 8) & 0xff);
6202 /* Output a binary NaN bit pattern in the target machine's format. */
6204 /* If special NaN bit patterns are required, define them in tm.h
6205 as arrays of unsigned 16-bit shorts. Otherwise, use the default
6211 static const UEMUSHORT TFbignan[8] =
6212 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6213 static const UEMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
6221 static const UEMUSHORT XFbignan[6] =
6222 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6223 static const UEMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
6231 static const UEMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
6232 static const UEMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
6240 static const UEMUSHORT SFbignan[2] = {0x7fff, 0xffff};
6241 static const UEMUSHORT SFlittlenan[2] = {0, 0xffc0};
6248 make_nan (nan, sign, mode)
6251 enum machine_mode mode;
6257 size = GET_MODE_BITSIZE (mode);
6258 if (LARGEST_EXPONENT_IS_NORMAL (size))
6260 warning ("%d-bit floats cannot hold NaNs", size);
6261 saturate (nan, sign, size, 0);
6266 /* Possibly the `reserved operand' patterns on a VAX can be
6267 used like NaN's, but probably not in the same way as IEEE. */
6268 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
6270 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6272 if (REAL_WORDS_BIG_ENDIAN)
6282 if (REAL_WORDS_BIG_ENDIAN)
6290 if (REAL_WORDS_BIG_ENDIAN)
6299 if (REAL_WORDS_BIG_ENDIAN)
6309 if (REAL_WORDS_BIG_ENDIAN)
6310 *nan++ = (sign << 15) | (*p++ & 0x7fff);
6313 if (! REAL_WORDS_BIG_ENDIAN)
6314 *nan = (sign << 15) | (*p & 0x7fff);
6319 /* Create a saturation value for a SIZE-bit float, assuming that
6320 LARGEST_EXPONENT_IS_NORMAL (SIZE).
6322 If SIGN is true, fill X with the most negative value, otherwise fill
6323 it with the most positive value. WARN is true if the function should
6324 warn about overflow. */
6327 saturate (x, sign, size, warn)
6329 int sign, size, warn;
6333 if (warn && extra_warnings)
6334 warning ("value exceeds the range of a %d-bit float", size);
6336 /* Create the most negative value. */
6337 for (i = 0; i < size / EMUSHORT_SIZE; i++)
6340 /* Make it positive, if necessary. */
6342 x[REAL_WORDS_BIG_ENDIAN? 0 : i - 1] = 0x7fff;
6346 /* This is the inverse of the function `etarsingle' invoked by
6347 REAL_VALUE_TO_TARGET_SINGLE. */
6350 ereal_unto_float (f)
6357 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6358 This is the inverse operation to what the function `endian' does. */
6359 if (REAL_WORDS_BIG_ENDIAN)
6361 s[0] = (UEMUSHORT) (f >> 16);
6362 s[1] = (UEMUSHORT) f;
6366 s[0] = (UEMUSHORT) f;
6367 s[1] = (UEMUSHORT) (f >> 16);
6369 /* Convert and promote the target float to E-type. */
6371 /* Output E-type to REAL_VALUE_TYPE. */
6377 /* This is the inverse of the function `etardouble' invoked by
6378 REAL_VALUE_TO_TARGET_DOUBLE. */
6381 ereal_unto_double (d)
6388 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
6389 if (REAL_WORDS_BIG_ENDIAN)
6391 s[0] = (UEMUSHORT) (d[0] >> 16);
6392 s[1] = (UEMUSHORT) d[0];
6393 s[2] = (UEMUSHORT) (d[1] >> 16);
6394 s[3] = (UEMUSHORT) d[1];
6398 /* Target float words are little-endian. */
6399 s[0] = (UEMUSHORT) d[0];
6400 s[1] = (UEMUSHORT) (d[0] >> 16);
6401 s[2] = (UEMUSHORT) d[1];
6402 s[3] = (UEMUSHORT) (d[1] >> 16);
6404 /* Convert target double to E-type. */
6406 /* Output E-type to REAL_VALUE_TYPE. */
6412 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
6413 This is somewhat like ereal_unto_float, but the input types
6414 for these are different. */
6417 ereal_from_float (f)
6424 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6425 This is the inverse operation to what the function `endian' does. */
6426 if (REAL_WORDS_BIG_ENDIAN)
6428 s[0] = (UEMUSHORT) (f >> 16);
6429 s[1] = (UEMUSHORT) f;
6433 s[0] = (UEMUSHORT) f;
6434 s[1] = (UEMUSHORT) (f >> 16);
6436 /* Convert and promote the target float to E-type. */
6438 /* Output E-type to REAL_VALUE_TYPE. */
6444 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
6445 This is somewhat like ereal_unto_double, but the input types
6446 for these are different.
6448 The DFmode is stored as an array of HOST_WIDE_INT in the target's
6449 data format, with no holes in the bit packing. The first element
6450 of the input array holds the bits that would come first in the
6451 target computer's memory. */
6454 ereal_from_double (d)
6461 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
6462 if (REAL_WORDS_BIG_ENDIAN)
6464 #if HOST_BITS_PER_WIDE_INT == 32
6465 s[0] = (UEMUSHORT) (d[0] >> 16);
6466 s[1] = (UEMUSHORT) d[0];
6467 s[2] = (UEMUSHORT) (d[1] >> 16);
6468 s[3] = (UEMUSHORT) d[1];
6470 /* In this case the entire target double is contained in the
6471 first array element. The second element of the input is
6473 s[0] = (UEMUSHORT) (d[0] >> 48);
6474 s[1] = (UEMUSHORT) (d[0] >> 32);
6475 s[2] = (UEMUSHORT) (d[0] >> 16);
6476 s[3] = (UEMUSHORT) d[0];
6481 /* Target float words are little-endian. */
6482 s[0] = (UEMUSHORT) d[0];
6483 s[1] = (UEMUSHORT) (d[0] >> 16);
6484 #if HOST_BITS_PER_WIDE_INT == 32
6485 s[2] = (UEMUSHORT) d[1];
6486 s[3] = (UEMUSHORT) (d[1] >> 16);
6488 s[2] = (UEMUSHORT) (d[0] >> 32);
6489 s[3] = (UEMUSHORT) (d[0] >> 48);
6492 /* Convert target double to E-type. */
6494 /* Output E-type to REAL_VALUE_TYPE. */
6501 /* Convert target computer unsigned 64-bit integer to e-type.
6502 The endian-ness of DImode follows the convention for integers,
6503 so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN. */
6507 const UEMUSHORT *di; /* Address of the 64-bit int. */
6514 if (WORDS_BIG_ENDIAN)
6516 for (k = M; k < M + 4; k++)
6521 for (k = M + 3; k >= M; k--)
6524 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
6525 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6526 ecleaz (yi); /* it was zero */
6528 yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
6532 /* Convert target computer signed 64-bit integer to e-type. */
6536 const UEMUSHORT *di; /* Address of the 64-bit int. */
6539 unsigned EMULONG acc;
6545 if (WORDS_BIG_ENDIAN)
6547 for (k = M; k < M + 4; k++)
6552 for (k = M + 3; k >= M; k--)
6555 /* Take absolute value */
6561 for (k = M + 3; k >= M; k--)
6563 acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
6570 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
6571 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6572 ecleaz (yi); /* it was zero */
6574 yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
6581 /* Convert e-type to unsigned 64-bit int. */
6597 k = (int) xi[E] - (EXONE - 1);
6600 for (j = 0; j < 4; j++)
6606 for (j = 0; j < 4; j++)
6609 warning ("overflow on truncation to integer");
6614 /* Shift more than 16 bits: first shift up k-16 mod 16,
6615 then shift up by 16's. */
6616 j = k - ((k >> 4) << 4);
6620 if (WORDS_BIG_ENDIAN)
6631 if (WORDS_BIG_ENDIAN)
6636 while ((k -= 16) > 0);
6640 /* shift not more than 16 bits */
6645 if (WORDS_BIG_ENDIAN)
6664 /* Convert e-type to signed 64-bit int. */
6671 unsigned EMULONG acc;
6678 k = (int) xi[E] - (EXONE - 1);
6681 for (j = 0; j < 4; j++)
6687 for (j = 0; j < 4; j++)
6690 warning ("overflow on truncation to integer");
6696 /* Shift more than 16 bits: first shift up k-16 mod 16,
6697 then shift up by 16's. */
6698 j = k - ((k >> 4) << 4);
6702 if (WORDS_BIG_ENDIAN)
6713 if (WORDS_BIG_ENDIAN)
6718 while ((k -= 16) > 0);
6722 /* shift not more than 16 bits */
6725 if (WORDS_BIG_ENDIAN)
6741 /* Negate if negative */
6745 if (WORDS_BIG_ENDIAN)
6747 for (k = 0; k < 4; k++)
6749 acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
6750 if (WORDS_BIG_ENDIAN)
6762 /* Longhand square root routine. */
6765 static int esqinited = 0;
6766 static unsigned short sqrndbit[NI];
6773 UEMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
6775 int i, j, k, n, nlups;
6780 sqrndbit[NI - 2] = 1;
6783 /* Check for arg <= 0 */
6784 i = ecmp (x, ezero);
6789 mtherr ("esqrt", DOMAIN);
6805 /* Bring in the arg and renormalize if it is denormal. */
6807 m = (EMULONG) xx[1]; /* local long word exponent */
6811 /* Divide exponent by 2 */
6813 exp = (unsigned short) ((m / 2) + 0x3ffe);
6815 /* Adjust if exponent odd */
6825 n = 8; /* get 8 bits of result per inner loop */
6831 /* bring in next word of arg */
6833 num[NI - 1] = xx[j + 3];
6834 /* Do additional bit on last outer loop, for roundoff. */
6837 for (i = 0; i < n; i++)
6839 /* Next 2 bits of arg */
6842 /* Shift up answer */
6844 /* Make trial divisor */
6845 for (k = 0; k < NI; k++)
6848 eaddm (sqrndbit, temp);
6849 /* Subtract and insert answer bit if it goes in */
6850 if (ecmpm (temp, num) <= 0)
6860 /* Adjust for extra, roundoff loop done. */
6861 exp += (NBITS - 1) - rndprc;
6863 /* Sticky bit = 1 if the remainder is nonzero. */
6865 for (i = 3; i < NI; i++)
6868 /* Renormalize and round off. */
6869 emdnorm (sq, k, 0, exp, !ROUND_TOWARDS_ZERO);
6874 /* Return the binary precision of the significand for a given
6875 floating point mode. The mode can hold an integer value
6876 that many bits wide, without losing any bits. */
6879 significand_size (mode)
6880 enum machine_mode mode;
6883 /* Don't test the modes, but their sizes, lest this
6884 code won't work for BITS_PER_UNIT != 8 . */
6886 switch (GET_MODE_BITSIZE (mode))
6890 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
6897 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
6900 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
6903 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
6906 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
6919 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)