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 /* You will have to modify this program to have a smaller unit size. */
169 #define EMU_NON_COMPILE
175 /* If no 16-bit type has been found and the compiler is GCC, try HImode. */
176 #if defined(__GNUC__) && EMUSHORT_SIZE != 16
177 typedef int HItype __attribute__ ((mode (HI)));
178 typedef unsigned int UHItype __attribute__ ((mode (HI)));
182 #define EMUSHORT HItype
183 #define UEMUSHORT UHItype
184 #define EMUSHORT_SIZE 16
185 #define EMULONG_SIZE 32
187 #define UEMUSHORT unsigned EMUSHORT
190 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
191 #define EMULONG short
193 #if HOST_BITS_PER_INT >= EMULONG_SIZE
196 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
199 #if HOST_BITS_PER_LONGLONG >= EMULONG_SIZE
200 #define EMULONG long long int
202 /* You will have to modify this program to have a smaller unit size. */
203 #define EMU_NON_COMPILE
210 /* The host interface doesn't work if no 16-bit size exists. */
211 #if EMUSHORT_SIZE != 16
212 #define EMU_NON_COMPILE
215 /* OK to continue compilation. */
216 #ifndef EMU_NON_COMPILE
218 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
219 In GET_REAL and PUT_REAL, r and e are pointers.
220 A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
221 in memory, with no holes. */
223 #if MAX_LONG_DOUBLE_TYPE_SIZE == 96 || \
224 ((INTEL_EXTENDED_IEEE_FORMAT != 0) && MAX_LONG_DOUBLE_TYPE_SIZE == 128)
225 /* Number of 16 bit words in external e type format */
227 # define MAXDECEXP 4932
228 # define MINDECEXP -4956
229 # define GET_REAL(r,e) memcpy ((e), (r), 2*NE)
230 # define PUT_REAL(e,r) \
232 memcpy ((r), (e), 2*NE); \
233 if (2*NE < sizeof (*r)) \
234 memset ((char *) (r) + 2*NE, 0, sizeof (*r) - 2*NE); \
236 # else /* no XFmode */
237 # if MAX_LONG_DOUBLE_TYPE_SIZE == 128
239 # define MAXDECEXP 4932
240 # define MINDECEXP -4977
241 # define GET_REAL(r,e) memcpy ((e), (r), 2*NE)
242 # define PUT_REAL(e,r) \
244 memcpy ((r), (e), 2*NE); \
245 if (2*NE < sizeof (*r)) \
246 memset ((char *) (r) + 2*NE, 0, sizeof (*r) - 2*NE); \
250 #define MAXDECEXP 4932
251 #define MINDECEXP -4956
252 /* Emulator uses target format internally
253 but host stores it in host endian-ness. */
255 #define GET_REAL(r,e) \
257 if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
258 e53toe ((const UEMUSHORT *) (r), (e)); \
262 memcpy (&w[3], ((const EMUSHORT *) r), sizeof (EMUSHORT)); \
263 memcpy (&w[2], ((const EMUSHORT *) r) + 1, sizeof (EMUSHORT)); \
264 memcpy (&w[1], ((const EMUSHORT *) r) + 2, sizeof (EMUSHORT)); \
265 memcpy (&w[0], ((const EMUSHORT *) r) + 3, sizeof (EMUSHORT)); \
270 #define PUT_REAL(e,r) \
272 if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
273 etoe53 ((e), (UEMUSHORT *) (r)); \
278 memcpy (((EMUSHORT *) r), &w[3], sizeof (EMUSHORT)); \
279 memcpy (((EMUSHORT *) r) + 1, &w[2], sizeof (EMUSHORT)); \
280 memcpy (((EMUSHORT *) r) + 2, &w[1], sizeof (EMUSHORT)); \
281 memcpy (((EMUSHORT *) r) + 3, &w[0], sizeof (EMUSHORT)); \
285 #endif /* not TFmode */
286 #endif /* not XFmode */
289 /* Number of 16 bit words in internal format */
292 /* Array offset to exponent */
295 /* Array offset to high guard word */
298 /* Number of bits of precision */
299 #define NBITS ((NI-4)*16)
301 /* Maximum number of decimal digits in ASCII conversion
304 #define NDEC (NBITS*8/27)
306 /* The exponent of 1.0 */
307 #define EXONE (0x3fff)
309 #if defined(HOST_EBCDIC)
310 /* bit 8 is significant in EBCDIC */
311 #define CHARMASK 0xff
313 #define CHARMASK 0x7f
316 extern int extra_warnings;
317 extern const UEMUSHORT ezero[NE], ehalf[NE], eone[NE], etwo[NE];
318 extern const UEMUSHORT elog2[NE], esqrt2[NE];
320 static void endian PARAMS ((const UEMUSHORT *, long *,
322 static void eclear PARAMS ((UEMUSHORT *));
323 static void emov PARAMS ((const UEMUSHORT *, UEMUSHORT *));
325 static void eabs PARAMS ((UEMUSHORT *));
327 static void eneg PARAMS ((UEMUSHORT *));
328 static int eisneg PARAMS ((const UEMUSHORT *));
329 static int eisinf PARAMS ((const UEMUSHORT *));
330 static int eisnan PARAMS ((const UEMUSHORT *));
331 static void einfin PARAMS ((UEMUSHORT *));
333 static void enan PARAMS ((UEMUSHORT *, int));
334 static void einan PARAMS ((UEMUSHORT *));
335 static int eiisnan PARAMS ((const UEMUSHORT *));
336 static int eiisneg PARAMS ((const UEMUSHORT *));
337 static void make_nan PARAMS ((UEMUSHORT *, int, enum machine_mode));
339 static void saturate PARAMS ((UEMUSHORT *, int, int, int));
340 static void emovi PARAMS ((const UEMUSHORT *, UEMUSHORT *));
341 static void emovo PARAMS ((const UEMUSHORT *, UEMUSHORT *));
342 static void ecleaz PARAMS ((UEMUSHORT *));
343 static void ecleazs PARAMS ((UEMUSHORT *));
344 static void emovz PARAMS ((const UEMUSHORT *, UEMUSHORT *));
346 static void eiinfin PARAMS ((UEMUSHORT *));
349 static int eiisinf PARAMS ((const UEMUSHORT *));
351 static int ecmpm PARAMS ((const UEMUSHORT *, const UEMUSHORT *));
352 static void eshdn1 PARAMS ((UEMUSHORT *));
353 static void eshup1 PARAMS ((UEMUSHORT *));
354 static void eshdn8 PARAMS ((UEMUSHORT *));
355 static void eshup8 PARAMS ((UEMUSHORT *));
356 static void eshup6 PARAMS ((UEMUSHORT *));
357 static void eshdn6 PARAMS ((UEMUSHORT *));
358 static void eaddm PARAMS ((const UEMUSHORT *, UEMUSHORT *));
\f
359 static void esubm PARAMS ((const UEMUSHORT *, UEMUSHORT *));
360 static void m16m PARAMS ((unsigned int, const UEMUSHORT *, UEMUSHORT *));
361 static int edivm PARAMS ((const UEMUSHORT *, UEMUSHORT *));
362 static int emulm PARAMS ((const UEMUSHORT *, UEMUSHORT *));
363 static void emdnorm PARAMS ((UEMUSHORT *, int, int, EMULONG, int));
364 static void esub PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
366 static void eadd PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
368 static void eadd1 PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
370 static void ediv PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
372 static void emul PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
374 static void e53toe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
375 static void e64toe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
376 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
377 static void e113toe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
379 static void e24toe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
380 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
381 static void etoe113 PARAMS ((const UEMUSHORT *, UEMUSHORT *));
382 static void toe113 PARAMS ((UEMUSHORT *, UEMUSHORT *));
384 static void etoe64 PARAMS ((const UEMUSHORT *, UEMUSHORT *));
385 static void toe64 PARAMS ((UEMUSHORT *, UEMUSHORT *));
386 static void etoe53 PARAMS ((const UEMUSHORT *, UEMUSHORT *));
387 static void toe53 PARAMS ((UEMUSHORT *, UEMUSHORT *));
388 static void etoe24 PARAMS ((const UEMUSHORT *, UEMUSHORT *));
389 static void toe24 PARAMS ((UEMUSHORT *, UEMUSHORT *));
390 static int ecmp PARAMS ((const UEMUSHORT *, const UEMUSHORT *));
392 static void eround PARAMS ((const UEMUSHORT *, UEMUSHORT *));
394 static void ltoe PARAMS ((const HOST_WIDE_INT *, UEMUSHORT *));
395 static void ultoe PARAMS ((const unsigned HOST_WIDE_INT *, UEMUSHORT *));
396 static void eifrac PARAMS ((const UEMUSHORT *, HOST_WIDE_INT *,
398 static void euifrac PARAMS ((const UEMUSHORT *, unsigned HOST_WIDE_INT *,
400 static int eshift PARAMS ((UEMUSHORT *, int));
401 static int enormlz PARAMS ((UEMUSHORT *));
403 static void e24toasc PARAMS ((const UEMUSHORT *, char *, int));
404 static void e53toasc PARAMS ((const UEMUSHORT *, char *, int));
405 static void e64toasc PARAMS ((const UEMUSHORT *, char *, int));
406 static void e113toasc PARAMS ((const UEMUSHORT *, char *, int));
408 static void etoasc PARAMS ((const UEMUSHORT *, char *, int));
409 static void asctoe24 PARAMS ((const char *, UEMUSHORT *));
410 static void asctoe53 PARAMS ((const char *, UEMUSHORT *));
411 static void asctoe64 PARAMS ((const char *, UEMUSHORT *));
412 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
413 static void asctoe113 PARAMS ((const char *, UEMUSHORT *));
415 static void asctoe PARAMS ((const char *, UEMUSHORT *));
416 static void asctoeg PARAMS ((const char *, UEMUSHORT *, int));
417 static void efloor PARAMS ((const UEMUSHORT *, UEMUSHORT *));
419 static void efrexp PARAMS ((const UEMUSHORT *, int *,
422 static void eldexp PARAMS ((const UEMUSHORT *, int, UEMUSHORT *));
424 static void eremain PARAMS ((const UEMUSHORT *, const UEMUSHORT *,
427 static void eiremain PARAMS ((UEMUSHORT *, UEMUSHORT *));
428 static void mtherr PARAMS ((const char *, int));
430 static void dectoe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
431 static void etodec PARAMS ((const UEMUSHORT *, UEMUSHORT *));
432 static void todec PARAMS ((UEMUSHORT *, UEMUSHORT *));
435 static void ibmtoe PARAMS ((const UEMUSHORT *, UEMUSHORT *,
437 static void etoibm PARAMS ((const UEMUSHORT *, UEMUSHORT *,
439 static void toibm PARAMS ((UEMUSHORT *, UEMUSHORT *,
443 static void c4xtoe PARAMS ((const UEMUSHORT *, UEMUSHORT *,
445 static void etoc4x PARAMS ((const UEMUSHORT *, UEMUSHORT *,
447 static void toc4x PARAMS ((UEMUSHORT *, UEMUSHORT *,
451 static void uditoe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
452 static void ditoe PARAMS ((const UEMUSHORT *, UEMUSHORT *));
453 static void etoudi PARAMS ((const UEMUSHORT *, UEMUSHORT *));
454 static void etodi PARAMS ((const UEMUSHORT *, UEMUSHORT *));
455 static void esqrt PARAMS ((const UEMUSHORT *, UEMUSHORT *));
458 /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
459 swapping ends if required, into output array of longs. The
460 result is normally passed to fprintf by the ASM_OUTPUT_ macros. */
466 enum machine_mode mode;
470 if (REAL_WORDS_BIG_ENDIAN)
475 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
476 /* Swap halfwords in the fourth long. */
477 th = (unsigned long) e[6] & 0xffff;
478 t = (unsigned long) e[7] & 0xffff;
487 /* Swap halfwords in the third long. */
488 th = (unsigned long) e[4] & 0xffff;
489 t = (unsigned long) e[5] & 0xffff;
495 /* Swap halfwords in the second word. */
496 th = (unsigned long) e[2] & 0xffff;
497 t = (unsigned long) e[3] & 0xffff;
504 /* Swap halfwords in the first word. */
505 th = (unsigned long) e[0] & 0xffff;
506 t = (unsigned long) e[1] & 0xffff;
517 /* Pack the output array without swapping. */
522 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
523 /* Pack the fourth long. */
524 th = (unsigned long) e[7] & 0xffff;
525 t = (unsigned long) e[6] & 0xffff;
534 /* Pack the third long.
535 Each element of the input REAL_VALUE_TYPE array has 16 useful bits
537 th = (unsigned long) e[5] & 0xffff;
538 t = (unsigned long) e[4] & 0xffff;
544 /* Pack the second long */
545 th = (unsigned long) e[3] & 0xffff;
546 t = (unsigned long) e[2] & 0xffff;
553 /* Pack the first long */
554 th = (unsigned long) e[1] & 0xffff;
555 t = (unsigned long) e[0] & 0xffff;
567 /* This is the implementation of the REAL_ARITHMETIC macro. */
570 earith (value, icode, r1, r2)
571 REAL_VALUE_TYPE *value;
576 UEMUSHORT d1[NE], d2[NE], v[NE];
582 /* Return NaN input back to the caller. */
585 PUT_REAL (d1, value);
590 PUT_REAL (d2, value);
594 code = (enum tree_code) icode;
602 esub (d2, d1, v); /* d1 - d2 */
610 #ifndef REAL_INFINITY
611 if (ecmp (d2, ezero) == 0)
614 enan (v, eisneg (d1) ^ eisneg (d2));
621 ediv (d2, d1, v); /* d1/d2 */
624 case MIN_EXPR: /* min (d1,d2) */
625 if (ecmp (d1, d2) < 0)
631 case MAX_EXPR: /* max (d1,d2) */
632 if (ecmp (d1, d2) > 0)
645 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
646 implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */
652 UEMUSHORT f[NE], g[NE];
668 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
669 implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */
675 UEMUSHORT f[NE], g[NE];
677 unsigned HOST_WIDE_INT l;
691 /* This is the REAL_VALUE_ATOF function. It converts a decimal or hexadecimal
692 string to binary, rounding off as indicated by the machine_mode argument.
693 Then it promotes the rounded value to REAL_VALUE_TYPE. */
700 UEMUSHORT tem[NE], e[NE];
726 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
746 /* Expansion of REAL_NEGATE. */
762 /* Round real toward zero to HOST_WIDE_INT;
763 implements REAL_VALUE_FIX (x). */
769 UEMUSHORT f[NE], g[NE];
776 warning ("conversion from NaN to int");
784 /* Round real toward zero to unsigned HOST_WIDE_INT
785 implements REAL_VALUE_UNSIGNED_FIX (x).
786 Negative input returns zero. */
788 unsigned HOST_WIDE_INT
792 UEMUSHORT f[NE], g[NE];
793 unsigned HOST_WIDE_INT l;
799 warning ("conversion from NaN to unsigned int");
808 /* REAL_VALUE_FROM_INT macro. */
811 ereal_from_int (d, i, j, mode)
814 enum machine_mode mode;
816 UEMUSHORT df[NE], dg[NE];
817 HOST_WIDE_INT low, high;
820 if (GET_MODE_CLASS (mode) != MODE_FLOAT)
827 /* complement and add 1 */
834 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
835 ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
837 ultoe ((unsigned HOST_WIDE_INT *) &low, df);
842 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
843 Avoid double-rounding errors later by rounding off now from the
844 extra-wide internal format to the requested precision. */
845 switch (GET_MODE_BITSIZE (mode))
863 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
880 /* REAL_VALUE_FROM_UNSIGNED_INT macro. */
883 ereal_from_uint (d, i, j, mode)
885 unsigned HOST_WIDE_INT i, j;
886 enum machine_mode mode;
888 UEMUSHORT df[NE], dg[NE];
889 unsigned HOST_WIDE_INT low, high;
891 if (GET_MODE_CLASS (mode) != MODE_FLOAT)
895 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
901 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
902 Avoid double-rounding errors later by rounding off now from the
903 extra-wide internal format to the requested precision. */
904 switch (GET_MODE_BITSIZE (mode))
922 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
939 /* REAL_VALUE_TO_INT macro. */
942 ereal_to_int (low, high, rr)
943 HOST_WIDE_INT *low, *high;
946 UEMUSHORT d[NE], df[NE], dg[NE], dh[NE];
953 warning ("conversion from NaN to int");
959 /* convert positive value */
966 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
967 ediv (df, d, dg); /* dg = d / 2^32 is the high word */
968 euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
969 emul (df, dh, dg); /* fractional part is the low word */
970 euifrac (dg, (unsigned HOST_WIDE_INT *) low, dh);
973 /* complement and add 1 */
983 /* REAL_VALUE_LDEXP macro. */
990 UEMUSHORT e[NE], y[NE];
1003 /* Check for infinity in a REAL_VALUE_TYPE. */
1007 REAL_VALUE_TYPE x ATTRIBUTE_UNUSED;
1013 return (eisinf (e));
1019 /* Check whether a REAL_VALUE_TYPE item is a NaN. */
1023 REAL_VALUE_TYPE x ATTRIBUTE_UNUSED;
1029 return (eisnan (e));
1036 /* Check for a negative REAL_VALUE_TYPE number.
1037 This just checks the sign bit, so that -0 counts as negative. */
1043 return ereal_isneg (x);
1046 /* Expansion of REAL_VALUE_TRUNCATE.
1047 The result is in floating point, rounded to nearest or even. */
1050 real_value_truncate (mode, arg)
1051 enum machine_mode mode;
1052 REAL_VALUE_TYPE arg;
1054 UEMUSHORT e[NE], t[NE];
1066 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
1103 /* If an unsupported type was requested, presume that
1104 the machine files know something useful to do with
1105 the unmodified value. */
1114 /* Try to change R into its exact multiplicative inverse in machine mode
1115 MODE. Return nonzero function value if successful. */
1118 exact_real_inverse (mode, r)
1119 enum machine_mode mode;
1122 UEMUSHORT e[NE], einv[NE];
1123 REAL_VALUE_TYPE rinv;
1128 /* Test for input in range. Don't transform IEEE special values. */
1129 if (eisinf (e) || eisnan (e) || (ecmp (e, ezero) == 0))
1132 /* Test for a power of 2: all significand bits zero except the MSB.
1133 We are assuming the target has binary (or hex) arithmetic. */
1134 if (e[NE - 2] != 0x8000)
1137 for (i = 0; i < NE - 2; i++)
1143 /* Compute the inverse and truncate it to the required mode. */
1144 ediv (e, eone, einv);
1145 PUT_REAL (einv, &rinv);
1146 rinv = real_value_truncate (mode, rinv);
1148 #ifdef CHECK_FLOAT_VALUE
1149 /* This check is not redundant. It may, for example, flush
1150 a supposedly IEEE denormal value to zero. */
1152 if (CHECK_FLOAT_VALUE (mode, rinv, i))
1155 GET_REAL (&rinv, einv);
1157 /* Check the bits again, because the truncation might have
1158 generated an arbitrary saturation value on overflow. */
1159 if (einv[NE - 2] != 0x8000)
1162 for (i = 0; i < NE - 2; i++)
1168 /* Fail if the computed inverse is out of range. */
1169 if (eisinf (einv) || eisnan (einv) || (ecmp (einv, ezero) == 0))
1172 /* Output the reciprocal and return success flag. */
1177 /* Used for debugging--print the value of R in human-readable format
1186 REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
1187 fprintf (stderr, "%s", dstr);
1191 /* The following routines convert REAL_VALUE_TYPE to the various floating
1192 point formats that are meaningful to supported computers.
1194 The results are returned in 32-bit pieces, each piece stored in a `long'.
1195 This is so they can be printed by statements like
1197 fprintf (file, "%lx, %lx", L[0], L[1]);
1199 that will work on both narrow- and wide-word host computers. */
1201 /* Convert R to a 128-bit long double precision value. The output array L
1202 contains four 32-bit pieces of the result, in the order they would appear
1213 #if INTEL_EXTENDED_IEEE_FORMAT == 0
1218 endian (e, l, TFmode);
1221 /* Convert R to a double extended precision value. The output array L
1222 contains three 32-bit pieces of the result, in the order they would
1223 appear in memory. */
1234 endian (e, l, XFmode);
1237 /* Convert R to a double precision value. The output array L contains two
1238 32-bit pieces of the result, in the order they would appear in memory. */
1249 endian (e, l, DFmode);
1252 /* Convert R to a single precision float value stored in the least-significant
1253 bits of a `long'. */
1264 endian (e, &l, SFmode);
1268 /* Convert X to a decimal ASCII string S for output to an assembly
1269 language file. Note, there is no standard way to spell infinity or
1270 a NaN, so these values may require special treatment in the tm.h
1274 ereal_to_decimal (x, s)
1284 /* Compare X and Y. Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
1285 or -2 if either is a NaN. */
1289 REAL_VALUE_TYPE x, y;
1291 UEMUSHORT ex[NE], ey[NE];
1295 return (ecmp (ex, ey));
1298 /* Return 1 if the sign bit of X is set, else return 0. */
1307 return (eisneg (ex));
1312 Extended precision IEEE binary floating point arithmetic routines
1314 Numbers are stored in C language as arrays of 16-bit unsigned
1315 short integers. The arguments of the routines are pointers to
1318 External e type data structure, similar to Intel 8087 chip
1319 temporary real format but possibly with a larger significand:
1321 NE-1 significand words (least significant word first,
1322 most significant bit is normally set)
1323 exponent (value = EXONE for 1.0,
1324 top bit is the sign)
1327 Internal exploded e-type data structure of a number (a "word" is 16 bits):
1329 ei[0] sign word (0 for positive, 0xffff for negative)
1330 ei[1] biased exponent (value = EXONE for the number 1.0)
1331 ei[2] high guard word (always zero after normalization)
1333 to ei[NI-2] significand (NI-4 significand words,
1334 most significant word first,
1335 most significant bit is set)
1336 ei[NI-1] low guard word (0x8000 bit is rounding place)
1340 Routines for external format e-type numbers
1342 asctoe (string, e) ASCII string to extended double e type
1343 asctoe64 (string, &d) ASCII string to long double
1344 asctoe53 (string, &d) ASCII string to double
1345 asctoe24 (string, &f) ASCII string to single
1346 asctoeg (string, e, prec) ASCII string to specified precision
1347 e24toe (&f, e) IEEE single precision to e type
1348 e53toe (&d, e) IEEE double precision to e type
1349 e64toe (&d, e) IEEE long double precision to e type
1350 e113toe (&d, e) 128-bit long double precision to e type
1352 eabs (e) absolute value
1354 eadd (a, b, c) c = b + a
1356 ecmp (a, b) Returns 1 if a > b, 0 if a == b,
1357 -1 if a < b, -2 if either a or b is a NaN.
1358 ediv (a, b, c) c = b / a
1359 efloor (a, b) truncate to integer, toward -infinity
1360 efrexp (a, exp, s) extract exponent and significand
1361 eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
1362 euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
1363 einfin (e) set e to infinity, leaving its sign alone
1364 eldexp (a, n, b) multiply by 2**n
1366 emul (a, b, c) c = b * a
1369 eround (a, b) b = nearest integer value to a
1371 esub (a, b, c) c = b - a
1373 e24toasc (&f, str, n) single to ASCII string, n digits after decimal
1374 e53toasc (&d, str, n) double to ASCII string, n digits after decimal
1375 e64toasc (&d, str, n) 80-bit long double to ASCII string
1376 e113toasc (&d, str, n) 128-bit long double to ASCII string
1378 etoasc (e, str, n) e to ASCII string, n digits after decimal
1379 etoe24 (e, &f) convert e type to IEEE single precision
1380 etoe53 (e, &d) convert e type to IEEE double precision
1381 etoe64 (e, &d) convert e type to IEEE long double precision
1382 ltoe (&l, e) HOST_WIDE_INT to e type
1383 ultoe (&l, e) unsigned HOST_WIDE_INT to e type
1384 eisneg (e) 1 if sign bit of e != 0, else 0
1385 eisinf (e) 1 if e has maximum exponent (non-IEEE)
1386 or is infinite (IEEE)
1387 eisnan (e) 1 if e is a NaN
1390 Routines for internal format exploded e-type numbers
1392 eaddm (ai, bi) add significands, bi = bi + ai
1394 ecleazs (ei) set ei = 0 but leave its sign alone
1395 ecmpm (ai, bi) compare significands, return 1, 0, or -1
1396 edivm (ai, bi) divide significands, bi = bi / ai
1397 emdnorm (ai,l,s,exp) normalize and round off
1398 emovi (a, ai) convert external a to internal ai
1399 emovo (ai, a) convert internal ai to external a
1400 emovz (ai, bi) bi = ai, low guard word of bi = 0
1401 emulm (ai, bi) multiply significands, bi = bi * ai
1402 enormlz (ei) left-justify the significand
1403 eshdn1 (ai) shift significand and guards down 1 bit
1404 eshdn8 (ai) shift down 8 bits
1405 eshdn6 (ai) shift down 16 bits
1406 eshift (ai, n) shift ai n bits up (or down if n < 0)
1407 eshup1 (ai) shift significand and guards up 1 bit
1408 eshup8 (ai) shift up 8 bits
1409 eshup6 (ai) shift up 16 bits
1410 esubm (ai, bi) subtract significands, bi = bi - ai
1411 eiisinf (ai) 1 if infinite
1412 eiisnan (ai) 1 if a NaN
1413 eiisneg (ai) 1 if sign bit of ai != 0, else 0
1414 einan (ai) set ai = NaN
1416 eiinfin (ai) set ai = infinity
1419 The result is always normalized and rounded to NI-4 word precision
1420 after each arithmetic operation.
1422 Exception flags are NOT fully supported.
1424 Signaling NaN's are NOT supported; they are treated the same
1427 Define INFINITY for support of infinity; otherwise a
1428 saturation arithmetic is implemented.
1430 Define NANS for support of Not-a-Number items; otherwise the
1431 arithmetic will never produce a NaN output, and might be confused
1433 If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1434 either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1435 may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1438 Denormals are always supported here where appropriate (e.g., not
1439 for conversion to DEC numbers). */
1441 /* Definitions for error codes that are passed to the common error handling
1444 For Digital Equipment PDP-11 and VAX computers, certain
1445 IBM systems, and others that use numbers with a 56-bit
1446 significand, the symbol DEC should be defined. In this
1447 mode, most floating point constants are given as arrays
1448 of octal integers to eliminate decimal to binary conversion
1449 errors that might be introduced by the compiler.
1451 For computers, such as IBM PC, that follow the IEEE
1452 Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1453 Std 754-1985), the symbol IEEE should be defined.
1454 These numbers have 53-bit significands. In this mode, constants
1455 are provided as arrays of hexadecimal 16 bit integers.
1456 The endian-ness of generated values is controlled by
1457 REAL_WORDS_BIG_ENDIAN.
1459 To accommodate other types of computer arithmetic, all
1460 constants are also provided in a normal decimal radix
1461 which one can hope are correctly converted to a suitable
1462 format by the available C language compiler. To invoke
1463 this mode, the symbol UNK is defined.
1465 An important difference among these modes is a predefined
1466 set of machine arithmetic constants for each. The numbers
1467 MACHEP (the machine roundoff error), MAXNUM (largest number
1468 represented), and several other parameters are preset by
1469 the configuration symbol. Check the file const.c to
1470 ensure that these values are correct for your computer.
1472 For ANSI C compatibility, define ANSIC equal to 1. Currently
1473 this affects only the atan2 function and others that use it. */
1475 /* Constant definitions for math error conditions. */
1477 #define DOMAIN 1 /* argument domain error */
1478 #define SING 2 /* argument singularity */
1479 #define OVERFLOW 3 /* overflow range error */
1480 #define UNDERFLOW 4 /* underflow range error */
1481 #define TLOSS 5 /* total loss of precision */
1482 #define PLOSS 6 /* partial loss of precision */
1483 #define INVALID 7 /* NaN-producing operation */
1485 /* e type constants used by high precision check routines */
1487 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
1489 const UEMUSHORT ezero[NE] =
1490 {0x0000, 0x0000, 0x0000, 0x0000,
1491 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
1494 const UEMUSHORT ehalf[NE] =
1495 {0x0000, 0x0000, 0x0000, 0x0000,
1496 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
1499 const UEMUSHORT eone[NE] =
1500 {0x0000, 0x0000, 0x0000, 0x0000,
1501 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
1504 const UEMUSHORT etwo[NE] =
1505 {0x0000, 0x0000, 0x0000, 0x0000,
1506 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
1509 const UEMUSHORT e32[NE] =
1510 {0x0000, 0x0000, 0x0000, 0x0000,
1511 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
1513 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1514 const UEMUSHORT elog2[NE] =
1515 {0x40f3, 0xf6af, 0x03f2, 0xb398,
1516 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1518 /* 1.41421356237309504880168872420969807856967187537695E0 */
1519 const UEMUSHORT esqrt2[NE] =
1520 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1521 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1523 /* 3.14159265358979323846264338327950288419716939937511E0 */
1524 const UEMUSHORT epi[NE] =
1525 {0x2902, 0x1cd1, 0x80dc, 0x628b,
1526 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1529 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1530 const UEMUSHORT ezero[NE] =
1531 {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1532 const UEMUSHORT ehalf[NE] =
1533 {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1534 const UEMUSHORT eone[NE] =
1535 {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1536 const UEMUSHORT etwo[NE] =
1537 {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1538 const UEMUSHORT e32[NE] =
1539 {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1540 const UEMUSHORT elog2[NE] =
1541 {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1542 const UEMUSHORT esqrt2[NE] =
1543 {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1544 const UEMUSHORT epi[NE] =
1545 {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1548 /* Control register for rounding precision.
1549 This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits. */
1554 /* Clear out entire e-type number X. */
1562 for (i = 0; i < NE; i++)
1566 /* Move e-type number from A to B. */
1575 for (i = 0; i < NE; i++)
1581 /* Absolute value of e-type X. */
1587 /* sign is top bit of last word of external format */
1588 x[NE - 1] &= 0x7fff;
1592 /* Negate the e-type number X. */
1599 x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
1602 /* Return 1 if sign bit of e-type number X is nonzero, else zero. */
1606 const UEMUSHORT x[];
1609 if (x[NE - 1] & 0x8000)
1615 /* Return 1 if e-type number X is infinity, else return zero. */
1619 const UEMUSHORT x[];
1626 if ((x[NE - 1] & 0x7fff) == 0x7fff)
1632 /* Check if e-type number is not a number. The bit pattern is one that we
1633 defined, so we know for sure how to detect it. */
1637 const UEMUSHORT x[] ATTRIBUTE_UNUSED;
1642 /* NaN has maximum exponent */
1643 if ((x[NE - 1] & 0x7fff) != 0x7fff)
1645 /* ... and non-zero significand field. */
1646 for (i = 0; i < NE - 1; i++)
1656 /* Fill e-type number X with infinity pattern (IEEE)
1657 or largest possible number (non-IEEE). */
1666 for (i = 0; i < NE - 1; i++)
1670 for (i = 0; i < NE - 1; i++)
1698 /* Output an e-type NaN.
1699 This generates Intel's quiet NaN pattern for extended real.
1700 The exponent is 7fff, the leading mantissa word is c000. */
1710 for (i = 0; i < NE - 2; i++)
1713 *x = (sign << 15) | 0x7fff;
1717 /* Move in an e-type number A, converting it to exploded e-type B. */
1729 p = a + (NE - 1); /* point to last word of external number */
1730 /* get the sign bit */
1735 /* get the exponent */
1737 *q++ &= 0x7fff; /* delete the sign bit */
1739 if ((*(q - 1) & 0x7fff) == 0x7fff)
1745 for (i = 3; i < NI; i++)
1751 for (i = 2; i < NI; i++)
1757 /* clear high guard word */
1759 /* move in the significand */
1760 for (i = 0; i < NE - 1; i++)
1762 /* clear low guard word */
1766 /* Move out exploded e-type number A, converting it to e type B. */
1779 q = b + (NE - 1); /* point to output exponent */
1780 /* combine sign and exponent */
1783 *q-- = *p++ | 0x8000;
1787 if (*(p - 1) == 0x7fff)
1792 enan (b, eiisneg (a));
1800 /* skip over guard word */
1802 /* move the significand */
1803 for (j = 0; j < NE - 1; j++)
1807 /* Clear out exploded e-type number XI. */
1815 for (i = 0; i < NI; i++)
1819 /* Clear out exploded e-type XI, but don't touch the sign. */
1828 for (i = 0; i < NI - 1; i++)
1832 /* Move exploded e-type number from A to B. */
1841 for (i = 0; i < NI - 1; i++)
1843 /* clear low guard word */
1847 /* Generate exploded e-type NaN.
1848 The explicit pattern for this is maximum exponent and
1849 top two significant bits set. */
1863 /* Return nonzero if exploded e-type X is a NaN. */
1868 const UEMUSHORT x[];
1872 if ((x[E] & 0x7fff) == 0x7fff)
1874 for (i = M + 1; i < NI; i++)
1884 /* Return nonzero if sign of exploded e-type X is nonzero. */
1889 const UEMUSHORT x[];
1897 /* Fill exploded e-type X with infinity pattern.
1898 This has maximum exponent and significand all zeros. */
1910 /* Return nonzero if exploded e-type X is infinite. */
1915 const UEMUSHORT x[];
1922 if ((x[E] & 0x7fff) == 0x7fff)
1926 #endif /* INFINITY */
1928 /* Compare significands of numbers in internal exploded e-type format.
1929 Guard words are included in the comparison.
1937 const UEMUSHORT *a, *b;
1941 a += M; /* skip up to significand area */
1943 for (i = M; i < NI; i++)
1951 if (*(--a) > *(--b))
1957 /* Shift significand of exploded e-type X down by 1 bit. */
1966 x += M; /* point to significand area */
1969 for (i = M; i < NI; i++)
1981 /* Shift significand of exploded e-type X up by 1 bit. */
1993 for (i = M; i < NI; i++)
2006 /* Shift significand of exploded e-type X down by 8 bits. */
2012 UEMUSHORT newbyt, oldbyt;
2017 for (i = M; i < NI; i++)
2027 /* Shift significand of exploded e-type X up by 8 bits. */
2034 UEMUSHORT newbyt, oldbyt;
2039 for (i = M; i < NI; i++)
2049 /* Shift significand of exploded e-type X up by 16 bits. */
2061 for (i = M; i < NI - 1; i++)
2067 /* Shift significand of exploded e-type X down by 16 bits. */
2079 for (i = M; i < NI - 1; i++)
2085 /* Add significands of exploded e-type X and Y. X + Y replaces Y. */
2099 for (i = M; i < NI; i++)
2101 a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
2112 /* Subtract significands of exploded e-type X and Y. Y - X replaces Y. */
2126 for (i = M; i < NI; i++)
2128 a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
2140 static UEMUSHORT equot[NI];
2144 /* Radix 2 shift-and-add versions of multiply and divide */
2147 /* Divide significands */
2151 UEMUSHORT den[], num[];
2161 for (i = M; i < NI; i++)
2166 /* Use faster compare and subtraction if denominator has only 15 bits of
2172 for (i = M + 3; i < NI; i++)
2177 if ((den[M + 1] & 1) != 0)
2185 for (i = 0; i < NBITS + 2; i++)
2203 /* The number of quotient bits to calculate is NBITS + 1 scaling guard
2204 bit + 1 roundoff bit. */
2209 for (i = 0; i < NBITS + 2; i++)
2211 if (ecmpm (den, num) <= 0)
2214 j = 1; /* quotient bit = 1 */
2228 /* test for nonzero remainder after roundoff bit */
2231 for (i = M; i < NI; i++)
2239 for (i = 0; i < NI; i++)
2245 /* Multiply significands */
2256 for (i = M; i < NI; i++)
2261 while (*p == 0) /* significand is not supposed to be zero */
2266 if ((*p & 0xff) == 0)
2274 for (i = 0; i < k; i++)
2278 /* remember if there were any nonzero bits shifted out */
2285 for (i = 0; i < NI; i++)
2288 /* return flag for lost nonzero bits */
2294 /* Radix 65536 versions of multiply and divide. */
2296 /* Multiply significand of e-type number B
2297 by 16-bit quantity A, return e-type result to C. */
2302 const UEMUSHORT b[];
2306 unsigned EMULONG carry;
2307 const UEMUSHORT *ps;
2309 unsigned EMULONG aa, m;
2318 for (i=M+1; i<NI; i++)
2328 m = (unsigned EMULONG) aa * *ps--;
2329 carry = (m & 0xffff) + *pp;
2330 *pp-- = (UEMUSHORT) carry;
2331 carry = (carry >> 16) + (m >> 16) + *pp;
2332 *pp = (UEMUSHORT) carry;
2333 *(pp-1) = carry >> 16;
2336 for (i=M; i<NI; i++)
2340 /* Divide significands of exploded e-types NUM / DEN. Neither the
2341 numerator NUM nor the denominator DEN is permitted to have its high guard
2346 const UEMUSHORT den[];
2351 unsigned EMULONG tnum;
2352 UEMUSHORT j, tdenm, tquot;
2353 UEMUSHORT tprod[NI+1];
2359 for (i=M; i<NI; i++)
2365 for (i=M; i<NI; i++)
2367 /* Find trial quotient digit (the radix is 65536). */
2368 tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
2370 /* Do not execute the divide instruction if it will overflow. */
2371 if ((tdenm * (unsigned long) 0xffff) < tnum)
2374 tquot = tnum / tdenm;
2375 /* Multiply denominator by trial quotient digit. */
2376 m16m ((unsigned int) tquot, den, tprod);
2377 /* The quotient digit may have been overestimated. */
2378 if (ecmpm (tprod, num) > 0)
2382 if (ecmpm (tprod, num) > 0)
2392 /* test for nonzero remainder after roundoff bit */
2395 for (i=M; i<NI; i++)
2402 for (i=0; i<NI; i++)
2408 /* Multiply significands of exploded e-type A and B, result in B. */
2412 const UEMUSHORT a[];
2417 UEMUSHORT pprod[NI];
2423 for (i=M; i<NI; i++)
2429 for (i=M+1; i<NI; i++)
2437 m16m ((unsigned int) *p--, b, pprod);
2438 eaddm (pprod, equot);
2444 for (i=0; i<NI; i++)
2447 /* return flag for lost nonzero bits */
2453 /* Normalize and round off.
2455 The internal format number to be rounded is S.
2456 Input LOST is 0 if the value is exact. This is the so-called sticky bit.
2458 Input SUBFLG indicates whether the number was obtained
2459 by a subtraction operation. In that case if LOST is nonzero
2460 then the number is slightly smaller than indicated.
2462 Input EXP is the biased exponent, which may be negative.
2463 the exponent field of S is ignored but is replaced by
2464 EXP as adjusted by normalization and rounding.
2466 Input RCNTRL is the rounding control. If it is nonzero, the
2467 returned value will be rounded to RNDPRC bits.
2469 For future reference: In order for emdnorm to round off denormal
2470 significands at the right point, the input exponent must be
2471 adjusted to be the actual value it would have after conversion to
2472 the final floating point type. This adjustment has been
2473 implemented for all type conversions (etoe53, etc.) and decimal
2474 conversions, but not for the arithmetic functions (eadd, etc.).
2475 Data types having standard 15-bit exponents are not affected by
2476 this, but SFmode and DFmode are affected. For example, ediv with
2477 rndprc = 24 will not round correctly to 24-bit precision if the
2478 result is denormal. */
2480 static int rlast = -1;
2482 static UEMUSHORT rmsk = 0;
2483 static UEMUSHORT rmbit = 0;
2484 static UEMUSHORT rebit = 0;
2486 static UEMUSHORT rbit[NI];
2489 emdnorm (s, lost, subflg, exp, rcntrl)
2502 /* a blank significand could mean either zero or infinity. */
2515 if ((j > NBITS) && (exp < 32767))
2523 if (exp > (EMULONG) (-NBITS - 1))
2536 /* Round off, unless told not to by rcntrl. */
2539 /* Set up rounding parameters if the control register changed. */
2540 if (rndprc != rlast)
2547 rw = NI - 1; /* low guard word */
2570 /* For DEC or IBM arithmetic */
2587 /* For C4x arithmetic */
2608 /* Shift down 1 temporarily if the data structure has an implied
2609 most significant bit and the number is denormal.
2610 Intel long double denormals also lose one bit of precision. */
2611 if ((exp <= 0) && (rndprc != NBITS)
2612 && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2614 lost |= s[NI - 1] & 1;
2617 /* Clear out all bits below the rounding bit,
2618 remembering in r if any were nonzero. */
2632 if ((r & rmbit) != 0)
2638 { /* round to even */
2639 if ((s[re] & rebit) == 0)
2652 /* Undo the temporary shift for denormal values. */
2653 if ((exp <= 0) && (rndprc != NBITS)
2654 && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2659 { /* overflow on roundoff */
2672 for (i = 2; i < NI - 1; i++)
2675 warning ("floating point overflow");
2679 for (i = M + 1; i < NI - 1; i++)
2682 if ((rndprc < 64) || (rndprc == 113))
2697 s[1] = (UEMUSHORT) exp;
2700 /* Subtract. C = B - A, all e type numbers. */
2702 static int subflg = 0;
2706 const UEMUSHORT *a, *b;
2721 /* Infinity minus infinity is a NaN.
2722 Test for subtracting infinities of the same sign. */
2723 if (eisinf (a) && eisinf (b)
2724 && ((eisneg (a) ^ eisneg (b)) == 0))
2726 mtherr ("esub", INVALID);
2735 /* Add. C = A + B, all e type. */
2739 const UEMUSHORT *a, *b;
2744 /* NaN plus anything is a NaN. */
2755 /* Infinity minus infinity is a NaN.
2756 Test for adding infinities of opposite signs. */
2757 if (eisinf (a) && eisinf (b)
2758 && ((eisneg (a) ^ eisneg (b)) != 0))
2760 mtherr ("esub", INVALID);
2769 /* Arithmetic common to both addition and subtraction. */
2773 const UEMUSHORT *a, *b;
2776 UEMUSHORT ai[NI], bi[NI], ci[NI];
2778 EMULONG lt, lta, ltb;
2799 /* compare exponents */
2804 { /* put the larger number in bi */
2814 if (lt < (EMULONG) (-NBITS - 1))
2815 goto done; /* answer same as larger addend */
2817 lost = eshift (ai, k); /* shift the smaller number down */
2821 /* exponents were the same, so must compare significands */
2824 { /* the numbers are identical in magnitude */
2825 /* if different signs, result is zero */
2831 /* if same sign, result is double */
2832 /* double denormalized tiny number */
2833 if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2838 /* add 1 to exponent unless both are zero! */
2839 for (j = 1; j < NI - 1; j++)
2855 bi[E] = (UEMUSHORT) ltb;
2859 { /* put the larger number in bi */
2875 emdnorm (bi, lost, subflg, ltb, !ROUND_TOWARDS_ZERO);
2881 /* Divide: C = B/A, all e type. */
2885 const UEMUSHORT *a, *b;
2888 UEMUSHORT ai[NI], bi[NI];
2890 EMULONG lt, lta, ltb;
2892 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2893 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
2894 sign = eisneg (a) ^ eisneg (b);
2897 /* Return any NaN input. */
2908 /* Zero over zero, or infinity over infinity, is a NaN. */
2909 if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
2910 || (eisinf (a) && eisinf (b)))
2912 mtherr ("ediv", INVALID);
2917 /* Infinity over anything else is infinity. */
2924 /* Anything else over infinity is zero. */
2936 { /* See if numerator is zero. */
2937 for (i = 1; i < NI - 1; i++)
2941 ltb -= enormlz (bi);
2951 { /* possible divide by zero */
2952 for (i = 1; i < NI - 1; i++)
2956 lta -= enormlz (ai);
2960 /* Divide by zero is not an invalid operation.
2961 It is a divide-by-zero operation! */
2963 mtherr ("ediv", SING);
2969 /* calculate exponent */
2970 lt = ltb - lta + EXONE;
2971 emdnorm (bi, i, 0, lt, !ROUND_TOWARDS_ZERO);
2978 && (ecmp (c, ezero) != 0)
2981 *(c+(NE-1)) |= 0x8000;
2983 *(c+(NE-1)) &= ~0x8000;
2986 /* Multiply e-types A and B, return e-type product C. */
2990 const UEMUSHORT *a, *b;
2993 UEMUSHORT ai[NI], bi[NI];
2995 EMULONG lt, lta, ltb;
2997 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2998 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
2999 sign = eisneg (a) ^ eisneg (b);
3002 /* NaN times anything is the same NaN. */
3013 /* Zero times infinity is a NaN. */
3014 if ((eisinf (a) && (ecmp (b, ezero) == 0))
3015 || (eisinf (b) && (ecmp (a, ezero) == 0)))
3017 mtherr ("emul", INVALID);
3022 /* Infinity times anything else is infinity. */
3024 if (eisinf (a) || eisinf (b))
3036 for (i = 1; i < NI - 1; i++)
3040 lta -= enormlz (ai);
3051 for (i = 1; i < NI - 1; i++)
3055 ltb -= enormlz (bi);
3064 /* Multiply significands */
3066 /* calculate exponent */
3067 lt = lta + ltb - (EXONE - 1);
3068 emdnorm (bi, j, 0, lt, !ROUND_TOWARDS_ZERO);
3075 && (ecmp (c, ezero) != 0)
3078 *(c+(NE-1)) |= 0x8000;
3080 *(c+(NE-1)) &= ~0x8000;
3083 /* Convert double precision PE to e-type Y. */
3087 const UEMUSHORT *pe;
3097 ibmtoe (pe, y, DFmode);
3102 c4xtoe (pe, y, HFmode);
3112 denorm = 0; /* flag if denormalized number */
3114 if (! REAL_WORDS_BIG_ENDIAN)
3120 yy[M] = (r & 0x0f) | 0x10;
3121 r &= ~0x800f; /* strip sign and 4 significand bits */
3126 if (! REAL_WORDS_BIG_ENDIAN)
3128 if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
3129 || (pe[1] != 0) || (pe[0] != 0))
3131 enan (y, yy[0] != 0);
3137 if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
3138 || (pe[2] != 0) || (pe[3] != 0))
3140 enan (y, yy[0] != 0);
3151 #endif /* INFINITY */
3153 /* If zero exponent, then the significand is denormalized.
3154 So take back the understood high significand bit. */
3165 if (! REAL_WORDS_BIG_ENDIAN)
3182 /* If zero exponent, then normalize the significand. */
3183 if ((k = enormlz (yy)) > NBITS)
3186 yy[E] -= (UEMUSHORT) (k - 1);
3189 #endif /* not C4X */
3190 #endif /* not IBM */
3191 #endif /* not DEC */
3194 /* Convert double extended precision float PE to e type Y. */
3198 const UEMUSHORT *pe;
3208 for (i = 0; i < NE - 5; i++)
3210 /* This precision is not ordinarily supported on DEC or IBM. */
3212 for (i = 0; i < 5; i++)
3216 p = &yy[0] + (NE - 1);
3219 for (i = 0; i < 5; i++)
3223 if (! REAL_WORDS_BIG_ENDIAN)
3225 for (i = 0; i < 5; i++)
3228 /* For denormal long double Intel format, shift significand up one
3229 -- but only if the top significand bit is zero. A top bit of 1
3230 is "pseudodenormal" when the exponent is zero. */
3231 if ((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
3243 p = &yy[0] + (NE - 1);
3244 #ifdef ARM_EXTENDED_IEEE_FORMAT
3245 /* For ARMs, the exponent is in the lowest 15 bits of the word. */
3246 *p-- = (e[0] & 0x8000) | (e[1] & 0x7ffff);
3252 for (i = 0; i < 4; i++)
3257 /* Point to the exponent field and check max exponent cases. */
3259 if ((*p & 0x7fff) == 0x7fff)
3262 if (! REAL_WORDS_BIG_ENDIAN)
3264 for (i = 0; i < 4; i++)
3266 if ((i != 3 && pe[i] != 0)
3267 /* Anything but 0x8000 here, including 0, is a NaN. */
3268 || (i == 3 && pe[i] != 0x8000))
3270 enan (y, (*p & 0x8000) != 0);
3277 #ifdef ARM_EXTENDED_IEEE_FORMAT
3278 for (i = 2; i <= 5; i++)
3282 enan (y, (*p & 0x8000) != 0);
3287 /* In Motorola extended precision format, the most significant
3288 bit of an infinity mantissa could be either 1 or 0. It is
3289 the lower order bits that tell whether the value is a NaN. */
3290 if ((pe[2] & 0x7fff) != 0)
3293 for (i = 3; i <= 5; i++)
3298 enan (y, (*p & 0x8000) != 0);
3302 #endif /* not ARM */
3311 #endif /* INFINITY */
3314 for (i = 0; i < NE; i++)
3318 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3319 /* Convert 128-bit long double precision float PE to e type Y. */
3323 const UEMUSHORT *pe;
3336 if (! REAL_WORDS_BIG_ENDIAN)
3348 if (! REAL_WORDS_BIG_ENDIAN)
3350 for (i = 0; i < 7; i++)
3354 enan (y, yy[0] != 0);
3361 for (i = 1; i < 8; i++)
3365 enan (y, yy[0] != 0);
3377 #endif /* INFINITY */
3381 if (! REAL_WORDS_BIG_ENDIAN)
3383 for (i = 0; i < 7; i++)
3389 for (i = 0; i < 7; i++)
3393 /* If denormal, remove the implied bit; else shift down 1. */
3407 /* Convert single precision float PE to e type Y. */
3411 const UEMUSHORT *pe;
3416 ibmtoe (pe, y, SFmode);
3422 c4xtoe (pe, y, QFmode);
3433 denorm = 0; /* flag if denormalized number */
3436 if (! REAL_WORDS_BIG_ENDIAN)
3446 yy[M] = (r & 0x7f) | 0200;
3447 r &= ~0x807f; /* strip sign and 7 significand bits */
3449 if (!LARGEST_EXPONENT_IS_NORMAL (32) && r == 0x7f80)
3452 if (REAL_WORDS_BIG_ENDIAN)
3454 if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
3456 enan (y, yy[0] != 0);
3462 if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
3464 enan (y, yy[0] != 0);
3475 #endif /* INFINITY */
3477 /* If zero exponent, then the significand is denormalized.
3478 So take back the understood high significand bit. */
3491 if (! REAL_WORDS_BIG_ENDIAN)
3501 { /* if zero exponent, then normalize the significand */
3502 if ((k = enormlz (yy)) > NBITS)
3505 yy[E] -= (UEMUSHORT) (k - 1);
3508 #endif /* not C4X */
3509 #endif /* not IBM */
3512 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3513 /* Convert e-type X to IEEE 128-bit long double format E. */
3527 make_nan (e, eisneg (x), TFmode);
3532 exp = (EMULONG) xi[E];
3537 /* round off to nearest or even */
3540 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
3548 /* Convert exploded e-type X, that has already been rounded to
3549 113-bit precision, to IEEE 128-bit long double format Y. */
3561 make_nan (b, eiisneg (a), TFmode);
3566 if (REAL_WORDS_BIG_ENDIAN)
3569 q = b + 7; /* point to output exponent */
3571 /* If not denormal, delete the implied bit. */
3576 /* combine sign and exponent */
3578 if (REAL_WORDS_BIG_ENDIAN)
3581 *q++ = *p++ | 0x8000;
3588 *q-- = *p++ | 0x8000;
3592 /* skip over guard word */
3594 /* move the significand */
3595 if (REAL_WORDS_BIG_ENDIAN)
3597 for (i = 0; i < 7; i++)
3602 for (i = 0; i < 7; i++)
3608 /* Convert e-type X to IEEE double extended format E. */
3622 make_nan (e, eisneg (x), XFmode);
3627 /* adjust exponent for offset */
3628 exp = (EMULONG) xi[E];
3633 /* round off to nearest or even */
3636 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
3644 /* Convert exploded e-type X, that has already been rounded to
3645 64-bit precision, to IEEE double extended format Y. */
3657 make_nan (b, eiisneg (a), XFmode);
3661 /* Shift denormal long double Intel format significand down one bit. */
3662 if ((a[E] == 0) && ! REAL_WORDS_BIG_ENDIAN)
3672 if (REAL_WORDS_BIG_ENDIAN)
3676 q = b + 4; /* point to output exponent */
3677 /* Clear the last two bytes of 12-byte Intel format. q is pointing
3678 into an array of size 6 (e.g. x[NE]), so the last two bytes are
3679 always there, and there are never more bytes, even when we are using
3680 INTEL_EXTENDED_IEEE_FORMAT. */
3685 /* combine sign and exponent */
3689 *q++ = *p++ | 0x8000;
3696 *q-- = *p++ | 0x8000;
3701 if (REAL_WORDS_BIG_ENDIAN)
3703 #ifdef ARM_EXTENDED_IEEE_FORMAT
3704 /* The exponent is in the lowest 15 bits of the first word. */
3705 *q++ = i ? 0x8000 : 0;
3709 *q++ = *p++ | 0x8000;
3718 *q-- = *p++ | 0x8000;
3723 /* skip over guard word */
3725 /* move the significand */
3727 for (i = 0; i < 4; i++)
3731 for (i = 0; i < 4; i++)
3735 if (REAL_WORDS_BIG_ENDIAN)
3737 for (i = 0; i < 4; i++)
3745 /* Intel long double infinity significand. */
3753 for (i = 0; i < 4; i++)
3759 /* e type to double precision. */
3762 /* Convert e-type X to DEC-format double E. */
3769 etodec (x, e); /* see etodec.c */
3772 /* Convert exploded e-type X, that has already been rounded to
3773 56-bit double precision, to DEC double Y. */
3784 /* Convert e-type X to IBM 370-format double E. */
3791 etoibm (x, e, DFmode);
3794 /* Convert exploded e-type X, that has already been rounded to
3795 56-bit precision, to IBM 370 double Y. */
3801 toibm (x, y, DFmode);
3804 #else /* it's neither DEC nor IBM */
3806 /* Convert e-type X to C4X-format long double E. */
3813 etoc4x (x, e, HFmode);
3816 /* Convert exploded e-type X, that has already been rounded to
3817 56-bit precision, to IBM 370 double Y. */
3823 toc4x (x, y, HFmode);
3826 #else /* it's neither DEC nor IBM nor C4X */
3828 /* Convert e-type X to IEEE double E. */
3842 make_nan (e, eisneg (x), DFmode);
3847 /* adjust exponent for offsets */
3848 exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
3853 /* round off to nearest or even */
3856 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
3864 /* Convert exploded e-type X, that has already been rounded to
3865 53-bit precision, to IEEE double Y. */
3877 make_nan (y, eiisneg (x), DFmode);
3881 if (LARGEST_EXPONENT_IS_NORMAL (64) && x[1] > 2047)
3883 saturate (y, eiisneg (x), 64, 1);
3888 if (! REAL_WORDS_BIG_ENDIAN)
3891 *y = 0; /* output high order */
3893 *y = 0x8000; /* output sign bit */
3896 if (i >= (unsigned int) 2047)
3898 /* Saturate at largest number less than infinity. */
3901 if (! REAL_WORDS_BIG_ENDIAN)
3915 *y |= (UEMUSHORT) 0x7fef;
3916 if (! REAL_WORDS_BIG_ENDIAN)
3941 i |= *p++ & (UEMUSHORT) 0x0f; /* *p = xi[M] */
3942 *y |= (UEMUSHORT) i; /* high order output already has sign bit set */
3943 if (! REAL_WORDS_BIG_ENDIAN)
3958 #endif /* not C4X */
3959 #endif /* not IBM */
3960 #endif /* not DEC */
3964 /* e type to single precision. */
3967 /* Convert e-type X to IBM 370 float E. */
3974 etoibm (x, e, SFmode);
3977 /* Convert exploded e-type X, that has already been rounded to
3978 float precision, to IBM 370 float Y. */
3984 toibm (x, y, SFmode);
3990 /* Convert e-type X to C4X float E. */
3997 etoc4x (x, e, QFmode);
4000 /* Convert exploded e-type X, that has already been rounded to
4001 float precision, to IBM 370 float Y. */
4007 toc4x (x, y, QFmode);
4012 /* Convert e-type X to IEEE float E. DEC float is the same as IEEE float. */
4026 make_nan (e, eisneg (x), SFmode);
4031 /* adjust exponent for offsets */
4032 exp = (EMULONG) xi[E] - (EXONE - 0177);
4037 /* round off to nearest or even */
4040 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
4048 /* Convert exploded e-type X, that has already been rounded to
4049 float precision, to IEEE float Y. */
4061 make_nan (y, eiisneg (x), SFmode);
4065 if (LARGEST_EXPONENT_IS_NORMAL (32) && x[1] > 255)
4067 saturate (y, eiisneg (x), 32, 1);
4072 if (! REAL_WORDS_BIG_ENDIAN)
4078 *y = 0; /* output high order */
4080 *y = 0x8000; /* output sign bit */
4083 /* Handle overflow cases. */
4084 if (!LARGEST_EXPONENT_IS_NORMAL (32) && i >= 255)
4087 *y |= (UEMUSHORT) 0x7f80;
4092 if (! REAL_WORDS_BIG_ENDIAN)
4100 #else /* no INFINITY */
4101 *y |= (UEMUSHORT) 0x7f7f;
4106 if (! REAL_WORDS_BIG_ENDIAN)
4117 #endif /* no INFINITY */
4129 i |= *p++ & (UEMUSHORT) 0x7f; /* *p = xi[M] */
4130 /* High order output already has sign bit set. */
4136 if (! REAL_WORDS_BIG_ENDIAN)
4145 #endif /* not C4X */
4146 #endif /* not IBM */
4148 /* Compare two e type numbers.
4152 -2 if either a or b is a NaN. */
4156 const UEMUSHORT *a, *b;
4158 UEMUSHORT ai[NI], bi[NI];
4164 if (eisnan (a) || eisnan (b))
4173 { /* the signs are different */
4175 for (i = 1; i < NI - 1; i++)
4189 /* both are the same sign */
4204 return (0); /* equality */
4208 if (*(--p) > *(--q))
4209 return (msign); /* p is bigger */
4211 return (-msign); /* p is littler */
4215 /* Find e-type nearest integer to X, as floor (X + 0.5). */
4227 /* Convert HOST_WIDE_INT LP to e type Y. */
4231 const HOST_WIDE_INT *lp;
4235 unsigned HOST_WIDE_INT ll;
4241 /* make it positive */
4242 ll = (unsigned HOST_WIDE_INT) (-(*lp));
4243 yi[0] = 0xffff; /* put correct sign in the e type number */
4247 ll = (unsigned HOST_WIDE_INT) (*lp);
4249 /* move the long integer to yi significand area */
4250 #if HOST_BITS_PER_WIDE_INT == 64
4251 yi[M] = (UEMUSHORT) (ll >> 48);
4252 yi[M + 1] = (UEMUSHORT) (ll >> 32);
4253 yi[M + 2] = (UEMUSHORT) (ll >> 16);
4254 yi[M + 3] = (UEMUSHORT) ll;
4255 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
4257 yi[M] = (UEMUSHORT) (ll >> 16);
4258 yi[M + 1] = (UEMUSHORT) ll;
4259 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
4262 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4263 ecleaz (yi); /* it was zero */
4265 yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
4266 emovo (yi, y); /* output the answer */
4269 /* Convert unsigned HOST_WIDE_INT LP to e type Y. */
4273 const unsigned HOST_WIDE_INT *lp;
4277 unsigned HOST_WIDE_INT ll;
4283 /* move the long integer to ayi significand area */
4284 #if HOST_BITS_PER_WIDE_INT == 64
4285 yi[M] = (UEMUSHORT) (ll >> 48);
4286 yi[M + 1] = (UEMUSHORT) (ll >> 32);
4287 yi[M + 2] = (UEMUSHORT) (ll >> 16);
4288 yi[M + 3] = (UEMUSHORT) ll;
4289 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
4291 yi[M] = (UEMUSHORT) (ll >> 16);
4292 yi[M + 1] = (UEMUSHORT) ll;
4293 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
4296 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4297 ecleaz (yi); /* it was zero */
4299 yi[E] -= (UEMUSHORT) k; /* subtract shift count from exponent */
4300 emovo (yi, y); /* output the answer */
4304 /* Find signed HOST_WIDE_INT integer I and floating point fractional
4305 part FRAC of e-type (packed internal format) floating point input X.
4306 The integer output I has the sign of the input, except that
4307 positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
4308 The output e-type fraction FRAC is the positive fractional
4319 unsigned HOST_WIDE_INT ll;
4322 k = (int) xi[E] - (EXONE - 1);
4325 /* if exponent <= 0, integer = 0 and real output is fraction */
4330 if (k > (HOST_BITS_PER_WIDE_INT - 1))
4332 /* long integer overflow: output large integer
4333 and correct fraction */
4335 *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
4338 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
4339 /* In this case, let it overflow and convert as if unsigned. */
4340 euifrac (x, &ll, frac);
4341 *i = (HOST_WIDE_INT) ll;
4344 /* In other cases, return the largest positive integer. */
4345 *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
4350 warning ("overflow on truncation to integer");
4354 /* Shift more than 16 bits: first shift up k-16 mod 16,
4355 then shift up by 16's. */
4356 j = k - ((k >> 4) << 4);
4363 ll = (ll << 16) | xi[M];
4365 while ((k -= 16) > 0);
4372 /* shift not more than 16 bits */
4374 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4381 if ((k = enormlz (xi)) > NBITS)
4384 xi[E] -= (UEMUSHORT) k;
4390 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
4391 FRAC of e-type X. A negative input yields integer output = 0 but
4392 correct fraction. */
4395 euifrac (x, i, frac)
4397 unsigned HOST_WIDE_INT *i;
4400 unsigned HOST_WIDE_INT ll;
4405 k = (int) xi[E] - (EXONE - 1);
4408 /* if exponent <= 0, integer = 0 and argument is fraction */
4413 if (k > HOST_BITS_PER_WIDE_INT)
4415 /* Long integer overflow: output large integer
4416 and correct fraction.
4417 Note, the BSD MicroVAX compiler says that ~(0UL)
4418 is a syntax error. */
4422 warning ("overflow on truncation to unsigned integer");
4426 /* Shift more than 16 bits: first shift up k-16 mod 16,
4427 then shift up by 16's. */
4428 j = k - ((k >> 4) << 4);
4435 ll = (ll << 16) | xi[M];
4437 while ((k -= 16) > 0);
4442 /* shift not more than 16 bits */
4444 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4447 if (xi[0]) /* A negative value yields unsigned integer 0. */
4453 if ((k = enormlz (xi)) > NBITS)
4456 xi[E] -= (UEMUSHORT) k;
4461 /* Shift the significand of exploded e-type X up or down by SC bits. */
4482 lost |= *p; /* remember lost bits */
4523 return ((int) lost);
4526 /* Shift normalize the significand area of exploded e-type X.
4527 Return the shift count (up = positive). */
4542 return (0); /* already normalized */
4548 /* With guard word, there are NBITS+16 bits available.
4549 Return true if all are zero. */
4553 /* see if high byte is zero */
4554 while ((*p & 0xff00) == 0)
4559 /* now shift 1 bit at a time */
4560 while ((*p & 0x8000) == 0)
4566 mtherr ("enormlz", UNDERFLOW);
4572 /* Normalize by shifting down out of the high guard word
4573 of the significand */
4588 mtherr ("enormlz", OVERFLOW);
4595 /* Powers of ten used in decimal <-> binary conversions. */
4600 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
4601 static const UEMUSHORT etens[NTEN + 1][NE] =
4603 {0x6576, 0x4a92, 0x804a, 0x153f,
4604 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4605 {0x6a32, 0xce52, 0x329a, 0x28ce,
4606 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4607 {0x526c, 0x50ce, 0xf18b, 0x3d28,
4608 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4609 {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4610 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4611 {0x851e, 0xeab7, 0x98fe, 0x901b,
4612 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4613 {0x0235, 0x0137, 0x36b1, 0x336c,
4614 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4615 {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4616 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4617 {0x0000, 0x0000, 0x0000, 0x0000,
4618 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4619 {0x0000, 0x0000, 0x0000, 0x0000,
4620 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4621 {0x0000, 0x0000, 0x0000, 0x0000,
4622 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4623 {0x0000, 0x0000, 0x0000, 0x0000,
4624 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4625 {0x0000, 0x0000, 0x0000, 0x0000,
4626 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4627 {0x0000, 0x0000, 0x0000, 0x0000,
4628 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4631 static const UEMUSHORT emtens[NTEN + 1][NE] =
4633 {0x2030, 0xcffc, 0xa1c3, 0x8123,
4634 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4635 {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4636 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4637 {0xf53f, 0xf698, 0x6bd3, 0x0158,
4638 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4639 {0xe731, 0x04d4, 0xe3f2, 0xd332,
4640 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4641 {0xa23e, 0x5308, 0xfefb, 0x1155,
4642 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4643 {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4644 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4645 {0x2a20, 0x6224, 0x47b3, 0x98d7,
4646 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4647 {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4648 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4649 {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4650 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4651 {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4652 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4653 {0xc155, 0xa4a8, 0x404e, 0x6113,
4654 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4655 {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4656 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4657 {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4658 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4661 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4662 static const UEMUSHORT etens[NTEN + 1][NE] =
4664 {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4665 {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4666 {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4667 {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4668 {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4669 {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4670 {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4671 {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4672 {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4673 {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4674 {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4675 {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4676 {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4679 static const UEMUSHORT emtens[NTEN + 1][NE] =
4681 {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4682 {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4683 {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4684 {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4685 {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4686 {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4687 {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4688 {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4689 {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4690 {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4691 {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4692 {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4693 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4698 /* Convert float value X to ASCII string STRING with NDIG digits after
4699 the decimal point. */
4702 e24toasc (x, string, ndigs)
4703 const UEMUSHORT x[];
4710 etoasc (w, string, ndigs);
4713 /* Convert double value X to ASCII string STRING with NDIG digits after
4714 the decimal point. */
4717 e53toasc (x, string, ndigs)
4718 const UEMUSHORT x[];
4725 etoasc (w, string, ndigs);
4728 /* Convert double extended value X to ASCII string STRING with NDIG digits
4729 after the decimal point. */
4732 e64toasc (x, string, ndigs)
4733 const UEMUSHORT x[];
4740 etoasc (w, string, ndigs);
4743 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4744 after the decimal point. */
4747 e113toasc (x, string, ndigs)
4748 const UEMUSHORT x[];
4755 etoasc (w, string, ndigs);
4759 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4760 the decimal point. */
4762 static char wstring[80]; /* working storage for ASCII output */
4765 etoasc (x, string, ndigs)
4766 const UEMUSHORT x[];
4771 UEMUSHORT y[NI], t[NI], u[NI], w[NI];
4772 const UEMUSHORT *p, *r, *ten;
4774 int i, j, k, expon, rndsav;
4787 sprintf (wstring, " NaN ");
4791 rndprc = NBITS; /* set to full precision */
4792 emov (x, y); /* retain external format */
4793 if (y[NE - 1] & 0x8000)
4796 y[NE - 1] &= 0x7fff;
4803 ten = &etens[NTEN][0];
4805 /* Test for zero exponent */
4808 for (k = 0; k < NE - 1; k++)
4811 goto tnzro; /* denormalized number */
4813 goto isone; /* valid all zeros */
4817 /* Test for infinity. */
4818 if (y[NE - 1] == 0x7fff)
4821 sprintf (wstring, " -Infinity ");
4823 sprintf (wstring, " Infinity ");
4827 /* Test for exponent nonzero but significand denormalized.
4828 * This is an error condition.
4830 if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4832 mtherr ("etoasc", DOMAIN);
4833 sprintf (wstring, "NaN");
4837 /* Compare to 1.0 */
4846 { /* Number is greater than 1 */
4847 /* Convert significand to an integer and strip trailing decimal zeros. */
4849 u[NE - 1] = EXONE + NBITS - 1;
4851 p = &etens[NTEN - 4][0];
4857 for (j = 0; j < NE - 1; j++)
4870 /* Rescale from integer significand */
4871 u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
4873 /* Find power of 10 */
4877 /* An unordered compare result shouldn't happen here. */
4878 while (ecmp (ten, u) <= 0)
4880 if (ecmp (p, u) <= 0)
4893 { /* Number is less than 1.0 */
4894 /* Pad significand with trailing decimal zeros. */
4897 while ((y[NE - 2] & 0x8000) == 0)
4906 for (i = 0; i < NDEC + 1; i++)
4908 if ((w[NI - 1] & 0x7) != 0)
4910 /* multiply by 10 */
4923 if (eone[NE - 1] <= u[1])
4935 while (ecmp (eone, w) > 0)
4937 if (ecmp (p, w) >= 0)
4952 /* Find the first (leading) digit. */
4958 digit = equot[NI - 1];
4959 while ((digit == 0) && (ecmp (y, ezero) != 0))
4967 digit = equot[NI - 1];
4975 /* Examine number of digits requested by caller. */
4993 *s++ = (char) digit + '0';
4996 /* Generate digits after the decimal point. */
4997 for (k = 0; k <= ndigs; k++)
4999 /* multiply current number by 10, without normalizing */
5006 *s++ = (char) equot[NI - 1] + '0';
5008 digit = equot[NI - 1];
5011 /* round off the ASCII string */
5014 /* Test for critical rounding case in ASCII output. */
5018 if (ecmp (t, ezero) != 0)
5019 goto roun; /* round to nearest */
5021 if ((*(s - 1) & 1) == 0)
5022 goto doexp; /* round to even */
5025 /* Round up and propagate carry-outs */
5029 /* Carry out to most significant digit? */
5036 /* Most significant digit carries to 10? */
5044 /* Round up and carry out from less significant digits */
5056 sprintf (ss, "e+%d", expon);
5058 sprintf (ss, "e%d", expon);
5060 sprintf (ss, "e%d", expon);
5063 /* copy out the working string */
5066 while (*ss == ' ') /* strip possible leading space */
5068 while ((*s++ = *ss++) != '\0')
5073 /* Convert ASCII string to floating point.
5075 Numeric input is a free format decimal number of any length, with
5076 or without decimal point. Entering E after the number followed by an
5077 integer number causes the second number to be interpreted as a power of
5078 10 to be multiplied by the first number (i.e., "scientific" notation). */
5080 /* Convert ASCII string S to single precision float value Y. */
5091 /* Convert ASCII string S to double precision value Y. */
5098 #if defined(DEC) || defined(IBM)
5110 /* Convert ASCII string S to double extended value Y. */
5120 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
5121 /* Convert ASCII string S to 128-bit long double Y. */
5128 asctoeg (s, y, 113);
5132 /* Convert ASCII string S to e type Y. */
5139 asctoeg (s, y, NBITS);
5142 /* Convert ASCII string SS to e type Y, with a specified rounding precision
5143 of OPREC bits. BASE is 16 for C99 hexadecimal floating constants. */
5146 asctoeg (ss, y, oprec)
5151 UEMUSHORT yy[NI], xt[NI], tt[NI];
5152 int esign, decflg, sgnflg, nexp, exp, prec, lost;
5153 int i, k, trail, c, rndsav;
5156 char *sp, *s, *lstr;
5159 /* Copy the input string. */
5160 lstr = (char *) alloca (strlen (ss) + 1);
5162 while (*ss == ' ') /* skip leading spaces */
5166 while ((*sp++ = *ss++) != '\0')
5170 if (s[0] == '0' && (s[1] == 'x' || s[1] == 'X'))
5177 rndprc = NBITS; /* Set to full precision */
5190 if ((k >= 0) && (k < base))
5192 /* Ignore leading zeros */
5193 if ((prec == 0) && (decflg == 0) && (k == 0))
5195 /* Identify and strip trailing zeros after the decimal point. */
5196 if ((trail == 0) && (decflg != 0))
5199 while (ISDIGIT (*sp) || (base == 16 && ISXDIGIT (*sp)))
5201 /* Check for syntax error */
5203 if ((base != 10 || ((c != 'e') && (c != 'E')))
5204 && (base != 16 || ((c != 'p') && (c != 'P')))
5206 && (c != '\n') && (c != '\r') && (c != ' ')
5208 goto unexpected_char_error;
5217 /* If enough digits were given to more than fill up the yy register,
5218 continuing until overflow into the high guard word yy[2]
5219 guarantees that there will be a roundoff bit at the top
5220 of the low guard word after normalization. */
5227 nexp += 4; /* count digits after decimal point */
5229 eshup1 (yy); /* multiply current number by 16 */
5237 nexp += 1; /* count digits after decimal point */
5239 eshup1 (yy); /* multiply current number by 10 */
5245 /* Insert the current digit. */
5247 xt[NI - 2] = (UEMUSHORT) k;
5252 /* Mark any lost non-zero digit. */
5254 /* Count lost digits before the decimal point. */
5276 case '.': /* decimal point */
5278 goto unexpected_char_error;
5284 goto unexpected_char_error;
5289 goto unexpected_char_error;
5302 unexpected_char_error:
5306 mtherr ("asctoe", DOMAIN);
5315 /* Exponent interpretation */
5317 /* 0.0eXXX is zero, regardless of XXX. Check for the 0.0. */
5318 for (k = 0; k < NI; k++)
5329 /* check for + or - */
5337 while (ISDIGIT (*s))
5346 if ((exp > MAXDECEXP) && (base == 10))
5350 yy[E] = 0x7fff; /* infinity */
5353 if ((exp < MINDECEXP) && (base == 10))
5363 /* Base 16 hexadecimal floating constant. */
5364 if ((k = enormlz (yy)) > NBITS)
5369 /* Adjust the exponent. NEXP is the number of hex digits,
5370 EXP is a power of 2. */
5371 lexp = (EXONE - 1 + NBITS) - k + yy[E] + exp - nexp;
5381 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
5382 while ((nexp > 0) && (yy[2] == 0))
5394 if ((k = enormlz (yy)) > NBITS)
5399 lexp = (EXONE - 1 + NBITS) - k;
5400 emdnorm (yy, lost, 0, lexp, 64);
5403 /* Convert to external format:
5405 Multiply by 10**nexp. If precision is 64 bits,
5406 the maximum relative error incurred in forming 10**n
5407 for 0 <= n <= 324 is 8.2e-20, at 10**180.
5408 For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5409 For 0 >= n >= -999, it is -1.55e-19 at 10**-435. */
5424 /* Punt. Can't handle this without 2 divides. */
5425 emovi (etens[0], tt);
5438 emul (etens[i], xt, xt);
5442 while (exp <= MAXP);
5461 /* Round and convert directly to the destination type */
5463 lexp -= EXONE - 0x3ff;
5465 else if (oprec == 24 || oprec == 32)
5466 lexp -= (EXONE - 0x7f);
5469 else if (oprec == 24 || oprec == 56)
5470 lexp -= EXONE - (0x41 << 2);
5472 else if (oprec == 24)
5473 lexp -= EXONE - 0177;
5477 else if (oprec == 56)
5478 lexp -= EXONE - 0201;
5481 emdnorm (yy, lost, 0, lexp, 64);
5491 todec (yy, y); /* see etodec.c */
5496 toibm (yy, y, DFmode);
5501 toc4x (yy, y, HFmode);
5514 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
5527 /* Return Y = largest integer not greater than X (truncated toward minus
5530 static const UEMUSHORT bmask[] =
5553 const UEMUSHORT x[];
5560 emov (x, f); /* leave in external format */
5561 expon = (int) f[NE - 1];
5562 e = (expon & 0x7fff) - (EXONE - 1);
5568 /* number of bits to clear out */
5580 /* clear the remaining bits */
5582 /* truncate negatives toward minus infinity */
5585 if ((UEMUSHORT) expon & (UEMUSHORT) 0x8000)
5587 for (i = 0; i < NE - 1; i++)
5600 /* Return S and EXP such that S * 2^EXP = X and .5 <= S < 1.
5601 For example, 1.1 = 0.55 * 2^1. */
5605 const UEMUSHORT x[];
5613 /* Handle denormalized numbers properly using long integer exponent. */
5614 li = (EMULONG) ((EMUSHORT) xi[1]);
5622 *exp = (int) (li - 0x3ffe);
5626 /* Return e type Y = X * 2^PWR2. */
5630 const UEMUSHORT x[];
5642 emdnorm (xi, i, i, li, !ROUND_TOWARDS_ZERO);
5648 /* C = remainder after dividing B by A, all e type values.
5649 Least significant integer quotient bits left in EQUOT. */
5653 const UEMUSHORT a[], b[];
5656 UEMUSHORT den[NI], num[NI];
5660 || (ecmp (a, ezero) == 0)
5668 if (ecmp (a, ezero) == 0)
5670 mtherr ("eremain", SING);
5676 eiremain (den, num);
5677 /* Sign of remainder = sign of quotient */
5686 /* Return quotient of exploded e-types NUM / DEN in EQUOT,
5687 remainder in NUM. */
5691 UEMUSHORT den[], num[];
5697 ld -= enormlz (den);
5699 ln -= enormlz (num);
5703 if (ecmpm (den, num) <= 0)
5715 emdnorm (num, 0, 0, ln, 0);
5718 /* Report an error condition CODE encountered in function NAME.
5720 Mnemonic Value Significance
5722 DOMAIN 1 argument domain error
5723 SING 2 function singularity
5724 OVERFLOW 3 overflow range error
5725 UNDERFLOW 4 underflow range error
5726 TLOSS 5 total loss of precision
5727 PLOSS 6 partial loss of precision
5728 INVALID 7 NaN - producing operation
5729 EDOM 33 Unix domain error code
5730 ERANGE 34 Unix range error code
5732 The order of appearance of the following messages is bound to the
5733 error codes defined above. */
5743 /* The string passed by the calling program is supposed to be the
5744 name of the function in which the error occurred.
5745 The code argument selects which error message string will be printed. */
5747 if (strcmp (name, "esub") == 0)
5748 name = "subtraction";
5749 else if (strcmp (name, "ediv") == 0)
5751 else if (strcmp (name, "emul") == 0)
5752 name = "multiplication";
5753 else if (strcmp (name, "enormlz") == 0)
5754 name = "normalization";
5755 else if (strcmp (name, "etoasc") == 0)
5756 name = "conversion to text";
5757 else if (strcmp (name, "asctoe") == 0)
5759 else if (strcmp (name, "eremain") == 0)
5761 else if (strcmp (name, "esqrt") == 0)
5762 name = "square root";
5767 case DOMAIN: warning ("%s: argument domain error" , name); break;
5768 case SING: warning ("%s: function singularity" , name); break;
5769 case OVERFLOW: warning ("%s: overflow range error" , name); break;
5770 case UNDERFLOW: warning ("%s: underflow range error" , name); break;
5771 case TLOSS: warning ("%s: total loss of precision" , name); break;
5772 case PLOSS: warning ("%s: partial loss of precision", name); break;
5773 case INVALID: warning ("%s: NaN - producing operation", name); break;
5778 /* Set global error message word */
5783 /* Convert DEC double precision D to e type E. */
5793 ecleaz (y); /* start with a zero */
5794 p = y; /* point to our number */
5795 r = *d; /* get DEC exponent word */
5796 if (*d & (unsigned int) 0x8000)
5797 *p = 0xffff; /* fill in our sign */
5798 ++p; /* bump pointer to our exponent word */
5799 r &= 0x7fff; /* strip the sign bit */
5800 if (r == 0) /* answer = 0 if high order DEC word = 0 */
5804 r >>= 7; /* shift exponent word down 7 bits */
5805 r += EXONE - 0201; /* subtract DEC exponent offset */
5806 /* add our e type exponent offset */
5807 *p++ = r; /* to form our exponent */
5809 r = *d++; /* now do the high order mantissa */
5810 r &= 0177; /* strip off the DEC exponent and sign bits */
5811 r |= 0200; /* the DEC understood high order mantissa bit */
5812 *p++ = r; /* put result in our high guard word */
5814 *p++ = *d++; /* fill in the rest of our mantissa */
5818 eshdn8 (y); /* shift our mantissa down 8 bits */
5823 /* Convert e type X to DEC double precision D. */
5835 /* Adjust exponent for offsets. */
5836 exp = (EMULONG) xi[E] - (EXONE - 0201);
5837 /* Round off to nearest or even. */
5840 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
5845 /* Convert exploded e-type X, that has already been rounded to
5846 56-bit precision, to DEC format double Y. */
5892 /* Convert IBM single/double precision to e type. */
5898 enum machine_mode mode;
5903 ecleaz (y); /* start with a zero */
5904 p = y; /* point to our number */
5905 r = *d; /* get IBM exponent word */
5906 if (*d & (unsigned int) 0x8000)
5907 *p = 0xffff; /* fill in our sign */
5908 ++p; /* bump pointer to our exponent word */
5909 r &= 0x7f00; /* strip the sign bit */
5910 r >>= 6; /* shift exponent word down 6 bits */
5911 /* in fact shift by 8 right and 2 left */
5912 r += EXONE - (0x41 << 2); /* subtract IBM exponent offset */
5913 /* add our e type exponent offset */
5914 *p++ = r; /* to form our exponent */
5916 *p++ = *d++ & 0xff; /* now do the high order mantissa */
5917 /* strip off the IBM exponent and sign bits */
5918 if (mode != SFmode) /* there are only 2 words in SFmode */
5920 *p++ = *d++; /* fill in the rest of our mantissa */
5925 if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5928 y[E] -= 5 + enormlz (y); /* now normalise the mantissa */
5929 /* handle change in RADIX */
5935 /* Convert e type to IBM single/double precision. */
5941 enum machine_mode mode;
5948 exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2)); /* adjust exponent for offsets */
5949 /* round off to nearest or even */
5952 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
5954 toibm (xi, d, mode);
5960 enum machine_mode mode;
6013 /* Convert C4X single/double precision to e type. */
6019 enum machine_mode mode;
6037 /* Short-circuit the zero case. */
6038 if ((dn[0] == 0x8000)
6039 && (dn[1] == 0x0000)
6040 && ((mode == QFmode) || ((dn[2] == 0x0000) && (dn[3] == 0x0000))))
6051 ecleaz (y); /* start with a zero */
6052 r = dn[0]; /* get sign/exponent part */
6053 if (r & (unsigned int) 0x0080)
6055 y[0] = 0xffff; /* fill in our sign */
6061 r >>= 8; /* Shift exponent word down 8 bits. */
6062 if (r & 0x80) /* Make the exponent negative if it is. */
6063 r = r | (~0 & ~0xff);
6067 /* Now do the high order mantissa. We don't "or" on the high bit
6068 because it is 2 (not 1) and is handled a little differently
6070 y[M] = dn[0] & 0x7f;
6073 if (mode != QFmode) /* There are only 2 words in QFmode. */
6075 y[M+2] = dn[2]; /* Fill in the rest of our mantissa. */
6083 /* Now do the two's complement on the data. */
6085 carry = 1; /* Initially add 1 for the two's complement. */
6086 for (i=size + M; i > M; i--)
6088 if (carry && (y[i] == 0x0000))
6089 /* We overflowed into the next word, carry is the same. */
6090 y[i] = carry ? 0x0000 : 0xffff;
6093 /* No overflow, just invert and add carry. */
6094 y[i] = ((~y[i]) + carry) & 0xffff;
6109 /* Add our e type exponent offset to form our exponent. */
6113 /* Now do the high order mantissa strip off the exponent and sign
6114 bits and add the high 1 bit. */
6115 y[M] = (dn[0] & 0x7f) | 0x80;
6118 if (mode != QFmode) /* There are only 2 words in QFmode. */
6120 y[M+2] = dn[2]; /* Fill in the rest of our mantissa. */
6130 /* Convert e type to C4X single/double precision. */
6136 enum machine_mode mode;
6144 /* Adjust exponent for offsets. */
6145 exp = (EMULONG) xi[E] - (EXONE - 0x7f);
6147 /* Round off to nearest or even. */
6149 rndprc = mode == QFmode ? 24 : 32;
6150 emdnorm (xi, 0, 0, exp, !ROUND_TOWARDS_ZERO);
6152 toc4x (xi, d, mode);
6158 enum machine_mode mode;
6164 /* Short-circuit the zero case */
6165 if ((x[0] == 0) /* Zero exponent and sign */
6167 && (x[M] == 0) /* The rest is for zero mantissa */
6169 /* Only check for double if necessary */
6170 && ((mode == QFmode) || ((x[M+2] == 0) && (x[M+3] == 0))))
6172 /* We have a zero. Put it into the output and return. */
6185 /* Negative number require a two's complement conversion of the
6191 i = ((int) x[1]) - 0x7f;
6193 /* Now add 1 to the inverted data to do the two's complement. */
6202 x[v] = carry ? 0x0000 : 0xffff;
6205 x[v] = ((~x[v]) + carry) & 0xffff;
6211 /* The following is a special case. The C4X negative float requires
6212 a zero in the high bit (because the format is (2 - x) x 2^m), so
6213 if a one is in that bit, we have to shift left one to get rid
6214 of it. This only occurs if the number is -1 x 2^m. */
6215 if (x[M+1] & 0x8000)
6217 /* This is the case of -1 x 2^m, we have to rid ourselves of the
6218 high sign bit and shift the exponent. */
6224 i = ((int) x[1]) - 0x7f;
6226 if ((i < -128) || (i > 127))
6234 y[3] = (y[1] << 8) | ((y[2] >> 8) & 0xff);
6235 y[2] = (y[0] << 8) | ((y[1] >> 8) & 0xff);
6243 y[0] |= ((i & 0xff) << 8);
6247 y[0] |= x[M] & 0x7f;
6253 y[3] = (y[1] << 8) | ((y[2] >> 8) & 0xff);
6254 y[2] = (y[0] << 8) | ((y[1] >> 8) & 0xff);
6259 /* Output a binary NaN bit pattern in the target machine's format. */
6261 /* If special NaN bit patterns are required, define them in tm.h
6262 as arrays of unsigned 16-bit shorts. Otherwise, use the default
6268 static const UEMUSHORT TFbignan[8] =
6269 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6270 static const UEMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
6278 static const UEMUSHORT XFbignan[6] =
6279 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6280 static const UEMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
6288 static const UEMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
6289 static const UEMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
6297 static const UEMUSHORT SFbignan[2] = {0x7fff, 0xffff};
6298 static const UEMUSHORT SFlittlenan[2] = {0, 0xffc0};
6305 make_nan (nan, sign, mode)
6308 enum machine_mode mode;
6314 size = GET_MODE_BITSIZE (mode);
6315 if (LARGEST_EXPONENT_IS_NORMAL (size))
6317 warning ("%d-bit floats cannot hold NaNs", size);
6318 saturate (nan, sign, size, 0);
6323 /* Possibly the `reserved operand' patterns on a VAX can be
6324 used like NaN's, but probably not in the same way as IEEE. */
6325 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
6327 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6329 if (REAL_WORDS_BIG_ENDIAN)
6339 if (REAL_WORDS_BIG_ENDIAN)
6347 if (REAL_WORDS_BIG_ENDIAN)
6356 if (REAL_WORDS_BIG_ENDIAN)
6366 if (REAL_WORDS_BIG_ENDIAN)
6367 *nan++ = (sign << 15) | (*p++ & 0x7fff);
6370 if (! REAL_WORDS_BIG_ENDIAN)
6371 *nan = (sign << 15) | (*p & 0x7fff);
6376 /* Create a saturation value for a SIZE-bit float, assuming that
6377 LARGEST_EXPONENT_IS_NORMAL (SIZE).
6379 If SIGN is true, fill X with the most negative value, otherwise fill
6380 it with the most positive value. WARN is true if the function should
6381 warn about overflow. */
6384 saturate (x, sign, size, warn)
6386 int sign, size, warn;
6390 if (warn && extra_warnings)
6391 warning ("value exceeds the range of a %d-bit float", size);
6393 /* Create the most negative value. */
6394 for (i = 0; i < size / EMUSHORT_SIZE; i++)
6397 /* Make it positive, if necessary. */
6399 x[REAL_WORDS_BIG_ENDIAN? 0 : i - 1] = 0x7fff;
6403 /* This is the inverse of the function `etarsingle' invoked by
6404 REAL_VALUE_TO_TARGET_SINGLE. */
6407 ereal_unto_float (f)
6414 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6415 This is the inverse operation to what the function `endian' does. */
6416 if (REAL_WORDS_BIG_ENDIAN)
6418 s[0] = (UEMUSHORT) (f >> 16);
6419 s[1] = (UEMUSHORT) f;
6423 s[0] = (UEMUSHORT) f;
6424 s[1] = (UEMUSHORT) (f >> 16);
6426 /* Convert and promote the target float to E-type. */
6428 /* Output E-type to REAL_VALUE_TYPE. */
6434 /* This is the inverse of the function `etardouble' invoked by
6435 REAL_VALUE_TO_TARGET_DOUBLE. */
6438 ereal_unto_double (d)
6445 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
6446 if (REAL_WORDS_BIG_ENDIAN)
6448 s[0] = (UEMUSHORT) (d[0] >> 16);
6449 s[1] = (UEMUSHORT) d[0];
6450 s[2] = (UEMUSHORT) (d[1] >> 16);
6451 s[3] = (UEMUSHORT) d[1];
6455 /* Target float words are little-endian. */
6456 s[0] = (UEMUSHORT) d[0];
6457 s[1] = (UEMUSHORT) (d[0] >> 16);
6458 s[2] = (UEMUSHORT) d[1];
6459 s[3] = (UEMUSHORT) (d[1] >> 16);
6461 /* Convert target double to E-type. */
6463 /* Output E-type to REAL_VALUE_TYPE. */
6469 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
6470 This is somewhat like ereal_unto_float, but the input types
6471 for these are different. */
6474 ereal_from_float (f)
6481 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6482 This is the inverse operation to what the function `endian' does. */
6483 if (REAL_WORDS_BIG_ENDIAN)
6485 s[0] = (UEMUSHORT) (f >> 16);
6486 s[1] = (UEMUSHORT) f;
6490 s[0] = (UEMUSHORT) f;
6491 s[1] = (UEMUSHORT) (f >> 16);
6493 /* Convert and promote the target float to E-type. */
6495 /* Output E-type to REAL_VALUE_TYPE. */
6501 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
6502 This is somewhat like ereal_unto_double, but the input types
6503 for these are different.
6505 The DFmode is stored as an array of HOST_WIDE_INT in the target's
6506 data format, with no holes in the bit packing. The first element
6507 of the input array holds the bits that would come first in the
6508 target computer's memory. */
6511 ereal_from_double (d)
6518 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
6519 if (REAL_WORDS_BIG_ENDIAN)
6521 #if HOST_BITS_PER_WIDE_INT == 32
6522 s[0] = (UEMUSHORT) (d[0] >> 16);
6523 s[1] = (UEMUSHORT) d[0];
6524 s[2] = (UEMUSHORT) (d[1] >> 16);
6525 s[3] = (UEMUSHORT) d[1];
6527 /* In this case the entire target double is contained in the
6528 first array element. The second element of the input is
6530 s[0] = (UEMUSHORT) (d[0] >> 48);
6531 s[1] = (UEMUSHORT) (d[0] >> 32);
6532 s[2] = (UEMUSHORT) (d[0] >> 16);
6533 s[3] = (UEMUSHORT) d[0];
6538 /* Target float words are little-endian. */
6539 s[0] = (UEMUSHORT) d[0];
6540 s[1] = (UEMUSHORT) (d[0] >> 16);
6541 #if HOST_BITS_PER_WIDE_INT == 32
6542 s[2] = (UEMUSHORT) d[1];
6543 s[3] = (UEMUSHORT) (d[1] >> 16);
6545 s[2] = (UEMUSHORT) (d[0] >> 32);
6546 s[3] = (UEMUSHORT) (d[0] >> 48);
6549 /* Convert target double to E-type. */
6551 /* Output E-type to REAL_VALUE_TYPE. */
6558 /* Convert target computer unsigned 64-bit integer to e-type.
6559 The endian-ness of DImode follows the convention for integers,
6560 so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN. */
6564 const UEMUSHORT *di; /* Address of the 64-bit int. */
6571 if (WORDS_BIG_ENDIAN)
6573 for (k = M; k < M + 4; k++)
6578 for (k = M + 3; k >= M; k--)
6581 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
6582 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6583 ecleaz (yi); /* it was zero */
6585 yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
6589 /* Convert target computer signed 64-bit integer to e-type. */
6593 const UEMUSHORT *di; /* Address of the 64-bit int. */
6596 unsigned EMULONG acc;
6602 if (WORDS_BIG_ENDIAN)
6604 for (k = M; k < M + 4; k++)
6609 for (k = M + 3; k >= M; k--)
6612 /* Take absolute value */
6618 for (k = M + 3; k >= M; k--)
6620 acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
6627 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
6628 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
6629 ecleaz (yi); /* it was zero */
6631 yi[E] -= (UEMUSHORT) k;/* subtract shift count from exponent */
6638 /* Convert e-type to unsigned 64-bit int. */
6654 k = (int) xi[E] - (EXONE - 1);
6657 for (j = 0; j < 4; j++)
6663 for (j = 0; j < 4; j++)
6666 warning ("overflow on truncation to integer");
6671 /* Shift more than 16 bits: first shift up k-16 mod 16,
6672 then shift up by 16's. */
6673 j = k - ((k >> 4) << 4);
6677 if (WORDS_BIG_ENDIAN)
6688 if (WORDS_BIG_ENDIAN)
6693 while ((k -= 16) > 0);
6697 /* shift not more than 16 bits */
6702 if (WORDS_BIG_ENDIAN)
6721 /* Convert e-type to signed 64-bit int. */
6728 unsigned EMULONG acc;
6735 k = (int) xi[E] - (EXONE - 1);
6738 for (j = 0; j < 4; j++)
6744 for (j = 0; j < 4; j++)
6747 warning ("overflow on truncation to integer");
6753 /* Shift more than 16 bits: first shift up k-16 mod 16,
6754 then shift up by 16's. */
6755 j = k - ((k >> 4) << 4);
6759 if (WORDS_BIG_ENDIAN)
6770 if (WORDS_BIG_ENDIAN)
6775 while ((k -= 16) > 0);
6779 /* shift not more than 16 bits */
6782 if (WORDS_BIG_ENDIAN)
6798 /* Negate if negative */
6802 if (WORDS_BIG_ENDIAN)
6804 for (k = 0; k < 4; k++)
6806 acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
6807 if (WORDS_BIG_ENDIAN)
6819 /* Longhand square root routine. */
6822 static int esqinited = 0;
6823 static unsigned short sqrndbit[NI];
6830 UEMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
6832 int i, j, k, n, nlups;
6837 sqrndbit[NI - 2] = 1;
6840 /* Check for arg <= 0 */
6841 i = ecmp (x, ezero);
6846 mtherr ("esqrt", DOMAIN);
6862 /* Bring in the arg and renormalize if it is denormal. */
6864 m = (EMULONG) xx[1]; /* local long word exponent */
6868 /* Divide exponent by 2 */
6870 exp = (unsigned short) ((m / 2) + 0x3ffe);
6872 /* Adjust if exponent odd */
6882 n = 8; /* get 8 bits of result per inner loop */
6888 /* bring in next word of arg */
6890 num[NI - 1] = xx[j + 3];
6891 /* Do additional bit on last outer loop, for roundoff. */
6894 for (i = 0; i < n; i++)
6896 /* Next 2 bits of arg */
6899 /* Shift up answer */
6901 /* Make trial divisor */
6902 for (k = 0; k < NI; k++)
6905 eaddm (sqrndbit, temp);
6906 /* Subtract and insert answer bit if it goes in */
6907 if (ecmpm (temp, num) <= 0)
6917 /* Adjust for extra, roundoff loop done. */
6918 exp += (NBITS - 1) - rndprc;
6920 /* Sticky bit = 1 if the remainder is nonzero. */
6922 for (i = 3; i < NI; i++)
6925 /* Renormalize and round off. */
6926 emdnorm (sq, k, 0, exp, !ROUND_TOWARDS_ZERO);
6930 #endif /* EMU_NON_COMPILE not defined */
6932 /* Return the binary precision of the significand for a given
6933 floating point mode. The mode can hold an integer value
6934 that many bits wide, without losing any bits. */
6937 significand_size (mode)
6938 enum machine_mode mode;
6941 /* Don't test the modes, but their sizes, lest this
6942 code won't work for BITS_PER_UNIT != 8 . */
6944 switch (GET_MODE_BITSIZE (mode))
6948 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
6955 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
6958 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
6961 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
6964 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
6977 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)