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 Free Software Foundation, Inc.
4 Contributed by Stephen L. Moshier (moshier@world.std.com).
6 This file is part of GNU CC.
8 GNU CC is free software; you can redistribute it and/or modify
9 it under the terms of the GNU General Public License as published by
10 the Free Software Foundation; either version 2, or (at your option)
13 GNU CC is distributed in the hope that it will be useful,
14 but WITHOUT ANY WARRANTY; without even the implied warranty of
15 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 GNU General Public License for more details.
18 You should have received a copy of the GNU General Public License
19 along with GNU CC; see the file COPYING. If not, write to
20 the Free Software Foundation, 59 Temple Place - Suite 330,
21 Boston, MA 02111-1307, USA. */
32 /* To enable support of XFmode extended real floating point, define
33 LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
35 To support cross compilation between IEEE, VAX and IBM floating
36 point formats, define REAL_ARITHMETIC in the tm.h file.
38 In either case the machine files (tm.h) must not contain any code
39 that tries to use host floating point arithmetic to convert
40 REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
41 etc. In cross-compile situations a REAL_VALUE_TYPE may not
42 be intelligible to the host computer's native arithmetic.
44 The emulator defaults to the host's floating point format so that
45 its decimal conversion functions can be used if desired (see
48 The first part of this file interfaces gcc to a floating point
49 arithmetic suite that was not written with gcc in mind. Avoid
50 changing the low-level arithmetic routines unless you have suitable
51 test programs available. A special version of the PARANOIA floating
52 point arithmetic tester, modified for this purpose, can be found on
53 usc.edu: /pub/C-numanal/ieeetest.zoo. Other tests, and libraries of
54 XFmode and TFmode transcendental functions, can be obtained by ftp from
55 netlib.att.com: netlib/cephes. */
57 /* Type of computer arithmetic.
58 Only one of DEC, IBM, IEEE, or UNK should get defined.
60 `IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
61 to big-endian IEEE floating-point data structure. This definition
62 should work in SFmode `float' type and DFmode `double' type on
63 virtually all big-endian IEEE machines. If LONG_DOUBLE_TYPE_SIZE
64 has been defined to be 96, then IEEE also invokes the particular
65 XFmode (`long double' type) data structure used by the Motorola
66 680x0 series processors.
68 `IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
69 little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
70 has been defined to be 96, then IEEE also invokes the particular
71 XFmode `long double' data structure used by the Intel 80x86 series
74 `DEC' refers specifically to the Digital Equipment Corp PDP-11
75 and VAX floating point data structure. This model currently
76 supports no type wider than DFmode.
78 `IBM' refers specifically to the IBM System/370 and compatible
79 floating point data structure. This model currently supports
80 no type wider than DFmode. The IBM conversions were contributed by
81 frank@atom.ansto.gov.au (Frank Crawford).
83 If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
84 then `long double' and `double' are both implemented, but they
85 both mean DFmode. In this case, the software floating-point
86 support available here is activated by writing
87 #define REAL_ARITHMETIC
90 The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
91 and may deactivate XFmode since `long double' is used to refer
94 The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
95 contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
96 separate the floating point unit's endian-ness from that of
97 the integer addressing. This permits one to define a big-endian
98 FPU on a little-endian machine (e.g., ARM). An extension to
99 BYTES_BIG_ENDIAN may be required for some machines in the future.
100 These optional macros may be defined in tm.h. In real.h, they
101 default to WORDS_BIG_ENDIAN, etc., so there is no need to define
102 them for any normal host or target machine on which the floats
103 and the integers have the same endian-ness. */
106 /* The following converts gcc macros into the ones used by this file. */
108 /* REAL_ARITHMETIC defined means that macros in real.h are
109 defined to call emulator functions. */
110 #ifdef REAL_ARITHMETIC
112 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
113 /* PDP-11, Pro350, VAX: */
115 #else /* it's not VAX */
116 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
117 /* IBM System/370 style */
119 #else /* it's also not an IBM */
120 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
122 #else /* it's not IEEE either */
123 /* UNKnown arithmetic. We don't support this and can't go on. */
124 unknown arithmetic type
126 #endif /* not IEEE */
130 #define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
133 /* REAL_ARITHMETIC not defined means that the *host's* data
134 structure will be used. It may differ by endian-ness from the
135 target machine's structure and will get its ends swapped
136 accordingly (but not here). Probably only the decimal <-> binary
137 functions in this file will actually be used in this case. */
139 #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
141 #else /* it's not VAX */
142 #if HOST_FLOAT_FORMAT == IBM_FLOAT_FORMAT
143 /* IBM System/370 style */
145 #else /* it's also not an IBM */
146 #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
148 #else /* it's not IEEE either */
149 unknown arithmetic type
151 #endif /* not IEEE */
155 #define REAL_WORDS_BIG_ENDIAN HOST_FLOAT_WORDS_BIG_ENDIAN
157 #endif /* REAL_ARITHMETIC not defined */
159 /* Define INFINITY for support of infinity.
160 Define NANS for support of Not-a-Number's (NaN's). */
161 #if !defined(DEC) && !defined(IBM)
166 /* Support of NaNs requires support of infinity. */
173 /* Find a host integer type that is at least 16 bits wide,
174 and another type at least twice whatever that size is. */
176 #if HOST_BITS_PER_CHAR >= 16
177 #define EMUSHORT char
178 #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
179 #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
181 #if HOST_BITS_PER_SHORT >= 16
182 #define EMUSHORT short
183 #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
184 #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
186 #if HOST_BITS_PER_INT >= 16
188 #define EMUSHORT_SIZE HOST_BITS_PER_INT
189 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
191 #if HOST_BITS_PER_LONG >= 16
192 #define EMUSHORT long
193 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
194 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
196 /* You will have to modify this program to have a smaller unit size. */
197 #define EMU_NON_COMPILE
203 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
204 #define EMULONG short
206 #if HOST_BITS_PER_INT >= EMULONG_SIZE
209 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
212 #if HOST_BITS_PER_LONG_LONG >= EMULONG_SIZE
213 #define EMULONG long long int
215 /* You will have to modify this program to have a smaller unit size. */
216 #define EMU_NON_COMPILE
223 /* The host interface doesn't work if no 16-bit size exists. */
224 #if EMUSHORT_SIZE != 16
225 #define EMU_NON_COMPILE
228 /* OK to continue compilation. */
229 #ifndef EMU_NON_COMPILE
231 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
232 In GET_REAL and PUT_REAL, r and e are pointers.
233 A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
234 in memory, with no holes. */
236 #if LONG_DOUBLE_TYPE_SIZE == 96
237 /* Number of 16 bit words in external e type format */
239 #define MAXDECEXP 4932
240 #define MINDECEXP -4956
241 #define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
242 #define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
243 #else /* no XFmode */
244 #if LONG_DOUBLE_TYPE_SIZE == 128
246 #define MAXDECEXP 4932
247 #define MINDECEXP -4977
248 #define GET_REAL(r,e) bcopy ((char *) r, (char *) e, 2*NE)
249 #define PUT_REAL(e,r) bcopy ((char *) e, (char *) r, 2*NE)
252 #define MAXDECEXP 4932
253 #define MINDECEXP -4956
254 #ifdef REAL_ARITHMETIC
255 /* Emulator uses target format internally
256 but host stores it in host endian-ness. */
258 #define GET_REAL(r,e) \
260 if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
261 e53toe ((unsigned EMUSHORT*) (r), (e)); \
264 unsigned EMUSHORT w[4]; \
265 w[3] = ((EMUSHORT *) r)[0]; \
266 w[2] = ((EMUSHORT *) r)[1]; \
267 w[1] = ((EMUSHORT *) r)[2]; \
268 w[0] = ((EMUSHORT *) r)[3]; \
273 #define PUT_REAL(e,r) \
275 if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
276 etoe53 ((e), (unsigned EMUSHORT *) (r)); \
279 unsigned EMUSHORT w[4]; \
281 *((EMUSHORT *) r) = w[3]; \
282 *((EMUSHORT *) r + 1) = w[2]; \
283 *((EMUSHORT *) r + 2) = w[1]; \
284 *((EMUSHORT *) r + 3) = w[0]; \
288 #else /* not REAL_ARITHMETIC */
290 /* emulator uses host format */
291 #define GET_REAL(r,e) e53toe ((unsigned EMUSHORT *) (r), (e))
292 #define PUT_REAL(e,r) etoe53 ((e), (unsigned EMUSHORT *) (r))
294 #endif /* not REAL_ARITHMETIC */
295 #endif /* not TFmode */
296 #endif /* no XFmode */
299 /* Number of 16 bit words in internal format */
302 /* Array offset to exponent */
305 /* Array offset to high guard word */
308 /* Number of bits of precision */
309 #define NBITS ((NI-4)*16)
311 /* Maximum number of decimal digits in ASCII conversion
314 #define NDEC (NBITS*8/27)
316 /* The exponent of 1.0 */
317 #define EXONE (0x3fff)
319 extern int extra_warnings;
320 extern unsigned EMUSHORT ezero[], ehalf[], eone[], etwo[];
321 extern unsigned EMUSHORT elog2[], esqrt2[];
323 static void endian PROTO((unsigned EMUSHORT *, long *,
325 static void eclear PROTO((unsigned EMUSHORT *));
326 static void emov PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
327 static void eabs PROTO((unsigned EMUSHORT *));
328 static void eneg PROTO((unsigned EMUSHORT *));
329 static int eisneg PROTO((unsigned EMUSHORT *));
330 static int eisinf PROTO((unsigned EMUSHORT *));
331 static int eisnan PROTO((unsigned EMUSHORT *));
332 static void einfin PROTO((unsigned EMUSHORT *));
333 static void enan PROTO((unsigned EMUSHORT *, int));
334 static void emovi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
335 static void emovo PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
336 static void ecleaz PROTO((unsigned EMUSHORT *));
337 static void ecleazs PROTO((unsigned EMUSHORT *));
338 static void emovz PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
339 static void einan PROTO((unsigned EMUSHORT *));
340 static int eiisnan PROTO((unsigned EMUSHORT *));
341 static int eiisneg PROTO((unsigned EMUSHORT *));
342 static void eiinfin PROTO((unsigned EMUSHORT *));
343 static int eiisinf PROTO((unsigned EMUSHORT *));
344 static int ecmpm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
345 static void eshdn1 PROTO((unsigned EMUSHORT *));
346 static void eshup1 PROTO((unsigned EMUSHORT *));
347 static void eshdn8 PROTO((unsigned EMUSHORT *));
348 static void eshup8 PROTO((unsigned EMUSHORT *));
349 static void eshup6 PROTO((unsigned EMUSHORT *));
350 static void eshdn6 PROTO((unsigned EMUSHORT *));
351 static void eaddm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
\f
352 static void esubm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
353 static void m16m PROTO((unsigned int, unsigned short *,
355 static int edivm PROTO((unsigned short *, unsigned short *));
356 static int emulm PROTO((unsigned short *, unsigned short *));
357 static void emdnorm PROTO((unsigned EMUSHORT *, int, int, EMULONG, int));
358 static void esub PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
359 unsigned EMUSHORT *));
360 static void eadd PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
361 unsigned EMUSHORT *));
362 static void eadd1 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
363 unsigned EMUSHORT *));
364 static void ediv PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
365 unsigned EMUSHORT *));
366 static void emul PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
367 unsigned EMUSHORT *));
368 static void e53toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
369 static void e64toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
370 static void e113toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
371 static void e24toe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
372 static void etoe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
373 static void toe113 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
374 static void etoe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
375 static void toe64 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
376 static void etoe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
377 static void toe53 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
378 static void etoe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
379 static void toe24 PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
380 static int ecmp PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
381 static void eround PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
382 static void ltoe PROTO((HOST_WIDE_INT *, unsigned EMUSHORT *));
383 static void ultoe PROTO((unsigned HOST_WIDE_INT *, unsigned EMUSHORT *));
384 static void eifrac PROTO((unsigned EMUSHORT *, HOST_WIDE_INT *,
385 unsigned EMUSHORT *));
386 static void euifrac PROTO((unsigned EMUSHORT *, unsigned HOST_WIDE_INT *,
387 unsigned EMUSHORT *));
388 static int eshift PROTO((unsigned EMUSHORT *, int));
389 static int enormlz PROTO((unsigned EMUSHORT *));
390 static void e24toasc PROTO((unsigned EMUSHORT *, char *, int));
391 static void e53toasc PROTO((unsigned EMUSHORT *, char *, int));
392 static void e64toasc PROTO((unsigned EMUSHORT *, char *, int));
393 static void e113toasc PROTO((unsigned EMUSHORT *, char *, int));
394 static void etoasc PROTO((unsigned EMUSHORT *, char *, int));
395 static void asctoe24 PROTO((char *, unsigned EMUSHORT *));
396 static void asctoe53 PROTO((char *, unsigned EMUSHORT *));
397 static void asctoe64 PROTO((char *, unsigned EMUSHORT *));
398 static void asctoe113 PROTO((char *, unsigned EMUSHORT *));
399 static void asctoe PROTO((char *, unsigned EMUSHORT *));
400 static void asctoeg PROTO((char *, unsigned EMUSHORT *, int));
401 static void efloor PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
402 static void efrexp PROTO((unsigned EMUSHORT *, int *,
403 unsigned EMUSHORT *));
404 static void eldexp PROTO((unsigned EMUSHORT *, int, unsigned EMUSHORT *));
405 static void eremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
406 unsigned EMUSHORT *));
407 static void eiremain PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
408 static void mtherr PROTO((char *, int));
409 static void dectoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
410 static void etodec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
411 static void todec PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
412 static void ibmtoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
414 static void etoibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
416 static void toibm PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *,
418 static void make_nan PROTO((unsigned EMUSHORT *, int, enum machine_mode));
419 static void uditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
420 static void ditoe PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
421 static void etoudi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
422 static void etodi PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
423 static void esqrt PROTO((unsigned EMUSHORT *, unsigned EMUSHORT *));
425 /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
426 swapping ends if required, into output array of longs. The
427 result is normally passed to fprintf by the ASM_OUTPUT_ macros. */
431 unsigned EMUSHORT e[];
433 enum machine_mode mode;
437 if (REAL_WORDS_BIG_ENDIAN)
443 /* Swap halfwords in the fourth long. */
444 th = (unsigned long) e[6] & 0xffff;
445 t = (unsigned long) e[7] & 0xffff;
451 /* Swap halfwords in the third long. */
452 th = (unsigned long) e[4] & 0xffff;
453 t = (unsigned long) e[5] & 0xffff;
456 /* fall into the double case */
460 /* swap halfwords in the second word */
461 th = (unsigned long) e[2] & 0xffff;
462 t = (unsigned long) e[3] & 0xffff;
465 /* fall into the float case */
470 /* swap halfwords in the first word */
471 th = (unsigned long) e[0] & 0xffff;
472 t = (unsigned long) e[1] & 0xffff;
483 /* Pack the output array without swapping. */
490 /* Pack the fourth long. */
491 th = (unsigned long) e[7] & 0xffff;
492 t = (unsigned long) e[6] & 0xffff;
498 /* Pack the third long.
499 Each element of the input REAL_VALUE_TYPE array has 16 useful bits
501 th = (unsigned long) e[5] & 0xffff;
502 t = (unsigned long) e[4] & 0xffff;
505 /* fall into the double case */
509 /* pack the second long */
510 th = (unsigned long) e[3] & 0xffff;
511 t = (unsigned long) e[2] & 0xffff;
514 /* fall into the float case */
519 /* pack the first long */
520 th = (unsigned long) e[1] & 0xffff;
521 t = (unsigned long) e[0] & 0xffff;
533 /* This is the implementation of the REAL_ARITHMETIC macro. */
536 earith (value, icode, r1, r2)
537 REAL_VALUE_TYPE *value;
542 unsigned EMUSHORT d1[NE], d2[NE], v[NE];
548 /* Return NaN input back to the caller. */
551 PUT_REAL (d1, value);
556 PUT_REAL (d2, value);
560 code = (enum tree_code) icode;
568 esub (d2, d1, v); /* d1 - d2 */
576 #ifndef REAL_INFINITY
577 if (ecmp (d2, ezero) == 0)
580 enan (v, eisneg (d1) ^ eisneg (d2));
587 ediv (d2, d1, v); /* d1/d2 */
590 case MIN_EXPR: /* min (d1,d2) */
591 if (ecmp (d1, d2) < 0)
597 case MAX_EXPR: /* max (d1,d2) */
598 if (ecmp (d1, d2) > 0)
611 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
612 implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */
618 unsigned EMUSHORT f[NE], g[NE];
634 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
635 implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */
641 unsigned EMUSHORT f[NE], g[NE];
643 unsigned HOST_WIDE_INT l;
657 /* This is the REAL_VALUE_ATOF function. It converts a decimal string to
658 binary, rounding off as indicated by the machine_mode argument. Then it
659 promotes the rounded value to REAL_VALUE_TYPE. */
666 unsigned EMUSHORT tem[NE], e[NE];
696 /* Expansion of REAL_NEGATE. */
702 unsigned EMUSHORT e[NE];
712 /* Round real toward zero to HOST_WIDE_INT;
713 implements REAL_VALUE_FIX (x). */
719 unsigned EMUSHORT f[NE], g[NE];
726 warning ("conversion from NaN to int");
734 /* Round real toward zero to unsigned HOST_WIDE_INT
735 implements REAL_VALUE_UNSIGNED_FIX (x).
736 Negative input returns zero. */
738 unsigned HOST_WIDE_INT
742 unsigned EMUSHORT f[NE], g[NE];
743 unsigned HOST_WIDE_INT l;
749 warning ("conversion from NaN to unsigned int");
758 /* REAL_VALUE_FROM_INT macro. */
761 ereal_from_int (d, i, j, mode)
764 enum machine_mode mode;
766 unsigned EMUSHORT df[NE], dg[NE];
767 HOST_WIDE_INT low, high;
770 if (GET_MODE_CLASS (mode) != MODE_FLOAT)
777 /* complement and add 1 */
784 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
785 ultoe ((unsigned HOST_WIDE_INT *) &high, dg);
787 ultoe ((unsigned HOST_WIDE_INT *) &low, df);
792 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
793 Avoid double-rounding errors later by rounding off now from the
794 extra-wide internal format to the requested precision. */
795 switch (GET_MODE_BITSIZE (mode))
825 /* REAL_VALUE_FROM_UNSIGNED_INT macro. */
828 ereal_from_uint (d, i, j, mode)
830 unsigned HOST_WIDE_INT i, j;
831 enum machine_mode mode;
833 unsigned EMUSHORT df[NE], dg[NE];
834 unsigned HOST_WIDE_INT low, high;
836 if (GET_MODE_CLASS (mode) != MODE_FLOAT)
840 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
846 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
847 Avoid double-rounding errors later by rounding off now from the
848 extra-wide internal format to the requested precision. */
849 switch (GET_MODE_BITSIZE (mode))
879 /* REAL_VALUE_TO_INT macro. */
882 ereal_to_int (low, high, rr)
883 HOST_WIDE_INT *low, *high;
886 unsigned EMUSHORT d[NE], df[NE], dg[NE], dh[NE];
893 warning ("conversion from NaN to int");
899 /* convert positive value */
906 eldexp (eone, HOST_BITS_PER_WIDE_INT, df);
907 ediv (df, d, dg); /* dg = d / 2^32 is the high word */
908 euifrac (dg, (unsigned HOST_WIDE_INT *) high, dh);
909 emul (df, dh, dg); /* fractional part is the low word */
910 euifrac (dg, (unsigned HOST_WIDE_INT *)low, dh);
913 /* complement and add 1 */
923 /* REAL_VALUE_LDEXP macro. */
930 unsigned EMUSHORT e[NE], y[NE];
943 /* These routines are conditionally compiled because functions
944 of the same names may be defined in fold-const.c. */
946 #ifdef REAL_ARITHMETIC
948 /* Check for infinity in a REAL_VALUE_TYPE. */
954 unsigned EMUSHORT e[NE];
964 /* Check whether a REAL_VALUE_TYPE item is a NaN. */
970 unsigned EMUSHORT e[NE];
981 /* Check for a negative REAL_VALUE_TYPE number.
982 This just checks the sign bit, so that -0 counts as negative. */
988 return ereal_isneg (x);
991 /* Expansion of REAL_VALUE_TRUNCATE.
992 The result is in floating point, rounded to nearest or even. */
995 real_value_truncate (mode, arg)
996 enum machine_mode mode;
999 unsigned EMUSHORT e[NE], t[NE];
1035 /* If an unsupported type was requested, presume that
1036 the machine files know something useful to do with
1037 the unmodified value. */
1046 #endif /* REAL_ARITHMETIC defined */
1048 /* Used for debugging--print the value of R in human-readable format
1057 REAL_VALUE_TO_DECIMAL (r, "%.20g", dstr);
1058 fprintf (stderr, "%s", dstr);
1062 /* The following routines convert REAL_VALUE_TYPE to the various floating
1063 point formats that are meaningful to supported computers.
1065 The results are returned in 32-bit pieces, each piece stored in a `long'.
1066 This is so they can be printed by statements like
1068 fprintf (file, "%lx, %lx", L[0], L[1]);
1070 that will work on both narrow- and wide-word host computers. */
1072 /* Convert R to a 128-bit long double precision value. The output array L
1073 contains four 32-bit pieces of the result, in the order they would appear
1081 unsigned EMUSHORT e[NE];
1085 endian (e, l, TFmode);
1088 /* Convert R to a double extended precision value. The output array L
1089 contains three 32-bit pieces of the result, in the order they would
1090 appear in memory. */
1097 unsigned EMUSHORT e[NE];
1101 endian (e, l, XFmode);
1104 /* Convert R to a double precision value. The output array L contains two
1105 32-bit pieces of the result, in the order they would appear in memory. */
1112 unsigned EMUSHORT e[NE];
1116 endian (e, l, DFmode);
1119 /* Convert R to a single precision float value stored in the least-significant
1120 bits of a `long'. */
1126 unsigned EMUSHORT e[NE];
1131 endian (e, &l, SFmode);
1135 /* Convert X to a decimal ASCII string S for output to an assembly
1136 language file. Note, there is no standard way to spell infinity or
1137 a NaN, so these values may require special treatment in the tm.h
1141 ereal_to_decimal (x, s)
1145 unsigned EMUSHORT e[NE];
1151 /* Compare X and Y. Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
1152 or -2 if either is a NaN. */
1156 REAL_VALUE_TYPE x, y;
1158 unsigned EMUSHORT ex[NE], ey[NE];
1162 return (ecmp (ex, ey));
1165 /* Return 1 if the sign bit of X is set, else return 0. */
1171 unsigned EMUSHORT ex[NE];
1174 return (eisneg (ex));
1177 /* End of REAL_ARITHMETIC interface */
1180 Extended precision IEEE binary floating point arithmetic routines
1182 Numbers are stored in C language as arrays of 16-bit unsigned
1183 short integers. The arguments of the routines are pointers to
1186 External e type data structure, similar to Intel 8087 chip
1187 temporary real format but possibly with a larger significand:
1189 NE-1 significand words (least significant word first,
1190 most significant bit is normally set)
1191 exponent (value = EXONE for 1.0,
1192 top bit is the sign)
1195 Internal exploded e-type data structure of a number (a "word" is 16 bits):
1197 ei[0] sign word (0 for positive, 0xffff for negative)
1198 ei[1] biased exponent (value = EXONE for the number 1.0)
1199 ei[2] high guard word (always zero after normalization)
1201 to ei[NI-2] significand (NI-4 significand words,
1202 most significant word first,
1203 most significant bit is set)
1204 ei[NI-1] low guard word (0x8000 bit is rounding place)
1208 Routines for external format e-type numbers
1210 asctoe (string, e) ASCII string to extended double e type
1211 asctoe64 (string, &d) ASCII string to long double
1212 asctoe53 (string, &d) ASCII string to double
1213 asctoe24 (string, &f) ASCII string to single
1214 asctoeg (string, e, prec) ASCII string to specified precision
1215 e24toe (&f, e) IEEE single precision to e type
1216 e53toe (&d, e) IEEE double precision to e type
1217 e64toe (&d, e) IEEE long double precision to e type
1218 e113toe (&d, e) 128-bit long double precision to e type
1219 eabs (e) absolute value
1220 eadd (a, b, c) c = b + a
1222 ecmp (a, b) Returns 1 if a > b, 0 if a == b,
1223 -1 if a < b, -2 if either a or b is a NaN.
1224 ediv (a, b, c) c = b / a
1225 efloor (a, b) truncate to integer, toward -infinity
1226 efrexp (a, exp, s) extract exponent and significand
1227 eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
1228 euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
1229 einfin (e) set e to infinity, leaving its sign alone
1230 eldexp (a, n, b) multiply by 2**n
1232 emul (a, b, c) c = b * a
1234 eround (a, b) b = nearest integer value to a
1235 esub (a, b, c) c = b - a
1236 e24toasc (&f, str, n) single to ASCII string, n digits after decimal
1237 e53toasc (&d, str, n) double to ASCII string, n digits after decimal
1238 e64toasc (&d, str, n) 80-bit long double to ASCII string
1239 e113toasc (&d, str, n) 128-bit long double to ASCII string
1240 etoasc (e, str, n) e to ASCII string, n digits after decimal
1241 etoe24 (e, &f) convert e type to IEEE single precision
1242 etoe53 (e, &d) convert e type to IEEE double precision
1243 etoe64 (e, &d) convert e type to IEEE long double precision
1244 ltoe (&l, e) HOST_WIDE_INT to e type
1245 ultoe (&l, e) unsigned HOST_WIDE_INT to e type
1246 eisneg (e) 1 if sign bit of e != 0, else 0
1247 eisinf (e) 1 if e has maximum exponent (non-IEEE)
1248 or is infinite (IEEE)
1249 eisnan (e) 1 if e is a NaN
1252 Routines for internal format exploded e-type numbers
1254 eaddm (ai, bi) add significands, bi = bi + ai
1256 ecleazs (ei) set ei = 0 but leave its sign alone
1257 ecmpm (ai, bi) compare significands, return 1, 0, or -1
1258 edivm (ai, bi) divide significands, bi = bi / ai
1259 emdnorm (ai,l,s,exp) normalize and round off
1260 emovi (a, ai) convert external a to internal ai
1261 emovo (ai, a) convert internal ai to external a
1262 emovz (ai, bi) bi = ai, low guard word of bi = 0
1263 emulm (ai, bi) multiply significands, bi = bi * ai
1264 enormlz (ei) left-justify the significand
1265 eshdn1 (ai) shift significand and guards down 1 bit
1266 eshdn8 (ai) shift down 8 bits
1267 eshdn6 (ai) shift down 16 bits
1268 eshift (ai, n) shift ai n bits up (or down if n < 0)
1269 eshup1 (ai) shift significand and guards up 1 bit
1270 eshup8 (ai) shift up 8 bits
1271 eshup6 (ai) shift up 16 bits
1272 esubm (ai, bi) subtract significands, bi = bi - ai
1273 eiisinf (ai) 1 if infinite
1274 eiisnan (ai) 1 if a NaN
1275 eiisneg (ai) 1 if sign bit of ai != 0, else 0
1276 einan (ai) set ai = NaN
1277 eiinfin (ai) set ai = infinity
1279 The result is always normalized and rounded to NI-4 word precision
1280 after each arithmetic operation.
1282 Exception flags are NOT fully supported.
1284 Signaling NaN's are NOT supported; they are treated the same
1287 Define INFINITY for support of infinity; otherwise a
1288 saturation arithmetic is implemented.
1290 Define NANS for support of Not-a-Number items; otherwise the
1291 arithmetic will never produce a NaN output, and might be confused
1293 If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1294 either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1295 may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1298 Denormals are always supported here where appropriate (e.g., not
1299 for conversion to DEC numbers). */
1301 /* Definitions for error codes that are passed to the common error handling
1304 For Digital Equipment PDP-11 and VAX computers, certain
1305 IBM systems, and others that use numbers with a 56-bit
1306 significand, the symbol DEC should be defined. In this
1307 mode, most floating point constants are given as arrays
1308 of octal integers to eliminate decimal to binary conversion
1309 errors that might be introduced by the compiler.
1311 For computers, such as IBM PC, that follow the IEEE
1312 Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1313 Std 754-1985), the symbol IEEE should be defined.
1314 These numbers have 53-bit significands. In this mode, constants
1315 are provided as arrays of hexadecimal 16 bit integers.
1316 The endian-ness of generated values is controlled by
1317 REAL_WORDS_BIG_ENDIAN.
1319 To accommodate other types of computer arithmetic, all
1320 constants are also provided in a normal decimal radix
1321 which one can hope are correctly converted to a suitable
1322 format by the available C language compiler. To invoke
1323 this mode, the symbol UNK is defined.
1325 An important difference among these modes is a predefined
1326 set of machine arithmetic constants for each. The numbers
1327 MACHEP (the machine roundoff error), MAXNUM (largest number
1328 represented), and several other parameters are preset by
1329 the configuration symbol. Check the file const.c to
1330 ensure that these values are correct for your computer.
1332 For ANSI C compatibility, define ANSIC equal to 1. Currently
1333 this affects only the atan2 function and others that use it. */
1335 /* Constant definitions for math error conditions. */
1337 #define DOMAIN 1 /* argument domain error */
1338 #define SING 2 /* argument singularity */
1339 #define OVERFLOW 3 /* overflow range error */
1340 #define UNDERFLOW 4 /* underflow range error */
1341 #define TLOSS 5 /* total loss of precision */
1342 #define PLOSS 6 /* partial loss of precision */
1343 #define INVALID 7 /* NaN-producing operation */
1345 /* e type constants used by high precision check routines */
1347 #if LONG_DOUBLE_TYPE_SIZE == 128
1349 unsigned EMUSHORT ezero[NE] =
1350 {0x0000, 0x0000, 0x0000, 0x0000,
1351 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
1352 extern unsigned EMUSHORT ezero[];
1355 unsigned EMUSHORT ehalf[NE] =
1356 {0x0000, 0x0000, 0x0000, 0x0000,
1357 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
1358 extern unsigned EMUSHORT ehalf[];
1361 unsigned EMUSHORT eone[NE] =
1362 {0x0000, 0x0000, 0x0000, 0x0000,
1363 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
1364 extern unsigned EMUSHORT eone[];
1367 unsigned EMUSHORT etwo[NE] =
1368 {0x0000, 0x0000, 0x0000, 0x0000,
1369 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
1370 extern unsigned EMUSHORT etwo[];
1373 unsigned EMUSHORT e32[NE] =
1374 {0x0000, 0x0000, 0x0000, 0x0000,
1375 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
1376 extern unsigned EMUSHORT e32[];
1378 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1379 unsigned EMUSHORT elog2[NE] =
1380 {0x40f3, 0xf6af, 0x03f2, 0xb398,
1381 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1382 extern unsigned EMUSHORT elog2[];
1384 /* 1.41421356237309504880168872420969807856967187537695E0 */
1385 unsigned EMUSHORT esqrt2[NE] =
1386 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1387 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1388 extern unsigned EMUSHORT esqrt2[];
1390 /* 3.14159265358979323846264338327950288419716939937511E0 */
1391 unsigned EMUSHORT epi[NE] =
1392 {0x2902, 0x1cd1, 0x80dc, 0x628b,
1393 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1394 extern unsigned EMUSHORT epi[];
1397 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1398 unsigned EMUSHORT ezero[NE] =
1399 {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1400 unsigned EMUSHORT ehalf[NE] =
1401 {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1402 unsigned EMUSHORT eone[NE] =
1403 {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1404 unsigned EMUSHORT etwo[NE] =
1405 {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1406 unsigned EMUSHORT e32[NE] =
1407 {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1408 unsigned EMUSHORT elog2[NE] =
1409 {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1410 unsigned EMUSHORT esqrt2[NE] =
1411 {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1412 unsigned EMUSHORT epi[NE] =
1413 {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1416 /* Control register for rounding precision.
1417 This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits. */
1422 /* Clear out entire e-type number X. */
1426 register unsigned EMUSHORT *x;
1430 for (i = 0; i < NE; i++)
1434 /* Move e-type number from A to B. */
1438 register unsigned EMUSHORT *a, *b;
1442 for (i = 0; i < NE; i++)
1447 /* Absolute value of e-type X. */
1451 unsigned EMUSHORT x[];
1453 /* sign is top bit of last word of external format */
1454 x[NE - 1] &= 0x7fff;
1457 /* Negate the e-type number X. */
1461 unsigned EMUSHORT x[];
1464 x[NE - 1] ^= 0x8000; /* Toggle the sign bit */
1467 /* Return 1 if sign bit of e-type number X is nonzero, else zero. */
1471 unsigned EMUSHORT x[];
1474 if (x[NE - 1] & 0x8000)
1480 /* Return 1 if e-type number X is infinity, else return zero. */
1484 unsigned EMUSHORT x[];
1491 if ((x[NE - 1] & 0x7fff) == 0x7fff)
1497 /* Check if e-type number is not a number. The bit pattern is one that we
1498 defined, so we know for sure how to detect it. */
1502 unsigned EMUSHORT x[];
1507 /* NaN has maximum exponent */
1508 if ((x[NE - 1] & 0x7fff) != 0x7fff)
1510 /* ... and non-zero significand field. */
1511 for (i = 0; i < NE - 1; i++)
1521 /* Fill e-type number X with infinity pattern (IEEE)
1522 or largest possible number (non-IEEE). */
1526 register unsigned EMUSHORT *x;
1531 for (i = 0; i < NE - 1; i++)
1535 for (i = 0; i < NE - 1; i++)
1563 /* Output an e-type NaN.
1564 This generates Intel's quiet NaN pattern for extended real.
1565 The exponent is 7fff, the leading mantissa word is c000. */
1569 register unsigned EMUSHORT *x;
1574 for (i = 0; i < NE - 2; i++)
1577 *x = (sign << 15) | 0x7fff;
1580 /* Move in an e-type number A, converting it to exploded e-type B. */
1584 unsigned EMUSHORT *a, *b;
1586 register unsigned EMUSHORT *p, *q;
1590 p = a + (NE - 1); /* point to last word of external number */
1591 /* get the sign bit */
1596 /* get the exponent */
1598 *q++ &= 0x7fff; /* delete the sign bit */
1600 if ((*(q - 1) & 0x7fff) == 0x7fff)
1606 for (i = 3; i < NI; i++)
1612 for (i = 2; i < NI; i++)
1618 /* clear high guard word */
1620 /* move in the significand */
1621 for (i = 0; i < NE - 1; i++)
1623 /* clear low guard word */
1627 /* Move out exploded e-type number A, converting it to e type B. */
1631 unsigned EMUSHORT *a, *b;
1633 register unsigned EMUSHORT *p, *q;
1634 unsigned EMUSHORT i;
1638 q = b + (NE - 1); /* point to output exponent */
1639 /* combine sign and exponent */
1642 *q-- = *p++ | 0x8000;
1646 if (*(p - 1) == 0x7fff)
1651 enan (b, eiisneg (a));
1659 /* skip over guard word */
1661 /* move the significand */
1662 for (j = 0; j < NE - 1; j++)
1666 /* Clear out exploded e-type number XI. */
1670 register unsigned EMUSHORT *xi;
1674 for (i = 0; i < NI; i++)
1678 /* Clear out exploded e-type XI, but don't touch the sign. */
1682 register unsigned EMUSHORT *xi;
1687 for (i = 0; i < NI - 1; i++)
1691 /* Move exploded e-type number from A to B. */
1695 register unsigned EMUSHORT *a, *b;
1699 for (i = 0; i < NI - 1; i++)
1701 /* clear low guard word */
1705 /* Generate exploded e-type NaN.
1706 The explicit pattern for this is maximum exponent and
1707 top two significant bits set. */
1711 unsigned EMUSHORT x[];
1719 /* Return nonzero if exploded e-type X is a NaN. */
1723 unsigned EMUSHORT x[];
1727 if ((x[E] & 0x7fff) == 0x7fff)
1729 for (i = M + 1; i < NI; i++)
1738 /* Return nonzero if sign of exploded e-type X is nonzero. */
1742 unsigned EMUSHORT x[];
1748 /* Fill exploded e-type X with infinity pattern.
1749 This has maximum exponent and significand all zeros. */
1753 unsigned EMUSHORT x[];
1760 /* Return nonzero if exploded e-type X is infinite. */
1764 unsigned EMUSHORT x[];
1771 if ((x[E] & 0x7fff) == 0x7fff)
1777 /* Compare significands of numbers in internal exploded e-type format.
1778 Guard words are included in the comparison.
1786 register unsigned EMUSHORT *a, *b;
1790 a += M; /* skip up to significand area */
1792 for (i = M; i < NI; i++)
1800 if (*(--a) > *(--b))
1806 /* Shift significand of exploded e-type X down by 1 bit. */
1810 register unsigned EMUSHORT *x;
1812 register unsigned EMUSHORT bits;
1815 x += M; /* point to significand area */
1818 for (i = M; i < NI; i++)
1830 /* Shift significand of exploded e-type X up by 1 bit. */
1834 register unsigned EMUSHORT *x;
1836 register unsigned EMUSHORT bits;
1842 for (i = M; i < NI; i++)
1855 /* Shift significand of exploded e-type X down by 8 bits. */
1859 register unsigned EMUSHORT *x;
1861 register unsigned EMUSHORT newbyt, oldbyt;
1866 for (i = M; i < NI; i++)
1876 /* Shift significand of exploded e-type X up by 8 bits. */
1880 register unsigned EMUSHORT *x;
1883 register unsigned EMUSHORT newbyt, oldbyt;
1888 for (i = M; i < NI; i++)
1898 /* Shift significand of exploded e-type X up by 16 bits. */
1902 register unsigned EMUSHORT *x;
1905 register unsigned EMUSHORT *p;
1910 for (i = M; i < NI - 1; i++)
1916 /* Shift significand of exploded e-type X down by 16 bits. */
1920 register unsigned EMUSHORT *x;
1923 register unsigned EMUSHORT *p;
1928 for (i = M; i < NI - 1; i++)
1934 /* Add significands of exploded e-type X and Y. X + Y replaces Y. */
1938 unsigned EMUSHORT *x, *y;
1940 register unsigned EMULONG a;
1947 for (i = M; i < NI; i++)
1949 a = (unsigned EMULONG) (*x) + (unsigned EMULONG) (*y) + carry;
1954 *y = (unsigned EMUSHORT) a;
1960 /* Subtract significands of exploded e-type X and Y. Y - X replaces Y. */
1964 unsigned EMUSHORT *x, *y;
1973 for (i = M; i < NI; i++)
1975 a = (unsigned EMULONG) (*y) - (unsigned EMULONG) (*x) - carry;
1980 *y = (unsigned EMUSHORT) a;
1987 static unsigned EMUSHORT equot[NI];
1991 /* Radix 2 shift-and-add versions of multiply and divide */
1994 /* Divide significands */
1998 unsigned EMUSHORT den[], num[];
2001 register unsigned EMUSHORT *p, *q;
2002 unsigned EMUSHORT j;
2008 for (i = M; i < NI; i++)
2013 /* Use faster compare and subtraction if denominator has only 15 bits of
2019 for (i = M + 3; i < NI; i++)
2024 if ((den[M + 1] & 1) != 0)
2032 for (i = 0; i < NBITS + 2; i++)
2050 /* The number of quotient bits to calculate is NBITS + 1 scaling guard
2051 bit + 1 roundoff bit. */
2056 for (i = 0; i < NBITS + 2; i++)
2058 if (ecmpm (den, num) <= 0)
2061 j = 1; /* quotient bit = 1 */
2075 /* test for nonzero remainder after roundoff bit */
2078 for (i = M; i < NI; i++)
2086 for (i = 0; i < NI; i++)
2092 /* Multiply significands */
2095 unsigned EMUSHORT a[], b[];
2097 unsigned EMUSHORT *p, *q;
2102 for (i = M; i < NI; i++)
2107 while (*p == 0) /* significand is not supposed to be zero */
2112 if ((*p & 0xff) == 0)
2120 for (i = 0; i < k; i++)
2124 /* remember if there were any nonzero bits shifted out */
2131 for (i = 0; i < NI; i++)
2134 /* return flag for lost nonzero bits */
2140 /* Radix 65536 versions of multiply and divide. */
2142 /* Multiply significand of e-type number B
2143 by 16-bit quantity A, return e-type result to C. */
2148 unsigned EMUSHORT b[], c[];
2150 register unsigned EMUSHORT *pp;
2151 register unsigned EMULONG carry;
2152 unsigned EMUSHORT *ps;
2153 unsigned EMUSHORT p[NI];
2154 unsigned EMULONG aa, m;
2163 for (i=M+1; i<NI; i++)
2173 m = (unsigned EMULONG) aa * *ps--;
2174 carry = (m & 0xffff) + *pp;
2175 *pp-- = (unsigned EMUSHORT)carry;
2176 carry = (carry >> 16) + (m >> 16) + *pp;
2177 *pp = (unsigned EMUSHORT)carry;
2178 *(pp-1) = carry >> 16;
2181 for (i=M; i<NI; i++)
2185 /* Divide significands of exploded e-types NUM / DEN. Neither the
2186 numerator NUM nor the denominator DEN is permitted to have its high guard
2191 unsigned EMUSHORT den[], num[];
2194 register unsigned EMUSHORT *p;
2195 unsigned EMULONG tnum;
2196 unsigned EMUSHORT j, tdenm, tquot;
2197 unsigned EMUSHORT tprod[NI+1];
2203 for (i=M; i<NI; i++)
2209 for (i=M; i<NI; i++)
2211 /* Find trial quotient digit (the radix is 65536). */
2212 tnum = (((unsigned EMULONG) num[M]) << 16) + num[M+1];
2214 /* Do not execute the divide instruction if it will overflow. */
2215 if ((tdenm * 0xffffL) < tnum)
2218 tquot = tnum / tdenm;
2219 /* Multiply denominator by trial quotient digit. */
2220 m16m ((unsigned int)tquot, den, tprod);
2221 /* The quotient digit may have been overestimated. */
2222 if (ecmpm (tprod, num) > 0)
2226 if (ecmpm (tprod, num) > 0)
2236 /* test for nonzero remainder after roundoff bit */
2239 for (i=M; i<NI; i++)
2246 for (i=0; i<NI; i++)
2252 /* Multiply significands of exploded e-type A and B, result in B. */
2256 unsigned EMUSHORT a[], b[];
2258 unsigned EMUSHORT *p, *q;
2259 unsigned EMUSHORT pprod[NI];
2260 unsigned EMUSHORT j;
2265 for (i=M; i<NI; i++)
2271 for (i=M+1; i<NI; i++)
2279 m16m ((unsigned int) *p--, b, pprod);
2280 eaddm(pprod, equot);
2286 for (i=0; i<NI; i++)
2289 /* return flag for lost nonzero bits */
2295 /* Normalize and round off.
2297 The internal format number to be rounded is S.
2298 Input LOST is 0 if the value is exact. This is the so-called sticky bit.
2300 Input SUBFLG indicates whether the number was obtained
2301 by a subtraction operation. In that case if LOST is nonzero
2302 then the number is slightly smaller than indicated.
2304 Input EXP is the biased exponent, which may be negative.
2305 the exponent field of S is ignored but is replaced by
2306 EXP as adjusted by normalization and rounding.
2308 Input RCNTRL is the rounding control. If it is nonzero, the
2309 returned value will be rounded to RNDPRC bits.
2311 For future reference: In order for emdnorm to round off denormal
2312 significands at the right point, the input exponent must be
2313 adjusted to be the actual value it would have after conversion to
2314 the final floating point type. This adjustment has been
2315 implemented for all type conversions (etoe53, etc.) and decimal
2316 conversions, but not for the arithmetic functions (eadd, etc.).
2317 Data types having standard 15-bit exponents are not affected by
2318 this, but SFmode and DFmode are affected. For example, ediv with
2319 rndprc = 24 will not round correctly to 24-bit precision if the
2320 result is denormal. */
2322 static int rlast = -1;
2324 static unsigned EMUSHORT rmsk = 0;
2325 static unsigned EMUSHORT rmbit = 0;
2326 static unsigned EMUSHORT rebit = 0;
2328 static unsigned EMUSHORT rbit[NI];
2331 emdnorm (s, lost, subflg, exp, rcntrl)
2332 unsigned EMUSHORT s[];
2339 unsigned EMUSHORT r;
2344 /* a blank significand could mean either zero or infinity. */
2357 if ((j > NBITS) && (exp < 32767))
2365 if (exp > (EMULONG) (-NBITS - 1))
2378 /* Round off, unless told not to by rcntrl. */
2381 /* Set up rounding parameters if the control register changed. */
2382 if (rndprc != rlast)
2389 rw = NI - 1; /* low guard word */
2409 /* For DEC or IBM arithmetic */
2436 /* Shift down 1 temporarily if the data structure has an implied
2437 most significant bit and the number is denormal.
2438 Intel long double denormals also lose one bit of precision. */
2439 if ((exp <= 0) && (rndprc != NBITS)
2440 && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2442 lost |= s[NI - 1] & 1;
2445 /* Clear out all bits below the rounding bit,
2446 remembering in r if any were nonzero. */
2460 if ((r & rmbit) != 0)
2465 { /* round to even */
2466 if ((s[re] & rebit) == 0)
2478 /* Undo the temporary shift for denormal values. */
2479 if ((exp <= 0) && (rndprc != NBITS)
2480 && ((rndprc != 64) || ((rndprc == 64) && ! REAL_WORDS_BIG_ENDIAN)))
2485 { /* overflow on roundoff */
2498 for (i = 2; i < NI - 1; i++)
2501 warning ("floating point overflow");
2505 for (i = M + 1; i < NI - 1; i++)
2508 if ((rndprc < 64) || (rndprc == 113))
2523 s[1] = (unsigned EMUSHORT) exp;
2526 /* Subtract. C = B - A, all e type numbers. */
2528 static int subflg = 0;
2532 unsigned EMUSHORT *a, *b, *c;
2546 /* Infinity minus infinity is a NaN.
2547 Test for subtracting infinities of the same sign. */
2548 if (eisinf (a) && eisinf (b)
2549 && ((eisneg (a) ^ eisneg (b)) == 0))
2551 mtherr ("esub", INVALID);
2560 /* Add. C = A + B, all e type. */
2564 unsigned EMUSHORT *a, *b, *c;
2568 /* NaN plus anything is a NaN. */
2579 /* Infinity minus infinity is a NaN.
2580 Test for adding infinities of opposite signs. */
2581 if (eisinf (a) && eisinf (b)
2582 && ((eisneg (a) ^ eisneg (b)) != 0))
2584 mtherr ("esub", INVALID);
2593 /* Arithmetic common to both addition and subtraction. */
2597 unsigned EMUSHORT *a, *b, *c;
2599 unsigned EMUSHORT ai[NI], bi[NI], ci[NI];
2601 EMULONG lt, lta, ltb;
2622 /* compare exponents */
2627 { /* put the larger number in bi */
2637 if (lt < (EMULONG) (-NBITS - 1))
2638 goto done; /* answer same as larger addend */
2640 lost = eshift (ai, k); /* shift the smaller number down */
2644 /* exponents were the same, so must compare significands */
2647 { /* the numbers are identical in magnitude */
2648 /* if different signs, result is zero */
2654 /* if same sign, result is double */
2655 /* double denormalized tiny number */
2656 if ((bi[E] == 0) && ((bi[3] & 0x8000) == 0))
2661 /* add 1 to exponent unless both are zero! */
2662 for (j = 1; j < NI - 1; j++)
2678 bi[E] = (unsigned EMUSHORT) ltb;
2682 { /* put the larger number in bi */
2698 emdnorm (bi, lost, subflg, ltb, 64);
2704 /* Divide: C = B/A, all e type. */
2708 unsigned EMUSHORT *a, *b, *c;
2710 unsigned EMUSHORT ai[NI], bi[NI];
2712 EMULONG lt, lta, ltb;
2714 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2715 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
2716 sign = eisneg(a) ^ eisneg(b);
2719 /* Return any NaN input. */
2730 /* Zero over zero, or infinity over infinity, is a NaN. */
2731 if (((ecmp (a, ezero) == 0) && (ecmp (b, ezero) == 0))
2732 || (eisinf (a) && eisinf (b)))
2734 mtherr ("ediv", INVALID);
2739 /* Infinity over anything else is infinity. */
2746 /* Anything else over infinity is zero. */
2758 { /* See if numerator is zero. */
2759 for (i = 1; i < NI - 1; i++)
2763 ltb -= enormlz (bi);
2773 { /* possible divide by zero */
2774 for (i = 1; i < NI - 1; i++)
2778 lta -= enormlz (ai);
2782 /* Divide by zero is not an invalid operation.
2783 It is a divide-by-zero operation! */
2785 mtherr ("ediv", SING);
2791 /* calculate exponent */
2792 lt = ltb - lta + EXONE;
2793 emdnorm (bi, i, 0, lt, 64);
2800 && (ecmp (c, ezero) != 0)
2803 *(c+(NE-1)) |= 0x8000;
2805 *(c+(NE-1)) &= ~0x8000;
2808 /* Multiply e-types A and B, return e-type product C. */
2812 unsigned EMUSHORT *a, *b, *c;
2814 unsigned EMUSHORT ai[NI], bi[NI];
2816 EMULONG lt, lta, ltb;
2818 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2819 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
2820 sign = eisneg(a) ^ eisneg(b);
2823 /* NaN times anything is the same NaN. */
2834 /* Zero times infinity is a NaN. */
2835 if ((eisinf (a) && (ecmp (b, ezero) == 0))
2836 || (eisinf (b) && (ecmp (a, ezero) == 0)))
2838 mtherr ("emul", INVALID);
2843 /* Infinity times anything else is infinity. */
2845 if (eisinf (a) || eisinf (b))
2857 for (i = 1; i < NI - 1; i++)
2861 lta -= enormlz (ai);
2872 for (i = 1; i < NI - 1; i++)
2876 ltb -= enormlz (bi);
2885 /* Multiply significands */
2887 /* calculate exponent */
2888 lt = lta + ltb - (EXONE - 1);
2889 emdnorm (bi, j, 0, lt, 64);
2896 && (ecmp (c, ezero) != 0)
2899 *(c+(NE-1)) |= 0x8000;
2901 *(c+(NE-1)) &= ~0x8000;
2904 /* Convert double precision PE to e-type Y. */
2908 unsigned EMUSHORT *pe, *y;
2917 ibmtoe (pe, y, DFmode);
2920 register unsigned EMUSHORT r;
2921 register unsigned EMUSHORT *e, *p;
2922 unsigned EMUSHORT yy[NI];
2926 denorm = 0; /* flag if denormalized number */
2928 if (! REAL_WORDS_BIG_ENDIAN)
2934 yy[M] = (r & 0x0f) | 0x10;
2935 r &= ~0x800f; /* strip sign and 4 significand bits */
2940 if (! REAL_WORDS_BIG_ENDIAN)
2942 if (((pe[3] & 0xf) != 0) || (pe[2] != 0)
2943 || (pe[1] != 0) || (pe[0] != 0))
2945 enan (y, yy[0] != 0);
2951 if (((pe[0] & 0xf) != 0) || (pe[1] != 0)
2952 || (pe[2] != 0) || (pe[3] != 0))
2954 enan (y, yy[0] != 0);
2965 #endif /* INFINITY */
2967 /* If zero exponent, then the significand is denormalized.
2968 So take back the understood high significand bit. */
2979 if (! REAL_WORDS_BIG_ENDIAN)
2995 { /* if zero exponent, then normalize the significand */
2996 if ((k = enormlz (yy)) > NBITS)
2999 yy[E] -= (unsigned EMUSHORT) (k - 1);
3002 #endif /* not IBM */
3003 #endif /* not DEC */
3006 /* Convert double extended precision float PE to e type Y. */
3010 unsigned EMUSHORT *pe, *y;
3012 unsigned EMUSHORT yy[NI];
3013 unsigned EMUSHORT *e, *p, *q;
3018 for (i = 0; i < NE - 5; i++)
3020 /* This precision is not ordinarily supported on DEC or IBM. */
3022 for (i = 0; i < 5; i++)
3026 p = &yy[0] + (NE - 1);
3029 for (i = 0; i < 5; i++)
3033 if (! REAL_WORDS_BIG_ENDIAN)
3035 for (i = 0; i < 5; i++)
3038 /* For denormal long double Intel format, shift significand up one
3039 -- but only if the top significand bit is zero. A top bit of 1
3040 is "pseudodenormal" when the exponent is zero. */
3041 if((yy[NE-1] & 0x7fff) == 0 && (yy[NE-2] & 0x8000) == 0)
3043 unsigned EMUSHORT temp[NI];
3053 p = &yy[0] + (NE - 1);
3054 #ifdef ARM_EXTENDED_IEEE_FORMAT
3055 /* For ARMs, the exponent is in the lowest 15 bits of the word. */
3056 *p-- = (e[0] & 0x8000) | (e[1] & 0x7ffff);
3062 for (i = 0; i < 4; i++)
3067 /* Point to the exponent field and check max exponent cases. */
3069 if ((*p & 0x7fff) == 0x7fff)
3072 if (! REAL_WORDS_BIG_ENDIAN)
3074 for (i = 0; i < 4; i++)
3076 if ((i != 3 && pe[i] != 0)
3077 /* Anything but 0x8000 here, including 0, is a NaN. */
3078 || (i == 3 && pe[i] != 0x8000))
3080 enan (y, (*p & 0x8000) != 0);
3087 #ifdef ARM_EXTENDED_IEEE_FORMAT
3088 for (i = 2; i <= 5; i++)
3092 enan (y, (*p & 0x8000) != 0);
3097 /* In Motorola extended precision format, the most significant
3098 bit of an infinity mantissa could be either 1 or 0. It is
3099 the lower order bits that tell whether the value is a NaN. */
3100 if ((pe[2] & 0x7fff) != 0)
3103 for (i = 3; i <= 5; i++)
3108 enan (y, (*p & 0x8000) != 0);
3112 #endif /* not ARM */
3121 #endif /* INFINITY */
3124 for (i = 0; i < NE; i++)
3128 /* Convert 128-bit long double precision float PE to e type Y. */
3132 unsigned EMUSHORT *pe, *y;
3134 register unsigned EMUSHORT r;
3135 unsigned EMUSHORT *e, *p;
3136 unsigned EMUSHORT yy[NI];
3143 if (! REAL_WORDS_BIG_ENDIAN)
3155 if (! REAL_WORDS_BIG_ENDIAN)
3157 for (i = 0; i < 7; i++)
3161 enan (y, yy[0] != 0);
3168 for (i = 1; i < 8; i++)
3172 enan (y, yy[0] != 0);
3184 #endif /* INFINITY */
3188 if (! REAL_WORDS_BIG_ENDIAN)
3190 for (i = 0; i < 7; i++)
3196 for (i = 0; i < 7; i++)
3200 /* If denormal, remove the implied bit; else shift down 1. */
3213 /* Convert single precision float PE to e type Y. */
3217 unsigned EMUSHORT *pe, *y;
3221 ibmtoe (pe, y, SFmode);
3224 register unsigned EMUSHORT r;
3225 register unsigned EMUSHORT *e, *p;
3226 unsigned EMUSHORT yy[NI];
3230 denorm = 0; /* flag if denormalized number */
3233 if (! REAL_WORDS_BIG_ENDIAN)
3243 yy[M] = (r & 0x7f) | 0200;
3244 r &= ~0x807f; /* strip sign and 7 significand bits */
3249 if (REAL_WORDS_BIG_ENDIAN)
3251 if (((pe[0] & 0x7f) != 0) || (pe[1] != 0))
3253 enan (y, yy[0] != 0);
3259 if (((pe[1] & 0x7f) != 0) || (pe[0] != 0))
3261 enan (y, yy[0] != 0);
3272 #endif /* INFINITY */
3274 /* If zero exponent, then the significand is denormalized.
3275 So take back the understood high significand bit. */
3288 if (! REAL_WORDS_BIG_ENDIAN)
3298 { /* if zero exponent, then normalize the significand */
3299 if ((k = enormlz (yy)) > NBITS)
3302 yy[E] -= (unsigned EMUSHORT) (k - 1);
3305 #endif /* not IBM */
3308 /* Convert e-type X to IEEE 128-bit long double format E. */
3312 unsigned EMUSHORT *x, *e;
3314 unsigned EMUSHORT xi[NI];
3321 make_nan (e, eisneg (x), TFmode);
3326 exp = (EMULONG) xi[E];
3331 /* round off to nearest or even */
3334 emdnorm (xi, 0, 0, exp, 64);
3340 /* Convert exploded e-type X, that has already been rounded to
3341 113-bit precision, to IEEE 128-bit long double format Y. */
3345 unsigned EMUSHORT *a, *b;
3347 register unsigned EMUSHORT *p, *q;
3348 unsigned EMUSHORT i;
3353 make_nan (b, eiisneg (a), TFmode);
3358 if (REAL_WORDS_BIG_ENDIAN)
3361 q = b + 7; /* point to output exponent */
3363 /* If not denormal, delete the implied bit. */
3368 /* combine sign and exponent */
3370 if (REAL_WORDS_BIG_ENDIAN)
3373 *q++ = *p++ | 0x8000;
3380 *q-- = *p++ | 0x8000;
3384 /* skip over guard word */
3386 /* move the significand */
3387 if (REAL_WORDS_BIG_ENDIAN)
3389 for (i = 0; i < 7; i++)
3394 for (i = 0; i < 7; i++)
3399 /* Convert e-type X to IEEE double extended format E. */
3403 unsigned EMUSHORT *x, *e;
3405 unsigned EMUSHORT xi[NI];
3412 make_nan (e, eisneg (x), XFmode);
3417 /* adjust exponent for offset */
3418 exp = (EMULONG) xi[E];
3423 /* round off to nearest or even */
3426 emdnorm (xi, 0, 0, exp, 64);
3432 /* Convert exploded e-type X, that has already been rounded to
3433 64-bit precision, to IEEE double extended format Y. */
3437 unsigned EMUSHORT *a, *b;
3439 register unsigned EMUSHORT *p, *q;
3440 unsigned EMUSHORT i;
3445 make_nan (b, eiisneg (a), XFmode);
3449 /* Shift denormal long double Intel format significand down one bit. */
3450 if ((a[E] == 0) && ! REAL_WORDS_BIG_ENDIAN)
3460 if (REAL_WORDS_BIG_ENDIAN)
3464 q = b + 4; /* point to output exponent */
3465 #if LONG_DOUBLE_TYPE_SIZE == 96
3466 /* Clear the last two bytes of 12-byte Intel format */
3472 /* combine sign and exponent */
3476 *q++ = *p++ | 0x8000;
3483 *q-- = *p++ | 0x8000;
3488 if (REAL_WORDS_BIG_ENDIAN)
3490 #ifdef ARM_EXTENDED_IEEE_FORMAT
3491 /* The exponent is in the lowest 15 bits of the first word. */
3492 *q++ = i ? 0x8000 : 0;
3496 *q++ = *p++ | 0x8000;
3505 *q-- = *p++ | 0x8000;
3510 /* skip over guard word */
3512 /* move the significand */
3514 for (i = 0; i < 4; i++)
3518 for (i = 0; i < 4; i++)
3522 if (REAL_WORDS_BIG_ENDIAN)
3524 for (i = 0; i < 4; i++)
3532 /* Intel long double infinity significand. */
3540 for (i = 0; i < 4; i++)
3546 /* e type to double precision. */
3549 /* Convert e-type X to DEC-format double E. */
3553 unsigned EMUSHORT *x, *e;
3555 etodec (x, e); /* see etodec.c */
3558 /* Convert exploded e-type X, that has already been rounded to
3559 56-bit double precision, to DEC double Y. */
3563 unsigned EMUSHORT *x, *y;
3570 /* Convert e-type X to IBM 370-format double E. */
3574 unsigned EMUSHORT *x, *e;
3576 etoibm (x, e, DFmode);
3579 /* Convert exploded e-type X, that has already been rounded to
3580 56-bit precision, to IBM 370 double Y. */
3584 unsigned EMUSHORT *x, *y;
3586 toibm (x, y, DFmode);
3589 #else /* it's neither DEC nor IBM */
3591 /* Convert e-type X to IEEE double E. */
3595 unsigned EMUSHORT *x, *e;
3597 unsigned EMUSHORT xi[NI];
3604 make_nan (e, eisneg (x), DFmode);
3609 /* adjust exponent for offsets */
3610 exp = (EMULONG) xi[E] - (EXONE - 0x3ff);
3615 /* round off to nearest or even */
3618 emdnorm (xi, 0, 0, exp, 64);
3624 /* Convert exploded e-type X, that has already been rounded to
3625 53-bit precision, to IEEE double Y. */
3629 unsigned EMUSHORT *x, *y;
3631 unsigned EMUSHORT i;
3632 unsigned EMUSHORT *p;
3637 make_nan (y, eiisneg (x), DFmode);
3643 if (! REAL_WORDS_BIG_ENDIAN)
3646 *y = 0; /* output high order */
3648 *y = 0x8000; /* output sign bit */
3651 if (i >= (unsigned int) 2047)
3652 { /* Saturate at largest number less than infinity. */
3655 if (! REAL_WORDS_BIG_ENDIAN)
3669 *y |= (unsigned EMUSHORT) 0x7fef;
3670 if (! REAL_WORDS_BIG_ENDIAN)
3695 i |= *p++ & (unsigned EMUSHORT) 0x0f; /* *p = xi[M] */
3696 *y |= (unsigned EMUSHORT) i; /* high order output already has sign bit set */
3697 if (! REAL_WORDS_BIG_ENDIAN)
3712 #endif /* not IBM */
3713 #endif /* not DEC */
3717 /* e type to single precision. */
3720 /* Convert e-type X to IBM 370 float E. */
3724 unsigned EMUSHORT *x, *e;
3726 etoibm (x, e, SFmode);
3729 /* Convert exploded e-type X, that has already been rounded to
3730 float precision, to IBM 370 float Y. */
3734 unsigned EMUSHORT *x, *y;
3736 toibm (x, y, SFmode);
3740 /* Convert e-type X to IEEE float E. DEC float is the same as IEEE float. */
3744 unsigned EMUSHORT *x, *e;
3747 unsigned EMUSHORT xi[NI];
3753 make_nan (e, eisneg (x), SFmode);
3758 /* adjust exponent for offsets */
3759 exp = (EMULONG) xi[E] - (EXONE - 0177);
3764 /* round off to nearest or even */
3767 emdnorm (xi, 0, 0, exp, 64);
3773 /* Convert exploded e-type X, that has already been rounded to
3774 float precision, to IEEE float Y. */
3778 unsigned EMUSHORT *x, *y;
3780 unsigned EMUSHORT i;
3781 unsigned EMUSHORT *p;
3786 make_nan (y, eiisneg (x), SFmode);
3792 if (! REAL_WORDS_BIG_ENDIAN)
3798 *y = 0; /* output high order */
3800 *y = 0x8000; /* output sign bit */
3803 /* Handle overflow cases. */
3807 *y |= (unsigned EMUSHORT) 0x7f80;
3812 if (! REAL_WORDS_BIG_ENDIAN)
3820 #else /* no INFINITY */
3821 *y |= (unsigned EMUSHORT) 0x7f7f;
3826 if (! REAL_WORDS_BIG_ENDIAN)
3837 #endif /* no INFINITY */
3849 i |= *p++ & (unsigned EMUSHORT) 0x7f; /* *p = xi[M] */
3850 /* High order output already has sign bit set. */
3856 if (! REAL_WORDS_BIG_ENDIAN)
3865 #endif /* not IBM */
3867 /* Compare two e type numbers.
3871 -2 if either a or b is a NaN. */
3875 unsigned EMUSHORT *a, *b;
3877 unsigned EMUSHORT ai[NI], bi[NI];
3878 register unsigned EMUSHORT *p, *q;
3883 if (eisnan (a) || eisnan (b))
3892 { /* the signs are different */
3894 for (i = 1; i < NI - 1; i++)
3908 /* both are the same sign */
3923 return (0); /* equality */
3927 if (*(--p) > *(--q))
3928 return (msign); /* p is bigger */
3930 return (-msign); /* p is littler */
3933 /* Find e-type nearest integer to X, as floor (X + 0.5). */
3937 unsigned EMUSHORT *x, *y;
3943 /* Convert HOST_WIDE_INT LP to e type Y. */
3948 unsigned EMUSHORT *y;
3950 unsigned EMUSHORT yi[NI];
3951 unsigned HOST_WIDE_INT ll;
3957 /* make it positive */
3958 ll = (unsigned HOST_WIDE_INT) (-(*lp));
3959 yi[0] = 0xffff; /* put correct sign in the e type number */
3963 ll = (unsigned HOST_WIDE_INT) (*lp);
3965 /* move the long integer to yi significand area */
3966 #if HOST_BITS_PER_WIDE_INT == 64
3967 yi[M] = (unsigned EMUSHORT) (ll >> 48);
3968 yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
3969 yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
3970 yi[M + 3] = (unsigned EMUSHORT) ll;
3971 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
3973 yi[M] = (unsigned EMUSHORT) (ll >> 16);
3974 yi[M + 1] = (unsigned EMUSHORT) ll;
3975 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
3978 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
3979 ecleaz (yi); /* it was zero */
3981 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
3982 emovo (yi, y); /* output the answer */
3985 /* Convert unsigned HOST_WIDE_INT LP to e type Y. */
3989 unsigned HOST_WIDE_INT *lp;
3990 unsigned EMUSHORT *y;
3992 unsigned EMUSHORT yi[NI];
3993 unsigned HOST_WIDE_INT ll;
3999 /* move the long integer to ayi significand area */
4000 #if HOST_BITS_PER_WIDE_INT == 64
4001 yi[M] = (unsigned EMUSHORT) (ll >> 48);
4002 yi[M + 1] = (unsigned EMUSHORT) (ll >> 32);
4003 yi[M + 2] = (unsigned EMUSHORT) (ll >> 16);
4004 yi[M + 3] = (unsigned EMUSHORT) ll;
4005 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
4007 yi[M] = (unsigned EMUSHORT) (ll >> 16);
4008 yi[M + 1] = (unsigned EMUSHORT) ll;
4009 yi[E] = EXONE + 15; /* exponent if normalize shift count were 0 */
4012 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
4013 ecleaz (yi); /* it was zero */
4015 yi[E] -= (unsigned EMUSHORT) k; /* subtract shift count from exponent */
4016 emovo (yi, y); /* output the answer */
4020 /* Find signed HOST_WIDE_INT integer I and floating point fractional
4021 part FRAC of e-type (packed internal format) floating point input X.
4022 The integer output I has the sign of the input, except that
4023 positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
4024 The output e-type fraction FRAC is the positive fractional
4029 unsigned EMUSHORT *x;
4031 unsigned EMUSHORT *frac;
4033 unsigned EMUSHORT xi[NI];
4035 unsigned HOST_WIDE_INT ll;
4038 k = (int) xi[E] - (EXONE - 1);
4041 /* if exponent <= 0, integer = 0 and real output is fraction */
4046 if (k > (HOST_BITS_PER_WIDE_INT - 1))
4048 /* long integer overflow: output large integer
4049 and correct fraction */
4051 *i = ((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1);
4054 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
4055 /* In this case, let it overflow and convert as if unsigned. */
4056 euifrac (x, &ll, frac);
4057 *i = (HOST_WIDE_INT) ll;
4060 /* In other cases, return the largest positive integer. */
4061 *i = (((unsigned HOST_WIDE_INT) 1) << (HOST_BITS_PER_WIDE_INT - 1)) - 1;
4066 warning ("overflow on truncation to integer");
4070 /* Shift more than 16 bits: first shift up k-16 mod 16,
4071 then shift up by 16's. */
4072 j = k - ((k >> 4) << 4);
4079 ll = (ll << 16) | xi[M];
4081 while ((k -= 16) > 0);
4088 /* shift not more than 16 bits */
4090 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4097 if ((k = enormlz (xi)) > NBITS)
4100 xi[E] -= (unsigned EMUSHORT) k;
4106 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
4107 FRAC of e-type X. A negative input yields integer output = 0 but
4108 correct fraction. */
4111 euifrac (x, i, frac)
4112 unsigned EMUSHORT *x;
4113 unsigned HOST_WIDE_INT *i;
4114 unsigned EMUSHORT *frac;
4116 unsigned HOST_WIDE_INT ll;
4117 unsigned EMUSHORT xi[NI];
4121 k = (int) xi[E] - (EXONE - 1);
4124 /* if exponent <= 0, integer = 0 and argument is fraction */
4129 if (k > HOST_BITS_PER_WIDE_INT)
4131 /* Long integer overflow: output large integer
4132 and correct fraction.
4133 Note, the BSD microvax compiler says that ~(0UL)
4134 is a syntax error. */
4138 warning ("overflow on truncation to unsigned integer");
4142 /* Shift more than 16 bits: first shift up k-16 mod 16,
4143 then shift up by 16's. */
4144 j = k - ((k >> 4) << 4);
4151 ll = (ll << 16) | xi[M];
4153 while ((k -= 16) > 0);
4158 /* shift not more than 16 bits */
4160 *i = (HOST_WIDE_INT) xi[M] & 0xffff;
4163 if (xi[0]) /* A negative value yields unsigned integer 0. */
4169 if ((k = enormlz (xi)) > NBITS)
4172 xi[E] -= (unsigned EMUSHORT) k;
4177 /* Shift the significand of exploded e-type X up or down by SC bits. */
4181 unsigned EMUSHORT *x;
4184 unsigned EMUSHORT lost;
4185 unsigned EMUSHORT *p;
4198 lost |= *p; /* remember lost bits */
4239 return ((int) lost);
4242 /* Shift normalize the significand area of exploded e-type X.
4243 Return the shift count (up = positive). */
4247 unsigned EMUSHORT x[];
4249 register unsigned EMUSHORT *p;
4258 return (0); /* already normalized */
4264 /* With guard word, there are NBITS+16 bits available.
4265 Return true if all are zero. */
4269 /* see if high byte is zero */
4270 while ((*p & 0xff00) == 0)
4275 /* now shift 1 bit at a time */
4276 while ((*p & 0x8000) == 0)
4282 mtherr ("enormlz", UNDERFLOW);
4288 /* Normalize by shifting down out of the high guard word
4289 of the significand */
4304 mtherr ("enormlz", OVERFLOW);
4311 /* Powers of ten used in decimal <-> binary conversions. */
4316 #if LONG_DOUBLE_TYPE_SIZE == 128
4317 static unsigned EMUSHORT etens[NTEN + 1][NE] =
4319 {0x6576, 0x4a92, 0x804a, 0x153f,
4320 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4321 {0x6a32, 0xce52, 0x329a, 0x28ce,
4322 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4323 {0x526c, 0x50ce, 0xf18b, 0x3d28,
4324 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4325 {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4326 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4327 {0x851e, 0xeab7, 0x98fe, 0x901b,
4328 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4329 {0x0235, 0x0137, 0x36b1, 0x336c,
4330 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4331 {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4332 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4333 {0x0000, 0x0000, 0x0000, 0x0000,
4334 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4335 {0x0000, 0x0000, 0x0000, 0x0000,
4336 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4337 {0x0000, 0x0000, 0x0000, 0x0000,
4338 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4339 {0x0000, 0x0000, 0x0000, 0x0000,
4340 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4341 {0x0000, 0x0000, 0x0000, 0x0000,
4342 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4343 {0x0000, 0x0000, 0x0000, 0x0000,
4344 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4347 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4349 {0x2030, 0xcffc, 0xa1c3, 0x8123,
4350 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4351 {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4352 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4353 {0xf53f, 0xf698, 0x6bd3, 0x0158,
4354 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4355 {0xe731, 0x04d4, 0xe3f2, 0xd332,
4356 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4357 {0xa23e, 0x5308, 0xfefb, 0x1155,
4358 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4359 {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4360 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4361 {0x2a20, 0x6224, 0x47b3, 0x98d7,
4362 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4363 {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4364 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4365 {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4366 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4367 {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4368 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4369 {0xc155, 0xa4a8, 0x404e, 0x6113,
4370 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4371 {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4372 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4373 {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4374 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4377 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4378 static unsigned EMUSHORT etens[NTEN + 1][NE] =
4380 {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4381 {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4382 {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4383 {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4384 {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4385 {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4386 {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4387 {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4388 {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4389 {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4390 {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4391 {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4392 {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4395 static unsigned EMUSHORT emtens[NTEN + 1][NE] =
4397 {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4398 {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4399 {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4400 {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4401 {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4402 {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4403 {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4404 {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4405 {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4406 {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4407 {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4408 {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4409 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4413 /* Convert float value X to ASCII string STRING with NDIG digits after
4414 the decimal point. */
4417 e24toasc (x, string, ndigs)
4418 unsigned EMUSHORT x[];
4422 unsigned EMUSHORT w[NI];
4425 etoasc (w, string, ndigs);
4428 /* Convert double value X to ASCII string STRING with NDIG digits after
4429 the decimal point. */
4432 e53toasc (x, string, ndigs)
4433 unsigned EMUSHORT x[];
4437 unsigned EMUSHORT w[NI];
4440 etoasc (w, string, ndigs);
4443 /* Convert double extended value X to ASCII string STRING with NDIG digits
4444 after the decimal point. */
4447 e64toasc (x, string, ndigs)
4448 unsigned EMUSHORT x[];
4452 unsigned EMUSHORT w[NI];
4455 etoasc (w, string, ndigs);
4458 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4459 after the decimal point. */
4462 e113toasc (x, string, ndigs)
4463 unsigned EMUSHORT x[];
4467 unsigned EMUSHORT w[NI];
4470 etoasc (w, string, ndigs);
4473 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4474 the decimal point. */
4476 static char wstring[80]; /* working storage for ASCII output */
4479 etoasc (x, string, ndigs)
4480 unsigned EMUSHORT x[];
4485 unsigned EMUSHORT y[NI], t[NI], u[NI], w[NI];
4486 unsigned EMUSHORT *p, *r, *ten;
4487 unsigned EMUSHORT sign;
4488 int i, j, k, expon, rndsav;
4490 unsigned EMUSHORT m;
4501 sprintf (wstring, " NaN ");
4505 rndprc = NBITS; /* set to full precision */
4506 emov (x, y); /* retain external format */
4507 if (y[NE - 1] & 0x8000)
4510 y[NE - 1] &= 0x7fff;
4517 ten = &etens[NTEN][0];
4519 /* Test for zero exponent */
4522 for (k = 0; k < NE - 1; k++)
4525 goto tnzro; /* denormalized number */
4527 goto isone; /* valid all zeros */
4531 /* Test for infinity. */
4532 if (y[NE - 1] == 0x7fff)
4535 sprintf (wstring, " -Infinity ");
4537 sprintf (wstring, " Infinity ");
4541 /* Test for exponent nonzero but significand denormalized.
4542 * This is an error condition.
4544 if ((y[NE - 1] != 0) && ((y[NE - 2] & 0x8000) == 0))
4546 mtherr ("etoasc", DOMAIN);
4547 sprintf (wstring, "NaN");
4551 /* Compare to 1.0 */
4560 { /* Number is greater than 1 */
4561 /* Convert significand to an integer and strip trailing decimal zeros. */
4563 u[NE - 1] = EXONE + NBITS - 1;
4565 p = &etens[NTEN - 4][0];
4571 for (j = 0; j < NE - 1; j++)
4584 /* Rescale from integer significand */
4585 u[NE - 1] += y[NE - 1] - (unsigned int) (EXONE + NBITS - 1);
4587 /* Find power of 10 */
4591 /* An unordered compare result shouldn't happen here. */
4592 while (ecmp (ten, u) <= 0)
4594 if (ecmp (p, u) <= 0)
4607 { /* Number is less than 1.0 */
4608 /* Pad significand with trailing decimal zeros. */
4611 while ((y[NE - 2] & 0x8000) == 0)
4620 for (i = 0; i < NDEC + 1; i++)
4622 if ((w[NI - 1] & 0x7) != 0)
4624 /* multiply by 10 */
4637 if (eone[NE - 1] <= u[1])
4649 while (ecmp (eone, w) > 0)
4651 if (ecmp (p, w) >= 0)
4666 /* Find the first (leading) digit. */
4672 digit = equot[NI - 1];
4673 while ((digit == 0) && (ecmp (y, ezero) != 0))
4681 digit = equot[NI - 1];
4689 /* Examine number of digits requested by caller. */
4707 *s++ = (char)digit + '0';
4710 /* Generate digits after the decimal point. */
4711 for (k = 0; k <= ndigs; k++)
4713 /* multiply current number by 10, without normalizing */
4720 *s++ = (char) equot[NI - 1] + '0';
4722 digit = equot[NI - 1];
4725 /* round off the ASCII string */
4728 /* Test for critical rounding case in ASCII output. */
4732 if (ecmp (t, ezero) != 0)
4733 goto roun; /* round to nearest */
4734 if ((*(s - 1) & 1) == 0)
4735 goto doexp; /* round to even */
4737 /* Round up and propagate carry-outs */
4741 /* Carry out to most significant digit? */
4748 /* Most significant digit carries to 10? */
4756 /* Round up and carry out from less significant digits */
4768 sprintf (ss, "e+%d", expon);
4770 sprintf (ss, "e%d", expon);
4772 sprintf (ss, "e%d", expon);
4775 /* copy out the working string */
4778 while (*ss == ' ') /* strip possible leading space */
4780 while ((*s++ = *ss++) != '\0')
4785 /* Convert ASCII string to floating point.
4787 Numeric input is a free format decimal number of any length, with
4788 or without decimal point. Entering E after the number followed by an
4789 integer number causes the second number to be interpreted as a power of
4790 10 to be multiplied by the first number (i.e., "scientific" notation). */
4792 /* Convert ASCII string S to single precision float value Y. */
4797 unsigned EMUSHORT *y;
4803 /* Convert ASCII string S to double precision value Y. */
4808 unsigned EMUSHORT *y;
4810 #if defined(DEC) || defined(IBM)
4818 /* Convert ASCII string S to double extended value Y. */
4823 unsigned EMUSHORT *y;
4828 /* Convert ASCII string S to 128-bit long double Y. */
4833 unsigned EMUSHORT *y;
4835 asctoeg (s, y, 113);
4838 /* Convert ASCII string S to e type Y. */
4843 unsigned EMUSHORT *y;
4845 asctoeg (s, y, NBITS);
4848 /* Convert ASCII string SS to e type Y, with a specified rounding precision
4852 asctoeg (ss, y, oprec)
4854 unsigned EMUSHORT *y;
4857 unsigned EMUSHORT yy[NI], xt[NI], tt[NI];
4858 int esign, decflg, sgnflg, nexp, exp, prec, lost;
4859 int k, trail, c, rndsav;
4861 unsigned EMUSHORT nsign, *p;
4862 char *sp, *s, *lstr;
4864 /* Copy the input string. */
4865 lstr = (char *) alloca (strlen (ss) + 1);
4867 while (*s == ' ') /* skip leading spaces */
4870 while ((*sp++ = *s++) != '\0')
4875 rndprc = NBITS; /* Set to full precision */
4888 if ((k >= 0) && (k <= 9))
4890 /* Ignore leading zeros */
4891 if ((prec == 0) && (decflg == 0) && (k == 0))
4893 /* Identify and strip trailing zeros after the decimal point. */
4894 if ((trail == 0) && (decflg != 0))
4897 while ((*sp >= '0') && (*sp <= '9'))
4899 /* Check for syntax error */
4901 if ((c != 'e') && (c != 'E') && (c != '\0')
4902 && (c != '\n') && (c != '\r') && (c != ' ')
4913 /* If enough digits were given to more than fill up the yy register,
4914 continuing until overflow into the high guard word yy[2]
4915 guarantees that there will be a roundoff bit at the top
4916 of the low guard word after normalization. */
4921 nexp += 1; /* count digits after decimal point */
4922 eshup1 (yy); /* multiply current number by 10 */
4928 xt[NI - 2] = (unsigned EMUSHORT) k;
4933 /* Mark any lost non-zero digit. */
4935 /* Count lost digits before the decimal point. */
4950 case '.': /* decimal point */
4980 mtherr ("asctoe", DOMAIN);
4989 /* Exponent interpretation */
4991 /* 0.0eXXX is zero, regardless of XXX. Check for the 0.0. */
4992 for (k = 0; k < NI; k++)
5003 /* check for + or - */
5011 while ((*s >= '0') && (*s <= '9'))
5015 if (exp > -(MINDECEXP))
5025 if (exp > MAXDECEXP)
5029 yy[E] = 0x7fff; /* infinity */
5032 if (exp < MINDECEXP)
5041 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
5042 while ((nexp > 0) && (yy[2] == 0))
5054 if ((k = enormlz (yy)) > NBITS)
5059 lexp = (EXONE - 1 + NBITS) - k;
5060 emdnorm (yy, lost, 0, lexp, 64);
5062 /* Convert to external format:
5064 Multiply by 10**nexp. If precision is 64 bits,
5065 the maximum relative error incurred in forming 10**n
5066 for 0 <= n <= 324 is 8.2e-20, at 10**180.
5067 For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5068 For 0 >= n >= -999, it is -1.55e-19 at 10**-435. */
5083 /* Punt. Can't handle this without 2 divides. */
5084 emovi (etens[0], tt);
5091 p = &etens[NTEN][0];
5101 while (exp <= MAXP);
5119 /* Round and convert directly to the destination type */
5121 lexp -= EXONE - 0x3ff;
5123 else if (oprec == 24 || oprec == 56)
5124 lexp -= EXONE - (0x41 << 2);
5126 else if (oprec == 24)
5127 lexp -= EXONE - 0177;
5130 else if (oprec == 56)
5131 lexp -= EXONE - 0201;
5134 emdnorm (yy, k, 0, lexp, 64);
5144 todec (yy, y); /* see etodec.c */
5149 toibm (yy, y, DFmode);
5172 /* Return Y = largest integer not greater than X (truncated toward minus
5175 static unsigned EMUSHORT bmask[] =
5198 unsigned EMUSHORT x[], y[];
5200 register unsigned EMUSHORT *p;
5202 unsigned EMUSHORT f[NE];
5204 emov (x, f); /* leave in external format */
5205 expon = (int) f[NE - 1];
5206 e = (expon & 0x7fff) - (EXONE - 1);
5212 /* number of bits to clear out */
5224 /* clear the remaining bits */
5226 /* truncate negatives toward minus infinity */
5229 if ((unsigned EMUSHORT) expon & (unsigned EMUSHORT) 0x8000)
5231 for (i = 0; i < NE - 1; i++)
5243 /* Return S and EXP such that S * 2^EXP = X and .5 <= S < 1.
5244 For example, 1.1 = 0.55 * 2^1. */
5248 unsigned EMUSHORT x[];
5250 unsigned EMUSHORT s[];
5252 unsigned EMUSHORT xi[NI];
5256 /* Handle denormalized numbers properly using long integer exponent. */
5257 li = (EMULONG) ((EMUSHORT) xi[1]);
5265 *exp = (int) (li - 0x3ffe);
5268 /* Return e type Y = X * 2^PWR2. */
5272 unsigned EMUSHORT x[];
5274 unsigned EMUSHORT y[];
5276 unsigned EMUSHORT xi[NI];
5284 emdnorm (xi, i, i, li, 64);
5289 /* C = remainder after dividing B by A, all e type values.
5290 Least significant integer quotient bits left in EQUOT. */
5294 unsigned EMUSHORT a[], b[], c[];
5296 unsigned EMUSHORT den[NI], num[NI];
5300 || (ecmp (a, ezero) == 0)
5308 if (ecmp (a, ezero) == 0)
5310 mtherr ("eremain", SING);
5316 eiremain (den, num);
5317 /* Sign of remainder = sign of quotient */
5325 /* Return quotient of exploded e-types NUM / DEN in EQUOT,
5326 remainder in NUM. */
5330 unsigned EMUSHORT den[], num[];
5333 unsigned EMUSHORT j;
5336 ld -= enormlz (den);
5338 ln -= enormlz (num);
5342 if (ecmpm (den, num) <= 0)
5354 emdnorm (num, 0, 0, ln, 0);
5357 /* Report an error condition CODE encountered in function NAME.
5358 CODE is one of the following:
5360 Mnemonic Value Significance
5362 DOMAIN 1 argument domain error
5363 SING 2 function singularity
5364 OVERFLOW 3 overflow range error
5365 UNDERFLOW 4 underflow range error
5366 TLOSS 5 total loss of precision
5367 PLOSS 6 partial loss of precision
5368 INVALID 7 NaN - producing operation
5369 EDOM 33 Unix domain error code
5370 ERANGE 34 Unix range error code
5372 The order of appearance of the following messages is bound to the
5373 error codes defined above. */
5376 static char *ermsg[NMSGS] =
5378 "unknown", /* error code 0 */
5379 "domain", /* error code 1 */
5380 "singularity", /* et seq. */
5383 "total loss of precision",
5384 "partial loss of precision",
5398 /* The string passed by the calling program is supposed to be the
5399 name of the function in which the error occurred.
5400 The code argument selects which error message string will be printed. */
5402 if ((code <= 0) || (code >= NMSGS))
5404 sprintf (errstr, " %s %s error", name, ermsg[code]);
5407 /* Set global error message word */
5412 /* Convert DEC double precision D to e type E. */
5416 unsigned EMUSHORT *d;
5417 unsigned EMUSHORT *e;
5419 unsigned EMUSHORT y[NI];
5420 register unsigned EMUSHORT r, *p;
5422 ecleaz (y); /* start with a zero */
5423 p = y; /* point to our number */
5424 r = *d; /* get DEC exponent word */
5425 if (*d & (unsigned int) 0x8000)
5426 *p = 0xffff; /* fill in our sign */
5427 ++p; /* bump pointer to our exponent word */
5428 r &= 0x7fff; /* strip the sign bit */
5429 if (r == 0) /* answer = 0 if high order DEC word = 0 */
5433 r >>= 7; /* shift exponent word down 7 bits */
5434 r += EXONE - 0201; /* subtract DEC exponent offset */
5435 /* add our e type exponent offset */
5436 *p++ = r; /* to form our exponent */
5438 r = *d++; /* now do the high order mantissa */
5439 r &= 0177; /* strip off the DEC exponent and sign bits */
5440 r |= 0200; /* the DEC understood high order mantissa bit */
5441 *p++ = r; /* put result in our high guard word */
5443 *p++ = *d++; /* fill in the rest of our mantissa */
5447 eshdn8 (y); /* shift our mantissa down 8 bits */
5452 /* Convert e type X to DEC double precision D. */
5456 unsigned EMUSHORT *x, *d;
5458 unsigned EMUSHORT xi[NI];
5463 /* Adjust exponent for offsets. */
5464 exp = (EMULONG) xi[E] - (EXONE - 0201);
5465 /* Round off to nearest or even. */
5468 emdnorm (xi, 0, 0, exp, 64);
5473 /* Convert exploded e-type X, that has already been rounded to
5474 56-bit precision, to DEC format double Y. */
5478 unsigned EMUSHORT *x, *y;
5480 unsigned EMUSHORT i;
5481 unsigned EMUSHORT *p;
5520 /* Convert IBM single/double precision to e type. */
5524 unsigned EMUSHORT *d;
5525 unsigned EMUSHORT *e;
5526 enum machine_mode mode;
5528 unsigned EMUSHORT y[NI];
5529 register unsigned EMUSHORT r, *p;
5532 ecleaz (y); /* start with a zero */
5533 p = y; /* point to our number */
5534 r = *d; /* get IBM exponent word */
5535 if (*d & (unsigned int) 0x8000)
5536 *p = 0xffff; /* fill in our sign */
5537 ++p; /* bump pointer to our exponent word */
5538 r &= 0x7f00; /* strip the sign bit */
5539 r >>= 6; /* shift exponent word down 6 bits */
5540 /* in fact shift by 8 right and 2 left */
5541 r += EXONE - (0x41 << 2); /* subtract IBM exponent offset */
5542 /* add our e type exponent offset */
5543 *p++ = r; /* to form our exponent */
5545 *p++ = *d++ & 0xff; /* now do the high order mantissa */
5546 /* strip off the IBM exponent and sign bits */
5547 if (mode != SFmode) /* there are only 2 words in SFmode */
5549 *p++ = *d++; /* fill in the rest of our mantissa */
5554 if (y[M] == 0 && y[M+1] == 0 && y[M+2] == 0 && y[M+3] == 0)
5557 y[E] -= 5 + enormlz (y); /* now normalise the mantissa */
5558 /* handle change in RADIX */
5564 /* Convert e type to IBM single/double precision. */
5568 unsigned EMUSHORT *x, *d;
5569 enum machine_mode mode;
5571 unsigned EMUSHORT xi[NI];
5576 exp = (EMULONG) xi[E] - (EXONE - (0x41 << 2)); /* adjust exponent for offsets */
5577 /* round off to nearest or even */
5580 emdnorm (xi, 0, 0, exp, 64);
5582 toibm (xi, d, mode);
5587 unsigned EMUSHORT *x, *y;
5588 enum machine_mode mode;
5590 unsigned EMUSHORT i;
5591 unsigned EMUSHORT *p;
5639 /* Output a binary NaN bit pattern in the target machine's format. */
5641 /* If special NaN bit patterns are required, define them in tm.h
5642 as arrays of unsigned 16-bit shorts. Otherwise, use the default
5648 unsigned EMUSHORT TFbignan[8] =
5649 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5650 unsigned EMUSHORT TFlittlenan[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
5658 unsigned EMUSHORT XFbignan[6] =
5659 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5660 unsigned EMUSHORT XFlittlenan[6] = {0, 0, 0, 0xc000, 0xffff, 0};
5668 unsigned EMUSHORT DFbignan[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
5669 unsigned EMUSHORT DFlittlenan[4] = {0, 0, 0, 0xfff8};
5677 unsigned EMUSHORT SFbignan[2] = {0x7fff, 0xffff};
5678 unsigned EMUSHORT SFlittlenan[2] = {0, 0xffc0};
5684 make_nan (nan, sign, mode)
5685 unsigned EMUSHORT *nan;
5687 enum machine_mode mode;
5690 unsigned EMUSHORT *p;
5694 /* Possibly the `reserved operand' patterns on a VAX can be
5695 used like NaN's, but probably not in the same way as IEEE. */
5696 #if !defined(DEC) && !defined(IBM)
5699 if (REAL_WORDS_BIG_ENDIAN)
5706 if (REAL_WORDS_BIG_ENDIAN)
5713 if (REAL_WORDS_BIG_ENDIAN)
5721 if (REAL_WORDS_BIG_ENDIAN)
5730 if (REAL_WORDS_BIG_ENDIAN)
5731 *nan++ = (sign << 15) | *p++;
5734 if (! REAL_WORDS_BIG_ENDIAN)
5735 *nan = (sign << 15) | *p;
5738 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
5739 This is the inverse of the function `etarsingle' invoked by
5740 REAL_VALUE_TO_TARGET_SINGLE. */
5743 ereal_from_float (f)
5747 unsigned EMUSHORT s[2];
5748 unsigned EMUSHORT e[NE];
5750 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
5751 This is the inverse operation to what the function `endian' does. */
5752 if (REAL_WORDS_BIG_ENDIAN)
5754 s[0] = (unsigned EMUSHORT) (f >> 16);
5755 s[1] = (unsigned EMUSHORT) f;
5759 s[0] = (unsigned EMUSHORT) f;
5760 s[1] = (unsigned EMUSHORT) (f >> 16);
5762 /* Convert and promote the target float to E-type. */
5764 /* Output E-type to REAL_VALUE_TYPE. */
5770 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
5771 This is the inverse of the function `etardouble' invoked by
5772 REAL_VALUE_TO_TARGET_DOUBLE.
5774 The DFmode is stored as an array of HOST_WIDE_INT in the target's
5775 data format, with no holes in the bit packing. The first element
5776 of the input array holds the bits that would come first in the
5777 target computer's memory. */
5780 ereal_from_double (d)
5784 unsigned EMUSHORT s[4];
5785 unsigned EMUSHORT e[NE];
5787 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
5788 if (REAL_WORDS_BIG_ENDIAN)
5790 s[0] = (unsigned EMUSHORT) (d[0] >> 16);
5791 s[1] = (unsigned EMUSHORT) d[0];
5792 #if HOST_BITS_PER_WIDE_INT == 32
5793 s[2] = (unsigned EMUSHORT) (d[1] >> 16);
5794 s[3] = (unsigned EMUSHORT) d[1];
5796 /* In this case the entire target double is contained in the
5797 first array element. The second element of the input is
5799 s[2] = (unsigned EMUSHORT) (d[0] >> 48);
5800 s[3] = (unsigned EMUSHORT) (d[0] >> 32);
5805 /* Target float words are little-endian. */
5806 s[0] = (unsigned EMUSHORT) d[0];
5807 s[1] = (unsigned EMUSHORT) (d[0] >> 16);
5808 #if HOST_BITS_PER_WIDE_INT == 32
5809 s[2] = (unsigned EMUSHORT) d[1];
5810 s[3] = (unsigned EMUSHORT) (d[1] >> 16);
5812 s[2] = (unsigned EMUSHORT) (d[0] >> 32);
5813 s[3] = (unsigned EMUSHORT) (d[0] >> 48);
5816 /* Convert target double to E-type. */
5818 /* Output E-type to REAL_VALUE_TYPE. */
5824 /* Convert target computer unsigned 64-bit integer to e-type.
5825 The endian-ness of DImode follows the convention for integers,
5826 so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN. */
5830 unsigned EMUSHORT *di; /* Address of the 64-bit int. */
5831 unsigned EMUSHORT *e;
5833 unsigned EMUSHORT yi[NI];
5837 if (WORDS_BIG_ENDIAN)
5839 for (k = M; k < M + 4; k++)
5844 for (k = M + 3; k >= M; k--)
5847 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
5848 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
5849 ecleaz (yi); /* it was zero */
5851 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
5855 /* Convert target computer signed 64-bit integer to e-type. */
5859 unsigned EMUSHORT *di; /* Address of the 64-bit int. */
5860 unsigned EMUSHORT *e;
5862 unsigned EMULONG acc;
5863 unsigned EMUSHORT yi[NI];
5864 unsigned EMUSHORT carry;
5868 if (WORDS_BIG_ENDIAN)
5870 for (k = M; k < M + 4; k++)
5875 for (k = M + 3; k >= M; k--)
5878 /* Take absolute value */
5884 for (k = M + 3; k >= M; k--)
5886 acc = (unsigned EMULONG) (~yi[k] & 0xffff) + carry;
5893 yi[E] = EXONE + 47; /* exponent if normalize shift count were 0 */
5894 if ((k = enormlz (yi)) > NBITS)/* normalize the significand */
5895 ecleaz (yi); /* it was zero */
5897 yi[E] -= (unsigned EMUSHORT) k;/* subtract shift count from exponent */
5904 /* Convert e-type to unsigned 64-bit int. */
5908 unsigned EMUSHORT *x;
5909 unsigned EMUSHORT *i;
5911 unsigned EMUSHORT xi[NI];
5920 k = (int) xi[E] - (EXONE - 1);
5923 for (j = 0; j < 4; j++)
5929 for (j = 0; j < 4; j++)
5932 warning ("overflow on truncation to integer");
5937 /* Shift more than 16 bits: first shift up k-16 mod 16,
5938 then shift up by 16's. */
5939 j = k - ((k >> 4) << 4);
5943 if (WORDS_BIG_ENDIAN)
5954 if (WORDS_BIG_ENDIAN)
5959 while ((k -= 16) > 0);
5963 /* shift not more than 16 bits */
5968 if (WORDS_BIG_ENDIAN)
5987 /* Convert e-type to signed 64-bit int. */
5991 unsigned EMUSHORT *x;
5992 unsigned EMUSHORT *i;
5994 unsigned EMULONG acc;
5995 unsigned EMUSHORT xi[NI];
5996 unsigned EMUSHORT carry;
5997 unsigned EMUSHORT *isave;
6001 k = (int) xi[E] - (EXONE - 1);
6004 for (j = 0; j < 4; j++)
6010 for (j = 0; j < 4; j++)
6013 warning ("overflow on truncation to integer");
6019 /* Shift more than 16 bits: first shift up k-16 mod 16,
6020 then shift up by 16's. */
6021 j = k - ((k >> 4) << 4);
6025 if (WORDS_BIG_ENDIAN)
6036 if (WORDS_BIG_ENDIAN)
6041 while ((k -= 16) > 0);
6045 /* shift not more than 16 bits */
6048 if (WORDS_BIG_ENDIAN)
6064 /* Negate if negative */
6068 if (WORDS_BIG_ENDIAN)
6070 for (k = 0; k < 4; k++)
6072 acc = (unsigned EMULONG) (~(*isave) & 0xffff) + carry;
6073 if (WORDS_BIG_ENDIAN)
6085 /* Longhand square root routine. */
6088 static int esqinited = 0;
6089 static unsigned short sqrndbit[NI];
6093 unsigned EMUSHORT *x, *y;
6095 unsigned EMUSHORT temp[NI], num[NI], sq[NI], xx[NI];
6097 int i, j, k, n, nlups;
6102 sqrndbit[NI - 2] = 1;
6105 /* Check for arg <= 0 */
6106 i = ecmp (x, ezero);
6111 mtherr ("esqrt", DOMAIN);
6127 /* Bring in the arg and renormalize if it is denormal. */
6129 m = (EMULONG) xx[1]; /* local long word exponent */
6133 /* Divide exponent by 2 */
6135 exp = (unsigned short) ((m / 2) + 0x3ffe);
6137 /* Adjust if exponent odd */
6147 n = 8; /* get 8 bits of result per inner loop */
6153 /* bring in next word of arg */
6155 num[NI - 1] = xx[j + 3];
6156 /* Do additional bit on last outer loop, for roundoff. */
6159 for (i = 0; i < n; i++)
6161 /* Next 2 bits of arg */
6164 /* Shift up answer */
6166 /* Make trial divisor */
6167 for (k = 0; k < NI; k++)
6170 eaddm (sqrndbit, temp);
6171 /* Subtract and insert answer bit if it goes in */
6172 if (ecmpm (temp, num) <= 0)
6182 /* Adjust for extra, roundoff loop done. */
6183 exp += (NBITS - 1) - rndprc;
6185 /* Sticky bit = 1 if the remainder is nonzero. */
6187 for (i = 3; i < NI; i++)
6190 /* Renormalize and round off. */
6191 emdnorm (sq, k, 0, exp, 64);
6194 #endif /* EMU_NON_COMPILE not defined */
6196 /* Return the binary precision of the significand for a given
6197 floating point mode. The mode can hold an integer value
6198 that many bits wide, without losing any bits. */
6201 significand_size (mode)
6202 enum machine_mode mode;
6205 /* Don't test the modes, but their sizes, lest this
6206 code won't work for BITS_PER_UNIT != 8 . */
6208 switch (GET_MODE_BITSIZE (mode))
6214 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
6217 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
6220 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT