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
= r
->exp
- shift
;
508 get_inf (r
, r
->sign
);
509 else if (exp
< -MAX_EXP
)
510 get_zero (r
, r
->sign
);
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
= a
->exp
- b
->exp
;
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
;
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
= (a
->exp
- (2*SIGSZ
-1-i
)*(HOST_BITS_PER_LONG
/2)
752 + (b
->exp
- (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
;
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
= a
->exp
- b
->exp
+ 1;
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
;
940 else if (a
->exp
< b
->exp
)
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
)
964 get_zero (r
, r
->sign
);
965 else if (r
->exp
< SIGNIFICAND_BITS
)
966 clear_significand_below (r
, SIGNIFICAND_BITS
- r
->exp
);
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;
1108 /* R = OP0 * 2**EXP. */
1111 real_ldexp (REAL_VALUE_TYPE
*r
, const REAL_VALUE_TYPE
*op0
, int exp
)
1124 get_inf (r
, r
->sign
);
1125 else if (exp
< -MAX_EXP
)
1126 get_zero (r
, r
->sign
);
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 (a
->exp
!= b
->exp
)
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);
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 (r
->exp
> 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
- r
->exp
;
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
= num
->exp
, expd
= den
->exp
;
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
);
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
= r
.exp
* 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 u
.exp
= 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 u
.exp
+= r
.exp
- (SIGNIFICAND_BITS
- 1);
1514 /* Find power of 10. Do this by dividing out 10**2**M when
1515 this is larger than the current remainder. Fill PTEN with
1516 the power of 10 that we compute. */
1519 m
= floor_log2 ((int)(r
.exp
* M_LOG10_2
)) + 1;
1522 const REAL_VALUE_TYPE
*ptentwo
= ten_to_ptwo (m
);
1523 if (do_compare (&u
, ptentwo
, 0) >= 0)
1525 do_divide (&u
, &u
, ptentwo
);
1526 do_multiply (&pten
, &pten
, ptentwo
);
1533 /* We managed to divide off enough tens in the above reduction
1534 loop that we've now got a negative exponent. Fall into the
1535 less-than-one code to compute the proper value for PTEN. */
1542 /* Number is less than one. Pad significand with leading
1548 /* Stop if we'd shift bits off the bottom. */
1552 do_multiply (&u
, &v
, ten
);
1554 /* Stop if we're now >= 1. */
1563 /* Find power of 10. Do this by multiplying in P=10**2**M when
1564 the current remainder is smaller than 1/P. Fill PTEN with the
1565 power of 10 that we compute. */
1566 m
= floor_log2 ((int)(-r
.exp
* M_LOG10_2
)) + 1;
1569 const REAL_VALUE_TYPE
*ptentwo
= ten_to_ptwo (m
);
1570 const REAL_VALUE_TYPE
*ptenmtwo
= ten_to_mptwo (m
);
1572 if (do_compare (&v
, ptenmtwo
, 0) <= 0)
1574 do_multiply (&v
, &v
, ptentwo
);
1575 do_multiply (&pten
, &pten
, ptentwo
);
1581 /* Invert the positive power of 10 that we've collected so far. */
1582 do_divide (&pten
, one
, &pten
);
1590 /* At this point, PTEN should contain the nearest power of 10 smaller
1591 than R, such that this division produces the first digit.
1593 Using a divide-step primitive that returns the complete integral
1594 remainder avoids the rounding error that would be produced if
1595 we were to use do_divide here and then simply multiply by 10 for
1596 each subsequent digit. */
1598 digit
= rtd_divmod (&r
, &pten
);
1600 /* Be prepared for error in that division via underflow ... */
1601 if (digit
== 0 && cmp_significand_0 (&r
))
1603 /* Multiply by 10 and try again. */
1604 do_multiply (&r
, &r
, ten
);
1605 digit
= rtd_divmod (&r
, &pten
);
1611 /* ... or overflow. */
1619 else if (digit
> 10)
1624 /* Generate subsequent digits. */
1625 while (--digits
> 0)
1627 do_multiply (&r
, &r
, ten
);
1628 digit
= rtd_divmod (&r
, &pten
);
1633 /* Generate one more digit with which to do rounding. */
1634 do_multiply (&r
, &r
, ten
);
1635 digit
= rtd_divmod (&r
, &pten
);
1637 /* Round the result. */
1640 /* Round to nearest. If R is nonzero there are additional
1641 nonzero digits to be extracted. */
1642 if (cmp_significand_0 (&r
))
1644 /* Round to even. */
1645 else if ((p
[-1] - '0') & 1)
1662 /* Carry out of the first digit. This means we had all 9's and
1663 now have all 0's. "Prepend" a 1 by overwriting the first 0. */
1671 /* Insert the decimal point. */
1672 first
[0] = first
[1];
1675 /* If requested, drop trailing zeros. Never crop past "1.0". */
1676 if (crop_trailing_zeros
)
1677 while (last
> first
+ 3 && last
[-1] == '0')
1680 /* Append the exponent. */
1681 sprintf (last
, "e%+d", dec_exp
);
1684 /* Render R as a hexadecimal floating point constant. Emit DIGITS
1685 significant digits in the result, bounded by BUF_SIZE. If DIGITS is 0,
1686 choose the maximum for the representation. If CROP_TRAILING_ZEROS,
1687 strip trailing zeros. */
1690 real_to_hexadecimal (char *str
, const REAL_VALUE_TYPE
*r
, size_t buf_size
,
1691 size_t digits
, int crop_trailing_zeros
)
1693 int i
, j
, exp
= r
->exp
;
1706 strcpy (str
, (r
->sign
? "-Inf" : "+Inf"));
1709 /* ??? Print the significand as well, if not canonical? */
1710 strcpy (str
, (r
->sign
? "-NaN" : "+NaN"));
1717 digits
= SIGNIFICAND_BITS
/ 4;
1719 /* Bound the number of digits printed by the size of the output buffer. */
1721 sprintf (exp_buf
, "p%+d", exp
);
1722 max_digits
= buf_size
- strlen (exp_buf
) - r
->sign
- 4 - 1;
1723 if (max_digits
> buf_size
)
1725 if (digits
> max_digits
)
1726 digits
= max_digits
;
1737 for (i
= SIGSZ
- 1; i
>= 0; --i
)
1738 for (j
= HOST_BITS_PER_LONG
- 4; j
>= 0; j
-= 4)
1740 *p
++ = "0123456789abcdef"[(r
->sig
[i
] >> j
) & 15];
1746 if (crop_trailing_zeros
)
1747 while (p
> first
+ 1 && p
[-1] == '0')
1750 sprintf (p
, "p%+d", exp
);
1753 /* Initialize R from a decimal or hexadecimal string. The string is
1754 assumed to have been syntax checked already. */
1757 real_from_string (REAL_VALUE_TYPE
*r
, const char *str
)
1769 else if (*str
== '+')
1772 if (str
[0] == '0' && (str
[1] == 'x' || str
[1] == 'X'))
1774 /* Hexadecimal floating point. */
1775 int pos
= SIGNIFICAND_BITS
- 4, d
;
1783 d
= hex_value (*str
);
1788 r
->sig
[pos
/ HOST_BITS_PER_LONG
]
1789 |= (unsigned long) d
<< (pos
% HOST_BITS_PER_LONG
);
1798 if (pos
== SIGNIFICAND_BITS
- 4)
1805 d
= hex_value (*str
);
1810 r
->sig
[pos
/ HOST_BITS_PER_LONG
]
1811 |= (unsigned long) d
<< (pos
% HOST_BITS_PER_LONG
);
1817 if (*str
== 'p' || *str
== 'P')
1819 bool exp_neg
= false;
1827 else if (*str
== '+')
1831 while (ISDIGIT (*str
))
1837 /* Overflowed the exponent. */
1851 r
->class = rvc_normal
;
1858 /* Decimal floating point. */
1859 const REAL_VALUE_TYPE
*ten
= ten_to_ptwo (0);
1864 while (ISDIGIT (*str
))
1867 do_multiply (r
, r
, ten
);
1869 do_add (r
, r
, real_digit (d
), 0);
1874 if (r
->class == rvc_zero
)
1879 while (ISDIGIT (*str
))
1882 do_multiply (r
, r
, ten
);
1884 do_add (r
, r
, real_digit (d
), 0);
1889 if (*str
== 'e' || *str
== 'E')
1891 bool exp_neg
= false;
1899 else if (*str
== '+')
1903 while (ISDIGIT (*str
))
1909 /* Overflowed the exponent. */
1923 times_pten (r
, exp
);
1938 /* Legacy. Similar, but return the result directly. */
1941 real_from_string2 (const char *s
, enum machine_mode mode
)
1945 real_from_string (&r
, s
);
1946 if (mode
!= VOIDmode
)
1947 real_convert (&r
, mode
, &r
);
1952 /* Initialize R from the integer pair HIGH+LOW. */
1955 real_from_integer (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
1956 unsigned HOST_WIDE_INT low
, HOST_WIDE_INT high
,
1959 if (low
== 0 && high
== 0)
1963 r
->class = rvc_normal
;
1964 r
->sign
= high
< 0 && !unsigned_p
;
1965 r
->exp
= 2 * HOST_BITS_PER_WIDE_INT
;
1976 if (HOST_BITS_PER_LONG
== HOST_BITS_PER_WIDE_INT
)
1978 r
->sig
[SIGSZ
-1] = high
;
1979 r
->sig
[SIGSZ
-2] = low
;
1980 memset (r
->sig
, 0, sizeof(long)*(SIGSZ
-2));
1982 else if (HOST_BITS_PER_LONG
*2 == HOST_BITS_PER_WIDE_INT
)
1984 r
->sig
[SIGSZ
-1] = high
>> (HOST_BITS_PER_LONG
- 1) >> 1;
1985 r
->sig
[SIGSZ
-2] = high
;
1986 r
->sig
[SIGSZ
-3] = low
>> (HOST_BITS_PER_LONG
- 1) >> 1;
1987 r
->sig
[SIGSZ
-4] = low
;
1989 memset (r
->sig
, 0, sizeof(long)*(SIGSZ
-4));
1997 if (mode
!= VOIDmode
)
1998 real_convert (r
, mode
, r
);
2001 /* Returns 10**2**N. */
2003 static const REAL_VALUE_TYPE
*
2006 static REAL_VALUE_TYPE tens
[EXP_BITS
];
2008 if (n
< 0 || n
>= EXP_BITS
)
2011 if (tens
[n
].class == rvc_zero
)
2013 if (n
< (HOST_BITS_PER_WIDE_INT
== 64 ? 5 : 4))
2015 HOST_WIDE_INT t
= 10;
2018 for (i
= 0; i
< n
; ++i
)
2021 real_from_integer (&tens
[n
], VOIDmode
, t
, 0, 1);
2025 const REAL_VALUE_TYPE
*t
= ten_to_ptwo (n
- 1);
2026 do_multiply (&tens
[n
], t
, t
);
2033 /* Returns 10**(-2**N). */
2035 static const REAL_VALUE_TYPE
*
2036 ten_to_mptwo (int n
)
2038 static REAL_VALUE_TYPE tens
[EXP_BITS
];
2040 if (n
< 0 || n
>= EXP_BITS
)
2043 if (tens
[n
].class == rvc_zero
)
2044 do_divide (&tens
[n
], real_digit (1), ten_to_ptwo (n
));
2051 static const REAL_VALUE_TYPE
*
2054 static REAL_VALUE_TYPE num
[10];
2059 if (n
> 0 && num
[n
].class == rvc_zero
)
2060 real_from_integer (&num
[n
], VOIDmode
, n
, 0, 1);
2065 /* Multiply R by 10**EXP. */
2068 times_pten (REAL_VALUE_TYPE
*r
, int exp
)
2070 REAL_VALUE_TYPE pten
, *rr
;
2071 bool negative
= (exp
< 0);
2077 pten
= *real_digit (1);
2083 for (i
= 0; exp
> 0; ++i
, exp
>>= 1)
2085 do_multiply (rr
, rr
, ten_to_ptwo (i
));
2088 do_divide (r
, r
, &pten
);
2091 /* Fills R with +Inf. */
2094 real_inf (REAL_VALUE_TYPE
*r
)
2099 /* Fills R with a NaN whose significand is described by STR. If QUIET,
2100 we force a QNaN, else we force an SNaN. The string, if not empty,
2101 is parsed as a number and placed in the significand. Return true
2102 if the string was successfully parsed. */
2105 real_nan (REAL_VALUE_TYPE
*r
, const char *str
, int quiet
,
2106 enum machine_mode mode
)
2108 const struct real_format
*fmt
;
2110 fmt
= REAL_MODE_FORMAT (mode
);
2117 get_canonical_qnan (r
, 0);
2119 get_canonical_snan (r
, 0);
2126 memset (r
, 0, sizeof (*r
));
2129 /* Parse akin to strtol into the significand of R. */
2131 while (ISSPACE (*str
))
2135 else if (*str
== '+')
2145 while ((d
= hex_value (*str
)) < base
)
2152 lshift_significand (r
, r
, 3);
2155 lshift_significand (r
, r
, 4);
2158 lshift_significand_1 (&u
, r
);
2159 lshift_significand (r
, r
, 3);
2160 add_significands (r
, r
, &u
);
2168 add_significands (r
, r
, &u
);
2173 /* Must have consumed the entire string for success. */
2177 /* Shift the significand into place such that the bits
2178 are in the most significant bits for the format. */
2179 lshift_significand (r
, r
, SIGNIFICAND_BITS
- fmt
->pnan
);
2181 /* Our MSB is always unset for NaNs. */
2182 r
->sig
[SIGSZ
-1] &= ~SIG_MSB
;
2184 /* Force quiet or signalling NaN. */
2185 r
->signalling
= !quiet
;
2191 /* Fills R with the largest finite value representable in mode MODE.
2192 If SIGN is nonzero, R is set to the most negative finite value. */
2195 real_maxval (REAL_VALUE_TYPE
*r
, int sign
, enum machine_mode mode
)
2197 const struct real_format
*fmt
;
2200 fmt
= REAL_MODE_FORMAT (mode
);
2204 r
->class = rvc_normal
;
2208 r
->exp
= fmt
->emax
* fmt
->log2_b
;
2210 np2
= SIGNIFICAND_BITS
- fmt
->p
* fmt
->log2_b
;
2211 memset (r
->sig
, -1, SIGSZ
* sizeof (unsigned long));
2212 clear_significand_below (r
, np2
);
2215 /* Fills R with 2**N. */
2218 real_2expN (REAL_VALUE_TYPE
*r
, int n
)
2220 memset (r
, 0, sizeof (*r
));
2225 else if (n
< -MAX_EXP
)
2229 r
->class = rvc_normal
;
2231 r
->sig
[SIGSZ
-1] = SIG_MSB
;
2237 round_for_format (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
)
2240 unsigned long sticky
;
2244 p2
= fmt
->p
* fmt
->log2_b
;
2245 emin2m1
= (fmt
->emin
- 1) * fmt
->log2_b
;
2246 emax2
= fmt
->emax
* fmt
->log2_b
;
2248 np2
= SIGNIFICAND_BITS
- p2
;
2252 get_zero (r
, r
->sign
);
2254 if (!fmt
->has_signed_zero
)
2259 get_inf (r
, r
->sign
);
2264 clear_significand_below (r
, np2
);
2274 /* If we're not base2, normalize the exponent to a multiple of
2276 if (fmt
->log2_b
!= 1)
2278 int shift
= r
->exp
& (fmt
->log2_b
- 1);
2281 shift
= fmt
->log2_b
- shift
;
2282 r
->sig
[0] |= sticky_rshift_significand (r
, r
, shift
);
2287 /* Check the range of the exponent. If we're out of range,
2288 either underflow or overflow. */
2291 else if (r
->exp
<= emin2m1
)
2295 if (!fmt
->has_denorm
)
2297 /* Don't underflow completely until we've had a chance to round. */
2298 if (r
->exp
< emin2m1
)
2303 diff
= emin2m1
- r
->exp
+ 1;
2307 /* De-normalize the significand. */
2308 r
->sig
[0] |= sticky_rshift_significand (r
, r
, diff
);
2313 /* There are P2 true significand bits, followed by one guard bit,
2314 followed by one sticky bit, followed by stuff. Fold nonzero
2315 stuff into the sticky bit. */
2318 for (i
= 0, w
= (np2
- 1) / HOST_BITS_PER_LONG
; i
< w
; ++i
)
2319 sticky
|= r
->sig
[i
];
2321 r
->sig
[w
] & (((unsigned long)1 << ((np2
- 1) % HOST_BITS_PER_LONG
)) - 1);
2323 guard
= test_significand_bit (r
, np2
- 1);
2324 lsb
= test_significand_bit (r
, np2
);
2326 /* Round to even. */
2327 if (guard
&& (sticky
|| lsb
))
2331 set_significand_bit (&u
, np2
);
2333 if (add_significands (r
, r
, &u
))
2335 /* Overflow. Means the significand had been all ones, and
2336 is now all zeros. Need to increase the exponent, and
2337 possibly re-normalize it. */
2338 if (++r
->exp
> emax2
)
2340 r
->sig
[SIGSZ
-1] = SIG_MSB
;
2342 if (fmt
->log2_b
!= 1)
2344 int shift
= r
->exp
& (fmt
->log2_b
- 1);
2347 shift
= fmt
->log2_b
- shift
;
2348 rshift_significand (r
, r
, shift
);
2357 /* Catch underflow that we deferred until after rounding. */
2358 if (r
->exp
<= emin2m1
)
2361 /* Clear out trailing garbage. */
2362 clear_significand_below (r
, np2
);
2365 /* Extend or truncate to a new mode. */
2368 real_convert (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
2369 const REAL_VALUE_TYPE
*a
)
2371 const struct real_format
*fmt
;
2373 fmt
= REAL_MODE_FORMAT (mode
);
2378 round_for_format (fmt
, r
);
2380 /* round_for_format de-normalizes denormals. Undo just that part. */
2381 if (r
->class == rvc_normal
)
2385 /* Legacy. Likewise, except return the struct directly. */
2388 real_value_truncate (enum machine_mode mode
, REAL_VALUE_TYPE a
)
2391 real_convert (&r
, mode
, &a
);
2395 /* Return true if truncating to MODE is exact. */
2398 exact_real_truncate (enum machine_mode mode
, const REAL_VALUE_TYPE
*a
)
2401 real_convert (&t
, mode
, a
);
2402 return real_identical (&t
, a
);
2405 /* Write R to the given target format. Place the words of the result
2406 in target word order in BUF. There are always 32 bits in each
2407 long, no matter the size of the host long.
2409 Legacy: return word 0 for implementing REAL_VALUE_TO_TARGET_SINGLE. */
2412 real_to_target_fmt (long *buf
, const REAL_VALUE_TYPE
*r_orig
,
2413 const struct real_format
*fmt
)
2419 round_for_format (fmt
, &r
);
2423 (*fmt
->encode
) (fmt
, buf
, &r
);
2428 /* Similar, but look up the format from MODE. */
2431 real_to_target (long *buf
, const REAL_VALUE_TYPE
*r
, enum machine_mode mode
)
2433 const struct real_format
*fmt
;
2435 fmt
= REAL_MODE_FORMAT (mode
);
2439 return real_to_target_fmt (buf
, r
, fmt
);
2442 /* Read R from the given target format. Read the words of the result
2443 in target word order in BUF. There are always 32 bits in each
2444 long, no matter the size of the host long. */
2447 real_from_target_fmt (REAL_VALUE_TYPE
*r
, const long *buf
,
2448 const struct real_format
*fmt
)
2450 (*fmt
->decode
) (fmt
, r
, buf
);
2453 /* Similar, but look up the format from MODE. */
2456 real_from_target (REAL_VALUE_TYPE
*r
, const long *buf
, enum machine_mode mode
)
2458 const struct real_format
*fmt
;
2460 fmt
= REAL_MODE_FORMAT (mode
);
2464 (*fmt
->decode
) (fmt
, r
, buf
);
2467 /* Return the number of bits in the significand for MODE. */
2468 /* ??? Legacy. Should get access to real_format directly. */
2471 significand_size (enum machine_mode mode
)
2473 const struct real_format
*fmt
;
2475 fmt
= REAL_MODE_FORMAT (mode
);
2479 return fmt
->p
* fmt
->log2_b
;
2482 /* Return a hash value for the given real value. */
2483 /* ??? The "unsigned int" return value is intended to be hashval_t,
2484 but I didn't want to pull hashtab.h into real.h. */
2487 real_hash (const REAL_VALUE_TYPE
*r
)
2492 h
= r
->class | (r
->sign
<< 2);
2505 h
^= (unsigned int)-1;
2514 if (sizeof(unsigned long) > sizeof(unsigned int))
2515 for (i
= 0; i
< SIGSZ
; ++i
)
2517 unsigned long s
= r
->sig
[i
];
2518 h
^= s
^ (s
>> (HOST_BITS_PER_LONG
/ 2));
2521 for (i
= 0; i
< SIGSZ
; ++i
)
2527 /* IEEE single-precision format. */
2529 static void encode_ieee_single (const struct real_format
*fmt
,
2530 long *, const REAL_VALUE_TYPE
*);
2531 static void decode_ieee_single (const struct real_format
*,
2532 REAL_VALUE_TYPE
*, const long *);
2535 encode_ieee_single (const struct real_format
*fmt
, long *buf
,
2536 const REAL_VALUE_TYPE
*r
)
2538 unsigned long image
, sig
, exp
;
2539 unsigned long sign
= r
->sign
;
2540 bool denormal
= (r
->sig
[SIGSZ
-1] & SIG_MSB
) == 0;
2543 sig
= (r
->sig
[SIGSZ
-1] >> (HOST_BITS_PER_LONG
- 24)) & 0x7fffff;
2554 image
|= 0x7fffffff;
2562 if (r
->signalling
== fmt
->qnan_msb_set
)
2566 /* We overload qnan_msb_set here: it's only clear for
2567 mips_ieee_single, which wants all mantissa bits but the
2568 quiet/signalling one set in canonical NaNs (at least
2570 if (r
->canonical
&& !fmt
->qnan_msb_set
)
2571 sig
|= (1 << 22) - 1;
2579 image
|= 0x7fffffff;
2583 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2584 whereas the intermediate representation is 0.F x 2**exp.
2585 Which means we're off by one. */
2589 exp
= r
->exp
+ 127 - 1;
2602 decode_ieee_single (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
,
2605 unsigned long image
= buf
[0] & 0xffffffff;
2606 bool sign
= (image
>> 31) & 1;
2607 int exp
= (image
>> 23) & 0xff;
2609 memset (r
, 0, sizeof (*r
));
2610 image
<<= HOST_BITS_PER_LONG
- 24;
2615 if (image
&& fmt
->has_denorm
)
2617 r
->class = rvc_normal
;
2620 r
->sig
[SIGSZ
-1] = image
<< 1;
2623 else if (fmt
->has_signed_zero
)
2626 else if (exp
== 255 && (fmt
->has_nans
|| fmt
->has_inf
))
2632 r
->signalling
= (((image
>> (HOST_BITS_PER_LONG
- 2)) & 1)
2633 ^ fmt
->qnan_msb_set
);
2634 r
->sig
[SIGSZ
-1] = image
;
2644 r
->class = rvc_normal
;
2646 r
->exp
= exp
- 127 + 1;
2647 r
->sig
[SIGSZ
-1] = image
| SIG_MSB
;
2651 const struct real_format ieee_single_format
=
2669 const struct real_format mips_single_format
=
2688 /* IEEE double-precision format. */
2690 static void encode_ieee_double (const struct real_format
*fmt
,
2691 long *, const REAL_VALUE_TYPE
*);
2692 static void decode_ieee_double (const struct real_format
*,
2693 REAL_VALUE_TYPE
*, const long *);
2696 encode_ieee_double (const struct real_format
*fmt
, long *buf
,
2697 const REAL_VALUE_TYPE
*r
)
2699 unsigned long image_lo
, image_hi
, sig_lo
, sig_hi
, exp
;
2700 bool denormal
= (r
->sig
[SIGSZ
-1] & SIG_MSB
) == 0;
2702 image_hi
= r
->sign
<< 31;
2705 if (HOST_BITS_PER_LONG
== 64)
2707 sig_hi
= r
->sig
[SIGSZ
-1];
2708 sig_lo
= (sig_hi
>> (64 - 53)) & 0xffffffff;
2709 sig_hi
= (sig_hi
>> (64 - 53 + 1) >> 31) & 0xfffff;
2713 sig_hi
= r
->sig
[SIGSZ
-1];
2714 sig_lo
= r
->sig
[SIGSZ
-2];
2715 sig_lo
= (sig_hi
<< 21) | (sig_lo
>> 11);
2716 sig_hi
= (sig_hi
>> 11) & 0xfffff;
2726 image_hi
|= 2047 << 20;
2729 image_hi
|= 0x7fffffff;
2730 image_lo
= 0xffffffff;
2738 sig_hi
= sig_lo
= 0;
2739 if (r
->signalling
== fmt
->qnan_msb_set
)
2740 sig_hi
&= ~(1 << 19);
2743 /* We overload qnan_msb_set here: it's only clear for
2744 mips_ieee_single, which wants all mantissa bits but the
2745 quiet/signalling one set in canonical NaNs (at least
2747 if (r
->canonical
&& !fmt
->qnan_msb_set
)
2749 sig_hi
|= (1 << 19) - 1;
2750 sig_lo
= 0xffffffff;
2752 else if (sig_hi
== 0 && sig_lo
== 0)
2755 image_hi
|= 2047 << 20;
2761 image_hi
|= 0x7fffffff;
2762 image_lo
= 0xffffffff;
2767 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2768 whereas the intermediate representation is 0.F x 2**exp.
2769 Which means we're off by one. */
2773 exp
= r
->exp
+ 1023 - 1;
2774 image_hi
|= exp
<< 20;
2783 if (FLOAT_WORDS_BIG_ENDIAN
)
2784 buf
[0] = image_hi
, buf
[1] = image_lo
;
2786 buf
[0] = image_lo
, buf
[1] = image_hi
;
2790 decode_ieee_double (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
,
2793 unsigned long image_hi
, image_lo
;
2797 if (FLOAT_WORDS_BIG_ENDIAN
)
2798 image_hi
= buf
[0], image_lo
= buf
[1];
2800 image_lo
= buf
[0], image_hi
= buf
[1];
2801 image_lo
&= 0xffffffff;
2802 image_hi
&= 0xffffffff;
2804 sign
= (image_hi
>> 31) & 1;
2805 exp
= (image_hi
>> 20) & 0x7ff;
2807 memset (r
, 0, sizeof (*r
));
2809 image_hi
<<= 32 - 21;
2810 image_hi
|= image_lo
>> 21;
2811 image_hi
&= 0x7fffffff;
2812 image_lo
<<= 32 - 21;
2816 if ((image_hi
|| image_lo
) && fmt
->has_denorm
)
2818 r
->class = rvc_normal
;
2821 if (HOST_BITS_PER_LONG
== 32)
2823 image_hi
= (image_hi
<< 1) | (image_lo
>> 31);
2825 r
->sig
[SIGSZ
-1] = image_hi
;
2826 r
->sig
[SIGSZ
-2] = image_lo
;
2830 image_hi
= (image_hi
<< 31 << 2) | (image_lo
<< 1);
2831 r
->sig
[SIGSZ
-1] = image_hi
;
2835 else if (fmt
->has_signed_zero
)
2838 else if (exp
== 2047 && (fmt
->has_nans
|| fmt
->has_inf
))
2840 if (image_hi
|| image_lo
)
2844 r
->signalling
= ((image_hi
>> 30) & 1) ^ fmt
->qnan_msb_set
;
2845 if (HOST_BITS_PER_LONG
== 32)
2847 r
->sig
[SIGSZ
-1] = image_hi
;
2848 r
->sig
[SIGSZ
-2] = image_lo
;
2851 r
->sig
[SIGSZ
-1] = (image_hi
<< 31 << 1) | image_lo
;
2861 r
->class = rvc_normal
;
2863 r
->exp
= exp
- 1023 + 1;
2864 if (HOST_BITS_PER_LONG
== 32)
2866 r
->sig
[SIGSZ
-1] = image_hi
| SIG_MSB
;
2867 r
->sig
[SIGSZ
-2] = image_lo
;
2870 r
->sig
[SIGSZ
-1] = (image_hi
<< 31 << 1) | image_lo
| SIG_MSB
;
2874 const struct real_format ieee_double_format
=
2892 const struct real_format mips_double_format
=
2911 /* IEEE extended double precision format. This comes in three
2912 flavors: Intel's as a 12 byte image, Intel's as a 16 byte image,
2915 static void encode_ieee_extended (const struct real_format
*fmt
,
2916 long *, const REAL_VALUE_TYPE
*);
2917 static void decode_ieee_extended (const struct real_format
*,
2918 REAL_VALUE_TYPE
*, const long *);
2920 static void encode_ieee_extended_128 (const struct real_format
*fmt
,
2921 long *, const REAL_VALUE_TYPE
*);
2922 static void decode_ieee_extended_128 (const struct real_format
*,
2923 REAL_VALUE_TYPE
*, const long *);
2926 encode_ieee_extended (const struct real_format
*fmt
, long *buf
,
2927 const REAL_VALUE_TYPE
*r
)
2929 unsigned long image_hi
, sig_hi
, sig_lo
;
2930 bool denormal
= (r
->sig
[SIGSZ
-1] & SIG_MSB
) == 0;
2932 image_hi
= r
->sign
<< 15;
2933 sig_hi
= sig_lo
= 0;
2945 /* Intel requires the explicit integer bit to be set, otherwise
2946 it considers the value a "pseudo-infinity". Motorola docs
2947 say it doesn't care. */
2948 sig_hi
= 0x80000000;
2953 sig_lo
= sig_hi
= 0xffffffff;
2961 if (HOST_BITS_PER_LONG
== 32)
2963 sig_hi
= r
->sig
[SIGSZ
-1];
2964 sig_lo
= r
->sig
[SIGSZ
-2];
2968 sig_lo
= r
->sig
[SIGSZ
-1];
2969 sig_hi
= sig_lo
>> 31 >> 1;
2970 sig_lo
&= 0xffffffff;
2972 if (r
->signalling
== fmt
->qnan_msb_set
)
2973 sig_hi
&= ~(1 << 30);
2976 if ((sig_hi
& 0x7fffffff) == 0 && sig_lo
== 0)
2979 /* Intel requires the explicit integer bit to be set, otherwise
2980 it considers the value a "pseudo-nan". Motorola docs say it
2982 sig_hi
|= 0x80000000;
2987 sig_lo
= sig_hi
= 0xffffffff;
2995 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2996 whereas the intermediate representation is 0.F x 2**exp.
2997 Which means we're off by one.
2999 Except for Motorola, which consider exp=0 and explicit
3000 integer bit set to continue to be normalized. In theory
3001 this discrepancy has been taken care of by the difference
3002 in fmt->emin in round_for_format. */
3014 if (HOST_BITS_PER_LONG
== 32)
3016 sig_hi
= r
->sig
[SIGSZ
-1];
3017 sig_lo
= r
->sig
[SIGSZ
-2];
3021 sig_lo
= r
->sig
[SIGSZ
-1];
3022 sig_hi
= sig_lo
>> 31 >> 1;
3023 sig_lo
&= 0xffffffff;
3032 if (FLOAT_WORDS_BIG_ENDIAN
)
3033 buf
[0] = image_hi
<< 16, buf
[1] = sig_hi
, buf
[2] = sig_lo
;
3035 buf
[0] = sig_lo
, buf
[1] = sig_hi
, buf
[2] = image_hi
;
3039 encode_ieee_extended_128 (const struct real_format
*fmt
, long *buf
,
3040 const REAL_VALUE_TYPE
*r
)
3042 buf
[3 * !FLOAT_WORDS_BIG_ENDIAN
] = 0;
3043 encode_ieee_extended (fmt
, buf
+!!FLOAT_WORDS_BIG_ENDIAN
, r
);
3047 decode_ieee_extended (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
,
3050 unsigned long image_hi
, sig_hi
, sig_lo
;
3054 if (FLOAT_WORDS_BIG_ENDIAN
)
3055 image_hi
= buf
[0] >> 16, sig_hi
= buf
[1], sig_lo
= buf
[2];
3057 sig_lo
= buf
[0], sig_hi
= buf
[1], image_hi
= buf
[2];
3058 sig_lo
&= 0xffffffff;
3059 sig_hi
&= 0xffffffff;
3060 image_hi
&= 0xffffffff;
3062 sign
= (image_hi
>> 15) & 1;
3063 exp
= image_hi
& 0x7fff;
3065 memset (r
, 0, sizeof (*r
));
3069 if ((sig_hi
|| sig_lo
) && fmt
->has_denorm
)
3071 r
->class = rvc_normal
;
3074 /* When the IEEE format contains a hidden bit, we know that
3075 it's zero at this point, and so shift up the significand
3076 and decrease the exponent to match. In this case, Motorola
3077 defines the explicit integer bit to be valid, so we don't
3078 know whether the msb is set or not. */
3080 if (HOST_BITS_PER_LONG
== 32)
3082 r
->sig
[SIGSZ
-1] = sig_hi
;
3083 r
->sig
[SIGSZ
-2] = sig_lo
;
3086 r
->sig
[SIGSZ
-1] = (sig_hi
<< 31 << 1) | sig_lo
;
3090 else if (fmt
->has_signed_zero
)
3093 else if (exp
== 32767 && (fmt
->has_nans
|| fmt
->has_inf
))
3095 /* See above re "pseudo-infinities" and "pseudo-nans".
3096 Short summary is that the MSB will likely always be
3097 set, and that we don't care about it. */
3098 sig_hi
&= 0x7fffffff;
3100 if (sig_hi
|| sig_lo
)
3104 r
->signalling
= ((sig_hi
>> 30) & 1) ^ fmt
->qnan_msb_set
;
3105 if (HOST_BITS_PER_LONG
== 32)
3107 r
->sig
[SIGSZ
-1] = sig_hi
;
3108 r
->sig
[SIGSZ
-2] = sig_lo
;
3111 r
->sig
[SIGSZ
-1] = (sig_hi
<< 31 << 1) | sig_lo
;
3121 r
->class = rvc_normal
;
3123 r
->exp
= exp
- 16383 + 1;
3124 if (HOST_BITS_PER_LONG
== 32)
3126 r
->sig
[SIGSZ
-1] = sig_hi
;
3127 r
->sig
[SIGSZ
-2] = sig_lo
;
3130 r
->sig
[SIGSZ
-1] = (sig_hi
<< 31 << 1) | sig_lo
;
3135 decode_ieee_extended_128 (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
,
3138 decode_ieee_extended (fmt
, r
, buf
+!!FLOAT_WORDS_BIG_ENDIAN
);
3141 const struct real_format ieee_extended_motorola_format
=
3143 encode_ieee_extended
,
3144 decode_ieee_extended
,
3159 const struct real_format ieee_extended_intel_96_format
=
3161 encode_ieee_extended
,
3162 decode_ieee_extended
,
3177 const struct real_format ieee_extended_intel_128_format
=
3179 encode_ieee_extended_128
,
3180 decode_ieee_extended_128
,
3195 /* The following caters to i386 systems that set the rounding precision
3196 to 53 bits instead of 64, e.g. FreeBSD. */
3197 const struct real_format ieee_extended_intel_96_round_53_format
=
3199 encode_ieee_extended
,
3200 decode_ieee_extended
,
3215 /* IBM 128-bit extended precision format: a pair of IEEE double precision
3216 numbers whose sum is equal to the extended precision value. The number
3217 with greater magnitude is first. This format has the same magnitude
3218 range as an IEEE double precision value, but effectively 106 bits of
3219 significand precision. Infinity and NaN are represented by their IEEE
3220 double precision value stored in the first number, the second number is
3221 ignored. Zeroes, Infinities, and NaNs are set in both doubles
3222 due to precedent. */
3224 static void encode_ibm_extended (const struct real_format
*fmt
,
3225 long *, const REAL_VALUE_TYPE
*);
3226 static void decode_ibm_extended (const struct real_format
*,
3227 REAL_VALUE_TYPE
*, const long *);
3230 encode_ibm_extended (const struct real_format
*fmt
, long *buf
,
3231 const REAL_VALUE_TYPE
*r
)
3233 REAL_VALUE_TYPE u
, normr
, v
;
3234 const struct real_format
*base_fmt
;
3236 base_fmt
= fmt
->qnan_msb_set
? &ieee_double_format
: &mips_double_format
;
3238 /* Renormlize R before doing any arithmetic on it. */
3240 if (normr
.class == rvc_normal
)
3243 /* u = IEEE double precision portion of significand. */
3245 round_for_format (base_fmt
, &u
);
3246 encode_ieee_double (base_fmt
, &buf
[0], &u
);
3248 if (u
.class == rvc_normal
)
3250 do_add (&v
, &normr
, &u
, 1);
3251 /* The low double won't need rounding, since we round to a 106 bit
3252 mantissa before calling this function, and we've just
3253 subtracted off the top 54 bits. (53+1 because u is rounded.) */
3254 encode_ieee_double (base_fmt
, &buf
[2], &v
);
3258 /* Inf, NaN, 0 are all representable as doubles, so the
3259 least-significant part can be 0.0. */
3266 decode_ibm_extended (const struct real_format
*fmt ATTRIBUTE_UNUSED
, REAL_VALUE_TYPE
*r
,
3269 REAL_VALUE_TYPE u
, v
;
3270 const struct real_format
*base_fmt
;
3272 base_fmt
= fmt
->qnan_msb_set
? &ieee_double_format
: &mips_double_format
;
3273 decode_ieee_double (base_fmt
, &u
, &buf
[0]);
3275 if (u
.class != rvc_zero
&& u
.class != rvc_inf
&& u
.class != rvc_nan
)
3277 decode_ieee_double (base_fmt
, &v
, &buf
[2]);
3278 do_add (r
, &u
, &v
, 0);
3284 const struct real_format ibm_extended_format
=
3286 encode_ibm_extended
,
3287 decode_ibm_extended
,
3302 const struct real_format mips_extended_format
=
3304 encode_ibm_extended
,
3305 decode_ibm_extended
,
3321 /* IEEE quad precision format. */
3323 static void encode_ieee_quad (const struct real_format
*fmt
,
3324 long *, const REAL_VALUE_TYPE
*);
3325 static void decode_ieee_quad (const struct real_format
*,
3326 REAL_VALUE_TYPE
*, const long *);
3329 encode_ieee_quad (const struct real_format
*fmt
, long *buf
,
3330 const REAL_VALUE_TYPE
*r
)
3332 unsigned long image3
, image2
, image1
, image0
, exp
;
3333 bool denormal
= (r
->sig
[SIGSZ
-1] & SIG_MSB
) == 0;
3336 image3
= r
->sign
<< 31;
3341 rshift_significand (&u
, r
, SIGNIFICAND_BITS
- 113);
3350 image3
|= 32767 << 16;
3353 image3
|= 0x7fffffff;
3354 image2
= 0xffffffff;
3355 image1
= 0xffffffff;
3356 image0
= 0xffffffff;
3363 image3
|= 32767 << 16;
3367 /* Don't use bits from the significand. The
3368 initialization above is right. */
3370 else if (HOST_BITS_PER_LONG
== 32)
3375 image3
|= u
.sig
[3] & 0xffff;
3380 image1
= image0
>> 31 >> 1;
3382 image3
|= (image2
>> 31 >> 1) & 0xffff;
3383 image0
&= 0xffffffff;
3384 image2
&= 0xffffffff;
3386 if (r
->signalling
== fmt
->qnan_msb_set
)
3390 /* We overload qnan_msb_set here: it's only clear for
3391 mips_ieee_single, which wants all mantissa bits but the
3392 quiet/signalling one set in canonical NaNs (at least
3394 if (r
->canonical
&& !fmt
->qnan_msb_set
)
3397 image2
= image1
= image0
= 0xffffffff;
3399 else if (((image3
& 0xffff) | image2
| image1
| image0
) == 0)
3404 image3
|= 0x7fffffff;
3405 image2
= 0xffffffff;
3406 image1
= 0xffffffff;
3407 image0
= 0xffffffff;
3412 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3413 whereas the intermediate representation is 0.F x 2**exp.
3414 Which means we're off by one. */
3418 exp
= r
->exp
+ 16383 - 1;
3419 image3
|= exp
<< 16;
3421 if (HOST_BITS_PER_LONG
== 32)
3426 image3
|= u
.sig
[3] & 0xffff;
3431 image1
= image0
>> 31 >> 1;
3433 image3
|= (image2
>> 31 >> 1) & 0xffff;
3434 image0
&= 0xffffffff;
3435 image2
&= 0xffffffff;
3443 if (FLOAT_WORDS_BIG_ENDIAN
)
3460 decode_ieee_quad (const struct real_format
*fmt
, REAL_VALUE_TYPE
*r
,
3463 unsigned long image3
, image2
, image1
, image0
;
3467 if (FLOAT_WORDS_BIG_ENDIAN
)
3481 image0
&= 0xffffffff;
3482 image1
&= 0xffffffff;
3483 image2
&= 0xffffffff;
3485 sign
= (image3
>> 31) & 1;
3486 exp
= (image3
>> 16) & 0x7fff;
3489 memset (r
, 0, sizeof (*r
));
3493 if ((image3
| image2
| image1
| image0
) && fmt
->has_denorm
)
3495 r
->class = rvc_normal
;
3498 r
->exp
= -16382 + (SIGNIFICAND_BITS
- 112);
3499 if (HOST_BITS_PER_LONG
== 32)
3508 r
->sig
[0] = (image1
<< 31 << 1) | image0
;
3509 r
->sig
[1] = (image3
<< 31 << 1) | image2
;
3514 else if (fmt
->has_signed_zero
)
3517 else if (exp
== 32767 && (fmt
->has_nans
|| fmt
->has_inf
))
3519 if (image3
| image2
| image1
| image0
)
3523 r
->signalling
= ((image3
>> 15) & 1) ^ fmt
->qnan_msb_set
;
3525 if (HOST_BITS_PER_LONG
== 32)
3534 r
->sig
[0] = (image1
<< 31 << 1) | image0
;
3535 r
->sig
[1] = (image3
<< 31 << 1) | image2
;
3537 lshift_significand (r
, r
, SIGNIFICAND_BITS
- 113);
3547 r
->class = rvc_normal
;
3549 r
->exp
= exp
- 16383 + 1;
3551 if (HOST_BITS_PER_LONG
== 32)
3560 r
->sig
[0] = (image1
<< 31 << 1) | image0
;
3561 r
->sig
[1] = (image3
<< 31 << 1) | image2
;
3563 lshift_significand (r
, r
, SIGNIFICAND_BITS
- 113);
3564 r
->sig
[SIGSZ
-1] |= SIG_MSB
;
3568 const struct real_format ieee_quad_format
=
3586 const struct real_format mips_quad_format
=
3604 /* Descriptions of VAX floating point formats can be found beginning at
3606 http://h71000.www7.hp.com/doc/73FINAL/4515/4515pro_013.html#f_floating_point_format
3608 The thing to remember is that they're almost IEEE, except for word
3609 order, exponent bias, and the lack of infinities, nans, and denormals.
3611 We don't implement the H_floating format here, simply because neither
3612 the VAX or Alpha ports use it. */
3614 static void encode_vax_f (const struct real_format
*fmt
,
3615 long *, const REAL_VALUE_TYPE
*);
3616 static void decode_vax_f (const struct real_format
*,
3617 REAL_VALUE_TYPE
*, const long *);
3618 static void encode_vax_d (const struct real_format
*fmt
,
3619 long *, const REAL_VALUE_TYPE
*);
3620 static void decode_vax_d (const struct real_format
*,
3621 REAL_VALUE_TYPE
*, const long *);
3622 static void encode_vax_g (const struct real_format
*fmt
,
3623 long *, const REAL_VALUE_TYPE
*);
3624 static void decode_vax_g (const struct real_format
*,
3625 REAL_VALUE_TYPE
*, const long *);
3628 encode_vax_f (const struct real_format
*fmt ATTRIBUTE_UNUSED
, long *buf
,
3629 const REAL_VALUE_TYPE
*r
)
3631 unsigned long sign
, exp
, sig
, image
;
3633 sign
= r
->sign
<< 15;
3643 image
= 0xffff7fff | sign
;
3647 sig
= (r
->sig
[SIGSZ
-1] >> (HOST_BITS_PER_LONG
- 24)) & 0x7fffff;
3650 image
= (sig
<< 16) & 0xffff0000;
3664 decode_vax_f (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
3665 REAL_VALUE_TYPE
*r
, const long *buf
)
3667 unsigned long image
= buf
[0] & 0xffffffff;
3668 int exp
= (image
>> 7) & 0xff;
3670 memset (r
, 0, sizeof (*r
));
3674 r
->class = rvc_normal
;
3675 r
->sign
= (image
>> 15) & 1;
3678 image
= ((image
& 0x7f) << 16) | ((image
>> 16) & 0xffff);
3679 r
->sig
[SIGSZ
-1] = (image
<< (HOST_BITS_PER_LONG
- 24)) | SIG_MSB
;
3684 encode_vax_d (const struct real_format
*fmt ATTRIBUTE_UNUSED
, long *buf
,
3685 const REAL_VALUE_TYPE
*r
)
3687 unsigned long image0
, image1
, sign
= r
->sign
<< 15;
3692 image0
= image1
= 0;
3697 image0
= 0xffff7fff | sign
;
3698 image1
= 0xffffffff;
3702 /* Extract the significand into straight hi:lo. */
3703 if (HOST_BITS_PER_LONG
== 64)
3705 image0
= r
->sig
[SIGSZ
-1];
3706 image1
= (image0
>> (64 - 56)) & 0xffffffff;
3707 image0
= (image0
>> (64 - 56 + 1) >> 31) & 0x7fffff;
3711 image0
= r
->sig
[SIGSZ
-1];
3712 image1
= r
->sig
[SIGSZ
-2];
3713 image1
= (image0
<< 24) | (image1
>> 8);
3714 image0
= (image0
>> 8) & 0xffffff;
3717 /* Rearrange the half-words of the significand to match the
3719 image0
= ((image0
<< 16) | (image0
>> 16)) & 0xffff007f;
3720 image1
= ((image1
<< 16) | (image1
>> 16)) & 0xffffffff;
3722 /* Add the sign and exponent. */
3724 image0
|= (r
->exp
+ 128) << 7;
3731 if (FLOAT_WORDS_BIG_ENDIAN
)
3732 buf
[0] = image1
, buf
[1] = image0
;
3734 buf
[0] = image0
, buf
[1] = image1
;
3738 decode_vax_d (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
3739 REAL_VALUE_TYPE
*r
, const long *buf
)
3741 unsigned long image0
, image1
;
3744 if (FLOAT_WORDS_BIG_ENDIAN
)
3745 image1
= buf
[0], image0
= buf
[1];
3747 image0
= buf
[0], image1
= buf
[1];
3748 image0
&= 0xffffffff;
3749 image1
&= 0xffffffff;
3751 exp
= (image0
>> 7) & 0xff;
3753 memset (r
, 0, sizeof (*r
));
3757 r
->class = rvc_normal
;
3758 r
->sign
= (image0
>> 15) & 1;
3761 /* Rearrange the half-words of the external format into
3762 proper ascending order. */
3763 image0
= ((image0
& 0x7f) << 16) | ((image0
>> 16) & 0xffff);
3764 image1
= ((image1
& 0xffff) << 16) | ((image1
>> 16) & 0xffff);
3766 if (HOST_BITS_PER_LONG
== 64)
3768 image0
= (image0
<< 31 << 1) | image1
;
3771 r
->sig
[SIGSZ
-1] = image0
;
3775 r
->sig
[SIGSZ
-1] = image0
;
3776 r
->sig
[SIGSZ
-2] = image1
;
3777 lshift_significand (r
, r
, 2*HOST_BITS_PER_LONG
- 56);
3778 r
->sig
[SIGSZ
-1] |= SIG_MSB
;
3784 encode_vax_g (const struct real_format
*fmt ATTRIBUTE_UNUSED
, long *buf
,
3785 const REAL_VALUE_TYPE
*r
)
3787 unsigned long image0
, image1
, sign
= r
->sign
<< 15;
3792 image0
= image1
= 0;
3797 image0
= 0xffff7fff | sign
;
3798 image1
= 0xffffffff;
3802 /* Extract the significand into straight hi:lo. */
3803 if (HOST_BITS_PER_LONG
== 64)
3805 image0
= r
->sig
[SIGSZ
-1];
3806 image1
= (image0
>> (64 - 53)) & 0xffffffff;
3807 image0
= (image0
>> (64 - 53 + 1) >> 31) & 0xfffff;
3811 image0
= r
->sig
[SIGSZ
-1];
3812 image1
= r
->sig
[SIGSZ
-2];
3813 image1
= (image0
<< 21) | (image1
>> 11);
3814 image0
= (image0
>> 11) & 0xfffff;
3817 /* Rearrange the half-words of the significand to match the
3819 image0
= ((image0
<< 16) | (image0
>> 16)) & 0xffff000f;
3820 image1
= ((image1
<< 16) | (image1
>> 16)) & 0xffffffff;
3822 /* Add the sign and exponent. */
3824 image0
|= (r
->exp
+ 1024) << 4;
3831 if (FLOAT_WORDS_BIG_ENDIAN
)
3832 buf
[0] = image1
, buf
[1] = image0
;
3834 buf
[0] = image0
, buf
[1] = image1
;
3838 decode_vax_g (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
3839 REAL_VALUE_TYPE
*r
, const long *buf
)
3841 unsigned long image0
, image1
;
3844 if (FLOAT_WORDS_BIG_ENDIAN
)
3845 image1
= buf
[0], image0
= buf
[1];
3847 image0
= buf
[0], image1
= buf
[1];
3848 image0
&= 0xffffffff;
3849 image1
&= 0xffffffff;
3851 exp
= (image0
>> 4) & 0x7ff;
3853 memset (r
, 0, sizeof (*r
));
3857 r
->class = rvc_normal
;
3858 r
->sign
= (image0
>> 15) & 1;
3859 r
->exp
= exp
- 1024;
3861 /* Rearrange the half-words of the external format into
3862 proper ascending order. */
3863 image0
= ((image0
& 0xf) << 16) | ((image0
>> 16) & 0xffff);
3864 image1
= ((image1
& 0xffff) << 16) | ((image1
>> 16) & 0xffff);
3866 if (HOST_BITS_PER_LONG
== 64)
3868 image0
= (image0
<< 31 << 1) | image1
;
3871 r
->sig
[SIGSZ
-1] = image0
;
3875 r
->sig
[SIGSZ
-1] = image0
;
3876 r
->sig
[SIGSZ
-2] = image1
;
3877 lshift_significand (r
, r
, 64 - 53);
3878 r
->sig
[SIGSZ
-1] |= SIG_MSB
;
3883 const struct real_format vax_f_format
=
3901 const struct real_format vax_d_format
=
3919 const struct real_format vax_g_format
=
3937 /* A good reference for these can be found in chapter 9 of
3938 "ESA/390 Principles of Operation", IBM document number SA22-7201-01.
3939 An on-line version can be found here:
3941 http://publibz.boulder.ibm.com/cgi-bin/bookmgr_OS390/BOOKS/DZ9AR001/9.1?DT=19930923083613
3944 static void encode_i370_single (const struct real_format
*fmt
,
3945 long *, const REAL_VALUE_TYPE
*);
3946 static void decode_i370_single (const struct real_format
*,
3947 REAL_VALUE_TYPE
*, const long *);
3948 static void encode_i370_double (const struct real_format
*fmt
,
3949 long *, const REAL_VALUE_TYPE
*);
3950 static void decode_i370_double (const struct real_format
*,
3951 REAL_VALUE_TYPE
*, const long *);
3954 encode_i370_single (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
3955 long *buf
, const REAL_VALUE_TYPE
*r
)
3957 unsigned long sign
, exp
, sig
, image
;
3959 sign
= r
->sign
<< 31;
3969 image
= 0x7fffffff | sign
;
3973 sig
= (r
->sig
[SIGSZ
-1] >> (HOST_BITS_PER_LONG
- 24)) & 0xffffff;
3974 exp
= ((r
->exp
/ 4) + 64) << 24;
3975 image
= sign
| exp
| sig
;
3986 decode_i370_single (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
3987 REAL_VALUE_TYPE
*r
, const long *buf
)
3989 unsigned long sign
, sig
, image
= buf
[0];
3992 sign
= (image
>> 31) & 1;
3993 exp
= (image
>> 24) & 0x7f;
3994 sig
= image
& 0xffffff;
3996 memset (r
, 0, sizeof (*r
));
4000 r
->class = rvc_normal
;
4002 r
->exp
= (exp
- 64) * 4;
4003 r
->sig
[SIGSZ
-1] = sig
<< (HOST_BITS_PER_LONG
- 24);
4009 encode_i370_double (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4010 long *buf
, const REAL_VALUE_TYPE
*r
)
4012 unsigned long sign
, exp
, image_hi
, image_lo
;
4014 sign
= r
->sign
<< 31;
4019 image_hi
= image_lo
= 0;
4024 image_hi
= 0x7fffffff | sign
;
4025 image_lo
= 0xffffffff;
4029 if (HOST_BITS_PER_LONG
== 64)
4031 image_hi
= r
->sig
[SIGSZ
-1];
4032 image_lo
= (image_hi
>> (64 - 56)) & 0xffffffff;
4033 image_hi
= (image_hi
>> (64 - 56 + 1) >> 31) & 0xffffff;
4037 image_hi
= r
->sig
[SIGSZ
-1];
4038 image_lo
= r
->sig
[SIGSZ
-2];
4039 image_lo
= (image_lo
>> 8) | (image_hi
<< 24);
4043 exp
= ((r
->exp
/ 4) + 64) << 24;
4044 image_hi
|= sign
| exp
;
4051 if (FLOAT_WORDS_BIG_ENDIAN
)
4052 buf
[0] = image_hi
, buf
[1] = image_lo
;
4054 buf
[0] = image_lo
, buf
[1] = image_hi
;
4058 decode_i370_double (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4059 REAL_VALUE_TYPE
*r
, const long *buf
)
4061 unsigned long sign
, image_hi
, image_lo
;
4064 if (FLOAT_WORDS_BIG_ENDIAN
)
4065 image_hi
= buf
[0], image_lo
= buf
[1];
4067 image_lo
= buf
[0], image_hi
= buf
[1];
4069 sign
= (image_hi
>> 31) & 1;
4070 exp
= (image_hi
>> 24) & 0x7f;
4071 image_hi
&= 0xffffff;
4072 image_lo
&= 0xffffffff;
4074 memset (r
, 0, sizeof (*r
));
4076 if (exp
|| image_hi
|| image_lo
)
4078 r
->class = rvc_normal
;
4080 r
->exp
= (exp
- 64) * 4 + (SIGNIFICAND_BITS
- 56);
4082 if (HOST_BITS_PER_LONG
== 32)
4084 r
->sig
[0] = image_lo
;
4085 r
->sig
[1] = image_hi
;
4088 r
->sig
[0] = image_lo
| (image_hi
<< 31 << 1);
4094 const struct real_format i370_single_format
=
4107 false, /* ??? The encoding does allow for "unnormals". */
4108 false, /* ??? The encoding does allow for "unnormals". */
4112 const struct real_format i370_double_format
=
4125 false, /* ??? The encoding does allow for "unnormals". */
4126 false, /* ??? The encoding does allow for "unnormals". */
4130 /* The "twos-complement" c4x format is officially defined as
4134 This is rather misleading. One must remember that F is signed.
4135 A better description would be
4137 x = -1**s * ((s + 1 + .f) * 2**e
4139 So if we have a (4 bit) fraction of .1000 with a sign bit of 1,
4140 that's -1 * (1+1+(-.5)) == -1.5. I think.
4142 The constructions here are taken from Tables 5-1 and 5-2 of the
4143 TMS320C4x User's Guide wherein step-by-step instructions for
4144 conversion from IEEE are presented. That's close enough to our
4145 internal representation so as to make things easy.
4147 See http://www-s.ti.com/sc/psheets/spru063c/spru063c.pdf */
4149 static void encode_c4x_single (const struct real_format
*fmt
,
4150 long *, const REAL_VALUE_TYPE
*);
4151 static void decode_c4x_single (const struct real_format
*,
4152 REAL_VALUE_TYPE
*, const long *);
4153 static void encode_c4x_extended (const struct real_format
*fmt
,
4154 long *, const REAL_VALUE_TYPE
*);
4155 static void decode_c4x_extended (const struct real_format
*,
4156 REAL_VALUE_TYPE
*, const long *);
4159 encode_c4x_single (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4160 long *buf
, const REAL_VALUE_TYPE
*r
)
4162 unsigned long image
, exp
, sig
;
4174 sig
= 0x800000 - r
->sign
;
4179 sig
= (r
->sig
[SIGSZ
-1] >> (HOST_BITS_PER_LONG
- 24)) & 0x7fffff;
4194 image
= ((exp
& 0xff) << 24) | (sig
& 0xffffff);
4199 decode_c4x_single (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4200 REAL_VALUE_TYPE
*r
, const long *buf
)
4202 unsigned long image
= buf
[0];
4206 exp
= (((image
>> 24) & 0xff) ^ 0x80) - 0x80;
4207 sf
= ((image
& 0xffffff) ^ 0x800000) - 0x800000;
4209 memset (r
, 0, sizeof (*r
));
4213 r
->class = rvc_normal
;
4215 sig
= sf
& 0x7fffff;
4224 sig
= (sig
<< (HOST_BITS_PER_LONG
- 24)) | SIG_MSB
;
4227 r
->sig
[SIGSZ
-1] = sig
;
4232 encode_c4x_extended (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4233 long *buf
, const REAL_VALUE_TYPE
*r
)
4235 unsigned long exp
, sig
;
4247 sig
= 0x80000000 - r
->sign
;
4253 sig
= r
->sig
[SIGSZ
-1];
4254 if (HOST_BITS_PER_LONG
== 64)
4255 sig
= sig
>> 1 >> 31;
4272 exp
= (exp
& 0xff) << 24;
4275 if (FLOAT_WORDS_BIG_ENDIAN
)
4276 buf
[0] = exp
, buf
[1] = sig
;
4278 buf
[0] = sig
, buf
[0] = exp
;
4282 decode_c4x_extended (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4283 REAL_VALUE_TYPE
*r
, const long *buf
)
4288 if (FLOAT_WORDS_BIG_ENDIAN
)
4289 exp
= buf
[0], sf
= buf
[1];
4291 sf
= buf
[0], exp
= buf
[1];
4293 exp
= (((exp
>> 24) & 0xff) & 0x80) - 0x80;
4294 sf
= ((sf
& 0xffffffff) ^ 0x80000000) - 0x80000000;
4296 memset (r
, 0, sizeof (*r
));
4300 r
->class = rvc_normal
;
4302 sig
= sf
& 0x7fffffff;
4311 if (HOST_BITS_PER_LONG
== 64)
4312 sig
= sig
<< 1 << 31;
4316 r
->sig
[SIGSZ
-1] = sig
;
4320 const struct real_format c4x_single_format
=
4338 const struct real_format c4x_extended_format
=
4340 encode_c4x_extended
,
4341 decode_c4x_extended
,
4357 /* A synthetic "format" for internal arithmetic. It's the size of the
4358 internal significand minus the two bits needed for proper rounding.
4359 The encode and decode routines exist only to satisfy our paranoia
4362 static void encode_internal (const struct real_format
*fmt
,
4363 long *, const REAL_VALUE_TYPE
*);
4364 static void decode_internal (const struct real_format
*,
4365 REAL_VALUE_TYPE
*, const long *);
4368 encode_internal (const struct real_format
*fmt ATTRIBUTE_UNUSED
, long *buf
,
4369 const REAL_VALUE_TYPE
*r
)
4371 memcpy (buf
, r
, sizeof (*r
));
4375 decode_internal (const struct real_format
*fmt ATTRIBUTE_UNUSED
,
4376 REAL_VALUE_TYPE
*r
, const long *buf
)
4378 memcpy (r
, buf
, sizeof (*r
));
4381 const struct real_format real_internal_format
=
4387 SIGNIFICAND_BITS
- 2,
4388 SIGNIFICAND_BITS
- 2,
4399 /* Calculate the square root of X in mode MODE, and store the result
4400 in R. Return TRUE if the operation does not raise an exception.
4401 For details see "High Precision Division and Square Root",
4402 Alan H. Karp and Peter Markstein, HP Lab Report 93-93-42, June
4403 1993. http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf. */
4406 real_sqrt (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4407 const REAL_VALUE_TYPE
*x
)
4409 static REAL_VALUE_TYPE halfthree
;
4410 static bool init
= false;
4411 REAL_VALUE_TYPE h
, t
, i
;
4414 /* sqrt(-0.0) is -0.0. */
4415 if (real_isnegzero (x
))
4421 /* Negative arguments return NaN. */
4424 get_canonical_qnan (r
, 0);
4428 /* Infinity and NaN return themselves. */
4429 if (real_isinf (x
) || real_isnan (x
))
4437 do_add (&halfthree
, &dconst1
, &dconsthalf
, 0);
4441 /* Initial guess for reciprocal sqrt, i. */
4442 exp
= real_exponent (x
);
4443 real_ldexp (&i
, &dconst1
, -exp
/2);
4445 /* Newton's iteration for reciprocal sqrt, i. */
4446 for (iter
= 0; iter
< 16; iter
++)
4448 /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x). */
4449 do_multiply (&t
, x
, &i
);
4450 do_multiply (&h
, &t
, &i
);
4451 do_multiply (&t
, &h
, &dconsthalf
);
4452 do_add (&h
, &halfthree
, &t
, 1);
4453 do_multiply (&t
, &i
, &h
);
4455 /* Check for early convergence. */
4456 if (iter
>= 6 && real_identical (&i
, &t
))
4459 /* ??? Unroll loop to avoid copying. */
4463 /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)). */
4464 do_multiply (&t
, x
, &i
);
4465 do_multiply (&h
, &t
, &i
);
4466 do_add (&i
, &dconst1
, &h
, 1);
4467 do_multiply (&h
, &t
, &i
);
4468 do_multiply (&i
, &dconsthalf
, &h
);
4469 do_add (&h
, &t
, &i
, 0);
4471 /* ??? We need a Tuckerman test to get the last bit. */
4473 real_convert (r
, mode
, &h
);
4477 /* Calculate X raised to the integer exponent N in mode MODE and store
4478 the result in R. Return true if the result may be inexact due to
4479 loss of precision. The algorithm is the classic "left-to-right binary
4480 method" described in section 4.6.3 of Donald Knuth's "Seminumerical
4481 Algorithms", "The Art of Computer Programming", Volume 2. */
4484 real_powi (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4485 const REAL_VALUE_TYPE
*x
, HOST_WIDE_INT n
)
4487 unsigned HOST_WIDE_INT bit
;
4489 bool inexact
= false;
4501 /* Don't worry about overflow, from now on n is unsigned. */
4509 bit
= (unsigned HOST_WIDE_INT
) 1 << (HOST_BITS_PER_WIDE_INT
- 1);
4510 for (i
= 0; i
< HOST_BITS_PER_WIDE_INT
; i
++)
4514 inexact
|= do_multiply (&t
, &t
, &t
);
4516 inexact
|= do_multiply (&t
, &t
, x
);
4524 inexact
|= do_divide (&t
, &dconst1
, &t
);
4526 real_convert (r
, mode
, &t
);
4530 /* Round X to the nearest integer not larger in absolute value, i.e.
4531 towards zero, placing the result in R in mode MODE. */
4534 real_trunc (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4535 const REAL_VALUE_TYPE
*x
)
4537 do_fix_trunc (r
, x
);
4538 if (mode
!= VOIDmode
)
4539 real_convert (r
, mode
, r
);
4542 /* Round X to the largest integer not greater in value, i.e. round
4543 down, placing the result in R in mode MODE. */
4546 real_floor (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4547 const REAL_VALUE_TYPE
*x
)
4551 do_fix_trunc (&t
, x
);
4552 if (! real_identical (&t
, x
) && x
->sign
)
4553 do_add (&t
, &t
, &dconstm1
, 0);
4554 if (mode
!= VOIDmode
)
4555 real_convert (r
, mode
, &t
);
4558 /* Round X to the smallest integer not less then argument, i.e. round
4559 up, placing the result in R in mode MODE. */
4562 real_ceil (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4563 const REAL_VALUE_TYPE
*x
)
4567 do_fix_trunc (&t
, x
);
4568 if (! real_identical (&t
, x
) && ! x
->sign
)
4569 do_add (&t
, &t
, &dconst1
, 0);
4570 if (mode
!= VOIDmode
)
4571 real_convert (r
, mode
, &t
);
4574 /* Round X to the nearest integer, but round halfway cases away from
4578 real_round (REAL_VALUE_TYPE
*r
, enum machine_mode mode
,
4579 const REAL_VALUE_TYPE
*x
)
4581 do_add (r
, x
, &dconsthalf
, x
->sign
);
4582 do_fix_trunc (r
, r
);
4583 if (mode
!= VOIDmode
)
4584 real_convert (r
, mode
, r
);