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 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 To support cross compilation between IEEE, VAX and IBM floating
34 point formats, define REAL_ARITHMETIC in the tm.h file.
36 In either case the machine files (tm.h) must not contain any code
37 that tries to use host floating point arithmetic to convert
38 REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
39 etc. In cross-compile situations a REAL_VALUE_TYPE may not
40 be intelligible to the host computer's native arithmetic.
42 The emulator defaults to the host's floating point format so that
43 its decimal conversion functions can be used if desired (see
46 The first part of this file interfaces gcc to a floating point
47 arithmetic suite that was not written with gcc in mind. Avoid
48 changing the low-level arithmetic routines unless you have suitable
49 test programs available. A special version of the PARANOIA floating
50 point arithmetic tester, modified for this purpose, can be found on
51 usc.edu: /pub/C-numanal/ieeetest.zoo. Other tests, and libraries of
52 XFmode and TFmode transcendental functions, can be obtained by ftp from
53 netlib.att.com: netlib/cephes. */
55 /* Type of computer arithmetic.
56 Only one of DEC, IBM, IEEE, C4X, or UNK should get defined.
58 `IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
59 to big-endian IEEE floating-point data structure. This definition
60 should work in SFmode `float' type and DFmode `double' type on
61 virtually all big-endian IEEE machines. If LONG_DOUBLE_TYPE_SIZE
62 has been defined to be 96, then IEEE also invokes the particular
63 XFmode (`long double' type) data structure used by the Motorola
64 680x0 series processors.
66 `IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
67 little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
68 has been defined to be 96, then IEEE also invokes the particular
69 XFmode `long double' data structure used by the Intel 80x86 series
72 `DEC' refers specifically to the Digital Equipment Corp PDP-11
73 and VAX floating point data structure. This model currently
74 supports no type wider than DFmode.
76 `IBM' refers specifically to the IBM System/370 and compatible
77 floating point data structure. This model currently supports
78 no type wider than DFmode. The IBM conversions were contributed by
79 frank@atom.ansto.gov.au (Frank Crawford).
81 `C4X' refers specifically to the floating point format used on
82 Texas Instruments TMS320C3x and TMS320C4x digital signal
83 processors. This supports QFmode (32-bit float, double) and HFmode
84 (40-bit long double) where BITS_PER_BYTE is 32. Unlike IEEE
85 floats, C4x floats are not rounded to be even. The C4x conversions
86 were contributed by m.hayes@elec.canterbury.ac.nz (Michael Hayes) and
87 Haj.Ten.Brugge@net.HCC.nl (Herman ten Brugge).
89 If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
90 then `long double' and `double' are both implemented, but they
91 both mean DFmode. In this case, the software floating-point
92 support available here is activated by writing
93 #define REAL_ARITHMETIC
96 The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
97 and may deactivate XFmode since `long double' is used to refer
98 to both modes. Defining INTEL_EXTENDED_IEEE_FORMAT to non-zero
99 at the same time enables 80387-style 80-bit floats in a 128-bit
100 padded image, as seen on IA-64.
102 The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
103 contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
104 separate the floating point unit's endian-ness from that of
105 the integer addressing. This permits one to define a big-endian
106 FPU on a little-endian machine (e.g., ARM). An extension to
107 BYTES_BIG_ENDIAN may be required for some machines in the future.
108 These optional macros may be defined in tm.h. In real.h, they
109 default to WORDS_BIG_ENDIAN, etc., so there is no need to define
110 them for any normal host or target machine on which the floats
111 and the integers have the same endian-ness. */
114 /* The following converts gcc macros into the ones used by this file. */
116 /* REAL_ARITHMETIC defined means that macros in real.h are
117 defined to call emulator functions. */
118 #ifdef REAL_ARITHMETIC
120 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
121 /* PDP-11, Pro350, VAX: */
123 #else /* it's not VAX */
124 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
125 /* IBM System/370 style */
127 #else /* it's also not an IBM */
128 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
129 /* TMS320C3x/C4x style */
131 #else /* it's also not a C4X */
132 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
134 #else /* it's not IEEE either */
135 /* UNKnown arithmetic. We don't support this and can't go on. */
136 unknown arithmetic type
138 #endif /* not IEEE */
143 #define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
146 /* REAL_ARITHMETIC not defined means that the *host's* data
147 structure will be used. It may differ by endian-ness from the
148 target machine's structure and will get its ends swapped
149 accordingly (but not here). Probably only the decimal <-> binary
150 functions in this file will actually be used in this case. */
152 #if HOST_FLOAT_FORMAT == VAX_FLOAT_FORMAT
154 #else /* it's not VAX */
155 #if HOST_FLOAT_FORMAT == IBM_FLOAT_FORMAT
156 /* IBM System/370 style */
158 #else /* it's also not an IBM */
159 #if HOST_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
161 #else /* it's not IEEE either */
162 unknown arithmetic type
164 #endif /* not IEEE */
168 #define REAL_WORDS_BIG_ENDIAN HOST_FLOAT_WORDS_BIG_ENDIAN
170 #endif /* REAL_ARITHMETIC not defined */
172 /* Define INFINITY for support of infinity.
173 Define NANS for support of Not-a-Number's (NaN's). */
174 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
179 /* Support of NaNs requires support of infinity. */
186 /* Find a host integer type that is at least 16 bits wide,
187 and another type at least twice whatever that size is. */
189 #if HOST_BITS_PER_CHAR >= 16
190 #define EMUSHORT char
191 #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
192 #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
194 #if HOST_BITS_PER_SHORT >= 16
195 #define EMUSHORT short
196 #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
197 #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
199 #if HOST_BITS_PER_INT >= 16
201 #define EMUSHORT_SIZE HOST_BITS_PER_INT
202 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
204 #if HOST_BITS_PER_LONG >= 16
205 #define EMUSHORT long
206 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
207 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
209 /* You will have to modify this program to have a smaller unit size. */
210 #define EMU_NON_COMPILE
216 /* If no 16-bit type has been found and the compiler is GCC, try HImode. */
217 #if defined(__GNUC__) && EMUSHORT_SIZE != 16
218 typedef int HItype
__attribute__ ((mode (HI
)));
219 typedef unsigned int UHItype
__attribute__ ((mode (HI
)));
223 #define EMUSHORT HItype
224 #define UEMUSHORT UHItype
225 #define EMUSHORT_SIZE 16
226 #define EMULONG_SIZE 32
228 #define UEMUSHORT unsigned EMUSHORT
231 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
232 #define EMULONG short
234 #if HOST_BITS_PER_INT >= EMULONG_SIZE
237 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
240 #if HOST_BITS_PER_LONGLONG >= EMULONG_SIZE
241 #define EMULONG long long int
243 /* You will have to modify this program to have a smaller unit size. */
244 #define EMU_NON_COMPILE
251 /* The host interface doesn't work if no 16-bit size exists. */
252 #if EMUSHORT_SIZE != 16
253 #define EMU_NON_COMPILE
256 /* OK to continue compilation. */
257 #ifndef EMU_NON_COMPILE
259 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
260 In GET_REAL and PUT_REAL, r and e are pointers.
261 A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
262 in memory, with no holes. */
264 #if MAX_LONG_DOUBLE_TYPE_SIZE == 96 || \
265 ((INTEL_EXTENDED_IEEE_FORMAT != 0) && MAX_LONG_DOUBLE_TYPE_SIZE == 128)
266 /* Number of 16 bit words in external e type format */
268 # define MAXDECEXP 4932
269 # define MINDECEXP -4956
270 # define GET_REAL(r,e) memcpy ((char *)(e), (char *)(r), 2*NE)
271 # define PUT_REAL(e,r) \
273 memcpy ((char *)(r), (char *)(e), 2*NE); \
274 if (2*NE < sizeof(*r)) \
275 memset ((char *)(r) + 2*NE, 0, sizeof(*r) - 2*NE); \
277 # else /* no XFmode */
278 # if MAX_LONG_DOUBLE_TYPE_SIZE == 128
280 # define MAXDECEXP 4932
281 # define MINDECEXP -4977
282 # define GET_REAL(r,e) memcpy ((char *)(e), (char *)(r), 2*NE)
283 # define PUT_REAL(e,r) \
285 memcpy ((char *)(r), (char *)(e), 2*NE); \
286 if (2*NE < sizeof(*r)) \
287 memset ((char *)(r) + 2*NE, 0, sizeof(*r) - 2*NE); \
291 #define MAXDECEXP 4932
292 #define MINDECEXP -4956
293 #ifdef REAL_ARITHMETIC
294 /* Emulator uses target format internally
295 but host stores it in host endian-ness. */
297 #define GET_REAL(r,e) \
299 if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
300 e53toe ((UEMUSHORT *) (r), (e)); \
304 memcpy (&w[3], ((EMUSHORT *) r), sizeof (EMUSHORT)); \
305 memcpy (&w[2], ((EMUSHORT *) r) + 1, sizeof (EMUSHORT)); \
306 memcpy (&w[1], ((EMUSHORT *) r) + 2, sizeof (EMUSHORT)); \
307 memcpy (&w[0], ((EMUSHORT *) r) + 3, sizeof (EMUSHORT)); \
312 #define PUT_REAL(e,r) \
314 if (HOST_FLOAT_WORDS_BIG_ENDIAN == REAL_WORDS_BIG_ENDIAN) \
315 etoe53 ((e), (UEMUSHORT *) (r)); \
320 memcpy (((EMUSHORT *) r), &w[3], sizeof (EMUSHORT)); \
321 memcpy (((EMUSHORT *) r) + 1, &w[2], sizeof (EMUSHORT)); \
322 memcpy (((EMUSHORT *) r) + 2, &w[1], sizeof (EMUSHORT)); \
323 memcpy (((EMUSHORT *) r) + 3, &w[0], sizeof (EMUSHORT)); \
327 #else /* not REAL_ARITHMETIC */
329 /* emulator uses host format */
330 #define GET_REAL(r,e) e53toe ((UEMUSHORT *) (r), (e))
331 #define PUT_REAL(e,r) etoe53 ((e), (UEMUSHORT *) (r))
333 #endif /* not REAL_ARITHMETIC */
334 #endif /* not TFmode */
335 #endif /* not XFmode */
338 /* Number of 16 bit words in internal format */
341 /* Array offset to exponent */
344 /* Array offset to high guard word */
347 /* Number of bits of precision */
348 #define NBITS ((NI-4)*16)
350 /* Maximum number of decimal digits in ASCII conversion
353 #define NDEC (NBITS*8/27)
355 /* The exponent of 1.0 */
356 #define EXONE (0x3fff)
358 #if defined(HOST_EBCDIC)
359 /* bit 8 is significant in EBCDIC */
360 #define CHARMASK 0xff
362 #define CHARMASK 0x7f
365 extern int extra_warnings
;
366 extern UEMUSHORT ezero
[], ehalf
[], eone
[], etwo
[];
367 extern UEMUSHORT elog2
[], esqrt2
[];
369 static void endian
PARAMS ((UEMUSHORT
*, long *,
371 static void eclear
PARAMS ((UEMUSHORT
*));
372 static void emov
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
374 static void eabs
PARAMS ((UEMUSHORT
*));
376 static void eneg
PARAMS ((UEMUSHORT
*));
377 static int eisneg
PARAMS ((UEMUSHORT
*));
378 static int eisinf
PARAMS ((UEMUSHORT
*));
379 static int eisnan
PARAMS ((UEMUSHORT
*));
380 static void einfin
PARAMS ((UEMUSHORT
*));
382 static void enan
PARAMS ((UEMUSHORT
*, int));
383 static void einan
PARAMS ((UEMUSHORT
*));
384 static int eiisnan
PARAMS ((UEMUSHORT
*));
385 static int eiisneg
PARAMS ((UEMUSHORT
*));
386 static void make_nan
PARAMS ((UEMUSHORT
*, int, enum machine_mode
));
388 static void emovi
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
389 static void emovo
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
390 static void ecleaz
PARAMS ((UEMUSHORT
*));
391 static void ecleazs
PARAMS ((UEMUSHORT
*));
392 static void emovz
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
394 static void eiinfin
PARAMS ((UEMUSHORT
*));
397 static int eiisinf
PARAMS ((UEMUSHORT
*));
399 static int ecmpm
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
400 static void eshdn1
PARAMS ((UEMUSHORT
*));
401 static void eshup1
PARAMS ((UEMUSHORT
*));
402 static void eshdn8
PARAMS ((UEMUSHORT
*));
403 static void eshup8
PARAMS ((UEMUSHORT
*));
404 static void eshup6
PARAMS ((UEMUSHORT
*));
405 static void eshdn6
PARAMS ((UEMUSHORT
*));
406 static void eaddm
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));\f
407 static void esubm
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
408 static void m16m
PARAMS ((unsigned int, UEMUSHORT
*, UEMUSHORT
*));
409 static int edivm
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
410 static int emulm
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
411 static void emdnorm
PARAMS ((UEMUSHORT
*, int, int, EMULONG
, int));
412 static void esub
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
414 static void eadd
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
416 static void eadd1
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
418 static void ediv
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
420 static void emul
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
422 static void e53toe
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
423 static void e64toe
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
424 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
425 static void e113toe
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
427 static void e24toe
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
428 static void etoe113
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
429 static void toe113
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
430 static void etoe64
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
431 static void toe64
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
432 static void etoe53
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
433 static void toe53
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
434 static void etoe24
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
435 static void toe24
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
436 static int ecmp
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
438 static void eround
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
440 static void ltoe
PARAMS ((HOST_WIDE_INT
*, UEMUSHORT
*));
441 static void ultoe
PARAMS ((unsigned HOST_WIDE_INT
*, UEMUSHORT
*));
442 static void eifrac
PARAMS ((UEMUSHORT
*, HOST_WIDE_INT
*,
444 static void euifrac
PARAMS ((UEMUSHORT
*, unsigned HOST_WIDE_INT
*,
446 static int eshift
PARAMS ((UEMUSHORT
*, int));
447 static int enormlz
PARAMS ((UEMUSHORT
*));
449 static void e24toasc
PARAMS ((UEMUSHORT
*, char *, int));
450 static void e53toasc
PARAMS ((UEMUSHORT
*, char *, int));
451 static void e64toasc
PARAMS ((UEMUSHORT
*, char *, int));
452 static void e113toasc
PARAMS ((UEMUSHORT
*, char *, int));
454 static void etoasc
PARAMS ((UEMUSHORT
*, char *, int));
455 static void asctoe24
PARAMS ((const char *, UEMUSHORT
*));
456 static void asctoe53
PARAMS ((const char *, UEMUSHORT
*));
457 static void asctoe64
PARAMS ((const char *, UEMUSHORT
*));
458 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
459 static void asctoe113
PARAMS ((const char *, UEMUSHORT
*));
461 static void asctoe
PARAMS ((const char *, UEMUSHORT
*));
462 static void asctoeg
PARAMS ((const char *, UEMUSHORT
*, int));
463 static void efloor
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
465 static void efrexp
PARAMS ((UEMUSHORT
*, int *,
468 static void eldexp
PARAMS ((UEMUSHORT
*, int, UEMUSHORT
*));
470 static void eremain
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
473 static void eiremain
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
474 static void mtherr
PARAMS ((const char *, int));
476 static void dectoe
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
477 static void etodec
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
478 static void todec
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
481 static void ibmtoe
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
483 static void etoibm
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
485 static void toibm
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
489 static void c4xtoe
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
491 static void etoc4x
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
493 static void toc4x
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
497 static void uditoe
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
498 static void ditoe
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
499 static void etoudi
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
500 static void etodi
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
501 static void esqrt
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
504 /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
505 swapping ends if required, into output array of longs. The
506 result is normally passed to fprintf by the ASM_OUTPUT_ macros. */
512 enum machine_mode mode
;
516 if (REAL_WORDS_BIG_ENDIAN
)
521 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
522 /* Swap halfwords in the fourth long. */
523 th
= (unsigned long) e
[6] & 0xffff;
524 t
= (unsigned long) e
[7] & 0xffff;
530 /* Swap halfwords in the third long. */
531 th
= (unsigned long) e
[4] & 0xffff;
532 t
= (unsigned long) e
[5] & 0xffff;
535 /* fall into the double case */
538 /* Swap halfwords in the second word. */
539 th
= (unsigned long) e
[2] & 0xffff;
540 t
= (unsigned long) e
[3] & 0xffff;
543 /* fall into the float case */
547 /* Swap halfwords in the first word. */
548 th
= (unsigned long) e
[0] & 0xffff;
549 t
= (unsigned long) e
[1] & 0xffff;
560 /* Pack the output array without swapping. */
565 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
566 /* Pack the fourth long. */
567 th
= (unsigned long) e
[7] & 0xffff;
568 t
= (unsigned long) e
[6] & 0xffff;
574 /* Pack the third long.
575 Each element of the input REAL_VALUE_TYPE array has 16 useful bits
577 th
= (unsigned long) e
[5] & 0xffff;
578 t
= (unsigned long) e
[4] & 0xffff;
581 /* fall into the double case */
584 /* Pack the second long */
585 th
= (unsigned long) e
[3] & 0xffff;
586 t
= (unsigned long) e
[2] & 0xffff;
589 /* fall into the float case */
593 /* Pack the first long */
594 th
= (unsigned long) e
[1] & 0xffff;
595 t
= (unsigned long) e
[0] & 0xffff;
607 /* This is the implementation of the REAL_ARITHMETIC macro. */
610 earith (value
, icode
, r1
, r2
)
611 REAL_VALUE_TYPE
*value
;
616 UEMUSHORT d1
[NE
], d2
[NE
], v
[NE
];
622 /* Return NaN input back to the caller. */
625 PUT_REAL (d1
, value
);
630 PUT_REAL (d2
, value
);
634 code
= (enum tree_code
) icode
;
642 esub (d2
, d1
, v
); /* d1 - d2 */
650 #ifndef REAL_INFINITY
651 if (ecmp (d2
, ezero
) == 0)
654 enan (v
, eisneg (d1
) ^ eisneg (d2
));
661 ediv (d2
, d1
, v
); /* d1/d2 */
664 case MIN_EXPR
: /* min (d1,d2) */
665 if (ecmp (d1
, d2
) < 0)
671 case MAX_EXPR
: /* max (d1,d2) */
672 if (ecmp (d1
, d2
) > 0)
685 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
686 implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */
692 UEMUSHORT f
[NE
], g
[NE
];
708 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
709 implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */
715 UEMUSHORT f
[NE
], g
[NE
];
717 unsigned HOST_WIDE_INT l
;
731 /* This is the REAL_VALUE_ATOF function. It converts a decimal or hexadecimal
732 string to binary, rounding off as indicated by the machine_mode argument.
733 Then it promotes the rounded value to REAL_VALUE_TYPE. */
740 UEMUSHORT tem
[NE
], e
[NE
];
766 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
786 /* Expansion of REAL_NEGATE. */
802 /* Round real toward zero to HOST_WIDE_INT;
803 implements REAL_VALUE_FIX (x). */
809 UEMUSHORT f
[NE
], g
[NE
];
816 warning ("conversion from NaN to int");
824 /* Round real toward zero to unsigned HOST_WIDE_INT
825 implements REAL_VALUE_UNSIGNED_FIX (x).
826 Negative input returns zero. */
828 unsigned HOST_WIDE_INT
832 UEMUSHORT f
[NE
], g
[NE
];
833 unsigned HOST_WIDE_INT l
;
839 warning ("conversion from NaN to unsigned int");
848 /* REAL_VALUE_FROM_INT macro. */
851 ereal_from_int (d
, i
, j
, mode
)
854 enum machine_mode mode
;
856 UEMUSHORT df
[NE
], dg
[NE
];
857 HOST_WIDE_INT low
, high
;
860 if (GET_MODE_CLASS (mode
) != MODE_FLOAT
)
867 /* complement and add 1 */
874 eldexp (eone
, HOST_BITS_PER_WIDE_INT
, df
);
875 ultoe ((unsigned HOST_WIDE_INT
*) &high
, dg
);
877 ultoe ((unsigned HOST_WIDE_INT
*) &low
, df
);
882 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
883 Avoid double-rounding errors later by rounding off now from the
884 extra-wide internal format to the requested precision. */
885 switch (GET_MODE_BITSIZE (mode
))
903 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
920 /* REAL_VALUE_FROM_UNSIGNED_INT macro. */
923 ereal_from_uint (d
, i
, j
, mode
)
925 unsigned HOST_WIDE_INT i
, j
;
926 enum machine_mode mode
;
928 UEMUSHORT df
[NE
], dg
[NE
];
929 unsigned HOST_WIDE_INT low
, high
;
931 if (GET_MODE_CLASS (mode
) != MODE_FLOAT
)
935 eldexp (eone
, HOST_BITS_PER_WIDE_INT
, df
);
941 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
942 Avoid double-rounding errors later by rounding off now from the
943 extra-wide internal format to the requested precision. */
944 switch (GET_MODE_BITSIZE (mode
))
962 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
979 /* REAL_VALUE_TO_INT macro. */
982 ereal_to_int (low
, high
, rr
)
983 HOST_WIDE_INT
*low
, *high
;
986 UEMUSHORT d
[NE
], df
[NE
], dg
[NE
], dh
[NE
];
993 warning ("conversion from NaN to int");
999 /* convert positive value */
1006 eldexp (eone
, HOST_BITS_PER_WIDE_INT
, df
);
1007 ediv (df
, d
, dg
); /* dg = d / 2^32 is the high word */
1008 euifrac (dg
, (unsigned HOST_WIDE_INT
*) high
, dh
);
1009 emul (df
, dh
, dg
); /* fractional part is the low word */
1010 euifrac (dg
, (unsigned HOST_WIDE_INT
*)low
, dh
);
1013 /* complement and add 1 */
1023 /* REAL_VALUE_LDEXP macro. */
1030 UEMUSHORT e
[NE
], y
[NE
];
1043 /* These routines are conditionally compiled because functions
1044 of the same names may be defined in fold-const.c. */
1046 #ifdef REAL_ARITHMETIC
1048 /* Check for infinity in a REAL_VALUE_TYPE. */
1052 REAL_VALUE_TYPE x ATTRIBUTE_UNUSED
;
1058 return (eisinf (e
));
1064 /* Check whether a REAL_VALUE_TYPE item is a NaN. */
1068 REAL_VALUE_TYPE x ATTRIBUTE_UNUSED
;
1074 return (eisnan (e
));
1081 /* Check for a negative REAL_VALUE_TYPE number.
1082 This just checks the sign bit, so that -0 counts as negative. */
1088 return ereal_isneg (x
);
1091 /* Expansion of REAL_VALUE_TRUNCATE.
1092 The result is in floating point, rounded to nearest or even. */
1095 real_value_truncate (mode
, arg
)
1096 enum machine_mode mode
;
1097 REAL_VALUE_TYPE arg
;
1099 UEMUSHORT e
[NE
], t
[NE
];
1111 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
1148 /* If an unsupported type was requested, presume that
1149 the machine files know something useful to do with
1150 the unmodified value. */
1159 /* Try to change R into its exact multiplicative inverse in machine mode
1160 MODE. Return nonzero function value if successful. */
1163 exact_real_inverse (mode
, r
)
1164 enum machine_mode mode
;
1167 UEMUSHORT e
[NE
], einv
[NE
];
1168 REAL_VALUE_TYPE rinv
;
1173 /* Test for input in range. Don't transform IEEE special values. */
1174 if (eisinf (e
) || eisnan (e
) || (ecmp (e
, ezero
) == 0))
1177 /* Test for a power of 2: all significand bits zero except the MSB.
1178 We are assuming the target has binary (or hex) arithmetic. */
1179 if (e
[NE
- 2] != 0x8000)
1182 for (i
= 0; i
< NE
- 2; i
++)
1188 /* Compute the inverse and truncate it to the required mode. */
1189 ediv (e
, eone
, einv
);
1190 PUT_REAL (einv
, &rinv
);
1191 rinv
= real_value_truncate (mode
, rinv
);
1193 #ifdef CHECK_FLOAT_VALUE
1194 /* This check is not redundant. It may, for example, flush
1195 a supposedly IEEE denormal value to zero. */
1197 if (CHECK_FLOAT_VALUE (mode
, rinv
, i
))
1200 GET_REAL (&rinv
, einv
);
1202 /* Check the bits again, because the truncation might have
1203 generated an arbitrary saturation value on overflow. */
1204 if (einv
[NE
- 2] != 0x8000)
1207 for (i
= 0; i
< NE
- 2; i
++)
1213 /* Fail if the computed inverse is out of range. */
1214 if (eisinf (einv
) || eisnan (einv
) || (ecmp (einv
, ezero
) == 0))
1217 /* Output the reciprocal and return success flag. */
1221 #endif /* REAL_ARITHMETIC defined */
1223 /* Used for debugging--print the value of R in human-readable format
1232 REAL_VALUE_TO_DECIMAL (r
, "%.20g", dstr
);
1233 fprintf (stderr
, "%s", dstr
);
1237 /* The following routines convert REAL_VALUE_TYPE to the various floating
1238 point formats that are meaningful to supported computers.
1240 The results are returned in 32-bit pieces, each piece stored in a `long'.
1241 This is so they can be printed by statements like
1243 fprintf (file, "%lx, %lx", L[0], L[1]);
1245 that will work on both narrow- and wide-word host computers. */
1247 /* Convert R to a 128-bit long double precision value. The output array L
1248 contains four 32-bit pieces of the result, in the order they would appear
1260 endian (e
, l
, TFmode
);
1263 /* Convert R to a double extended precision value. The output array L
1264 contains three 32-bit pieces of the result, in the order they would
1265 appear in memory. */
1276 endian (e
, l
, XFmode
);
1279 /* Convert R to a double precision value. The output array L contains two
1280 32-bit pieces of the result, in the order they would appear in memory. */
1291 endian (e
, l
, DFmode
);
1294 /* Convert R to a single precision float value stored in the least-significant
1295 bits of a `long'. */
1306 endian (e
, &l
, SFmode
);
1310 /* Convert X to a decimal ASCII string S for output to an assembly
1311 language file. Note, there is no standard way to spell infinity or
1312 a NaN, so these values may require special treatment in the tm.h
1316 ereal_to_decimal (x
, s
)
1326 /* Compare X and Y. Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
1327 or -2 if either is a NaN. */
1331 REAL_VALUE_TYPE x
, y
;
1333 UEMUSHORT ex
[NE
], ey
[NE
];
1337 return (ecmp (ex
, ey
));
1340 /* Return 1 if the sign bit of X is set, else return 0. */
1349 return (eisneg (ex
));
1352 /* End of REAL_ARITHMETIC interface */
1355 Extended precision IEEE binary floating point arithmetic routines
1357 Numbers are stored in C language as arrays of 16-bit unsigned
1358 short integers. The arguments of the routines are pointers to
1361 External e type data structure, similar to Intel 8087 chip
1362 temporary real format but possibly with a larger significand:
1364 NE-1 significand words (least significant word first,
1365 most significant bit is normally set)
1366 exponent (value = EXONE for 1.0,
1367 top bit is the sign)
1370 Internal exploded e-type data structure of a number (a "word" is 16 bits):
1372 ei[0] sign word (0 for positive, 0xffff for negative)
1373 ei[1] biased exponent (value = EXONE for the number 1.0)
1374 ei[2] high guard word (always zero after normalization)
1376 to ei[NI-2] significand (NI-4 significand words,
1377 most significant word first,
1378 most significant bit is set)
1379 ei[NI-1] low guard word (0x8000 bit is rounding place)
1383 Routines for external format e-type numbers
1385 asctoe (string, e) ASCII string to extended double e type
1386 asctoe64 (string, &d) ASCII string to long double
1387 asctoe53 (string, &d) ASCII string to double
1388 asctoe24 (string, &f) ASCII string to single
1389 asctoeg (string, e, prec) ASCII string to specified precision
1390 e24toe (&f, e) IEEE single precision to e type
1391 e53toe (&d, e) IEEE double precision to e type
1392 e64toe (&d, e) IEEE long double precision to e type
1393 e113toe (&d, e) 128-bit long double precision to e type
1395 eabs (e) absolute value
1397 eadd (a, b, c) c = b + a
1399 ecmp (a, b) Returns 1 if a > b, 0 if a == b,
1400 -1 if a < b, -2 if either a or b is a NaN.
1401 ediv (a, b, c) c = b / a
1402 efloor (a, b) truncate to integer, toward -infinity
1403 efrexp (a, exp, s) extract exponent and significand
1404 eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
1405 euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
1406 einfin (e) set e to infinity, leaving its sign alone
1407 eldexp (a, n, b) multiply by 2**n
1409 emul (a, b, c) c = b * a
1412 eround (a, b) b = nearest integer value to a
1414 esub (a, b, c) c = b - a
1416 e24toasc (&f, str, n) single to ASCII string, n digits after decimal
1417 e53toasc (&d, str, n) double to ASCII string, n digits after decimal
1418 e64toasc (&d, str, n) 80-bit long double to ASCII string
1419 e113toasc (&d, str, n) 128-bit long double to ASCII string
1421 etoasc (e, str, n) e to ASCII string, n digits after decimal
1422 etoe24 (e, &f) convert e type to IEEE single precision
1423 etoe53 (e, &d) convert e type to IEEE double precision
1424 etoe64 (e, &d) convert e type to IEEE long double precision
1425 ltoe (&l, e) HOST_WIDE_INT to e type
1426 ultoe (&l, e) unsigned HOST_WIDE_INT to e type
1427 eisneg (e) 1 if sign bit of e != 0, else 0
1428 eisinf (e) 1 if e has maximum exponent (non-IEEE)
1429 or is infinite (IEEE)
1430 eisnan (e) 1 if e is a NaN
1433 Routines for internal format exploded e-type numbers
1435 eaddm (ai, bi) add significands, bi = bi + ai
1437 ecleazs (ei) set ei = 0 but leave its sign alone
1438 ecmpm (ai, bi) compare significands, return 1, 0, or -1
1439 edivm (ai, bi) divide significands, bi = bi / ai
1440 emdnorm (ai,l,s,exp) normalize and round off
1441 emovi (a, ai) convert external a to internal ai
1442 emovo (ai, a) convert internal ai to external a
1443 emovz (ai, bi) bi = ai, low guard word of bi = 0
1444 emulm (ai, bi) multiply significands, bi = bi * ai
1445 enormlz (ei) left-justify the significand
1446 eshdn1 (ai) shift significand and guards down 1 bit
1447 eshdn8 (ai) shift down 8 bits
1448 eshdn6 (ai) shift down 16 bits
1449 eshift (ai, n) shift ai n bits up (or down if n < 0)
1450 eshup1 (ai) shift significand and guards up 1 bit
1451 eshup8 (ai) shift up 8 bits
1452 eshup6 (ai) shift up 16 bits
1453 esubm (ai, bi) subtract significands, bi = bi - ai
1454 eiisinf (ai) 1 if infinite
1455 eiisnan (ai) 1 if a NaN
1456 eiisneg (ai) 1 if sign bit of ai != 0, else 0
1457 einan (ai) set ai = NaN
1459 eiinfin (ai) set ai = infinity
1462 The result is always normalized and rounded to NI-4 word precision
1463 after each arithmetic operation.
1465 Exception flags are NOT fully supported.
1467 Signaling NaN's are NOT supported; they are treated the same
1470 Define INFINITY for support of infinity; otherwise a
1471 saturation arithmetic is implemented.
1473 Define NANS for support of Not-a-Number items; otherwise the
1474 arithmetic will never produce a NaN output, and might be confused
1476 If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1477 either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1478 may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1481 Denormals are always supported here where appropriate (e.g., not
1482 for conversion to DEC numbers). */
1484 /* Definitions for error codes that are passed to the common error handling
1487 For Digital Equipment PDP-11 and VAX computers, certain
1488 IBM systems, and others that use numbers with a 56-bit
1489 significand, the symbol DEC should be defined. In this
1490 mode, most floating point constants are given as arrays
1491 of octal integers to eliminate decimal to binary conversion
1492 errors that might be introduced by the compiler.
1494 For computers, such as IBM PC, that follow the IEEE
1495 Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1496 Std 754-1985), the symbol IEEE should be defined.
1497 These numbers have 53-bit significands. In this mode, constants
1498 are provided as arrays of hexadecimal 16 bit integers.
1499 The endian-ness of generated values is controlled by
1500 REAL_WORDS_BIG_ENDIAN.
1502 To accommodate other types of computer arithmetic, all
1503 constants are also provided in a normal decimal radix
1504 which one can hope are correctly converted to a suitable
1505 format by the available C language compiler. To invoke
1506 this mode, the symbol UNK is defined.
1508 An important difference among these modes is a predefined
1509 set of machine arithmetic constants for each. The numbers
1510 MACHEP (the machine roundoff error), MAXNUM (largest number
1511 represented), and several other parameters are preset by
1512 the configuration symbol. Check the file const.c to
1513 ensure that these values are correct for your computer.
1515 For ANSI C compatibility, define ANSIC equal to 1. Currently
1516 this affects only the atan2 function and others that use it. */
1518 /* Constant definitions for math error conditions. */
1520 #define DOMAIN 1 /* argument domain error */
1521 #define SING 2 /* argument singularity */
1522 #define OVERFLOW 3 /* overflow range error */
1523 #define UNDERFLOW 4 /* underflow range error */
1524 #define TLOSS 5 /* total loss of precision */
1525 #define PLOSS 6 /* partial loss of precision */
1526 #define INVALID 7 /* NaN-producing operation */
1528 /* e type constants used by high precision check routines */
1530 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
1532 UEMUSHORT ezero
[NE
] =
1533 {0x0000, 0x0000, 0x0000, 0x0000,
1534 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
1535 extern UEMUSHORT ezero
[];
1538 UEMUSHORT ehalf
[NE
] =
1539 {0x0000, 0x0000, 0x0000, 0x0000,
1540 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
1541 extern UEMUSHORT ehalf
[];
1544 UEMUSHORT eone
[NE
] =
1545 {0x0000, 0x0000, 0x0000, 0x0000,
1546 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
1547 extern UEMUSHORT eone
[];
1550 UEMUSHORT etwo
[NE
] =
1551 {0x0000, 0x0000, 0x0000, 0x0000,
1552 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
1553 extern UEMUSHORT etwo
[];
1557 {0x0000, 0x0000, 0x0000, 0x0000,
1558 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
1559 extern UEMUSHORT e32
[];
1561 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1562 UEMUSHORT elog2
[NE
] =
1563 {0x40f3, 0xf6af, 0x03f2, 0xb398,
1564 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1565 extern UEMUSHORT elog2
[];
1567 /* 1.41421356237309504880168872420969807856967187537695E0 */
1568 UEMUSHORT esqrt2
[NE
] =
1569 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1570 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1571 extern UEMUSHORT esqrt2
[];
1573 /* 3.14159265358979323846264338327950288419716939937511E0 */
1575 {0x2902, 0x1cd1, 0x80dc, 0x628b,
1576 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1577 extern UEMUSHORT epi
[];
1580 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1581 UEMUSHORT ezero
[NE
] =
1582 {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1583 UEMUSHORT ehalf
[NE
] =
1584 {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1585 UEMUSHORT eone
[NE
] =
1586 {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1587 UEMUSHORT etwo
[NE
] =
1588 {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1590 {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1591 UEMUSHORT elog2
[NE
] =
1592 {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1593 UEMUSHORT esqrt2
[NE
] =
1594 {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1596 {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1599 /* Control register for rounding precision.
1600 This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits. */
1605 /* Clear out entire e-type number X. */
1609 register UEMUSHORT
*x
;
1613 for (i
= 0; i
< NE
; i
++)
1617 /* Move e-type number from A to B. */
1621 register UEMUSHORT
*a
, *b
;
1625 for (i
= 0; i
< NE
; i
++)
1631 /* Absolute value of e-type X. */
1637 /* sign is top bit of last word of external format */
1638 x
[NE
- 1] &= 0x7fff;
1642 /* Negate the e-type number X. */
1649 x
[NE
- 1] ^= 0x8000; /* Toggle the sign bit */
1652 /* Return 1 if sign bit of e-type number X is nonzero, else zero. */
1659 if (x
[NE
- 1] & 0x8000)
1665 /* Return 1 if e-type number X is infinity, else return zero. */
1676 if ((x
[NE
- 1] & 0x7fff) == 0x7fff)
1682 /* Check if e-type number is not a number. The bit pattern is one that we
1683 defined, so we know for sure how to detect it. */
1687 UEMUSHORT x
[] ATTRIBUTE_UNUSED
;
1692 /* NaN has maximum exponent */
1693 if ((x
[NE
- 1] & 0x7fff) != 0x7fff)
1695 /* ... and non-zero significand field. */
1696 for (i
= 0; i
< NE
- 1; i
++)
1706 /* Fill e-type number X with infinity pattern (IEEE)
1707 or largest possible number (non-IEEE). */
1711 register UEMUSHORT
*x
;
1716 for (i
= 0; i
< NE
- 1; i
++)
1720 for (i
= 0; i
< NE
- 1; i
++)
1748 /* Output an e-type NaN.
1749 This generates Intel's quiet NaN pattern for extended real.
1750 The exponent is 7fff, the leading mantissa word is c000. */
1755 register UEMUSHORT
*x
;
1760 for (i
= 0; i
< NE
- 2; i
++)
1763 *x
= (sign
<< 15) | 0x7fff;
1767 /* Move in an e-type number A, converting it to exploded e-type B. */
1773 register UEMUSHORT
*p
, *q
;
1777 p
= a
+ (NE
- 1); /* point to last word of external number */
1778 /* get the sign bit */
1783 /* get the exponent */
1785 *q
++ &= 0x7fff; /* delete the sign bit */
1787 if ((*(q
- 1) & 0x7fff) == 0x7fff)
1793 for (i
= 3; i
< NI
; i
++)
1799 for (i
= 2; i
< NI
; i
++)
1805 /* clear high guard word */
1807 /* move in the significand */
1808 for (i
= 0; i
< NE
- 1; i
++)
1810 /* clear low guard word */
1814 /* Move out exploded e-type number A, converting it to e type B. */
1820 register UEMUSHORT
*p
, *q
;
1825 q
= b
+ (NE
- 1); /* point to output exponent */
1826 /* combine sign and exponent */
1829 *q
-- = *p
++ | 0x8000;
1833 if (*(p
- 1) == 0x7fff)
1838 enan (b
, eiisneg (a
));
1846 /* skip over guard word */
1848 /* move the significand */
1849 for (j
= 0; j
< NE
- 1; j
++)
1853 /* Clear out exploded e-type number XI. */
1857 register UEMUSHORT
*xi
;
1861 for (i
= 0; i
< NI
; i
++)
1865 /* Clear out exploded e-type XI, but don't touch the sign. */
1869 register UEMUSHORT
*xi
;
1874 for (i
= 0; i
< NI
- 1; i
++)
1878 /* Move exploded e-type number from A to B. */
1882 register UEMUSHORT
*a
, *b
;
1886 for (i
= 0; i
< NI
- 1; i
++)
1888 /* clear low guard word */
1892 /* Generate exploded e-type NaN.
1893 The explicit pattern for this is maximum exponent and
1894 top two significant bits set. */
1908 /* Return nonzero if exploded e-type X is a NaN. */
1917 if ((x
[E
] & 0x7fff) == 0x7fff)
1919 for (i
= M
+ 1; i
< NI
; i
++)
1929 /* Return nonzero if sign of exploded e-type X is nonzero. */
1942 /* Fill exploded e-type X with infinity pattern.
1943 This has maximum exponent and significand all zeros. */
1955 /* Return nonzero if exploded e-type X is infinite. */
1967 if ((x
[E
] & 0x7fff) == 0x7fff)
1971 #endif /* INFINITY */
1973 /* Compare significands of numbers in internal exploded e-type format.
1974 Guard words are included in the comparison.
1982 register UEMUSHORT
*a
, *b
;
1986 a
+= M
; /* skip up to significand area */
1988 for (i
= M
; i
< NI
; i
++)
1996 if (*(--a
) > *(--b
))
2002 /* Shift significand of exploded e-type X down by 1 bit. */
2006 register UEMUSHORT
*x
;
2008 register UEMUSHORT bits
;
2011 x
+= M
; /* point to significand area */
2014 for (i
= M
; i
< NI
; i
++)
2026 /* Shift significand of exploded e-type X up by 1 bit. */
2030 register UEMUSHORT
*x
;
2032 register UEMUSHORT bits
;
2038 for (i
= M
; i
< NI
; i
++)
2051 /* Shift significand of exploded e-type X down by 8 bits. */
2055 register UEMUSHORT
*x
;
2057 register UEMUSHORT newbyt
, oldbyt
;
2062 for (i
= M
; i
< NI
; i
++)
2072 /* Shift significand of exploded e-type X up by 8 bits. */
2076 register UEMUSHORT
*x
;
2079 register UEMUSHORT newbyt
, oldbyt
;
2084 for (i
= M
; i
< NI
; i
++)
2094 /* Shift significand of exploded e-type X up by 16 bits. */
2098 register UEMUSHORT
*x
;
2101 register UEMUSHORT
*p
;
2106 for (i
= M
; i
< NI
- 1; i
++)
2112 /* Shift significand of exploded e-type X down by 16 bits. */
2116 register UEMUSHORT
*x
;
2119 register UEMUSHORT
*p
;
2124 for (i
= M
; i
< NI
- 1; i
++)
2130 /* Add significands of exploded e-type X and Y. X + Y replaces Y. */
2136 register unsigned EMULONG a
;
2143 for (i
= M
; i
< NI
; i
++)
2145 a
= (unsigned EMULONG
) (*x
) + (unsigned EMULONG
) (*y
) + carry
;
2156 /* Subtract significands of exploded e-type X and Y. Y - X replaces Y. */
2169 for (i
= M
; i
< NI
; i
++)
2171 a
= (unsigned EMULONG
) (*y
) - (unsigned EMULONG
) (*x
) - carry
;
2183 static UEMUSHORT equot
[NI
];
2187 /* Radix 2 shift-and-add versions of multiply and divide */
2190 /* Divide significands */
2194 UEMUSHORT den
[], num
[];
2197 register UEMUSHORT
*p
, *q
;
2204 for (i
= M
; i
< NI
; i
++)
2209 /* Use faster compare and subtraction if denominator has only 15 bits of
2215 for (i
= M
+ 3; i
< NI
; i
++)
2220 if ((den
[M
+ 1] & 1) != 0)
2228 for (i
= 0; i
< NBITS
+ 2; i
++)
2246 /* The number of quotient bits to calculate is NBITS + 1 scaling guard
2247 bit + 1 roundoff bit. */
2252 for (i
= 0; i
< NBITS
+ 2; i
++)
2254 if (ecmpm (den
, num
) <= 0)
2257 j
= 1; /* quotient bit = 1 */
2271 /* test for nonzero remainder after roundoff bit */
2274 for (i
= M
; i
< NI
; i
++)
2282 for (i
= 0; i
< NI
; i
++)
2288 /* Multiply significands */
2299 for (i
= M
; i
< NI
; i
++)
2304 while (*p
== 0) /* significand is not supposed to be zero */
2309 if ((*p
& 0xff) == 0)
2317 for (i
= 0; i
< k
; i
++)
2321 /* remember if there were any nonzero bits shifted out */
2328 for (i
= 0; i
< NI
; i
++)
2331 /* return flag for lost nonzero bits */
2337 /* Radix 65536 versions of multiply and divide. */
2339 /* Multiply significand of e-type number B
2340 by 16-bit quantity A, return e-type result to C. */
2347 register UEMUSHORT
*pp
;
2348 register unsigned EMULONG carry
;
2351 unsigned EMULONG aa
, m
;
2360 for (i
=M
+1; i
<NI
; i
++)
2370 m
= (unsigned EMULONG
) aa
* *ps
--;
2371 carry
= (m
& 0xffff) + *pp
;
2372 *pp
-- = (UEMUSHORT
)carry
;
2373 carry
= (carry
>> 16) + (m
>> 16) + *pp
;
2374 *pp
= (UEMUSHORT
)carry
;
2375 *(pp
-1) = carry
>> 16;
2378 for (i
=M
; i
<NI
; i
++)
2382 /* Divide significands of exploded e-types NUM / DEN. Neither the
2383 numerator NUM nor the denominator DEN is permitted to have its high guard
2388 UEMUSHORT den
[], num
[];
2391 register UEMUSHORT
*p
;
2392 unsigned EMULONG tnum
;
2393 UEMUSHORT j
, tdenm
, tquot
;
2394 UEMUSHORT tprod
[NI
+1];
2400 for (i
=M
; i
<NI
; i
++)
2406 for (i
=M
; i
<NI
; i
++)
2408 /* Find trial quotient digit (the radix is 65536). */
2409 tnum
= (((unsigned EMULONG
) num
[M
]) << 16) + num
[M
+1];
2411 /* Do not execute the divide instruction if it will overflow. */
2412 if ((tdenm
* (unsigned long)0xffff) < tnum
)
2415 tquot
= tnum
/ tdenm
;
2416 /* Multiply denominator by trial quotient digit. */
2417 m16m ((unsigned int)tquot
, den
, tprod
);
2418 /* The quotient digit may have been overestimated. */
2419 if (ecmpm (tprod
, num
) > 0)
2423 if (ecmpm (tprod
, num
) > 0)
2433 /* test for nonzero remainder after roundoff bit */
2436 for (i
=M
; i
<NI
; i
++)
2443 for (i
=0; i
<NI
; i
++)
2449 /* Multiply significands of exploded e-type A and B, result in B. */
2456 UEMUSHORT pprod
[NI
];
2462 for (i
=M
; i
<NI
; i
++)
2468 for (i
=M
+1; i
<NI
; i
++)
2476 m16m ((unsigned int) *p
--, b
, pprod
);
2477 eaddm(pprod
, equot
);
2483 for (i
=0; i
<NI
; i
++)
2486 /* return flag for lost nonzero bits */
2492 /* Normalize and round off.
2494 The internal format number to be rounded is S.
2495 Input LOST is 0 if the value is exact. This is the so-called sticky bit.
2497 Input SUBFLG indicates whether the number was obtained
2498 by a subtraction operation. In that case if LOST is nonzero
2499 then the number is slightly smaller than indicated.
2501 Input EXP is the biased exponent, which may be negative.
2502 the exponent field of S is ignored but is replaced by
2503 EXP as adjusted by normalization and rounding.
2505 Input RCNTRL is the rounding control. If it is nonzero, the
2506 returned value will be rounded to RNDPRC bits.
2508 For future reference: In order for emdnorm to round off denormal
2509 significands at the right point, the input exponent must be
2510 adjusted to be the actual value it would have after conversion to
2511 the final floating point type. This adjustment has been
2512 implemented for all type conversions (etoe53, etc.) and decimal
2513 conversions, but not for the arithmetic functions (eadd, etc.).
2514 Data types having standard 15-bit exponents are not affected by
2515 this, but SFmode and DFmode are affected. For example, ediv with
2516 rndprc = 24 will not round correctly to 24-bit precision if the
2517 result is denormal. */
2519 static int rlast
= -1;
2521 static UEMUSHORT rmsk
= 0;
2522 static UEMUSHORT rmbit
= 0;
2523 static UEMUSHORT rebit
= 0;
2525 static UEMUSHORT rbit
[NI
];
2528 emdnorm (s
, lost
, subflg
, exp
, rcntrl
)
2541 /* a blank significand could mean either zero or infinity. */
2554 if ((j
> NBITS
) && (exp
< 32767))
2562 if (exp
> (EMULONG
) (-NBITS
- 1))
2575 /* Round off, unless told not to by rcntrl. */
2578 /* Set up rounding parameters if the control register changed. */
2579 if (rndprc
!= rlast
)
2586 rw
= NI
- 1; /* low guard word */
2609 /* For DEC or IBM arithmetic */
2626 /* For C4x arithmetic */
2647 /* Shift down 1 temporarily if the data structure has an implied
2648 most significant bit and the number is denormal.
2649 Intel long double denormals also lose one bit of precision. */
2650 if ((exp
<= 0) && (rndprc
!= NBITS
)
2651 && ((rndprc
!= 64) || ((rndprc
== 64) && ! REAL_WORDS_BIG_ENDIAN
)))
2653 lost
|= s
[NI
- 1] & 1;
2656 /* Clear out all bits below the rounding bit,
2657 remembering in r if any were nonzero. */
2671 if ((r
& rmbit
) != 0)
2677 { /* round to even */
2678 if ((s
[re
] & rebit
) == 0)
2691 /* Undo the temporary shift for denormal values. */
2692 if ((exp
<= 0) && (rndprc
!= NBITS
)
2693 && ((rndprc
!= 64) || ((rndprc
== 64) && ! REAL_WORDS_BIG_ENDIAN
)))
2698 { /* overflow on roundoff */
2711 for (i
= 2; i
< NI
- 1; i
++)
2714 warning ("floating point overflow");
2718 for (i
= M
+ 1; i
< NI
- 1; i
++)
2721 if ((rndprc
< 64) || (rndprc
== 113))
2736 s
[1] = (UEMUSHORT
) exp
;
2739 /* Subtract. C = B - A, all e type numbers. */
2741 static int subflg
= 0;
2745 UEMUSHORT
*a
, *b
, *c
;
2759 /* Infinity minus infinity is a NaN.
2760 Test for subtracting infinities of the same sign. */
2761 if (eisinf (a
) && eisinf (b
)
2762 && ((eisneg (a
) ^ eisneg (b
)) == 0))
2764 mtherr ("esub", INVALID
);
2773 /* Add. C = A + B, all e type. */
2777 UEMUSHORT
*a
, *b
, *c
;
2781 /* NaN plus anything is a NaN. */
2792 /* Infinity minus infinity is a NaN.
2793 Test for adding infinities of opposite signs. */
2794 if (eisinf (a
) && eisinf (b
)
2795 && ((eisneg (a
) ^ eisneg (b
)) != 0))
2797 mtherr ("esub", INVALID
);
2806 /* Arithmetic common to both addition and subtraction. */
2810 UEMUSHORT
*a
, *b
, *c
;
2812 UEMUSHORT ai
[NI
], bi
[NI
], ci
[NI
];
2814 EMULONG lt
, lta
, ltb
;
2835 /* compare exponents */
2840 { /* put the larger number in bi */
2850 if (lt
< (EMULONG
) (-NBITS
- 1))
2851 goto done
; /* answer same as larger addend */
2853 lost
= eshift (ai
, k
); /* shift the smaller number down */
2857 /* exponents were the same, so must compare significands */
2860 { /* the numbers are identical in magnitude */
2861 /* if different signs, result is zero */
2867 /* if same sign, result is double */
2868 /* double denormalized tiny number */
2869 if ((bi
[E
] == 0) && ((bi
[3] & 0x8000) == 0))
2874 /* add 1 to exponent unless both are zero! */
2875 for (j
= 1; j
< NI
- 1; j
++)
2891 bi
[E
] = (UEMUSHORT
) ltb
;
2895 { /* put the larger number in bi */
2911 emdnorm (bi
, lost
, subflg
, ltb
, 64);
2917 /* Divide: C = B/A, all e type. */
2921 UEMUSHORT
*a
, *b
, *c
;
2923 UEMUSHORT ai
[NI
], bi
[NI
];
2925 EMULONG lt
, lta
, ltb
;
2927 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2928 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
2929 sign
= eisneg(a
) ^ eisneg(b
);
2932 /* Return any NaN input. */
2943 /* Zero over zero, or infinity over infinity, is a NaN. */
2944 if (((ecmp (a
, ezero
) == 0) && (ecmp (b
, ezero
) == 0))
2945 || (eisinf (a
) && eisinf (b
)))
2947 mtherr ("ediv", INVALID
);
2952 /* Infinity over anything else is infinity. */
2959 /* Anything else over infinity is zero. */
2971 { /* See if numerator is zero. */
2972 for (i
= 1; i
< NI
- 1; i
++)
2976 ltb
-= enormlz (bi
);
2986 { /* possible divide by zero */
2987 for (i
= 1; i
< NI
- 1; i
++)
2991 lta
-= enormlz (ai
);
2995 /* Divide by zero is not an invalid operation.
2996 It is a divide-by-zero operation! */
2998 mtherr ("ediv", SING
);
3004 /* calculate exponent */
3005 lt
= ltb
- lta
+ EXONE
;
3006 emdnorm (bi
, i
, 0, lt
, 64);
3013 && (ecmp (c
, ezero
) != 0)
3016 *(c
+(NE
-1)) |= 0x8000;
3018 *(c
+(NE
-1)) &= ~0x8000;
3021 /* Multiply e-types A and B, return e-type product C. */
3025 UEMUSHORT
*a
, *b
, *c
;
3027 UEMUSHORT ai
[NI
], bi
[NI
];
3029 EMULONG lt
, lta
, ltb
;
3031 /* IEEE says if result is not a NaN, the sign is "-" if and only if
3032 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
3033 sign
= eisneg(a
) ^ eisneg(b
);
3036 /* NaN times anything is the same NaN. */
3047 /* Zero times infinity is a NaN. */
3048 if ((eisinf (a
) && (ecmp (b
, ezero
) == 0))
3049 || (eisinf (b
) && (ecmp (a
, ezero
) == 0)))
3051 mtherr ("emul", INVALID
);
3056 /* Infinity times anything else is infinity. */
3058 if (eisinf (a
) || eisinf (b
))
3070 for (i
= 1; i
< NI
- 1; i
++)
3074 lta
-= enormlz (ai
);
3085 for (i
= 1; i
< NI
- 1; i
++)
3089 ltb
-= enormlz (bi
);
3098 /* Multiply significands */
3100 /* calculate exponent */
3101 lt
= lta
+ ltb
- (EXONE
- 1);
3102 emdnorm (bi
, j
, 0, lt
, 64);
3109 && (ecmp (c
, ezero
) != 0)
3112 *(c
+(NE
-1)) |= 0x8000;
3114 *(c
+(NE
-1)) &= ~0x8000;
3117 /* Convert double precision PE to e-type Y. */
3130 ibmtoe (pe
, y
, DFmode
);
3135 c4xtoe (pe
, y
, HFmode
);
3138 register UEMUSHORT r
;
3139 register UEMUSHORT
*e
, *p
;
3144 denorm
= 0; /* flag if denormalized number */
3146 if (! REAL_WORDS_BIG_ENDIAN
)
3152 yy
[M
] = (r
& 0x0f) | 0x10;
3153 r
&= ~0x800f; /* strip sign and 4 significand bits */
3158 if (! REAL_WORDS_BIG_ENDIAN
)
3160 if (((pe
[3] & 0xf) != 0) || (pe
[2] != 0)
3161 || (pe
[1] != 0) || (pe
[0] != 0))
3163 enan (y
, yy
[0] != 0);
3169 if (((pe
[0] & 0xf) != 0) || (pe
[1] != 0)
3170 || (pe
[2] != 0) || (pe
[3] != 0))
3172 enan (y
, yy
[0] != 0);
3183 #endif /* INFINITY */
3185 /* If zero exponent, then the significand is denormalized.
3186 So take back the understood high significand bit. */
3197 if (! REAL_WORDS_BIG_ENDIAN
)
3214 /* If zero exponent, then normalize the significand. */
3215 if ((k
= enormlz (yy
)) > NBITS
)
3218 yy
[E
] -= (UEMUSHORT
) (k
- 1);
3221 #endif /* not C4X */
3222 #endif /* not IBM */
3223 #endif /* not DEC */
3226 /* Convert double extended precision float PE to e type Y. */
3233 UEMUSHORT
*e
, *p
, *q
;
3238 for (i
= 0; i
< NE
- 5; i
++)
3240 /* This precision is not ordinarily supported on DEC or IBM. */
3242 for (i
= 0; i
< 5; i
++)
3246 p
= &yy
[0] + (NE
- 1);
3249 for (i
= 0; i
< 5; i
++)
3253 if (! REAL_WORDS_BIG_ENDIAN
)
3255 for (i
= 0; i
< 5; i
++)
3258 /* For denormal long double Intel format, shift significand up one
3259 -- but only if the top significand bit is zero. A top bit of 1
3260 is "pseudodenormal" when the exponent is zero. */
3261 if((yy
[NE
-1] & 0x7fff) == 0 && (yy
[NE
-2] & 0x8000) == 0)
3273 p
= &yy
[0] + (NE
- 1);
3274 #ifdef ARM_EXTENDED_IEEE_FORMAT
3275 /* For ARMs, the exponent is in the lowest 15 bits of the word. */
3276 *p
-- = (e
[0] & 0x8000) | (e
[1] & 0x7ffff);
3282 for (i
= 0; i
< 4; i
++)
3287 /* Point to the exponent field and check max exponent cases. */
3289 if ((*p
& 0x7fff) == 0x7fff)
3292 if (! REAL_WORDS_BIG_ENDIAN
)
3294 for (i
= 0; i
< 4; i
++)
3296 if ((i
!= 3 && pe
[i
] != 0)
3297 /* Anything but 0x8000 here, including 0, is a NaN. */
3298 || (i
== 3 && pe
[i
] != 0x8000))
3300 enan (y
, (*p
& 0x8000) != 0);
3307 #ifdef ARM_EXTENDED_IEEE_FORMAT
3308 for (i
= 2; i
<= 5; i
++)
3312 enan (y
, (*p
& 0x8000) != 0);
3317 /* In Motorola extended precision format, the most significant
3318 bit of an infinity mantissa could be either 1 or 0. It is
3319 the lower order bits that tell whether the value is a NaN. */
3320 if ((pe
[2] & 0x7fff) != 0)
3323 for (i
= 3; i
<= 5; i
++)
3328 enan (y
, (*p
& 0x8000) != 0);
3332 #endif /* not ARM */
3341 #endif /* INFINITY */
3344 for (i
= 0; i
< NE
; i
++)
3348 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3349 /* Convert 128-bit long double precision float PE to e type Y. */
3355 register UEMUSHORT r
;
3364 if (! REAL_WORDS_BIG_ENDIAN
)
3376 if (! REAL_WORDS_BIG_ENDIAN
)
3378 for (i
= 0; i
< 7; i
++)
3382 enan (y
, yy
[0] != 0);
3389 for (i
= 1; i
< 8; i
++)
3393 enan (y
, yy
[0] != 0);
3405 #endif /* INFINITY */
3409 if (! REAL_WORDS_BIG_ENDIAN
)
3411 for (i
= 0; i
< 7; i
++)
3417 for (i
= 0; i
< 7; i
++)
3421 /* If denormal, remove the implied bit; else shift down 1. */
3435 /* Convert single precision float PE to e type Y. */
3443 ibmtoe (pe
, y
, SFmode
);
3449 c4xtoe (pe
, y
, QFmode
);
3453 register UEMUSHORT r
;
3454 register UEMUSHORT
*e
, *p
;
3459 denorm
= 0; /* flag if denormalized number */
3462 if (! REAL_WORDS_BIG_ENDIAN
)
3472 yy
[M
] = (r
& 0x7f) | 0200;
3473 r
&= ~0x807f; /* strip sign and 7 significand bits */
3478 if (REAL_WORDS_BIG_ENDIAN
)
3480 if (((pe
[0] & 0x7f) != 0) || (pe
[1] != 0))
3482 enan (y
, yy
[0] != 0);
3488 if (((pe
[1] & 0x7f) != 0) || (pe
[0] != 0))
3490 enan (y
, yy
[0] != 0);
3501 #endif /* INFINITY */
3503 /* If zero exponent, then the significand is denormalized.
3504 So take back the understood high significand bit. */
3517 if (! REAL_WORDS_BIG_ENDIAN
)
3527 { /* if zero exponent, then normalize the significand */
3528 if ((k
= enormlz (yy
)) > NBITS
)
3531 yy
[E
] -= (UEMUSHORT
) (k
- 1);
3534 #endif /* not C4X */
3535 #endif /* not IBM */
3538 /* Convert e-type X to IEEE 128-bit long double format E. */
3551 make_nan (e
, eisneg (x
), TFmode
);
3556 exp
= (EMULONG
) xi
[E
];
3561 /* round off to nearest or even */
3564 emdnorm (xi
, 0, 0, exp
, 64);
3572 /* Convert exploded e-type X, that has already been rounded to
3573 113-bit precision, to IEEE 128-bit long double format Y. */
3579 register UEMUSHORT
*p
, *q
;
3585 make_nan (b
, eiisneg (a
), TFmode
);
3590 if (REAL_WORDS_BIG_ENDIAN
)
3593 q
= b
+ 7; /* point to output exponent */
3595 /* If not denormal, delete the implied bit. */
3600 /* combine sign and exponent */
3602 if (REAL_WORDS_BIG_ENDIAN
)
3605 *q
++ = *p
++ | 0x8000;
3612 *q
-- = *p
++ | 0x8000;
3616 /* skip over guard word */
3618 /* move the significand */
3619 if (REAL_WORDS_BIG_ENDIAN
)
3621 for (i
= 0; i
< 7; i
++)
3626 for (i
= 0; i
< 7; i
++)
3631 /* Convert e-type X to IEEE double extended format E. */
3644 make_nan (e
, eisneg (x
), XFmode
);
3649 /* adjust exponent for offset */
3650 exp
= (EMULONG
) xi
[E
];
3655 /* round off to nearest or even */
3658 emdnorm (xi
, 0, 0, exp
, 64);
3666 /* Convert exploded e-type X, that has already been rounded to
3667 64-bit precision, to IEEE double extended format Y. */
3673 register UEMUSHORT
*p
, *q
;
3679 make_nan (b
, eiisneg (a
), XFmode
);
3683 /* Shift denormal long double Intel format significand down one bit. */
3684 if ((a
[E
] == 0) && ! REAL_WORDS_BIG_ENDIAN
)
3694 if (REAL_WORDS_BIG_ENDIAN
)
3698 q
= b
+ 4; /* point to output exponent */
3699 /* Clear the last two bytes of 12-byte Intel format. q is pointing
3700 into an array of size 6 (e.g. x[NE]), so the last two bytes are
3701 always there, and there are never more bytes, even when we are using
3702 INTEL_EXTENDED_IEEE_FORMAT. */
3707 /* combine sign and exponent */
3711 *q
++ = *p
++ | 0x8000;
3718 *q
-- = *p
++ | 0x8000;
3723 if (REAL_WORDS_BIG_ENDIAN
)
3725 #ifdef ARM_EXTENDED_IEEE_FORMAT
3726 /* The exponent is in the lowest 15 bits of the first word. */
3727 *q
++ = i
? 0x8000 : 0;
3731 *q
++ = *p
++ | 0x8000;
3740 *q
-- = *p
++ | 0x8000;
3745 /* skip over guard word */
3747 /* move the significand */
3749 for (i
= 0; i
< 4; i
++)
3753 for (i
= 0; i
< 4; i
++)
3757 if (REAL_WORDS_BIG_ENDIAN
)
3759 for (i
= 0; i
< 4; i
++)
3767 /* Intel long double infinity significand. */
3775 for (i
= 0; i
< 4; i
++)
3781 /* e type to double precision. */
3784 /* Convert e-type X to DEC-format double E. */
3790 etodec (x
, e
); /* see etodec.c */
3793 /* Convert exploded e-type X, that has already been rounded to
3794 56-bit double precision, to DEC double Y. */
3805 /* Convert e-type X to IBM 370-format double E. */
3811 etoibm (x
, e
, DFmode
);
3814 /* Convert exploded e-type X, that has already been rounded to
3815 56-bit precision, to IBM 370 double Y. */
3821 toibm (x
, y
, DFmode
);
3824 #else /* it's neither DEC nor IBM */
3826 /* Convert e-type X to C4X-format long double E. */
3832 etoc4x (x
, e
, HFmode
);
3835 /* Convert exploded e-type X, that has already been rounded to
3836 56-bit precision, to IBM 370 double Y. */
3842 toc4x (x
, y
, HFmode
);
3845 #else /* it's neither DEC nor IBM nor C4X */
3847 /* Convert e-type X to IEEE double E. */
3860 make_nan (e
, eisneg (x
), DFmode
);
3865 /* adjust exponent for offsets */
3866 exp
= (EMULONG
) xi
[E
] - (EXONE
- 0x3ff);
3871 /* round off to nearest or even */
3874 emdnorm (xi
, 0, 0, exp
, 64);
3882 /* Convert exploded e-type X, that has already been rounded to
3883 53-bit precision, to IEEE double Y. */
3895 make_nan (y
, eiisneg (x
), DFmode
);
3901 if (! REAL_WORDS_BIG_ENDIAN
)
3904 *y
= 0; /* output high order */
3906 *y
= 0x8000; /* output sign bit */
3909 if (i
>= (unsigned int) 2047)
3911 /* Saturate at largest number less than infinity. */
3914 if (! REAL_WORDS_BIG_ENDIAN
)
3928 *y
|= (UEMUSHORT
) 0x7fef;
3929 if (! REAL_WORDS_BIG_ENDIAN
)
3954 i
|= *p
++ & (UEMUSHORT
) 0x0f; /* *p = xi[M] */
3955 *y
|= (UEMUSHORT
) i
; /* high order output already has sign bit set */
3956 if (! REAL_WORDS_BIG_ENDIAN
)
3971 #endif /* not C4X */
3972 #endif /* not IBM */
3973 #endif /* not DEC */
3977 /* e type to single precision. */
3980 /* Convert e-type X to IBM 370 float E. */
3986 etoibm (x
, e
, SFmode
);
3989 /* Convert exploded e-type X, that has already been rounded to
3990 float precision, to IBM 370 float Y. */
3996 toibm (x
, y
, SFmode
);
4002 /* Convert e-type X to C4X float E. */
4008 etoc4x (x
, e
, QFmode
);
4011 /* Convert exploded e-type X, that has already been rounded to
4012 float precision, to IBM 370 float Y. */
4018 toc4x (x
, y
, QFmode
);
4023 /* Convert e-type X to IEEE float E. DEC float is the same as IEEE float. */
4036 make_nan (e
, eisneg (x
), SFmode
);
4041 /* adjust exponent for offsets */
4042 exp
= (EMULONG
) xi
[E
] - (EXONE
- 0177);
4047 /* round off to nearest or even */
4050 emdnorm (xi
, 0, 0, exp
, 64);
4058 /* Convert exploded e-type X, that has already been rounded to
4059 float precision, to IEEE float Y. */
4071 make_nan (y
, eiisneg (x
), SFmode
);
4077 if (! REAL_WORDS_BIG_ENDIAN
)
4083 *y
= 0; /* output high order */
4085 *y
= 0x8000; /* output sign bit */
4088 /* Handle overflow cases. */
4092 *y
|= (UEMUSHORT
) 0x7f80;
4097 if (! REAL_WORDS_BIG_ENDIAN
)
4105 #else /* no INFINITY */
4106 *y
|= (UEMUSHORT
) 0x7f7f;
4111 if (! REAL_WORDS_BIG_ENDIAN
)
4122 #endif /* no INFINITY */
4134 i
|= *p
++ & (UEMUSHORT
) 0x7f; /* *p = xi[M] */
4135 /* High order output already has sign bit set. */
4141 if (! REAL_WORDS_BIG_ENDIAN
)
4150 #endif /* not C4X */
4151 #endif /* not IBM */
4153 /* Compare two e type numbers.
4157 -2 if either a or b is a NaN. */
4163 UEMUSHORT ai
[NI
], bi
[NI
];
4164 register UEMUSHORT
*p
, *q
;
4169 if (eisnan (a
) || eisnan (b
))
4178 { /* the signs are different */
4180 for (i
= 1; i
< NI
- 1; i
++)
4194 /* both are the same sign */
4209 return (0); /* equality */
4213 if (*(--p
) > *(--q
))
4214 return (msign
); /* p is bigger */
4216 return (-msign
); /* p is littler */
4220 /* Find e-type nearest integer to X, as floor (X + 0.5). */
4231 /* Convert HOST_WIDE_INT LP to e type Y. */
4239 unsigned HOST_WIDE_INT ll
;
4245 /* make it positive */
4246 ll
= (unsigned HOST_WIDE_INT
) (-(*lp
));
4247 yi
[0] = 0xffff; /* put correct sign in the e type number */
4251 ll
= (unsigned HOST_WIDE_INT
) (*lp
);
4253 /* move the long integer to yi significand area */
4254 #if HOST_BITS_PER_WIDE_INT == 64
4255 yi
[M
] = (UEMUSHORT
) (ll
>> 48);
4256 yi
[M
+ 1] = (UEMUSHORT
) (ll
>> 32);
4257 yi
[M
+ 2] = (UEMUSHORT
) (ll
>> 16);
4258 yi
[M
+ 3] = (UEMUSHORT
) ll
;
4259 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
4261 yi
[M
] = (UEMUSHORT
) (ll
>> 16);
4262 yi
[M
+ 1] = (UEMUSHORT
) ll
;
4263 yi
[E
] = EXONE
+ 15; /* exponent if normalize shift count were 0 */
4266 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
4267 ecleaz (yi
); /* it was zero */
4269 yi
[E
] -= (UEMUSHORT
) k
;/* subtract shift count from exponent */
4270 emovo (yi
, y
); /* output the answer */
4273 /* Convert unsigned HOST_WIDE_INT LP to e type Y. */
4277 unsigned HOST_WIDE_INT
*lp
;
4281 unsigned HOST_WIDE_INT ll
;
4287 /* move the long integer to ayi significand area */
4288 #if HOST_BITS_PER_WIDE_INT == 64
4289 yi
[M
] = (UEMUSHORT
) (ll
>> 48);
4290 yi
[M
+ 1] = (UEMUSHORT
) (ll
>> 32);
4291 yi
[M
+ 2] = (UEMUSHORT
) (ll
>> 16);
4292 yi
[M
+ 3] = (UEMUSHORT
) ll
;
4293 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
4295 yi
[M
] = (UEMUSHORT
) (ll
>> 16);
4296 yi
[M
+ 1] = (UEMUSHORT
) ll
;
4297 yi
[E
] = EXONE
+ 15; /* exponent if normalize shift count were 0 */
4300 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
4301 ecleaz (yi
); /* it was zero */
4303 yi
[E
] -= (UEMUSHORT
) k
; /* subtract shift count from exponent */
4304 emovo (yi
, y
); /* output the answer */
4308 /* Find signed HOST_WIDE_INT integer I and floating point fractional
4309 part FRAC of e-type (packed internal format) floating point input X.
4310 The integer output I has the sign of the input, except that
4311 positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
4312 The output e-type fraction FRAC is the positive fractional
4323 unsigned HOST_WIDE_INT ll
;
4326 k
= (int) xi
[E
] - (EXONE
- 1);
4329 /* if exponent <= 0, integer = 0 and real output is fraction */
4334 if (k
> (HOST_BITS_PER_WIDE_INT
- 1))
4336 /* long integer overflow: output large integer
4337 and correct fraction */
4339 *i
= ((unsigned HOST_WIDE_INT
) 1) << (HOST_BITS_PER_WIDE_INT
- 1);
4342 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
4343 /* In this case, let it overflow and convert as if unsigned. */
4344 euifrac (x
, &ll
, frac
);
4345 *i
= (HOST_WIDE_INT
) ll
;
4348 /* In other cases, return the largest positive integer. */
4349 *i
= (((unsigned HOST_WIDE_INT
) 1) << (HOST_BITS_PER_WIDE_INT
- 1)) - 1;
4354 warning ("overflow on truncation to integer");
4358 /* Shift more than 16 bits: first shift up k-16 mod 16,
4359 then shift up by 16's. */
4360 j
= k
- ((k
>> 4) << 4);
4367 ll
= (ll
<< 16) | xi
[M
];
4369 while ((k
-= 16) > 0);
4376 /* shift not more than 16 bits */
4378 *i
= (HOST_WIDE_INT
) xi
[M
] & 0xffff;
4385 if ((k
= enormlz (xi
)) > NBITS
)
4388 xi
[E
] -= (UEMUSHORT
) k
;
4394 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
4395 FRAC of e-type X. A negative input yields integer output = 0 but
4396 correct fraction. */
4399 euifrac (x
, i
, frac
)
4401 unsigned HOST_WIDE_INT
*i
;
4404 unsigned HOST_WIDE_INT ll
;
4409 k
= (int) xi
[E
] - (EXONE
- 1);
4412 /* if exponent <= 0, integer = 0 and argument is fraction */
4417 if (k
> HOST_BITS_PER_WIDE_INT
)
4419 /* Long integer overflow: output large integer
4420 and correct fraction.
4421 Note, the BSD MicroVAX compiler says that ~(0UL)
4422 is a syntax error. */
4426 warning ("overflow on truncation to unsigned integer");
4430 /* Shift more than 16 bits: first shift up k-16 mod 16,
4431 then shift up by 16's. */
4432 j
= k
- ((k
>> 4) << 4);
4439 ll
= (ll
<< 16) | xi
[M
];
4441 while ((k
-= 16) > 0);
4446 /* shift not more than 16 bits */
4448 *i
= (HOST_WIDE_INT
) xi
[M
] & 0xffff;
4451 if (xi
[0]) /* A negative value yields unsigned integer 0. */
4457 if ((k
= enormlz (xi
)) > NBITS
)
4460 xi
[E
] -= (UEMUSHORT
) k
;
4465 /* Shift the significand of exploded e-type X up or down by SC bits. */
4486 lost
|= *p
; /* remember lost bits */
4527 return ((int) lost
);
4530 /* Shift normalize the significand area of exploded e-type X.
4531 Return the shift count (up = positive). */
4537 register UEMUSHORT
*p
;
4546 return (0); /* already normalized */
4552 /* With guard word, there are NBITS+16 bits available.
4553 Return true if all are zero. */
4557 /* see if high byte is zero */
4558 while ((*p
& 0xff00) == 0)
4563 /* now shift 1 bit at a time */
4564 while ((*p
& 0x8000) == 0)
4570 mtherr ("enormlz", UNDERFLOW
);
4576 /* Normalize by shifting down out of the high guard word
4577 of the significand */
4592 mtherr ("enormlz", OVERFLOW
);
4599 /* Powers of ten used in decimal <-> binary conversions. */
4604 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
4605 static UEMUSHORT etens
[NTEN
+ 1][NE
] =
4607 {0x6576, 0x4a92, 0x804a, 0x153f,
4608 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4609 {0x6a32, 0xce52, 0x329a, 0x28ce,
4610 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4611 {0x526c, 0x50ce, 0xf18b, 0x3d28,
4612 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4613 {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4614 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4615 {0x851e, 0xeab7, 0x98fe, 0x901b,
4616 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4617 {0x0235, 0x0137, 0x36b1, 0x336c,
4618 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4619 {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4620 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4621 {0x0000, 0x0000, 0x0000, 0x0000,
4622 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4623 {0x0000, 0x0000, 0x0000, 0x0000,
4624 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4625 {0x0000, 0x0000, 0x0000, 0x0000,
4626 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4627 {0x0000, 0x0000, 0x0000, 0x0000,
4628 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4629 {0x0000, 0x0000, 0x0000, 0x0000,
4630 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4631 {0x0000, 0x0000, 0x0000, 0x0000,
4632 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4635 static UEMUSHORT emtens
[NTEN
+ 1][NE
] =
4637 {0x2030, 0xcffc, 0xa1c3, 0x8123,
4638 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4639 {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4640 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4641 {0xf53f, 0xf698, 0x6bd3, 0x0158,
4642 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4643 {0xe731, 0x04d4, 0xe3f2, 0xd332,
4644 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4645 {0xa23e, 0x5308, 0xfefb, 0x1155,
4646 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4647 {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4648 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4649 {0x2a20, 0x6224, 0x47b3, 0x98d7,
4650 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4651 {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4652 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4653 {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4654 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4655 {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4656 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4657 {0xc155, 0xa4a8, 0x404e, 0x6113,
4658 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4659 {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4660 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4661 {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4662 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4665 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4666 static UEMUSHORT etens
[NTEN
+ 1][NE
] =
4668 {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4669 {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4670 {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4671 {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4672 {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4673 {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4674 {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4675 {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4676 {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4677 {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4678 {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4679 {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4680 {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4683 static UEMUSHORT emtens
[NTEN
+ 1][NE
] =
4685 {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4686 {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4687 {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4688 {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4689 {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4690 {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4691 {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4692 {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4693 {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4694 {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4695 {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4696 {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4697 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4702 /* Convert float value X to ASCII string STRING with NDIG digits after
4703 the decimal point. */
4706 e24toasc (x
, string
, ndigs
)
4714 etoasc (w
, string
, ndigs
);
4717 /* Convert double value X to ASCII string STRING with NDIG digits after
4718 the decimal point. */
4721 e53toasc (x
, string
, ndigs
)
4729 etoasc (w
, string
, ndigs
);
4732 /* Convert double extended value X to ASCII string STRING with NDIG digits
4733 after the decimal point. */
4736 e64toasc (x
, string
, ndigs
)
4744 etoasc (w
, string
, ndigs
);
4747 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4748 after the decimal point. */
4751 e113toasc (x
, string
, ndigs
)
4759 etoasc (w
, string
, ndigs
);
4763 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4764 the decimal point. */
4766 static char wstring
[80]; /* working storage for ASCII output */
4769 etoasc (x
, string
, ndigs
)
4775 UEMUSHORT y
[NI
], t
[NI
], u
[NI
], w
[NI
];
4776 UEMUSHORT
*p
, *r
, *ten
;
4778 int i
, j
, k
, expon
, rndsav
;
4791 sprintf (wstring
, " NaN ");
4795 rndprc
= NBITS
; /* set to full precision */
4796 emov (x
, y
); /* retain external format */
4797 if (y
[NE
- 1] & 0x8000)
4800 y
[NE
- 1] &= 0x7fff;
4807 ten
= &etens
[NTEN
][0];
4809 /* Test for zero exponent */
4812 for (k
= 0; k
< NE
- 1; k
++)
4815 goto tnzro
; /* denormalized number */
4817 goto isone
; /* valid all zeros */
4821 /* Test for infinity. */
4822 if (y
[NE
- 1] == 0x7fff)
4825 sprintf (wstring
, " -Infinity ");
4827 sprintf (wstring
, " Infinity ");
4831 /* Test for exponent nonzero but significand denormalized.
4832 * This is an error condition.
4834 if ((y
[NE
- 1] != 0) && ((y
[NE
- 2] & 0x8000) == 0))
4836 mtherr ("etoasc", DOMAIN
);
4837 sprintf (wstring
, "NaN");
4841 /* Compare to 1.0 */
4850 { /* Number is greater than 1 */
4851 /* Convert significand to an integer and strip trailing decimal zeros. */
4853 u
[NE
- 1] = EXONE
+ NBITS
- 1;
4855 p
= &etens
[NTEN
- 4][0];
4861 for (j
= 0; j
< NE
- 1; j
++)
4874 /* Rescale from integer significand */
4875 u
[NE
- 1] += y
[NE
- 1] - (unsigned int) (EXONE
+ NBITS
- 1);
4877 /* Find power of 10 */
4881 /* An unordered compare result shouldn't happen here. */
4882 while (ecmp (ten
, u
) <= 0)
4884 if (ecmp (p
, u
) <= 0)
4897 { /* Number is less than 1.0 */
4898 /* Pad significand with trailing decimal zeros. */
4901 while ((y
[NE
- 2] & 0x8000) == 0)
4910 for (i
= 0; i
< NDEC
+ 1; i
++)
4912 if ((w
[NI
- 1] & 0x7) != 0)
4914 /* multiply by 10 */
4927 if (eone
[NE
- 1] <= u
[1])
4939 while (ecmp (eone
, w
) > 0)
4941 if (ecmp (p
, w
) >= 0)
4956 /* Find the first (leading) digit. */
4962 digit
= equot
[NI
- 1];
4963 while ((digit
== 0) && (ecmp (y
, ezero
) != 0))
4971 digit
= equot
[NI
- 1];
4979 /* Examine number of digits requested by caller. */
4997 *s
++ = (char)digit
+ '0';
5000 /* Generate digits after the decimal point. */
5001 for (k
= 0; k
<= ndigs
; k
++)
5003 /* multiply current number by 10, without normalizing */
5010 *s
++ = (char) equot
[NI
- 1] + '0';
5012 digit
= equot
[NI
- 1];
5015 /* round off the ASCII string */
5018 /* Test for critical rounding case in ASCII output. */
5022 if (ecmp (t
, ezero
) != 0)
5023 goto roun
; /* round to nearest */
5025 if ((*(s
- 1) & 1) == 0)
5026 goto doexp
; /* round to even */
5029 /* Round up and propagate carry-outs */
5033 /* Carry out to most significant digit? */
5040 /* Most significant digit carries to 10? */
5048 /* Round up and carry out from less significant digits */
5060 sprintf (ss, "e+%d", expon);
5062 sprintf (ss, "e%d", expon);
5064 sprintf (ss
, "e%d", expon
);
5067 /* copy out the working string */
5070 while (*ss
== ' ') /* strip possible leading space */
5072 while ((*s
++ = *ss
++) != '\0')
5077 /* Convert ASCII string to floating point.
5079 Numeric input is a free format decimal number of any length, with
5080 or without decimal point. Entering E after the number followed by an
5081 integer number causes the second number to be interpreted as a power of
5082 10 to be multiplied by the first number (i.e., "scientific" notation). */
5084 /* Convert ASCII string S to single precision float value Y. */
5095 /* Convert ASCII string S to double precision value Y. */
5102 #if defined(DEC) || defined(IBM)
5114 /* Convert ASCII string S to double extended value Y. */
5124 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
5125 /* Convert ASCII string S to 128-bit long double Y. */
5132 asctoeg (s
, y
, 113);
5136 /* Convert ASCII string S to e type Y. */
5143 asctoeg (s
, y
, NBITS
);
5146 /* Convert ASCII string SS to e type Y, with a specified rounding precision
5147 of OPREC bits. BASE is 16 for C99 hexadecimal floating constants. */
5150 asctoeg (ss
, y
, oprec
)
5155 UEMUSHORT yy
[NI
], xt
[NI
], tt
[NI
];
5156 int esign
, decflg
, sgnflg
, nexp
, exp
, prec
, lost
;
5157 int i
, k
, trail
, c
, rndsav
;
5160 char *sp
, *s
, *lstr
;
5163 /* Copy the input string. */
5164 lstr
= (char *) alloca (strlen (ss
) + 1);
5166 while (*ss
== ' ') /* skip leading spaces */
5170 while ((*sp
++ = *ss
++) != '\0')
5174 if (s
[0] == '0' && (s
[1] == 'x' || s
[1] == 'X'))
5181 rndprc
= NBITS
; /* Set to full precision */
5193 if (*s
>= '0' && *s
<= '9')
5195 else if (*s
>= 'a' && *s
<= 'f')
5199 if ((k
>= 0) && (k
< base
))
5201 /* Ignore leading zeros */
5202 if ((prec
== 0) && (decflg
== 0) && (k
== 0))
5204 /* Identify and strip trailing zeros after the decimal point. */
5205 if ((trail
== 0) && (decflg
!= 0))
5208 while ((*sp
>= '0' && *sp
<= '9')
5209 || (base
== 16 && ((*sp
>= 'a' && *sp
<= 'f')
5210 || (*sp
>= 'A' && *sp
<= 'F'))))
5212 /* Check for syntax error */
5214 if ((base
!= 10 || ((c
!= 'e') && (c
!= 'E')))
5215 && (base
!= 16 || ((c
!= 'p') && (c
!= 'P')))
5217 && (c
!= '\n') && (c
!= '\r') && (c
!= ' ')
5219 goto unexpected_char_error
;
5228 /* If enough digits were given to more than fill up the yy register,
5229 continuing until overflow into the high guard word yy[2]
5230 guarantees that there will be a roundoff bit at the top
5231 of the low guard word after normalization. */
5238 nexp
+= 4; /* count digits after decimal point */
5240 eshup1 (yy
); /* multiply current number by 16 */
5248 nexp
+= 1; /* count digits after decimal point */
5250 eshup1 (yy
); /* multiply current number by 10 */
5256 /* Insert the current digit. */
5258 xt
[NI
- 2] = (UEMUSHORT
) k
;
5263 /* Mark any lost non-zero digit. */
5265 /* Count lost digits before the decimal point. */
5287 case '.': /* decimal point */
5289 goto unexpected_char_error
;
5295 goto unexpected_char_error
;
5300 goto unexpected_char_error
;
5313 unexpected_char_error
:
5317 mtherr ("asctoe", DOMAIN
);
5326 /* Exponent interpretation */
5328 /* 0.0eXXX is zero, regardless of XXX. Check for the 0.0. */
5329 for (k
= 0; k
< NI
; k
++)
5340 /* check for + or - */
5348 while ((*s
>= '0') && (*s
<= '9'))
5357 if ((exp
> MAXDECEXP
) && (base
== 10))
5361 yy
[E
] = 0x7fff; /* infinity */
5364 if ((exp
< MINDECEXP
) && (base
== 10))
5374 /* Base 16 hexadecimal floating constant. */
5375 if ((k
= enormlz (yy
)) > NBITS
)
5380 /* Adjust the exponent. NEXP is the number of hex digits,
5381 EXP is a power of 2. */
5382 lexp
= (EXONE
- 1 + NBITS
) - k
+ yy
[E
] + exp
- nexp
;
5392 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
5393 while ((nexp
> 0) && (yy
[2] == 0))
5405 if ((k
= enormlz (yy
)) > NBITS
)
5410 lexp
= (EXONE
- 1 + NBITS
) - k
;
5411 emdnorm (yy
, lost
, 0, lexp
, 64);
5414 /* Convert to external format:
5416 Multiply by 10**nexp. If precision is 64 bits,
5417 the maximum relative error incurred in forming 10**n
5418 for 0 <= n <= 324 is 8.2e-20, at 10**180.
5419 For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5420 For 0 >= n >= -999, it is -1.55e-19 at 10**-435. */
5435 /* Punt. Can't handle this without 2 divides. */
5436 emovi (etens
[0], tt
);
5449 emul (etens
[i
], xt
, xt
);
5453 while (exp
<= MAXP
);
5472 /* Round and convert directly to the destination type */
5474 lexp
-= EXONE
- 0x3ff;
5476 else if (oprec
== 24 || oprec
== 32)
5477 lexp
-= (EXONE
- 0x7f);
5480 else if (oprec
== 24 || oprec
== 56)
5481 lexp
-= EXONE
- (0x41 << 2);
5483 else if (oprec
== 24)
5484 lexp
-= EXONE
- 0177;
5488 else if (oprec
== 56)
5489 lexp
-= EXONE
- 0201;
5492 emdnorm (yy
, lost
, 0, lexp
, 64);
5502 todec (yy
, y
); /* see etodec.c */
5507 toibm (yy
, y
, DFmode
);
5512 toc4x (yy
, y
, HFmode
);
5536 /* Return Y = largest integer not greater than X (truncated toward minus
5539 static UEMUSHORT bmask
[] =
5564 register UEMUSHORT
*p
;
5568 emov (x
, f
); /* leave in external format */
5569 expon
= (int) f
[NE
- 1];
5570 e
= (expon
& 0x7fff) - (EXONE
- 1);
5576 /* number of bits to clear out */
5588 /* clear the remaining bits */
5590 /* truncate negatives toward minus infinity */
5593 if ((UEMUSHORT
) expon
& (UEMUSHORT
) 0x8000)
5595 for (i
= 0; i
< NE
- 1; i
++)
5608 /* Return S and EXP such that S * 2^EXP = X and .5 <= S < 1.
5609 For example, 1.1 = 0.55 * 2^1. */
5621 /* Handle denormalized numbers properly using long integer exponent. */
5622 li
= (EMULONG
) ((EMUSHORT
) xi
[1]);
5630 *exp
= (int) (li
- 0x3ffe);
5634 /* Return e type Y = X * 2^PWR2. */
5650 emdnorm (xi
, i
, i
, li
, 64);
5656 /* C = remainder after dividing B by A, all e type values.
5657 Least significant integer quotient bits left in EQUOT. */
5661 UEMUSHORT a
[], b
[], c
[];
5663 UEMUSHORT den
[NI
], num
[NI
];
5667 || (ecmp (a
, ezero
) == 0)
5675 if (ecmp (a
, ezero
) == 0)
5677 mtherr ("eremain", SING
);
5683 eiremain (den
, num
);
5684 /* Sign of remainder = sign of quotient */
5693 /* Return quotient of exploded e-types NUM / DEN in EQUOT,
5694 remainder in NUM. */
5698 UEMUSHORT den
[], num
[];
5704 ld
-= enormlz (den
);
5706 ln
-= enormlz (num
);
5710 if (ecmpm (den
, num
) <= 0)
5722 emdnorm (num
, 0, 0, ln
, 0);
5725 /* Report an error condition CODE encountered in function NAME.
5727 Mnemonic Value Significance
5729 DOMAIN 1 argument domain error
5730 SING 2 function singularity
5731 OVERFLOW 3 overflow range error
5732 UNDERFLOW 4 underflow range error
5733 TLOSS 5 total loss of precision
5734 PLOSS 6 partial loss of precision
5735 INVALID 7 NaN - producing operation
5736 EDOM 33 Unix domain error code
5737 ERANGE 34 Unix range error code
5739 The order of appearance of the following messages is bound to the
5740 error codes defined above. */
5750 /* The string passed by the calling program is supposed to be the
5751 name of the function in which the error occurred.
5752 The code argument selects which error message string will be printed. */
5754 if (strcmp (name
, "esub") == 0)
5755 name
= "subtraction";
5756 else if (strcmp (name
, "ediv") == 0)
5758 else if (strcmp (name
, "emul") == 0)
5759 name
= "multiplication";
5760 else if (strcmp (name
, "enormlz") == 0)
5761 name
= "normalization";
5762 else if (strcmp (name
, "etoasc") == 0)
5763 name
= "conversion to text";
5764 else if (strcmp (name
, "asctoe") == 0)
5766 else if (strcmp (name
, "eremain") == 0)
5768 else if (strcmp (name
, "esqrt") == 0)
5769 name
= "square root";
5774 case DOMAIN
: warning ("%s: argument domain error" , name
); break;
5775 case SING
: warning ("%s: function singularity" , name
); break;
5776 case OVERFLOW
: warning ("%s: overflow range error" , name
); break;
5777 case UNDERFLOW
: warning ("%s: underflow range error" , name
); break;
5778 case TLOSS
: warning ("%s: total loss of precision" , name
); break;
5779 case PLOSS
: warning ("%s: partial loss of precision", name
); break;
5780 case INVALID
: warning ("%s: NaN - producing operation", name
); break;
5785 /* Set global error message word */
5790 /* Convert DEC double precision D to e type E. */
5798 register UEMUSHORT r
, *p
;
5800 ecleaz (y
); /* start with a zero */
5801 p
= y
; /* point to our number */
5802 r
= *d
; /* get DEC exponent word */
5803 if (*d
& (unsigned int) 0x8000)
5804 *p
= 0xffff; /* fill in our sign */
5805 ++p
; /* bump pointer to our exponent word */
5806 r
&= 0x7fff; /* strip the sign bit */
5807 if (r
== 0) /* answer = 0 if high order DEC word = 0 */
5811 r
>>= 7; /* shift exponent word down 7 bits */
5812 r
+= EXONE
- 0201; /* subtract DEC exponent offset */
5813 /* add our e type exponent offset */
5814 *p
++ = r
; /* to form our exponent */
5816 r
= *d
++; /* now do the high order mantissa */
5817 r
&= 0177; /* strip off the DEC exponent and sign bits */
5818 r
|= 0200; /* the DEC understood high order mantissa bit */
5819 *p
++ = r
; /* put result in our high guard word */
5821 *p
++ = *d
++; /* fill in the rest of our mantissa */
5825 eshdn8 (y
); /* shift our mantissa down 8 bits */
5830 /* Convert e type X to DEC double precision D. */
5841 /* Adjust exponent for offsets. */
5842 exp
= (EMULONG
) xi
[E
] - (EXONE
- 0201);
5843 /* Round off to nearest or even. */
5846 emdnorm (xi
, 0, 0, exp
, 64);
5851 /* Convert exploded e-type X, that has already been rounded to
5852 56-bit precision, to DEC format double Y. */
5898 /* Convert IBM single/double precision to e type. */
5904 enum machine_mode mode
;
5907 register UEMUSHORT r
, *p
;
5909 ecleaz (y
); /* start with a zero */
5910 p
= y
; /* point to our number */
5911 r
= *d
; /* get IBM exponent word */
5912 if (*d
& (unsigned int) 0x8000)
5913 *p
= 0xffff; /* fill in our sign */
5914 ++p
; /* bump pointer to our exponent word */
5915 r
&= 0x7f00; /* strip the sign bit */
5916 r
>>= 6; /* shift exponent word down 6 bits */
5917 /* in fact shift by 8 right and 2 left */
5918 r
+= EXONE
- (0x41 << 2); /* subtract IBM exponent offset */
5919 /* add our e type exponent offset */
5920 *p
++ = r
; /* to form our exponent */
5922 *p
++ = *d
++ & 0xff; /* now do the high order mantissa */
5923 /* strip off the IBM exponent and sign bits */
5924 if (mode
!= SFmode
) /* there are only 2 words in SFmode */
5926 *p
++ = *d
++; /* fill in the rest of our mantissa */
5931 if (y
[M
] == 0 && y
[M
+1] == 0 && y
[M
+2] == 0 && y
[M
+3] == 0)
5934 y
[E
] -= 5 + enormlz (y
); /* now normalise the mantissa */
5935 /* handle change in RADIX */
5941 /* Convert e type to IBM single/double precision. */
5946 enum machine_mode mode
;
5953 exp
= (EMULONG
) xi
[E
] - (EXONE
- (0x41 << 2)); /* adjust exponent for offsets */
5954 /* round off to nearest or even */
5957 emdnorm (xi
, 0, 0, exp
, 64);
5959 toibm (xi
, d
, mode
);
5965 enum machine_mode mode
;
6018 /* Convert C4X single/double precision to e type. */
6024 enum machine_mode mode
;
6033 /* Short-circuit the zero case. */
6034 if ((d
[0] == 0x8000)
6036 && ((mode
== QFmode
) || ((d
[2] == 0x0000) && (d
[3] == 0x0000))))
6047 ecleaz (y
); /* start with a zero */
6048 r
= d
[0]; /* get sign/exponent part */
6049 if (r
& (unsigned int) 0x0080)
6051 y
[0] = 0xffff; /* fill in our sign */
6059 r
>>= 8; /* Shift exponent word down 8 bits. */
6060 if (r
& 0x80) /* Make the exponent negative if it is. */
6062 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
6073 if (mode
!= QFmode
) /* There are only 2 words in QFmode. */
6075 y
[M
+2] = d
[2]; /* Fill in the rest of our mantissa. */
6085 /* Now do the two's complement on the data. */
6087 carry
= 1; /* Initially add 1 for the two's complement. */
6088 for (i
=size
+ M
; i
> M
; i
--)
6090 if (carry
&& (y
[i
] == 0x0000))
6092 /* We overflowed into the next word, carry is the same. */
6093 y
[i
] = carry
? 0x0000 : 0xffff;
6097 /* No overflow, just invert and add carry. */
6098 y
[i
] = ((~y
[i
]) + carry
) & 0xffff;
6113 /* Add our e type exponent offset to form our exponent. */
6117 /* Now do the high order mantissa strip off the exponent and sign
6118 bits and add the high 1 bit. */
6119 y
[M
] = (d
[0] & 0x7f) | 0x80;
6122 if (mode
!= QFmode
) /* There are only 2 words in QFmode. */
6124 y
[M
+2] = d
[2]; /* Fill in the rest of our mantissa. */
6134 /* Convert e type to C4X single/double precision. */
6139 enum machine_mode mode
;
6147 /* Adjust exponent for offsets. */
6148 exp
= (EMULONG
) xi
[E
] - (EXONE
- 0x7f);
6150 /* Round off to nearest or even. */
6152 rndprc
= mode
== QFmode
? 24 : 32;
6153 emdnorm (xi
, 0, 0, exp
, 64);
6155 toc4x (xi
, d
, mode
);
6161 enum machine_mode mode
;
6167 /* Short-circuit the zero case */
6168 if ((x
[0] == 0) /* Zero exponent and sign */
6170 && (x
[M
] == 0) /* The rest is for zero mantissa */
6172 /* Only check for double if necessary */
6173 && ((mode
== QFmode
) || ((x
[M
+2] == 0) && (x
[M
+3] == 0))))
6175 /* We have a zero. Put it into the output and return. */
6188 /* Negative number require a two's complement conversion of the
6194 i
= ((int) x
[1]) - 0x7f;
6196 /* Now add 1 to the inverted data to do the two's complement. */
6206 x
[v
] = carry
? 0x0000 : 0xffff;
6210 x
[v
] = ((~x
[v
]) + carry
) & 0xffff;
6216 /* The following is a special case. The C4X negative float requires
6217 a zero in the high bit (because the format is (2 - x) x 2^m), so
6218 if a one is in that bit, we have to shift left one to get rid
6219 of it. This only occurs if the number is -1 x 2^m. */
6220 if (x
[M
+1] & 0x8000)
6222 /* This is the case of -1 x 2^m, we have to rid ourselves of the
6223 high sign bit and shift the exponent. */
6230 i
= ((int) x
[1]) - 0x7f;
6233 if ((i
< -128) || (i
> 127))
6248 y
[0] |= ((i
& 0xff) << 8);
6252 y
[0] |= x
[M
] & 0x7f;
6262 /* Output a binary NaN bit pattern in the target machine's format. */
6264 /* If special NaN bit patterns are required, define them in tm.h
6265 as arrays of unsigned 16-bit shorts. Otherwise, use the default
6271 UEMUSHORT TFbignan
[8] =
6272 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6273 UEMUSHORT TFlittlenan
[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
6281 UEMUSHORT XFbignan
[6] =
6282 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
6283 UEMUSHORT XFlittlenan
[6] = {0, 0, 0, 0xc000, 0xffff, 0};
6291 UEMUSHORT DFbignan
[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
6292 UEMUSHORT DFlittlenan
[4] = {0, 0, 0, 0xfff8};
6300 UEMUSHORT SFbignan
[2] = {0x7fff, 0xffff};
6301 UEMUSHORT SFlittlenan
[2] = {0, 0xffc0};
6308 make_nan (nan
, sign
, mode
)
6311 enum machine_mode mode
;
6318 /* Possibly the `reserved operand' patterns on a VAX can be
6319 used like NaN's, but probably not in the same way as IEEE. */
6320 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
6322 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6324 if (REAL_WORDS_BIG_ENDIAN
)
6334 if (REAL_WORDS_BIG_ENDIAN
)
6342 if (REAL_WORDS_BIG_ENDIAN
)
6351 if (REAL_WORDS_BIG_ENDIAN
)
6361 if (REAL_WORDS_BIG_ENDIAN
)
6362 *nan
++ = (sign
<< 15) | (*p
++ & 0x7fff);
6365 if (! REAL_WORDS_BIG_ENDIAN
)
6366 *nan
= (sign
<< 15) | (*p
& 0x7fff);
6370 /* This is the inverse of the function `etarsingle' invoked by
6371 REAL_VALUE_TO_TARGET_SINGLE. */
6374 ereal_unto_float (f
)
6381 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6382 This is the inverse operation to what the function `endian' does. */
6383 if (REAL_WORDS_BIG_ENDIAN
)
6385 s
[0] = (UEMUSHORT
) (f
>> 16);
6386 s
[1] = (UEMUSHORT
) f
;
6390 s
[0] = (UEMUSHORT
) f
;
6391 s
[1] = (UEMUSHORT
) (f
>> 16);
6393 /* Convert and promote the target float to E-type. */
6395 /* Output E-type to REAL_VALUE_TYPE. */
6401 /* This is the inverse of the function `etardouble' invoked by
6402 REAL_VALUE_TO_TARGET_DOUBLE. */
6405 ereal_unto_double (d
)
6412 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
6413 if (REAL_WORDS_BIG_ENDIAN
)
6415 s
[0] = (UEMUSHORT
) (d
[0] >> 16);
6416 s
[1] = (UEMUSHORT
) d
[0];
6417 s
[2] = (UEMUSHORT
) (d
[1] >> 16);
6418 s
[3] = (UEMUSHORT
) d
[1];
6422 /* Target float words are little-endian. */
6423 s
[0] = (UEMUSHORT
) d
[0];
6424 s
[1] = (UEMUSHORT
) (d
[0] >> 16);
6425 s
[2] = (UEMUSHORT
) d
[1];
6426 s
[3] = (UEMUSHORT
) (d
[1] >> 16);
6428 /* Convert target double to E-type. */
6430 /* Output E-type to REAL_VALUE_TYPE. */
6436 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
6437 This is somewhat like ereal_unto_float, but the input types
6438 for these are different. */
6441 ereal_from_float (f
)
6448 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6449 This is the inverse operation to what the function `endian' does. */
6450 if (REAL_WORDS_BIG_ENDIAN
)
6452 s
[0] = (UEMUSHORT
) (f
>> 16);
6453 s
[1] = (UEMUSHORT
) f
;
6457 s
[0] = (UEMUSHORT
) f
;
6458 s
[1] = (UEMUSHORT
) (f
>> 16);
6460 /* Convert and promote the target float to E-type. */
6462 /* Output E-type to REAL_VALUE_TYPE. */
6468 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
6469 This is somewhat like ereal_unto_double, but the input types
6470 for these are different.
6472 The DFmode is stored as an array of HOST_WIDE_INT in the target's
6473 data format, with no holes in the bit packing. The first element
6474 of the input array holds the bits that would come first in the
6475 target computer's memory. */
6478 ereal_from_double (d
)
6485 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
6486 if (REAL_WORDS_BIG_ENDIAN
)
6488 #if HOST_BITS_PER_WIDE_INT == 32
6489 s
[0] = (UEMUSHORT
) (d
[0] >> 16);
6490 s
[1] = (UEMUSHORT
) d
[0];
6491 s
[2] = (UEMUSHORT
) (d
[1] >> 16);
6492 s
[3] = (UEMUSHORT
) d
[1];
6494 /* In this case the entire target double is contained in the
6495 first array element. The second element of the input is
6497 s
[0] = (UEMUSHORT
) (d
[0] >> 48);
6498 s
[1] = (UEMUSHORT
) (d
[0] >> 32);
6499 s
[2] = (UEMUSHORT
) (d
[0] >> 16);
6500 s
[3] = (UEMUSHORT
) d
[0];
6505 /* Target float words are little-endian. */
6506 s
[0] = (UEMUSHORT
) d
[0];
6507 s
[1] = (UEMUSHORT
) (d
[0] >> 16);
6508 #if HOST_BITS_PER_WIDE_INT == 32
6509 s
[2] = (UEMUSHORT
) d
[1];
6510 s
[3] = (UEMUSHORT
) (d
[1] >> 16);
6512 s
[2] = (UEMUSHORT
) (d
[0] >> 32);
6513 s
[3] = (UEMUSHORT
) (d
[0] >> 48);
6516 /* Convert target double to E-type. */
6518 /* Output E-type to REAL_VALUE_TYPE. */
6525 /* Convert target computer unsigned 64-bit integer to e-type.
6526 The endian-ness of DImode follows the convention for integers,
6527 so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN. */
6531 UEMUSHORT
*di
; /* Address of the 64-bit int. */
6538 if (WORDS_BIG_ENDIAN
)
6540 for (k
= M
; k
< M
+ 4; k
++)
6545 for (k
= M
+ 3; k
>= M
; k
--)
6548 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
6549 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
6550 ecleaz (yi
); /* it was zero */
6552 yi
[E
] -= (UEMUSHORT
) k
;/* subtract shift count from exponent */
6556 /* Convert target computer signed 64-bit integer to e-type. */
6560 UEMUSHORT
*di
; /* Address of the 64-bit int. */
6563 unsigned EMULONG acc
;
6569 if (WORDS_BIG_ENDIAN
)
6571 for (k
= M
; k
< M
+ 4; k
++)
6576 for (k
= M
+ 3; k
>= M
; k
--)
6579 /* Take absolute value */
6585 for (k
= M
+ 3; k
>= M
; k
--)
6587 acc
= (unsigned EMULONG
) (~yi
[k
] & 0xffff) + carry
;
6594 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
6595 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
6596 ecleaz (yi
); /* it was zero */
6598 yi
[E
] -= (UEMUSHORT
) k
;/* subtract shift count from exponent */
6605 /* Convert e-type to unsigned 64-bit int. */
6621 k
= (int) xi
[E
] - (EXONE
- 1);
6624 for (j
= 0; j
< 4; j
++)
6630 for (j
= 0; j
< 4; j
++)
6633 warning ("overflow on truncation to integer");
6638 /* Shift more than 16 bits: first shift up k-16 mod 16,
6639 then shift up by 16's. */
6640 j
= k
- ((k
>> 4) << 4);
6644 if (WORDS_BIG_ENDIAN
)
6655 if (WORDS_BIG_ENDIAN
)
6660 while ((k
-= 16) > 0);
6664 /* shift not more than 16 bits */
6669 if (WORDS_BIG_ENDIAN
)
6688 /* Convert e-type to signed 64-bit int. */
6695 unsigned EMULONG acc
;
6702 k
= (int) xi
[E
] - (EXONE
- 1);
6705 for (j
= 0; j
< 4; j
++)
6711 for (j
= 0; j
< 4; j
++)
6714 warning ("overflow on truncation to integer");
6720 /* Shift more than 16 bits: first shift up k-16 mod 16,
6721 then shift up by 16's. */
6722 j
= k
- ((k
>> 4) << 4);
6726 if (WORDS_BIG_ENDIAN
)
6737 if (WORDS_BIG_ENDIAN
)
6742 while ((k
-= 16) > 0);
6746 /* shift not more than 16 bits */
6749 if (WORDS_BIG_ENDIAN
)
6765 /* Negate if negative */
6769 if (WORDS_BIG_ENDIAN
)
6771 for (k
= 0; k
< 4; k
++)
6773 acc
= (unsigned EMULONG
) (~(*isave
) & 0xffff) + carry
;
6774 if (WORDS_BIG_ENDIAN
)
6786 /* Longhand square root routine. */
6789 static int esqinited
= 0;
6790 static unsigned short sqrndbit
[NI
];
6796 UEMUSHORT temp
[NI
], num
[NI
], sq
[NI
], xx
[NI
];
6798 int i
, j
, k
, n
, nlups
;
6803 sqrndbit
[NI
- 2] = 1;
6806 /* Check for arg <= 0 */
6807 i
= ecmp (x
, ezero
);
6812 mtherr ("esqrt", DOMAIN
);
6828 /* Bring in the arg and renormalize if it is denormal. */
6830 m
= (EMULONG
) xx
[1]; /* local long word exponent */
6834 /* Divide exponent by 2 */
6836 exp
= (unsigned short) ((m
/ 2) + 0x3ffe);
6838 /* Adjust if exponent odd */
6848 n
= 8; /* get 8 bits of result per inner loop */
6854 /* bring in next word of arg */
6856 num
[NI
- 1] = xx
[j
+ 3];
6857 /* Do additional bit on last outer loop, for roundoff. */
6860 for (i
= 0; i
< n
; i
++)
6862 /* Next 2 bits of arg */
6865 /* Shift up answer */
6867 /* Make trial divisor */
6868 for (k
= 0; k
< NI
; k
++)
6871 eaddm (sqrndbit
, temp
);
6872 /* Subtract and insert answer bit if it goes in */
6873 if (ecmpm (temp
, num
) <= 0)
6883 /* Adjust for extra, roundoff loop done. */
6884 exp
+= (NBITS
- 1) - rndprc
;
6886 /* Sticky bit = 1 if the remainder is nonzero. */
6888 for (i
= 3; i
< NI
; i
++)
6891 /* Renormalize and round off. */
6892 emdnorm (sq
, k
, 0, exp
, 64);
6896 #endif /* EMU_NON_COMPILE not defined */
6898 /* Return the binary precision of the significand for a given
6899 floating point mode. The mode can hold an integer value
6900 that many bits wide, without losing any bits. */
6903 significand_size (mode
)
6904 enum machine_mode mode
;
6907 /* Don't test the modes, but their sizes, lest this
6908 code won't work for BITS_PER_UNIT != 8 . */
6910 switch (GET_MODE_BITSIZE (mode
))
6914 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
6921 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
6924 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
6927 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
6930 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
6943 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)