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 27.
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;
1084 return do_compare (op0, op1, 0) != 0;
1091 /* Return floor log2(R). */
1094 real_exponent (const REAL_VALUE_TYPE *r)
1102 return (unsigned int)-1 >> 1;
1104 return REAL_EXP (r);
1110 /* R = OP0 * 2**EXP. */
1113 real_ldexp (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *op0, int exp)
1124 exp += REAL_EXP (op0);
1126 get_inf (r, r->sign);
1127 else if (exp < -MAX_EXP)
1128 get_zero (r, r->sign);
1130 SET_REAL_EXP (r, exp);
1138 /* Determine whether a floating-point value X is infinite. */
1141 real_isinf (const REAL_VALUE_TYPE *r)
1143 return (r->class == rvc_inf);
1146 /* Determine whether a floating-point value X is a NaN. */
1149 real_isnan (const REAL_VALUE_TYPE *r)
1151 return (r->class == rvc_nan);
1154 /* Determine whether a floating-point value X is negative. */
1157 real_isneg (const REAL_VALUE_TYPE *r)
1162 /* Determine whether a floating-point value X is minus zero. */
1165 real_isnegzero (const REAL_VALUE_TYPE *r)
1167 return r->sign && r->class == rvc_zero;
1170 /* Compare two floating-point objects for bitwise identity. */
1173 real_identical (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
1177 if (a->class != b->class)
1179 if (a->sign != b->sign)
1189 if (REAL_EXP (a) != REAL_EXP (b))
1194 if (a->signalling != b->signalling)
1196 /* The significand is ignored for canonical NaNs. */
1197 if (a->canonical || b->canonical)
1198 return a->canonical == b->canonical;
1205 for (i = 0; i < SIGSZ; ++i)
1206 if (a->sig[i] != b->sig[i])
1212 /* Try to change R into its exact multiplicative inverse in machine
1213 mode MODE. Return true if successful. */
1216 exact_real_inverse (enum machine_mode mode, REAL_VALUE_TYPE *r)
1218 const REAL_VALUE_TYPE *one = real_digit (1);
1222 if (r->class != rvc_normal)
1225 /* Check for a power of two: all significand bits zero except the MSB. */
1226 for (i = 0; i < SIGSZ-1; ++i)
1229 if (r->sig[SIGSZ-1] != SIG_MSB)
1232 /* Find the inverse and truncate to the required mode. */
1233 do_divide (&u, one, r);
1234 real_convert (&u, mode, &u);
1236 /* The rounding may have overflowed. */
1237 if (u.class != rvc_normal)
1239 for (i = 0; i < SIGSZ-1; ++i)
1242 if (u.sig[SIGSZ-1] != SIG_MSB)
1249 /* Render R as an integer. */
1252 real_to_integer (const REAL_VALUE_TYPE *r)
1254 unsigned HOST_WIDE_INT i;
1265 i = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1271 if (REAL_EXP (r) <= 0)
1273 /* Only force overflow for unsigned overflow. Signed overflow is
1274 undefined, so it doesn't matter what we return, and some callers
1275 expect to be able to use this routine for both signed and
1276 unsigned conversions. */
1277 if (REAL_EXP (r) > HOST_BITS_PER_WIDE_INT)
1280 if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1281 i = r->sig[SIGSZ-1];
1282 else if (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG)
1284 i = r->sig[SIGSZ-1];
1285 i = i << (HOST_BITS_PER_LONG - 1) << 1;
1286 i |= r->sig[SIGSZ-2];
1291 i >>= HOST_BITS_PER_WIDE_INT - REAL_EXP (r);
1302 /* Likewise, but to an integer pair, HI+LOW. */
1305 real_to_integer2 (HOST_WIDE_INT *plow, HOST_WIDE_INT *phigh,
1306 const REAL_VALUE_TYPE *r)
1309 HOST_WIDE_INT low, high;
1322 high = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1336 /* Only force overflow for unsigned overflow. Signed overflow is
1337 undefined, so it doesn't matter what we return, and some callers
1338 expect to be able to use this routine for both signed and
1339 unsigned conversions. */
1340 if (exp > 2*HOST_BITS_PER_WIDE_INT)
1343 rshift_significand (&t, r, 2*HOST_BITS_PER_WIDE_INT - exp);
1344 if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1346 high = t.sig[SIGSZ-1];
1347 low = t.sig[SIGSZ-2];
1349 else if (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG)
1351 high = t.sig[SIGSZ-1];
1352 high = high << (HOST_BITS_PER_LONG - 1) << 1;
1353 high |= t.sig[SIGSZ-2];
1355 low = t.sig[SIGSZ-3];
1356 low = low << (HOST_BITS_PER_LONG - 1) << 1;
1357 low |= t.sig[SIGSZ-4];
1367 low = -low, high = ~high;
1379 /* A subroutine of real_to_decimal. Compute the quotient and remainder
1380 of NUM / DEN. Return the quotient and place the remainder in NUM.
1381 It is expected that NUM / DEN are close enough that the quotient is
1384 static unsigned long
1385 rtd_divmod (REAL_VALUE_TYPE *num, REAL_VALUE_TYPE *den)
1387 unsigned long q, msb;
1388 int expn = REAL_EXP (num), expd = REAL_EXP (den);
1397 msb = num->sig[SIGSZ-1] & SIG_MSB;
1399 lshift_significand_1 (num, num);
1401 if (msb || cmp_significands (num, den) >= 0)
1403 sub_significands (num, num, den, 0);
1407 while (--expn >= expd);
1409 SET_REAL_EXP (num, expd);
1415 /* Render R as a decimal floating point constant. Emit DIGITS significant
1416 digits in the result, bounded by BUF_SIZE. If DIGITS is 0, choose the
1417 maximum for the representation. If CROP_TRAILING_ZEROS, strip trailing
1420 #define M_LOG10_2 0.30102999566398119521
1423 real_to_decimal (char *str, const REAL_VALUE_TYPE *r_orig, size_t buf_size,
1424 size_t digits, int crop_trailing_zeros)
1426 const REAL_VALUE_TYPE *one, *ten;
1427 REAL_VALUE_TYPE r, pten, u, v;
1428 int dec_exp, cmp_one, digit;
1430 char *p, *first, *last;
1437 strcpy (str, (r.sign ? "-0.0" : "0.0"));
1442 strcpy (str, (r.sign ? "-Inf" : "+Inf"));
1445 /* ??? Print the significand as well, if not canonical? */
1446 strcpy (str, (r.sign ? "-NaN" : "+NaN"));
1452 /* Bound the number of digits printed by the size of the representation. */
1453 max_digits = SIGNIFICAND_BITS * M_LOG10_2;
1454 if (digits == 0 || digits > max_digits)
1455 digits = max_digits;
1457 /* Estimate the decimal exponent, and compute the length of the string it
1458 will print as. Be conservative and add one to account for possible
1459 overflow or rounding error. */
1460 dec_exp = REAL_EXP (&r) * M_LOG10_2;
1461 for (max_digits = 1; dec_exp ; max_digits++)
1464 /* Bound the number of digits printed by the size of the output buffer. */
1465 max_digits = buf_size - 1 - 1 - 2 - max_digits - 1;
1466 if (max_digits > buf_size)
1468 if (digits > max_digits)
1469 digits = max_digits;
1471 one = real_digit (1);
1472 ten = ten_to_ptwo (0);
1480 cmp_one = do_compare (&r, one, 0);
1485 /* Number is greater than one. Convert significand to an integer
1486 and strip trailing decimal zeros. */
1489 SET_REAL_EXP (&u, SIGNIFICAND_BITS - 1);
1491 /* Largest M, such that 10**2**M fits within SIGNIFICAND_BITS. */
1492 m = floor_log2 (max_digits);
1494 /* Iterate over the bits of the possible powers of 10 that might
1495 be present in U and eliminate them. That is, if we find that
1496 10**2**M divides U evenly, keep the division and increase
1502 do_divide (&t, &u, ten_to_ptwo (m));
1503 do_fix_trunc (&v, &t);
1504 if (cmp_significands (&v, &t) == 0)
1512 /* Revert the scaling to integer that we performed earlier. */
1513 SET_REAL_EXP (&u, REAL_EXP (&u) + REAL_EXP (&r)
1514 - (SIGNIFICAND_BITS - 1));
1517 /* Find power of 10. Do this by dividing out 10**2**M when
1518 this is larger than the current remainder. Fill PTEN with
1519 the power of 10 that we compute. */
1520 if (REAL_EXP (&r) > 0)
1522 m = floor_log2 ((int)(REAL_EXP (&r) * M_LOG10_2)) + 1;
1525 const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1526 if (do_compare (&u, ptentwo, 0) >= 0)
1528 do_divide (&u, &u, ptentwo);
1529 do_multiply (&pten, &pten, ptentwo);
1536 /* We managed to divide off enough tens in the above reduction
1537 loop that we've now got a negative exponent. Fall into the
1538 less-than-one code to compute the proper value for PTEN. */
1545 /* Number is less than one. Pad significand with leading
1551 /* Stop if we'd shift bits off the bottom. */
1555 do_multiply (&u, &v, ten);
1557 /* Stop if we're now >= 1. */
1558 if (REAL_EXP (&u) > 0)
1566 /* Find power of 10. Do this by multiplying in P=10**2**M when
1567 the current remainder is smaller than 1/P. Fill PTEN with the
1568 power of 10 that we compute. */
1569 m = floor_log2 ((int)(-REAL_EXP (&r) * M_LOG10_2)) + 1;
1572 const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1573 const REAL_VALUE_TYPE *ptenmtwo = ten_to_mptwo (m);
1575 if (do_compare (&v, ptenmtwo, 0) <= 0)
1577 do_multiply (&v, &v, ptentwo);
1578 do_multiply (&pten, &pten, ptentwo);
1584 /* Invert the positive power of 10 that we've collected so far. */
1585 do_divide (&pten, one, &pten);
1593 /* At this point, PTEN should contain the nearest power of 10 smaller
1594 than R, such that this division produces the first digit.
1596 Using a divide-step primitive that returns the complete integral
1597 remainder avoids the rounding error that would be produced if
1598 we were to use do_divide here and then simply multiply by 10 for
1599 each subsequent digit. */
1601 digit = rtd_divmod (&r, &pten);
1603 /* Be prepared for error in that division via underflow ... */
1604 if (digit == 0 && cmp_significand_0 (&r))
1606 /* Multiply by 10 and try again. */
1607 do_multiply (&r, &r, ten);
1608 digit = rtd_divmod (&r, &pten);
1614 /* ... or overflow. */
1622 else if (digit > 10)
1627 /* Generate subsequent digits. */
1628 while (--digits > 0)
1630 do_multiply (&r, &r, ten);
1631 digit = rtd_divmod (&r, &pten);
1636 /* Generate one more digit with which to do rounding. */
1637 do_multiply (&r, &r, ten);
1638 digit = rtd_divmod (&r, &pten);
1640 /* Round the result. */
1643 /* Round to nearest. If R is nonzero there are additional
1644 nonzero digits to be extracted. */
1645 if (cmp_significand_0 (&r))
1647 /* Round to even. */
1648 else if ((p[-1] - '0') & 1)
1665 /* Carry out of the first digit. This means we had all 9's and
1666 now have all 0's. "Prepend" a 1 by overwriting the first 0. */
1674 /* Insert the decimal point. */
1675 first[0] = first[1];
1678 /* If requested, drop trailing zeros. Never crop past "1.0". */
1679 if (crop_trailing_zeros)
1680 while (last > first + 3 && last[-1] == '0')
1683 /* Append the exponent. */
1684 sprintf (last, "e%+d", dec_exp);
1687 /* Render R as a hexadecimal floating point constant. Emit DIGITS
1688 significant digits in the result, bounded by BUF_SIZE. If DIGITS is 0,
1689 choose the maximum for the representation. If CROP_TRAILING_ZEROS,
1690 strip trailing zeros. */
1693 real_to_hexadecimal (char *str, const REAL_VALUE_TYPE *r, size_t buf_size,
1694 size_t digits, int crop_trailing_zeros)
1696 int i, j, exp = REAL_EXP (r);
1709 strcpy (str, (r->sign ? "-Inf" : "+Inf"));
1712 /* ??? Print the significand as well, if not canonical? */
1713 strcpy (str, (r->sign ? "-NaN" : "+NaN"));
1720 digits = SIGNIFICAND_BITS / 4;
1722 /* Bound the number of digits printed by the size of the output buffer. */
1724 sprintf (exp_buf, "p%+d", exp);
1725 max_digits = buf_size - strlen (exp_buf) - r->sign - 4 - 1;
1726 if (max_digits > buf_size)
1728 if (digits > max_digits)
1729 digits = max_digits;
1740 for (i = SIGSZ - 1; i >= 0; --i)
1741 for (j = HOST_BITS_PER_LONG - 4; j >= 0; j -= 4)
1743 *p++ = "0123456789abcdef"[(r->sig[i] >> j) & 15];
1749 if (crop_trailing_zeros)
1750 while (p > first + 1 && p[-1] == '0')
1753 sprintf (p, "p%+d", exp);
1756 /* Initialize R from a decimal or hexadecimal string. The string is
1757 assumed to have been syntax checked already. */
1760 real_from_string (REAL_VALUE_TYPE *r, const char *str)
1772 else if (*str == '+')
1775 if (str[0] == '0' && (str[1] == 'x' || str[1] == 'X'))
1777 /* Hexadecimal floating point. */
1778 int pos = SIGNIFICAND_BITS - 4, d;
1786 d = hex_value (*str);
1791 r->sig[pos / HOST_BITS_PER_LONG]
1792 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1801 if (pos == SIGNIFICAND_BITS - 4)
1808 d = hex_value (*str);
1813 r->sig[pos / HOST_BITS_PER_LONG]
1814 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1820 if (*str == 'p' || *str == 'P')
1822 bool exp_neg = false;
1830 else if (*str == '+')
1834 while (ISDIGIT (*str))
1840 /* Overflowed the exponent. */
1854 r->class = rvc_normal;
1855 SET_REAL_EXP (r, exp);
1861 /* Decimal floating point. */
1862 const REAL_VALUE_TYPE *ten = ten_to_ptwo (0);
1867 while (ISDIGIT (*str))
1870 do_multiply (r, r, ten);
1872 do_add (r, r, real_digit (d), 0);
1877 if (r->class == rvc_zero)
1882 while (ISDIGIT (*str))
1885 do_multiply (r, r, ten);
1887 do_add (r, r, real_digit (d), 0);
1892 if (*str == 'e' || *str == 'E')
1894 bool exp_neg = false;
1902 else if (*str == '+')
1906 while (ISDIGIT (*str))
1912 /* Overflowed the exponent. */
1926 times_pten (r, exp);
1941 /* Legacy. Similar, but return the result directly. */
1944 real_from_string2 (const char *s, enum machine_mode mode)
1948 real_from_string (&r, s);
1949 if (mode != VOIDmode)
1950 real_convert (&r, mode, &r);
1955 /* Initialize R from the integer pair HIGH+LOW. */
1958 real_from_integer (REAL_VALUE_TYPE *r, enum machine_mode mode,
1959 unsigned HOST_WIDE_INT low, HOST_WIDE_INT high,
1962 if (low == 0 && high == 0)
1966 r->class = rvc_normal;
1967 r->sign = high < 0 && !unsigned_p;
1968 SET_REAL_EXP (r, 2 * HOST_BITS_PER_WIDE_INT);
1979 if (HOST_BITS_PER_LONG == HOST_BITS_PER_WIDE_INT)
1981 r->sig[SIGSZ-1] = high;
1982 r->sig[SIGSZ-2] = low;
1983 memset (r->sig, 0, sizeof(long)*(SIGSZ-2));
1985 else if (HOST_BITS_PER_LONG*2 == HOST_BITS_PER_WIDE_INT)
1987 r->sig[SIGSZ-1] = high >> (HOST_BITS_PER_LONG - 1) >> 1;
1988 r->sig[SIGSZ-2] = high;
1989 r->sig[SIGSZ-3] = low >> (HOST_BITS_PER_LONG - 1) >> 1;
1990 r->sig[SIGSZ-4] = low;
1992 memset (r->sig, 0, sizeof(long)*(SIGSZ-4));
2000 if (mode != VOIDmode)
2001 real_convert (r, mode, r);
2004 /* Returns 10**2**N. */
2006 static const REAL_VALUE_TYPE *
2009 static REAL_VALUE_TYPE tens[EXP_BITS];
2011 if (n < 0 || n >= EXP_BITS)
2014 if (tens[n].class == rvc_zero)
2016 if (n < (HOST_BITS_PER_WIDE_INT == 64 ? 5 : 4))
2018 HOST_WIDE_INT t = 10;
2021 for (i = 0; i < n; ++i)
2024 real_from_integer (&tens[n], VOIDmode, t, 0, 1);
2028 const REAL_VALUE_TYPE *t = ten_to_ptwo (n - 1);
2029 do_multiply (&tens[n], t, t);
2036 /* Returns 10**(-2**N). */
2038 static const REAL_VALUE_TYPE *
2039 ten_to_mptwo (int n)
2041 static REAL_VALUE_TYPE tens[EXP_BITS];
2043 if (n < 0 || n >= EXP_BITS)
2046 if (tens[n].class == rvc_zero)
2047 do_divide (&tens[n], real_digit (1), ten_to_ptwo (n));
2054 static const REAL_VALUE_TYPE *
2057 static REAL_VALUE_TYPE num[10];
2062 if (n > 0 && num[n].class == rvc_zero)
2063 real_from_integer (&num[n], VOIDmode, n, 0, 1);
2068 /* Multiply R by 10**EXP. */
2071 times_pten (REAL_VALUE_TYPE *r, int exp)
2073 REAL_VALUE_TYPE pten, *rr;
2074 bool negative = (exp < 0);
2080 pten = *real_digit (1);
2086 for (i = 0; exp > 0; ++i, exp >>= 1)
2088 do_multiply (rr, rr, ten_to_ptwo (i));
2091 do_divide (r, r, &pten);
2094 /* Fills R with +Inf. */
2097 real_inf (REAL_VALUE_TYPE *r)
2102 /* Fills R with a NaN whose significand is described by STR. If QUIET,
2103 we force a QNaN, else we force an SNaN. The string, if not empty,
2104 is parsed as a number and placed in the significand. Return true
2105 if the string was successfully parsed. */
2108 real_nan (REAL_VALUE_TYPE *r, const char *str, int quiet,
2109 enum machine_mode mode)
2111 const struct real_format *fmt;
2113 fmt = REAL_MODE_FORMAT (mode);
2120 get_canonical_qnan (r, 0);
2122 get_canonical_snan (r, 0);
2129 memset (r, 0, sizeof (*r));
2132 /* Parse akin to strtol into the significand of R. */
2134 while (ISSPACE (*str))
2138 else if (*str == '+')
2148 while ((d = hex_value (*str)) < base)
2155 lshift_significand (r, r, 3);
2158 lshift_significand (r, r, 4);
2161 lshift_significand_1 (&u, r);
2162 lshift_significand (r, r, 3);
2163 add_significands (r, r, &u);
2171 add_significands (r, r, &u);
2176 /* Must have consumed the entire string for success. */
2180 /* Shift the significand into place such that the bits
2181 are in the most significant bits for the format. */
2182 lshift_significand (r, r, SIGNIFICAND_BITS - fmt->pnan);
2184 /* Our MSB is always unset for NaNs. */
2185 r->sig[SIGSZ-1] &= ~SIG_MSB;
2187 /* Force quiet or signalling NaN. */
2188 r->signalling = !quiet;
2194 /* Fills R with the largest finite value representable in mode MODE.
2195 If SIGN is nonzero, R is set to the most negative finite value. */
2198 real_maxval (REAL_VALUE_TYPE *r, int sign, enum machine_mode mode)
2200 const struct real_format *fmt;
2203 fmt = REAL_MODE_FORMAT (mode);
2207 r->class = rvc_normal;
2211 SET_REAL_EXP (r, fmt->emax * fmt->log2_b);
2213 np2 = SIGNIFICAND_BITS - fmt->p * fmt->log2_b;
2214 memset (r->sig, -1, SIGSZ * sizeof (unsigned long));
2215 clear_significand_below (r, np2);
2218 /* Fills R with 2**N. */
2221 real_2expN (REAL_VALUE_TYPE *r, int n)
2223 memset (r, 0, sizeof (*r));
2228 else if (n < -MAX_EXP)
2232 r->class = rvc_normal;
2233 SET_REAL_EXP (r, n);
2234 r->sig[SIGSZ-1] = SIG_MSB;
2240 round_for_format (const struct real_format *fmt, REAL_VALUE_TYPE *r)
2243 unsigned long sticky;
2247 p2 = fmt->p * fmt->log2_b;
2248 emin2m1 = (fmt->emin - 1) * fmt->log2_b;
2249 emax2 = fmt->emax * fmt->log2_b;
2251 np2 = SIGNIFICAND_BITS - p2;
2255 get_zero (r, r->sign);
2257 if (!fmt->has_signed_zero)
2262 get_inf (r, r->sign);
2267 clear_significand_below (r, np2);
2277 /* If we're not base2, normalize the exponent to a multiple of
2279 if (fmt->log2_b != 1)
2281 int shift = REAL_EXP (r) & (fmt->log2_b - 1);
2284 shift = fmt->log2_b - shift;
2285 r->sig[0] |= sticky_rshift_significand (r, r, shift);
2286 SET_REAL_EXP (r, REAL_EXP (r) + shift);
2290 /* Check the range of the exponent. If we're out of range,
2291 either underflow or overflow. */
2292 if (REAL_EXP (r) > emax2)
2294 else if (REAL_EXP (r) <= emin2m1)
2298 if (!fmt->has_denorm)
2300 /* Don't underflow completely until we've had a chance to round. */
2301 if (REAL_EXP (r) < emin2m1)
2306 diff = emin2m1 - REAL_EXP (r) + 1;
2310 /* De-normalize the significand. */
2311 r->sig[0] |= sticky_rshift_significand (r, r, diff);
2312 SET_REAL_EXP (r, REAL_EXP (r) + diff);
2316 /* There are P2 true significand bits, followed by one guard bit,
2317 followed by one sticky bit, followed by stuff. Fold nonzero
2318 stuff into the sticky bit. */
2321 for (i = 0, w = (np2 - 1) / HOST_BITS_PER_LONG; i < w; ++i)
2322 sticky |= r->sig[i];
2324 r->sig[w] & (((unsigned long)1 << ((np2 - 1) % HOST_BITS_PER_LONG)) - 1);
2326 guard = test_significand_bit (r, np2 - 1);
2327 lsb = test_significand_bit (r, np2);
2329 /* Round to even. */
2330 if (guard && (sticky || lsb))
2334 set_significand_bit (&u, np2);
2336 if (add_significands (r, r, &u))
2338 /* Overflow. Means the significand had been all ones, and
2339 is now all zeros. Need to increase the exponent, and
2340 possibly re-normalize it. */
2341 SET_REAL_EXP (r, REAL_EXP (r) + 1);
2342 if (REAL_EXP (r) > emax2)
2344 r->sig[SIGSZ-1] = SIG_MSB;
2346 if (fmt->log2_b != 1)
2348 int shift = REAL_EXP (r) & (fmt->log2_b - 1);
2351 shift = fmt->log2_b - shift;
2352 rshift_significand (r, r, shift);
2353 SET_REAL_EXP (r, REAL_EXP (r) + shift);
2354 if (REAL_EXP (r) > emax2)
2361 /* Catch underflow that we deferred until after rounding. */
2362 if (REAL_EXP (r) <= emin2m1)
2365 /* Clear out trailing garbage. */
2366 clear_significand_below (r, np2);
2369 /* Extend or truncate to a new mode. */
2372 real_convert (REAL_VALUE_TYPE *r, enum machine_mode mode,
2373 const REAL_VALUE_TYPE *a)
2375 const struct real_format *fmt;
2377 fmt = REAL_MODE_FORMAT (mode);
2382 round_for_format (fmt, r);
2384 /* round_for_format de-normalizes denormals. Undo just that part. */
2385 if (r->class == rvc_normal)
2389 /* Legacy. Likewise, except return the struct directly. */
2392 real_value_truncate (enum machine_mode mode, REAL_VALUE_TYPE a)
2395 real_convert (&r, mode, &a);
2399 /* Return true if truncating to MODE is exact. */
2402 exact_real_truncate (enum machine_mode mode, const REAL_VALUE_TYPE *a)
2405 real_convert (&t, mode, a);
2406 return real_identical (&t, a);
2409 /* Write R to the given target format. Place the words of the result
2410 in target word order in BUF. There are always 32 bits in each
2411 long, no matter the size of the host long.
2413 Legacy: return word 0 for implementing REAL_VALUE_TO_TARGET_SINGLE. */
2416 real_to_target_fmt (long *buf, const REAL_VALUE_TYPE *r_orig,
2417 const struct real_format *fmt)
2423 round_for_format (fmt, &r);
2427 (*fmt->encode) (fmt, buf, &r);
2432 /* Similar, but look up the format from MODE. */
2435 real_to_target (long *buf, const REAL_VALUE_TYPE *r, enum machine_mode mode)
2437 const struct real_format *fmt;
2439 fmt = REAL_MODE_FORMAT (mode);
2443 return real_to_target_fmt (buf, r, fmt);
2446 /* Read R from the given target format. Read the words of the result
2447 in target word order in BUF. There are always 32 bits in each
2448 long, no matter the size of the host long. */
2451 real_from_target_fmt (REAL_VALUE_TYPE *r, const long *buf,
2452 const struct real_format *fmt)
2454 (*fmt->decode) (fmt, r, buf);
2457 /* Similar, but look up the format from MODE. */
2460 real_from_target (REAL_VALUE_TYPE *r, const long *buf, enum machine_mode mode)
2462 const struct real_format *fmt;
2464 fmt = REAL_MODE_FORMAT (mode);
2468 (*fmt->decode) (fmt, r, buf);
2471 /* Return the number of bits in the significand for MODE. */
2472 /* ??? Legacy. Should get access to real_format directly. */
2475 significand_size (enum machine_mode mode)
2477 const struct real_format *fmt;
2479 fmt = REAL_MODE_FORMAT (mode);
2483 return fmt->p * fmt->log2_b;
2486 /* Return a hash value for the given real value. */
2487 /* ??? The "unsigned int" return value is intended to be hashval_t,
2488 but I didn't want to pull hashtab.h into real.h. */
2491 real_hash (const REAL_VALUE_TYPE *r)
2496 h = r->class | (r->sign << 2);
2504 h |= REAL_EXP (r) << 3;
2509 h ^= (unsigned int)-1;
2518 if (sizeof(unsigned long) > sizeof(unsigned int))
2519 for (i = 0; i < SIGSZ; ++i)
2521 unsigned long s = r->sig[i];
2522 h ^= s ^ (s >> (HOST_BITS_PER_LONG / 2));
2525 for (i = 0; i < SIGSZ; ++i)
2531 /* IEEE single-precision format. */
2533 static void encode_ieee_single (const struct real_format *fmt,
2534 long *, const REAL_VALUE_TYPE *);
2535 static void decode_ieee_single (const struct real_format *,
2536 REAL_VALUE_TYPE *, const long *);
2539 encode_ieee_single (const struct real_format *fmt, long *buf,
2540 const REAL_VALUE_TYPE *r)
2542 unsigned long image, sig, exp;
2543 unsigned long sign = r->sign;
2544 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2547 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
2558 image |= 0x7fffffff;
2566 if (r->signalling == fmt->qnan_msb_set)
2570 /* We overload qnan_msb_set here: it's only clear for
2571 mips_ieee_single, which wants all mantissa bits but the
2572 quiet/signalling one set in canonical NaNs (at least
2574 if (r->canonical && !fmt->qnan_msb_set)
2575 sig |= (1 << 22) - 1;
2583 image |= 0x7fffffff;
2587 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2588 whereas the intermediate representation is 0.F x 2**exp.
2589 Which means we're off by one. */
2593 exp = REAL_EXP (r) + 127 - 1;
2606 decode_ieee_single (const struct real_format *fmt, REAL_VALUE_TYPE *r,
2609 unsigned long image = buf[0] & 0xffffffff;
2610 bool sign = (image >> 31) & 1;
2611 int exp = (image >> 23) & 0xff;
2613 memset (r, 0, sizeof (*r));
2614 image <<= HOST_BITS_PER_LONG - 24;
2619 if (image && fmt->has_denorm)
2621 r->class = rvc_normal;
2623 SET_REAL_EXP (r, -126);
2624 r->sig[SIGSZ-1] = image << 1;
2627 else if (fmt->has_signed_zero)
2630 else if (exp == 255 && (fmt->has_nans || fmt->has_inf))
2636 r->signalling = (((image >> (HOST_BITS_PER_LONG - 2)) & 1)
2637 ^ fmt->qnan_msb_set);
2638 r->sig[SIGSZ-1] = image;
2648 r->class = rvc_normal;
2650 SET_REAL_EXP (r, exp - 127 + 1);
2651 r->sig[SIGSZ-1] = image | SIG_MSB;
2655 const struct real_format ieee_single_format =
2673 const struct real_format mips_single_format =
2692 /* IEEE double-precision format. */
2694 static void encode_ieee_double (const struct real_format *fmt,
2695 long *, const REAL_VALUE_TYPE *);
2696 static void decode_ieee_double (const struct real_format *,
2697 REAL_VALUE_TYPE *, const long *);
2700 encode_ieee_double (const struct real_format *fmt, long *buf,
2701 const REAL_VALUE_TYPE *r)
2703 unsigned long image_lo, image_hi, sig_lo, sig_hi, exp;
2704 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2706 image_hi = r->sign << 31;
2709 if (HOST_BITS_PER_LONG == 64)
2711 sig_hi = r->sig[SIGSZ-1];
2712 sig_lo = (sig_hi >> (64 - 53)) & 0xffffffff;
2713 sig_hi = (sig_hi >> (64 - 53 + 1) >> 31) & 0xfffff;
2717 sig_hi = r->sig[SIGSZ-1];
2718 sig_lo = r->sig[SIGSZ-2];
2719 sig_lo = (sig_hi << 21) | (sig_lo >> 11);
2720 sig_hi = (sig_hi >> 11) & 0xfffff;
2730 image_hi |= 2047 << 20;
2733 image_hi |= 0x7fffffff;
2734 image_lo = 0xffffffff;
2742 sig_hi = sig_lo = 0;
2743 if (r->signalling == fmt->qnan_msb_set)
2744 sig_hi &= ~(1 << 19);
2747 /* We overload qnan_msb_set here: it's only clear for
2748 mips_ieee_single, which wants all mantissa bits but the
2749 quiet/signalling one set in canonical NaNs (at least
2751 if (r->canonical && !fmt->qnan_msb_set)
2753 sig_hi |= (1 << 19) - 1;
2754 sig_lo = 0xffffffff;
2756 else if (sig_hi == 0 && sig_lo == 0)
2759 image_hi |= 2047 << 20;
2765 image_hi |= 0x7fffffff;
2766 image_lo = 0xffffffff;
2771 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2772 whereas the intermediate representation is 0.F x 2**exp.
2773 Which means we're off by one. */
2777 exp = REAL_EXP (r) + 1023 - 1;
2778 image_hi |= exp << 20;
2787 if (FLOAT_WORDS_BIG_ENDIAN)
2788 buf[0] = image_hi, buf[1] = image_lo;
2790 buf[0] = image_lo, buf[1] = image_hi;
2794 decode_ieee_double (const struct real_format *fmt, REAL_VALUE_TYPE *r,
2797 unsigned long image_hi, image_lo;
2801 if (FLOAT_WORDS_BIG_ENDIAN)
2802 image_hi = buf[0], image_lo = buf[1];
2804 image_lo = buf[0], image_hi = buf[1];
2805 image_lo &= 0xffffffff;
2806 image_hi &= 0xffffffff;
2808 sign = (image_hi >> 31) & 1;
2809 exp = (image_hi >> 20) & 0x7ff;
2811 memset (r, 0, sizeof (*r));
2813 image_hi <<= 32 - 21;
2814 image_hi |= image_lo >> 21;
2815 image_hi &= 0x7fffffff;
2816 image_lo <<= 32 - 21;
2820 if ((image_hi || image_lo) && fmt->has_denorm)
2822 r->class = rvc_normal;
2824 SET_REAL_EXP (r, -1022);
2825 if (HOST_BITS_PER_LONG == 32)
2827 image_hi = (image_hi << 1) | (image_lo >> 31);
2829 r->sig[SIGSZ-1] = image_hi;
2830 r->sig[SIGSZ-2] = image_lo;
2834 image_hi = (image_hi << 31 << 2) | (image_lo << 1);
2835 r->sig[SIGSZ-1] = image_hi;
2839 else if (fmt->has_signed_zero)
2842 else if (exp == 2047 && (fmt->has_nans || fmt->has_inf))
2844 if (image_hi || image_lo)
2848 r->signalling = ((image_hi >> 30) & 1) ^ fmt->qnan_msb_set;
2849 if (HOST_BITS_PER_LONG == 32)
2851 r->sig[SIGSZ-1] = image_hi;
2852 r->sig[SIGSZ-2] = image_lo;
2855 r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo;
2865 r->class = rvc_normal;
2867 SET_REAL_EXP (r, exp - 1023 + 1);
2868 if (HOST_BITS_PER_LONG == 32)
2870 r->sig[SIGSZ-1] = image_hi | SIG_MSB;
2871 r->sig[SIGSZ-2] = image_lo;
2874 r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo | SIG_MSB;
2878 const struct real_format ieee_double_format =
2896 const struct real_format mips_double_format =
2915 /* IEEE extended double precision format. This comes in three
2916 flavors: Intel's as a 12 byte image, Intel's as a 16 byte image,
2919 static void encode_ieee_extended (const struct real_format *fmt,
2920 long *, const REAL_VALUE_TYPE *);
2921 static void decode_ieee_extended (const struct real_format *,
2922 REAL_VALUE_TYPE *, const long *);
2924 static void encode_ieee_extended_128 (const struct real_format *fmt,
2925 long *, const REAL_VALUE_TYPE *);
2926 static void decode_ieee_extended_128 (const struct real_format *,
2927 REAL_VALUE_TYPE *, const long *);
2930 encode_ieee_extended (const struct real_format *fmt, long *buf,
2931 const REAL_VALUE_TYPE *r)
2933 unsigned long image_hi, sig_hi, sig_lo;
2934 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2936 image_hi = r->sign << 15;
2937 sig_hi = sig_lo = 0;
2949 /* Intel requires the explicit integer bit to be set, otherwise
2950 it considers the value a "pseudo-infinity". Motorola docs
2951 say it doesn't care. */
2952 sig_hi = 0x80000000;
2957 sig_lo = sig_hi = 0xffffffff;
2965 if (HOST_BITS_PER_LONG == 32)
2967 sig_hi = r->sig[SIGSZ-1];
2968 sig_lo = r->sig[SIGSZ-2];
2972 sig_lo = r->sig[SIGSZ-1];
2973 sig_hi = sig_lo >> 31 >> 1;
2974 sig_lo &= 0xffffffff;
2976 if (r->signalling == fmt->qnan_msb_set)
2977 sig_hi &= ~(1 << 30);
2980 if ((sig_hi & 0x7fffffff) == 0 && sig_lo == 0)
2983 /* Intel requires the explicit integer bit to be set, otherwise
2984 it considers the value a "pseudo-nan". Motorola docs say it
2986 sig_hi |= 0x80000000;
2991 sig_lo = sig_hi = 0xffffffff;
2997 int exp = REAL_EXP (r);
2999 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3000 whereas the intermediate representation is 0.F x 2**exp.
3001 Which means we're off by one.
3003 Except for Motorola, which consider exp=0 and explicit
3004 integer bit set to continue to be normalized. In theory
3005 this discrepancy has been taken care of by the difference
3006 in fmt->emin in round_for_format. */
3018 if (HOST_BITS_PER_LONG == 32)
3020 sig_hi = r->sig[SIGSZ-1];
3021 sig_lo = r->sig[SIGSZ-2];
3025 sig_lo = r->sig[SIGSZ-1];
3026 sig_hi = sig_lo >> 31 >> 1;
3027 sig_lo &= 0xffffffff;
3036 if (FLOAT_WORDS_BIG_ENDIAN)
3037 buf[0] = image_hi << 16, buf[1] = sig_hi, buf[2] = sig_lo;
3039 buf[0] = sig_lo, buf[1] = sig_hi, buf[2] = image_hi;
3043 encode_ieee_extended_128 (const struct real_format *fmt, long *buf,
3044 const REAL_VALUE_TYPE *r)
3046 buf[3 * !FLOAT_WORDS_BIG_ENDIAN] = 0;
3047 encode_ieee_extended (fmt, buf+!!FLOAT_WORDS_BIG_ENDIAN, r);
3051 decode_ieee_extended (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3054 unsigned long image_hi, sig_hi, sig_lo;
3058 if (FLOAT_WORDS_BIG_ENDIAN)
3059 image_hi = buf[0] >> 16, sig_hi = buf[1], sig_lo = buf[2];
3061 sig_lo = buf[0], sig_hi = buf[1], image_hi = buf[2];
3062 sig_lo &= 0xffffffff;
3063 sig_hi &= 0xffffffff;
3064 image_hi &= 0xffffffff;
3066 sign = (image_hi >> 15) & 1;
3067 exp = image_hi & 0x7fff;
3069 memset (r, 0, sizeof (*r));
3073 if ((sig_hi || sig_lo) && fmt->has_denorm)
3075 r->class = rvc_normal;
3078 /* When the IEEE format contains a hidden bit, we know that
3079 it's zero at this point, and so shift up the significand
3080 and decrease the exponent to match. In this case, Motorola
3081 defines the explicit integer bit to be valid, so we don't
3082 know whether the msb is set or not. */
3083 SET_REAL_EXP (r, fmt->emin);
3084 if (HOST_BITS_PER_LONG == 32)
3086 r->sig[SIGSZ-1] = sig_hi;
3087 r->sig[SIGSZ-2] = sig_lo;
3090 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3094 else if (fmt->has_signed_zero)
3097 else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3099 /* See above re "pseudo-infinities" and "pseudo-nans".
3100 Short summary is that the MSB will likely always be
3101 set, and that we don't care about it. */
3102 sig_hi &= 0x7fffffff;
3104 if (sig_hi || sig_lo)
3108 r->signalling = ((sig_hi >> 30) & 1) ^ fmt->qnan_msb_set;
3109 if (HOST_BITS_PER_LONG == 32)
3111 r->sig[SIGSZ-1] = sig_hi;
3112 r->sig[SIGSZ-2] = sig_lo;
3115 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3125 r->class = rvc_normal;
3127 SET_REAL_EXP (r, exp - 16383 + 1);
3128 if (HOST_BITS_PER_LONG == 32)
3130 r->sig[SIGSZ-1] = sig_hi;
3131 r->sig[SIGSZ-2] = sig_lo;
3134 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3139 decode_ieee_extended_128 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3142 decode_ieee_extended (fmt, r, buf+!!FLOAT_WORDS_BIG_ENDIAN);
3145 const struct real_format ieee_extended_motorola_format =
3147 encode_ieee_extended,
3148 decode_ieee_extended,
3163 const struct real_format ieee_extended_intel_96_format =
3165 encode_ieee_extended,
3166 decode_ieee_extended,
3181 const struct real_format ieee_extended_intel_128_format =
3183 encode_ieee_extended_128,
3184 decode_ieee_extended_128,
3199 /* The following caters to i386 systems that set the rounding precision
3200 to 53 bits instead of 64, e.g. FreeBSD. */
3201 const struct real_format ieee_extended_intel_96_round_53_format =
3203 encode_ieee_extended,
3204 decode_ieee_extended,
3219 /* IBM 128-bit extended precision format: a pair of IEEE double precision
3220 numbers whose sum is equal to the extended precision value. The number
3221 with greater magnitude is first. This format has the same magnitude
3222 range as an IEEE double precision value, but effectively 106 bits of
3223 significand precision. Infinity and NaN are represented by their IEEE
3224 double precision value stored in the first number, the second number is
3225 ignored. Zeroes, Infinities, and NaNs are set in both doubles
3226 due to precedent. */
3228 static void encode_ibm_extended (const struct real_format *fmt,
3229 long *, const REAL_VALUE_TYPE *);
3230 static void decode_ibm_extended (const struct real_format *,
3231 REAL_VALUE_TYPE *, const long *);
3234 encode_ibm_extended (const struct real_format *fmt, long *buf,
3235 const REAL_VALUE_TYPE *r)
3237 REAL_VALUE_TYPE u, normr, v;
3238 const struct real_format *base_fmt;
3240 base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3242 /* Renormlize R before doing any arithmetic on it. */
3244 if (normr.class == rvc_normal)
3247 /* u = IEEE double precision portion of significand. */
3249 round_for_format (base_fmt, &u);
3250 encode_ieee_double (base_fmt, &buf[0], &u);
3252 if (u.class == rvc_normal)
3254 do_add (&v, &normr, &u, 1);
3255 /* Call round_for_format since we might need to denormalize. */
3256 round_for_format (base_fmt, &v);
3257 encode_ieee_double (base_fmt, &buf[2], &v);
3261 /* Inf, NaN, 0 are all representable as doubles, so the
3262 least-significant part can be 0.0. */
3269 decode_ibm_extended (const struct real_format *fmt ATTRIBUTE_UNUSED, REAL_VALUE_TYPE *r,
3272 REAL_VALUE_TYPE u, v;
3273 const struct real_format *base_fmt;
3275 base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3276 decode_ieee_double (base_fmt, &u, &buf[0]);
3278 if (u.class != rvc_zero && u.class != rvc_inf && u.class != rvc_nan)
3280 decode_ieee_double (base_fmt, &v, &buf[2]);
3281 do_add (r, &u, &v, 0);
3287 const struct real_format ibm_extended_format =
3289 encode_ibm_extended,
3290 decode_ibm_extended,
3305 const struct real_format mips_extended_format =
3307 encode_ibm_extended,
3308 decode_ibm_extended,
3324 /* IEEE quad precision format. */
3326 static void encode_ieee_quad (const struct real_format *fmt,
3327 long *, const REAL_VALUE_TYPE *);
3328 static void decode_ieee_quad (const struct real_format *,
3329 REAL_VALUE_TYPE *, const long *);
3332 encode_ieee_quad (const struct real_format *fmt, long *buf,
3333 const REAL_VALUE_TYPE *r)
3335 unsigned long image3, image2, image1, image0, exp;
3336 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3339 image3 = r->sign << 31;
3344 rshift_significand (&u, r, SIGNIFICAND_BITS - 113);
3353 image3 |= 32767 << 16;
3356 image3 |= 0x7fffffff;
3357 image2 = 0xffffffff;
3358 image1 = 0xffffffff;
3359 image0 = 0xffffffff;
3366 image3 |= 32767 << 16;
3370 /* Don't use bits from the significand. The
3371 initialization above is right. */
3373 else if (HOST_BITS_PER_LONG == 32)
3378 image3 |= u.sig[3] & 0xffff;
3383 image1 = image0 >> 31 >> 1;
3385 image3 |= (image2 >> 31 >> 1) & 0xffff;
3386 image0 &= 0xffffffff;
3387 image2 &= 0xffffffff;
3389 if (r->signalling == fmt->qnan_msb_set)
3393 /* We overload qnan_msb_set here: it's only clear for
3394 mips_ieee_single, which wants all mantissa bits but the
3395 quiet/signalling one set in canonical NaNs (at least
3397 if (r->canonical && !fmt->qnan_msb_set)
3400 image2 = image1 = image0 = 0xffffffff;
3402 else if (((image3 & 0xffff) | image2 | image1 | image0) == 0)
3407 image3 |= 0x7fffffff;
3408 image2 = 0xffffffff;
3409 image1 = 0xffffffff;
3410 image0 = 0xffffffff;
3415 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3416 whereas the intermediate representation is 0.F x 2**exp.
3417 Which means we're off by one. */
3421 exp = REAL_EXP (r) + 16383 - 1;
3422 image3 |= exp << 16;
3424 if (HOST_BITS_PER_LONG == 32)
3429 image3 |= u.sig[3] & 0xffff;
3434 image1 = image0 >> 31 >> 1;
3436 image3 |= (image2 >> 31 >> 1) & 0xffff;
3437 image0 &= 0xffffffff;
3438 image2 &= 0xffffffff;
3446 if (FLOAT_WORDS_BIG_ENDIAN)
3463 decode_ieee_quad (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3466 unsigned long image3, image2, image1, image0;
3470 if (FLOAT_WORDS_BIG_ENDIAN)
3484 image0 &= 0xffffffff;
3485 image1 &= 0xffffffff;
3486 image2 &= 0xffffffff;
3488 sign = (image3 >> 31) & 1;
3489 exp = (image3 >> 16) & 0x7fff;
3492 memset (r, 0, sizeof (*r));
3496 if ((image3 | image2 | image1 | image0) && fmt->has_denorm)
3498 r->class = rvc_normal;
3501 SET_REAL_EXP (r, -16382 + (SIGNIFICAND_BITS - 112));
3502 if (HOST_BITS_PER_LONG == 32)
3511 r->sig[0] = (image1 << 31 << 1) | image0;
3512 r->sig[1] = (image3 << 31 << 1) | image2;
3517 else if (fmt->has_signed_zero)
3520 else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3522 if (image3 | image2 | image1 | image0)
3526 r->signalling = ((image3 >> 15) & 1) ^ fmt->qnan_msb_set;
3528 if (HOST_BITS_PER_LONG == 32)
3537 r->sig[0] = (image1 << 31 << 1) | image0;
3538 r->sig[1] = (image3 << 31 << 1) | image2;
3540 lshift_significand (r, r, SIGNIFICAND_BITS - 113);
3550 r->class = rvc_normal;
3552 SET_REAL_EXP (r, exp - 16383 + 1);
3554 if (HOST_BITS_PER_LONG == 32)
3563 r->sig[0] = (image1 << 31 << 1) | image0;
3564 r->sig[1] = (image3 << 31 << 1) | image2;
3566 lshift_significand (r, r, SIGNIFICAND_BITS - 113);
3567 r->sig[SIGSZ-1] |= SIG_MSB;
3571 const struct real_format ieee_quad_format =
3589 const struct real_format mips_quad_format =
3607 /* Descriptions of VAX floating point formats can be found beginning at
3609 http://h71000.www7.hp.com/doc/73FINAL/4515/4515pro_013.html#f_floating_point_format
3611 The thing to remember is that they're almost IEEE, except for word
3612 order, exponent bias, and the lack of infinities, nans, and denormals.
3614 We don't implement the H_floating format here, simply because neither
3615 the VAX or Alpha ports use it. */
3617 static void encode_vax_f (const struct real_format *fmt,
3618 long *, const REAL_VALUE_TYPE *);
3619 static void decode_vax_f (const struct real_format *,
3620 REAL_VALUE_TYPE *, const long *);
3621 static void encode_vax_d (const struct real_format *fmt,
3622 long *, const REAL_VALUE_TYPE *);
3623 static void decode_vax_d (const struct real_format *,
3624 REAL_VALUE_TYPE *, const long *);
3625 static void encode_vax_g (const struct real_format *fmt,
3626 long *, const REAL_VALUE_TYPE *);
3627 static void decode_vax_g (const struct real_format *,
3628 REAL_VALUE_TYPE *, const long *);
3631 encode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3632 const REAL_VALUE_TYPE *r)
3634 unsigned long sign, exp, sig, image;
3636 sign = r->sign << 15;
3646 image = 0xffff7fff | sign;
3650 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
3651 exp = REAL_EXP (r) + 128;
3653 image = (sig << 16) & 0xffff0000;
3667 decode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED,
3668 REAL_VALUE_TYPE *r, const long *buf)
3670 unsigned long image = buf[0] & 0xffffffff;
3671 int exp = (image >> 7) & 0xff;
3673 memset (r, 0, sizeof (*r));
3677 r->class = rvc_normal;
3678 r->sign = (image >> 15) & 1;
3679 SET_REAL_EXP (r, exp - 128);
3681 image = ((image & 0x7f) << 16) | ((image >> 16) & 0xffff);
3682 r->sig[SIGSZ-1] = (image << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
3687 encode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3688 const REAL_VALUE_TYPE *r)
3690 unsigned long image0, image1, sign = r->sign << 15;
3695 image0 = image1 = 0;
3700 image0 = 0xffff7fff | sign;
3701 image1 = 0xffffffff;
3705 /* Extract the significand into straight hi:lo. */
3706 if (HOST_BITS_PER_LONG == 64)
3708 image0 = r->sig[SIGSZ-1];
3709 image1 = (image0 >> (64 - 56)) & 0xffffffff;
3710 image0 = (image0 >> (64 - 56 + 1) >> 31) & 0x7fffff;
3714 image0 = r->sig[SIGSZ-1];
3715 image1 = r->sig[SIGSZ-2];
3716 image1 = (image0 << 24) | (image1 >> 8);
3717 image0 = (image0 >> 8) & 0xffffff;
3720 /* Rearrange the half-words of the significand to match the
3722 image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff007f;
3723 image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
3725 /* Add the sign and exponent. */
3727 image0 |= (REAL_EXP (r) + 128) << 7;
3734 if (FLOAT_WORDS_BIG_ENDIAN)
3735 buf[0] = image1, buf[1] = image0;
3737 buf[0] = image0, buf[1] = image1;
3741 decode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED,
3742 REAL_VALUE_TYPE *r, const long *buf)
3744 unsigned long image0, image1;
3747 if (FLOAT_WORDS_BIG_ENDIAN)
3748 image1 = buf[0], image0 = buf[1];
3750 image0 = buf[0], image1 = buf[1];
3751 image0 &= 0xffffffff;
3752 image1 &= 0xffffffff;
3754 exp = (image0 >> 7) & 0xff;
3756 memset (r, 0, sizeof (*r));
3760 r->class = rvc_normal;
3761 r->sign = (image0 >> 15) & 1;
3762 SET_REAL_EXP (r, exp - 128);
3764 /* Rearrange the half-words of the external format into
3765 proper ascending order. */
3766 image0 = ((image0 & 0x7f) << 16) | ((image0 >> 16) & 0xffff);
3767 image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
3769 if (HOST_BITS_PER_LONG == 64)
3771 image0 = (image0 << 31 << 1) | image1;
3774 r->sig[SIGSZ-1] = image0;
3778 r->sig[SIGSZ-1] = image0;
3779 r->sig[SIGSZ-2] = image1;
3780 lshift_significand (r, r, 2*HOST_BITS_PER_LONG - 56);
3781 r->sig[SIGSZ-1] |= SIG_MSB;
3787 encode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3788 const REAL_VALUE_TYPE *r)
3790 unsigned long image0, image1, sign = r->sign << 15;
3795 image0 = image1 = 0;
3800 image0 = 0xffff7fff | sign;
3801 image1 = 0xffffffff;
3805 /* Extract the significand into straight hi:lo. */
3806 if (HOST_BITS_PER_LONG == 64)
3808 image0 = r->sig[SIGSZ-1];
3809 image1 = (image0 >> (64 - 53)) & 0xffffffff;
3810 image0 = (image0 >> (64 - 53 + 1) >> 31) & 0xfffff;
3814 image0 = r->sig[SIGSZ-1];
3815 image1 = r->sig[SIGSZ-2];
3816 image1 = (image0 << 21) | (image1 >> 11);
3817 image0 = (image0 >> 11) & 0xfffff;
3820 /* Rearrange the half-words of the significand to match the
3822 image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff000f;
3823 image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
3825 /* Add the sign and exponent. */
3827 image0 |= (REAL_EXP (r) + 1024) << 4;
3834 if (FLOAT_WORDS_BIG_ENDIAN)
3835 buf[0] = image1, buf[1] = image0;
3837 buf[0] = image0, buf[1] = image1;
3841 decode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED,
3842 REAL_VALUE_TYPE *r, const long *buf)
3844 unsigned long image0, image1;
3847 if (FLOAT_WORDS_BIG_ENDIAN)
3848 image1 = buf[0], image0 = buf[1];
3850 image0 = buf[0], image1 = buf[1];
3851 image0 &= 0xffffffff;
3852 image1 &= 0xffffffff;
3854 exp = (image0 >> 4) & 0x7ff;
3856 memset (r, 0, sizeof (*r));
3860 r->class = rvc_normal;
3861 r->sign = (image0 >> 15) & 1;
3862 SET_REAL_EXP (r, exp - 1024);
3864 /* Rearrange the half-words of the external format into
3865 proper ascending order. */
3866 image0 = ((image0 & 0xf) << 16) | ((image0 >> 16) & 0xffff);
3867 image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
3869 if (HOST_BITS_PER_LONG == 64)
3871 image0 = (image0 << 31 << 1) | image1;
3874 r->sig[SIGSZ-1] = image0;
3878 r->sig[SIGSZ-1] = image0;
3879 r->sig[SIGSZ-2] = image1;
3880 lshift_significand (r, r, 64 - 53);
3881 r->sig[SIGSZ-1] |= SIG_MSB;
3886 const struct real_format vax_f_format =
3904 const struct real_format vax_d_format =
3922 const struct real_format vax_g_format =
3940 /* A good reference for these can be found in chapter 9 of
3941 "ESA/390 Principles of Operation", IBM document number SA22-7201-01.
3942 An on-line version can be found here:
3944 http://publibz.boulder.ibm.com/cgi-bin/bookmgr_OS390/BOOKS/DZ9AR001/9.1?DT=19930923083613
3947 static void encode_i370_single (const struct real_format *fmt,
3948 long *, const REAL_VALUE_TYPE *);
3949 static void decode_i370_single (const struct real_format *,
3950 REAL_VALUE_TYPE *, const long *);
3951 static void encode_i370_double (const struct real_format *fmt,
3952 long *, const REAL_VALUE_TYPE *);
3953 static void decode_i370_double (const struct real_format *,
3954 REAL_VALUE_TYPE *, const long *);
3957 encode_i370_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
3958 long *buf, const REAL_VALUE_TYPE *r)
3960 unsigned long sign, exp, sig, image;
3962 sign = r->sign << 31;
3972 image = 0x7fffffff | sign;
3976 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0xffffff;
3977 exp = ((REAL_EXP (r) / 4) + 64) << 24;
3978 image = sign | exp | sig;
3989 decode_i370_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
3990 REAL_VALUE_TYPE *r, const long *buf)
3992 unsigned long sign, sig, image = buf[0];
3995 sign = (image >> 31) & 1;
3996 exp = (image >> 24) & 0x7f;
3997 sig = image & 0xffffff;
3999 memset (r, 0, sizeof (*r));
4003 r->class = rvc_normal;
4005 SET_REAL_EXP (r, (exp - 64) * 4);
4006 r->sig[SIGSZ-1] = sig << (HOST_BITS_PER_LONG - 24);
4012 encode_i370_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4013 long *buf, const REAL_VALUE_TYPE *r)
4015 unsigned long sign, exp, image_hi, image_lo;
4017 sign = r->sign << 31;
4022 image_hi = image_lo = 0;
4027 image_hi = 0x7fffffff | sign;
4028 image_lo = 0xffffffff;
4032 if (HOST_BITS_PER_LONG == 64)
4034 image_hi = r->sig[SIGSZ-1];
4035 image_lo = (image_hi >> (64 - 56)) & 0xffffffff;
4036 image_hi = (image_hi >> (64 - 56 + 1) >> 31) & 0xffffff;
4040 image_hi = r->sig[SIGSZ-1];
4041 image_lo = r->sig[SIGSZ-2];
4042 image_lo = (image_lo >> 8) | (image_hi << 24);
4046 exp = ((REAL_EXP (r) / 4) + 64) << 24;
4047 image_hi |= sign | exp;
4054 if (FLOAT_WORDS_BIG_ENDIAN)
4055 buf[0] = image_hi, buf[1] = image_lo;
4057 buf[0] = image_lo, buf[1] = image_hi;
4061 decode_i370_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4062 REAL_VALUE_TYPE *r, const long *buf)
4064 unsigned long sign, image_hi, image_lo;
4067 if (FLOAT_WORDS_BIG_ENDIAN)
4068 image_hi = buf[0], image_lo = buf[1];
4070 image_lo = buf[0], image_hi = buf[1];
4072 sign = (image_hi >> 31) & 1;
4073 exp = (image_hi >> 24) & 0x7f;
4074 image_hi &= 0xffffff;
4075 image_lo &= 0xffffffff;
4077 memset (r, 0, sizeof (*r));
4079 if (exp || image_hi || image_lo)
4081 r->class = rvc_normal;
4083 SET_REAL_EXP (r, (exp - 64) * 4 + (SIGNIFICAND_BITS - 56));
4085 if (HOST_BITS_PER_LONG == 32)
4087 r->sig[0] = image_lo;
4088 r->sig[1] = image_hi;
4091 r->sig[0] = image_lo | (image_hi << 31 << 1);
4097 const struct real_format i370_single_format =
4110 false, /* ??? The encoding does allow for "unnormals". */
4111 false, /* ??? The encoding does allow for "unnormals". */
4115 const struct real_format i370_double_format =
4128 false, /* ??? The encoding does allow for "unnormals". */
4129 false, /* ??? The encoding does allow for "unnormals". */
4133 /* The "twos-complement" c4x format is officially defined as
4137 This is rather misleading. One must remember that F is signed.
4138 A better description would be
4140 x = -1**s * ((s + 1 + .f) * 2**e
4142 So if we have a (4 bit) fraction of .1000 with a sign bit of 1,
4143 that's -1 * (1+1+(-.5)) == -1.5. I think.
4145 The constructions here are taken from Tables 5-1 and 5-2 of the
4146 TMS320C4x User's Guide wherein step-by-step instructions for
4147 conversion from IEEE are presented. That's close enough to our
4148 internal representation so as to make things easy.
4150 See http://www-s.ti.com/sc/psheets/spru063c/spru063c.pdf */
4152 static void encode_c4x_single (const struct real_format *fmt,
4153 long *, const REAL_VALUE_TYPE *);
4154 static void decode_c4x_single (const struct real_format *,
4155 REAL_VALUE_TYPE *, const long *);
4156 static void encode_c4x_extended (const struct real_format *fmt,
4157 long *, const REAL_VALUE_TYPE *);
4158 static void decode_c4x_extended (const struct real_format *,
4159 REAL_VALUE_TYPE *, const long *);
4162 encode_c4x_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4163 long *buf, const REAL_VALUE_TYPE *r)
4165 unsigned long image, exp, sig;
4177 sig = 0x800000 - r->sign;
4181 exp = REAL_EXP (r) - 1;
4182 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
4197 image = ((exp & 0xff) << 24) | (sig & 0xffffff);
4202 decode_c4x_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4203 REAL_VALUE_TYPE *r, const long *buf)
4205 unsigned long image = buf[0];
4209 exp = (((image >> 24) & 0xff) ^ 0x80) - 0x80;
4210 sf = ((image & 0xffffff) ^ 0x800000) - 0x800000;
4212 memset (r, 0, sizeof (*r));
4216 r->class = rvc_normal;
4218 sig = sf & 0x7fffff;
4227 sig = (sig << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
4229 SET_REAL_EXP (r, exp + 1);
4230 r->sig[SIGSZ-1] = sig;
4235 encode_c4x_extended (const struct real_format *fmt ATTRIBUTE_UNUSED,
4236 long *buf, const REAL_VALUE_TYPE *r)
4238 unsigned long exp, sig;
4250 sig = 0x80000000 - r->sign;
4254 exp = REAL_EXP (r) - 1;
4256 sig = r->sig[SIGSZ-1];
4257 if (HOST_BITS_PER_LONG == 64)
4258 sig = sig >> 1 >> 31;
4275 exp = (exp & 0xff) << 24;
4278 if (FLOAT_WORDS_BIG_ENDIAN)
4279 buf[0] = exp, buf[1] = sig;
4281 buf[0] = sig, buf[0] = exp;
4285 decode_c4x_extended (const struct real_format *fmt ATTRIBUTE_UNUSED,
4286 REAL_VALUE_TYPE *r, const long *buf)
4291 if (FLOAT_WORDS_BIG_ENDIAN)
4292 exp = buf[0], sf = buf[1];
4294 sf = buf[0], exp = buf[1];
4296 exp = (((exp >> 24) & 0xff) & 0x80) - 0x80;
4297 sf = ((sf & 0xffffffff) ^ 0x80000000) - 0x80000000;
4299 memset (r, 0, sizeof (*r));
4303 r->class = rvc_normal;
4305 sig = sf & 0x7fffffff;
4314 if (HOST_BITS_PER_LONG == 64)
4315 sig = sig << 1 << 31;
4318 SET_REAL_EXP (r, exp + 1);
4319 r->sig[SIGSZ-1] = sig;
4323 const struct real_format c4x_single_format =
4341 const struct real_format c4x_extended_format =
4343 encode_c4x_extended,
4344 decode_c4x_extended,
4360 /* A synthetic "format" for internal arithmetic. It's the size of the
4361 internal significand minus the two bits needed for proper rounding.
4362 The encode and decode routines exist only to satisfy our paranoia
4365 static void encode_internal (const struct real_format *fmt,
4366 long *, const REAL_VALUE_TYPE *);
4367 static void decode_internal (const struct real_format *,
4368 REAL_VALUE_TYPE *, const long *);
4371 encode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4372 const REAL_VALUE_TYPE *r)
4374 memcpy (buf, r, sizeof (*r));
4378 decode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED,
4379 REAL_VALUE_TYPE *r, const long *buf)
4381 memcpy (r, buf, sizeof (*r));
4384 const struct real_format real_internal_format =
4390 SIGNIFICAND_BITS - 2,
4391 SIGNIFICAND_BITS - 2,
4402 /* Calculate the square root of X in mode MODE, and store the result
4403 in R. Return TRUE if the operation does not raise an exception.
4404 For details see "High Precision Division and Square Root",
4405 Alan H. Karp and Peter Markstein, HP Lab Report 93-93-42, June
4406 1993. http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf. */
4409 real_sqrt (REAL_VALUE_TYPE *r, enum machine_mode mode,
4410 const REAL_VALUE_TYPE *x)
4412 static REAL_VALUE_TYPE halfthree;
4413 static bool init = false;
4414 REAL_VALUE_TYPE h, t, i;
4417 /* sqrt(-0.0) is -0.0. */
4418 if (real_isnegzero (x))
4424 /* Negative arguments return NaN. */
4427 get_canonical_qnan (r, 0);
4431 /* Infinity and NaN return themselves. */
4432 if (real_isinf (x) || real_isnan (x))
4440 do_add (&halfthree, &dconst1, &dconsthalf, 0);
4444 /* Initial guess for reciprocal sqrt, i. */
4445 exp = real_exponent (x);
4446 real_ldexp (&i, &dconst1, -exp/2);
4448 /* Newton's iteration for reciprocal sqrt, i. */
4449 for (iter = 0; iter < 16; iter++)
4451 /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x). */
4452 do_multiply (&t, x, &i);
4453 do_multiply (&h, &t, &i);
4454 do_multiply (&t, &h, &dconsthalf);
4455 do_add (&h, &halfthree, &t, 1);
4456 do_multiply (&t, &i, &h);
4458 /* Check for early convergence. */
4459 if (iter >= 6 && real_identical (&i, &t))
4462 /* ??? Unroll loop to avoid copying. */
4466 /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)). */
4467 do_multiply (&t, x, &i);
4468 do_multiply (&h, &t, &i);
4469 do_add (&i, &dconst1, &h, 1);
4470 do_multiply (&h, &t, &i);
4471 do_multiply (&i, &dconsthalf, &h);
4472 do_add (&h, &t, &i, 0);
4474 /* ??? We need a Tuckerman test to get the last bit. */
4476 real_convert (r, mode, &h);
4480 /* Calculate X raised to the integer exponent N in mode MODE and store
4481 the result in R. Return true if the result may be inexact due to
4482 loss of precision. The algorithm is the classic "left-to-right binary
4483 method" described in section 4.6.3 of Donald Knuth's "Seminumerical
4484 Algorithms", "The Art of Computer Programming", Volume 2. */
4487 real_powi (REAL_VALUE_TYPE *r, enum machine_mode mode,
4488 const REAL_VALUE_TYPE *x, HOST_WIDE_INT n)
4490 unsigned HOST_WIDE_INT bit;
4492 bool inexact = false;
4504 /* Don't worry about overflow, from now on n is unsigned. */
4512 bit = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
4513 for (i = 0; i < HOST_BITS_PER_WIDE_INT; i++)
4517 inexact |= do_multiply (&t, &t, &t);
4519 inexact |= do_multiply (&t, &t, x);
4527 inexact |= do_divide (&t, &dconst1, &t);
4529 real_convert (r, mode, &t);
4533 /* Round X to the nearest integer not larger in absolute value, i.e.
4534 towards zero, placing the result in R in mode MODE. */
4537 real_trunc (REAL_VALUE_TYPE *r, enum machine_mode mode,
4538 const REAL_VALUE_TYPE *x)
4540 do_fix_trunc (r, x);
4541 if (mode != VOIDmode)
4542 real_convert (r, mode, r);
4545 /* Round X to the largest integer not greater in value, i.e. round
4546 down, placing the result in R in mode MODE. */
4549 real_floor (REAL_VALUE_TYPE *r, enum machine_mode mode,
4550 const REAL_VALUE_TYPE *x)
4554 do_fix_trunc (&t, x);
4555 if (! real_identical (&t, x) && x->sign)
4556 do_add (&t, &t, &dconstm1, 0);
4557 if (mode != VOIDmode)
4558 real_convert (r, mode, &t);
4561 /* Round X to the smallest integer not less then argument, i.e. round
4562 up, placing the result in R in mode MODE. */
4565 real_ceil (REAL_VALUE_TYPE *r, enum machine_mode mode,
4566 const REAL_VALUE_TYPE *x)
4570 do_fix_trunc (&t, x);
4571 if (! real_identical (&t, x) && ! x->sign)
4572 do_add (&t, &t, &dconst1, 0);
4573 if (mode != VOIDmode)
4574 real_convert (r, mode, &t);
4577 /* Round X to the nearest integer, but round halfway cases away from
4581 real_round (REAL_VALUE_TYPE *r, enum machine_mode mode,
4582 const REAL_VALUE_TYPE *x)
4584 do_add (r, x, &dconsthalf, x->sign);
4585 do_fix_trunc (r, r);
4586 if (mode != VOIDmode)
4587 real_convert (r, mode, r);
4590 /* Set the sign of R to the sign of X. */
4593 real_copysign (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *x)