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. */
324 /* IEEE float (24 bits). */
325 static const struct ieee_format ieee_24
=
334 /* IEEE double (53 bits). */
335 static const struct ieee_format ieee_53
=
344 /* IEEE extended double (64 bits). */
345 static const struct ieee_format ieee_64
=
354 /* IEEE long double (113 bits). */
355 static const struct ieee_format ieee_113
=
364 /* DEC F float (24 bits). */
365 static const struct ieee_format dec_f
=
374 /* DEC D float (56 bits). */
375 static const struct ieee_format dec_d
=
384 /* DEC G float (53 bits). */
385 static const struct ieee_format dec_g
=
394 /* DEC H float (113 bits). (not yet used) */
395 static const struct ieee_format dec_h
=
404 extern int extra_warnings
;
405 extern const UEMUSHORT ezero
[NE
], ehalf
[NE
], eone
[NE
], etwo
[NE
];
406 extern const UEMUSHORT elog2
[NE
], esqrt2
[NE
];
408 static void endian
PARAMS ((const UEMUSHORT
*, long *,
410 static void eclear
PARAMS ((UEMUSHORT
*));
411 static void emov
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
413 static void eabs
PARAMS ((UEMUSHORT
*));
415 static void eneg
PARAMS ((UEMUSHORT
*));
416 static int eisneg
PARAMS ((const UEMUSHORT
*));
417 static int eisinf
PARAMS ((const UEMUSHORT
*));
418 static int eisnan
PARAMS ((const UEMUSHORT
*));
419 static void einfin
PARAMS ((UEMUSHORT
*));
421 static void enan
PARAMS ((UEMUSHORT
*, int));
422 static void einan
PARAMS ((UEMUSHORT
*));
423 static int eiisnan
PARAMS ((const UEMUSHORT
*));
424 static void make_nan
PARAMS ((UEMUSHORT
*, int, enum machine_mode
));
426 static int eiisneg
PARAMS ((const UEMUSHORT
*));
427 static void saturate
PARAMS ((UEMUSHORT
*, int, int, int));
428 static void emovi
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
429 static void emovo
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
430 static void ecleaz
PARAMS ((UEMUSHORT
*));
431 static void ecleazs
PARAMS ((UEMUSHORT
*));
432 static void emovz
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
434 static void eiinfin
PARAMS ((UEMUSHORT
*));
437 static int eiisinf
PARAMS ((const UEMUSHORT
*));
439 static int ecmpm
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*));
440 static void eshdn1
PARAMS ((UEMUSHORT
*));
441 static void eshup1
PARAMS ((UEMUSHORT
*));
442 static void eshdn8
PARAMS ((UEMUSHORT
*));
443 static void eshup8
PARAMS ((UEMUSHORT
*));
444 static void eshup6
PARAMS ((UEMUSHORT
*));
445 static void eshdn6
PARAMS ((UEMUSHORT
*));
446 static void eaddm
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));\f
447 static void esubm
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
448 static void m16m
PARAMS ((unsigned int, const UEMUSHORT
*, UEMUSHORT
*));
449 static int edivm
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
450 static int emulm
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
451 static void emdnorm
PARAMS ((UEMUSHORT
*, int, int, EMULONG
, int));
452 static void esub
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
454 static void eadd
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
456 static void eadd1
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
458 static void ediv
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
460 static void emul
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
462 static void e53toe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
463 static void e64toe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
464 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
465 static void e113toe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
467 static void e24toe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
468 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
469 static void etoe113
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
470 static void toe113
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
472 static void etoe64
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
473 static void toe64
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
474 static void etoe53
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
475 static void toe53
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
476 static void etoe24
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
477 static void toe24
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
478 static void ieeetoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
479 const struct ieee_format
*));
480 static void etoieee
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
481 const struct ieee_format
*));
482 static void toieee
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
483 const struct ieee_format
*));
484 static int ecmp
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*));
486 static void eround
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
488 static void ltoe
PARAMS ((const HOST_WIDE_INT
*, UEMUSHORT
*));
489 static void ultoe
PARAMS ((const unsigned HOST_WIDE_INT
*, UEMUSHORT
*));
490 static void eifrac
PARAMS ((const UEMUSHORT
*, HOST_WIDE_INT
*,
492 static void euifrac
PARAMS ((const UEMUSHORT
*, unsigned HOST_WIDE_INT
*,
494 static int eshift
PARAMS ((UEMUSHORT
*, int));
495 static int enormlz
PARAMS ((UEMUSHORT
*));
497 static void e24toasc
PARAMS ((const UEMUSHORT
*, char *, int));
498 static void e53toasc
PARAMS ((const UEMUSHORT
*, char *, int));
499 static void e64toasc
PARAMS ((const UEMUSHORT
*, char *, int));
500 static void e113toasc
PARAMS ((const UEMUSHORT
*, char *, int));
502 static void etoasc
PARAMS ((const UEMUSHORT
*, char *, int));
503 static void asctoe24
PARAMS ((const char *, UEMUSHORT
*));
504 static void asctoe53
PARAMS ((const char *, UEMUSHORT
*));
505 static void asctoe64
PARAMS ((const char *, UEMUSHORT
*));
506 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
507 static void asctoe113
PARAMS ((const char *, UEMUSHORT
*));
509 static void asctoe
PARAMS ((const char *, UEMUSHORT
*));
510 static void asctoeg
PARAMS ((const char *, UEMUSHORT
*, int));
511 static void efloor
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
513 static void efrexp
PARAMS ((const UEMUSHORT
*, int *,
516 static void eldexp
PARAMS ((const UEMUSHORT
*, int, UEMUSHORT
*));
518 static void eremain
PARAMS ((const UEMUSHORT
*, const UEMUSHORT
*,
521 static void eiremain
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
522 static void mtherr
PARAMS ((const char *, int));
524 static void dectoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
525 static void etodec
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
526 static void todec
PARAMS ((UEMUSHORT
*, UEMUSHORT
*));
529 static void ibmtoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
531 static void etoibm
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
533 static void toibm
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
537 static void c4xtoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
539 static void etoc4x
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*,
541 static void toc4x
PARAMS ((UEMUSHORT
*, UEMUSHORT
*,
545 static void uditoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
546 static void ditoe
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
547 static void etoudi
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
548 static void etodi
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
549 static void esqrt
PARAMS ((const UEMUSHORT
*, UEMUSHORT
*));
552 /* Copy 32-bit numbers obtained from array containing 16-bit numbers,
553 swapping ends if required, into output array of longs. The
554 result is normally passed to fprintf by the ASM_OUTPUT_ macros. */
560 enum machine_mode mode
;
564 if (REAL_WORDS_BIG_ENDIAN
&& !VAX_HALFWORD_ORDER
)
569 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
570 /* Swap halfwords in the fourth long. */
571 th
= (unsigned long) e
[6] & 0xffff;
572 t
= (unsigned long) e
[7] & 0xffff;
581 /* Swap halfwords in the third long. */
582 th
= (unsigned long) e
[4] & 0xffff;
583 t
= (unsigned long) e
[5] & 0xffff;
589 /* Swap halfwords in the second word. */
590 th
= (unsigned long) e
[2] & 0xffff;
591 t
= (unsigned long) e
[3] & 0xffff;
598 /* Swap halfwords in the first word. */
599 th
= (unsigned long) e
[0] & 0xffff;
600 t
= (unsigned long) e
[1] & 0xffff;
611 /* Pack the output array without swapping. */
616 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
617 /* Pack the fourth long. */
618 th
= (unsigned long) e
[7] & 0xffff;
619 t
= (unsigned long) e
[6] & 0xffff;
628 /* Pack the third long.
629 Each element of the input REAL_VALUE_TYPE array has 16 useful bits
631 th
= (unsigned long) e
[5] & 0xffff;
632 t
= (unsigned long) e
[4] & 0xffff;
638 /* Pack the second long */
639 th
= (unsigned long) e
[3] & 0xffff;
640 t
= (unsigned long) e
[2] & 0xffff;
647 /* Pack the first long */
648 th
= (unsigned long) e
[1] & 0xffff;
649 t
= (unsigned long) e
[0] & 0xffff;
661 /* This is the implementation of the REAL_ARITHMETIC macro. */
664 earith (value
, icode
, r1
, r2
)
665 REAL_VALUE_TYPE
*value
;
670 UEMUSHORT d1
[NE
], d2
[NE
], v
[NE
];
676 /* Return NaN input back to the caller. */
679 PUT_REAL (d1
, value
);
684 PUT_REAL (d2
, value
);
688 code
= (enum tree_code
) icode
;
696 esub (d2
, d1
, v
); /* d1 - d2 */
705 if (ecmp (d2
, ezero
) == 0)
708 ediv (d2
, d1
, v
); /* d1/d2 */
711 case MIN_EXPR
: /* min (d1,d2) */
712 if (ecmp (d1
, d2
) < 0)
718 case MAX_EXPR
: /* max (d1,d2) */
719 if (ecmp (d1
, d2
) > 0)
732 /* Truncate REAL_VALUE_TYPE toward zero to signed HOST_WIDE_INT.
733 implements REAL_VALUE_RNDZINT (x) (etrunci (x)). */
739 UEMUSHORT f
[NE
], g
[NE
];
755 /* Truncate REAL_VALUE_TYPE toward zero to unsigned HOST_WIDE_INT;
756 implements REAL_VALUE_UNSIGNED_RNDZINT (x) (etruncui (x)). */
762 UEMUSHORT f
[NE
], g
[NE
];
764 unsigned HOST_WIDE_INT l
;
778 /* This is the REAL_VALUE_ATOF function. It converts a decimal or hexadecimal
779 string to binary, rounding off as indicated by the machine_mode argument.
780 Then it promotes the rounded value to REAL_VALUE_TYPE. */
787 UEMUSHORT tem
[NE
], e
[NE
];
813 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
833 /* Expansion of REAL_NEGATE. */
849 /* Round real toward zero to HOST_WIDE_INT;
850 implements REAL_VALUE_FIX (x). */
856 UEMUSHORT f
[NE
], g
[NE
];
863 warning ("conversion from NaN to int");
871 /* Round real toward zero to unsigned HOST_WIDE_INT
872 implements REAL_VALUE_UNSIGNED_FIX (x).
873 Negative input returns zero. */
875 unsigned HOST_WIDE_INT
879 UEMUSHORT f
[NE
], g
[NE
];
880 unsigned HOST_WIDE_INT l
;
886 warning ("conversion from NaN to unsigned int");
895 /* REAL_VALUE_FROM_INT macro. */
898 ereal_from_int (d
, i
, j
, mode
)
901 enum machine_mode mode
;
903 UEMUSHORT df
[NE
], dg
[NE
];
904 HOST_WIDE_INT low
, high
;
907 if (GET_MODE_CLASS (mode
) != MODE_FLOAT
)
914 /* complement and add 1 */
921 eldexp (eone
, HOST_BITS_PER_WIDE_INT
, df
);
922 ultoe ((unsigned HOST_WIDE_INT
*) &high
, dg
);
924 ultoe ((unsigned HOST_WIDE_INT
*) &low
, df
);
929 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
930 Avoid double-rounding errors later by rounding off now from the
931 extra-wide internal format to the requested precision. */
932 switch (GET_MODE_BITSIZE (mode
))
950 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
967 /* REAL_VALUE_FROM_UNSIGNED_INT macro. */
970 ereal_from_uint (d
, i
, j
, mode
)
972 unsigned HOST_WIDE_INT i
, j
;
973 enum machine_mode mode
;
975 UEMUSHORT df
[NE
], dg
[NE
];
976 unsigned HOST_WIDE_INT low
, high
;
978 if (GET_MODE_CLASS (mode
) != MODE_FLOAT
)
982 eldexp (eone
, HOST_BITS_PER_WIDE_INT
, df
);
988 /* A REAL_VALUE_TYPE may not be wide enough to hold the two HOST_WIDE_INTS.
989 Avoid double-rounding errors later by rounding off now from the
990 extra-wide internal format to the requested precision. */
991 switch (GET_MODE_BITSIZE (mode
))
1009 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
1026 /* REAL_VALUE_TO_INT macro. */
1029 ereal_to_int (low
, high
, rr
)
1030 HOST_WIDE_INT
*low
, *high
;
1033 UEMUSHORT d
[NE
], df
[NE
], dg
[NE
], dh
[NE
];
1040 warning ("conversion from NaN to int");
1046 /* convert positive value */
1053 eldexp (eone
, HOST_BITS_PER_WIDE_INT
, df
);
1054 ediv (df
, d
, dg
); /* dg = d / 2^32 is the high word */
1055 euifrac (dg
, (unsigned HOST_WIDE_INT
*) high
, dh
);
1056 emul (df
, dh
, dg
); /* fractional part is the low word */
1057 euifrac (dg
, (unsigned HOST_WIDE_INT
*) low
, dh
);
1060 /* complement and add 1 */
1070 /* REAL_VALUE_LDEXP macro. */
1077 UEMUSHORT e
[NE
], y
[NE
];
1090 /* Check for infinity in a REAL_VALUE_TYPE. */
1094 REAL_VALUE_TYPE x ATTRIBUTE_UNUSED
;
1100 return (eisinf (e
));
1106 /* Check whether a REAL_VALUE_TYPE item is a NaN. */
1110 REAL_VALUE_TYPE x ATTRIBUTE_UNUSED
;
1116 return (eisnan (e
));
1123 /* Check for a negative REAL_VALUE_TYPE number.
1124 This just checks the sign bit, so that -0 counts as negative. */
1130 return ereal_isneg (x
);
1133 /* Expansion of REAL_VALUE_TRUNCATE.
1134 The result is in floating point, rounded to nearest or even. */
1137 real_value_truncate (mode
, arg
)
1138 enum machine_mode mode
;
1139 REAL_VALUE_TYPE arg
;
1141 UEMUSHORT e
[NE
], t
[NE
];
1153 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
1190 /* If an unsupported type was requested, presume that
1191 the machine files know something useful to do with
1192 the unmodified value. */
1201 /* Return true if ARG can be represented exactly in MODE. */
1204 exact_real_truncate (mode
, arg
)
1205 enum machine_mode mode
;
1206 REAL_VALUE_TYPE
*arg
;
1208 REAL_VALUE_TYPE trunc
;
1210 if (target_isnan (*arg
))
1213 trunc
= real_value_truncate (mode
, *arg
);
1214 return ereal_cmp (*arg
, trunc
) == 0;
1217 /* Try to change R into its exact multiplicative inverse in machine mode
1218 MODE. Return nonzero function value if successful. */
1221 exact_real_inverse (mode
, r
)
1222 enum machine_mode mode
;
1225 UEMUSHORT e
[NE
], einv
[NE
];
1226 REAL_VALUE_TYPE rinv
;
1231 /* Test for input in range. Don't transform IEEE special values. */
1232 if (eisinf (e
) || eisnan (e
) || (ecmp (e
, ezero
) == 0))
1235 /* Test for a power of 2: all significand bits zero except the MSB.
1236 We are assuming the target has binary (or hex) arithmetic. */
1237 if (e
[NE
- 2] != 0x8000)
1240 for (i
= 0; i
< NE
- 2; i
++)
1246 /* Compute the inverse and truncate it to the required mode. */
1247 ediv (e
, eone
, einv
);
1248 PUT_REAL (einv
, &rinv
);
1249 rinv
= real_value_truncate (mode
, rinv
);
1251 #ifdef CHECK_FLOAT_VALUE
1252 /* This check is not redundant. It may, for example, flush
1253 a supposedly IEEE denormal value to zero. */
1255 if (CHECK_FLOAT_VALUE (mode
, rinv
, i
))
1258 GET_REAL (&rinv
, einv
);
1260 /* Check the bits again, because the truncation might have
1261 generated an arbitrary saturation value on overflow. */
1262 if (einv
[NE
- 2] != 0x8000)
1265 for (i
= 0; i
< NE
- 2; i
++)
1271 /* Fail if the computed inverse is out of range. */
1272 if (eisinf (einv
) || eisnan (einv
) || (ecmp (einv
, ezero
) == 0))
1275 /* Output the reciprocal and return success flag. */
1280 /* Used for debugging--print the value of R in human-readable format
1289 REAL_VALUE_TO_DECIMAL (r
, "%.20g", dstr
);
1290 fprintf (stderr
, "%s", dstr
);
1294 /* The following routines convert REAL_VALUE_TYPE to the various floating
1295 point formats that are meaningful to supported computers.
1297 The results are returned in 32-bit pieces, each piece stored in a `long'.
1298 This is so they can be printed by statements like
1300 fprintf (file, "%lx, %lx", L[0], L[1]);
1302 that will work on both narrow- and wide-word host computers. */
1304 /* Convert R to a 128-bit long double precision value. The output array L
1305 contains four 32-bit pieces of the result, in the order they would appear
1316 #if INTEL_EXTENDED_IEEE_FORMAT == 0
1321 endian (e
, l
, TFmode
);
1324 /* Convert R to a double extended precision value. The output array L
1325 contains three 32-bit pieces of the result, in the order they would
1326 appear in memory. */
1337 endian (e
, l
, XFmode
);
1340 /* Convert R to a double precision value. The output array L contains two
1341 32-bit pieces of the result, in the order they would appear in memory. */
1352 endian (e
, l
, DFmode
);
1355 /* Convert R to a single precision float value stored in the least-significant
1356 bits of a `long'. */
1367 endian (e
, &l
, SFmode
);
1371 /* Convert X to a decimal ASCII string S for output to an assembly
1372 language file. Note, there is no standard way to spell infinity or
1373 a NaN, so these values may require special treatment in the tm.h
1377 ereal_to_decimal (x
, s
)
1387 /* Compare X and Y. Return 1 if X > Y, 0 if X == Y, -1 if X < Y,
1388 or -2 if either is a NaN. */
1392 REAL_VALUE_TYPE x
, y
;
1394 UEMUSHORT ex
[NE
], ey
[NE
];
1398 return (ecmp (ex
, ey
));
1401 /* Return 1 if the sign bit of X is set, else return 0. */
1410 return (eisneg (ex
));
1415 Extended precision IEEE binary floating point arithmetic routines
1417 Numbers are stored in C language as arrays of 16-bit unsigned
1418 short integers. The arguments of the routines are pointers to
1421 External e type data structure, similar to Intel 8087 chip
1422 temporary real format but possibly with a larger significand:
1424 NE-1 significand words (least significant word first,
1425 most significant bit is normally set)
1426 exponent (value = EXONE for 1.0,
1427 top bit is the sign)
1430 Internal exploded e-type data structure of a number (a "word" is 16 bits):
1432 ei[0] sign word (0 for positive, 0xffff for negative)
1433 ei[1] biased exponent (value = EXONE for the number 1.0)
1434 ei[2] high guard word (always zero after normalization)
1436 to ei[NI-2] significand (NI-4 significand words,
1437 most significant word first,
1438 most significant bit is set)
1439 ei[NI-1] low guard word (0x8000 bit is rounding place)
1443 Routines for external format e-type numbers
1445 asctoe (string, e) ASCII string to extended double e type
1446 asctoe64 (string, &d) ASCII string to long double
1447 asctoe53 (string, &d) ASCII string to double
1448 asctoe24 (string, &f) ASCII string to single
1449 asctoeg (string, e, prec) ASCII string to specified precision
1450 e24toe (&f, e) IEEE single precision to e type
1451 e53toe (&d, e) IEEE double precision to e type
1452 e64toe (&d, e) IEEE long double precision to e type
1453 e113toe (&d, e) 128-bit long double precision to e type
1455 eabs (e) absolute value
1457 eadd (a, b, c) c = b + a
1459 ecmp (a, b) Returns 1 if a > b, 0 if a == b,
1460 -1 if a < b, -2 if either a or b is a NaN.
1461 ediv (a, b, c) c = b / a
1462 efloor (a, b) truncate to integer, toward -infinity
1463 efrexp (a, exp, s) extract exponent and significand
1464 eifrac (e, &l, frac) e to HOST_WIDE_INT and e type fraction
1465 euifrac (e, &l, frac) e to unsigned HOST_WIDE_INT and e type fraction
1466 einfin (e) set e to infinity, leaving its sign alone
1467 eldexp (a, n, b) multiply by 2**n
1469 emul (a, b, c) c = b * a
1472 eround (a, b) b = nearest integer value to a
1474 esub (a, b, c) c = b - a
1476 e24toasc (&f, str, n) single to ASCII string, n digits after decimal
1477 e53toasc (&d, str, n) double to ASCII string, n digits after decimal
1478 e64toasc (&d, str, n) 80-bit long double to ASCII string
1479 e113toasc (&d, str, n) 128-bit long double to ASCII string
1481 etoasc (e, str, n) e to ASCII string, n digits after decimal
1482 etoe24 (e, &f) convert e type to IEEE single precision
1483 etoe53 (e, &d) convert e type to IEEE double precision
1484 etoe64 (e, &d) convert e type to IEEE long double precision
1485 ltoe (&l, e) HOST_WIDE_INT to e type
1486 ultoe (&l, e) unsigned HOST_WIDE_INT to e type
1487 eisneg (e) 1 if sign bit of e != 0, else 0
1488 eisinf (e) 1 if e has maximum exponent (non-IEEE)
1489 or is infinite (IEEE)
1490 eisnan (e) 1 if e is a NaN
1493 Routines for internal format exploded e-type numbers
1495 eaddm (ai, bi) add significands, bi = bi + ai
1497 ecleazs (ei) set ei = 0 but leave its sign alone
1498 ecmpm (ai, bi) compare significands, return 1, 0, or -1
1499 edivm (ai, bi) divide significands, bi = bi / ai
1500 emdnorm (ai,l,s,exp) normalize and round off
1501 emovi (a, ai) convert external a to internal ai
1502 emovo (ai, a) convert internal ai to external a
1503 emovz (ai, bi) bi = ai, low guard word of bi = 0
1504 emulm (ai, bi) multiply significands, bi = bi * ai
1505 enormlz (ei) left-justify the significand
1506 eshdn1 (ai) shift significand and guards down 1 bit
1507 eshdn8 (ai) shift down 8 bits
1508 eshdn6 (ai) shift down 16 bits
1509 eshift (ai, n) shift ai n bits up (or down if n < 0)
1510 eshup1 (ai) shift significand and guards up 1 bit
1511 eshup8 (ai) shift up 8 bits
1512 eshup6 (ai) shift up 16 bits
1513 esubm (ai, bi) subtract significands, bi = bi - ai
1514 eiisinf (ai) 1 if infinite
1515 eiisnan (ai) 1 if a NaN
1516 eiisneg (ai) 1 if sign bit of ai != 0, else 0
1517 einan (ai) set ai = NaN
1519 eiinfin (ai) set ai = infinity
1522 The result is always normalized and rounded to NI-4 word precision
1523 after each arithmetic operation.
1525 Exception flags are NOT fully supported.
1527 Signaling NaN's are NOT supported; they are treated the same
1530 Define INFINITY for support of infinity; otherwise a
1531 saturation arithmetic is implemented.
1533 Define NANS for support of Not-a-Number items; otherwise the
1534 arithmetic will never produce a NaN output, and might be confused
1536 If NaN's are supported, the output of `ecmp (a,b)' is -2 if
1537 either a or b is a NaN. This means asking `if (ecmp (a,b) < 0)'
1538 may not be legitimate. Use `if (ecmp (a,b) == -1)' for `less than'
1541 Denormals are always supported here where appropriate (e.g., not
1542 for conversion to DEC numbers). */
1544 /* Definitions for error codes that are passed to the common error handling
1547 For Digital Equipment PDP-11 and VAX computers, certain
1548 IBM systems, and others that use numbers with a 56-bit
1549 significand, the symbol DEC should be defined. In this
1550 mode, most floating point constants are given as arrays
1551 of octal integers to eliminate decimal to binary conversion
1552 errors that might be introduced by the compiler.
1554 For computers, such as IBM PC, that follow the IEEE
1555 Standard for Binary Floating Point Arithmetic (ANSI/IEEE
1556 Std 754-1985), the symbol IEEE should be defined.
1557 These numbers have 53-bit significands. In this mode, constants
1558 are provided as arrays of hexadecimal 16 bit integers.
1559 The endian-ness of generated values is controlled by
1560 REAL_WORDS_BIG_ENDIAN.
1562 To accommodate other types of computer arithmetic, all
1563 constants are also provided in a normal decimal radix
1564 which one can hope are correctly converted to a suitable
1565 format by the available C language compiler. To invoke
1566 this mode, the symbol UNK is defined.
1568 An important difference among these modes is a predefined
1569 set of machine arithmetic constants for each. The numbers
1570 MACHEP (the machine roundoff error), MAXNUM (largest number
1571 represented), and several other parameters are preset by
1572 the configuration symbol. Check the file const.c to
1573 ensure that these values are correct for your computer.
1575 For ANSI C compatibility, define ANSIC equal to 1. Currently
1576 this affects only the atan2 function and others that use it. */
1578 /* Constant definitions for math error conditions. */
1580 #define DOMAIN 1 /* argument domain error */
1581 #define SING 2 /* argument singularity */
1582 #define OVERFLOW 3 /* overflow range error */
1583 #define UNDERFLOW 4 /* underflow range error */
1584 #define TLOSS 5 /* total loss of precision */
1585 #define PLOSS 6 /* partial loss of precision */
1586 #define INVALID 7 /* NaN-producing operation */
1588 /* e type constants used by high precision check routines */
1590 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
1592 const UEMUSHORT ezero
[NE
] =
1593 {0x0000, 0x0000, 0x0000, 0x0000,
1594 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000,};
1597 const UEMUSHORT ehalf
[NE
] =
1598 {0x0000, 0x0000, 0x0000, 0x0000,
1599 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3ffe,};
1602 const UEMUSHORT eone
[NE
] =
1603 {0x0000, 0x0000, 0x0000, 0x0000,
1604 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x3fff,};
1607 const UEMUSHORT etwo
[NE
] =
1608 {0x0000, 0x0000, 0x0000, 0x0000,
1609 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4000,};
1612 const UEMUSHORT e32
[NE
] =
1613 {0x0000, 0x0000, 0x0000, 0x0000,
1614 0x0000, 0x0000, 0x0000, 0x0000, 0x8000, 0x4004,};
1616 /* 6.93147180559945309417232121458176568075500134360255E-1 */
1617 const UEMUSHORT elog2
[NE
] =
1618 {0x40f3, 0xf6af, 0x03f2, 0xb398,
1619 0xc9e3, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1621 /* 1.41421356237309504880168872420969807856967187537695E0 */
1622 const UEMUSHORT esqrt2
[NE
] =
1623 {0x1d6f, 0xbe9f, 0x754a, 0x89b3,
1624 0x597d, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1626 /* 3.14159265358979323846264338327950288419716939937511E0 */
1627 const UEMUSHORT epi
[NE
] =
1628 {0x2902, 0x1cd1, 0x80dc, 0x628b,
1629 0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1632 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
1633 const UEMUSHORT ezero
[NE
] =
1634 {0, 0000000, 0000000, 0000000, 0000000, 0000000,};
1635 const UEMUSHORT ehalf
[NE
] =
1636 {0, 0000000, 0000000, 0000000, 0100000, 0x3ffe,};
1637 const UEMUSHORT eone
[NE
] =
1638 {0, 0000000, 0000000, 0000000, 0100000, 0x3fff,};
1639 const UEMUSHORT etwo
[NE
] =
1640 {0, 0000000, 0000000, 0000000, 0100000, 0040000,};
1641 const UEMUSHORT e32
[NE
] =
1642 {0, 0000000, 0000000, 0000000, 0100000, 0040004,};
1643 const UEMUSHORT elog2
[NE
] =
1644 {0xc9e4, 0x79ab, 0150717, 0013767, 0130562, 0x3ffe,};
1645 const UEMUSHORT esqrt2
[NE
] =
1646 {0x597e, 0x6484, 0174736, 0171463, 0132404, 0x3fff,};
1647 const UEMUSHORT epi
[NE
] =
1648 {0xc4c6, 0xc234, 0020550, 0155242, 0144417, 0040000,};
1651 /* Control register for rounding precision.
1652 This can be set to 113 (if NE=10), 80 (if NE=6), 64, 56, 53, or 24 bits. */
1657 /* Clear out entire e-type number X. */
1665 for (i
= 0; i
< NE
; i
++)
1669 /* Move e-type number from A to B. */
1678 for (i
= 0; i
< NE
; i
++)
1684 /* Absolute value of e-type X. */
1690 /* sign is top bit of last word of external format */
1691 x
[NE
- 1] &= 0x7fff;
1695 /* Negate the e-type number X. */
1702 x
[NE
- 1] ^= 0x8000; /* Toggle the sign bit */
1705 /* Return 1 if sign bit of e-type number X is nonzero, else zero. */
1709 const UEMUSHORT x
[];
1712 if (x
[NE
- 1] & 0x8000)
1718 /* Return 1 if e-type number X is infinity, else return zero. */
1722 const UEMUSHORT x
[];
1729 if ((x
[NE
- 1] & 0x7fff) == 0x7fff)
1735 /* Check if e-type number is not a number. The bit pattern is one that we
1736 defined, so we know for sure how to detect it. */
1740 const UEMUSHORT x
[] ATTRIBUTE_UNUSED
;
1745 /* NaN has maximum exponent */
1746 if ((x
[NE
- 1] & 0x7fff) != 0x7fff)
1748 /* ... and non-zero significand field. */
1749 for (i
= 0; i
< NE
- 1; i
++)
1759 /* Fill e-type number X with infinity pattern (IEEE)
1760 or largest possible number (non-IEEE). */
1769 for (i
= 0; i
< NE
- 1; i
++)
1773 for (i
= 0; i
< NE
- 1; i
++)
1801 /* Output an e-type NaN.
1802 This generates Intel's quiet NaN pattern for extended real.
1803 The exponent is 7fff, the leading mantissa word is c000. */
1813 for (i
= 0; i
< NE
- 2; i
++)
1816 *x
= (sign
<< 15) | 0x7fff;
1820 /* Move in an e-type number A, converting it to exploded e-type B. */
1832 p
= a
+ (NE
- 1); /* point to last word of external number */
1833 /* get the sign bit */
1838 /* get the exponent */
1840 *q
++ &= 0x7fff; /* delete the sign bit */
1842 if ((*(q
- 1) & 0x7fff) == 0x7fff)
1848 for (i
= 3; i
< NI
; i
++)
1854 for (i
= 2; i
< NI
; i
++)
1860 /* clear high guard word */
1862 /* move in the significand */
1863 for (i
= 0; i
< NE
- 1; i
++)
1865 /* clear low guard word */
1869 /* Move out exploded e-type number A, converting it to e type B. */
1882 q
= b
+ (NE
- 1); /* point to output exponent */
1883 /* combine sign and exponent */
1886 *q
-- = *p
++ | 0x8000;
1890 if (*(p
- 1) == 0x7fff)
1895 enan (b
, eiisneg (a
));
1903 /* skip over guard word */
1905 /* move the significand */
1906 for (j
= 0; j
< NE
- 1; j
++)
1910 /* Clear out exploded e-type number XI. */
1918 for (i
= 0; i
< NI
; i
++)
1922 /* Clear out exploded e-type XI, but don't touch the sign. */
1931 for (i
= 0; i
< NI
- 1; i
++)
1935 /* Move exploded e-type number from A to B. */
1944 for (i
= 0; i
< NI
- 1; i
++)
1946 /* clear low guard word */
1950 /* Generate exploded e-type NaN.
1951 The explicit pattern for this is maximum exponent and
1952 top two significant bits set. */
1966 /* Return nonzero if exploded e-type X is a NaN. */
1971 const UEMUSHORT x
[];
1975 if ((x
[E
] & 0x7fff) == 0x7fff)
1977 for (i
= M
+ 1; i
< NI
; i
++)
1987 /* Return nonzero if sign of exploded e-type X is nonzero. */
1991 const UEMUSHORT x
[];
1998 /* Fill exploded e-type X with infinity pattern.
1999 This has maximum exponent and significand all zeros. */
2011 /* Return nonzero if exploded e-type X is infinite. */
2016 const UEMUSHORT x
[];
2023 if ((x
[E
] & 0x7fff) == 0x7fff)
2027 #endif /* INFINITY */
2029 /* Compare significands of numbers in internal exploded e-type format.
2030 Guard words are included in the comparison.
2038 const UEMUSHORT
*a
, *b
;
2042 a
+= M
; /* skip up to significand area */
2044 for (i
= M
; i
< NI
; i
++)
2052 if (*(--a
) > *(--b
))
2058 /* Shift significand of exploded e-type X down by 1 bit. */
2067 x
+= M
; /* point to significand area */
2070 for (i
= M
; i
< NI
; i
++)
2082 /* Shift significand of exploded e-type X up by 1 bit. */
2094 for (i
= M
; i
< NI
; i
++)
2107 /* Shift significand of exploded e-type X down by 8 bits. */
2113 UEMUSHORT newbyt
, oldbyt
;
2118 for (i
= M
; i
< NI
; i
++)
2128 /* Shift significand of exploded e-type X up by 8 bits. */
2135 UEMUSHORT newbyt
, oldbyt
;
2140 for (i
= M
; i
< NI
; i
++)
2150 /* Shift significand of exploded e-type X up by 16 bits. */
2162 for (i
= M
; i
< NI
- 1; i
++)
2168 /* Shift significand of exploded e-type X down by 16 bits. */
2180 for (i
= M
; i
< NI
- 1; i
++)
2186 /* Add significands of exploded e-type X and Y. X + Y replaces Y. */
2200 for (i
= M
; i
< NI
; i
++)
2202 a
= (unsigned EMULONG
) (*x
) + (unsigned EMULONG
) (*y
) + carry
;
2213 /* Subtract significands of exploded e-type X and Y. Y - X replaces Y. */
2227 for (i
= M
; i
< NI
; i
++)
2229 a
= (unsigned EMULONG
) (*y
) - (unsigned EMULONG
) (*x
) - carry
;
2241 static UEMUSHORT equot
[NI
];
2245 /* Radix 2 shift-and-add versions of multiply and divide */
2248 /* Divide significands */
2252 UEMUSHORT den
[], num
[];
2262 for (i
= M
; i
< NI
; i
++)
2267 /* Use faster compare and subtraction if denominator has only 15 bits of
2273 for (i
= M
+ 3; i
< NI
; i
++)
2278 if ((den
[M
+ 1] & 1) != 0)
2286 for (i
= 0; i
< NBITS
+ 2; i
++)
2304 /* The number of quotient bits to calculate is NBITS + 1 scaling guard
2305 bit + 1 roundoff bit. */
2310 for (i
= 0; i
< NBITS
+ 2; i
++)
2312 if (ecmpm (den
, num
) <= 0)
2315 j
= 1; /* quotient bit = 1 */
2329 /* test for nonzero remainder after roundoff bit */
2332 for (i
= M
; i
< NI
; i
++)
2340 for (i
= 0; i
< NI
; i
++)
2346 /* Multiply significands */
2357 for (i
= M
; i
< NI
; i
++)
2362 while (*p
== 0) /* significand is not supposed to be zero */
2367 if ((*p
& 0xff) == 0)
2375 for (i
= 0; i
< k
; i
++)
2379 /* remember if there were any nonzero bits shifted out */
2386 for (i
= 0; i
< NI
; i
++)
2389 /* return flag for lost nonzero bits */
2395 /* Radix 65536 versions of multiply and divide. */
2397 /* Multiply significand of e-type number B
2398 by 16-bit quantity A, return e-type result to C. */
2403 const UEMUSHORT b
[];
2407 unsigned EMULONG carry
;
2408 const UEMUSHORT
*ps
;
2410 unsigned EMULONG aa
, m
;
2419 for (i
=M
+1; i
<NI
; i
++)
2429 m
= (unsigned EMULONG
) aa
* *ps
--;
2430 carry
= (m
& 0xffff) + *pp
;
2431 *pp
-- = (UEMUSHORT
) carry
;
2432 carry
= (carry
>> 16) + (m
>> 16) + *pp
;
2433 *pp
= (UEMUSHORT
) carry
;
2434 *(pp
-1) = carry
>> 16;
2437 for (i
=M
; i
<NI
; i
++)
2441 /* Divide significands of exploded e-types NUM / DEN. Neither the
2442 numerator NUM nor the denominator DEN is permitted to have its high guard
2447 const UEMUSHORT den
[];
2452 unsigned EMULONG tnum
;
2453 UEMUSHORT j
, tdenm
, tquot
;
2454 UEMUSHORT tprod
[NI
+1];
2460 for (i
=M
; i
<NI
; i
++)
2466 for (i
=M
; i
<NI
; i
++)
2468 /* Find trial quotient digit (the radix is 65536). */
2469 tnum
= (((unsigned EMULONG
) num
[M
]) << 16) + num
[M
+1];
2471 /* Do not execute the divide instruction if it will overflow. */
2472 if ((tdenm
* (unsigned long) 0xffff) < tnum
)
2475 tquot
= tnum
/ tdenm
;
2476 /* Multiply denominator by trial quotient digit. */
2477 m16m ((unsigned int) tquot
, den
, tprod
);
2478 /* The quotient digit may have been overestimated. */
2479 if (ecmpm (tprod
, num
) > 0)
2483 if (ecmpm (tprod
, num
) > 0)
2493 /* test for nonzero remainder after roundoff bit */
2496 for (i
=M
; i
<NI
; i
++)
2503 for (i
=0; i
<NI
; i
++)
2509 /* Multiply significands of exploded e-type A and B, result in B. */
2513 const UEMUSHORT a
[];
2518 UEMUSHORT pprod
[NI
];
2524 for (i
=M
; i
<NI
; i
++)
2530 for (i
=M
+1; i
<NI
; i
++)
2538 m16m ((unsigned int) *p
--, b
, pprod
);
2539 eaddm (pprod
, equot
);
2545 for (i
=0; i
<NI
; i
++)
2548 /* return flag for lost nonzero bits */
2554 /* Normalize and round off.
2556 The internal format number to be rounded is S.
2557 Input LOST is 0 if the value is exact. This is the so-called sticky bit.
2559 Input SUBFLG indicates whether the number was obtained
2560 by a subtraction operation. In that case if LOST is nonzero
2561 then the number is slightly smaller than indicated.
2563 Input EXP is the biased exponent, which may be negative.
2564 the exponent field of S is ignored but is replaced by
2565 EXP as adjusted by normalization and rounding.
2567 Input RCNTRL is the rounding control. If it is nonzero, the
2568 returned value will be rounded to RNDPRC bits.
2570 For future reference: In order for emdnorm to round off denormal
2571 significands at the right point, the input exponent must be
2572 adjusted to be the actual value it would have after conversion to
2573 the final floating point type. This adjustment has been
2574 implemented for all type conversions (etoe53, etc.) and decimal
2575 conversions, but not for the arithmetic functions (eadd, etc.).
2576 Data types having standard 15-bit exponents are not affected by
2577 this, but SFmode and DFmode are affected. For example, ediv with
2578 rndprc = 24 will not round correctly to 24-bit precision if the
2579 result is denormal. */
2581 static int rlast
= -1;
2583 static UEMUSHORT rmsk
= 0;
2584 static UEMUSHORT rmbit
= 0;
2585 static UEMUSHORT rebit
= 0;
2587 static UEMUSHORT rbit
[NI
];
2590 emdnorm (s
, lost
, subflg
, exp
, rcntrl
)
2603 /* a blank significand could mean either zero or infinity. */
2616 if ((j
> NBITS
) && (exp
< 32767))
2624 if (exp
> (EMULONG
) (-NBITS
- 1))
2637 /* Round off, unless told not to by rcntrl. */
2640 /* Set up rounding parameters if the control register changed. */
2641 if (rndprc
!= rlast
)
2648 rw
= NI
- 1; /* low guard word */
2671 /* For DEC or IBM arithmetic */
2688 /* For C4x arithmetic */
2709 /* Shift down 1 temporarily if the data structure has an implied
2710 most significant bit and the number is denormal.
2711 Intel long double denormals also lose one bit of precision. */
2712 if ((exp
<= 0) && (rndprc
!= NBITS
)
2713 && ((rndprc
!= 64) || ((rndprc
== 64) && ! REAL_WORDS_BIG_ENDIAN
)))
2715 lost
|= s
[NI
- 1] & 1;
2718 /* Clear out all bits below the rounding bit,
2719 remembering in r if any were nonzero. */
2733 if ((r
& rmbit
) != 0)
2739 { /* round to even */
2740 if ((s
[re
] & rebit
) == 0)
2753 /* Undo the temporary shift for denormal values. */
2754 if ((exp
<= 0) && (rndprc
!= NBITS
)
2755 && ((rndprc
!= 64) || ((rndprc
== 64) && ! REAL_WORDS_BIG_ENDIAN
)))
2760 { /* overflow on roundoff */
2773 for (i
= 2; i
< NI
- 1; i
++)
2776 warning ("floating point overflow");
2780 for (i
= M
+ 1; i
< NI
- 1; i
++)
2783 if ((rndprc
< 64) || (rndprc
== 113))
2798 s
[1] = (UEMUSHORT
) exp
;
2801 /* Subtract. C = B - A, all e type numbers. */
2803 static int subflg
= 0;
2807 const UEMUSHORT
*a
, *b
;
2822 /* Infinity minus infinity is a NaN.
2823 Test for subtracting infinities of the same sign. */
2824 if (eisinf (a
) && eisinf (b
)
2825 && ((eisneg (a
) ^ eisneg (b
)) == 0))
2827 mtherr ("esub", INVALID
);
2836 /* Add. C = A + B, all e type. */
2840 const UEMUSHORT
*a
, *b
;
2845 /* NaN plus anything is a NaN. */
2856 /* Infinity minus infinity is a NaN.
2857 Test for adding infinities of opposite signs. */
2858 if (eisinf (a
) && eisinf (b
)
2859 && ((eisneg (a
) ^ eisneg (b
)) != 0))
2861 mtherr ("esub", INVALID
);
2870 /* Arithmetic common to both addition and subtraction. */
2874 const UEMUSHORT
*a
, *b
;
2877 UEMUSHORT ai
[NI
], bi
[NI
], ci
[NI
];
2879 EMULONG lt
, lta
, ltb
;
2900 /* compare exponents */
2905 { /* put the larger number in bi */
2915 if (lt
< (EMULONG
) (-NBITS
- 1))
2916 goto done
; /* answer same as larger addend */
2918 lost
= eshift (ai
, k
); /* shift the smaller number down */
2922 /* exponents were the same, so must compare significands */
2925 { /* the numbers are identical in magnitude */
2926 /* if different signs, result is zero */
2932 /* if same sign, result is double */
2933 /* double denormalized tiny number */
2934 if ((bi
[E
] == 0) && ((bi
[3] & 0x8000) == 0))
2939 /* add 1 to exponent unless both are zero! */
2940 for (j
= 1; j
< NI
- 1; j
++)
2956 bi
[E
] = (UEMUSHORT
) ltb
;
2960 { /* put the larger number in bi */
2976 emdnorm (bi
, lost
, subflg
, ltb
, !ROUND_TOWARDS_ZERO
);
2982 /* Divide: C = B/A, all e type. */
2986 const UEMUSHORT
*a
, *b
;
2989 UEMUSHORT ai
[NI
], bi
[NI
];
2991 EMULONG lt
, lta
, ltb
;
2993 /* IEEE says if result is not a NaN, the sign is "-" if and only if
2994 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
2995 sign
= eisneg (a
) ^ eisneg (b
);
2998 /* Return any NaN input. */
3009 /* Zero over zero, or infinity over infinity, is a NaN. */
3010 if (((ecmp (a
, ezero
) == 0) && (ecmp (b
, ezero
) == 0))
3011 || (eisinf (a
) && eisinf (b
)))
3013 mtherr ("ediv", INVALID
);
3018 /* Infinity over anything else is infinity. */
3025 /* Anything else over infinity is zero. */
3037 { /* See if numerator is zero. */
3038 for (i
= 1; i
< NI
- 1; i
++)
3042 ltb
-= enormlz (bi
);
3052 { /* possible divide by zero */
3053 for (i
= 1; i
< NI
- 1; i
++)
3057 lta
-= enormlz (ai
);
3061 /* Divide by zero is not an invalid operation.
3062 It is a divide-by-zero operation! */
3064 mtherr ("ediv", SING
);
3070 /* calculate exponent */
3071 lt
= ltb
- lta
+ EXONE
;
3072 emdnorm (bi
, i
, 0, lt
, !ROUND_TOWARDS_ZERO
);
3079 && (ecmp (c
, ezero
) != 0)
3082 *(c
+(NE
-1)) |= 0x8000;
3084 *(c
+(NE
-1)) &= ~0x8000;
3087 /* Multiply e-types A and B, return e-type product C. */
3091 const UEMUSHORT
*a
, *b
;
3094 UEMUSHORT ai
[NI
], bi
[NI
];
3096 EMULONG lt
, lta
, ltb
;
3098 /* IEEE says if result is not a NaN, the sign is "-" if and only if
3099 operands have opposite signs -- but flush -0 to 0 later if not IEEE. */
3100 sign
= eisneg (a
) ^ eisneg (b
);
3103 /* NaN times anything is the same NaN. */
3114 /* Zero times infinity is a NaN. */
3115 if ((eisinf (a
) && (ecmp (b
, ezero
) == 0))
3116 || (eisinf (b
) && (ecmp (a
, ezero
) == 0)))
3118 mtherr ("emul", INVALID
);
3123 /* Infinity times anything else is infinity. */
3125 if (eisinf (a
) || eisinf (b
))
3137 for (i
= 1; i
< NI
- 1; i
++)
3141 lta
-= enormlz (ai
);
3152 for (i
= 1; i
< NI
- 1; i
++)
3156 ltb
-= enormlz (bi
);
3165 /* Multiply significands */
3167 /* calculate exponent */
3168 lt
= lta
+ ltb
- (EXONE
- 1);
3169 emdnorm (bi
, j
, 0, lt
, !ROUND_TOWARDS_ZERO
);
3176 && (ecmp (c
, ezero
) != 0)
3179 *(c
+(NE
-1)) |= 0x8000;
3181 *(c
+(NE
-1)) &= ~0x8000;
3184 /* Convert double precision PE to e-type Y. */
3188 const UEMUSHORT
*pe
;
3198 ibmtoe (pe
, y
, DFmode
);
3203 c4xtoe (pe
, y
, HFmode
);
3207 ieeetoe (pe
, y
, &ieee_53
);
3209 #endif /* not C4X */
3210 #endif /* not IBM */
3211 #endif /* not DEC */
3214 /* Convert double extended precision float PE to e type Y. */
3218 const UEMUSHORT
*pe
;
3228 for (i
= 0; i
< NE
- 5; i
++)
3231 /* REAL_WORDS_BIG_ENDIAN is always 0 for DEC and 1 for IBM.
3232 This precision is not ordinarily supported on DEC or IBM. */
3233 if (! REAL_WORDS_BIG_ENDIAN
)
3235 for (i
= 0; i
< 5; i
++)
3239 /* For denormal long double Intel format, shift significand up one
3240 -- but only if the top significand bit is zero. A top bit of 1
3241 is "pseudodenormal" when the exponent is zero. */
3242 if ((yy
[NE
-1] & 0x7fff) == 0 && (yy
[NE
-2] & 0x8000) == 0)
3255 p
= &yy
[0] + (NE
- 1);
3258 for (i
= 0; i
< 4; i
++)
3261 #endif /* not C4X */
3263 /* Point to the exponent field and check max exponent cases. */
3265 if ((*p
& 0x7fff) == 0x7fff)
3268 if (! REAL_WORDS_BIG_ENDIAN
)
3270 for (i
= 0; i
< 4; i
++)
3272 if ((i
!= 3 && pe
[i
] != 0)
3273 /* Anything but 0x8000 here, including 0, is a NaN. */
3274 || (i
== 3 && pe
[i
] != 0x8000))
3276 enan (y
, (*p
& 0x8000) != 0);
3283 /* In Motorola extended precision format, the most significant
3284 bit of an infinity mantissa could be either 1 or 0. It is
3285 the lower order bits that tell whether the value is a NaN. */
3286 if ((pe
[2] & 0x7fff) != 0)
3289 for (i
= 3; i
<= 5; i
++)
3294 enan (y
, (*p
& 0x8000) != 0);
3306 #endif /* INFINITY */
3309 for (i
= 0; i
< NE
; i
++)
3313 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3314 /* Convert 128-bit long double precision float PE to e type Y. */
3318 const UEMUSHORT
*pe
;
3321 ieeetoe (pe
, y
, &ieee_113
);
3323 #endif /* INTEL_EXTENDED_IEEE_FORMAT == 0 */
3325 /* Convert single precision float PE to e type Y. */
3329 const UEMUSHORT
*pe
;
3334 ibmtoe (pe
, y
, SFmode
);
3340 c4xtoe (pe
, y
, QFmode
);
3345 ieeetoe (pe
, y
, &dec_f
);
3349 ieeetoe (pe
, y
, &ieee_24
);
3351 #endif /* not DEC */
3352 #endif /* not C4X */
3353 #endif /* not IBM */
3356 /* Convert machine format float of specified format PE to e type Y. */
3359 ieeetoe (pe
, y
, fmt
)
3360 const UEMUSHORT
*pe
;
3362 const struct ieee_format
*fmt
;
3369 int shortsm1
= fmt
->bits
/ 16 - 1;
3371 int expmask
= (1 << fmt
->expbits
) - 1;
3373 int expshift
= (fmt
->precision
- 1) & 0x0f;
3374 int highbit
= 1 << expshift
;
3379 if (! REAL_WORDS_BIG_ENDIAN
)
3385 yy
[M
] = (r
& (highbit
- 1)) | highbit
;
3386 r
= (r
& 0x7fff) >> expshift
;
3388 if (!LARGEST_EXPONENT_IS_NORMAL (fmt
->precision
) && r
== expmask
)
3391 /* First check the word where high order mantissa and exponent live */
3392 if ((*e
& (highbit
- 1)) != 0)
3394 enan (y
, yy
[0] != 0);
3397 if (! REAL_WORDS_BIG_ENDIAN
)
3399 for (i
= 0; i
< shortsm1
; i
++)
3403 enan (y
, yy
[0] != 0);
3410 for (i
= 1; i
< shortsm1
+ 1; i
++)
3414 enan (y
, yy
[0] != 0);
3426 #endif /* INFINITY */
3427 /* If zero exponent, then the significand is denormalized.
3428 So take back the understood high significand bit. */
3434 r
+= fmt
->adjustment
;
3437 if (! REAL_WORDS_BIG_ENDIAN
)
3439 for (i
= 0; i
< shortsm1
; i
++)
3445 for (i
= 0; i
< shortsm1
; i
++)
3448 if (fmt
->precision
== 113)
3450 /* denorm is left alone in 113 bit format */
3456 eshift (yy
, -(expshift
+ 1));
3458 { /* if zero exponent, then normalize the significand */
3459 if ((k
= enormlz (yy
)) > NBITS
)
3462 yy
[E
] -= (UEMUSHORT
) (k
- 1);
3468 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
3469 /* Convert e-type X to IEEE 128-bit long double format E. */
3476 etoieee (x
, e
, &ieee_113
);
3479 /* Convert exploded e-type X, that has already been rounded to
3480 113-bit precision, to IEEE 128-bit long double format Y. */
3486 toieee (x
, y
, &ieee_113
);
3489 #endif /* INTEL_EXTENDED_IEEE_FORMAT == 0 */
3491 /* Convert e-type X to IEEE double extended format E. */
3498 etoieee (x
, e
, &ieee_64
);
3501 /* Convert exploded e-type X, that has already been rounded to
3502 64-bit precision, to IEEE double extended format Y. */
3508 toieee (x
, y
, &ieee_64
);
3511 /* e type to double precision. */
3514 /* Convert e-type X to DEC-format double E. */
3524 /* Convert exploded e-type X, that has already been rounded to
3525 56-bit double precision, to DEC double Y. */
3536 /* Convert e-type X to IBM 370-format double E. */
3543 etoibm (x
, e
, DFmode
);
3546 /* Convert exploded e-type X, that has already been rounded to
3547 56-bit precision, to IBM 370 double Y. */
3553 toibm (x
, y
, DFmode
);
3556 #else /* it's neither DEC nor IBM */
3558 /* Convert e-type X to C4X-format long double E. */
3565 etoc4x (x
, e
, HFmode
);
3568 /* Convert exploded e-type X, that has already been rounded to
3569 56-bit precision, to IBM 370 double Y. */
3575 toc4x (x
, y
, HFmode
);
3578 #else /* it's neither DEC nor IBM nor C4X */
3580 /* Convert e-type X to IEEE double E. */
3587 etoieee (x
, e
, &ieee_53
);
3590 /* Convert exploded e-type X, that has already been rounded to
3591 53-bit precision, to IEEE double Y. */
3597 toieee (x
, y
, &ieee_53
);
3600 #endif /* not C4X */
3601 #endif /* not IBM */
3602 #endif /* not DEC */
3606 /* e type to single precision. */
3609 /* Convert e-type X to IBM 370 float E. */
3616 etoibm (x
, e
, SFmode
);
3619 /* Convert exploded e-type X, that has already been rounded to
3620 float precision, to IBM 370 float Y. */
3626 toibm (x
, y
, SFmode
);
3629 #else /* it's not IBM */
3632 /* Convert e-type X to C4X float E. */
3639 etoc4x (x
, e
, QFmode
);
3642 /* Convert exploded e-type X, that has already been rounded to
3643 float precision, to IBM 370 float Y. */
3649 toc4x (x
, y
, QFmode
);
3652 #else /* it's neither IBM nor C4X */
3656 /* Convert e-type X to DEC F-float E. */
3663 etoieee (x
, e
, &dec_f
);
3666 /* Convert exploded e-type X, that has already been rounded to
3667 float precision, to DEC F-float Y. */
3673 toieee (x
, y
, &dec_f
);
3678 /* Convert e-type X to IEEE float E. */
3685 etoieee (x
, e
, &ieee_24
);
3688 /* Convert exploded e-type X, that has already been rounded to
3689 float precision, to IEEE float Y. */
3695 toieee (x
, y
, &ieee_24
);
3698 #endif /* not DEC */
3699 #endif /* not C4X */
3700 #endif /* not IBM */
3703 /* Convert e-type X to the IEEE format described by FMT. */
3709 const struct ieee_format
*fmt
;
3718 make_nan (e
, eisneg (x
), fmt
->mode
);
3729 /* Adjust exponent for offset. */
3730 exp
= (EMULONG
) xi
[E
] - fmt
->adjustment
;
3732 /* Round off to nearest or even. */
3734 rndprc
= fmt
->precision
;
3735 emdnorm (xi
, 0, 0, exp
, !ROUND_TOWARDS_ZERO
);
3740 toieee (xi
, e
, fmt
);
3743 /* Convert exploded e-type X, that has already been rounded to
3744 the necessary precision, to the IEEE format described by FMT. */
3749 const struct ieee_format
*fmt
;
3756 maxexp
= (1 << fmt
->expbits
) - 1;
3757 words
= (fmt
->bits
- fmt
->expbits
) / EMUSHORT_SIZE
;
3762 make_nan (y
, eiisneg (x
), fmt
->mode
);
3767 if (fmt
->expbits
< 15
3768 && LARGEST_EXPONENT_IS_NORMAL (fmt
->bits
)
3771 saturate (y
, eiisneg (x
), fmt
->bits
, 1);
3775 /* Point to the exponent. */
3776 if (REAL_WORDS_BIG_ENDIAN
)
3781 /* Copy the sign. */
3787 if (fmt
->expbits
< 15
3788 && !LARGEST_EXPONENT_IS_NORMAL (fmt
->bits
)
3791 /* Saturate at largest number less that infinity. */
3794 *q
|= maxexp
<< (15 - fmt
->expbits
);
3797 *q
|= (maxexp
<< (15 - fmt
->expbits
)) - 1;
3801 if (!REAL_WORDS_BIG_ENDIAN
)
3803 for (i
= 0; i
< words
; i
++)
3808 for (i
= 0; i
< words
; i
++)
3811 #if defined(INFINITY) && defined(ERANGE)
3817 /* If denormal and DEC float, return zero (DEC has no denormals) */
3821 for (i
= 0; i
< fmt
->bits
/ EMUSHORT_SIZE
; i
++)
3827 /* Delete the implied bit unless denormal, except for
3828 64-bit precision. */
3829 if (fmt
->precision
!= 64 && x
[E
] != 0)
3834 /* Shift denormal double extended Intel format significand down
3836 if (fmt
->precision
== 64 && x
[E
] == 0 && ! REAL_WORDS_BIG_ENDIAN
)
3839 if (fmt
->expbits
< 15)
3841 /* Shift the significand. */
3842 eshift (x
, 15 - fmt
->expbits
);
3844 /* Combine the exponent and upper bits of the significand. */
3845 *q
|= x
[E
] << (15 - fmt
->expbits
);
3846 *q
|= x
[M
] & (UEMUSHORT
) ~((maxexp
<< (15 - fmt
->expbits
)) | 0x8000);
3850 /* Copy the exponent. */
3854 /* Add padding after the exponent. At the moment, this is only necessary for
3855 64-bit precision; in this case, the padding is 16 bits. */
3856 if (fmt
->precision
== 64)
3861 if (REAL_WORDS_BIG_ENDIAN
)
3865 /* Copy the significand. */
3866 if (REAL_WORDS_BIG_ENDIAN
)
3868 for (i
= 0; i
< words
; i
++)
3869 *(++q
) = x
[i
+ M
+ 1];
3872 else if (fmt
->precision
== 64 && eiisinf (x
))
3874 /* Intel double extended infinity significand. */
3883 for (i
= 0; i
< words
; i
++)
3884 *(--q
) = x
[i
+ M
+ 1];
3889 /* Compare two e type numbers.
3893 -2 if either a or b is a NaN. */
3897 const UEMUSHORT
*a
, *b
;
3899 UEMUSHORT ai
[NI
], bi
[NI
];
3905 if (eisnan (a
) || eisnan (b
))
3914 { /* the signs are different */
3916 for (i
= 1; i
< NI
- 1; i
++)
3930 /* both are the same sign */
3945 return (0); /* equality */
3949 if (*(--p
) > *(--q
))
3950 return (msign
); /* p is bigger */
3952 return (-msign
); /* p is littler */
3956 /* Find e-type nearest integer to X, as floor (X + 0.5). */
3968 /* Convert HOST_WIDE_INT LP to e type Y. */
3972 const HOST_WIDE_INT
*lp
;
3976 unsigned HOST_WIDE_INT ll
;
3982 /* make it positive */
3983 ll
= (unsigned HOST_WIDE_INT
) (-(*lp
));
3984 yi
[0] = 0xffff; /* put correct sign in the e type number */
3988 ll
= (unsigned HOST_WIDE_INT
) (*lp
);
3990 /* move the long integer to yi significand area */
3991 #if HOST_BITS_PER_WIDE_INT == 64
3992 yi
[M
] = (UEMUSHORT
) (ll
>> 48);
3993 yi
[M
+ 1] = (UEMUSHORT
) (ll
>> 32);
3994 yi
[M
+ 2] = (UEMUSHORT
) (ll
>> 16);
3995 yi
[M
+ 3] = (UEMUSHORT
) ll
;
3996 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
3998 yi
[M
] = (UEMUSHORT
) (ll
>> 16);
3999 yi
[M
+ 1] = (UEMUSHORT
) ll
;
4000 yi
[E
] = EXONE
+ 15; /* exponent if normalize shift count were 0 */
4003 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
4004 ecleaz (yi
); /* it was zero */
4006 yi
[E
] -= (UEMUSHORT
) k
;/* subtract shift count from exponent */
4007 emovo (yi
, y
); /* output the answer */
4010 /* Convert unsigned HOST_WIDE_INT LP to e type Y. */
4014 const unsigned HOST_WIDE_INT
*lp
;
4018 unsigned HOST_WIDE_INT ll
;
4024 /* move the long integer to ayi significand area */
4025 #if HOST_BITS_PER_WIDE_INT == 64
4026 yi
[M
] = (UEMUSHORT
) (ll
>> 48);
4027 yi
[M
+ 1] = (UEMUSHORT
) (ll
>> 32);
4028 yi
[M
+ 2] = (UEMUSHORT
) (ll
>> 16);
4029 yi
[M
+ 3] = (UEMUSHORT
) ll
;
4030 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
4032 yi
[M
] = (UEMUSHORT
) (ll
>> 16);
4033 yi
[M
+ 1] = (UEMUSHORT
) ll
;
4034 yi
[E
] = EXONE
+ 15; /* exponent if normalize shift count were 0 */
4037 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
4038 ecleaz (yi
); /* it was zero */
4040 yi
[E
] -= (UEMUSHORT
) k
; /* subtract shift count from exponent */
4041 emovo (yi
, y
); /* output the answer */
4045 /* Find signed HOST_WIDE_INT integer I and floating point fractional
4046 part FRAC of e-type (packed internal format) floating point input X.
4047 The integer output I has the sign of the input, except that
4048 positive overflow is permitted if FIXUNS_TRUNC_LIKE_FIX_TRUNC.
4049 The output e-type fraction FRAC is the positive fractional
4060 unsigned HOST_WIDE_INT ll
;
4063 k
= (int) xi
[E
] - (EXONE
- 1);
4066 /* if exponent <= 0, integer = 0 and real output is fraction */
4071 if (k
> (HOST_BITS_PER_WIDE_INT
- 1))
4073 /* long integer overflow: output large integer
4074 and correct fraction */
4076 *i
= ((unsigned HOST_WIDE_INT
) 1) << (HOST_BITS_PER_WIDE_INT
- 1);
4079 #ifdef FIXUNS_TRUNC_LIKE_FIX_TRUNC
4080 /* In this case, let it overflow and convert as if unsigned. */
4081 euifrac (x
, &ll
, frac
);
4082 *i
= (HOST_WIDE_INT
) ll
;
4085 /* In other cases, return the largest positive integer. */
4086 *i
= (((unsigned HOST_WIDE_INT
) 1) << (HOST_BITS_PER_WIDE_INT
- 1)) - 1;
4091 warning ("overflow on truncation to integer");
4095 /* Shift more than 16 bits: first shift up k-16 mod 16,
4096 then shift up by 16's. */
4097 j
= k
- ((k
>> 4) << 4);
4104 ll
= (ll
<< 16) | xi
[M
];
4106 while ((k
-= 16) > 0);
4113 /* shift not more than 16 bits */
4115 *i
= (HOST_WIDE_INT
) xi
[M
] & 0xffff;
4122 if ((k
= enormlz (xi
)) > NBITS
)
4125 xi
[E
] -= (UEMUSHORT
) k
;
4131 /* Find unsigned HOST_WIDE_INT integer I and floating point fractional part
4132 FRAC of e-type X. A negative input yields integer output = 0 but
4133 correct fraction. */
4136 euifrac (x
, i
, frac
)
4138 unsigned HOST_WIDE_INT
*i
;
4141 unsigned HOST_WIDE_INT ll
;
4146 k
= (int) xi
[E
] - (EXONE
- 1);
4149 /* if exponent <= 0, integer = 0 and argument is fraction */
4154 if (k
> HOST_BITS_PER_WIDE_INT
)
4156 /* Long integer overflow: output large integer
4157 and correct fraction.
4158 Note, the BSD MicroVAX compiler says that ~(0UL)
4159 is a syntax error. */
4163 warning ("overflow on truncation to unsigned integer");
4167 /* Shift more than 16 bits: first shift up k-16 mod 16,
4168 then shift up by 16's. */
4169 j
= k
- ((k
>> 4) << 4);
4176 ll
= (ll
<< 16) | xi
[M
];
4178 while ((k
-= 16) > 0);
4183 /* shift not more than 16 bits */
4185 *i
= (HOST_WIDE_INT
) xi
[M
] & 0xffff;
4188 if (xi
[0]) /* A negative value yields unsigned integer 0. */
4194 if ((k
= enormlz (xi
)) > NBITS
)
4197 xi
[E
] -= (UEMUSHORT
) k
;
4202 /* Shift the significand of exploded e-type X up or down by SC bits. */
4223 lost
|= *p
; /* remember lost bits */
4264 return ((int) lost
);
4267 /* Shift normalize the significand area of exploded e-type X.
4268 Return the shift count (up = positive). */
4283 return (0); /* already normalized */
4289 /* With guard word, there are NBITS+16 bits available.
4290 Return true if all are zero. */
4294 /* see if high byte is zero */
4295 while ((*p
& 0xff00) == 0)
4300 /* now shift 1 bit at a time */
4301 while ((*p
& 0x8000) == 0)
4307 mtherr ("enormlz", UNDERFLOW
);
4313 /* Normalize by shifting down out of the high guard word
4314 of the significand */
4329 mtherr ("enormlz", OVERFLOW
);
4336 /* Powers of ten used in decimal <-> binary conversions. */
4341 #if MAX_LONG_DOUBLE_TYPE_SIZE == 128 && (INTEL_EXTENDED_IEEE_FORMAT == 0)
4342 static const UEMUSHORT etens
[NTEN
+ 1][NE
] =
4344 {0x6576, 0x4a92, 0x804a, 0x153f,
4345 0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4346 {0x6a32, 0xce52, 0x329a, 0x28ce,
4347 0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4348 {0x526c, 0x50ce, 0xf18b, 0x3d28,
4349 0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4350 {0x9c66, 0x58f8, 0xbc50, 0x5c54,
4351 0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4352 {0x851e, 0xeab7, 0x98fe, 0x901b,
4353 0xddbb, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4354 {0x0235, 0x0137, 0x36b1, 0x336c,
4355 0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4356 {0x50f8, 0x25fb, 0xc76b, 0x6b71,
4357 0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4358 {0x0000, 0x0000, 0x0000, 0x0000,
4359 0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4360 {0x0000, 0x0000, 0x0000, 0x0000,
4361 0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4362 {0x0000, 0x0000, 0x0000, 0x0000,
4363 0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4364 {0x0000, 0x0000, 0x0000, 0x0000,
4365 0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4366 {0x0000, 0x0000, 0x0000, 0x0000,
4367 0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4368 {0x0000, 0x0000, 0x0000, 0x0000,
4369 0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4372 static const UEMUSHORT emtens
[NTEN
+ 1][NE
] =
4374 {0x2030, 0xcffc, 0xa1c3, 0x8123,
4375 0x2de3, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4376 {0x8264, 0xd2cb, 0xf2ea, 0x12d4,
4377 0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4378 {0xf53f, 0xf698, 0x6bd3, 0x0158,
4379 0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4380 {0xe731, 0x04d4, 0xe3f2, 0xd332,
4381 0x7132, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4382 {0xa23e, 0x5308, 0xfefb, 0x1155,
4383 0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4384 {0xe26d, 0xdbde, 0xd05d, 0xb3f6,
4385 0xac7c, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4386 {0x2a20, 0x6224, 0x47b3, 0x98d7,
4387 0x3f23, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4388 {0x0b5b, 0x4af2, 0xa581, 0x18ed,
4389 0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4390 {0xbf71, 0xa9b3, 0x7989, 0xbe68,
4391 0x4c2e, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4392 {0x3d4d, 0x7c3d, 0x36ba, 0x0d2b,
4393 0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4394 {0xc155, 0xa4a8, 0x404e, 0x6113,
4395 0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4396 {0xd70a, 0x70a3, 0x0a3d, 0xa3d7,
4397 0x3d70, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4398 {0xcccd, 0xcccc, 0xcccc, 0xcccc,
4399 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4402 /* LONG_DOUBLE_TYPE_SIZE is other than 128 */
4403 static const UEMUSHORT etens
[NTEN
+ 1][NE
] =
4405 {0xc94c, 0x979a, 0x8a20, 0x5202, 0xc460, 0x7525,}, /* 10**4096 */
4406 {0xa74d, 0x5de4, 0xc53d, 0x3b5d, 0x9e8b, 0x5a92,}, /* 10**2048 */
4407 {0x650d, 0x0c17, 0x8175, 0x7586, 0xc976, 0x4d48,},
4408 {0xcc65, 0x91c6, 0xa60e, 0xa0ae, 0xe319, 0x46a3,},
4409 {0xddbc, 0xde8d, 0x9df9, 0xebfb, 0xaa7e, 0x4351,},
4410 {0xc66f, 0x8cdf, 0x80e9, 0x47c9, 0x93ba, 0x41a8,},
4411 {0x3cbf, 0xa6d5, 0xffcf, 0x1f49, 0xc278, 0x40d3,},
4412 {0xf020, 0xb59d, 0x2b70, 0xada8, 0x9dc5, 0x4069,},
4413 {0x0000, 0x0000, 0x0400, 0xc9bf, 0x8e1b, 0x4034,},
4414 {0x0000, 0x0000, 0x0000, 0x2000, 0xbebc, 0x4019,},
4415 {0x0000, 0x0000, 0x0000, 0x0000, 0x9c40, 0x400c,},
4416 {0x0000, 0x0000, 0x0000, 0x0000, 0xc800, 0x4005,},
4417 {0x0000, 0x0000, 0x0000, 0x0000, 0xa000, 0x4002,}, /* 10**1 */
4420 static const UEMUSHORT emtens
[NTEN
+ 1][NE
] =
4422 {0x2de4, 0x9fde, 0xd2ce, 0x04c8, 0xa6dd, 0x0ad8,}, /* 10**-4096 */
4423 {0x4925, 0x2de4, 0x3436, 0x534f, 0xceae, 0x256b,}, /* 10**-2048 */
4424 {0x87a6, 0xc0bd, 0xda57, 0x82a5, 0xa2a6, 0x32b5,},
4425 {0x7133, 0xd21c, 0xdb23, 0xee32, 0x9049, 0x395a,},
4426 {0xfa91, 0x1939, 0x637a, 0x4325, 0xc031, 0x3cac,},
4427 {0xac7d, 0xe4a0, 0x64bc, 0x467c, 0xddd0, 0x3e55,},
4428 {0x3f24, 0xe9a5, 0xa539, 0xea27, 0xa87f, 0x3f2a,},
4429 {0x67de, 0x94ba, 0x4539, 0x1ead, 0xcfb1, 0x3f94,},
4430 {0x4c2f, 0xe15b, 0xc44d, 0x94be, 0xe695, 0x3fc9,},
4431 {0xfdc2, 0xcefc, 0x8461, 0x7711, 0xabcc, 0x3fe4,},
4432 {0xd3c3, 0x652b, 0xe219, 0x1758, 0xd1b7, 0x3ff1,},
4433 {0x3d71, 0xd70a, 0x70a3, 0x0a3d, 0xa3d7, 0x3ff8,},
4434 {0xcccd, 0xcccc, 0xcccc, 0xcccc, 0xcccc, 0x3ffb,}, /* 10**-1 */
4439 /* Convert float value X to ASCII string STRING with NDIG digits after
4440 the decimal point. */
4443 e24toasc (x
, string
, ndigs
)
4444 const UEMUSHORT x
[];
4451 etoasc (w
, string
, ndigs
);
4454 /* Convert double value X to ASCII string STRING with NDIG digits after
4455 the decimal point. */
4458 e53toasc (x
, string
, ndigs
)
4459 const UEMUSHORT x
[];
4466 etoasc (w
, string
, ndigs
);
4469 /* Convert double extended value X to ASCII string STRING with NDIG digits
4470 after the decimal point. */
4473 e64toasc (x
, string
, ndigs
)
4474 const UEMUSHORT x
[];
4481 etoasc (w
, string
, ndigs
);
4484 /* Convert 128-bit long double value X to ASCII string STRING with NDIG digits
4485 after the decimal point. */
4488 e113toasc (x
, string
, ndigs
)
4489 const UEMUSHORT x
[];
4496 etoasc (w
, string
, ndigs
);
4500 /* Convert e-type X to ASCII string STRING with NDIGS digits after
4501 the decimal point. */
4503 static char wstring
[80]; /* working storage for ASCII output */
4506 etoasc (x
, string
, ndigs
)
4507 const UEMUSHORT x
[];
4512 UEMUSHORT y
[NI
], t
[NI
], u
[NI
], w
[NI
];
4513 const UEMUSHORT
*p
, *r
, *ten
;
4515 int i
, j
, k
, expon
, rndsav
;
4528 sprintf (wstring
, " NaN ");
4532 rndprc
= NBITS
; /* set to full precision */
4533 emov (x
, y
); /* retain external format */
4534 if (y
[NE
- 1] & 0x8000)
4537 y
[NE
- 1] &= 0x7fff;
4544 ten
= &etens
[NTEN
][0];
4546 /* Test for zero exponent */
4549 for (k
= 0; k
< NE
- 1; k
++)
4552 goto tnzro
; /* denormalized number */
4554 goto isone
; /* valid all zeros */
4558 /* Test for infinity. */
4559 if (y
[NE
- 1] == 0x7fff)
4562 sprintf (wstring
, " -Infinity ");
4564 sprintf (wstring
, " Infinity ");
4568 /* Test for exponent nonzero but significand denormalized.
4569 * This is an error condition.
4571 if ((y
[NE
- 1] != 0) && ((y
[NE
- 2] & 0x8000) == 0))
4573 mtherr ("etoasc", DOMAIN
);
4574 sprintf (wstring
, "NaN");
4578 /* Compare to 1.0 */
4587 { /* Number is greater than 1 */
4588 /* Convert significand to an integer and strip trailing decimal zeros. */
4590 u
[NE
- 1] = EXONE
+ NBITS
- 1;
4592 p
= &etens
[NTEN
- 4][0];
4598 for (j
= 0; j
< NE
- 1; j
++)
4611 /* Rescale from integer significand */
4612 u
[NE
- 1] += y
[NE
- 1] - (unsigned int) (EXONE
+ NBITS
- 1);
4614 /* Find power of 10 */
4618 /* An unordered compare result shouldn't happen here. */
4619 while (ecmp (ten
, u
) <= 0)
4621 if (ecmp (p
, u
) <= 0)
4634 { /* Number is less than 1.0 */
4635 /* Pad significand with trailing decimal zeros. */
4638 while ((y
[NE
- 2] & 0x8000) == 0)
4647 for (i
= 0; i
< NDEC
+ 1; i
++)
4649 if ((w
[NI
- 1] & 0x7) != 0)
4651 /* multiply by 10 */
4664 if (eone
[NE
- 1] <= u
[1])
4676 while (ecmp (eone
, w
) > 0)
4678 if (ecmp (p
, w
) >= 0)
4693 /* Find the first (leading) digit. */
4699 digit
= equot
[NI
- 1];
4700 while ((digit
== 0) && (ecmp (y
, ezero
) != 0))
4708 digit
= equot
[NI
- 1];
4716 /* Examine number of digits requested by caller. */
4734 *s
++ = (char) digit
+ '0';
4737 /* Generate digits after the decimal point. */
4738 for (k
= 0; k
<= ndigs
; k
++)
4740 /* multiply current number by 10, without normalizing */
4747 *s
++ = (char) equot
[NI
- 1] + '0';
4749 digit
= equot
[NI
- 1];
4752 /* round off the ASCII string */
4755 /* Test for critical rounding case in ASCII output. */
4759 if (ecmp (t
, ezero
) != 0)
4760 goto roun
; /* round to nearest */
4762 if ((*(s
- 1) & 1) == 0)
4763 goto doexp
; /* round to even */
4766 /* Round up and propagate carry-outs */
4770 /* Carry out to most significant digit? */
4777 /* Most significant digit carries to 10? */
4785 /* Round up and carry out from less significant digits */
4795 /* Strip trailing zeros, but leave at least one. */
4796 while (ss
[-1] == '0' && ss
[-2] != '.')
4798 sprintf (ss
, "e%d", expon
);
4801 /* copy out the working string */
4804 while (*ss
== ' ') /* strip possible leading space */
4806 while ((*s
++ = *ss
++) != '\0')
4811 /* Convert ASCII string to floating point.
4813 Numeric input is a free format decimal number of any length, with
4814 or without decimal point. Entering E after the number followed by an
4815 integer number causes the second number to be interpreted as a power of
4816 10 to be multiplied by the first number (i.e., "scientific" notation). */
4818 /* Convert ASCII string S to single precision float value Y. */
4829 /* Convert ASCII string S to double precision value Y. */
4836 #if defined(DEC) || defined(IBM)
4848 /* Convert ASCII string S to double extended value Y. */
4858 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
4859 /* Convert ASCII string S to 128-bit long double Y. */
4866 asctoeg (s
, y
, 113);
4870 /* Convert ASCII string S to e type Y. */
4877 asctoeg (s
, y
, NBITS
);
4880 /* Convert ASCII string SS to e type Y, with a specified rounding precision
4881 of OPREC bits. BASE is 16 for C99 hexadecimal floating constants. */
4884 asctoeg (ss
, y
, oprec
)
4889 UEMUSHORT yy
[NI
], xt
[NI
], tt
[NI
];
4890 int esign
, decflg
, sgnflg
, nexp
, exp
, prec
, lost
;
4891 int i
, k
, trail
, c
, rndsav
;
4894 char *sp
, *s
, *lstr
;
4897 /* Copy the input string. */
4898 lstr
= (char *) alloca (strlen (ss
) + 1);
4900 while (*ss
== ' ') /* skip leading spaces */
4904 while ((*sp
++ = *ss
++) != '\0')
4908 if (s
[0] == '0' && (s
[1] == 'x' || s
[1] == 'X'))
4915 rndprc
= NBITS
; /* Set to full precision */
4928 if ((k
>= 0) && (k
< base
))
4930 /* Ignore leading zeros */
4931 if ((prec
== 0) && (decflg
== 0) && (k
== 0))
4933 /* Identify and strip trailing zeros after the decimal point. */
4934 if ((trail
== 0) && (decflg
!= 0))
4937 while (ISDIGIT (*sp
) || (base
== 16 && ISXDIGIT (*sp
)))
4939 /* Check for syntax error */
4941 if ((base
!= 10 || ((c
!= 'e') && (c
!= 'E')))
4942 && (base
!= 16 || ((c
!= 'p') && (c
!= 'P')))
4944 && (c
!= '\n') && (c
!= '\r') && (c
!= ' ')
4946 goto unexpected_char_error
;
4955 /* If enough digits were given to more than fill up the yy register,
4956 continuing until overflow into the high guard word yy[2]
4957 guarantees that there will be a roundoff bit at the top
4958 of the low guard word after normalization. */
4965 nexp
+= 4; /* count digits after decimal point */
4967 eshup1 (yy
); /* multiply current number by 16 */
4975 nexp
+= 1; /* count digits after decimal point */
4977 eshup1 (yy
); /* multiply current number by 10 */
4983 /* Insert the current digit. */
4985 xt
[NI
- 2] = (UEMUSHORT
) k
;
4990 /* Mark any lost non-zero digit. */
4992 /* Count lost digits before the decimal point. */
5014 case '.': /* decimal point */
5016 goto unexpected_char_error
;
5022 goto unexpected_char_error
;
5027 goto unexpected_char_error
;
5040 unexpected_char_error
:
5044 mtherr ("asctoe", DOMAIN
);
5053 /* Exponent interpretation */
5055 /* 0.0eXXX is zero, regardless of XXX. Check for the 0.0. */
5056 for (k
= 0; k
< NI
; k
++)
5067 /* check for + or - */
5075 while (ISDIGIT (*s
))
5084 if ((exp
> MAXDECEXP
) && (base
== 10))
5088 yy
[E
] = 0x7fff; /* infinity */
5091 if ((exp
< MINDECEXP
) && (base
== 10))
5101 /* Base 16 hexadecimal floating constant. */
5102 if ((k
= enormlz (yy
)) > NBITS
)
5107 /* Adjust the exponent. NEXP is the number of hex digits,
5108 EXP is a power of 2. */
5109 lexp
= (EXONE
- 1 + NBITS
) - k
+ yy
[E
] + exp
- nexp
;
5119 /* Pad trailing zeros to minimize power of 10, per IEEE spec. */
5120 while ((nexp
> 0) && (yy
[2] == 0))
5132 if ((k
= enormlz (yy
)) > NBITS
)
5137 lexp
= (EXONE
- 1 + NBITS
) - k
;
5138 emdnorm (yy
, lost
, 0, lexp
, 64);
5141 /* Convert to external format:
5143 Multiply by 10**nexp. If precision is 64 bits,
5144 the maximum relative error incurred in forming 10**n
5145 for 0 <= n <= 324 is 8.2e-20, at 10**180.
5146 For 0 <= n <= 999, the peak relative error is 1.4e-19 at 10**947.
5147 For 0 >= n >= -999, it is -1.55e-19 at 10**-435. */
5162 /* Punt. Can't handle this without 2 divides. */
5163 emovi (etens
[0], tt
);
5176 emul (etens
[i
], xt
, xt
);
5180 while (exp
<= MAXP
);
5199 /* Round and convert directly to the destination type */
5201 lexp
-= EXONE
- 0x3ff;
5203 else if (oprec
== 24 || oprec
== 32)
5204 lexp
-= (EXONE
- 0x7f);
5207 else if (oprec
== 24 || oprec
== 56)
5208 lexp
-= EXONE
- (0x41 << 2);
5211 else if (oprec
== 24)
5212 lexp
-= dec_f
.adjustment
;
5213 else if (oprec
== 56)
5216 lexp
-= dec_g
.adjustment
;
5218 lexp
-= dec_d
.adjustment
;
5221 else if (oprec
== 24)
5222 lexp
-= EXONE
- 0177;
5227 emdnorm (yy
, lost
, 0, lexp
, 64);
5242 toibm (yy
, y
, DFmode
);
5247 toc4x (yy
, y
, HFmode
);
5260 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
5273 /* Return Y = largest integer not greater than X (truncated toward minus
5276 static const UEMUSHORT bmask
[] =
5299 const UEMUSHORT x
[];
5306 emov (x
, f
); /* leave in external format */
5307 expon
= (int) f
[NE
- 1];
5308 e
= (expon
& 0x7fff) - (EXONE
- 1);
5314 /* number of bits to clear out */
5326 /* clear the remaining bits */
5328 /* truncate negatives toward minus infinity */
5331 if ((UEMUSHORT
) expon
& (UEMUSHORT
) 0x8000)
5333 for (i
= 0; i
< NE
- 1; i
++)
5346 /* Return S and EXP such that S * 2^EXP = X and .5 <= S < 1.
5347 For example, 1.1 = 0.55 * 2^1. */
5351 const UEMUSHORT x
[];
5359 /* Handle denormalized numbers properly using long integer exponent. */
5360 li
= (EMULONG
) ((EMUSHORT
) xi
[1]);
5368 *exp
= (int) (li
- 0x3ffe);
5372 /* Return e type Y = X * 2^PWR2. */
5376 const UEMUSHORT x
[];
5388 emdnorm (xi
, i
, i
, li
, !ROUND_TOWARDS_ZERO
);
5394 /* C = remainder after dividing B by A, all e type values.
5395 Least significant integer quotient bits left in EQUOT. */
5399 const UEMUSHORT a
[], b
[];
5402 UEMUSHORT den
[NI
], num
[NI
];
5406 || (ecmp (a
, ezero
) == 0)
5414 if (ecmp (a
, ezero
) == 0)
5416 mtherr ("eremain", SING
);
5422 eiremain (den
, num
);
5423 /* Sign of remainder = sign of quotient */
5432 /* Return quotient of exploded e-types NUM / DEN in EQUOT,
5433 remainder in NUM. */
5437 UEMUSHORT den
[], num
[];
5443 ld
-= enormlz (den
);
5445 ln
-= enormlz (num
);
5449 if (ecmpm (den
, num
) <= 0)
5461 emdnorm (num
, 0, 0, ln
, 0);
5464 /* Report an error condition CODE encountered in function NAME.
5466 Mnemonic Value Significance
5468 DOMAIN 1 argument domain error
5469 SING 2 function singularity
5470 OVERFLOW 3 overflow range error
5471 UNDERFLOW 4 underflow range error
5472 TLOSS 5 total loss of precision
5473 PLOSS 6 partial loss of precision
5474 INVALID 7 NaN - producing operation
5475 EDOM 33 Unix domain error code
5476 ERANGE 34 Unix range error code
5478 The order of appearance of the following messages is bound to the
5479 error codes defined above. */
5489 /* The string passed by the calling program is supposed to be the
5490 name of the function in which the error occurred.
5491 The code argument selects which error message string will be printed. */
5493 if (strcmp (name
, "esub") == 0)
5494 name
= "subtraction";
5495 else if (strcmp (name
, "ediv") == 0)
5497 else if (strcmp (name
, "emul") == 0)
5498 name
= "multiplication";
5499 else if (strcmp (name
, "enormlz") == 0)
5500 name
= "normalization";
5501 else if (strcmp (name
, "etoasc") == 0)
5502 name
= "conversion to text";
5503 else if (strcmp (name
, "asctoe") == 0)
5505 else if (strcmp (name
, "eremain") == 0)
5507 else if (strcmp (name
, "esqrt") == 0)
5508 name
= "square root";
5513 case DOMAIN
: warning ("%s: argument domain error" , name
); break;
5514 case SING
: warning ("%s: function singularity" , name
); break;
5515 case OVERFLOW
: warning ("%s: overflow range error" , name
); break;
5516 case UNDERFLOW
: warning ("%s: underflow range error" , name
); break;
5517 case TLOSS
: warning ("%s: total loss of precision" , name
); break;
5518 case PLOSS
: warning ("%s: partial loss of precision", name
); break;
5519 case INVALID
: warning ("%s: NaN - producing operation", name
); break;
5524 /* Set global error message word */
5529 /* Convert DEC double precision D to e type E. */
5537 ieeetoe (d
, e
, &dec_g
);
5539 ieeetoe (d
, e
, &dec_d
);
5542 /* Convert e type X to DEC double precision D. */
5552 const struct ieee_format
*fmt
;
5560 /* Adjust exponent for offsets. */
5561 exp
= (EMULONG
) xi
[E
] - fmt
->adjustment
;
5562 /* Round off to nearest or even. */
5564 rndprc
= fmt
->precision
;
5565 emdnorm (xi
, 0, 0, exp
, !ROUND_TOWARDS_ZERO
);
5570 /* Convert exploded e-type X, that has already been rounded to
5571 56-bit precision, to DEC format double Y. */
5578 toieee (x
, y
, &dec_g
);
5580 toieee (x
, y
, &dec_d
);
5585 /* Convert IBM single/double precision to e type. */
5591 enum machine_mode mode
;
5596 ecleaz (y
); /* start with a zero */
5597 p
= y
; /* point to our number */
5598 r
= *d
; /* get IBM exponent word */
5599 if (*d
& (unsigned int) 0x8000)
5600 *p
= 0xffff; /* fill in our sign */
5601 ++p
; /* bump pointer to our exponent word */
5602 r
&= 0x7f00; /* strip the sign bit */
5603 r
>>= 6; /* shift exponent word down 6 bits */
5604 /* in fact shift by 8 right and 2 left */
5605 r
+= EXONE
- (0x41 << 2); /* subtract IBM exponent offset */
5606 /* add our e type exponent offset */
5607 *p
++ = r
; /* to form our exponent */
5609 *p
++ = *d
++ & 0xff; /* now do the high order mantissa */
5610 /* strip off the IBM exponent and sign bits */
5611 if (mode
!= SFmode
) /* there are only 2 words in SFmode */
5613 *p
++ = *d
++; /* fill in the rest of our mantissa */
5618 if (y
[M
] == 0 && y
[M
+1] == 0 && y
[M
+2] == 0 && y
[M
+3] == 0)
5621 y
[E
] -= 5 + enormlz (y
); /* now normalise the mantissa */
5622 /* handle change in RADIX */
5628 /* Convert e type to IBM single/double precision. */
5634 enum machine_mode mode
;
5641 exp
= (EMULONG
) xi
[E
] - (EXONE
- (0x41 << 2)); /* adjust exponent for offsets */
5642 /* round off to nearest or even */
5645 emdnorm (xi
, 0, 0, exp
, !ROUND_TOWARDS_ZERO
);
5647 toibm (xi
, d
, mode
);
5653 enum machine_mode mode
;
5706 /* Convert C4X single/double precision to e type. */
5712 enum machine_mode mode
;
5730 /* Short-circuit the zero case. */
5731 if ((dn
[0] == 0x8000)
5732 && (dn
[1] == 0x0000)
5733 && ((mode
== QFmode
) || ((dn
[2] == 0x0000) && (dn
[3] == 0x0000))))
5744 ecleaz (y
); /* start with a zero */
5745 r
= dn
[0]; /* get sign/exponent part */
5746 if (r
& (unsigned int) 0x0080)
5748 y
[0] = 0xffff; /* fill in our sign */
5754 r
>>= 8; /* Shift exponent word down 8 bits. */
5755 if (r
& 0x80) /* Make the exponent negative if it is. */
5756 r
= r
| (~0 & ~0xff);
5760 /* Now do the high order mantissa. We don't "or" on the high bit
5761 because it is 2 (not 1) and is handled a little differently
5763 y
[M
] = dn
[0] & 0x7f;
5766 if (mode
!= QFmode
) /* There are only 2 words in QFmode. */
5768 y
[M
+2] = dn
[2]; /* Fill in the rest of our mantissa. */
5776 /* Now do the two's complement on the data. */
5778 carry
= 1; /* Initially add 1 for the two's complement. */
5779 for (i
=size
+ M
; i
> M
; i
--)
5781 if (carry
&& (y
[i
] == 0x0000))
5782 /* We overflowed into the next word, carry is the same. */
5783 y
[i
] = carry
? 0x0000 : 0xffff;
5786 /* No overflow, just invert and add carry. */
5787 y
[i
] = ((~y
[i
]) + carry
) & 0xffff;
5802 /* Add our e type exponent offset to form our exponent. */
5806 /* Now do the high order mantissa strip off the exponent and sign
5807 bits and add the high 1 bit. */
5808 y
[M
] = (dn
[0] & 0x7f) | 0x80;
5811 if (mode
!= QFmode
) /* There are only 2 words in QFmode. */
5813 y
[M
+2] = dn
[2]; /* Fill in the rest of our mantissa. */
5823 /* Convert e type to C4X single/double precision. */
5829 enum machine_mode mode
;
5837 /* Adjust exponent for offsets. */
5838 exp
= (EMULONG
) xi
[E
] - (EXONE
- 0x7f);
5840 /* Round off to nearest or even. */
5842 rndprc
= mode
== QFmode
? 24 : 32;
5843 emdnorm (xi
, 0, 0, exp
, !ROUND_TOWARDS_ZERO
);
5845 toc4x (xi
, d
, mode
);
5851 enum machine_mode mode
;
5857 /* Short-circuit the zero case */
5858 if ((x
[0] == 0) /* Zero exponent and sign */
5860 && (x
[M
] == 0) /* The rest is for zero mantissa */
5862 /* Only check for double if necessary */
5863 && ((mode
== QFmode
) || ((x
[M
+2] == 0) && (x
[M
+3] == 0))))
5865 /* We have a zero. Put it into the output and return. */
5878 /* Negative number require a two's complement conversion of the
5884 i
= ((int) x
[1]) - 0x7f;
5886 /* Now add 1 to the inverted data to do the two's complement. */
5895 x
[v
] = carry
? 0x0000 : 0xffff;
5898 x
[v
] = ((~x
[v
]) + carry
) & 0xffff;
5904 /* The following is a special case. The C4X negative float requires
5905 a zero in the high bit (because the format is (2 - x) x 2^m), so
5906 if a one is in that bit, we have to shift left one to get rid
5907 of it. This only occurs if the number is -1 x 2^m. */
5908 if (x
[M
+1] & 0x8000)
5910 /* This is the case of -1 x 2^m, we have to rid ourselves of the
5911 high sign bit and shift the exponent. */
5917 i
= ((int) x
[1]) - 0x7f;
5919 if ((i
< -128) || (i
> 127))
5927 y
[3] = (y
[1] << 8) | ((y
[2] >> 8) & 0xff);
5928 y
[2] = (y
[0] << 8) | ((y
[1] >> 8) & 0xff);
5936 y
[0] |= ((i
& 0xff) << 8);
5940 y
[0] |= x
[M
] & 0x7f;
5946 y
[3] = (y
[1] << 8) | ((y
[2] >> 8) & 0xff);
5947 y
[2] = (y
[0] << 8) | ((y
[1] >> 8) & 0xff);
5952 /* Output a binary NaN bit pattern in the target machine's format. */
5954 /* If special NaN bit patterns are required, define them in tm.h
5955 as arrays of unsigned 16-bit shorts. Otherwise, use the default
5961 static const UEMUSHORT TFbignan
[8] =
5962 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5963 static const UEMUSHORT TFlittlenan
[8] = {0, 0, 0, 0, 0, 0, 0x8000, 0xffff};
5971 static const UEMUSHORT XFbignan
[6] =
5972 {0x7fff, 0xffff, 0xffff, 0xffff, 0xffff, 0xffff};
5973 static const UEMUSHORT XFlittlenan
[6] = {0, 0, 0, 0xc000, 0xffff, 0};
5981 static const UEMUSHORT DFbignan
[4] = {0x7fff, 0xffff, 0xffff, 0xffff};
5982 static const UEMUSHORT DFlittlenan
[4] = {0, 0, 0, 0xfff8};
5990 static const UEMUSHORT SFbignan
[2] = {0x7fff, 0xffff};
5991 static const UEMUSHORT SFlittlenan
[2] = {0, 0xffc0};
5998 make_nan (nan
, sign
, mode
)
6001 enum machine_mode mode
;
6007 size
= GET_MODE_BITSIZE (mode
);
6008 if (LARGEST_EXPONENT_IS_NORMAL (size
))
6010 warning ("%d-bit floats cannot hold NaNs", size
);
6011 saturate (nan
, sign
, size
, 0);
6016 /* Possibly the `reserved operand' patterns on a VAX can be
6017 used like NaN's, but probably not in the same way as IEEE. */
6018 #if !defined(DEC) && !defined(IBM) && !defined(C4X)
6020 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)
6022 if (REAL_WORDS_BIG_ENDIAN
)
6032 if (REAL_WORDS_BIG_ENDIAN
)
6040 if (REAL_WORDS_BIG_ENDIAN
)
6049 if (REAL_WORDS_BIG_ENDIAN
)
6059 if (REAL_WORDS_BIG_ENDIAN
)
6060 *nan
++ = (sign
<< 15) | (*p
++ & 0x7fff);
6063 if (! REAL_WORDS_BIG_ENDIAN
)
6064 *nan
= (sign
<< 15) | (*p
& 0x7fff);
6069 /* Create a saturation value for a SIZE-bit float, assuming that
6070 LARGEST_EXPONENT_IS_NORMAL (SIZE).
6072 If SIGN is true, fill X with the most negative value, otherwise fill
6073 it with the most positive value. WARN is true if the function should
6074 warn about overflow. */
6077 saturate (x
, sign
, size
, warn
)
6079 int sign
, size
, warn
;
6083 if (warn
&& extra_warnings
)
6084 warning ("value exceeds the range of a %d-bit float", size
);
6086 /* Create the most negative value. */
6087 for (i
= 0; i
< size
/ EMUSHORT_SIZE
; i
++)
6090 /* Make it positive, if necessary. */
6092 x
[REAL_WORDS_BIG_ENDIAN
? 0 : i
- 1] = 0x7fff;
6096 /* This is the inverse of the function `etarsingle' invoked by
6097 REAL_VALUE_TO_TARGET_SINGLE. */
6100 ereal_unto_float (f
)
6107 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6108 This is the inverse operation to what the function `endian' does. */
6109 if (REAL_WORDS_BIG_ENDIAN
)
6111 s
[0] = (UEMUSHORT
) (f
>> 16);
6112 s
[1] = (UEMUSHORT
) f
;
6116 s
[0] = (UEMUSHORT
) f
;
6117 s
[1] = (UEMUSHORT
) (f
>> 16);
6119 /* Convert and promote the target float to E-type. */
6121 /* Output E-type to REAL_VALUE_TYPE. */
6127 /* This is the inverse of the function `etardouble' invoked by
6128 REAL_VALUE_TO_TARGET_DOUBLE. */
6131 ereal_unto_double (d
)
6138 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
6139 if (REAL_WORDS_BIG_ENDIAN
)
6141 s
[0] = (UEMUSHORT
) (d
[0] >> 16);
6142 s
[1] = (UEMUSHORT
) d
[0];
6143 s
[2] = (UEMUSHORT
) (d
[1] >> 16);
6144 s
[3] = (UEMUSHORT
) d
[1];
6148 /* Target float words are little-endian. */
6149 s
[0] = (UEMUSHORT
) d
[0];
6150 s
[1] = (UEMUSHORT
) (d
[0] >> 16);
6151 s
[2] = (UEMUSHORT
) d
[1];
6152 s
[3] = (UEMUSHORT
) (d
[1] >> 16);
6154 /* Convert target double to E-type. */
6156 /* Output E-type to REAL_VALUE_TYPE. */
6162 /* Convert an SFmode target `float' value to a REAL_VALUE_TYPE.
6163 This is somewhat like ereal_unto_float, but the input types
6164 for these are different. */
6167 ereal_from_float (f
)
6174 /* Convert 32 bit integer to array of 16 bit pieces in target machine order.
6175 This is the inverse operation to what the function `endian' does. */
6176 if (REAL_WORDS_BIG_ENDIAN
)
6178 s
[0] = (UEMUSHORT
) (f
>> 16);
6179 s
[1] = (UEMUSHORT
) f
;
6183 s
[0] = (UEMUSHORT
) f
;
6184 s
[1] = (UEMUSHORT
) (f
>> 16);
6186 /* Convert and promote the target float to E-type. */
6188 /* Output E-type to REAL_VALUE_TYPE. */
6194 /* Convert a DFmode target `double' value to a REAL_VALUE_TYPE.
6195 This is somewhat like ereal_unto_double, but the input types
6196 for these are different.
6198 The DFmode is stored as an array of HOST_WIDE_INT in the target's
6199 data format, with no holes in the bit packing. The first element
6200 of the input array holds the bits that would come first in the
6201 target computer's memory. */
6204 ereal_from_double (d
)
6211 /* Convert array of HOST_WIDE_INT to equivalent array of 16-bit pieces. */
6212 if (REAL_WORDS_BIG_ENDIAN
)
6214 #if HOST_BITS_PER_WIDE_INT == 32
6215 s
[0] = (UEMUSHORT
) (d
[0] >> 16);
6216 s
[1] = (UEMUSHORT
) d
[0];
6217 s
[2] = (UEMUSHORT
) (d
[1] >> 16);
6218 s
[3] = (UEMUSHORT
) d
[1];
6220 /* In this case the entire target double is contained in the
6221 first array element. The second element of the input is
6223 s
[0] = (UEMUSHORT
) (d
[0] >> 48);
6224 s
[1] = (UEMUSHORT
) (d
[0] >> 32);
6225 s
[2] = (UEMUSHORT
) (d
[0] >> 16);
6226 s
[3] = (UEMUSHORT
) d
[0];
6231 /* Target float words are little-endian. */
6232 s
[0] = (UEMUSHORT
) d
[0];
6233 s
[1] = (UEMUSHORT
) (d
[0] >> 16);
6234 #if HOST_BITS_PER_WIDE_INT == 32
6235 s
[2] = (UEMUSHORT
) d
[1];
6236 s
[3] = (UEMUSHORT
) (d
[1] >> 16);
6238 s
[2] = (UEMUSHORT
) (d
[0] >> 32);
6239 s
[3] = (UEMUSHORT
) (d
[0] >> 48);
6242 /* Convert target double to E-type. */
6244 /* Output E-type to REAL_VALUE_TYPE. */
6251 /* Convert target computer unsigned 64-bit integer to e-type.
6252 The endian-ness of DImode follows the convention for integers,
6253 so we use WORDS_BIG_ENDIAN here, not REAL_WORDS_BIG_ENDIAN. */
6257 const UEMUSHORT
*di
; /* Address of the 64-bit int. */
6264 if (WORDS_BIG_ENDIAN
)
6266 for (k
= M
; k
< M
+ 4; k
++)
6271 for (k
= M
+ 3; k
>= M
; k
--)
6274 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
6275 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
6276 ecleaz (yi
); /* it was zero */
6278 yi
[E
] -= (UEMUSHORT
) k
;/* subtract shift count from exponent */
6282 /* Convert target computer signed 64-bit integer to e-type. */
6286 const UEMUSHORT
*di
; /* Address of the 64-bit int. */
6289 unsigned EMULONG acc
;
6295 if (WORDS_BIG_ENDIAN
)
6297 for (k
= M
; k
< M
+ 4; k
++)
6302 for (k
= M
+ 3; k
>= M
; k
--)
6305 /* Take absolute value */
6311 for (k
= M
+ 3; k
>= M
; k
--)
6313 acc
= (unsigned EMULONG
) (~yi
[k
] & 0xffff) + carry
;
6320 yi
[E
] = EXONE
+ 47; /* exponent if normalize shift count were 0 */
6321 if ((k
= enormlz (yi
)) > NBITS
)/* normalize the significand */
6322 ecleaz (yi
); /* it was zero */
6324 yi
[E
] -= (UEMUSHORT
) k
;/* subtract shift count from exponent */
6331 /* Convert e-type to unsigned 64-bit int. */
6347 k
= (int) xi
[E
] - (EXONE
- 1);
6350 for (j
= 0; j
< 4; j
++)
6356 for (j
= 0; j
< 4; j
++)
6359 warning ("overflow on truncation to integer");
6364 /* Shift more than 16 bits: first shift up k-16 mod 16,
6365 then shift up by 16's. */
6366 j
= k
- ((k
>> 4) << 4);
6370 if (WORDS_BIG_ENDIAN
)
6381 if (WORDS_BIG_ENDIAN
)
6386 while ((k
-= 16) > 0);
6390 /* shift not more than 16 bits */
6395 if (WORDS_BIG_ENDIAN
)
6414 /* Convert e-type to signed 64-bit int. */
6421 unsigned EMULONG acc
;
6428 k
= (int) xi
[E
] - (EXONE
- 1);
6431 for (j
= 0; j
< 4; j
++)
6437 for (j
= 0; j
< 4; j
++)
6440 warning ("overflow on truncation to integer");
6446 /* Shift more than 16 bits: first shift up k-16 mod 16,
6447 then shift up by 16's. */
6448 j
= k
- ((k
>> 4) << 4);
6452 if (WORDS_BIG_ENDIAN
)
6463 if (WORDS_BIG_ENDIAN
)
6468 while ((k
-= 16) > 0);
6472 /* shift not more than 16 bits */
6475 if (WORDS_BIG_ENDIAN
)
6491 /* Negate if negative */
6495 if (WORDS_BIG_ENDIAN
)
6497 for (k
= 0; k
< 4; k
++)
6499 acc
= (unsigned EMULONG
) (~(*isave
) & 0xffff) + carry
;
6500 if (WORDS_BIG_ENDIAN
)
6512 /* Longhand square root routine. */
6515 static int esqinited
= 0;
6516 static unsigned short sqrndbit
[NI
];
6523 UEMUSHORT temp
[NI
], num
[NI
], sq
[NI
], xx
[NI
];
6525 int i
, j
, k
, n
, nlups
;
6530 sqrndbit
[NI
- 2] = 1;
6533 /* Check for arg <= 0 */
6534 i
= ecmp (x
, ezero
);
6539 mtherr ("esqrt", DOMAIN
);
6555 /* Bring in the arg and renormalize if it is denormal. */
6557 m
= (EMULONG
) xx
[1]; /* local long word exponent */
6561 /* Divide exponent by 2 */
6563 exp
= (unsigned short) ((m
/ 2) + 0x3ffe);
6565 /* Adjust if exponent odd */
6575 n
= 8; /* get 8 bits of result per inner loop */
6581 /* bring in next word of arg */
6583 num
[NI
- 1] = xx
[j
+ 3];
6584 /* Do additional bit on last outer loop, for roundoff. */
6587 for (i
= 0; i
< n
; i
++)
6589 /* Next 2 bits of arg */
6592 /* Shift up answer */
6594 /* Make trial divisor */
6595 for (k
= 0; k
< NI
; k
++)
6598 eaddm (sqrndbit
, temp
);
6599 /* Subtract and insert answer bit if it goes in */
6600 if (ecmpm (temp
, num
) <= 0)
6610 /* Adjust for extra, roundoff loop done. */
6611 exp
+= (NBITS
- 1) - rndprc
;
6613 /* Sticky bit = 1 if the remainder is nonzero. */
6615 for (i
= 3; i
< NI
; i
++)
6618 /* Renormalize and round off. */
6619 emdnorm (sq
, k
, 0, exp
, !ROUND_TOWARDS_ZERO
);
6624 /* Return the binary precision of the significand for a given
6625 floating point mode. The mode can hold an integer value
6626 that many bits wide, without losing any bits. */
6629 significand_size (mode
)
6630 enum machine_mode mode
;
6633 /* Don't test the modes, but their sizes, lest this
6634 code won't work for BITS_PER_UNIT != 8 . */
6636 switch (GET_MODE_BITSIZE (mode
))
6657 #if (INTEL_EXTENDED_IEEE_FORMAT == 0)