1 /* real.c - software floating point emulation.
2 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 1999,
3 2000, 2002, 2003, 2004 Free Software Foundation, Inc.
4 Contributed by Stephen L. Moshier (moshier@world.std.com).
5 Re-written by Richard Henderson <rth@redhat.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
26 #include "coretypes.h"
33 /* The floating point model used internally is not exactly IEEE 754
34 compliant, and close to the description in the ISO C99 standard,
35 section 5.2.4.2.2 Characteristics of floating types.
39 x = s * b^e * \sum_{k=1}^p f_k * b^{-k}
43 b = base or radix, here always 2
45 p = precision (the number of base-b digits in the significand)
46 f_k = the digits of the significand.
48 We differ from typical IEEE 754 encodings in that the entire
49 significand is fractional. Normalized significands are in the
52 A requirement of the model is that P be larger than the largest
53 supported target floating-point type by at least 2 bits. This gives
54 us proper rounding when we truncate to the target type. In addition,
55 E must be large enough to hold the smallest supported denormal number
58 Both of these requirements are easily satisfied. The largest target
59 significand is 113 bits; we store at least 160. The smallest
60 denormal number fits in 17 exponent bits; we store 29.
62 Note that the decimal string conversion routines are sensitive to
63 rounding errors. Since the raw arithmetic routines do not themselves
64 have guard digits or rounding, the computation of 10**exp can
65 accumulate more than a few digits of error. The previous incarnation
66 of real.c successfully used a 144-bit fraction; given the current
67 layout of REAL_VALUE_TYPE we're forced to expand to at least 160 bits.
69 Target floating point models that use base 16 instead of base 2
70 (i.e. IBM 370), are handled during round_for_format, in which we
71 canonicalize the exponent to be a multiple of 4 (log2(16)), and
72 adjust the significand to match. */
75 /* Used to classify two numbers simultaneously. */
76 #define CLASS2(A, B) ((A) << 2 | (B))
78 #if HOST_BITS_PER_LONG != 64 && HOST_BITS_PER_LONG != 32
79 #error "Some constant folding done by hand to avoid shift count warnings"
82 static void get_zero (REAL_VALUE_TYPE
*, int);
83 static void get_canonical_qnan (REAL_VALUE_TYPE
*, int);
84 static void get_canonical_snan (REAL_VALUE_TYPE
*, int);
85 static void get_inf (REAL_VALUE_TYPE
*, int);
86 static bool sticky_rshift_significand (REAL_VALUE_TYPE
*,
87 const REAL_VALUE_TYPE
*, unsigned int);
88 static void rshift_significand (REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*,
90 static void lshift_significand (REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*,
92 static void lshift_significand_1 (REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*);
93 static bool add_significands (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*,
94 const REAL_VALUE_TYPE
*);
95 static bool sub_significands (REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*,
96 const REAL_VALUE_TYPE
*, int);
97 static void neg_significand (REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*);
98 static int cmp_significands (const REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*);
99 static int cmp_significand_0 (const REAL_VALUE_TYPE
*);
100 static void set_significand_bit (REAL_VALUE_TYPE
*, unsigned int);
101 static void clear_significand_bit (REAL_VALUE_TYPE
*, unsigned int);
102 static bool test_significand_bit (REAL_VALUE_TYPE
*, unsigned int);
103 static void clear_significand_below (REAL_VALUE_TYPE
*, unsigned int);
104 static bool div_significands (REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*,
105 const REAL_VALUE_TYPE
*);
106 static void normalize (REAL_VALUE_TYPE
*);
108 static bool do_add (REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*,
109 const REAL_VALUE_TYPE
*, int);
110 static bool do_multiply (REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*,
111 const REAL_VALUE_TYPE
*);
112 static bool do_divide (REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*,
113 const REAL_VALUE_TYPE
*);
114 static int do_compare (const REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*, int);
115 static void do_fix_trunc (REAL_VALUE_TYPE
*, const REAL_VALUE_TYPE
*);
117 static unsigned long rtd_divmod (REAL_VALUE_TYPE
*, REAL_VALUE_TYPE
*);
119 static const REAL_VALUE_TYPE
* ten_to_ptwo (int);
120 static const REAL_VALUE_TYPE
* ten_to_mptwo (int);
121 static const REAL_VALUE_TYPE
* real_digit (int);
122 static void times_pten (REAL_VALUE_TYPE
*, int);
124 static void round_for_format (const struct real_format
*, REAL_VALUE_TYPE
*);
126 /* Initialize R with a positive zero. */
129 get_zero (REAL_VALUE_TYPE
*r
, int sign
)
131 memset (r
, 0, sizeof (*r
));
135 /* Initialize R with the canonical quiet NaN. */
138 get_canonical_qnan (REAL_VALUE_TYPE
*r
, int sign
)
140 memset (r
, 0, sizeof (*r
));
147 get_canonical_snan (REAL_VALUE_TYPE
*r
, int sign
)
149 memset (r
, 0, sizeof (*r
));
157 get_inf (REAL_VALUE_TYPE
*r
, int sign
)
159 memset (r
, 0, sizeof (*r
));
165 /* Right-shift the significand of A by N bits; put the result in the
166 significand of R. If any one bits are shifted out, return true. */
169 sticky_rshift_significand (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
,
172 unsigned long sticky
= 0;
173 unsigned int i
, ofs
= 0;
175 if (n
>= HOST_BITS_PER_LONG
)
177 for (i
= 0, ofs
= n
/ HOST_BITS_PER_LONG
; i
< ofs
; ++i
)
179 n
&= HOST_BITS_PER_LONG
- 1;
184 sticky
|= a
->sig
[ofs
] & (((unsigned long)1 << n
) - 1);
185 for (i
= 0; i
< SIGSZ
; ++i
)
188 = (((ofs
+ i
>= SIGSZ
? 0 : a
->sig
[ofs
+ i
]) >> n
)
189 | ((ofs
+ i
+ 1 >= SIGSZ
? 0 : a
->sig
[ofs
+ i
+ 1])
190 << (HOST_BITS_PER_LONG
- n
)));
195 for (i
= 0; ofs
+ i
< SIGSZ
; ++i
)
196 r
->sig
[i
] = a
->sig
[ofs
+ i
];
197 for (; i
< SIGSZ
; ++i
)
204 /* Right-shift the significand of A by N bits; put the result in the
208 rshift_significand (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
,
211 unsigned int i
, ofs
= n
/ HOST_BITS_PER_LONG
;
213 n
&= HOST_BITS_PER_LONG
- 1;
216 for (i
= 0; i
< SIGSZ
; ++i
)
219 = (((ofs
+ i
>= SIGSZ
? 0 : a
->sig
[ofs
+ i
]) >> n
)
220 | ((ofs
+ i
+ 1 >= SIGSZ
? 0 : a
->sig
[ofs
+ i
+ 1])
221 << (HOST_BITS_PER_LONG
- n
)));
226 for (i
= 0; ofs
+ i
< SIGSZ
; ++i
)
227 r
->sig
[i
] = a
->sig
[ofs
+ i
];
228 for (; i
< SIGSZ
; ++i
)
233 /* Left-shift the significand of A by N bits; put the result in the
237 lshift_significand (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
,
240 unsigned int i
, ofs
= n
/ HOST_BITS_PER_LONG
;
242 n
&= HOST_BITS_PER_LONG
- 1;
245 for (i
= 0; ofs
+ i
< SIGSZ
; ++i
)
246 r
->sig
[SIGSZ
-1-i
] = a
->sig
[SIGSZ
-1-i
-ofs
];
247 for (; i
< SIGSZ
; ++i
)
248 r
->sig
[SIGSZ
-1-i
] = 0;
251 for (i
= 0; i
< SIGSZ
; ++i
)
254 = (((ofs
+ i
>= SIGSZ
? 0 : a
->sig
[SIGSZ
-1-i
-ofs
]) << n
)
255 | ((ofs
+ i
+ 1 >= SIGSZ
? 0 : a
->sig
[SIGSZ
-1-i
-ofs
-1])
256 >> (HOST_BITS_PER_LONG
- n
)));
260 /* Likewise, but N is specialized to 1. */
263 lshift_significand_1 (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
)
267 for (i
= SIGSZ
- 1; i
> 0; --i
)
268 r
->sig
[i
] = (a
->sig
[i
] << 1) | (a
->sig
[i
-1] >> (HOST_BITS_PER_LONG
- 1));
269 r
->sig
[0] = a
->sig
[0] << 1;
272 /* Add the significands of A and B, placing the result in R. Return
273 true if there was carry out of the most significant word. */
276 add_significands (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
,
277 const REAL_VALUE_TYPE
*b
)
282 for (i
= 0; i
< SIGSZ
; ++i
)
284 unsigned long ai
= a
->sig
[i
];
285 unsigned long ri
= ai
+ b
->sig
[i
];
301 /* Subtract the significands of A and B, placing the result in R. CARRY is
302 true if there's a borrow incoming to the least significant word.
303 Return true if there was borrow out of the most significant word. */
306 sub_significands (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
,
307 const REAL_VALUE_TYPE
*b
, int carry
)
311 for (i
= 0; i
< SIGSZ
; ++i
)
313 unsigned long ai
= a
->sig
[i
];
314 unsigned long ri
= ai
- b
->sig
[i
];
330 /* Negate the significand A, placing the result in R. */
333 neg_significand (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
)
338 for (i
= 0; i
< SIGSZ
; ++i
)
340 unsigned long ri
, ai
= a
->sig
[i
];
359 /* Compare significands. Return tri-state vs zero. */
362 cmp_significands (const REAL_VALUE_TYPE
*a
, const REAL_VALUE_TYPE
*b
)
366 for (i
= SIGSZ
- 1; i
>= 0; --i
)
368 unsigned long ai
= a
->sig
[i
];
369 unsigned long bi
= b
->sig
[i
];
380 /* Return true if A is nonzero. */
383 cmp_significand_0 (const REAL_VALUE_TYPE
*a
)
387 for (i
= SIGSZ
- 1; i
>= 0; --i
)
394 /* Set bit N of the significand of R. */
397 set_significand_bit (REAL_VALUE_TYPE
*r
, unsigned int n
)
399 r
->sig
[n
/ HOST_BITS_PER_LONG
]
400 |= (unsigned long)1 << (n
% HOST_BITS_PER_LONG
);
403 /* Clear bit N of the significand of R. */
406 clear_significand_bit (REAL_VALUE_TYPE
*r
, unsigned int n
)
408 r
->sig
[n
/ HOST_BITS_PER_LONG
]
409 &= ~((unsigned long)1 << (n
% HOST_BITS_PER_LONG
));
412 /* Test bit N of the significand of R. */
415 test_significand_bit (REAL_VALUE_TYPE
*r
, unsigned int n
)
417 /* ??? Compiler bug here if we return this expression directly.
418 The conversion to bool strips the "&1" and we wind up testing
419 e.g. 2 != 0 -> true. Seen in gcc version 3.2 20020520. */
420 int t
= (r
->sig
[n
/ HOST_BITS_PER_LONG
] >> (n
% HOST_BITS_PER_LONG
)) & 1;
424 /* Clear bits 0..N-1 of the significand of R. */
427 clear_significand_below (REAL_VALUE_TYPE
*r
, unsigned int n
)
429 int i
, w
= n
/ HOST_BITS_PER_LONG
;
431 for (i
= 0; i
< w
; ++i
)
434 r
->sig
[w
] &= ~(((unsigned long)1 << (n
% HOST_BITS_PER_LONG
)) - 1);
437 /* Divide the significands of A and B, placing the result in R. Return
438 true if the division was inexact. */
441 div_significands (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
,
442 const REAL_VALUE_TYPE
*b
)
445 int i
, bit
= SIGNIFICAND_BITS
- 1;
446 unsigned long msb
, inexact
;
449 memset (r
->sig
, 0, sizeof (r
->sig
));
455 msb
= u
.sig
[SIGSZ
-1] & SIG_MSB
;
456 lshift_significand_1 (&u
, &u
);
458 if (msb
|| cmp_significands (&u
, b
) >= 0)
460 sub_significands (&u
, &u
, b
, 0);
461 set_significand_bit (r
, bit
);
466 for (i
= 0, inexact
= 0; i
< SIGSZ
; i
++)
472 /* Adjust the exponent and significand of R such that the most
473 significant bit is set. We underflow to zero and overflow to
474 infinity here, without denormals. (The intermediate representation
475 exponent is large enough to handle target denormals normalized.) */
478 normalize (REAL_VALUE_TYPE
*r
)
483 /* Find the first word that is nonzero. */
484 for (i
= SIGSZ
- 1; i
>= 0; i
--)
486 shift
+= HOST_BITS_PER_LONG
;
490 /* Zero significand flushes to zero. */
498 /* Find the first bit that is nonzero. */
500 if (r
->sig
[i
] & ((unsigned long)1 << (HOST_BITS_PER_LONG
- 1 - j
)))
506 exp
= REAL_EXP (r
) - shift
;
508 get_inf (r
, r
->sign
);
509 else if (exp
< -MAX_EXP
)
510 get_zero (r
, r
->sign
);
513 SET_REAL_EXP (r
, exp
);
514 lshift_significand (r
, r
, shift
);
519 /* Calculate R = A + (SUBTRACT_P ? -B : B). Return true if the
520 result may be inexact due to a loss of precision. */
523 do_add (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
,
524 const REAL_VALUE_TYPE
*b
, int subtract_p
)
528 bool inexact
= false;
530 /* Determine if we need to add or subtract. */
532 subtract_p
= (sign
^ b
->sign
) ^ subtract_p
;
534 switch (CLASS2 (a
->class, b
->class))
536 case CLASS2 (rvc_zero
, rvc_zero
):
537 /* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0. */
538 get_zero (r
, sign
& !subtract_p
);
541 case CLASS2 (rvc_zero
, rvc_normal
):
542 case CLASS2 (rvc_zero
, rvc_inf
):
543 case CLASS2 (rvc_zero
, rvc_nan
):
545 case CLASS2 (rvc_normal
, rvc_nan
):
546 case CLASS2 (rvc_inf
, rvc_nan
):
547 case CLASS2 (rvc_nan
, rvc_nan
):
548 /* ANY + NaN = NaN. */
549 case CLASS2 (rvc_normal
, rvc_inf
):
552 r
->sign
= sign
^ subtract_p
;
555 case CLASS2 (rvc_normal
, rvc_zero
):
556 case CLASS2 (rvc_inf
, rvc_zero
):
557 case CLASS2 (rvc_nan
, rvc_zero
):
559 case CLASS2 (rvc_nan
, rvc_normal
):
560 case CLASS2 (rvc_nan
, rvc_inf
):
561 /* NaN + ANY = NaN. */
562 case CLASS2 (rvc_inf
, rvc_normal
):
567 case CLASS2 (rvc_inf
, rvc_inf
):
569 /* Inf - Inf = NaN. */
570 get_canonical_qnan (r
, 0);
572 /* Inf + Inf = Inf. */
576 case CLASS2 (rvc_normal
, rvc_normal
):
583 /* Swap the arguments such that A has the larger exponent. */
584 dexp
= REAL_EXP (a
) - REAL_EXP (b
);
587 const REAL_VALUE_TYPE
*t
;
594 /* If the exponents are not identical, we need to shift the
595 significand of B down. */
598 /* If the exponents are too far apart, the significands
599 do not overlap, which makes the subtraction a noop. */
600 if (dexp
>= SIGNIFICAND_BITS
)
607 inexact
|= sticky_rshift_significand (&t
, b
, dexp
);
613 if (sub_significands (r
, a
, b
, inexact
))
615 /* We got a borrow out of the subtraction. That means that
616 A and B had the same exponent, and B had the larger
617 significand. We need to swap the sign and negate the
620 neg_significand (r
, r
);
625 if (add_significands (r
, a
, b
))
627 /* We got carry out of the addition. This means we need to
628 shift the significand back down one bit and increase the
630 inexact
|= sticky_rshift_significand (r
, r
, 1);
631 r
->sig
[SIGSZ
-1] |= SIG_MSB
;
640 r
->class = rvc_normal
;
642 SET_REAL_EXP (r
, exp
);
644 /* Re-normalize the result. */
647 /* Special case: if the subtraction results in zero, the result
649 if (r
->class == rvc_zero
)
652 r
->sig
[0] |= inexact
;
657 /* Calculate R = A * B. Return true if the result may be inexact. */
660 do_multiply (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
,
661 const REAL_VALUE_TYPE
*b
)
663 REAL_VALUE_TYPE u
, t
, *rr
;
664 unsigned int i
, j
, k
;
665 int sign
= a
->sign
^ b
->sign
;
666 bool inexact
= false;
668 switch (CLASS2 (a
->class, b
->class))
670 case CLASS2 (rvc_zero
, rvc_zero
):
671 case CLASS2 (rvc_zero
, rvc_normal
):
672 case CLASS2 (rvc_normal
, rvc_zero
):
673 /* +-0 * ANY = 0 with appropriate sign. */
677 case CLASS2 (rvc_zero
, rvc_nan
):
678 case CLASS2 (rvc_normal
, rvc_nan
):
679 case CLASS2 (rvc_inf
, rvc_nan
):
680 case CLASS2 (rvc_nan
, rvc_nan
):
681 /* ANY * NaN = NaN. */
686 case CLASS2 (rvc_nan
, rvc_zero
):
687 case CLASS2 (rvc_nan
, rvc_normal
):
688 case CLASS2 (rvc_nan
, rvc_inf
):
689 /* NaN * ANY = NaN. */
694 case CLASS2 (rvc_zero
, rvc_inf
):
695 case CLASS2 (rvc_inf
, rvc_zero
):
697 get_canonical_qnan (r
, sign
);
700 case CLASS2 (rvc_inf
, rvc_inf
):
701 case CLASS2 (rvc_normal
, rvc_inf
):
702 case CLASS2 (rvc_inf
, rvc_normal
):
703 /* Inf * Inf = Inf, R * Inf = Inf */
707 case CLASS2 (rvc_normal
, rvc_normal
):
714 if (r
== a
|| r
== b
)
720 /* Collect all the partial products. Since we don't have sure access
721 to a widening multiply, we split each long into two half-words.
723 Consider the long-hand form of a four half-word multiplication:
733 We construct partial products of the widened half-word products
734 that are known to not overlap, e.g. DF+DH. Each such partial
735 product is given its proper exponent, which allows us to sum them
736 and obtain the finished product. */
738 for (i
= 0; i
< SIGSZ
* 2; ++i
)
740 unsigned long ai
= a
->sig
[i
/ 2];
742 ai
>>= HOST_BITS_PER_LONG
/ 2;
744 ai
&= ((unsigned long)1 << (HOST_BITS_PER_LONG
/ 2)) - 1;
749 for (j
= 0; j
< 2; ++j
)
751 int exp
= (REAL_EXP (a
) - (2*SIGSZ
-1-i
)*(HOST_BITS_PER_LONG
/2)
752 + (REAL_EXP (b
) - (1-j
)*(HOST_BITS_PER_LONG
/2)));
761 /* Would underflow to zero, which we shouldn't bother adding. */
766 memset (&u
, 0, sizeof (u
));
767 u
.class = rvc_normal
;
768 SET_REAL_EXP (&u
, exp
);
770 for (k
= j
; k
< SIGSZ
* 2; k
+= 2)
772 unsigned long bi
= b
->sig
[k
/ 2];
774 bi
>>= HOST_BITS_PER_LONG
/ 2;
776 bi
&= ((unsigned long)1 << (HOST_BITS_PER_LONG
/ 2)) - 1;
778 u
.sig
[k
/ 2] = ai
* bi
;
782 inexact
|= do_add (rr
, rr
, &u
, 0);
793 /* Calculate R = A / B. Return true if the result may be inexact. */
796 do_divide (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
,
797 const REAL_VALUE_TYPE
*b
)
799 int exp
, sign
= a
->sign
^ b
->sign
;
800 REAL_VALUE_TYPE t
, *rr
;
803 switch (CLASS2 (a
->class, b
->class))
805 case CLASS2 (rvc_zero
, rvc_zero
):
807 case CLASS2 (rvc_inf
, rvc_inf
):
808 /* Inf / Inf = NaN. */
809 get_canonical_qnan (r
, sign
);
812 case CLASS2 (rvc_zero
, rvc_normal
):
813 case CLASS2 (rvc_zero
, rvc_inf
):
815 case CLASS2 (rvc_normal
, rvc_inf
):
820 case CLASS2 (rvc_normal
, rvc_zero
):
822 case CLASS2 (rvc_inf
, rvc_zero
):
827 case CLASS2 (rvc_zero
, rvc_nan
):
828 case CLASS2 (rvc_normal
, rvc_nan
):
829 case CLASS2 (rvc_inf
, rvc_nan
):
830 case CLASS2 (rvc_nan
, rvc_nan
):
831 /* ANY / NaN = NaN. */
836 case CLASS2 (rvc_nan
, rvc_zero
):
837 case CLASS2 (rvc_nan
, rvc_normal
):
838 case CLASS2 (rvc_nan
, rvc_inf
):
839 /* NaN / ANY = NaN. */
844 case CLASS2 (rvc_inf
, rvc_normal
):
849 case CLASS2 (rvc_normal
, rvc_normal
):
856 if (r
== a
|| r
== b
)
861 /* Make sure all fields in the result are initialized. */
863 rr
->class = rvc_normal
;
866 exp
= REAL_EXP (a
) - REAL_EXP (b
) + 1;
877 SET_REAL_EXP (rr
, exp
);
879 inexact
= div_significands (rr
, a
, b
);
881 /* Re-normalize the result. */
883 rr
->sig
[0] |= inexact
;
891 /* Return a tri-state comparison of A vs B. Return NAN_RESULT if
892 one of the two operands is a NaN. */
895 do_compare (const REAL_VALUE_TYPE
*a
, const REAL_VALUE_TYPE
*b
,
900 switch (CLASS2 (a
->class, b
->class))
902 case CLASS2 (rvc_zero
, rvc_zero
):
903 /* Sign of zero doesn't matter for compares. */
906 case CLASS2 (rvc_inf
, rvc_zero
):
907 case CLASS2 (rvc_inf
, rvc_normal
):
908 case CLASS2 (rvc_normal
, rvc_zero
):
909 return (a
->sign
? -1 : 1);
911 case CLASS2 (rvc_inf
, rvc_inf
):
912 return -a
->sign
- -b
->sign
;
914 case CLASS2 (rvc_zero
, rvc_normal
):
915 case CLASS2 (rvc_zero
, rvc_inf
):
916 case CLASS2 (rvc_normal
, rvc_inf
):
917 return (b
->sign
? 1 : -1);
919 case CLASS2 (rvc_zero
, rvc_nan
):
920 case CLASS2 (rvc_normal
, rvc_nan
):
921 case CLASS2 (rvc_inf
, rvc_nan
):
922 case CLASS2 (rvc_nan
, rvc_nan
):
923 case CLASS2 (rvc_nan
, rvc_zero
):
924 case CLASS2 (rvc_nan
, rvc_normal
):
925 case CLASS2 (rvc_nan
, rvc_inf
):
928 case CLASS2 (rvc_normal
, rvc_normal
):
935 if (a
->sign
!= b
->sign
)
936 return -a
->sign
- -b
->sign
;
938 if (REAL_EXP (a
) > REAL_EXP (b
))
940 else if (REAL_EXP (a
) < REAL_EXP (b
))
943 ret
= cmp_significands (a
, b
);
945 return (a
->sign
? -ret
: ret
);
948 /* Return A truncated to an integral value toward zero. */
951 do_fix_trunc (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*a
)
963 if (REAL_EXP (r
) <= 0)
964 get_zero (r
, r
->sign
);
965 else if (REAL_EXP (r
) < SIGNIFICAND_BITS
)
966 clear_significand_below (r
, SIGNIFICAND_BITS
- REAL_EXP (r
));
974 /* Perform the binary or unary operation described by CODE.
975 For a unary operation, leave OP1 NULL. */
978 real_arithmetic (REAL_VALUE_TYPE
*r
, int icode
, const REAL_VALUE_TYPE
*op0
,
979 const REAL_VALUE_TYPE
*op1
)
981 enum tree_code code
= icode
;
986 do_add (r
, op0
, op1
, 0);
990 do_add (r
, op0
, op1
, 1);
994 do_multiply (r
, op0
, op1
);
998 do_divide (r
, op0
, op1
);
1002 if (op1
->class == rvc_nan
)
1004 else if (do_compare (op0
, op1
, -1) < 0)
1011 if (op1
->class == rvc_nan
)
1013 else if (do_compare (op0
, op1
, 1) < 0)
1029 case FIX_TRUNC_EXPR
:
1030 do_fix_trunc (r
, op0
);
1038 /* Legacy. Similar, but return the result directly. */
1041 real_arithmetic2 (int icode
, const REAL_VALUE_TYPE
*op0
,
1042 const REAL_VALUE_TYPE
*op1
)
1045 real_arithmetic (&r
, icode
, op0
, op1
);
1050 real_compare (int icode
, const REAL_VALUE_TYPE
*op0
,
1051 const REAL_VALUE_TYPE
*op1
)
1053 enum tree_code code
= icode
;
1058 return do_compare (op0
, op1
, 1) < 0;
1060 return do_compare (op0
, op1
, 1) <= 0;
1062 return do_compare (op0
, op1
, -1) > 0;
1064 return do_compare (op0
, op1
, -1) >= 0;
1066 return do_compare (op0
, op1
, -1) == 0;
1068 return do_compare (op0
, op1
, -1) != 0;
1069 case UNORDERED_EXPR
:
1070 return op0
->class == rvc_nan
|| op1
->class == rvc_nan
;
1072 return op0
->class != rvc_nan
&& op1
->class != rvc_nan
;
1074 return do_compare (op0
, op1
, -1) < 0;
1076 return do_compare (op0
, op1
, -1) <= 0;
1078 return do_compare (op0
, op1
, 1) > 0;
1080 return do_compare (op0
, op1
, 1) >= 0;
1082 return do_compare (op0
, op1
, 0) == 0;
1089 /* Return floor log2(R). */
1092 real_exponent (const REAL_VALUE_TYPE
*r
)
1100 return (unsigned int)-1 >> 1;
1102 return REAL_EXP (r
);
1108 /* R = OP0 * 2**EXP. */
1111 real_ldexp (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*op0
, int exp
)
1122 exp
+= REAL_EXP (op0
);
1124 get_inf (r
, r
->sign
);
1125 else if (exp
< -MAX_EXP
)
1126 get_zero (r
, r
->sign
);
1128 SET_REAL_EXP (r
, exp
);
1136 /* Determine whether a floating-point value X is infinite. */
1139 real_isinf (const REAL_VALUE_TYPE
*r
)
1141 return (r
->class == rvc_inf
);
1144 /* Determine whether a floating-point value X is a NaN. */
1147 real_isnan (const REAL_VALUE_TYPE
*r
)
1149 return (r
->class == rvc_nan
);
1152 /* Determine whether a floating-point value X is negative. */
1155 real_isneg (const REAL_VALUE_TYPE
*r
)
1160 /* Determine whether a floating-point value X is minus zero. */
1163 real_isnegzero (const REAL_VALUE_TYPE
*r
)
1165 return r
->sign
&& r
->class == rvc_zero
;
1168 /* Compare two floating-point objects for bitwise identity. */
1171 real_identical (const REAL_VALUE_TYPE
*a
, const REAL_VALUE_TYPE
*b
)
1175 if (a
->class != b
->class)
1177 if (a
->sign
!= b
->sign
)
1187 if (REAL_EXP (a
) != REAL_EXP (b
))
1192 if (a
->signalling
!= b
->signalling
)
1194 /* The significand is ignored for canonical NaNs. */
1195 if (a
->canonical
|| b
->canonical
)
1196 return a
->canonical
== b
->canonical
;
1203 for (i
= 0; i
< SIGSZ
; ++i
)
1204 if (a
->sig
[i
] != b
->sig
[i
])
1210 /* Try to change R into its exact multiplicative inverse in machine
1211 mode MODE. Return true if successful. */
1214 exact_real_inverse (enum machine_mode mode
, REAL_VALUE_TYPE
*r
)
1216 const REAL_VALUE_TYPE
*one
= real_digit (1);
1220 if (r
->class != rvc_normal
)
1223 /* Check for a power of two: all significand bits zero except the MSB. */
1224 for (i
= 0; i
< SIGSZ
-1; ++i
)
1227 if (r
->sig
[SIGSZ
-1] != SIG_MSB
)
1230 /* Find the inverse and truncate to the required mode. */
1231 do_divide (&u
, one
, r
);
1232 real_convert (&u
, mode
, &u
);
1234 /* The rounding may have overflowed. */
1235 if (u
.class != rvc_normal
)
1237 for (i
= 0; i
< SIGSZ
-1; ++i
)
1240 if (u
.sig
[SIGSZ
-1] != SIG_MSB
)
1247 /* Render R as an integer. */
1250 real_to_integer (const REAL_VALUE_TYPE
*r
)
1252 unsigned HOST_WIDE_INT i
;
1263 i
= (unsigned HOST_WIDE_INT
) 1 << (HOST_BITS_PER_WIDE_INT
- 1);
1269 if (REAL_EXP (r
) <= 0)
1271 /* Only force overflow for unsigned overflow. Signed overflow is
1272 undefined, so it doesn't matter what we return, and some callers
1273 expect to be able to use this routine for both signed and
1274 unsigned conversions. */
1275 if (REAL_EXP (r
) > HOST_BITS_PER_WIDE_INT
)
1278 if (HOST_BITS_PER_WIDE_INT
== HOST_BITS_PER_LONG
)
1279 i
= r
->sig
[SIGSZ
-1];
1280 else if (HOST_BITS_PER_WIDE_INT
== 2*HOST_BITS_PER_LONG
)
1282 i
= r
->sig
[SIGSZ
-1];
1283 i
= i
<< (HOST_BITS_PER_LONG
- 1) << 1;
1284 i
|= r
->sig
[SIGSZ
-2];
1289 i
>>= HOST_BITS_PER_WIDE_INT
- REAL_EXP (r
);
1300 /* Likewise, but to an integer pair, HI+LOW. */
1303 real_to_integer2 (HOST_WIDE_INT
*plow
, HOST_WIDE_INT
*phigh
,
1304 const REAL_VALUE_TYPE
*r
)
1307 HOST_WIDE_INT low
, high
;
1320 high
= (unsigned HOST_WIDE_INT
) 1 << (HOST_BITS_PER_WIDE_INT
- 1);
1334 /* Only force overflow for unsigned overflow. Signed overflow is
1335 undefined, so it doesn't matter what we return, and some callers
1336 expect to be able to use this routine for both signed and
1337 unsigned conversions. */
1338 if (exp
> 2*HOST_BITS_PER_WIDE_INT
)
1341 rshift_significand (&t
, r
, 2*HOST_BITS_PER_WIDE_INT
- exp
);
1342 if (HOST_BITS_PER_WIDE_INT
== HOST_BITS_PER_LONG
)
1344 high
= t
.sig
[SIGSZ
-1];
1345 low
= t
.sig
[SIGSZ
-2];
1347 else if (HOST_BITS_PER_WIDE_INT
== 2*HOST_BITS_PER_LONG
)
1349 high
= t
.sig
[SIGSZ
-1];
1350 high
= high
<< (HOST_BITS_PER_LONG
- 1) << 1;
1351 high
|= t
.sig
[SIGSZ
-2];
1353 low
= t
.sig
[SIGSZ
-3];
1354 low
= low
<< (HOST_BITS_PER_LONG
- 1) << 1;
1355 low
|= t
.sig
[SIGSZ
-4];
1365 low
= -low
, high
= ~high
;
1377 /* A subroutine of real_to_decimal. Compute the quotient and remainder
1378 of NUM / DEN. Return the quotient and place the remainder in NUM.
1379 It is expected that NUM / DEN are close enough that the quotient is
1382 static unsigned long
1383 rtd_divmod (REAL_VALUE_TYPE
*num
, REAL_VALUE_TYPE
*den
)
1385 unsigned long q
, msb
;
1386 int expn
= REAL_EXP (num
), expd
= REAL_EXP (den
);
1395 msb
= num
->sig
[SIGSZ
-1] & SIG_MSB
;
1397 lshift_significand_1 (num
, num
);
1399 if (msb
|| cmp_significands (num
, den
) >= 0)
1401 sub_significands (num
, num
, den
, 0);
1405 while (--expn
>= expd
);
1407 SET_REAL_EXP (num
, expd
);
1413 /* Render R as a decimal floating point constant. Emit DIGITS significant
1414 digits in the result, bounded by BUF_SIZE. If DIGITS is 0, choose the
1415 maximum for the representation. If CROP_TRAILING_ZEROS, strip trailing
1418 #define M_LOG10_2 0.30102999566398119521
1421 real_to_decimal (char *str
, const REAL_VALUE_TYPE
*r_orig
, size_t buf_size
,
1422 size_t digits
, int crop_trailing_zeros
)
1424 const REAL_VALUE_TYPE
*one
, *ten
;
1425 REAL_VALUE_TYPE r
, pten
, u
, v
;
1426 int dec_exp
, cmp_one
, digit
;
1428 char *p
, *first
, *last
;
1435 strcpy (str
, (r
.sign
? "-0.0" : "0.0"));
1440 strcpy (str
, (r
.sign
? "-Inf" : "+Inf"));
1443 /* ??? Print the significand as well, if not canonical? */
1444 strcpy (str
, (r
.sign
? "-NaN" : "+NaN"));
1450 /* Bound the number of digits printed by the size of the representation. */
1451 max_digits
= SIGNIFICAND_BITS
* M_LOG10_2
;
1452 if (digits
== 0 || digits
> max_digits
)
1453 digits
= max_digits
;
1455 /* Estimate the decimal exponent, and compute the length of the string it
1456 will print as. Be conservative and add one to account for possible
1457 overflow or rounding error. */
1458 dec_exp
= REAL_EXP (&r
) * M_LOG10_2
;
1459 for (max_digits
= 1; dec_exp
; max_digits
++)
1462 /* Bound the number of digits printed by the size of the output buffer. */
1463 max_digits
= buf_size
- 1 - 1 - 2 - max_digits
- 1;
1464 if (max_digits
> buf_size
)
1466 if (digits
> max_digits
)
1467 digits
= max_digits
;
1469 one
= real_digit (1);
1470 ten
= ten_to_ptwo (0);
1478 cmp_one
= do_compare (&r
, one
, 0);
1483 /* Number is greater than one. Convert significand to an integer
1484 and strip trailing decimal zeros. */
1487 SET_REAL_EXP (&u
, SIGNIFICAND_BITS
- 1);
1489 /* Largest M, such that 10**2**M fits within SIGNIFICAND_BITS. */
1490 m
= floor_log2 (max_digits
);
1492 /* Iterate over the bits of the possible powers of 10 that might
1493 be present in U and eliminate them. That is, if we find that
1494 10**2**M divides U evenly, keep the division and increase
1500 do_divide (&t
, &u
, ten_to_ptwo (m
));
1501 do_fix_trunc (&v
, &t
);
1502 if (cmp_significands (&v
, &t
) == 0)
1510 /* Revert the scaling to integer that we performed earlier. */
1511 SET_REAL_EXP (&u
, REAL_EXP (&u
) + REAL_EXP (&r
)
1512 - (SIGNIFICAND_BITS
- 1));
1515 /* Find power of 10. Do this by dividing out 10**2**M when
1516 this is larger than the current remainder. Fill PTEN with
1517 the power of 10 that we compute. */
1518 if (REAL_EXP (&r
) > 0)
1520 m
= floor_log2 ((int)(REAL_EXP (&r
) * M_LOG10_2
)) + 1;
1523 const REAL_VALUE_TYPE
*ptentwo
= ten_to_ptwo (m
);
1524 if (do_compare (&u
, ptentwo
, 0) >= 0)
1526 do_divide (&u
, &u
, ptentwo
);
1527 do_multiply (&pten
, &pten
, ptentwo
);
1534 /* We managed to divide off enough tens in the above reduction
1535 loop that we've now got a negative exponent. Fall into the
1536 less-than-one code to compute the proper value for PTEN. */
1543 /* Number is less than one. Pad significand with leading
1549 /* Stop if we'd shift bits off the bottom. */
1553 do_multiply (&u
, &v
, ten
);
1555 /* Stop if we're now >= 1. */
1556 if (REAL_EXP (&u
) > 0)
1564 /* Find power of 10. Do this by multiplying in P=10**2**M when
1565 the current remainder is smaller than 1/P. Fill PTEN with the
1566 power of 10 that we compute. */
1567 m
= floor_log2 ((int)(-REAL_EXP (&r
) * M_LOG10_2
)) + 1;
1570 const REAL_VALUE_TYPE
*ptentwo
= ten_to_ptwo (m
);
1571 const REAL_VALUE_TYPE
*ptenmtwo
= ten_to_mptwo (m
);
1573 if (do_compare (&v
, ptenmtwo
, 0) <= 0)
1575 do_multiply (&v
, &v
, ptentwo
);
1576 do_multiply (&pten
, &pten
, ptentwo
);
1582 /* Invert the positive power of 10 that we've collected so far. */
1583 do_divide (&pten
, one
, &pten
);
1591 /* At this point, PTEN should contain the nearest power of 10 smaller
1592 than R, such that this division produces the first digit.
1594 Using a divide-step primitive that returns the complete integral
1595 remainder avoids the rounding error that would be produced if
1596 we were to use do_divide here and then simply multiply by 10 for
1597 each subsequent digit. */
1599 digit
= rtd_divmod (&r
, &pten
);
1601 /* Be prepared for error in that division via underflow ... */
1602 if (digit
== 0 && cmp_significand_0 (&r
))
1604 /* Multiply by 10 and try again. */
1605 do_multiply (&r
, &r
, ten
);
1606 digit
= rtd_divmod (&r
, &pten
);
1612 /* ... or overflow. */
1620 else if (digit
> 10)
1625 /* Generate subsequent digits. */
1626 while (--digits
> 0)
1628 do_multiply (&r
, &r
, ten
);
1629 digit
= rtd_divmod (&r
, &pten
);
1634 /* Generate one more digit with which to do rounding. */
1635 do_multiply (&r
, &r
, ten
);
1636 digit
= rtd_divmod (&r
, &pten
);
1638 /* Round the result. */
1641 /* Round to nearest. If R is nonzero there are additional
1642 nonzero digits to be extracted. */
1643 if (cmp_significand_0 (&r
))
1645 /* Round to even. */
1646 else if ((p
[-1] - '0') & 1)
1663 /* Carry out of the first digit. This means we had all 9's and
1664 now have all 0's. "Prepend" a 1 by overwriting the first 0. */
1672 /* Insert the decimal point. */
1673 first
[0] = first
[1];
1676 /* If requested, drop trailing zeros. Never crop past "1.0". */
1677 if (crop_trailing_zeros
)
1678 while (last
> first
+ 3 && last
[-1] == '0')
1681 /* Append the exponent. */
1682 sprintf (last
, "e%+d", dec_exp
);
1685 /* Render R as a hexadecimal floating point constant. Emit DIGITS
1686 significant digits in the result, bounded by BUF_SIZE. If DIGITS is 0,
1687 choose the maximum for the representation. If CROP_TRAILING_ZEROS,
1688 strip trailing zeros. */
1691 real_to_hexadecimal (char *str
, const REAL_VALUE_TYPE
*r
, size_t buf_size
,
1692 size_t digits
, int crop_trailing_zeros
)
1694 int i
, j
, exp
= REAL_EXP (r
);
1707 strcpy (str
, (r
->sign
? "-Inf" : "+Inf"));
1710 /* ??? Print the significand as well, if not canonical? */
1711 strcpy (str
, (r
->sign
? "-NaN" : "+NaN"));
1718 digits
= SIGNIFICAND_BITS
/ 4;
1720 /* Bound the number of digits printed by the size of the output buffer. */
1722 sprintf (exp_buf
, "p%+d", exp
);
1723 max_digits
= buf_size
- strlen (exp_buf
) - r
->sign
- 4 - 1;
1724 if (max_digits
> buf_size
)
1726 if (digits
> max_digits
)
1727 digits
= max_digits
;
1738 for (i
= SIGSZ
- 1; i
>= 0; --i
)
1739 for (j
= HOST_BITS_PER_LONG
- 4; j
>= 0; j
-= 4)
1741 *p
++ = "0123456789abcdef"[(r
->sig
[i
] >> j
) & 15];
1747 if (crop_trailing_zeros
)
1748 while (p
> first
+ 1 && p
[-1] == '0')
1751 sprintf (p
, "p%+d", exp
);
1754 /* Initialize R from a decimal or hexadecimal string. The string is
1755 assumed to have been syntax checked already. */
1758 real_from_string (REAL_VALUE_TYPE
*r
, const char *str
)
1770 else if (*str
== '+')
1773 if (str
[0] == '0' && (str
[1] == 'x' || str
[1] == 'X'))
1775 /* Hexadecimal floating point. */
1776 int pos
= SIGNIFICAND_BITS
- 4, d
;
1784 d
= hex_value (*str
);
1789 r
->sig
[pos
/ HOST_BITS_PER_LONG
]
1790 |= (unsigned long) d
<< (pos
% HOST_BITS_PER_LONG
);
1799 if (pos
== SIGNIFICAND_BITS
- 4)
1806 d
= hex_value (*str
);
1811 r
->sig
[pos
/ HOST_BITS_PER_LONG
]
1812 |= (unsigned long) d
<< (pos
% HOST_BITS_PER_LONG
);
1818 if (*str
== 'p' || *str
== 'P')
1820 bool exp_neg
= false;
1828 else if (*str
== '+')
1832 while (ISDIGIT (*str
))
1838 /* Overflowed the exponent. */
1852 r
->class = rvc_normal
;
1853 SET_REAL_EXP (r
, exp
);
1859 /* Decimal floating point. */
1860 const REAL_VALUE_TYPE
*ten
= ten_to_ptwo (0);
1865 while (ISDIGIT (*str
))
1868 do_multiply (r
, r
, ten
);
1870 do_add (r
, r
, real_digit (d
), 0);
1875 if (r
->class == rvc_zero
)
1880 while (ISDIGIT (*str
))
1883 do_multiply (r
, r
, ten
);
1885 do_add (r
, r
, real_digit (d
), 0);
1890 if (*str
== 'e' || *str
== 'E')
1892 bool exp_neg
= false;
1900 else if (*str
== '+')
1904 while (ISDIGIT (*str
))
1910 /* Overflowed the exponent. */
1924 times_pten (r
, exp
);
1939 /* Legacy. Similar, but return the result directly. */
1942 real_from_string2 (const char *s
, enum machine_mode mode
)
1946 real_from_string (&r
, s
);
1947 if (mode
!= VOIDmode
)
1948 real_convert (&r
, mode
, &r
);
1953 /* Initialize R from the integer pair HIGH+LOW. */
1956 real_from_integer (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
1957 unsigned HOST_WIDE_INT low
, HOST_WIDE_INT high
,
1960 if (low
== 0 && high
== 0)
1964 r
->class = rvc_normal
;
1965 r
->sign
= high
< 0 && !unsigned_p
;
1966 SET_REAL_EXP (r
, 2 * HOST_BITS_PER_WIDE_INT
);
1977 if (HOST_BITS_PER_LONG
== HOST_BITS_PER_WIDE_INT
)
1979 r
->sig
[SIGSZ
-1] = high
;
1980 r
->sig
[SIGSZ
-2] = low
;
1981 memset (r
->sig
, 0, sizeof(long)*(SIGSZ
-2));
1983 else if (HOST_BITS_PER_LONG
*2 == HOST_BITS_PER_WIDE_INT
)
1985 r
->sig
[SIGSZ
-1] = high
>> (HOST_BITS_PER_LONG
- 1) >> 1;
1986 r
->sig
[SIGSZ
-2] = high
;
1987 r
->sig
[SIGSZ
-3] = low
>> (HOST_BITS_PER_LONG
- 1) >> 1;
1988 r
->sig
[SIGSZ
-4] = low
;
1990 memset (r
->sig
, 0, sizeof(long)*(SIGSZ
-4));
1998 if (mode
!= VOIDmode
)
1999 real_convert (r
, mode
, r
);
2002 /* Returns 10**2**N. */
2004 static const REAL_VALUE_TYPE
*
2007 static REAL_VALUE_TYPE tens
[EXP_BITS
];
2009 if (n
< 0 || n
>= EXP_BITS
)
2012 if (tens
[n
].class == rvc_zero
)
2014 if (n
< (HOST_BITS_PER_WIDE_INT
== 64 ? 5 : 4))
2016 HOST_WIDE_INT t
= 10;
2019 for (i
= 0; i
< n
; ++i
)
2022 real_from_integer (&tens
[n
], VOIDmode
, t
, 0, 1);
2026 const REAL_VALUE_TYPE
*t
= ten_to_ptwo (n
- 1);
2027 do_multiply (&tens
[n
], t
, t
);
2034 /* Returns 10**(-2**N). */
2036 static const REAL_VALUE_TYPE
*
2037 ten_to_mptwo (int n
)
2039 static REAL_VALUE_TYPE tens
[EXP_BITS
];
2041 if (n
< 0 || n
>= EXP_BITS
)
2044 if (tens
[n
].class == rvc_zero
)
2045 do_divide (&tens
[n
], real_digit (1), ten_to_ptwo (n
));
2052 static const REAL_VALUE_TYPE
*
2055 static REAL_VALUE_TYPE num
[10];
2060 if (n
> 0 && num
[n
].class == rvc_zero
)
2061 real_from_integer (&num
[n
], VOIDmode
, n
, 0, 1);
2066 /* Multiply R by 10**EXP. */
2069 times_pten (REAL_VALUE_TYPE
*r
, int exp
)
2071 REAL_VALUE_TYPE pten
, *rr
;
2072 bool negative
= (exp
< 0);
2078 pten
= *real_digit (1);
2084 for (i
= 0; exp
> 0; ++i
, exp
>>= 1)
2086 do_multiply (rr
, rr
, ten_to_ptwo (i
));
2089 do_divide (r
, r
, &pten
);
2092 /* Fills R with +Inf. */
2095 real_inf (REAL_VALUE_TYPE
*r
)
2100 /* Fills R with a NaN whose significand is described by STR. If QUIET,
2101 we force a QNaN, else we force an SNaN. The string, if not empty,
2102 is parsed as a number and placed in the significand. Return true
2103 if the string was successfully parsed. */
2106 real_nan (REAL_VALUE_TYPE
*r
, const char *str
, int quiet
,
2107 enum machine_mode mode
)
2109 const struct real_format
*fmt
;
2111 fmt
= REAL_MODE_FORMAT (mode
);
2118 get_canonical_qnan (r
, 0);
2120 get_canonical_snan (r
, 0);
2127 memset (r
, 0, sizeof (*r
));
2130 /* Parse akin to strtol into the significand of R. */
2132 while (ISSPACE (*str
))
2136 else if (*str
== '+')
2146 while ((d
= hex_value (*str
)) < base
)
2153 lshift_significand (r
, r
, 3);
2156 lshift_significand (r
, r
, 4);
2159 lshift_significand_1 (&u
, r
);
2160 lshift_significand (r
, r
, 3);
2161 add_significands (r
, r
, &u
);
2169 add_significands (r
, r
, &u
);
2174 /* Must have consumed the entire string for success. */
2178 /* Shift the significand into place such that the bits
2179 are in the most significant bits for the format. */
2180 lshift_significand (r
, r
, SIGNIFICAND_BITS
- fmt
->pnan
);
2182 /* Our MSB is always unset for NaNs. */
2183 r
->sig
[SIGSZ
-1] &= ~SIG_MSB
;
2185 /* Force quiet or signalling NaN. */
2186 r
->signalling
= !quiet
;
2192 /* Fills R with the largest finite value representable in mode MODE.
2193 If SIGN is nonzero, R is set to the most negative finite value. */
2196 real_maxval (REAL_VALUE_TYPE
*r
, int sign
, enum machine_mode mode
)
2198 const struct real_format
*fmt
;
2201 fmt
= REAL_MODE_FORMAT (mode
);
2205 r
->class = rvc_normal
;
2209 SET_REAL_EXP (r
, fmt
->emax
* fmt
->log2_b
);
2211 np2
= SIGNIFICAND_BITS
- fmt
->p
* fmt
->log2_b
;
2212 memset (r
->sig
, -1, SIGSZ
* sizeof (unsigned long));
2213 clear_significand_below (r
, np2
);
2216 /* Fills R with 2**N. */
2219 real_2expN (REAL_VALUE_TYPE
*r
, int n
)
2221 memset (r
, 0, sizeof (*r
));
2226 else if (n
< -MAX_EXP
)
2230 r
->class = rvc_normal
;
2231 SET_REAL_EXP (r
, n
);
2232 r
->sig
[SIGSZ
-1] = SIG_MSB
;
2238 round_for_format (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
)
2241 unsigned long sticky
;
2245 p2
= fmt
->p
* fmt
->log2_b
;
2246 emin2m1
= (fmt
->emin
- 1) * fmt
->log2_b
;
2247 emax2
= fmt
->emax
* fmt
->log2_b
;
2249 np2
= SIGNIFICAND_BITS
- p2
;
2253 get_zero (r
, r
->sign
);
2255 if (!fmt
->has_signed_zero
)
2260 get_inf (r
, r
->sign
);
2265 clear_significand_below (r
, np2
);
2275 /* If we're not base2, normalize the exponent to a multiple of
2277 if (fmt
->log2_b
!= 1)
2279 int shift
= REAL_EXP (r
) & (fmt
->log2_b
- 1);
2282 shift
= fmt
->log2_b
- shift
;
2283 r
->sig
[0] |= sticky_rshift_significand (r
, r
, shift
);
2284 SET_REAL_EXP (r
, REAL_EXP (r
) + shift
);
2288 /* Check the range of the exponent. If we're out of range,
2289 either underflow or overflow. */
2290 if (REAL_EXP (r
) > emax2
)
2292 else if (REAL_EXP (r
) <= emin2m1
)
2296 if (!fmt
->has_denorm
)
2298 /* Don't underflow completely until we've had a chance to round. */
2299 if (REAL_EXP (r
) < emin2m1
)
2304 diff
= emin2m1
- REAL_EXP (r
) + 1;
2308 /* De-normalize the significand. */
2309 r
->sig
[0] |= sticky_rshift_significand (r
, r
, diff
);
2310 SET_REAL_EXP (r
, REAL_EXP (r
) + diff
);
2314 /* There are P2 true significand bits, followed by one guard bit,
2315 followed by one sticky bit, followed by stuff. Fold nonzero
2316 stuff into the sticky bit. */
2319 for (i
= 0, w
= (np2
- 1) / HOST_BITS_PER_LONG
; i
< w
; ++i
)
2320 sticky
|= r
->sig
[i
];
2322 r
->sig
[w
] & (((unsigned long)1 << ((np2
- 1) % HOST_BITS_PER_LONG
)) - 1);
2324 guard
= test_significand_bit (r
, np2
- 1);
2325 lsb
= test_significand_bit (r
, np2
);
2327 /* Round to even. */
2328 if (guard
&& (sticky
|| lsb
))
2332 set_significand_bit (&u
, np2
);
2334 if (add_significands (r
, r
, &u
))
2336 /* Overflow. Means the significand had been all ones, and
2337 is now all zeros. Need to increase the exponent, and
2338 possibly re-normalize it. */
2339 SET_REAL_EXP (r
, REAL_EXP (r
) + 1);
2340 if (REAL_EXP (r
) > emax2
)
2342 r
->sig
[SIGSZ
-1] = SIG_MSB
;
2344 if (fmt
->log2_b
!= 1)
2346 int shift
= REAL_EXP (r
) & (fmt
->log2_b
- 1);
2349 shift
= fmt
->log2_b
- shift
;
2350 rshift_significand (r
, r
, shift
);
2351 SET_REAL_EXP (r
, REAL_EXP (r
) + shift
);
2352 if (REAL_EXP (r
) > emax2
)
2359 /* Catch underflow that we deferred until after rounding. */
2360 if (REAL_EXP (r
) <= emin2m1
)
2363 /* Clear out trailing garbage. */
2364 clear_significand_below (r
, np2
);
2367 /* Extend or truncate to a new mode. */
2370 real_convert (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
2371 const REAL_VALUE_TYPE
*a
)
2373 const struct real_format
*fmt
;
2375 fmt
= REAL_MODE_FORMAT (mode
);
2380 round_for_format (fmt
, r
);
2382 /* round_for_format de-normalizes denormals. Undo just that part. */
2383 if (r
->class == rvc_normal
)
2387 /* Legacy. Likewise, except return the struct directly. */
2390 real_value_truncate (enum machine_mode mode
, REAL_VALUE_TYPE a
)
2393 real_convert (&r
, mode
, &a
);
2397 /* Return true if truncating to MODE is exact. */
2400 exact_real_truncate (enum machine_mode mode
, const REAL_VALUE_TYPE
*a
)
2403 real_convert (&t
, mode
, a
);
2404 return real_identical (&t
, a
);
2407 /* Write R to the given target format. Place the words of the result
2408 in target word order in BUF. There are always 32 bits in each
2409 long, no matter the size of the host long.
2411 Legacy: return word 0 for implementing REAL_VALUE_TO_TARGET_SINGLE. */
2414 real_to_target_fmt (long *buf
, const REAL_VALUE_TYPE
*r_orig
,
2415 const struct real_format
*fmt
)
2421 round_for_format (fmt
, &r
);
2425 (*fmt
->encode
) (fmt
, buf
, &r
);
2430 /* Similar, but look up the format from MODE. */
2433 real_to_target (long *buf
, const REAL_VALUE_TYPE
*r
, enum machine_mode mode
)
2435 const struct real_format
*fmt
;
2437 fmt
= REAL_MODE_FORMAT (mode
);
2441 return real_to_target_fmt (buf
, r
, fmt
);
2444 /* Read R from the given target format. Read the words of the result
2445 in target word order in BUF. There are always 32 bits in each
2446 long, no matter the size of the host long. */
2449 real_from_target_fmt (REAL_VALUE_TYPE
*r
, const long *buf
,
2450 const struct real_format
*fmt
)
2452 (*fmt
->decode
) (fmt
, r
, buf
);
2455 /* Similar, but look up the format from MODE. */
2458 real_from_target (REAL_VALUE_TYPE
*r
, const long *buf
, enum machine_mode mode
)
2460 const struct real_format
*fmt
;
2462 fmt
= REAL_MODE_FORMAT (mode
);
2466 (*fmt
->decode
) (fmt
, r
, buf
);
2469 /* Return the number of bits in the significand for MODE. */
2470 /* ??? Legacy. Should get access to real_format directly. */
2473 significand_size (enum machine_mode mode
)
2475 const struct real_format
*fmt
;
2477 fmt
= REAL_MODE_FORMAT (mode
);
2481 return fmt
->p
* fmt
->log2_b
;
2484 /* Return a hash value for the given real value. */
2485 /* ??? The "unsigned int" return value is intended to be hashval_t,
2486 but I didn't want to pull hashtab.h into real.h. */
2489 real_hash (const REAL_VALUE_TYPE
*r
)
2494 h
= r
->class | (r
->sign
<< 2);
2502 h
|= REAL_EXP (r
) << 3;
2507 h
^= (unsigned int)-1;
2516 if (sizeof(unsigned long) > sizeof(unsigned int))
2517 for (i
= 0; i
< SIGSZ
; ++i
)
2519 unsigned long s
= r
->sig
[i
];
2520 h
^= s
^ (s
>> (HOST_BITS_PER_LONG
/ 2));
2523 for (i
= 0; i
< SIGSZ
; ++i
)
2529 /* IEEE single-precision format. */
2531 static void encode_ieee_single (const struct real_format
*fmt
,
2532 long *, const REAL_VALUE_TYPE
*);
2533 static void decode_ieee_single (const struct real_format
*,
2534 REAL_VALUE_TYPE
*, const long *);
2537 encode_ieee_single (const struct real_format
*fmt
, long *buf
,
2538 const REAL_VALUE_TYPE
*r
)
2540 unsigned long image
, sig
, exp
;
2541 unsigned long sign
= r
->sign
;
2542 bool denormal
= (r
->sig
[SIGSZ
-1] & SIG_MSB
) == 0;
2545 sig
= (r
->sig
[SIGSZ
-1] >> (HOST_BITS_PER_LONG
- 24)) & 0x7fffff;
2556 image
|= 0x7fffffff;
2564 if (r
->signalling
== fmt
->qnan_msb_set
)
2568 /* We overload qnan_msb_set here: it's only clear for
2569 mips_ieee_single, which wants all mantissa bits but the
2570 quiet/signalling one set in canonical NaNs (at least
2572 if (r
->canonical
&& !fmt
->qnan_msb_set
)
2573 sig
|= (1 << 22) - 1;
2581 image
|= 0x7fffffff;
2585 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2586 whereas the intermediate representation is 0.F x 2**exp.
2587 Which means we're off by one. */
2591 exp
= REAL_EXP (r
) + 127 - 1;
2604 decode_ieee_single (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
,
2607 unsigned long image
= buf
[0] & 0xffffffff;
2608 bool sign
= (image
>> 31) & 1;
2609 int exp
= (image
>> 23) & 0xff;
2611 memset (r
, 0, sizeof (*r
));
2612 image
<<= HOST_BITS_PER_LONG
- 24;
2617 if (image
&& fmt
->has_denorm
)
2619 r
->class = rvc_normal
;
2621 SET_REAL_EXP (r
, -126);
2622 r
->sig
[SIGSZ
-1] = image
<< 1;
2625 else if (fmt
->has_signed_zero
)
2628 else if (exp
== 255 && (fmt
->has_nans
|| fmt
->has_inf
))
2634 r
->signalling
= (((image
>> (HOST_BITS_PER_LONG
- 2)) & 1)
2635 ^ fmt
->qnan_msb_set
);
2636 r
->sig
[SIGSZ
-1] = image
;
2646 r
->class = rvc_normal
;
2648 SET_REAL_EXP (r
, exp
- 127 + 1);
2649 r
->sig
[SIGSZ
-1] = image
| SIG_MSB
;
2653 const struct real_format ieee_single_format
=
2671 const struct real_format mips_single_format
=
2690 /* IEEE double-precision format. */
2692 static void encode_ieee_double (const struct real_format
*fmt
,
2693 long *, const REAL_VALUE_TYPE
*);
2694 static void decode_ieee_double (const struct real_format
*,
2695 REAL_VALUE_TYPE
*, const long *);
2698 encode_ieee_double (const struct real_format
*fmt
, long *buf
,
2699 const REAL_VALUE_TYPE
*r
)
2701 unsigned long image_lo
, image_hi
, sig_lo
, sig_hi
, exp
;
2702 bool denormal
= (r
->sig
[SIGSZ
-1] & SIG_MSB
) == 0;
2704 image_hi
= r
->sign
<< 31;
2707 if (HOST_BITS_PER_LONG
== 64)
2709 sig_hi
= r
->sig
[SIGSZ
-1];
2710 sig_lo
= (sig_hi
>> (64 - 53)) & 0xffffffff;
2711 sig_hi
= (sig_hi
>> (64 - 53 + 1) >> 31) & 0xfffff;
2715 sig_hi
= r
->sig
[SIGSZ
-1];
2716 sig_lo
= r
->sig
[SIGSZ
-2];
2717 sig_lo
= (sig_hi
<< 21) | (sig_lo
>> 11);
2718 sig_hi
= (sig_hi
>> 11) & 0xfffff;
2728 image_hi
|= 2047 << 20;
2731 image_hi
|= 0x7fffffff;
2732 image_lo
= 0xffffffff;
2740 sig_hi
= sig_lo
= 0;
2741 if (r
->signalling
== fmt
->qnan_msb_set
)
2742 sig_hi
&= ~(1 << 19);
2745 /* We overload qnan_msb_set here: it's only clear for
2746 mips_ieee_single, which wants all mantissa bits but the
2747 quiet/signalling one set in canonical NaNs (at least
2749 if (r
->canonical
&& !fmt
->qnan_msb_set
)
2751 sig_hi
|= (1 << 19) - 1;
2752 sig_lo
= 0xffffffff;
2754 else if (sig_hi
== 0 && sig_lo
== 0)
2757 image_hi
|= 2047 << 20;
2763 image_hi
|= 0x7fffffff;
2764 image_lo
= 0xffffffff;
2769 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2770 whereas the intermediate representation is 0.F x 2**exp.
2771 Which means we're off by one. */
2775 exp
= REAL_EXP (r
) + 1023 - 1;
2776 image_hi
|= exp
<< 20;
2785 if (FLOAT_WORDS_BIG_ENDIAN
)
2786 buf
[0] = image_hi
, buf
[1] = image_lo
;
2788 buf
[0] = image_lo
, buf
[1] = image_hi
;
2792 decode_ieee_double (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
,
2795 unsigned long image_hi
, image_lo
;
2799 if (FLOAT_WORDS_BIG_ENDIAN
)
2800 image_hi
= buf
[0], image_lo
= buf
[1];
2802 image_lo
= buf
[0], image_hi
= buf
[1];
2803 image_lo
&= 0xffffffff;
2804 image_hi
&= 0xffffffff;
2806 sign
= (image_hi
>> 31) & 1;
2807 exp
= (image_hi
>> 20) & 0x7ff;
2809 memset (r
, 0, sizeof (*r
));
2811 image_hi
<<= 32 - 21;
2812 image_hi
|= image_lo
>> 21;
2813 image_hi
&= 0x7fffffff;
2814 image_lo
<<= 32 - 21;
2818 if ((image_hi
|| image_lo
) && fmt
->has_denorm
)
2820 r
->class = rvc_normal
;
2822 SET_REAL_EXP (r
, -1022);
2823 if (HOST_BITS_PER_LONG
== 32)
2825 image_hi
= (image_hi
<< 1) | (image_lo
>> 31);
2827 r
->sig
[SIGSZ
-1] = image_hi
;
2828 r
->sig
[SIGSZ
-2] = image_lo
;
2832 image_hi
= (image_hi
<< 31 << 2) | (image_lo
<< 1);
2833 r
->sig
[SIGSZ
-1] = image_hi
;
2837 else if (fmt
->has_signed_zero
)
2840 else if (exp
== 2047 && (fmt
->has_nans
|| fmt
->has_inf
))
2842 if (image_hi
|| image_lo
)
2846 r
->signalling
= ((image_hi
>> 30) & 1) ^ fmt
->qnan_msb_set
;
2847 if (HOST_BITS_PER_LONG
== 32)
2849 r
->sig
[SIGSZ
-1] = image_hi
;
2850 r
->sig
[SIGSZ
-2] = image_lo
;
2853 r
->sig
[SIGSZ
-1] = (image_hi
<< 31 << 1) | image_lo
;
2863 r
->class = rvc_normal
;
2865 SET_REAL_EXP (r
, exp
- 1023 + 1);
2866 if (HOST_BITS_PER_LONG
== 32)
2868 r
->sig
[SIGSZ
-1] = image_hi
| SIG_MSB
;
2869 r
->sig
[SIGSZ
-2] = image_lo
;
2872 r
->sig
[SIGSZ
-1] = (image_hi
<< 31 << 1) | image_lo
| SIG_MSB
;
2876 const struct real_format ieee_double_format
=
2894 const struct real_format mips_double_format
=
2913 /* IEEE extended double precision format. This comes in three
2914 flavors: Intel's as a 12 byte image, Intel's as a 16 byte image,
2917 static void encode_ieee_extended (const struct real_format
*fmt
,
2918 long *, const REAL_VALUE_TYPE
*);
2919 static void decode_ieee_extended (const struct real_format
*,
2920 REAL_VALUE_TYPE
*, const long *);
2922 static void encode_ieee_extended_128 (const struct real_format
*fmt
,
2923 long *, const REAL_VALUE_TYPE
*);
2924 static void decode_ieee_extended_128 (const struct real_format
*,
2925 REAL_VALUE_TYPE
*, const long *);
2928 encode_ieee_extended (const struct real_format
*fmt
, long *buf
,
2929 const REAL_VALUE_TYPE
*r
)
2931 unsigned long image_hi
, sig_hi
, sig_lo
;
2932 bool denormal
= (r
->sig
[SIGSZ
-1] & SIG_MSB
) == 0;
2934 image_hi
= r
->sign
<< 15;
2935 sig_hi
= sig_lo
= 0;
2947 /* Intel requires the explicit integer bit to be set, otherwise
2948 it considers the value a "pseudo-infinity". Motorola docs
2949 say it doesn't care. */
2950 sig_hi
= 0x80000000;
2955 sig_lo
= sig_hi
= 0xffffffff;
2963 if (HOST_BITS_PER_LONG
== 32)
2965 sig_hi
= r
->sig
[SIGSZ
-1];
2966 sig_lo
= r
->sig
[SIGSZ
-2];
2970 sig_lo
= r
->sig
[SIGSZ
-1];
2971 sig_hi
= sig_lo
>> 31 >> 1;
2972 sig_lo
&= 0xffffffff;
2974 if (r
->signalling
== fmt
->qnan_msb_set
)
2975 sig_hi
&= ~(1 << 30);
2978 if ((sig_hi
& 0x7fffffff) == 0 && sig_lo
== 0)
2981 /* Intel requires the explicit integer bit to be set, otherwise
2982 it considers the value a "pseudo-nan". Motorola docs say it
2984 sig_hi
|= 0x80000000;
2989 sig_lo
= sig_hi
= 0xffffffff;
2995 int exp
= REAL_EXP (r
);
2997 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2998 whereas the intermediate representation is 0.F x 2**exp.
2999 Which means we're off by one.
3001 Except for Motorola, which consider exp=0 and explicit
3002 integer bit set to continue to be normalized. In theory
3003 this discrepancy has been taken care of by the difference
3004 in fmt->emin in round_for_format. */
3016 if (HOST_BITS_PER_LONG
== 32)
3018 sig_hi
= r
->sig
[SIGSZ
-1];
3019 sig_lo
= r
->sig
[SIGSZ
-2];
3023 sig_lo
= r
->sig
[SIGSZ
-1];
3024 sig_hi
= sig_lo
>> 31 >> 1;
3025 sig_lo
&= 0xffffffff;
3034 if (FLOAT_WORDS_BIG_ENDIAN
)
3035 buf
[0] = image_hi
<< 16, buf
[1] = sig_hi
, buf
[2] = sig_lo
;
3037 buf
[0] = sig_lo
, buf
[1] = sig_hi
, buf
[2] = image_hi
;
3041 encode_ieee_extended_128 (const struct real_format
*fmt
, long *buf
,
3042 const REAL_VALUE_TYPE
*r
)
3044 buf
[3 * !FLOAT_WORDS_BIG_ENDIAN
] = 0;
3045 encode_ieee_extended (fmt
, buf
+!!FLOAT_WORDS_BIG_ENDIAN
, r
);
3049 decode_ieee_extended (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
,
3052 unsigned long image_hi
, sig_hi
, sig_lo
;
3056 if (FLOAT_WORDS_BIG_ENDIAN
)
3057 image_hi
= buf
[0] >> 16, sig_hi
= buf
[1], sig_lo
= buf
[2];
3059 sig_lo
= buf
[0], sig_hi
= buf
[1], image_hi
= buf
[2];
3060 sig_lo
&= 0xffffffff;
3061 sig_hi
&= 0xffffffff;
3062 image_hi
&= 0xffffffff;
3064 sign
= (image_hi
>> 15) & 1;
3065 exp
= image_hi
& 0x7fff;
3067 memset (r
, 0, sizeof (*r
));
3071 if ((sig_hi
|| sig_lo
) && fmt
->has_denorm
)
3073 r
->class = rvc_normal
;
3076 /* When the IEEE format contains a hidden bit, we know that
3077 it's zero at this point, and so shift up the significand
3078 and decrease the exponent to match. In this case, Motorola
3079 defines the explicit integer bit to be valid, so we don't
3080 know whether the msb is set or not. */
3081 SET_REAL_EXP (r
, fmt
->emin
);
3082 if (HOST_BITS_PER_LONG
== 32)
3084 r
->sig
[SIGSZ
-1] = sig_hi
;
3085 r
->sig
[SIGSZ
-2] = sig_lo
;
3088 r
->sig
[SIGSZ
-1] = (sig_hi
<< 31 << 1) | sig_lo
;
3092 else if (fmt
->has_signed_zero
)
3095 else if (exp
== 32767 && (fmt
->has_nans
|| fmt
->has_inf
))
3097 /* See above re "pseudo-infinities" and "pseudo-nans".
3098 Short summary is that the MSB will likely always be
3099 set, and that we don't care about it. */
3100 sig_hi
&= 0x7fffffff;
3102 if (sig_hi
|| sig_lo
)
3106 r
->signalling
= ((sig_hi
>> 30) & 1) ^ fmt
->qnan_msb_set
;
3107 if (HOST_BITS_PER_LONG
== 32)
3109 r
->sig
[SIGSZ
-1] = sig_hi
;
3110 r
->sig
[SIGSZ
-2] = sig_lo
;
3113 r
->sig
[SIGSZ
-1] = (sig_hi
<< 31 << 1) | sig_lo
;
3123 r
->class = rvc_normal
;
3125 SET_REAL_EXP (r
, exp
- 16383 + 1);
3126 if (HOST_BITS_PER_LONG
== 32)
3128 r
->sig
[SIGSZ
-1] = sig_hi
;
3129 r
->sig
[SIGSZ
-2] = sig_lo
;
3132 r
->sig
[SIGSZ
-1] = (sig_hi
<< 31 << 1) | sig_lo
;
3137 decode_ieee_extended_128 (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
,
3140 decode_ieee_extended (fmt
, r
, buf
+!!FLOAT_WORDS_BIG_ENDIAN
);
3143 const struct real_format ieee_extended_motorola_format
=
3145 encode_ieee_extended
,
3146 decode_ieee_extended
,
3161 const struct real_format ieee_extended_intel_96_format
=
3163 encode_ieee_extended
,
3164 decode_ieee_extended
,
3179 const struct real_format ieee_extended_intel_128_format
=
3181 encode_ieee_extended_128
,
3182 decode_ieee_extended_128
,
3197 /* The following caters to i386 systems that set the rounding precision
3198 to 53 bits instead of 64, e.g. FreeBSD. */
3199 const struct real_format ieee_extended_intel_96_round_53_format
=
3201 encode_ieee_extended
,
3202 decode_ieee_extended
,
3217 /* IBM 128-bit extended precision format: a pair of IEEE double precision
3218 numbers whose sum is equal to the extended precision value. The number
3219 with greater magnitude is first. This format has the same magnitude
3220 range as an IEEE double precision value, but effectively 106 bits of
3221 significand precision. Infinity and NaN are represented by their IEEE
3222 double precision value stored in the first number, the second number is
3223 ignored. Zeroes, Infinities, and NaNs are set in both doubles
3224 due to precedent. */
3226 static void encode_ibm_extended (const struct real_format
*fmt
,
3227 long *, const REAL_VALUE_TYPE
*);
3228 static void decode_ibm_extended (const struct real_format
*,
3229 REAL_VALUE_TYPE
*, const long *);
3232 encode_ibm_extended (const struct real_format
*fmt
, long *buf
,
3233 const REAL_VALUE_TYPE
*r
)
3235 REAL_VALUE_TYPE u
, normr
, v
;
3236 const struct real_format
*base_fmt
;
3238 base_fmt
= fmt
->qnan_msb_set
? &ieee_double_format
: &mips_double_format
;
3240 /* Renormlize R before doing any arithmetic on it. */
3242 if (normr
.class == rvc_normal
)
3245 /* u = IEEE double precision portion of significand. */
3247 round_for_format (base_fmt
, &u
);
3248 encode_ieee_double (base_fmt
, &buf
[0], &u
);
3250 if (u
.class == rvc_normal
)
3252 do_add (&v
, &normr
, &u
, 1);
3253 /* Call round_for_format since we might need to denormalize. */
3254 round_for_format (base_fmt
, &v
);
3255 encode_ieee_double (base_fmt
, &buf
[2], &v
);
3259 /* Inf, NaN, 0 are all representable as doubles, so the
3260 least-significant part can be 0.0. */
3267 decode_ibm_extended (const struct real_format
*fmt ATTRIBUTE_UNUSED
, REAL_VALUE_TYPE
*r
,
3270 REAL_VALUE_TYPE u
, v
;
3271 const struct real_format
*base_fmt
;
3273 base_fmt
= fmt
->qnan_msb_set
? &ieee_double_format
: &mips_double_format
;
3274 decode_ieee_double (base_fmt
, &u
, &buf
[0]);
3276 if (u
.class != rvc_zero
&& u
.class != rvc_inf
&& u
.class != rvc_nan
)
3278 decode_ieee_double (base_fmt
, &v
, &buf
[2]);
3279 do_add (r
, &u
, &v
, 0);
3285 const struct real_format ibm_extended_format
=
3287 encode_ibm_extended
,
3288 decode_ibm_extended
,
3303 const struct real_format mips_extended_format
=
3305 encode_ibm_extended
,
3306 decode_ibm_extended
,
3322 /* IEEE quad precision format. */
3324 static void encode_ieee_quad (const struct real_format
*fmt
,
3325 long *, const REAL_VALUE_TYPE
*);
3326 static void decode_ieee_quad (const struct real_format
*,
3327 REAL_VALUE_TYPE
*, const long *);
3330 encode_ieee_quad (const struct real_format
*fmt
, long *buf
,
3331 const REAL_VALUE_TYPE
*r
)
3333 unsigned long image3
, image2
, image1
, image0
, exp
;
3334 bool denormal
= (r
->sig
[SIGSZ
-1] & SIG_MSB
) == 0;
3337 image3
= r
->sign
<< 31;
3342 rshift_significand (&u
, r
, SIGNIFICAND_BITS
- 113);
3351 image3
|= 32767 << 16;
3354 image3
|= 0x7fffffff;
3355 image2
= 0xffffffff;
3356 image1
= 0xffffffff;
3357 image0
= 0xffffffff;
3364 image3
|= 32767 << 16;
3368 /* Don't use bits from the significand. The
3369 initialization above is right. */
3371 else if (HOST_BITS_PER_LONG
== 32)
3376 image3
|= u
.sig
[3] & 0xffff;
3381 image1
= image0
>> 31 >> 1;
3383 image3
|= (image2
>> 31 >> 1) & 0xffff;
3384 image0
&= 0xffffffff;
3385 image2
&= 0xffffffff;
3387 if (r
->signalling
== fmt
->qnan_msb_set
)
3391 /* We overload qnan_msb_set here: it's only clear for
3392 mips_ieee_single, which wants all mantissa bits but the
3393 quiet/signalling one set in canonical NaNs (at least
3395 if (r
->canonical
&& !fmt
->qnan_msb_set
)
3398 image2
= image1
= image0
= 0xffffffff;
3400 else if (((image3
& 0xffff) | image2
| image1
| image0
) == 0)
3405 image3
|= 0x7fffffff;
3406 image2
= 0xffffffff;
3407 image1
= 0xffffffff;
3408 image0
= 0xffffffff;
3413 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3414 whereas the intermediate representation is 0.F x 2**exp.
3415 Which means we're off by one. */
3419 exp
= REAL_EXP (r
) + 16383 - 1;
3420 image3
|= exp
<< 16;
3422 if (HOST_BITS_PER_LONG
== 32)
3427 image3
|= u
.sig
[3] & 0xffff;
3432 image1
= image0
>> 31 >> 1;
3434 image3
|= (image2
>> 31 >> 1) & 0xffff;
3435 image0
&= 0xffffffff;
3436 image2
&= 0xffffffff;
3444 if (FLOAT_WORDS_BIG_ENDIAN
)
3461 decode_ieee_quad (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
,
3464 unsigned long image3
, image2
, image1
, image0
;
3468 if (FLOAT_WORDS_BIG_ENDIAN
)
3482 image0
&= 0xffffffff;
3483 image1
&= 0xffffffff;
3484 image2
&= 0xffffffff;
3486 sign
= (image3
>> 31) & 1;
3487 exp
= (image3
>> 16) & 0x7fff;
3490 memset (r
, 0, sizeof (*r
));
3494 if ((image3
| image2
| image1
| image0
) && fmt
->has_denorm
)
3496 r
->class = rvc_normal
;
3499 SET_REAL_EXP (r
, -16382 + (SIGNIFICAND_BITS
- 112));
3500 if (HOST_BITS_PER_LONG
== 32)
3509 r
->sig
[0] = (image1
<< 31 << 1) | image0
;
3510 r
->sig
[1] = (image3
<< 31 << 1) | image2
;
3515 else if (fmt
->has_signed_zero
)
3518 else if (exp
== 32767 && (fmt
->has_nans
|| fmt
->has_inf
))
3520 if (image3
| image2
| image1
| image0
)
3524 r
->signalling
= ((image3
>> 15) & 1) ^ fmt
->qnan_msb_set
;
3526 if (HOST_BITS_PER_LONG
== 32)
3535 r
->sig
[0] = (image1
<< 31 << 1) | image0
;
3536 r
->sig
[1] = (image3
<< 31 << 1) | image2
;
3538 lshift_significand (r
, r
, SIGNIFICAND_BITS
- 113);
3548 r
->class = rvc_normal
;
3550 SET_REAL_EXP (r
, exp
- 16383 + 1);
3552 if (HOST_BITS_PER_LONG
== 32)
3561 r
->sig
[0] = (image1
<< 31 << 1) | image0
;
3562 r
->sig
[1] = (image3
<< 31 << 1) | image2
;
3564 lshift_significand (r
, r
, SIGNIFICAND_BITS
- 113);
3565 r
->sig
[SIGSZ
-1] |= SIG_MSB
;
3569 const struct real_format ieee_quad_format
=
3587 const struct real_format mips_quad_format
=
3605 /* Descriptions of VAX floating point formats can be found beginning at
3607 http://h71000.www7.hp.com/doc/73FINAL/4515/4515pro_013.html#f_floating_point_format
3609 The thing to remember is that they're almost IEEE, except for word
3610 order, exponent bias, and the lack of infinities, nans, and denormals.
3612 We don't implement the H_floating format here, simply because neither
3613 the VAX or Alpha ports use it. */
3615 static void encode_vax_f (const struct real_format
*fmt
,
3616 long *, const REAL_VALUE_TYPE
*);
3617 static void decode_vax_f (const struct real_format
*,
3618 REAL_VALUE_TYPE
*, const long *);
3619 static void encode_vax_d (const struct real_format
*fmt
,
3620 long *, const REAL_VALUE_TYPE
*);
3621 static void decode_vax_d (const struct real_format
*,
3622 REAL_VALUE_TYPE
*, const long *);
3623 static void encode_vax_g (const struct real_format
*fmt
,
3624 long *, const REAL_VALUE_TYPE
*);
3625 static void decode_vax_g (const struct real_format
*,
3626 REAL_VALUE_TYPE
*, const long *);
3629 encode_vax_f (const struct real_format
*fmt ATTRIBUTE_UNUSED
, long *buf
,
3630 const REAL_VALUE_TYPE
*r
)
3632 unsigned long sign
, exp
, sig
, image
;
3634 sign
= r
->sign
<< 15;
3644 image
= 0xffff7fff | sign
;
3648 sig
= (r
->sig
[SIGSZ
-1] >> (HOST_BITS_PER_LONG
- 24)) & 0x7fffff;
3649 exp
= REAL_EXP (r
) + 128;
3651 image
= (sig
<< 16) & 0xffff0000;
3665 decode_vax_f (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
3666 REAL_VALUE_TYPE
*r
, const long *buf
)
3668 unsigned long image
= buf
[0] & 0xffffffff;
3669 int exp
= (image
>> 7) & 0xff;
3671 memset (r
, 0, sizeof (*r
));
3675 r
->class = rvc_normal
;
3676 r
->sign
= (image
>> 15) & 1;
3677 SET_REAL_EXP (r
, exp
- 128);
3679 image
= ((image
& 0x7f) << 16) | ((image
>> 16) & 0xffff);
3680 r
->sig
[SIGSZ
-1] = (image
<< (HOST_BITS_PER_LONG
- 24)) | SIG_MSB
;
3685 encode_vax_d (const struct real_format
*fmt ATTRIBUTE_UNUSED
, long *buf
,
3686 const REAL_VALUE_TYPE
*r
)
3688 unsigned long image0
, image1
, sign
= r
->sign
<< 15;
3693 image0
= image1
= 0;
3698 image0
= 0xffff7fff | sign
;
3699 image1
= 0xffffffff;
3703 /* Extract the significand into straight hi:lo. */
3704 if (HOST_BITS_PER_LONG
== 64)
3706 image0
= r
->sig
[SIGSZ
-1];
3707 image1
= (image0
>> (64 - 56)) & 0xffffffff;
3708 image0
= (image0
>> (64 - 56 + 1) >> 31) & 0x7fffff;
3712 image0
= r
->sig
[SIGSZ
-1];
3713 image1
= r
->sig
[SIGSZ
-2];
3714 image1
= (image0
<< 24) | (image1
>> 8);
3715 image0
= (image0
>> 8) & 0xffffff;
3718 /* Rearrange the half-words of the significand to match the
3720 image0
= ((image0
<< 16) | (image0
>> 16)) & 0xffff007f;
3721 image1
= ((image1
<< 16) | (image1
>> 16)) & 0xffffffff;
3723 /* Add the sign and exponent. */
3725 image0
|= (REAL_EXP (r
) + 128) << 7;
3732 if (FLOAT_WORDS_BIG_ENDIAN
)
3733 buf
[0] = image1
, buf
[1] = image0
;
3735 buf
[0] = image0
, buf
[1] = image1
;
3739 decode_vax_d (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
3740 REAL_VALUE_TYPE
*r
, const long *buf
)
3742 unsigned long image0
, image1
;
3745 if (FLOAT_WORDS_BIG_ENDIAN
)
3746 image1
= buf
[0], image0
= buf
[1];
3748 image0
= buf
[0], image1
= buf
[1];
3749 image0
&= 0xffffffff;
3750 image1
&= 0xffffffff;
3752 exp
= (image0
>> 7) & 0xff;
3754 memset (r
, 0, sizeof (*r
));
3758 r
->class = rvc_normal
;
3759 r
->sign
= (image0
>> 15) & 1;
3760 SET_REAL_EXP (r
, exp
- 128);
3762 /* Rearrange the half-words of the external format into
3763 proper ascending order. */
3764 image0
= ((image0
& 0x7f) << 16) | ((image0
>> 16) & 0xffff);
3765 image1
= ((image1
& 0xffff) << 16) | ((image1
>> 16) & 0xffff);
3767 if (HOST_BITS_PER_LONG
== 64)
3769 image0
= (image0
<< 31 << 1) | image1
;
3772 r
->sig
[SIGSZ
-1] = image0
;
3776 r
->sig
[SIGSZ
-1] = image0
;
3777 r
->sig
[SIGSZ
-2] = image1
;
3778 lshift_significand (r
, r
, 2*HOST_BITS_PER_LONG
- 56);
3779 r
->sig
[SIGSZ
-1] |= SIG_MSB
;
3785 encode_vax_g (const struct real_format
*fmt ATTRIBUTE_UNUSED
, long *buf
,
3786 const REAL_VALUE_TYPE
*r
)
3788 unsigned long image0
, image1
, sign
= r
->sign
<< 15;
3793 image0
= image1
= 0;
3798 image0
= 0xffff7fff | sign
;
3799 image1
= 0xffffffff;
3803 /* Extract the significand into straight hi:lo. */
3804 if (HOST_BITS_PER_LONG
== 64)
3806 image0
= r
->sig
[SIGSZ
-1];
3807 image1
= (image0
>> (64 - 53)) & 0xffffffff;
3808 image0
= (image0
>> (64 - 53 + 1) >> 31) & 0xfffff;
3812 image0
= r
->sig
[SIGSZ
-1];
3813 image1
= r
->sig
[SIGSZ
-2];
3814 image1
= (image0
<< 21) | (image1
>> 11);
3815 image0
= (image0
>> 11) & 0xfffff;
3818 /* Rearrange the half-words of the significand to match the
3820 image0
= ((image0
<< 16) | (image0
>> 16)) & 0xffff000f;
3821 image1
= ((image1
<< 16) | (image1
>> 16)) & 0xffffffff;
3823 /* Add the sign and exponent. */
3825 image0
|= (REAL_EXP (r
) + 1024) << 4;
3832 if (FLOAT_WORDS_BIG_ENDIAN
)
3833 buf
[0] = image1
, buf
[1] = image0
;
3835 buf
[0] = image0
, buf
[1] = image1
;
3839 decode_vax_g (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
3840 REAL_VALUE_TYPE
*r
, const long *buf
)
3842 unsigned long image0
, image1
;
3845 if (FLOAT_WORDS_BIG_ENDIAN
)
3846 image1
= buf
[0], image0
= buf
[1];
3848 image0
= buf
[0], image1
= buf
[1];
3849 image0
&= 0xffffffff;
3850 image1
&= 0xffffffff;
3852 exp
= (image0
>> 4) & 0x7ff;
3854 memset (r
, 0, sizeof (*r
));
3858 r
->class = rvc_normal
;
3859 r
->sign
= (image0
>> 15) & 1;
3860 SET_REAL_EXP (r
, exp
- 1024);
3862 /* Rearrange the half-words of the external format into
3863 proper ascending order. */
3864 image0
= ((image0
& 0xf) << 16) | ((image0
>> 16) & 0xffff);
3865 image1
= ((image1
& 0xffff) << 16) | ((image1
>> 16) & 0xffff);
3867 if (HOST_BITS_PER_LONG
== 64)
3869 image0
= (image0
<< 31 << 1) | image1
;
3872 r
->sig
[SIGSZ
-1] = image0
;
3876 r
->sig
[SIGSZ
-1] = image0
;
3877 r
->sig
[SIGSZ
-2] = image1
;
3878 lshift_significand (r
, r
, 64 - 53);
3879 r
->sig
[SIGSZ
-1] |= SIG_MSB
;
3884 const struct real_format vax_f_format
=
3902 const struct real_format vax_d_format
=
3920 const struct real_format vax_g_format
=
3938 /* A good reference for these can be found in chapter 9 of
3939 "ESA/390 Principles of Operation", IBM document number SA22-7201-01.
3940 An on-line version can be found here:
3942 http://publibz.boulder.ibm.com/cgi-bin/bookmgr_OS390/BOOKS/DZ9AR001/9.1?DT=19930923083613
3945 static void encode_i370_single (const struct real_format
*fmt
,
3946 long *, const REAL_VALUE_TYPE
*);
3947 static void decode_i370_single (const struct real_format
*,
3948 REAL_VALUE_TYPE
*, const long *);
3949 static void encode_i370_double (const struct real_format
*fmt
,
3950 long *, const REAL_VALUE_TYPE
*);
3951 static void decode_i370_double (const struct real_format
*,
3952 REAL_VALUE_TYPE
*, const long *);
3955 encode_i370_single (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
3956 long *buf
, const REAL_VALUE_TYPE
*r
)
3958 unsigned long sign
, exp
, sig
, image
;
3960 sign
= r
->sign
<< 31;
3970 image
= 0x7fffffff | sign
;
3974 sig
= (r
->sig
[SIGSZ
-1] >> (HOST_BITS_PER_LONG
- 24)) & 0xffffff;
3975 exp
= ((REAL_EXP (r
) / 4) + 64) << 24;
3976 image
= sign
| exp
| sig
;
3987 decode_i370_single (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
3988 REAL_VALUE_TYPE
*r
, const long *buf
)
3990 unsigned long sign
, sig
, image
= buf
[0];
3993 sign
= (image
>> 31) & 1;
3994 exp
= (image
>> 24) & 0x7f;
3995 sig
= image
& 0xffffff;
3997 memset (r
, 0, sizeof (*r
));
4001 r
->class = rvc_normal
;
4003 SET_REAL_EXP (r
, (exp
- 64) * 4);
4004 r
->sig
[SIGSZ
-1] = sig
<< (HOST_BITS_PER_LONG
- 24);
4010 encode_i370_double (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4011 long *buf
, const REAL_VALUE_TYPE
*r
)
4013 unsigned long sign
, exp
, image_hi
, image_lo
;
4015 sign
= r
->sign
<< 31;
4020 image_hi
= image_lo
= 0;
4025 image_hi
= 0x7fffffff | sign
;
4026 image_lo
= 0xffffffff;
4030 if (HOST_BITS_PER_LONG
== 64)
4032 image_hi
= r
->sig
[SIGSZ
-1];
4033 image_lo
= (image_hi
>> (64 - 56)) & 0xffffffff;
4034 image_hi
= (image_hi
>> (64 - 56 + 1) >> 31) & 0xffffff;
4038 image_hi
= r
->sig
[SIGSZ
-1];
4039 image_lo
= r
->sig
[SIGSZ
-2];
4040 image_lo
= (image_lo
>> 8) | (image_hi
<< 24);
4044 exp
= ((REAL_EXP (r
) / 4) + 64) << 24;
4045 image_hi
|= sign
| exp
;
4052 if (FLOAT_WORDS_BIG_ENDIAN
)
4053 buf
[0] = image_hi
, buf
[1] = image_lo
;
4055 buf
[0] = image_lo
, buf
[1] = image_hi
;
4059 decode_i370_double (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4060 REAL_VALUE_TYPE
*r
, const long *buf
)
4062 unsigned long sign
, image_hi
, image_lo
;
4065 if (FLOAT_WORDS_BIG_ENDIAN
)
4066 image_hi
= buf
[0], image_lo
= buf
[1];
4068 image_lo
= buf
[0], image_hi
= buf
[1];
4070 sign
= (image_hi
>> 31) & 1;
4071 exp
= (image_hi
>> 24) & 0x7f;
4072 image_hi
&= 0xffffff;
4073 image_lo
&= 0xffffffff;
4075 memset (r
, 0, sizeof (*r
));
4077 if (exp
|| image_hi
|| image_lo
)
4079 r
->class = rvc_normal
;
4081 SET_REAL_EXP (r
, (exp
- 64) * 4 + (SIGNIFICAND_BITS
- 56));
4083 if (HOST_BITS_PER_LONG
== 32)
4085 r
->sig
[0] = image_lo
;
4086 r
->sig
[1] = image_hi
;
4089 r
->sig
[0] = image_lo
| (image_hi
<< 31 << 1);
4095 const struct real_format i370_single_format
=
4108 false, /* ??? The encoding does allow for "unnormals". */
4109 false, /* ??? The encoding does allow for "unnormals". */
4113 const struct real_format i370_double_format
=
4126 false, /* ??? The encoding does allow for "unnormals". */
4127 false, /* ??? The encoding does allow for "unnormals". */
4131 /* The "twos-complement" c4x format is officially defined as
4135 This is rather misleading. One must remember that F is signed.
4136 A better description would be
4138 x = -1**s * ((s + 1 + .f) * 2**e
4140 So if we have a (4 bit) fraction of .1000 with a sign bit of 1,
4141 that's -1 * (1+1+(-.5)) == -1.5. I think.
4143 The constructions here are taken from Tables 5-1 and 5-2 of the
4144 TMS320C4x User's Guide wherein step-by-step instructions for
4145 conversion from IEEE are presented. That's close enough to our
4146 internal representation so as to make things easy.
4148 See http://www-s.ti.com/sc/psheets/spru063c/spru063c.pdf */
4150 static void encode_c4x_single (const struct real_format
*fmt
,
4151 long *, const REAL_VALUE_TYPE
*);
4152 static void decode_c4x_single (const struct real_format
*,
4153 REAL_VALUE_TYPE
*, const long *);
4154 static void encode_c4x_extended (const struct real_format
*fmt
,
4155 long *, const REAL_VALUE_TYPE
*);
4156 static void decode_c4x_extended (const struct real_format
*,
4157 REAL_VALUE_TYPE
*, const long *);
4160 encode_c4x_single (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4161 long *buf
, const REAL_VALUE_TYPE
*r
)
4163 unsigned long image
, exp
, sig
;
4175 sig
= 0x800000 - r
->sign
;
4179 exp
= REAL_EXP (r
) - 1;
4180 sig
= (r
->sig
[SIGSZ
-1] >> (HOST_BITS_PER_LONG
- 24)) & 0x7fffff;
4195 image
= ((exp
& 0xff) << 24) | (sig
& 0xffffff);
4200 decode_c4x_single (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4201 REAL_VALUE_TYPE
*r
, const long *buf
)
4203 unsigned long image
= buf
[0];
4207 exp
= (((image
>> 24) & 0xff) ^ 0x80) - 0x80;
4208 sf
= ((image
& 0xffffff) ^ 0x800000) - 0x800000;
4210 memset (r
, 0, sizeof (*r
));
4214 r
->class = rvc_normal
;
4216 sig
= sf
& 0x7fffff;
4225 sig
= (sig
<< (HOST_BITS_PER_LONG
- 24)) | SIG_MSB
;
4227 SET_REAL_EXP (r
, exp
+ 1);
4228 r
->sig
[SIGSZ
-1] = sig
;
4233 encode_c4x_extended (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4234 long *buf
, const REAL_VALUE_TYPE
*r
)
4236 unsigned long exp
, sig
;
4248 sig
= 0x80000000 - r
->sign
;
4252 exp
= REAL_EXP (r
) - 1;
4254 sig
= r
->sig
[SIGSZ
-1];
4255 if (HOST_BITS_PER_LONG
== 64)
4256 sig
= sig
>> 1 >> 31;
4273 exp
= (exp
& 0xff) << 24;
4276 if (FLOAT_WORDS_BIG_ENDIAN
)
4277 buf
[0] = exp
, buf
[1] = sig
;
4279 buf
[0] = sig
, buf
[0] = exp
;
4283 decode_c4x_extended (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4284 REAL_VALUE_TYPE
*r
, const long *buf
)
4289 if (FLOAT_WORDS_BIG_ENDIAN
)
4290 exp
= buf
[0], sf
= buf
[1];
4292 sf
= buf
[0], exp
= buf
[1];
4294 exp
= (((exp
>> 24) & 0xff) & 0x80) - 0x80;
4295 sf
= ((sf
& 0xffffffff) ^ 0x80000000) - 0x80000000;
4297 memset (r
, 0, sizeof (*r
));
4301 r
->class = rvc_normal
;
4303 sig
= sf
& 0x7fffffff;
4312 if (HOST_BITS_PER_LONG
== 64)
4313 sig
= sig
<< 1 << 31;
4316 SET_REAL_EXP (r
, exp
+ 1);
4317 r
->sig
[SIGSZ
-1] = sig
;
4321 const struct real_format c4x_single_format
=
4339 const struct real_format c4x_extended_format
=
4341 encode_c4x_extended
,
4342 decode_c4x_extended
,
4358 /* A synthetic "format" for internal arithmetic. It's the size of the
4359 internal significand minus the two bits needed for proper rounding.
4360 The encode and decode routines exist only to satisfy our paranoia
4363 static void encode_internal (const struct real_format
*fmt
,
4364 long *, const REAL_VALUE_TYPE
*);
4365 static void decode_internal (const struct real_format
*,
4366 REAL_VALUE_TYPE
*, const long *);
4369 encode_internal (const struct real_format
*fmt ATTRIBUTE_UNUSED
, long *buf
,
4370 const REAL_VALUE_TYPE
*r
)
4372 memcpy (buf
, r
, sizeof (*r
));
4376 decode_internal (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4377 REAL_VALUE_TYPE
*r
, const long *buf
)
4379 memcpy (r
, buf
, sizeof (*r
));
4382 const struct real_format real_internal_format
=
4388 SIGNIFICAND_BITS
- 2,
4389 SIGNIFICAND_BITS
- 2,
4400 /* Calculate the square root of X in mode MODE, and store the result
4401 in R. Return TRUE if the operation does not raise an exception.
4402 For details see "High Precision Division and Square Root",
4403 Alan H. Karp and Peter Markstein, HP Lab Report 93-93-42, June
4404 1993. http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf. */
4407 real_sqrt (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4408 const REAL_VALUE_TYPE
*x
)
4410 static REAL_VALUE_TYPE halfthree
;
4411 static bool init
= false;
4412 REAL_VALUE_TYPE h
, t
, i
;
4415 /* sqrt(-0.0) is -0.0. */
4416 if (real_isnegzero (x
))
4422 /* Negative arguments return NaN. */
4425 get_canonical_qnan (r
, 0);
4429 /* Infinity and NaN return themselves. */
4430 if (real_isinf (x
) || real_isnan (x
))
4438 do_add (&halfthree
, &dconst1
, &dconsthalf
, 0);
4442 /* Initial guess for reciprocal sqrt, i. */
4443 exp
= real_exponent (x
);
4444 real_ldexp (&i
, &dconst1
, -exp
/2);
4446 /* Newton's iteration for reciprocal sqrt, i. */
4447 for (iter
= 0; iter
< 16; iter
++)
4449 /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x). */
4450 do_multiply (&t
, x
, &i
);
4451 do_multiply (&h
, &t
, &i
);
4452 do_multiply (&t
, &h
, &dconsthalf
);
4453 do_add (&h
, &halfthree
, &t
, 1);
4454 do_multiply (&t
, &i
, &h
);
4456 /* Check for early convergence. */
4457 if (iter
>= 6 && real_identical (&i
, &t
))
4460 /* ??? Unroll loop to avoid copying. */
4464 /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)). */
4465 do_multiply (&t
, x
, &i
);
4466 do_multiply (&h
, &t
, &i
);
4467 do_add (&i
, &dconst1
, &h
, 1);
4468 do_multiply (&h
, &t
, &i
);
4469 do_multiply (&i
, &dconsthalf
, &h
);
4470 do_add (&h
, &t
, &i
, 0);
4472 /* ??? We need a Tuckerman test to get the last bit. */
4474 real_convert (r
, mode
, &h
);
4478 /* Calculate X raised to the integer exponent N in mode MODE and store
4479 the result in R. Return true if the result may be inexact due to
4480 loss of precision. The algorithm is the classic "left-to-right binary
4481 method" described in section 4.6.3 of Donald Knuth's "Seminumerical
4482 Algorithms", "The Art of Computer Programming", Volume 2. */
4485 real_powi (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4486 const REAL_VALUE_TYPE
*x
, HOST_WIDE_INT n
)
4488 unsigned HOST_WIDE_INT bit
;
4490 bool inexact
= false;
4502 /* Don't worry about overflow, from now on n is unsigned. */
4510 bit
= (unsigned HOST_WIDE_INT
) 1 << (HOST_BITS_PER_WIDE_INT
- 1);
4511 for (i
= 0; i
< HOST_BITS_PER_WIDE_INT
; i
++)
4515 inexact
|= do_multiply (&t
, &t
, &t
);
4517 inexact
|= do_multiply (&t
, &t
, x
);
4525 inexact
|= do_divide (&t
, &dconst1
, &t
);
4527 real_convert (r
, mode
, &t
);
4531 /* Round X to the nearest integer not larger in absolute value, i.e.
4532 towards zero, placing the result in R in mode MODE. */
4535 real_trunc (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4536 const REAL_VALUE_TYPE
*x
)
4538 do_fix_trunc (r
, x
);
4539 if (mode
!= VOIDmode
)
4540 real_convert (r
, mode
, r
);
4543 /* Round X to the largest integer not greater in value, i.e. round
4544 down, placing the result in R in mode MODE. */
4547 real_floor (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4548 const REAL_VALUE_TYPE
*x
)
4552 do_fix_trunc (&t
, x
);
4553 if (! real_identical (&t
, x
) && x
->sign
)
4554 do_add (&t
, &t
, &dconstm1
, 0);
4555 if (mode
!= VOIDmode
)
4556 real_convert (r
, mode
, &t
);
4559 /* Round X to the smallest integer not less then argument, i.e. round
4560 up, placing the result in R in mode MODE. */
4563 real_ceil (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4564 const REAL_VALUE_TYPE
*x
)
4568 do_fix_trunc (&t
, x
);
4569 if (! real_identical (&t
, x
) && ! x
->sign
)
4570 do_add (&t
, &t
, &dconst1
, 0);
4571 if (mode
!= VOIDmode
)
4572 real_convert (r
, mode
, &t
);
4575 /* Round X to the nearest integer, but round halfway cases away from
4579 real_round (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4580 const REAL_VALUE_TYPE
*x
)
4582 do_add (r
, x
, &dconsthalf
, x
->sign
);
4583 do_fix_trunc (r
, r
);
4584 if (mode
!= VOIDmode
)
4585 real_convert (r
, mode
, r
);