1 /* real.c - implementation of REAL_ARITHMETIC, REAL_VALUE_ATOF,
2 and support for XFmode IEEE extended real floating point arithmetic.
3 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998,
4 1999, 2000, 2002 Free Software Foundation, Inc.
5 Contributed by Stephen L. Moshier (moshier@world.std.com).
7 This file is part of GCC.
9 GCC is free software; you can redistribute it and/or modify it under
10 the terms of the GNU General Public License as published by the Free
11 Software Foundation; either version 2, or (at your option) any later
14 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
15 WARRANTY; without even the implied warranty of MERCHANTABILITY or
16 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19 You should have received a copy of the GNU General Public License
20 along with GCC; see the file COPYING. If not, write to the Free
21 Software Foundation, 59 Temple Place - Suite 330, Boston, MA
31 /* To enable support of XFmode extended real floating point, define
32 LONG_DOUBLE_TYPE_SIZE 96 in the tm.h file (m68k.h or i386.h).
34 Machine files (tm.h etc) must not contain any code
35 that tries to use host floating point arithmetic to convert
36 REAL_VALUE_TYPEs from `double' to `float', pass them to fprintf,
37 etc. In cross-compile situations a REAL_VALUE_TYPE may not
38 be intelligible to the host computer's native arithmetic.
40 The first part of this file interfaces gcc to a floating point
41 arithmetic suite that was not written with gcc in mind. Avoid
42 changing the low-level arithmetic routines unless you have suitable
43 test programs available. A special version of the PARANOIA floating
44 point arithmetic tester, modified for this purpose, can be found on
45 usc.edu: /pub/C-numanal/ieeetest.zoo. Other tests, and libraries of
46 XFmode and TFmode transcendental functions, can be obtained by ftp from
47 netlib.att.com: netlib/cephes. */
49 /* Type of computer arithmetic.
50 Only one of DEC, IBM, IEEE, C4X, or UNK should get defined.
52 `IEEE', when REAL_WORDS_BIG_ENDIAN is non-zero, refers generically
53 to big-endian IEEE floating-point data structure. This definition
54 should work in SFmode `float' type and DFmode `double' type on
55 virtually all big-endian IEEE machines. If LONG_DOUBLE_TYPE_SIZE
56 has been defined to be 96, then IEEE also invokes the particular
57 XFmode (`long double' type) data structure used by the Motorola
58 680x0 series processors.
60 `IEEE', when REAL_WORDS_BIG_ENDIAN is zero, refers generally to
61 little-endian IEEE machines. In this case, if LONG_DOUBLE_TYPE_SIZE
62 has been defined to be 96, then IEEE also invokes the particular
63 XFmode `long double' data structure used by the Intel 80x86 series
66 `DEC' refers specifically to the Digital Equipment Corp PDP-11
67 and VAX floating point data structure. This model currently
68 supports no type wider than DFmode.
70 `IBM' refers specifically to the IBM System/370 and compatible
71 floating point data structure. This model currently supports
72 no type wider than DFmode. The IBM conversions were contributed by
73 frank@atom.ansto.gov.au (Frank Crawford).
75 `C4X' refers specifically to the floating point format used on
76 Texas Instruments TMS320C3x and TMS320C4x digital signal
77 processors. This supports QFmode (32-bit float, double) and HFmode
78 (40-bit long double) where BITS_PER_BYTE is 32. Unlike IEEE
79 floats, C4x floats are not rounded to be even. The C4x conversions
80 were contributed by m.hayes@elec.canterbury.ac.nz (Michael Hayes) and
81 Haj.Ten.Brugge@net.HCC.nl (Herman ten Brugge).
83 If LONG_DOUBLE_TYPE_SIZE = 64 (the default, unless tm.h defines it)
84 then `long double' and `double' are both implemented, but they
87 The case LONG_DOUBLE_TYPE_SIZE = 128 activates TFmode support
88 and may deactivate XFmode since `long double' is used to refer
89 to both modes. Defining INTEL_EXTENDED_IEEE_FORMAT to non-zero
90 at the same time enables 80387-style 80-bit floats in a 128-bit
91 padded image, as seen on IA-64.
93 The macros FLOAT_WORDS_BIG_ENDIAN, HOST_FLOAT_WORDS_BIG_ENDIAN,
94 contributed by Richard Earnshaw <Richard.Earnshaw@cl.cam.ac.uk>,
95 separate the floating point unit's endian-ness from that of
96 the integer addressing. This permits one to define a big-endian
97 FPU on a little-endian machine (e.g., ARM). An extension to
98 BYTES_BIG_ENDIAN may be required for some machines in the future.
99 These optional macros may be defined in tm.h. In real.h, they
100 default to WORDS_BIG_ENDIAN, etc., so there is no need to define
101 them for any normal host or target machine on which the floats
102 and the integers have the same endian-ness. */
105 /* The following converts gcc macros into the ones used by this file. */
107 #if TARGET_FLOAT_FORMAT == VAX_FLOAT_FORMAT
108 /* PDP-11, Pro350, VAX: */
110 #else /* it's not VAX */
111 #if TARGET_FLOAT_FORMAT == IBM_FLOAT_FORMAT
112 /* IBM System/370 style */
114 #else /* it's also not an IBM */
115 #if TARGET_FLOAT_FORMAT == C4X_FLOAT_FORMAT
116 /* TMS320C3x/C4x style */
118 #else /* it's also not a C4X */
119 #if TARGET_FLOAT_FORMAT == IEEE_FLOAT_FORMAT
121 #else /* it's not IEEE either */
122 /* UNKnown arithmetic. We don't support this and can't go on. */
123 unknown arithmetic type
125 #endif /* not IEEE */
130 #define REAL_WORDS_BIG_ENDIAN FLOAT_WORDS_BIG_ENDIAN
132 /* Make sure that the endianness is correct for IBM and DEC. */
134 #undef LARGEST_EXPONENT_IS_NORMAL
135 #define LARGEST_EXPONENT_IS_NORMAL(x) 1
136 #undef REAL_WORDS_BIG_ENDIAN
137 /* Strangely enough, DEC float most closely resembles big endian IEEE */
138 #define REAL_WORDS_BIG_ENDIAN 1
139 /* ... but the halfwords are reversed from IEEE big endian. */
140 #ifndef VAX_HALFWORD_ORDER
141 #define VAX_HALFWORD_ORDER 1
144 #if defined(IBM) && !REAL_WORDS_BIG_ENDIAN
145 #error "Little-endian representations are not supported for IBM."
149 #if defined(DEC) && !defined (TARGET_G_FLOAT)
150 #define TARGET_G_FLOAT 0
153 #ifndef VAX_HALFWORD_ORDER
154 #define VAX_HALFWORD_ORDER 0
157 /* Define INFINITY for support of infinity.
158 Define NANS for support of Not-a-Number's (NaN's). */
159 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
164 /* Support of NaNs requires support of infinity. */
171 /* Find a host integer type that is at least 16 bits wide,
172 and another type at least twice whatever that size is. */
174 #if HOST_BITS_PER_CHAR >= 16
175 #define EMUSHORT char
176 #define EMUSHORT_SIZE HOST_BITS_PER_CHAR
177 #define EMULONG_SIZE (2 * HOST_BITS_PER_CHAR)
179 #if HOST_BITS_PER_SHORT >= 16
180 #define EMUSHORT short
181 #define EMUSHORT_SIZE HOST_BITS_PER_SHORT
182 #define EMULONG_SIZE (2 * HOST_BITS_PER_SHORT)
184 #if HOST_BITS_PER_INT >= 16
186 #define EMUSHORT_SIZE HOST_BITS_PER_INT
187 #define EMULONG_SIZE (2 * HOST_BITS_PER_INT)
189 #if HOST_BITS_PER_LONG >= 16
190 #define EMUSHORT long
191 #define EMUSHORT_SIZE HOST_BITS_PER_LONG
192 #define EMULONG_SIZE (2 * HOST_BITS_PER_LONG)
194 #error "You will have to modify this program to have a smaller unit size."
200 /* If no 16-bit type has been found and the compiler is GCC, try HImode. */
201 #if defined(__GNUC__) && EMUSHORT_SIZE != 16
202 typedef int HItype
__attribute__ ((mode (HI
)));
203 typedef unsigned int UHItype
__attribute__ ((mode (HI
)));
207 #define EMUSHORT HItype
208 #define UEMUSHORT UHItype
209 #define EMUSHORT_SIZE 16
210 #define EMULONG_SIZE 32
212 #define UEMUSHORT unsigned EMUSHORT
215 #if HOST_BITS_PER_SHORT >= EMULONG_SIZE
216 #define EMULONG short
218 #if HOST_BITS_PER_INT >= EMULONG_SIZE
221 #if HOST_BITS_PER_LONG >= EMULONG_SIZE
224 #if HOST_BITS_PER_LONGLONG >= EMULONG_SIZE
225 #define EMULONG long long int
227 #error "You will have to modify this program to have a smaller unit size."
233 #if EMUSHORT_SIZE != 16
234 #error "The host interface doesn't work if no 16-bit size exists."
237 /* Calculate the size of the generic "e" type. This always has
238 identical in-memory size to REAL_VALUE_TYPE. The sizes are supposed
239 to be the same as well, but when REAL_VALUE_TYPE_SIZE is not evenly
240 divisible by HOST_BITS_PER_WIDE_INT we have some padding in
242 There are only two supported sizes: ten and six 16-bit words (160
245 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && !INTEL_EXTENDED_IEEE_FORMAT
248 # define MAXDECEXP 4932
249 # define MINDECEXP -4977
252 # define MAXDECEXP 4932
253 # define MINDECEXP -4956
256 /* Fail compilation if 2*NE is not the appropriate size.
257 If HOST_BITS_PER_WIDE_INT is 64, we're going to have padding
258 at the end of the array, because neither 96 nor 160 is
259 evenly divisible by 64. */
260 struct compile_test_dummy
{
261 char twice_NE_must_equal_sizeof_REAL_VALUE_TYPE
262 [(sizeof (REAL_VALUE_TYPE
) >= 2*NE
) ? 1 : -1];
265 /* Construct macros to translate between REAL_VALUE_TYPE and e type.
266 In GET_REAL and PUT_REAL, r and e are pointers.
267 A REAL_VALUE_TYPE is guaranteed to occupy contiguous locations
268 in memory, with no holes. */
269 #define GET_REAL(r, e) memcpy ((e), (r), 2*NE)
270 #define PUT_REAL(e, r) \
272 memcpy (r, e, 2*NE); \
273 if (2*NE < sizeof (*r)) \
274 memset ((char *) (r) + 2*NE, 0, sizeof (*r) - 2*NE); \
277 /* Number of 16 bit words in internal format */
280 /* Array offset to exponent */
283 /* Array offset to high guard word */
286 /* Number of bits of precision */
287 #define NBITS ((NI-4)*16)
289 /* Maximum number of decimal digits in ASCII conversion
292 #define NDEC (NBITS*8/27)
294 /* The exponent of 1.0 */
295 #define EXONE (0x3fff)
297 #if defined(HOST_EBCDIC)
298 /* bit 8 is significant in EBCDIC */
299 #define CHARMASK 0xff
301 #define CHARMASK 0x7f
304 /* Information about the various IEEE precisions. At the moment, we only
305 support exponents of 15 bits or less. */
311 /* Size of the exponent in bits. */
314 /* Overall size of the value in bits. */
317 /* Mode used for representing the value. */
318 enum machine_mode mode
;
320 /* Exponent adjustment for offsets. */
325 /* IEEE float (24 bits). */
326 static const struct ieee_format ieee_24
=
335 /* IEEE double (53 bits). */
336 static const struct ieee_format ieee_53
=
347 /* IEEE extended double (64 bits). */
348 static const struct ieee_format ieee_64
=
357 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
358 /* IEEE long double (113 bits). */
359 static const struct ieee_format ieee_113
=
367 #endif /* INTEL_EXTENDED_IEEE_FORMAT == 0 */
370 /* DEC F float (24 bits). */
371 static const struct ieee_format dec_f
=
380 /* DEC D float (56 bits). */
381 static const struct ieee_format dec_d
=
390 /* DEC G float (53 bits). */
391 static const struct ieee_format dec_g
=
401 /* DEC H float (113 bits). (not yet used) */
402 static const struct ieee_format dec_h
=
413 extern int extra_warnings
;
414 extern const UEMUSHORT ezero
[NE
], ehalf
[NE
], eone
[NE
], etwo
[NE
];
415 extern const UEMUSHORT elog2
[NE
], esqrt2
[NE
];
417 static void endian
PARAMS ((const UEMUSHORT
*, long *,
419 static void eclear
PARAMS ((UEMUSHORT
*));
420 static void emov
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
422 static void eabs
PARAMS ((UEMUSHORT
*));
424 static void eneg
PARAMS ((UEMUSHORT
*));
425 static int eisneg
PARAMS ((const UEMUSHORT
*));
426 static int eisinf
PARAMS ((const UEMUSHORT
*));
427 static int eisnan
PARAMS ((const UEMUSHORT
*));
428 static void einfin
PARAMS ((UEMUSHORT
*));
430 static void enan
PARAMS ((UEMUSHORT
*, int));
431 static void einan
PARAMS ((UEMUSHORT
*));
432 static int eiisnan
PARAMS ((const UEMUSHORT
*));
433 static void make_nan
PARAMS ((UEMUSHORT
*, int, enum machine_mode
));
435 static int eiisneg
PARAMS ((const UEMUSHORT
*));
436 static void saturate
PARAMS ((UEMUSHORT
*, int, int, int));
437 static void emovi
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
438 static void emovo
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
439 static void ecleaz
PARAMS ((UEMUSHORT
*));
440 static void ecleazs
PARAMS ((UEMUSHORT
*));
441 static void emovz
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
443 static void eiinfin
PARAMS ((UEMUSHORT
*));
446 static int eiisinf
PARAMS ((const UEMUSHORT
*));
448 static int ecmpm
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*));
449 static void eshdn1
PARAMS ((UEMUSHORT
*));
450 static void eshup1
PARAMS ((UEMUSHORT
*));
451 static void eshdn8
PARAMS ((UEMUSHORT
*));
452 static void eshup8
PARAMS ((UEMUSHORT
*));
453 static void eshup6
PARAMS ((UEMUSHORT
*));
454 static void eshdn6
PARAMS ((UEMUSHORT
*));
455 static void eaddm
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));\f
456 static void esubm
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
457 static void m16m
PARAMS ((unsigned int, const UEMUSHORT
*, UEMUSHORT
*));
458 static int edivm
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
459 static int emulm
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
460 static void emdnorm
PARAMS ((UEMUSHORT
*, int, int, EMULONG
, int));
461 static void esub
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
463 static void eadd
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
465 static void eadd1
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
467 static void ediv
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
469 static void emul
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
471 static void e53toe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
472 static void e64toe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
473 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
474 static void e113toe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
476 static void e24toe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
477 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
478 static void etoe113
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
479 static void toe113
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
481 static void etoe64
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
482 static void toe64
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
483 static void etoe53
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
484 static void toe53
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
485 static void etoe24
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
486 static void toe24
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
487 static void ieeetoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
488 const struct ieee_format
*));
489 static void etoieee
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
490 const struct ieee_format
*));
491 static void toieee
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
492 const struct ieee_format
*));
493 static int ecmp
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*));
495 static void eround
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
497 static void ltoe
PARAMS ((const HOST_WIDE_INT
*, UEMUSHORT
*));
498 static void ultoe
PARAMS ((const unsigned HOST_WIDE_INT
*, UEMUSHORT
*));
499 static void eifrac
PARAMS ((const UEMUSHORT
*, HOST_WIDE_INT
*,
501 static void euifrac
PARAMS ((const UEMUSHORT
*, unsigned HOST_WIDE_INT
*,
503 static int eshift
PARAMS ((UEMUSHORT
*, int));
504 static int enormlz
PARAMS ((UEMUSHORT
*));
506 static void e24toasc
PARAMS ((const UEMUSHORT
*, char *, int));
507 static void e53toasc
PARAMS ((const UEMUSHORT
*, char *, int));
508 static void e64toasc
PARAMS ((const UEMUSHORT
*, char *, int));
509 static void e113toasc
PARAMS ((const UEMUSHORT
*, char *, int));
511 static void etoasc
PARAMS ((const UEMUSHORT
*, char *, int));
512 static void asctoe24
PARAMS ((const char *, UEMUSHORT
*));
513 static void asctoe53
PARAMS ((const char *, UEMUSHORT
*));
514 static void asctoe64
PARAMS ((const char *, UEMUSHORT
*));
515 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
516 static void asctoe113
PARAMS ((const char *, UEMUSHORT
*));
518 static void asctoe
PARAMS ((const char *, UEMUSHORT
*));
519 static void asctoeg
PARAMS ((const char *, UEMUSHORT
*, int));
520 static void efloor
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
522 static void efrexp
PARAMS ((const UEMUSHORT
*, int *,
525 static void eldexp
PARAMS ((const UEMUSHORT
*, int, UEMUSHORT
*));
527 static void eremain
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
530 static void eiremain
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
531 static void mtherr
PARAMS ((const char *, int));
533 static void dectoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
534 static void etodec
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
535 static void todec
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
538 static void ibmtoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
540 static void etoibm
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
542 static void toibm
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
546 static void c4xtoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
548 static void etoc4x
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
550 static void toc4x
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
554 static void uditoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
555 static void ditoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
556 static void etoudi
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
557 static void etodi
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
558 static void esqrt
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
561 /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
562 swapping ends if required, into output array of longs. The
563 result is normally passed to fprintf by the ASM_OUTPUT_ macros. */
569 enum machine_mode mode
;
573 if (REAL_WORDS_BIG_ENDIAN
&& !VAX_HALFWORD_ORDER
)
578 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
579 /* Swap halfwords in the fourth long. */
580 th
= (unsigned long) e
[6] & 0xffff;
581 t
= (unsigned long) e
[7] & 0xffff;
590 /* Swap halfwords in the third long. */
591 th
= (unsigned long) e
[4] & 0xffff;
592 t
= (unsigned long) e
[5] & 0xffff;
598 /* Swap halfwords in the second word. */
599 th
= (unsigned long) e
[2] & 0xffff;
600 t
= (unsigned long) e
[3] & 0xffff;
607 /* Swap halfwords in the first word. */
608 th
= (unsigned long) e
[0] & 0xffff;
609 t
= (unsigned long) e
[1] & 0xffff;
620 /* Pack the output array without swapping. */
625 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
626 /* Pack the fourth long. */
627 th
= (unsigned long) e
[7] & 0xffff;
628 t
= (unsigned long) e
[6] & 0xffff;
637 /* Pack the third long.
638 Each element of the input REAL_VALUE_TYPE array has 16 useful bits
640 th
= (unsigned long) e
[5] & 0xffff;
641 t
= (unsigned long) e
[4] & 0xffff;
647 /* Pack the second long */
648 th
= (unsigned long) e
[3] & 0xffff;
649 t
= (unsigned long) e
[2] & 0xffff;
656 /* Pack the first long */
657 th
= (unsigned long) e
[1] & 0xffff;
658 t
= (unsigned long) e
[0] & 0xffff;
670 /* This is the implementation of the REAL_ARITHMETIC macro. */
673 earith (value
, icode
, r1
, r2
)
674 REAL_VALUE_TYPE
*value
;
679 UEMUSHORT d1
[NE
], d2
[NE
], v
[NE
];
685 /* Return NaN input back to the caller. */
688 PUT_REAL (d1
, value
);
693 PUT_REAL (d2
, value
);
697 code
= (enum tree_code
) icode
;
705 esub (d2
, d1
, v
); /* d1 - d2 */
714 if (ecmp (d2
, ezero
) == 0)
717 ediv (d2
, d1
, v
); /* d1/d2 */
720 case MIN_EXPR
: /* min (d1,d2) */
721 if (ecmp (d1
, d2
) < 0)
727 case MAX_EXPR
: /* max (d1,d2) */
728 if (ecmp (d1
, d2
) > 0)
741 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
742 implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */
748 UEMUSHORT f
[NE
], g
[NE
];
764 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
765 implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */
771 UEMUSHORT f
[NE
], g
[NE
];
773 unsigned HOST_WIDE_INT l
;
787 /* This is the REAL_VALUE_ATOF function. It converts a decimal or hexadecimal
788 string to binary, rounding off as indicated by the machine_mode argument.
789 Then it promotes the rounded value to REAL_VALUE_TYPE. */
796 UEMUSHORT tem
[NE
], e
[NE
];
822 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
842 /* Expansion of REAL_NEGATE. */
858 /* Round real toward zero to HOST_WIDE_INT;
859 implements REAL_VALUE_FIX (x). */
865 UEMUSHORT f
[NE
], g
[NE
];
872 warning ("conversion from NaN to int");
880 /* Round real toward zero to unsigned HOST_WIDE_INT
881 implements REAL_VALUE_UNSIGNED_FIX (x).
882 Negative input returns zero. */
884 unsigned HOST_WIDE_INT
888 UEMUSHORT f
[NE
], g
[NE
];
889 unsigned HOST_WIDE_INT l
;
895 warning ("conversion from NaN to unsigned int");
904 /* REAL_VALUE_FROM_INT macro. */
907 ereal_from_int (d
, i
, j
, mode
)
910 enum machine_mode mode
;
912 UEMUSHORT df
[NE
], dg
[NE
];
913 HOST_WIDE_INT low
, high
;
916 if (GET_MODE_CLASS (mode
) != MODE_FLOAT
)
923 /* complement and add 1 */
930 eldexp (eone
, HOST_BITS_PER_WIDE_INT
, df
);
931 ultoe ((unsigned HOST_WIDE_INT
*) &high
, dg
);
933 ultoe ((unsigned HOST_WIDE_INT
*) &low
, df
);
938 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
939 Avoid double-rounding errors later by rounding off now from the
940 extra-wide internal format to the requested precision. */
941 switch (GET_MODE_BITSIZE (mode
))
959 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
976 /* REAL_VALUE_FROM_UNSIGNED_INT macro. */
979 ereal_from_uint (d
, i
, j
, mode
)
981 unsigned HOST_WIDE_INT i
, j
;
982 enum machine_mode mode
;
984 UEMUSHORT df
[NE
], dg
[NE
];
985 unsigned HOST_WIDE_INT low
, high
;
987 if (GET_MODE_CLASS (mode
) != MODE_FLOAT
)
991 eldexp (eone
, HOST_BITS_PER_WIDE_INT
, df
);
997 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
998 Avoid double-rounding errors later by rounding off now from the
999 extra-wide internal format to the requested precision. */
1000 switch (GET_MODE_BITSIZE (mode
))
1018 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
1035 /* REAL_VALUE_TO_INT macro. */
1038 ereal_to_int (low
, high
, rr
)
1039 HOST_WIDE_INT
*low
, *high
;
1042 UEMUSHORT d
[NE
], df
[NE
], dg
[NE
], dh
[NE
];
1049 warning ("conversion from NaN to int");
1055 /* convert positive value */
1062 eldexp (eone
, HOST_BITS_PER_WIDE_INT
, df
);
1063 ediv (df
, d
, dg
); /* dg = d / 2^32 is the high word */
1064 euifrac (dg
, (unsigned HOST_WIDE_INT
*) high
, dh
);
1065 emul (df
, dh
, dg
); /* fractional part is the low word */
1066 euifrac (dg
, (unsigned HOST_WIDE_INT
*) low
, dh
);
1069 /* complement and add 1 */
1079 /* REAL_VALUE_LDEXP macro. */
1086 UEMUSHORT e
[NE
], y
[NE
];
1099 /* Check for infinity in a REAL_VALUE_TYPE. */
1103 REAL_VALUE_TYPE x ATTRIBUTE_UNUSED
;
1109 return (eisinf (e
));
1115 /* Check whether a REAL_VALUE_TYPE item is a NaN. */
1119 REAL_VALUE_TYPE x ATTRIBUTE_UNUSED
;
1125 return (eisnan (e
));
1132 /* Check for a negative REAL_VALUE_TYPE number.
1133 This just checks the sign bit, so that -0 counts as negative. */
1139 return ereal_isneg (x
);
1142 /* Expansion of REAL_VALUE_TRUNCATE.
1143 The result is in floating point, rounded to nearest or even. */
1146 real_value_truncate (mode
, arg
)
1147 enum machine_mode mode
;
1148 REAL_VALUE_TYPE arg
;
1150 UEMUSHORT e
[NE
], t
[NE
];
1162 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
1199 /* If an unsupported type was requested, presume that
1200 the machine files know something useful to do with
1201 the unmodified value. */
1210 /* Return true if ARG can be represented exactly in MODE. */
1213 exact_real_truncate (mode
, arg
)
1214 enum machine_mode mode
;
1215 REAL_VALUE_TYPE
*arg
;
1217 REAL_VALUE_TYPE trunc
;
1219 if (target_isnan (*arg
))
1222 trunc
= real_value_truncate (mode
, *arg
);
1223 return ereal_cmp (*arg
, trunc
) == 0;
1226 /* Try to change R into its exact multiplicative inverse in machine mode
1227 MODE. Return nonzero function value if successful. */
1230 exact_real_inverse (mode
, r
)
1231 enum machine_mode mode
;
1234 UEMUSHORT e
[NE
], einv
[NE
];
1235 REAL_VALUE_TYPE rinv
;
1240 /* Test for input in range. Don't transform IEEE special values. */
1241 if (eisinf (e
) || eisnan (e
) || (ecmp (e
, ezero
) == 0))
1244 /* Test for a power of 2: all significand bits zero except the MSB.
1245 We are assuming the target has binary (or hex) arithmetic. */
1246 if (e
[NE
- 2] != 0x8000)
1249 for (i
= 0; i
< NE
- 2; i
++)
1255 /* Compute the inverse and truncate it to the required mode. */
1256 ediv (e
, eone
, einv
);
1257 PUT_REAL (einv
, &rinv
);
1258 rinv
= real_value_truncate (mode
, rinv
);
1260 #ifdef CHECK_FLOAT_VALUE
1261 /* This check is not redundant. It may, for example, flush
1262 a supposedly IEEE denormal value to zero. */
1264 if (CHECK_FLOAT_VALUE (mode
, rinv
, i
))
1267 GET_REAL (&rinv
, einv
);
1269 /* Check the bits again, because the truncation might have
1270 generated an arbitrary saturation value on overflow. */
1271 if (einv
[NE
- 2] != 0x8000)
1274 for (i
= 0; i
< NE
- 2; i
++)
1280 /* Fail if the computed inverse is out of range. */
1281 if (eisinf (einv
) || eisnan (einv
) || (ecmp (einv
, ezero
) == 0))
1284 /* Output the reciprocal and return success flag. */
1289 /* Used for debugging--print the value of R in human-readable format
1298 REAL_VALUE_TO_DECIMAL (r
, "%.20g", dstr
);
1299 fprintf (stderr
, "%s", dstr
);
1303 /* The following routines convert REAL_VALUE_TYPE to the various floating
1304 point formats that are meaningful to supported computers.
1306 The results are returned in 32-bit pieces, each piece stored in a `long'.
1307 This is so they can be printed by statements like
1309 fprintf (file, "%lx, %lx", L[0], L[1]);
1311 that will work on both narrow- and wide-word host computers. */
1313 /* Convert R to a 128-bit long double precision value. The output array L
1314 contains four 32-bit pieces of the result, in the order they would appear
1325 #if INTEL_EXTENDED_IEEE_FORMAT == 0
1330 endian (e
, l
, TFmode
);
1333 /* Convert R to a double extended precision value. The output array L
1334 contains three 32-bit pieces of the result, in the order they would
1335 appear in memory. */
1346 endian (e
, l
, XFmode
);
1349 /* Convert R to a double precision value. The output array L contains two
1350 32-bit pieces of the result, in the order they would appear in memory. */
1361 endian (e
, l
, DFmode
);
1364 /* Convert R to a single precision float value stored in the least-significant
1365 bits of a `long'. */
1376 endian (e
, &l
, SFmode
);
1380 /* Convert X to a decimal ASCII string S for output to an assembly
1381 language file. Note, there is no standard way to spell infinity or
1382 a NaN, so these values may require special treatment in the tm.h
1386 ereal_to_decimal (x
, s
)
1396 /* Compare X and Y. Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
1397 or -2 if either is a NaN. */
1401 REAL_VALUE_TYPE x
, y
;
1403 UEMUSHORT ex
[NE
], ey
[NE
];
1407 return (ecmp (ex
, ey
));
1410 /* Return 1 if the sign bit of X is set, else return 0. */
1419 return (eisneg (ex
));
1424 Extended precision IEEE binary floating point arithmetic routines
1426 Numbers are stored in C language as arrays of 16-bit unsigned
1427 short integers. The arguments of the routines are pointers to
1430 External e type data structure, similar to Intel 8087 chip
1431 temporary real format but possibly with a larger significand:
1433 NE-1 significand words (least significant word first,
1434 most significant bit is normally set)
1435 exponent (value = EXONE for 1.0,
1436 top bit is the sign)
1439 Internal exploded e-type data structure of a number (a "word" is 16 bits):
1441 ei[0] sign word (0 for positive, 0xffff for negative)
1442 ei[1] biased exponent (value = EXONE for the number 1.0)
1443 ei[2] high guard word (always zero after normalization)
1445 to ei[NI-2] significand (NI-4 significand words,
1446 most significant word first,
1447 most significant bit is set)
1448 ei[NI-1] low guard word (0x8000 bit is rounding place)
1452 Routines for external format e-type numbers
1454 asctoe (string, e) ASCII string to extended double e type
1455 asctoe64 (string, &d) ASCII string to long double
1456 asctoe53 (string, &d) ASCII string to double
1457 asctoe24 (string, &f) ASCII string to single
1458 asctoeg (string, e, prec) ASCII string to specified precision
1459 e24toe (&f, e) IEEE single precision to e type
1460 e53toe (&d, e) IEEE double precision to e type
1461 e64toe (&d, e) IEEE long double precision to e type
1462 e113toe (&d, e) 128-bit long double precision to e type
1464 eabs (e) absolute value
1466 eadd (a, b, c) c = b + a
1468 ecmp (a, b) Returns 1 if a > b, 0 if a == b,
1469 -1 if a < b, -2 if either a or b is a NaN.
1470 ediv (a, b, c) c = b / a
1471 efloor (a, b) truncate to integer, toward -infinity
1472 efrexp (a, exp, s) extract exponent and significand
1473 eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
1474 euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
1475 einfin (e) set e to infinity, leaving its sign alone
1476 eldexp (a, n, b) multiply by 2**n
1478 emul (a, b, c) c = b * a
1481 eround (a, b) b = nearest integer value to a
1483 esub (a, b, c) c = b - a
1485 e24toasc (&f, str, n) single to ASCII string, n digits after decimal
1486 e53toasc (&d, str, n) double to ASCII string, n digits after decimal
1487 e64toasc (&d, str, n) 80-bit long double to ASCII string
1488 e113toasc (&d, str, n) 128-bit long double to ASCII string
1490 etoasc (e, str, n) e to ASCII string, n digits after decimal
1491 etoe24 (e, &f) convert e type to IEEE single precision
1492 etoe53 (e, &d) convert e type to IEEE double precision
1493 etoe64 (e, &d) convert e type to IEEE long double precision
1494 ltoe (&l, e) HOST_WIDE_INT to e type
1495 ultoe (&l, e) unsigned HOST_WIDE_INT to e type
1496 eisneg (e) 1 if sign bit of e != 0, else 0
1497 eisinf (e) 1 if e has maximum exponent (non-IEEE)
1498 or is infinite (IEEE)
1499 eisnan (e) 1 if e is a NaN
1502 Routines for internal format exploded e-type numbers
1504 eaddm (ai, bi) add significands, bi = bi + ai
1506 ecleazs (ei) set ei = 0 but leave its sign alone
1507 ecmpm (ai, bi) compare significands, return 1, 0, or -1
1508 edivm (ai, bi) divide significands, bi = bi / ai
1509 emdnorm (ai,l,s,exp) normalize and round off
1510 emovi (a, ai) convert external a to internal ai
1511 emovo (ai, a) convert internal ai to external a
1512 emovz (ai, bi) bi = ai, low guard word of bi = 0
1513 emulm (ai, bi) multiply significands, bi = bi * ai
1514 enormlz (ei) left-justify the significand
1515 eshdn1 (ai) shift significand and guards down 1 bit
1516 eshdn8 (ai) shift down 8 bits
1517 eshdn6 (ai) shift down 16 bits
1518 eshift (ai, n) shift ai n bits up (or down if n < 0)
1519 eshup1 (ai) shift significand and guards up 1 bit
1520 eshup8 (ai) shift up 8 bits
1521 eshup6 (ai) shift up 16 bits
1522 esubm (ai, bi) subtract significands, bi = bi - ai
1523 eiisinf (ai) 1 if infinite
1524 eiisnan (ai) 1 if a NaN
1525 eiisneg (ai) 1 if sign bit of ai != 0, else 0
1526 einan (ai) set ai = NaN
1528 eiinfin (ai) set ai = infinity
1531 The result is always normalized and rounded to NI-4 word precision
1532 after each arithmetic operation.
1534 Exception flags are NOT fully supported.
1536 Signaling NaN's are NOT supported; they are treated the same
1539 Define INFINITY for support of infinity; otherwise a
1540 saturation arithmetic is implemented.
1542 Define NANS for support of Not-a-Number items; otherwise the
1543 arithmetic will never produce a NaN output, and might be confused
1545 If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1546 either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1547 may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1550 Denormals are always supported here where appropriate (e.g., not
1551 for conversion to DEC numbers). */
1553 /* Definitions for error codes that are passed to the common error handling
1556 For Digital Equipment PDP-11 and VAX computers, certain
1557 IBM systems, and others that use numbers with a 56-bit
1558 significand, the symbol DEC should be defined. In this
1559 mode, most floating point constants are given as arrays
1560 of octal integers to eliminate decimal to binary conversion
1561 errors that might be introduced by the compiler.
1563 For computers, such as IBM PC, that follow the IEEE
1564 Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1565 Std 754-1985), the symbol IEEE should be defined.
1566 These numbers have 53-bit significands. In this mode, constants
1567 are provided as arrays of hexadecimal 16 bit integers.
1568 The endian-ness of generated values is controlled by
1569 REAL_WORDS_BIG_ENDIAN.
1571 To accommodate other types of computer arithmetic, all
1572 constants are also provided in a normal decimal radix
1573 which one can hope are correctly converted to a suitable
1574 format by the available C language compiler. To invoke
1575 this mode, the symbol UNK is defined.
1577 An important difference among these modes is a predefined
1578 set of machine arithmetic constants for each. The numbers
1579 MACHEP (the machine roundoff error), MAXNUM (largest number
1580 represented), and several other parameters are preset by
1581 the configuration symbol. Check the file const.c to
1582 ensure that these values are correct for your computer.
1584 For ANSI C compatibility, define ANSIC equal to 1. Currently
1585 this affects only the atan2 function and others that use it. */
1587 /* Constant definitions for math error conditions. */
1589 #define DOMAIN 1 /* argument domain error */
1590 #define SING 2 /* argument singularity */
1591 #define OVERFLOW 3 /* overflow range error */
1592 #define UNDERFLOW 4 /* underflow range error */
1593 #define TLOSS 5 /* total loss of precision */
1594 #define PLOSS 6 /* partial loss of precision */
1595 #define INVALID 7 /* NaN-producing operation */
1597 /* e type constants used by high precision check routines */
1599 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
1601 const UEMUSHORT ezero
[NE
] =
1602 {0x0000, 0x0000, 0x0000, 0x0000,
1603 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
1606 const UEMUSHORT ehalf
[NE
] =
1607 {0x0000, 0x0000, 0x0000, 0x0000,
1608 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
1611 const UEMUSHORT eone
[NE
] =
1612 {0x0000, 0x0000, 0x0000, 0x0000,
1613 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
1616 const UEMUSHORT etwo
[NE
] =
1617 {0x0000, 0x0000, 0x0000, 0x0000,
1618 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
1621 const UEMUSHORT e32
[NE
] =
1622 {0x0000, 0x0000, 0x0000, 0x0000,
1623 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
1625 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1626 const UEMUSHORT elog2
[NE
] =
1627 {0x40f3, 0xf6af, 0x03f2, 0xb398,
1628 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1630 /* 1.41421356237309504880168872420969807856967187537695E0 */
1631 const UEMUSHORT esqrt2
[NE
] =
1632 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1633 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1635 /* 3.14159265358979323846264338327950288419716939937511E0 */
1636 const UEMUSHORT epi
[NE
] =
1637 {0x2902, 0x1cd1, 0x80dc, 0x628b,
1638 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1641 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1642 const UEMUSHORT ezero
[NE
] =
1643 {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1644 const UEMUSHORT ehalf
[NE
] =
1645 {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1646 const UEMUSHORT eone
[NE
] =
1647 {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1648 const UEMUSHORT etwo
[NE
] =
1649 {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1650 const UEMUSHORT e32
[NE
] =
1651 {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1652 const UEMUSHORT elog2
[NE
] =
1653 {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1654 const UEMUSHORT esqrt2
[NE
] =
1655 {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1656 const UEMUSHORT epi
[NE
] =
1657 {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1660 /* Control register for rounding precision.
1661 This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits. */
1666 /* Clear out entire e-type number X. */
1674 for (i
= 0; i
< NE
; i
++)
1678 /* Move e-type number from A to B. */
1687 for (i
= 0; i
< NE
; i
++)
1693 /* Absolute value of e-type X. */
1699 /* sign is top bit of last word of external format */
1700 x
[NE
- 1] &= 0x7fff;
1704 /* Negate the e-type number X. */
1711 x
[NE
- 1] ^= 0x8000; /* Toggle the sign bit */
1714 /* Return 1 if sign bit of e-type number X is nonzero, else zero. */
1718 const UEMUSHORT x
[];
1721 if (x
[NE
- 1] & 0x8000)
1727 /* Return 1 if e-type number X is infinity, else return zero. */
1731 const UEMUSHORT x
[];
1738 if ((x
[NE
- 1] & 0x7fff) == 0x7fff)
1744 /* Check if e-type number is not a number. The bit pattern is one that we
1745 defined, so we know for sure how to detect it. */
1749 const UEMUSHORT x
[] ATTRIBUTE_UNUSED
;
1754 /* NaN has maximum exponent */
1755 if ((x
[NE
- 1] & 0x7fff) != 0x7fff)
1757 /* ... and non-zero significand field. */
1758 for (i
= 0; i
< NE
- 1; i
++)
1768 /* Fill e-type number X with infinity pattern (IEEE)
1769 or largest possible number (non-IEEE). */
1778 for (i
= 0; i
< NE
- 1; i
++)
1782 for (i
= 0; i
< NE
- 1; i
++)
1810 /* Output an e-type NaN.
1811 This generates Intel's quiet NaN pattern for extended real.
1812 The exponent is 7fff, the leading mantissa word is c000. */
1822 for (i
= 0; i
< NE
- 2; i
++)
1825 *x
= (sign
<< 15) | 0x7fff;
1829 /* Move in an e-type number A, converting it to exploded e-type B. */
1841 p
= a
+ (NE
- 1); /* point to last word of external number */
1842 /* get the sign bit */
1847 /* get the exponent */
1849 *q
++ &= 0x7fff; /* delete the sign bit */
1851 if ((*(q
- 1) & 0x7fff) == 0x7fff)
1857 for (i
= 3; i
< NI
; i
++)
1863 for (i
= 2; i
< NI
; i
++)
1869 /* clear high guard word */
1871 /* move in the significand */
1872 for (i
= 0; i
< NE
- 1; i
++)
1874 /* clear low guard word */
1878 /* Move out exploded e-type number A, converting it to e type B. */
1891 q
= b
+ (NE
- 1); /* point to output exponent */
1892 /* combine sign and exponent */
1895 *q
-- = *p
++ | 0x8000;
1899 if (*(p
- 1) == 0x7fff)
1904 enan (b
, eiisneg (a
));
1912 /* skip over guard word */
1914 /* move the significand */
1915 for (j
= 0; j
< NE
- 1; j
++)
1919 /* Clear out exploded e-type number XI. */
1927 for (i
= 0; i
< NI
; i
++)
1931 /* Clear out exploded e-type XI, but don't touch the sign. */
1940 for (i
= 0; i
< NI
- 1; i
++)
1944 /* Move exploded e-type number from A to B. */
1953 for (i
= 0; i
< NI
- 1; i
++)
1955 /* clear low guard word */
1959 /* Generate exploded e-type NaN.
1960 The explicit pattern for this is maximum exponent and
1961 top two significant bits set. */
1975 /* Return nonzero if exploded e-type X is a NaN. */
1980 const UEMUSHORT x
[];
1984 if ((x
[E
] & 0x7fff) == 0x7fff)
1986 for (i
= M
+ 1; i
< NI
; i
++)
1996 /* Return nonzero if sign of exploded e-type X is nonzero. */
2000 const UEMUSHORT x
[];
2007 /* Fill exploded e-type X with infinity pattern.
2008 This has maximum exponent and significand all zeros. */
2020 /* Return nonzero if exploded e-type X is infinite. */
2025 const UEMUSHORT x
[];
2032 if ((x
[E
] & 0x7fff) == 0x7fff)
2036 #endif /* INFINITY */
2038 /* Compare significands of numbers in internal exploded e-type format.
2039 Guard words are included in the comparison.
2047 const UEMUSHORT
*a
, *b
;
2051 a
+= M
; /* skip up to significand area */
2053 for (i
= M
; i
< NI
; i
++)
2061 if (*(--a
) > *(--b
))
2067 /* Shift significand of exploded e-type X down by 1 bit. */
2076 x
+= M
; /* point to significand area */
2079 for (i
= M
; i
< NI
; i
++)
2091 /* Shift significand of exploded e-type X up by 1 bit. */
2103 for (i
= M
; i
< NI
; i
++)
2116 /* Shift significand of exploded e-type X down by 8 bits. */
2122 UEMUSHORT newbyt
, oldbyt
;
2127 for (i
= M
; i
< NI
; i
++)
2137 /* Shift significand of exploded e-type X up by 8 bits. */
2144 UEMUSHORT newbyt
, oldbyt
;
2149 for (i
= M
; i
< NI
; i
++)
2159 /* Shift significand of exploded e-type X up by 16 bits. */
2171 for (i
= M
; i
< NI
- 1; i
++)
2177 /* Shift significand of exploded e-type X down by 16 bits. */
2189 for (i
= M
; i
< NI
- 1; i
++)
2195 /* Add significands of exploded e-type X and Y. X + Y replaces Y. */
2209 for (i
= M
; i
< NI
; i
++)
2211 a
= (unsigned EMULONG
) (*x
) + (unsigned EMULONG
) (*y
) + carry
;
2222 /* Subtract significands of exploded e-type X and Y. Y - X replaces Y. */
2236 for (i
= M
; i
< NI
; i
++)
2238 a
= (unsigned EMULONG
) (*y
) - (unsigned EMULONG
) (*x
) - carry
;
2250 static UEMUSHORT equot
[NI
];
2254 /* Radix 2 shift-and-add versions of multiply and divide */
2257 /* Divide significands */
2261 UEMUSHORT den
[], num
[];
2271 for (i
= M
; i
< NI
; i
++)
2276 /* Use faster compare and subtraction if denominator has only 15 bits of
2282 for (i
= M
+ 3; i
< NI
; i
++)
2287 if ((den
[M
+ 1] & 1) != 0)
2295 for (i
= 0; i
< NBITS
+ 2; i
++)
2313 /* The number of quotient bits to calculate is NBITS + 1 scaling guard
2314 bit + 1 roundoff bit. */
2319 for (i
= 0; i
< NBITS
+ 2; i
++)
2321 if (ecmpm (den
, num
) <= 0)
2324 j
= 1; /* quotient bit = 1 */
2338 /* test for nonzero remainder after roundoff bit */
2341 for (i
= M
; i
< NI
; i
++)
2349 for (i
= 0; i
< NI
; i
++)
2355 /* Multiply significands */
2366 for (i
= M
; i
< NI
; i
++)
2371 while (*p
== 0) /* significand is not supposed to be zero */
2376 if ((*p
& 0xff) == 0)
2384 for (i
= 0; i
< k
; i
++)
2388 /* remember if there were any nonzero bits shifted out */
2395 for (i
= 0; i
< NI
; i
++)
2398 /* return flag for lost nonzero bits */
2404 /* Radix 65536 versions of multiply and divide. */
2406 /* Multiply significand of e-type number B
2407 by 16-bit quantity A, return e-type result to C. */
2412 const UEMUSHORT b
[];
2416 unsigned EMULONG carry
;
2417 const UEMUSHORT
*ps
;
2419 unsigned EMULONG aa
, m
;
2428 for (i
=M
+1; i
<NI
; i
++)
2438 m
= (unsigned EMULONG
) aa
* *ps
--;
2439 carry
= (m
& 0xffff) + *pp
;
2440 *pp
-- = (UEMUSHORT
) carry
;
2441 carry
= (carry
>> 16) + (m
>> 16) + *pp
;
2442 *pp
= (UEMUSHORT
) carry
;
2443 *(pp
-1) = carry
>> 16;
2446 for (i
=M
; i
<NI
; i
++)
2450 /* Divide significands of exploded e-types NUM / DEN. Neither the
2451 numerator NUM nor the denominator DEN is permitted to have its high guard
2456 const UEMUSHORT den
[];
2461 unsigned EMULONG tnum
;
2462 UEMUSHORT j
, tdenm
, tquot
;
2463 UEMUSHORT tprod
[NI
+1];
2469 for (i
=M
; i
<NI
; i
++)
2475 for (i
=M
; i
<NI
; i
++)
2477 /* Find trial quotient digit (the radix is 65536). */
2478 tnum
= (((unsigned EMULONG
) num
[M
]) << 16) + num
[M
+1];
2480 /* Do not execute the divide instruction if it will overflow. */
2481 if ((tdenm
* (unsigned long) 0xffff) < tnum
)
2484 tquot
= tnum
/ tdenm
;
2485 /* Multiply denominator by trial quotient digit. */
2486 m16m ((unsigned int) tquot
, den
, tprod
);
2487 /* The quotient digit may have been overestimated. */
2488 if (ecmpm (tprod
, num
) > 0)
2492 if (ecmpm (tprod
, num
) > 0)
2502 /* test for nonzero remainder after roundoff bit */
2505 for (i
=M
; i
<NI
; i
++)
2512 for (i
=0; i
<NI
; i
++)
2518 /* Multiply significands of exploded e-type A and B, result in B. */
2522 const UEMUSHORT a
[];
2527 UEMUSHORT pprod
[NI
];
2533 for (i
=M
; i
<NI
; i
++)
2539 for (i
=M
+1; i
<NI
; i
++)
2547 m16m ((unsigned int) *p
--, b
, pprod
);
2548 eaddm (pprod
, equot
);
2554 for (i
=0; i
<NI
; i
++)
2557 /* return flag for lost nonzero bits */
2563 /* Normalize and round off.
2565 The internal format number to be rounded is S.
2566 Input LOST is 0 if the value is exact. This is the so-called sticky bit.
2568 Input SUBFLG indicates whether the number was obtained
2569 by a subtraction operation. In that case if LOST is nonzero
2570 then the number is slightly smaller than indicated.
2572 Input EXP is the biased exponent, which may be negative.
2573 the exponent field of S is ignored but is replaced by
2574 EXP as adjusted by normalization and rounding.
2576 Input RCNTRL is the rounding control. If it is nonzero, the
2577 returned value will be rounded to RNDPRC bits.
2579 For future reference: In order for emdnorm to round off denormal
2580 significands at the right point, the input exponent must be
2581 adjusted to be the actual value it would have after conversion to
2582 the final floating point type. This adjustment has been
2583 implemented for all type conversions (etoe53, etc.) and decimal
2584 conversions, but not for the arithmetic functions (eadd, etc.).
2585 Data types having standard 15-bit exponents are not affected by
2586 this, but SFmode and DFmode are affected. For example, ediv with
2587 rndprc = 24 will not round correctly to 24-bit precision if the
2588 result is denormal. */
2590 static int rlast
= -1;
2592 static UEMUSHORT rmsk
= 0;
2593 static UEMUSHORT rmbit
= 0;
2594 static UEMUSHORT rebit
= 0;
2596 static UEMUSHORT rbit
[NI
];
2599 emdnorm (s
, lost
, subflg
, exp
, rcntrl
)
2602 int subflg ATTRIBUTE_UNUSED
;
2612 /* a blank significand could mean either zero or infinity. */
2625 if ((j
> NBITS
) && (exp
< 32767))
2633 if (exp
> (EMULONG
) (-NBITS
- 1))
2646 /* Round off, unless told not to by rcntrl. */
2649 /* Set up rounding parameters if the control register changed. */
2650 if (rndprc
!= rlast
)
2657 rw
= NI
- 1; /* low guard word */
2680 /* For DEC or IBM arithmetic */
2697 /* For C4x arithmetic */
2718 /* Shift down 1 temporarily if the data structure has an implied
2719 most significant bit and the number is denormal.
2720 Intel long double denormals also lose one bit of precision. */
2721 if ((exp
<= 0) && (rndprc
!= NBITS
)
2722 && ((rndprc
!= 64) || ((rndprc
== 64) && ! REAL_WORDS_BIG_ENDIAN
)))
2724 lost
|= s
[NI
- 1] & 1;
2727 /* Clear out all bits below the rounding bit,
2728 remembering in r if any were nonzero. */
2742 if ((r
& rmbit
) != 0)
2748 { /* round to even */
2749 if ((s
[re
] & rebit
) == 0)
2764 /* Undo the temporary shift for denormal values. */
2765 if ((exp
<= 0) && (rndprc
!= NBITS
)
2766 && ((rndprc
!= 64) || ((rndprc
== 64) && ! REAL_WORDS_BIG_ENDIAN
)))
2771 { /* overflow on roundoff */
2784 for (i
= 2; i
< NI
- 1; i
++)
2787 warning ("floating point overflow");
2791 for (i
= M
+ 1; i
< NI
- 1; i
++)
2794 if ((rndprc
< 64) || (rndprc
== 113))
2809 s
[1] = (UEMUSHORT
) exp
;
2812 /* Subtract. C = B - A, all e type numbers. */
2814 static int subflg
= 0;
2818 const UEMUSHORT
*a
, *b
;
2833 /* Infinity minus infinity is a NaN.
2834 Test for subtracting infinities of the same sign. */
2835 if (eisinf (a
) && eisinf (b
)
2836 && ((eisneg (a
) ^ eisneg (b
)) == 0))
2838 mtherr ("esub", INVALID
);
2847 /* Add. C = A + B, all e type. */
2851 const UEMUSHORT
*a
, *b
;
2856 /* NaN plus anything is a NaN. */
2867 /* Infinity minus infinity is a NaN.
2868 Test for adding infinities of opposite signs. */
2869 if (eisinf (a
) && eisinf (b
)
2870 && ((eisneg (a
) ^ eisneg (b
)) != 0))
2872 mtherr ("esub", INVALID
);
2881 /* Arithmetic common to both addition and subtraction. */
2885 const UEMUSHORT
*a
, *b
;
2888 UEMUSHORT ai
[NI
], bi
[NI
], ci
[NI
];
2890 EMULONG lt
, lta
, ltb
;
2911 /* compare exponents */
2916 { /* put the larger number in bi */
2926 if (lt
< (EMULONG
) (-NBITS
- 1))
2927 goto done
; /* answer same as larger addend */
2929 lost
= eshift (ai
, k
); /* shift the smaller number down */
2933 /* exponents were the same, so must compare significands */
2936 { /* the numbers are identical in magnitude */
2937 /* if different signs, result is zero */
2943 /* if same sign, result is double */
2944 /* double denormalized tiny number */
2945 if ((bi
[E
] == 0) && ((bi
[3] & 0x8000) == 0))
2950 /* add 1 to exponent unless both are zero! */
2951 for (j
= 1; j
< NI
- 1; j
++)
2967 bi
[E
] = (UEMUSHORT
) ltb
;
2971 { /* put the larger number in bi */
2987 emdnorm (bi
, lost
, subflg
, ltb
, !ROUND_TOWARDS_ZERO
);
2993 /* Divide: C = B/A, all e type. */
2997 const UEMUSHORT
*a
, *b
;
3000 UEMUSHORT ai
[NI
], bi
[NI
];
3002 EMULONG lt
, lta
, ltb
;
3004 /* IEEE says if result is not a NaN, the sign is "-" if and only if
3005 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
3006 sign
= eisneg (a
) ^ eisneg (b
);
3009 /* Return any NaN input. */
3020 /* Zero over zero, or infinity over infinity, is a NaN. */
3021 if (((ecmp (a
, ezero
) == 0) && (ecmp (b
, ezero
) == 0))
3022 || (eisinf (a
) && eisinf (b
)))
3024 mtherr ("ediv", INVALID
);
3029 /* Infinity over anything else is infinity. */
3036 /* Anything else over infinity is zero. */
3048 { /* See if numerator is zero. */
3049 for (i
= 1; i
< NI
- 1; i
++)
3053 ltb
-= enormlz (bi
);
3063 { /* possible divide by zero */
3064 for (i
= 1; i
< NI
- 1; i
++)
3068 lta
-= enormlz (ai
);
3072 /* Divide by zero is not an invalid operation.
3073 It is a divide-by-zero operation! */
3075 mtherr ("ediv", SING
);
3081 /* calculate exponent */
3082 lt
= ltb
- lta
+ EXONE
;
3083 emdnorm (bi
, i
, 0, lt
, !ROUND_TOWARDS_ZERO
);
3090 && (ecmp (c
, ezero
) != 0)
3093 *(c
+(NE
-1)) |= 0x8000;
3095 *(c
+(NE
-1)) &= ~0x8000;
3098 /* Multiply e-types A and B, return e-type product C. */
3102 const UEMUSHORT
*a
, *b
;
3105 UEMUSHORT ai
[NI
], bi
[NI
];
3107 EMULONG lt
, lta
, ltb
;
3109 /* IEEE says if result is not a NaN, the sign is "-" if and only if
3110 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
3111 sign
= eisneg (a
) ^ eisneg (b
);
3114 /* NaN times anything is the same NaN. */
3125 /* Zero times infinity is a NaN. */
3126 if ((eisinf (a
) && (ecmp (b
, ezero
) == 0))
3127 || (eisinf (b
) && (ecmp (a
, ezero
) == 0)))
3129 mtherr ("emul", INVALID
);
3134 /* Infinity times anything else is infinity. */
3136 if (eisinf (a
) || eisinf (b
))
3148 for (i
= 1; i
< NI
- 1; i
++)
3152 lta
-= enormlz (ai
);
3163 for (i
= 1; i
< NI
- 1; i
++)
3167 ltb
-= enormlz (bi
);
3176 /* Multiply significands */
3178 /* calculate exponent */
3179 lt
= lta
+ ltb
- (EXONE
- 1);
3180 emdnorm (bi
, j
, 0, lt
, !ROUND_TOWARDS_ZERO
);
3187 && (ecmp (c
, ezero
) != 0)
3190 *(c
+(NE
-1)) |= 0x8000;
3192 *(c
+(NE
-1)) &= ~0x8000;
3195 /* Convert double precision PE to e-type Y. */
3199 const UEMUSHORT
*pe
;
3209 ibmtoe (pe
, y
, DFmode
);
3214 c4xtoe (pe
, y
, HFmode
);
3218 ieeetoe (pe
, y
, &ieee_53
);
3220 #endif /* not C4X */
3221 #endif /* not IBM */
3222 #endif /* not DEC */
3225 /* Convert double extended precision float PE to e type Y. */
3229 const UEMUSHORT
*pe
;
3239 for (i
= 0; i
< NE
- 5; i
++)
3242 /* REAL_WORDS_BIG_ENDIAN is always 0 for DEC and 1 for IBM.
3243 This precision is not ordinarily supported on DEC or IBM. */
3244 if (! REAL_WORDS_BIG_ENDIAN
)
3246 for (i
= 0; i
< 5; i
++)
3250 /* For denormal long double Intel format, shift significand up one
3251 -- but only if the top significand bit is zero. A top bit of 1
3252 is "pseudodenormal" when the exponent is zero. */
3253 if ((yy
[NE
-1] & 0x7fff) == 0 && (yy
[NE
-2] & 0x8000) == 0)
3266 p
= &yy
[0] + (NE
- 1);
3269 for (i
= 0; i
< 4; i
++)
3272 #endif /* not C4X */
3274 /* Point to the exponent field and check max exponent cases. */
3276 if ((*p
& 0x7fff) == 0x7fff)
3279 if (! REAL_WORDS_BIG_ENDIAN
)
3281 for (i
= 0; i
< 4; i
++)
3283 if ((i
!= 3 && pe
[i
] != 0)
3284 /* Anything but 0x8000 here, including 0, is a NaN. */
3285 || (i
== 3 && pe
[i
] != 0x8000))
3287 enan (y
, (*p
& 0x8000) != 0);
3294 /* In Motorola extended precision format, the most significant
3295 bit of an infinity mantissa could be either 1 or 0. It is
3296 the lower order bits that tell whether the value is a NaN. */
3297 if ((pe
[2] & 0x7fff) != 0)
3300 for (i
= 3; i
<= 5; i
++)
3305 enan (y
, (*p
& 0x8000) != 0);
3317 #endif /* INFINITY */
3320 for (i
= 0; i
< NE
; i
++)
3324 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3325 /* Convert 128-bit long double precision float PE to e type Y. */
3329 const UEMUSHORT
*pe
;
3332 ieeetoe (pe
, y
, &ieee_113
);
3334 #endif /* INTEL_EXTENDED_IEEE_FORMAT == 0 */
3336 /* Convert single precision float PE to e type Y. */
3340 const UEMUSHORT
*pe
;
3345 ibmtoe (pe
, y
, SFmode
);
3351 c4xtoe (pe
, y
, QFmode
);
3356 ieeetoe (pe
, y
, &dec_f
);
3360 ieeetoe (pe
, y
, &ieee_24
);
3362 #endif /* not DEC */
3363 #endif /* not C4X */
3364 #endif /* not IBM */
3367 /* Convert machine format float of specified format PE to e type Y. */
3370 ieeetoe (pe
, y
, fmt
)
3371 const UEMUSHORT
*pe
;
3373 const struct ieee_format
*fmt
;
3380 int shortsm1
= fmt
->bits
/ 16 - 1;
3382 int expmask
= (1 << fmt
->expbits
) - 1;
3384 int expshift
= (fmt
->precision
- 1) & 0x0f;
3385 int highbit
= 1 << expshift
;
3390 if (! REAL_WORDS_BIG_ENDIAN
)
3396 yy
[M
] = (r
& (highbit
- 1)) | highbit
;
3397 r
= (r
& 0x7fff) >> expshift
;
3399 if (!LARGEST_EXPONENT_IS_NORMAL (fmt
->precision
) && r
== expmask
)
3402 /* First check the word where high order mantissa and exponent live */
3403 if ((*e
& (highbit
- 1)) != 0)
3405 enan (y
, yy
[0] != 0);
3408 if (! REAL_WORDS_BIG_ENDIAN
)
3410 for (i
= 0; i
< shortsm1
; i
++)
3414 enan (y
, yy
[0] != 0);
3421 for (i
= 1; i
< shortsm1
+ 1; i
++)
3425 enan (y
, yy
[0] != 0);
3437 #endif /* INFINITY */
3438 /* If zero exponent, then the significand is denormalized.
3439 So take back the understood high significand bit. */
3445 r
+= fmt
->adjustment
;
3448 if (! REAL_WORDS_BIG_ENDIAN
)
3450 for (i
= 0; i
< shortsm1
; i
++)
3456 for (i
= 0; i
< shortsm1
; i
++)
3459 if (fmt
->precision
== 113)
3461 /* denorm is left alone in 113 bit format */
3467 eshift (yy
, -(expshift
+ 1));
3469 { /* if zero exponent, then normalize the significand */
3470 if ((k
= enormlz (yy
)) > NBITS
)
3473 yy
[E
] -= (UEMUSHORT
) (k
- 1);
3479 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3480 /* Convert e-type X to IEEE 128-bit long double format E. */
3487 etoieee (x
, e
, &ieee_113
);
3490 /* Convert exploded e-type X, that has already been rounded to
3491 113-bit precision, to IEEE 128-bit long double format Y. */
3497 toieee (x
, y
, &ieee_113
);
3500 #endif /* INTEL_EXTENDED_IEEE_FORMAT == 0 */
3502 /* Convert e-type X to IEEE double extended format E. */
3509 etoieee (x
, e
, &ieee_64
);
3512 /* Convert exploded e-type X, that has already been rounded to
3513 64-bit precision, to IEEE double extended format Y. */
3519 toieee (x
, y
, &ieee_64
);
3522 /* e type to double precision. */
3525 /* Convert e-type X to DEC-format double E. */
3535 /* Convert exploded e-type X, that has already been rounded to
3536 56-bit double precision, to DEC double Y. */
3547 /* Convert e-type X to IBM 370-format double E. */
3554 etoibm (x
, e
, DFmode
);
3557 /* Convert exploded e-type X, that has already been rounded to
3558 56-bit precision, to IBM 370 double Y. */
3564 toibm (x
, y
, DFmode
);
3567 #else /* it's neither DEC nor IBM */
3569 /* Convert e-type X to C4X-format long double E. */
3576 etoc4x (x
, e
, HFmode
);
3579 /* Convert exploded e-type X, that has already been rounded to
3580 56-bit precision, to IBM 370 double Y. */
3586 toc4x (x
, y
, HFmode
);
3589 #else /* it's neither DEC nor IBM nor C4X */
3591 /* Convert e-type X to IEEE double E. */
3598 etoieee (x
, e
, &ieee_53
);
3601 /* Convert exploded e-type X, that has already been rounded to
3602 53-bit precision, to IEEE double Y. */
3608 toieee (x
, y
, &ieee_53
);
3611 #endif /* not C4X */
3612 #endif /* not IBM */
3613 #endif /* not DEC */
3617 /* e type to single precision. */
3620 /* Convert e-type X to IBM 370 float E. */
3627 etoibm (x
, e
, SFmode
);
3630 /* Convert exploded e-type X, that has already been rounded to
3631 float precision, to IBM 370 float Y. */
3637 toibm (x
, y
, SFmode
);
3640 #else /* it's not IBM */
3643 /* Convert e-type X to C4X float E. */
3650 etoc4x (x
, e
, QFmode
);
3653 /* Convert exploded e-type X, that has already been rounded to
3654 float precision, to IBM 370 float Y. */
3660 toc4x (x
, y
, QFmode
);
3663 #else /* it's neither IBM nor C4X */
3667 /* Convert e-type X to DEC F-float E. */
3674 etoieee (x
, e
, &dec_f
);
3677 /* Convert exploded e-type X, that has already been rounded to
3678 float precision, to DEC F-float Y. */
3684 toieee (x
, y
, &dec_f
);
3689 /* Convert e-type X to IEEE float E. */
3696 etoieee (x
, e
, &ieee_24
);
3699 /* Convert exploded e-type X, that has already been rounded to
3700 float precision, to IEEE float Y. */
3706 toieee (x
, y
, &ieee_24
);
3709 #endif /* not DEC */
3710 #endif /* not C4X */
3711 #endif /* not IBM */
3714 /* Convert e-type X to the IEEE format described by FMT. */
3720 const struct ieee_format
*fmt
;
3729 make_nan (e
, eisneg (x
), fmt
->mode
);
3740 /* Adjust exponent for offset. */
3741 exp
= (EMULONG
) xi
[E
] - fmt
->adjustment
;
3743 /* Round off to nearest or even. */
3745 rndprc
= fmt
->precision
;
3746 emdnorm (xi
, 0, 0, exp
, !ROUND_TOWARDS_ZERO
);
3751 toieee (xi
, e
, fmt
);
3754 /* Convert exploded e-type X, that has already been rounded to
3755 the necessary precision, to the IEEE format described by FMT. */
3760 const struct ieee_format
*fmt
;
3767 maxexp
= (1 << fmt
->expbits
) - 1;
3768 words
= (fmt
->bits
- fmt
->expbits
) / EMUSHORT_SIZE
;
3773 make_nan (y
, eiisneg (x
), fmt
->mode
);
3778 if (fmt
->expbits
< 15
3779 && LARGEST_EXPONENT_IS_NORMAL (fmt
->bits
)
3782 saturate (y
, eiisneg (x
), fmt
->bits
, 1);
3786 /* Point to the exponent. */
3787 if (REAL_WORDS_BIG_ENDIAN
)
3792 /* Copy the sign. */
3798 if (fmt
->expbits
< 15
3799 && !LARGEST_EXPONENT_IS_NORMAL (fmt
->bits
)
3802 /* Saturate at largest number less that infinity. */
3805 *q
|= maxexp
<< (15 - fmt
->expbits
);
3808 *q
|= (maxexp
<< (15 - fmt
->expbits
)) - 1;
3812 if (!REAL_WORDS_BIG_ENDIAN
)
3814 for (i
= 0; i
< words
; i
++)
3819 for (i
= 0; i
< words
; i
++)
3822 #if defined(INFINITY) && defined(ERANGE)
3828 /* If denormal and DEC float, return zero (DEC has no denormals) */
3832 for (i
= 0; i
< fmt
->bits
/ EMUSHORT_SIZE
; i
++)
3838 /* Delete the implied bit unless denormal, except for
3839 64-bit precision. */
3840 if (fmt
->precision
!= 64 && x
[E
] != 0)
3845 /* Shift denormal double extended Intel format significand down
3847 if (fmt
->precision
== 64 && x
[E
] == 0 && ! REAL_WORDS_BIG_ENDIAN
)
3850 if (fmt
->expbits
< 15)
3852 /* Shift the significand. */
3853 eshift (x
, 15 - fmt
->expbits
);
3855 /* Combine the exponent and upper bits of the significand. */
3856 *q
|= x
[E
] << (15 - fmt
->expbits
);
3857 *q
|= x
[M
] & (UEMUSHORT
) ~((maxexp
<< (15 - fmt
->expbits
)) | 0x8000);
3861 /* Copy the exponent. */
3865 /* Add padding after the exponent. At the moment, this is only necessary for
3866 64-bit precision; in this case, the padding is 16 bits. */
3867 if (fmt
->precision
== 64)
3872 if (REAL_WORDS_BIG_ENDIAN
)
3876 /* Copy the significand. */
3877 if (REAL_WORDS_BIG_ENDIAN
)
3879 for (i
= 0; i
< words
; i
++)
3880 *(++q
) = x
[i
+ M
+ 1];
3883 else if (fmt
->precision
== 64 && eiisinf (x
))
3885 /* Intel double extended infinity significand. */
3894 for (i
= 0; i
< words
; i
++)
3895 *(--q
) = x
[i
+ M
+ 1];
3900 /* Compare two e type numbers.
3904 -2 if either a or b is a NaN. */
3908 const UEMUSHORT
*a
, *b
;
3910 UEMUSHORT ai
[NI
], bi
[NI
];
3916 if (eisnan (a
) || eisnan (b
))
3925 { /* the signs are different */
3927 for (i
= 1; i
< NI
- 1; i
++)
3941 /* both are the same sign */
3956 return (0); /* equality */
3960 if (*(--p
) > *(--q
))
3961 return (msign
); /* p is bigger */
3963 return (-msign
); /* p is littler */
3967 /* Find e-type nearest integer to X, as floor (X + 0.5). */
3979 /* Convert HOST_WIDE_INT LP to e type Y. */
3983 const HOST_WIDE_INT
*lp
;
3987 unsigned HOST_WIDE_INT ll
;
3993 /* make it positive */
3994 ll
= (unsigned HOST_WIDE_INT
) (-(*lp
));
3995 yi
[0] = 0xffff; /* put correct sign in the e type number */
3999 ll
= (unsigned HOST_WIDE_INT
) (*lp
);
4001 /* move the long integer to yi significand area */
4002 #if HOST_BITS_PER_WIDE_INT == 64
4003 yi
[M
] = (UEMUSHORT
) (ll
>> 48);
4004 yi
[M
+ 1] = (UEMUSHORT
) (ll
>> 32);
4005 yi
[M
+ 2] = (UEMUSHORT
) (ll
>> 16);
4006 yi
[M
+ 3] = (UEMUSHORT
) ll
;
4007 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
4009 yi
[M
] = (UEMUSHORT
) (ll
>> 16);
4010 yi
[M
+ 1] = (UEMUSHORT
) ll
;
4011 yi
[E
] = EXONE
+ 15; /* exponent if normalize shift count were 0 */
4014 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
4015 ecleaz (yi
); /* it was zero */
4017 yi
[E
] -= (UEMUSHORT
) k
;/* subtract shift count from exponent */
4018 emovo (yi
, y
); /* output the answer */
4021 /* Convert unsigned HOST_WIDE_INT LP to e type Y. */
4025 const unsigned HOST_WIDE_INT
*lp
;
4029 unsigned HOST_WIDE_INT ll
;
4035 /* move the long integer to ayi significand area */
4036 #if HOST_BITS_PER_WIDE_INT == 64
4037 yi
[M
] = (UEMUSHORT
) (ll
>> 48);
4038 yi
[M
+ 1] = (UEMUSHORT
) (ll
>> 32);
4039 yi
[M
+ 2] = (UEMUSHORT
) (ll
>> 16);
4040 yi
[M
+ 3] = (UEMUSHORT
) ll
;
4041 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
4043 yi
[M
] = (UEMUSHORT
) (ll
>> 16);
4044 yi
[M
+ 1] = (UEMUSHORT
) ll
;
4045 yi
[E
] = EXONE
+ 15; /* exponent if normalize shift count were 0 */
4048 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
4049 ecleaz (yi
); /* it was zero */
4051 yi
[E
] -= (UEMUSHORT
) k
; /* subtract shift count from exponent */
4052 emovo (yi
, y
); /* output the answer */
4056 /* Find signed HOST_WIDE_INT integer I and floating point fractional
4057 part FRAC of e-type (packed internal format) floating point input X.
4058 The integer output I has the sign of the input, except that
4059 positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
4060 The output e-type fraction FRAC is the positive fractional
4071 unsigned HOST_WIDE_INT ll
;
4074 k
= (int) xi
[E
] - (EXONE
- 1);
4077 /* if exponent <= 0, integer = 0 and real output is fraction */
4082 if (k
> (HOST_BITS_PER_WIDE_INT
- 1))
4084 /* long integer overflow: output large integer
4085 and correct fraction */
4087 *i
= ((unsigned HOST_WIDE_INT
) 1) << (HOST_BITS_PER_WIDE_INT
- 1);
4090 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
4091 /* In this case, let it overflow and convert as if unsigned. */
4092 euifrac (x
, &ll
, frac
);
4093 *i
= (HOST_WIDE_INT
) ll
;
4096 /* In other cases, return the largest positive integer. */
4097 *i
= (((unsigned HOST_WIDE_INT
) 1) << (HOST_BITS_PER_WIDE_INT
- 1)) - 1;
4102 warning ("overflow on truncation to integer");
4106 /* Shift more than 16 bits: first shift up k-16 mod 16,
4107 then shift up by 16's. */
4108 j
= k
- ((k
>> 4) << 4);
4115 ll
= (ll
<< 16) | xi
[M
];
4117 while ((k
-= 16) > 0);
4124 /* shift not more than 16 bits */
4126 *i
= (HOST_WIDE_INT
) xi
[M
] & 0xffff;
4133 if ((k
= enormlz (xi
)) > NBITS
)
4136 xi
[E
] -= (UEMUSHORT
) k
;
4142 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
4143 FRAC of e-type X. A negative input yields integer output = 0 but
4144 correct fraction. */
4147 euifrac (x
, i
, frac
)
4149 unsigned HOST_WIDE_INT
*i
;
4152 unsigned HOST_WIDE_INT ll
;
4157 k
= (int) xi
[E
] - (EXONE
- 1);
4160 /* if exponent <= 0, integer = 0 and argument is fraction */
4165 if (k
> HOST_BITS_PER_WIDE_INT
)
4167 /* Long integer overflow: output large integer
4168 and correct fraction.
4169 Note, the BSD MicroVAX compiler says that ~(0UL)
4170 is a syntax error. */
4174 warning ("overflow on truncation to unsigned integer");
4178 /* Shift more than 16 bits: first shift up k-16 mod 16,
4179 then shift up by 16's. */
4180 j
= k
- ((k
>> 4) << 4);
4187 ll
= (ll
<< 16) | xi
[M
];
4189 while ((k
-= 16) > 0);
4194 /* shift not more than 16 bits */
4196 *i
= (HOST_WIDE_INT
) xi
[M
] & 0xffff;
4199 if (xi
[0]) /* A negative value yields unsigned integer 0. */
4205 if ((k
= enormlz (xi
)) > NBITS
)
4208 xi
[E
] -= (UEMUSHORT
) k
;
4213 /* Shift the significand of exploded e-type X up or down by SC bits. */
4234 lost
|= *p
; /* remember lost bits */
4275 return ((int) lost
);
4278 /* Shift normalize the significand area of exploded e-type X.
4279 Return the shift count (up = positive). */
4294 return (0); /* already normalized */
4300 /* With guard word, there are NBITS+16 bits available.
4301 Return true if all are zero. */
4305 /* see if high byte is zero */
4306 while ((*p
& 0xff00) == 0)
4311 /* now shift 1 bit at a time */
4312 while ((*p
& 0x8000) == 0)
4318 mtherr ("enormlz", UNDERFLOW
);
4324 /* Normalize by shifting down out of the high guard word
4325 of the significand */
4340 mtherr ("enormlz", OVERFLOW
);
4347 /* Powers of ten used in decimal <-> binary conversions. */
4352 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
4353 static const UEMUSHORT etens
[NTEN
+ 1][NE
] =
4355 {0x6576, 0x4a92, 0x804a, 0x153f,
4356 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4357 {0x6a32, 0xce52, 0x329a, 0x28ce,
4358 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4359 {0x526c, 0x50ce, 0xf18b, 0x3d28,
4360 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4361 {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4362 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4363 {0x851e, 0xeab7, 0x98fe, 0x901b,
4364 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4365 {0x0235, 0x0137, 0x36b1, 0x336c,
4366 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4367 {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4368 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4369 {0x0000, 0x0000, 0x0000, 0x0000,
4370 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4371 {0x0000, 0x0000, 0x0000, 0x0000,
4372 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4373 {0x0000, 0x0000, 0x0000, 0x0000,
4374 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4375 {0x0000, 0x0000, 0x0000, 0x0000,
4376 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4377 {0x0000, 0x0000, 0x0000, 0x0000,
4378 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4379 {0x0000, 0x0000, 0x0000, 0x0000,
4380 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4383 static const UEMUSHORT emtens
[NTEN
+ 1][NE
] =
4385 {0x2030, 0xcffc, 0xa1c3, 0x8123,
4386 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4387 {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4388 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4389 {0xf53f, 0xf698, 0x6bd3, 0x0158,
4390 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4391 {0xe731, 0x04d4, 0xe3f2, 0xd332,
4392 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4393 {0xa23e, 0x5308, 0xfefb, 0x1155,
4394 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4395 {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4396 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4397 {0x2a20, 0x6224, 0x47b3, 0x98d7,
4398 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4399 {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4400 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4401 {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4402 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4403 {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4404 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4405 {0xc155, 0xa4a8, 0x404e, 0x6113,
4406 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4407 {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4408 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4409 {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4410 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4413 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4414 static const UEMUSHORT etens
[NTEN
+ 1][NE
] =
4416 {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4417 {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4418 {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4419 {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4420 {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4421 {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4422 {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4423 {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4424 {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4425 {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4426 {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4427 {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4428 {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4431 static const UEMUSHORT emtens
[NTEN
+ 1][NE
] =
4433 {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4434 {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4435 {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4436 {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4437 {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4438 {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4439 {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4440 {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4441 {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4442 {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4443 {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4444 {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4445 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4450 /* Convert float value X to ASCII string STRING with NDIG digits after
4451 the decimal point. */
4454 e24toasc (x
, string
, ndigs
)
4455 const UEMUSHORT x
[];
4462 etoasc (w
, string
, ndigs
);
4465 /* Convert double value X to ASCII string STRING with NDIG digits after
4466 the decimal point. */
4469 e53toasc (x
, string
, ndigs
)
4470 const UEMUSHORT x
[];
4477 etoasc (w
, string
, ndigs
);
4480 /* Convert double extended value X to ASCII string STRING with NDIG digits
4481 after the decimal point. */
4484 e64toasc (x
, string
, ndigs
)
4485 const UEMUSHORT x
[];
4492 etoasc (w
, string
, ndigs
);
4495 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4496 after the decimal point. */
4499 e113toasc (x
, string
, ndigs
)
4500 const UEMUSHORT x
[];
4507 etoasc (w
, string
, ndigs
);
4511 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4512 the decimal point. */
4514 static char wstring
[80]; /* working storage for ASCII output */
4517 etoasc (x
, string
, ndigs
)
4518 const UEMUSHORT x
[];
4523 UEMUSHORT y
[NI
], t
[NI
], u
[NI
], w
[NI
];
4524 const UEMUSHORT
*p
, *r
, *ten
;
4526 int i
, j
, k
, expon
, rndsav
;
4539 sprintf (wstring
, " NaN ");
4543 rndprc
= NBITS
; /* set to full precision */
4544 emov (x
, y
); /* retain external format */
4545 if (y
[NE
- 1] & 0x8000)
4548 y
[NE
- 1] &= 0x7fff;
4555 ten
= &etens
[NTEN
][0];
4557 /* Test for zero exponent */
4560 for (k
= 0; k
< NE
- 1; k
++)
4563 goto tnzro
; /* denormalized number */
4565 goto isone
; /* valid all zeros */
4569 /* Test for infinity. */
4570 if (y
[NE
- 1] == 0x7fff)
4573 sprintf (wstring
, " -Infinity ");
4575 sprintf (wstring
, " Infinity ");
4579 /* Test for exponent nonzero but significand denormalized.
4580 * This is an error condition.
4582 if ((y
[NE
- 1] != 0) && ((y
[NE
- 2] & 0x8000) == 0))
4584 mtherr ("etoasc", DOMAIN
);
4585 sprintf (wstring
, "NaN");
4589 /* Compare to 1.0 */
4598 { /* Number is greater than 1 */
4599 /* Convert significand to an integer and strip trailing decimal zeros. */
4601 u
[NE
- 1] = EXONE
+ NBITS
- 1;
4603 p
= &etens
[NTEN
- 4][0];
4609 for (j
= 0; j
< NE
- 1; j
++)
4622 /* Rescale from integer significand */
4623 u
[NE
- 1] += y
[NE
- 1] - (unsigned int) (EXONE
+ NBITS
- 1);
4625 /* Find power of 10 */
4629 /* An unordered compare result shouldn't happen here. */
4630 while (ecmp (ten
, u
) <= 0)
4632 if (ecmp (p
, u
) <= 0)
4645 { /* Number is less than 1.0 */
4646 /* Pad significand with trailing decimal zeros. */
4649 while ((y
[NE
- 2] & 0x8000) == 0)
4658 for (i
= 0; i
< NDEC
+ 1; i
++)
4660 if ((w
[NI
- 1] & 0x7) != 0)
4662 /* multiply by 10 */
4675 if (eone
[NE
- 1] <= u
[1])
4687 while (ecmp (eone
, w
) > 0)
4689 if (ecmp (p
, w
) >= 0)
4704 /* Find the first (leading) digit. */
4710 digit
= equot
[NI
- 1];
4711 while ((digit
== 0) && (ecmp (y
, ezero
) != 0))
4719 digit
= equot
[NI
- 1];
4727 /* Examine number of digits requested by caller. */
4745 *s
++ = (char) digit
+ '0';
4748 /* Generate digits after the decimal point. */
4749 for (k
= 0; k
<= ndigs
; k
++)
4751 /* multiply current number by 10, without normalizing */
4758 *s
++ = (char) equot
[NI
- 1] + '0';
4760 digit
= equot
[NI
- 1];
4763 /* round off the ASCII string */
4766 /* Test for critical rounding case in ASCII output. */
4770 if (ecmp (t
, ezero
) != 0)
4771 goto roun
; /* round to nearest */
4773 if ((*(s
- 1) & 1) == 0)
4774 goto doexp
; /* round to even */
4777 /* Round up and propagate carry-outs */
4781 /* Carry out to most significant digit? */
4788 /* Most significant digit carries to 10? */
4796 /* Round up and carry out from less significant digits */
4806 /* Strip trailing zeros, but leave at least one. */
4807 while (ss
[-1] == '0' && ss
[-2] != '.')
4809 sprintf (ss
, "e%d", expon
);
4812 /* copy out the working string */
4815 while (*ss
== ' ') /* strip possible leading space */
4817 while ((*s
++ = *ss
++) != '\0')
4822 /* Convert ASCII string to floating point.
4824 Numeric input is a free format decimal number of any length, with
4825 or without decimal point. Entering E after the number followed by an
4826 integer number causes the second number to be interpreted as a power of
4827 10 to be multiplied by the first number (i.e., "scientific" notation). */
4829 /* Convert ASCII string S to single precision float value Y. */
4840 /* Convert ASCII string S to double precision value Y. */
4847 #if defined(DEC) || defined(IBM)
4859 /* Convert ASCII string S to double extended value Y. */
4869 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
4870 /* Convert ASCII string S to 128-bit long double Y. */
4877 asctoeg (s
, y
, 113);
4881 /* Convert ASCII string S to e type Y. */
4888 asctoeg (s
, y
, NBITS
);
4891 /* Convert ASCII string SS to e type Y, with a specified rounding precision
4892 of OPREC bits. BASE is 16 for C99 hexadecimal floating constants. */
4895 asctoeg (ss
, y
, oprec
)
4900 UEMUSHORT yy
[NI
], xt
[NI
], tt
[NI
];
4901 int esign
, decflg
, sgnflg
, nexp
, exp
, prec
, lost
;
4902 int i
, k
, trail
, c
, rndsav
;
4905 char *sp
, *s
, *lstr
;
4908 /* Copy the input string. */
4909 lstr
= (char *) alloca (strlen (ss
) + 1);
4911 while (*ss
== ' ') /* skip leading spaces */
4915 while ((*sp
++ = *ss
++) != '\0')
4919 if (s
[0] == '0' && (s
[1] == 'x' || s
[1] == 'X'))
4926 rndprc
= NBITS
; /* Set to full precision */
4939 if ((k
>= 0) && (k
< base
))
4941 /* Ignore leading zeros */
4942 if ((prec
== 0) && (decflg
== 0) && (k
== 0))
4944 /* Identify and strip trailing zeros after the decimal point. */
4945 if ((trail
== 0) && (decflg
!= 0))
4948 while (ISDIGIT (*sp
) || (base
== 16 && ISXDIGIT (*sp
)))
4950 /* Check for syntax error */
4952 if ((base
!= 10 || ((c
!= 'e') && (c
!= 'E')))
4953 && (base
!= 16 || ((c
!= 'p') && (c
!= 'P')))
4955 && (c
!= '\n') && (c
!= '\r') && (c
!= ' ')
4957 goto unexpected_char_error
;
4966 /* If enough digits were given to more than fill up the yy register,
4967 continuing until overflow into the high guard word yy[2]
4968 guarantees that there will be a roundoff bit at the top
4969 of the low guard word after normalization. */
4976 nexp
+= 4; /* count digits after decimal point */
4978 eshup1 (yy
); /* multiply current number by 16 */
4986 nexp
+= 1; /* count digits after decimal point */
4988 eshup1 (yy
); /* multiply current number by 10 */
4994 /* Insert the current digit. */
4996 xt
[NI
- 2] = (UEMUSHORT
) k
;
5001 /* Mark any lost non-zero digit. */
5003 /* Count lost digits before the decimal point. */
5025 case '.': /* decimal point */
5027 goto unexpected_char_error
;
5033 goto unexpected_char_error
;
5038 goto unexpected_char_error
;
5051 unexpected_char_error
:
5055 mtherr ("asctoe", DOMAIN
);
5064 /* Exponent interpretation */
5066 /* 0.0eXXX is zero, regardless of XXX. Check for the 0.0. */
5067 for (k
= 0; k
< NI
; k
++)
5078 /* check for + or - */
5086 while (ISDIGIT (*s
))
5095 if ((exp
> MAXDECEXP
) && (base
== 10))
5099 yy
[E
] = 0x7fff; /* infinity */
5102 if ((exp
< MINDECEXP
) && (base
== 10))
5112 /* Base 16 hexadecimal floating constant. */
5113 if ((k
= enormlz (yy
)) > NBITS
)
5118 /* Adjust the exponent. NEXP is the number of hex digits,
5119 EXP is a power of 2. */
5120 lexp
= (EXONE
- 1 + NBITS
) - k
+ yy
[E
] + exp
- nexp
;
5130 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
5131 while ((nexp
> 0) && (yy
[2] == 0))
5143 if ((k
= enormlz (yy
)) > NBITS
)
5148 lexp
= (EXONE
- 1 + NBITS
) - k
;
5149 emdnorm (yy
, lost
, 0, lexp
, 64);
5152 /* Convert to external format:
5154 Multiply by 10**nexp. If precision is 64 bits,
5155 the maximum relative error incurred in forming 10**n
5156 for 0 <= n <= 324 is 8.2e-20, at 10**180.
5157 For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5158 For 0 >= n >= -999, it is -1.55e-19 at 10**-435. */
5173 /* Punt. Can't handle this without 2 divides. */
5174 emovi (etens
[0], tt
);
5187 emul (etens
[i
], xt
, xt
);
5191 while (exp
<= MAXP
);
5210 /* Round and convert directly to the destination type */
5212 lexp
-= EXONE
- 0x3ff;
5214 else if (oprec
== 24 || oprec
== 32)
5215 lexp
-= (EXONE
- 0x7f);
5218 else if (oprec
== 24 || oprec
== 56)
5219 lexp
-= EXONE
- (0x41 << 2);
5222 else if (oprec
== 24)
5223 lexp
-= dec_f
.adjustment
;
5224 else if (oprec
== 56)
5227 lexp
-= dec_g
.adjustment
;
5229 lexp
-= dec_d
.adjustment
;
5232 else if (oprec
== 24)
5233 lexp
-= EXONE
- 0177;
5238 emdnorm (yy
, lost
, 0, lexp
, 64);
5253 toibm (yy
, y
, DFmode
);
5258 toc4x (yy
, y
, HFmode
);
5271 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
5284 /* Return Y = largest integer not greater than X (truncated toward minus
5287 static const UEMUSHORT bmask
[] =
5310 const UEMUSHORT x
[];
5317 emov (x
, f
); /* leave in external format */
5318 expon
= (int) f
[NE
- 1];
5319 e
= (expon
& 0x7fff) - (EXONE
- 1);
5325 /* number of bits to clear out */
5337 /* clear the remaining bits */
5339 /* truncate negatives toward minus infinity */
5342 if ((UEMUSHORT
) expon
& (UEMUSHORT
) 0x8000)
5344 for (i
= 0; i
< NE
- 1; i
++)
5357 /* Return S and EXP such that S * 2^EXP = X and .5 <= S < 1.
5358 For example, 1.1 = 0.55 * 2^1. */
5362 const UEMUSHORT x
[];
5370 /* Handle denormalized numbers properly using long integer exponent. */
5371 li
= (EMULONG
) ((EMUSHORT
) xi
[1]);
5379 *exp
= (int) (li
- 0x3ffe);
5383 /* Return e type Y = X * 2^PWR2. */
5387 const UEMUSHORT x
[];
5399 emdnorm (xi
, i
, i
, li
, !ROUND_TOWARDS_ZERO
);
5405 /* C = remainder after dividing B by A, all e type values.
5406 Least significant integer quotient bits left in EQUOT. */
5410 const UEMUSHORT a
[], b
[];
5413 UEMUSHORT den
[NI
], num
[NI
];
5417 || (ecmp (a
, ezero
) == 0)
5425 if (ecmp (a
, ezero
) == 0)
5427 mtherr ("eremain", SING
);
5433 eiremain (den
, num
);
5434 /* Sign of remainder = sign of quotient */
5443 /* Return quotient of exploded e-types NUM / DEN in EQUOT,
5444 remainder in NUM. */
5448 UEMUSHORT den
[], num
[];
5454 ld
-= enormlz (den
);
5456 ln
-= enormlz (num
);
5460 if (ecmpm (den
, num
) <= 0)
5472 emdnorm (num
, 0, 0, ln
, 0);
5475 /* Report an error condition CODE encountered in function NAME.
5477 Mnemonic Value Significance
5479 DOMAIN 1 argument domain error
5480 SING 2 function singularity
5481 OVERFLOW 3 overflow range error
5482 UNDERFLOW 4 underflow range error
5483 TLOSS 5 total loss of precision
5484 PLOSS 6 partial loss of precision
5485 INVALID 7 NaN - producing operation
5486 EDOM 33 Unix domain error code
5487 ERANGE 34 Unix range error code
5489 The order of appearance of the following messages is bound to the
5490 error codes defined above. */
5500 /* The string passed by the calling program is supposed to be the
5501 name of the function in which the error occurred.
5502 The code argument selects which error message string will be printed. */
5504 if (strcmp (name
, "esub") == 0)
5505 name
= "subtraction";
5506 else if (strcmp (name
, "ediv") == 0)
5508 else if (strcmp (name
, "emul") == 0)
5509 name
= "multiplication";
5510 else if (strcmp (name
, "enormlz") == 0)
5511 name
= "normalization";
5512 else if (strcmp (name
, "etoasc") == 0)
5513 name
= "conversion to text";
5514 else if (strcmp (name
, "asctoe") == 0)
5516 else if (strcmp (name
, "eremain") == 0)
5518 else if (strcmp (name
, "esqrt") == 0)
5519 name
= "square root";
5524 case DOMAIN
: warning ("%s: argument domain error" , name
); break;
5525 case SING
: warning ("%s: function singularity" , name
); break;
5526 case OVERFLOW
: warning ("%s: overflow range error" , name
); break;
5527 case UNDERFLOW
: warning ("%s: underflow range error" , name
); break;
5528 case TLOSS
: warning ("%s: total loss of precision" , name
); break;
5529 case PLOSS
: warning ("%s: partial loss of precision", name
); break;
5530 case INVALID
: warning ("%s: NaN - producing operation", name
); break;
5535 /* Set global error message word */
5540 /* Convert DEC double precision D to e type E. */
5548 ieeetoe (d
, e
, &dec_g
);
5550 ieeetoe (d
, e
, &dec_d
);
5553 /* Convert e type X to DEC double precision D. */
5563 const struct ieee_format
*fmt
;
5571 /* Adjust exponent for offsets. */
5572 exp
= (EMULONG
) xi
[E
] - fmt
->adjustment
;
5573 /* Round off to nearest or even. */
5575 rndprc
= fmt
->precision
;
5576 emdnorm (xi
, 0, 0, exp
, !ROUND_TOWARDS_ZERO
);
5581 /* Convert exploded e-type X, that has already been rounded to
5582 56-bit precision, to DEC format double Y. */
5589 toieee (x
, y
, &dec_g
);
5591 toieee (x
, y
, &dec_d
);
5596 /* Convert IBM single/double precision to e type. */
5602 enum machine_mode mode
;
5607 ecleaz (y
); /* start with a zero */
5608 p
= y
; /* point to our number */
5609 r
= *d
; /* get IBM exponent word */
5610 if (*d
& (unsigned int) 0x8000)
5611 *p
= 0xffff; /* fill in our sign */
5612 ++p
; /* bump pointer to our exponent word */
5613 r
&= 0x7f00; /* strip the sign bit */
5614 r
>>= 6; /* shift exponent word down 6 bits */
5615 /* in fact shift by 8 right and 2 left */
5616 r
+= EXONE
- (0x41 << 2); /* subtract IBM exponent offset */
5617 /* add our e type exponent offset */
5618 *p
++ = r
; /* to form our exponent */
5620 *p
++ = *d
++ & 0xff; /* now do the high order mantissa */
5621 /* strip off the IBM exponent and sign bits */
5622 if (mode
!= SFmode
) /* there are only 2 words in SFmode */
5624 *p
++ = *d
++; /* fill in the rest of our mantissa */
5629 if (y
[M
] == 0 && y
[M
+1] == 0 && y
[M
+2] == 0 && y
[M
+3] == 0)
5632 y
[E
] -= 5 + enormlz (y
); /* now normalise the mantissa */
5633 /* handle change in RADIX */
5639 /* Convert e type to IBM single/double precision. */
5645 enum machine_mode mode
;
5652 exp
= (EMULONG
) xi
[E
] - (EXONE
- (0x41 << 2)); /* adjust exponent for offsets */
5653 /* round off to nearest or even */
5656 emdnorm (xi
, 0, 0, exp
, !ROUND_TOWARDS_ZERO
);
5658 toibm (xi
, d
, mode
);
5664 enum machine_mode mode
;
5717 /* Convert C4X single/double precision to e type. */
5723 enum machine_mode mode
;
5741 /* Short-circuit the zero case. */
5742 if ((dn
[0] == 0x8000)
5743 && (dn
[1] == 0x0000)
5744 && ((mode
== QFmode
) || ((dn
[2] == 0x0000) && (dn
[3] == 0x0000))))
5755 ecleaz (y
); /* start with a zero */
5756 r
= dn
[0]; /* get sign/exponent part */
5757 if (r
& (unsigned int) 0x0080)
5759 y
[0] = 0xffff; /* fill in our sign */
5765 r
>>= 8; /* Shift exponent word down 8 bits. */
5766 if (r
& 0x80) /* Make the exponent negative if it is. */
5767 r
= r
| (~0 & ~0xff);
5771 /* Now do the high order mantissa. We don't "or" on the high bit
5772 because it is 2 (not 1) and is handled a little differently
5774 y
[M
] = dn
[0] & 0x7f;
5777 if (mode
!= QFmode
) /* There are only 2 words in QFmode. */
5779 y
[M
+2] = dn
[2]; /* Fill in the rest of our mantissa. */
5787 /* Now do the two's complement on the data. */
5789 carry
= 1; /* Initially add 1 for the two's complement. */
5790 for (i
=size
+ M
; i
> M
; i
--)
5792 if (carry
&& (y
[i
] == 0x0000))
5793 /* We overflowed into the next word, carry is the same. */
5794 y
[i
] = carry
? 0x0000 : 0xffff;
5797 /* No overflow, just invert and add carry. */
5798 y
[i
] = ((~y
[i
]) + carry
) & 0xffff;
5813 /* Add our e type exponent offset to form our exponent. */
5817 /* Now do the high order mantissa strip off the exponent and sign
5818 bits and add the high 1 bit. */
5819 y
[M
] = (dn
[0] & 0x7f) | 0x80;
5822 if (mode
!= QFmode
) /* There are only 2 words in QFmode. */
5824 y
[M
+2] = dn
[2]; /* Fill in the rest of our mantissa. */
5834 /* Convert e type to C4X single/double precision. */
5840 enum machine_mode mode
;
5848 /* Adjust exponent for offsets. */
5849 exp
= (EMULONG
) xi
[E
] - (EXONE
- 0x7f);
5851 /* Round off to nearest or even. */
5853 rndprc
= mode
== QFmode
? 24 : 32;
5854 emdnorm (xi
, 0, 0, exp
, !ROUND_TOWARDS_ZERO
);
5856 toc4x (xi
, d
, mode
);
5862 enum machine_mode mode
;
5868 /* Short-circuit the zero case */
5869 if ((x
[0] == 0) /* Zero exponent and sign */
5871 && (x
[M
] == 0) /* The rest is for zero mantissa */
5873 /* Only check for double if necessary */
5874 && ((mode
== QFmode
) || ((x
[M
+2] == 0) && (x
[M
+3] == 0))))
5876 /* We have a zero. Put it into the output and return. */
5889 /* Negative number require a two's complement conversion of the
5895 i
= ((int) x
[1]) - 0x7f;
5897 /* Now add 1 to the inverted data to do the two's complement. */
5906 x
[v
] = carry
? 0x0000 : 0xffff;
5909 x
[v
] = ((~x
[v
]) + carry
) & 0xffff;
5915 /* The following is a special case. The C4X negative float requires
5916 a zero in the high bit (because the format is (2 - x) x 2^m), so
5917 if a one is in that bit, we have to shift left one to get rid
5918 of it. This only occurs if the number is -1 x 2^m. */
5919 if (x
[M
+1] & 0x8000)
5921 /* This is the case of -1 x 2^m, we have to rid ourselves of the
5922 high sign bit and shift the exponent. */
5928 i
= ((int) x
[1]) - 0x7f;
5930 if ((i
< -128) || (i
> 127))
5938 y
[3] = (y
[1] << 8) | ((y
[2] >> 8) & 0xff);
5939 y
[2] = (y
[0] << 8) | ((y
[1] >> 8) & 0xff);
5947 y
[0] |= ((i
& 0xff) << 8);
5951 y
[0] |= x
[M
] & 0x7f;
5957 y
[3] = (y
[1] << 8) | ((y
[2] >> 8) & 0xff);
5958 y
[2] = (y
[0] << 8) | ((y
[1] >> 8) & 0xff);
5963 /* Output a binary NaN bit pattern in the target machine's format. */
5965 /* If special NaN bit patterns are required, define them in tm.h
5966 as arrays of unsigned 16-bit shorts. Otherwise, use the default
5971 #if defined (IEEE) && (INTEL_EXTENDED_IEEE_FORMAT == 0)
5972 static const UEMUSHORT TFbignan
[8] =
5973 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5974 static const UEMUSHORT TFlittlenan
[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
5982 static const UEMUSHORT XFbignan
[6] =
5983 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5984 static const UEMUSHORT XFlittlenan
[6] = {0, 0, 0, 0xc000, 0xffff, 0};
5992 static const UEMUSHORT DFbignan
[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
5993 static const UEMUSHORT DFlittlenan
[4] = {0, 0, 0, 0xfff8};
6001 static const UEMUSHORT SFbignan
[2] = {0x7fff, 0xffff};
6002 static const UEMUSHORT SFlittlenan
[2] = {0, 0xffc0};
6009 make_nan (nan
, sign
, mode
)
6012 enum machine_mode mode
;
6018 size
= GET_MODE_BITSIZE (mode
);
6019 if (LARGEST_EXPONENT_IS_NORMAL (size
))
6021 warning ("%d-bit floats cannot hold NaNs", size
);
6022 saturate (nan
, sign
, size
, 0);
6027 /* Possibly the `reserved operand' patterns on a VAX can be
6028 used like NaN's, but probably not in the same way as IEEE. */
6029 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
6031 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6033 if (REAL_WORDS_BIG_ENDIAN
)
6043 if (REAL_WORDS_BIG_ENDIAN
)
6051 if (REAL_WORDS_BIG_ENDIAN
)
6060 if (REAL_WORDS_BIG_ENDIAN
)
6070 if (REAL_WORDS_BIG_ENDIAN
)
6071 *nan
++ = (sign
<< 15) | (*p
++ & 0x7fff);
6074 if (! REAL_WORDS_BIG_ENDIAN
)
6075 *nan
= (sign
<< 15) | (*p
& 0x7fff);
6080 /* Create a saturation value for a SIZE-bit float, assuming that
6081 LARGEST_EXPONENT_IS_NORMAL (SIZE).
6083 If SIGN is true, fill X with the most negative value, otherwise fill
6084 it with the most positive value. WARN is true if the function should
6085 warn about overflow. */
6088 saturate (x
, sign
, size
, warn
)
6090 int sign
, size
, warn
;
6094 if (warn
&& extra_warnings
)
6095 warning ("value exceeds the range of a %d-bit float", size
);
6097 /* Create the most negative value. */
6098 for (i
= 0; i
< size
/ EMUSHORT_SIZE
; i
++)
6101 /* Make it positive, if necessary. */
6103 x
[REAL_WORDS_BIG_ENDIAN
? 0 : i
- 1] = 0x7fff;
6107 /* This is the inverse of the function `etarsingle' invoked by
6108 REAL_VALUE_TO_TARGET_SINGLE. */
6111 ereal_unto_float (f
)
6118 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6119 This is the inverse operation to what the function `endian' does. */
6120 if (REAL_WORDS_BIG_ENDIAN
)
6122 s
[0] = (UEMUSHORT
) (f
>> 16);
6123 s
[1] = (UEMUSHORT
) f
;
6127 s
[0] = (UEMUSHORT
) f
;
6128 s
[1] = (UEMUSHORT
) (f
>> 16);
6130 /* Convert and promote the target float to E-type. */
6132 /* Output E-type to REAL_VALUE_TYPE. */
6138 /* This is the inverse of the function `etardouble' invoked by
6139 REAL_VALUE_TO_TARGET_DOUBLE. */
6142 ereal_unto_double (d
)
6149 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
6150 if (REAL_WORDS_BIG_ENDIAN
)
6152 s
[0] = (UEMUSHORT
) (d
[0] >> 16);
6153 s
[1] = (UEMUSHORT
) d
[0];
6154 s
[2] = (UEMUSHORT
) (d
[1] >> 16);
6155 s
[3] = (UEMUSHORT
) d
[1];
6159 /* Target float words are little-endian. */
6160 s
[0] = (UEMUSHORT
) d
[0];
6161 s
[1] = (UEMUSHORT
) (d
[0] >> 16);
6162 s
[2] = (UEMUSHORT
) d
[1];
6163 s
[3] = (UEMUSHORT
) (d
[1] >> 16);
6165 /* Convert target double to E-type. */
6167 /* Output E-type to REAL_VALUE_TYPE. */
6173 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
6174 This is somewhat like ereal_unto_float, but the input types
6175 for these are different. */
6178 ereal_from_float (f
)
6185 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6186 This is the inverse operation to what the function `endian' does. */
6187 if (REAL_WORDS_BIG_ENDIAN
)
6189 s
[0] = (UEMUSHORT
) (f
>> 16);
6190 s
[1] = (UEMUSHORT
) f
;
6194 s
[0] = (UEMUSHORT
) f
;
6195 s
[1] = (UEMUSHORT
) (f
>> 16);
6197 /* Convert and promote the target float to E-type. */
6199 /* Output E-type to REAL_VALUE_TYPE. */
6205 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
6206 This is somewhat like ereal_unto_double, but the input types
6207 for these are different.
6209 The DFmode is stored as an array of HOST_WIDE_INT in the target's
6210 data format, with no holes in the bit packing. The first element
6211 of the input array holds the bits that would come first in the
6212 target computer's memory. */
6215 ereal_from_double (d
)
6222 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
6223 if (REAL_WORDS_BIG_ENDIAN
)
6225 #if HOST_BITS_PER_WIDE_INT == 32
6226 s
[0] = (UEMUSHORT
) (d
[0] >> 16);
6227 s
[1] = (UEMUSHORT
) d
[0];
6228 s
[2] = (UEMUSHORT
) (d
[1] >> 16);
6229 s
[3] = (UEMUSHORT
) d
[1];
6231 /* In this case the entire target double is contained in the
6232 first array element. The second element of the input is
6234 s
[0] = (UEMUSHORT
) (d
[0] >> 48);
6235 s
[1] = (UEMUSHORT
) (d
[0] >> 32);
6236 s
[2] = (UEMUSHORT
) (d
[0] >> 16);
6237 s
[3] = (UEMUSHORT
) d
[0];
6242 /* Target float words are little-endian. */
6243 s
[0] = (UEMUSHORT
) d
[0];
6244 s
[1] = (UEMUSHORT
) (d
[0] >> 16);
6245 #if HOST_BITS_PER_WIDE_INT == 32
6246 s
[2] = (UEMUSHORT
) d
[1];
6247 s
[3] = (UEMUSHORT
) (d
[1] >> 16);
6249 s
[2] = (UEMUSHORT
) (d
[0] >> 32);
6250 s
[3] = (UEMUSHORT
) (d
[0] >> 48);
6253 /* Convert target double to E-type. */
6255 /* Output E-type to REAL_VALUE_TYPE. */
6262 /* Convert target computer unsigned 64-bit integer to e-type.
6263 The endian-ness of DImode follows the convention for integers,
6264 so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN. */
6268 const UEMUSHORT
*di
; /* Address of the 64-bit int. */
6275 if (WORDS_BIG_ENDIAN
)
6277 for (k
= M
; k
< M
+ 4; k
++)
6282 for (k
= M
+ 3; k
>= M
; k
--)
6285 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
6286 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
6287 ecleaz (yi
); /* it was zero */
6289 yi
[E
] -= (UEMUSHORT
) k
;/* subtract shift count from exponent */
6293 /* Convert target computer signed 64-bit integer to e-type. */
6297 const UEMUSHORT
*di
; /* Address of the 64-bit int. */
6300 unsigned EMULONG acc
;
6306 if (WORDS_BIG_ENDIAN
)
6308 for (k
= M
; k
< M
+ 4; k
++)
6313 for (k
= M
+ 3; k
>= M
; k
--)
6316 /* Take absolute value */
6322 for (k
= M
+ 3; k
>= M
; k
--)
6324 acc
= (unsigned EMULONG
) (~yi
[k
] & 0xffff) + carry
;
6331 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
6332 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
6333 ecleaz (yi
); /* it was zero */
6335 yi
[E
] -= (UEMUSHORT
) k
;/* subtract shift count from exponent */
6342 /* Convert e-type to unsigned 64-bit int. */
6358 k
= (int) xi
[E
] - (EXONE
- 1);
6361 for (j
= 0; j
< 4; j
++)
6367 for (j
= 0; j
< 4; j
++)
6370 warning ("overflow on truncation to integer");
6375 /* Shift more than 16 bits: first shift up k-16 mod 16,
6376 then shift up by 16's. */
6377 j
= k
- ((k
>> 4) << 4);
6381 if (WORDS_BIG_ENDIAN
)
6392 if (WORDS_BIG_ENDIAN
)
6397 while ((k
-= 16) > 0);
6401 /* shift not more than 16 bits */
6406 if (WORDS_BIG_ENDIAN
)
6425 /* Convert e-type to signed 64-bit int. */
6432 unsigned EMULONG acc
;
6439 k
= (int) xi
[E
] - (EXONE
- 1);
6442 for (j
= 0; j
< 4; j
++)
6448 for (j
= 0; j
< 4; j
++)
6451 warning ("overflow on truncation to integer");
6457 /* Shift more than 16 bits: first shift up k-16 mod 16,
6458 then shift up by 16's. */
6459 j
= k
- ((k
>> 4) << 4);
6463 if (WORDS_BIG_ENDIAN
)
6474 if (WORDS_BIG_ENDIAN
)
6479 while ((k
-= 16) > 0);
6483 /* shift not more than 16 bits */
6486 if (WORDS_BIG_ENDIAN
)
6502 /* Negate if negative */
6506 if (WORDS_BIG_ENDIAN
)
6508 for (k
= 0; k
< 4; k
++)
6510 acc
= (unsigned EMULONG
) (~(*isave
) & 0xffff) + carry
;
6511 if (WORDS_BIG_ENDIAN
)
6523 /* Longhand square root routine. */
6526 static int esqinited
= 0;
6527 static unsigned short sqrndbit
[NI
];
6534 UEMUSHORT temp
[NI
], num
[NI
], sq
[NI
], xx
[NI
];
6536 int i
, j
, k
, n
, nlups
;
6541 sqrndbit
[NI
- 2] = 1;
6544 /* Check for arg <= 0 */
6545 i
= ecmp (x
, ezero
);
6550 mtherr ("esqrt", DOMAIN
);
6566 /* Bring in the arg and renormalize if it is denormal. */
6568 m
= (EMULONG
) xx
[1]; /* local long word exponent */
6572 /* Divide exponent by 2 */
6574 exp
= (unsigned short) ((m
/ 2) + 0x3ffe);
6576 /* Adjust if exponent odd */
6586 n
= 8; /* get 8 bits of result per inner loop */
6592 /* bring in next word of arg */
6594 num
[NI
- 1] = xx
[j
+ 3];
6595 /* Do additional bit on last outer loop, for roundoff. */
6598 for (i
= 0; i
< n
; i
++)
6600 /* Next 2 bits of arg */
6603 /* Shift up answer */
6605 /* Make trial divisor */
6606 for (k
= 0; k
< NI
; k
++)
6609 eaddm (sqrndbit
, temp
);
6610 /* Subtract and insert answer bit if it goes in */
6611 if (ecmpm (temp
, num
) <= 0)
6621 /* Adjust for extra, roundoff loop done. */
6622 exp
+= (NBITS
- 1) - rndprc
;
6624 /* Sticky bit = 1 if the remainder is nonzero. */
6626 for (i
= 3; i
< NI
; i
++)
6629 /* Renormalize and round off. */
6630 emdnorm (sq
, k
, 0, exp
, !ROUND_TOWARDS_ZERO
);
6635 /* Return the binary precision of the significand for a given
6636 floating point mode. The mode can hold an integer value
6637 that many bits wide, without losing any bits. */
6640 significand_size (mode
)
6641 enum machine_mode mode
;
6644 /* Don't test the modes, but their sizes, lest this
6645 code won't work for BITS_PER_UNIT != 8 . */
6647 switch (GET_MODE_BITSIZE (mode
))
6668 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)