1 /* real.c - software floating point emulation.
2 Copyright (C) 1993-2013 Free Software Foundation, Inc.
3 Contributed by Stephen L. Moshier (moshier@world.std.com).
4 Re-written by Richard Henderson <rth@redhat.com>
6 This file is part of GCC.
8 GCC is free software; you can redistribute it and/or modify it under
9 the terms of the GNU General Public License as published by the Free
10 Software Foundation; either version 3, or (at your option) any later
13 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
14 WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 You should have received a copy of the GNU General Public License
19 along with GCC; see the file COPYING3. If not see
20 <http://www.gnu.org/licenses/>. */
24 #include "coretypes.h"
27 #include "diagnostic-core.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 26.
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. */
70 /* Used to classify two numbers simultaneously. */
71 #define CLASS2(A, B) ((A) << 2 | (B))
73 #if HOST_BITS_PER_LONG != 64 && HOST_BITS_PER_LONG != 32
74 #error "Some constant folding done by hand to avoid shift count warnings"
77 static void get_zero (REAL_VALUE_TYPE *, int);
78 static void get_canonical_qnan (REAL_VALUE_TYPE *, int);
79 static void get_canonical_snan (REAL_VALUE_TYPE *, int);
80 static void get_inf (REAL_VALUE_TYPE *, int);
81 static bool sticky_rshift_significand (REAL_VALUE_TYPE *,
82 const REAL_VALUE_TYPE *, unsigned int);
83 static void rshift_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
85 static void lshift_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
87 static void lshift_significand_1 (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
88 static bool add_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *,
89 const REAL_VALUE_TYPE *);
90 static bool sub_significands (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
91 const REAL_VALUE_TYPE *, int);
92 static void neg_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
93 static int cmp_significands (const REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
94 static int cmp_significand_0 (const REAL_VALUE_TYPE *);
95 static void set_significand_bit (REAL_VALUE_TYPE *, unsigned int);
96 static void clear_significand_bit (REAL_VALUE_TYPE *, unsigned int);
97 static bool test_significand_bit (REAL_VALUE_TYPE *, unsigned int);
98 static void clear_significand_below (REAL_VALUE_TYPE *, unsigned int);
99 static bool div_significands (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
100 const REAL_VALUE_TYPE *);
101 static void normalize (REAL_VALUE_TYPE *);
103 static bool do_add (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
104 const REAL_VALUE_TYPE *, int);
105 static bool do_multiply (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
106 const REAL_VALUE_TYPE *);
107 static bool do_divide (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
108 const REAL_VALUE_TYPE *);
109 static int do_compare (const REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *, int);
110 static void do_fix_trunc (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
112 static unsigned long rtd_divmod (REAL_VALUE_TYPE *, REAL_VALUE_TYPE *);
113 static void decimal_from_integer (REAL_VALUE_TYPE *);
114 static void decimal_integer_string (char *, const REAL_VALUE_TYPE *,
117 static const REAL_VALUE_TYPE * ten_to_ptwo (int);
118 static const REAL_VALUE_TYPE * ten_to_mptwo (int);
119 static const REAL_VALUE_TYPE * real_digit (int);
120 static void times_pten (REAL_VALUE_TYPE *, int);
122 static void round_for_format (const struct real_format *, REAL_VALUE_TYPE *);
124 /* Initialize R with a positive zero. */
127 get_zero (REAL_VALUE_TYPE *r, int sign)
129 memset (r, 0, sizeof (*r));
133 /* Initialize R with the canonical quiet NaN. */
136 get_canonical_qnan (REAL_VALUE_TYPE *r, int sign)
138 memset (r, 0, sizeof (*r));
145 get_canonical_snan (REAL_VALUE_TYPE *r, int sign)
147 memset (r, 0, sizeof (*r));
155 get_inf (REAL_VALUE_TYPE *r, int sign)
157 memset (r, 0, sizeof (*r));
163 /* Right-shift the significand of A by N bits; put the result in the
164 significand of R. If any one bits are shifted out, return true. */
167 sticky_rshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
170 unsigned long sticky = 0;
171 unsigned int i, ofs = 0;
173 if (n >= HOST_BITS_PER_LONG)
175 for (i = 0, ofs = n / HOST_BITS_PER_LONG; i < ofs; ++i)
177 n &= HOST_BITS_PER_LONG - 1;
182 sticky |= a->sig[ofs] & (((unsigned long)1 << n) - 1);
183 for (i = 0; i < SIGSZ; ++i)
186 = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
187 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
188 << (HOST_BITS_PER_LONG - n)));
193 for (i = 0; ofs + i < SIGSZ; ++i)
194 r->sig[i] = a->sig[ofs + i];
195 for (; i < SIGSZ; ++i)
202 /* Right-shift the significand of A by N bits; put the result in the
206 rshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
209 unsigned int i, ofs = n / HOST_BITS_PER_LONG;
211 n &= HOST_BITS_PER_LONG - 1;
214 for (i = 0; i < SIGSZ; ++i)
217 = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
218 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
219 << (HOST_BITS_PER_LONG - n)));
224 for (i = 0; ofs + i < SIGSZ; ++i)
225 r->sig[i] = a->sig[ofs + i];
226 for (; i < SIGSZ; ++i)
231 /* Left-shift the significand of A by N bits; put the result in the
235 lshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
238 unsigned int i, ofs = n / HOST_BITS_PER_LONG;
240 n &= HOST_BITS_PER_LONG - 1;
243 for (i = 0; ofs + i < SIGSZ; ++i)
244 r->sig[SIGSZ-1-i] = a->sig[SIGSZ-1-i-ofs];
245 for (; i < SIGSZ; ++i)
246 r->sig[SIGSZ-1-i] = 0;
249 for (i = 0; i < SIGSZ; ++i)
252 = (((ofs + i >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs]) << n)
253 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs-1])
254 >> (HOST_BITS_PER_LONG - n)));
258 /* Likewise, but N is specialized to 1. */
261 lshift_significand_1 (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
265 for (i = SIGSZ - 1; i > 0; --i)
266 r->sig[i] = (a->sig[i] << 1) | (a->sig[i-1] >> (HOST_BITS_PER_LONG - 1));
267 r->sig[0] = a->sig[0] << 1;
270 /* Add the significands of A and B, placing the result in R. Return
271 true if there was carry out of the most significant word. */
274 add_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
275 const REAL_VALUE_TYPE *b)
280 for (i = 0; i < SIGSZ; ++i)
282 unsigned long ai = a->sig[i];
283 unsigned long ri = ai + b->sig[i];
299 /* Subtract the significands of A and B, placing the result in R. CARRY is
300 true if there's a borrow incoming to the least significant word.
301 Return true if there was borrow out of the most significant word. */
304 sub_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
305 const REAL_VALUE_TYPE *b, int carry)
309 for (i = 0; i < SIGSZ; ++i)
311 unsigned long ai = a->sig[i];
312 unsigned long ri = ai - b->sig[i];
328 /* Negate the significand A, placing the result in R. */
331 neg_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
336 for (i = 0; i < SIGSZ; ++i)
338 unsigned long ri, ai = a->sig[i];
357 /* Compare significands. Return tri-state vs zero. */
360 cmp_significands (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
364 for (i = SIGSZ - 1; i >= 0; --i)
366 unsigned long ai = a->sig[i];
367 unsigned long bi = b->sig[i];
378 /* Return true if A is nonzero. */
381 cmp_significand_0 (const REAL_VALUE_TYPE *a)
385 for (i = SIGSZ - 1; i >= 0; --i)
392 /* Set bit N of the significand of R. */
395 set_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
397 r->sig[n / HOST_BITS_PER_LONG]
398 |= (unsigned long)1 << (n % HOST_BITS_PER_LONG);
401 /* Clear bit N of the significand of R. */
404 clear_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
406 r->sig[n / HOST_BITS_PER_LONG]
407 &= ~((unsigned long)1 << (n % HOST_BITS_PER_LONG));
410 /* Test bit N of the significand of R. */
413 test_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
415 /* ??? Compiler bug here if we return this expression directly.
416 The conversion to bool strips the "&1" and we wind up testing
417 e.g. 2 != 0 -> true. Seen in gcc version 3.2 20020520. */
418 int t = (r->sig[n / HOST_BITS_PER_LONG] >> (n % HOST_BITS_PER_LONG)) & 1;
422 /* Clear bits 0..N-1 of the significand of R. */
425 clear_significand_below (REAL_VALUE_TYPE *r, unsigned int n)
427 int i, w = n / HOST_BITS_PER_LONG;
429 for (i = 0; i < w; ++i)
432 r->sig[w] &= ~(((unsigned long)1 << (n % HOST_BITS_PER_LONG)) - 1);
435 /* Divide the significands of A and B, placing the result in R. Return
436 true if the division was inexact. */
439 div_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
440 const REAL_VALUE_TYPE *b)
443 int i, bit = SIGNIFICAND_BITS - 1;
444 unsigned long msb, inexact;
447 memset (r->sig, 0, sizeof (r->sig));
453 msb = u.sig[SIGSZ-1] & SIG_MSB;
454 lshift_significand_1 (&u, &u);
456 if (msb || cmp_significands (&u, b) >= 0)
458 sub_significands (&u, &u, b, 0);
459 set_significand_bit (r, bit);
464 for (i = 0, inexact = 0; i < SIGSZ; i++)
470 /* Adjust the exponent and significand of R such that the most
471 significant bit is set. We underflow to zero and overflow to
472 infinity here, without denormals. (The intermediate representation
473 exponent is large enough to handle target denormals normalized.) */
476 normalize (REAL_VALUE_TYPE *r)
484 /* Find the first word that is nonzero. */
485 for (i = SIGSZ - 1; i >= 0; i--)
487 shift += HOST_BITS_PER_LONG;
491 /* Zero significand flushes to zero. */
499 /* Find the first bit that is nonzero. */
501 if (r->sig[i] & ((unsigned long)1 << (HOST_BITS_PER_LONG - 1 - j)))
507 exp = REAL_EXP (r) - shift;
509 get_inf (r, r->sign);
510 else if (exp < -MAX_EXP)
511 get_zero (r, r->sign);
514 SET_REAL_EXP (r, exp);
515 lshift_significand (r, r, shift);
520 /* Calculate R = A + (SUBTRACT_P ? -B : B). Return true if the
521 result may be inexact due to a loss of precision. */
524 do_add (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
525 const REAL_VALUE_TYPE *b, int subtract_p)
529 bool inexact = false;
531 /* Determine if we need to add or subtract. */
533 subtract_p = (sign ^ b->sign) ^ subtract_p;
535 switch (CLASS2 (a->cl, b->cl))
537 case CLASS2 (rvc_zero, rvc_zero):
538 /* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0. */
539 get_zero (r, sign & !subtract_p);
542 case CLASS2 (rvc_zero, rvc_normal):
543 case CLASS2 (rvc_zero, rvc_inf):
544 case CLASS2 (rvc_zero, rvc_nan):
546 case CLASS2 (rvc_normal, rvc_nan):
547 case CLASS2 (rvc_inf, rvc_nan):
548 case CLASS2 (rvc_nan, rvc_nan):
549 /* ANY + NaN = NaN. */
550 case CLASS2 (rvc_normal, rvc_inf):
553 r->sign = sign ^ subtract_p;
556 case CLASS2 (rvc_normal, rvc_zero):
557 case CLASS2 (rvc_inf, rvc_zero):
558 case CLASS2 (rvc_nan, rvc_zero):
560 case CLASS2 (rvc_nan, rvc_normal):
561 case CLASS2 (rvc_nan, rvc_inf):
562 /* NaN + ANY = NaN. */
563 case CLASS2 (rvc_inf, rvc_normal):
568 case CLASS2 (rvc_inf, rvc_inf):
570 /* Inf - Inf = NaN. */
571 get_canonical_qnan (r, 0);
573 /* Inf + Inf = Inf. */
577 case CLASS2 (rvc_normal, rvc_normal):
584 /* Swap the arguments such that A has the larger exponent. */
585 dexp = REAL_EXP (a) - REAL_EXP (b);
588 const REAL_VALUE_TYPE *t;
595 /* If the exponents are not identical, we need to shift the
596 significand of B down. */
599 /* If the exponents are too far apart, the significands
600 do not overlap, which makes the subtraction a noop. */
601 if (dexp >= SIGNIFICAND_BITS)
608 inexact |= sticky_rshift_significand (&t, b, dexp);
614 if (sub_significands (r, a, b, inexact))
616 /* We got a borrow out of the subtraction. That means that
617 A and B had the same exponent, and B had the larger
618 significand. We need to swap the sign and negate the
621 neg_significand (r, r);
626 if (add_significands (r, a, b))
628 /* We got carry out of the addition. This means we need to
629 shift the significand back down one bit and increase the
631 inexact |= sticky_rshift_significand (r, r, 1);
632 r->sig[SIGSZ-1] |= SIG_MSB;
643 SET_REAL_EXP (r, exp);
644 /* Zero out the remaining fields. */
649 /* Re-normalize the result. */
652 /* Special case: if the subtraction results in zero, the result
654 if (r->cl == rvc_zero)
657 r->sig[0] |= inexact;
662 /* Calculate R = A * B. Return true if the result may be inexact. */
665 do_multiply (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
666 const REAL_VALUE_TYPE *b)
668 REAL_VALUE_TYPE u, t, *rr;
669 unsigned int i, j, k;
670 int sign = a->sign ^ b->sign;
671 bool inexact = false;
673 switch (CLASS2 (a->cl, b->cl))
675 case CLASS2 (rvc_zero, rvc_zero):
676 case CLASS2 (rvc_zero, rvc_normal):
677 case CLASS2 (rvc_normal, rvc_zero):
678 /* +-0 * ANY = 0 with appropriate sign. */
682 case CLASS2 (rvc_zero, rvc_nan):
683 case CLASS2 (rvc_normal, rvc_nan):
684 case CLASS2 (rvc_inf, rvc_nan):
685 case CLASS2 (rvc_nan, rvc_nan):
686 /* ANY * NaN = NaN. */
691 case CLASS2 (rvc_nan, rvc_zero):
692 case CLASS2 (rvc_nan, rvc_normal):
693 case CLASS2 (rvc_nan, rvc_inf):
694 /* NaN * ANY = NaN. */
699 case CLASS2 (rvc_zero, rvc_inf):
700 case CLASS2 (rvc_inf, rvc_zero):
702 get_canonical_qnan (r, sign);
705 case CLASS2 (rvc_inf, rvc_inf):
706 case CLASS2 (rvc_normal, rvc_inf):
707 case CLASS2 (rvc_inf, rvc_normal):
708 /* Inf * Inf = Inf, R * Inf = Inf */
712 case CLASS2 (rvc_normal, rvc_normal):
719 if (r == a || r == b)
725 /* Collect all the partial products. Since we don't have sure access
726 to a widening multiply, we split each long into two half-words.
728 Consider the long-hand form of a four half-word multiplication:
738 We construct partial products of the widened half-word products
739 that are known to not overlap, e.g. DF+DH. Each such partial
740 product is given its proper exponent, which allows us to sum them
741 and obtain the finished product. */
743 for (i = 0; i < SIGSZ * 2; ++i)
745 unsigned long ai = a->sig[i / 2];
747 ai >>= HOST_BITS_PER_LONG / 2;
749 ai &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
754 for (j = 0; j < 2; ++j)
756 int exp = (REAL_EXP (a) - (2*SIGSZ-1-i)*(HOST_BITS_PER_LONG/2)
757 + (REAL_EXP (b) - (1-j)*(HOST_BITS_PER_LONG/2)));
766 /* Would underflow to zero, which we shouldn't bother adding. */
771 memset (&u, 0, sizeof (u));
773 SET_REAL_EXP (&u, exp);
775 for (k = j; k < SIGSZ * 2; k += 2)
777 unsigned long bi = b->sig[k / 2];
779 bi >>= HOST_BITS_PER_LONG / 2;
781 bi &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
783 u.sig[k / 2] = ai * bi;
787 inexact |= do_add (rr, rr, &u, 0);
798 /* Calculate R = A / B. Return true if the result may be inexact. */
801 do_divide (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
802 const REAL_VALUE_TYPE *b)
804 int exp, sign = a->sign ^ b->sign;
805 REAL_VALUE_TYPE t, *rr;
808 switch (CLASS2 (a->cl, b->cl))
810 case CLASS2 (rvc_zero, rvc_zero):
812 case CLASS2 (rvc_inf, rvc_inf):
813 /* Inf / Inf = NaN. */
814 get_canonical_qnan (r, sign);
817 case CLASS2 (rvc_zero, rvc_normal):
818 case CLASS2 (rvc_zero, rvc_inf):
820 case CLASS2 (rvc_normal, rvc_inf):
825 case CLASS2 (rvc_normal, rvc_zero):
827 case CLASS2 (rvc_inf, rvc_zero):
832 case CLASS2 (rvc_zero, rvc_nan):
833 case CLASS2 (rvc_normal, rvc_nan):
834 case CLASS2 (rvc_inf, rvc_nan):
835 case CLASS2 (rvc_nan, rvc_nan):
836 /* ANY / NaN = NaN. */
841 case CLASS2 (rvc_nan, rvc_zero):
842 case CLASS2 (rvc_nan, rvc_normal):
843 case CLASS2 (rvc_nan, rvc_inf):
844 /* NaN / ANY = NaN. */
849 case CLASS2 (rvc_inf, rvc_normal):
854 case CLASS2 (rvc_normal, rvc_normal):
861 if (r == a || r == b)
866 /* Make sure all fields in the result are initialized. */
871 exp = REAL_EXP (a) - REAL_EXP (b) + 1;
882 SET_REAL_EXP (rr, exp);
884 inexact = div_significands (rr, a, b);
886 /* Re-normalize the result. */
888 rr->sig[0] |= inexact;
896 /* Return a tri-state comparison of A vs B. Return NAN_RESULT if
897 one of the two operands is a NaN. */
900 do_compare (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b,
905 switch (CLASS2 (a->cl, b->cl))
907 case CLASS2 (rvc_zero, rvc_zero):
908 /* Sign of zero doesn't matter for compares. */
911 case CLASS2 (rvc_normal, rvc_zero):
912 /* Decimal float zero is special and uses rvc_normal, not rvc_zero. */
914 return decimal_do_compare (a, b, nan_result);
916 case CLASS2 (rvc_inf, rvc_zero):
917 case CLASS2 (rvc_inf, rvc_normal):
918 return (a->sign ? -1 : 1);
920 case CLASS2 (rvc_inf, rvc_inf):
921 return -a->sign - -b->sign;
923 case CLASS2 (rvc_zero, rvc_normal):
924 /* Decimal float zero is special and uses rvc_normal, not rvc_zero. */
926 return decimal_do_compare (a, b, nan_result);
928 case CLASS2 (rvc_zero, rvc_inf):
929 case CLASS2 (rvc_normal, rvc_inf):
930 return (b->sign ? 1 : -1);
932 case CLASS2 (rvc_zero, rvc_nan):
933 case CLASS2 (rvc_normal, rvc_nan):
934 case CLASS2 (rvc_inf, rvc_nan):
935 case CLASS2 (rvc_nan, rvc_nan):
936 case CLASS2 (rvc_nan, rvc_zero):
937 case CLASS2 (rvc_nan, rvc_normal):
938 case CLASS2 (rvc_nan, rvc_inf):
941 case CLASS2 (rvc_normal, rvc_normal):
948 if (a->sign != b->sign)
949 return -a->sign - -b->sign;
951 if (a->decimal || b->decimal)
952 return decimal_do_compare (a, b, nan_result);
954 if (REAL_EXP (a) > REAL_EXP (b))
956 else if (REAL_EXP (a) < REAL_EXP (b))
959 ret = cmp_significands (a, b);
961 return (a->sign ? -ret : ret);
964 /* Return A truncated to an integral value toward zero. */
967 do_fix_trunc (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
981 decimal_do_fix_trunc (r, a);
984 if (REAL_EXP (r) <= 0)
985 get_zero (r, r->sign);
986 else if (REAL_EXP (r) < SIGNIFICAND_BITS)
987 clear_significand_below (r, SIGNIFICAND_BITS - REAL_EXP (r));
995 /* Perform the binary or unary operation described by CODE.
996 For a unary operation, leave OP1 NULL. This function returns
997 true if the result may be inexact due to loss of precision. */
1000 real_arithmetic (REAL_VALUE_TYPE *r, int icode, const REAL_VALUE_TYPE *op0,
1001 const REAL_VALUE_TYPE *op1)
1003 enum tree_code code = (enum tree_code) icode;
1005 if (op0->decimal || (op1 && op1->decimal))
1006 return decimal_real_arithmetic (r, code, op0, op1);
1011 /* Clear any padding areas in *r if it isn't equal to one of the
1012 operands so that we can later do bitwise comparisons later on. */
1013 if (r != op0 && r != op1)
1014 memset (r, '\0', sizeof (*r));
1015 return do_add (r, op0, op1, 0);
1018 if (r != op0 && r != op1)
1019 memset (r, '\0', sizeof (*r));
1020 return do_add (r, op0, op1, 1);
1023 if (r != op0 && r != op1)
1024 memset (r, '\0', sizeof (*r));
1025 return do_multiply (r, op0, op1);
1028 if (r != op0 && r != op1)
1029 memset (r, '\0', sizeof (*r));
1030 return do_divide (r, op0, op1);
1033 if (op1->cl == rvc_nan)
1035 else if (do_compare (op0, op1, -1) < 0)
1042 if (op1->cl == rvc_nan)
1044 else if (do_compare (op0, op1, 1) < 0)
1060 case FIX_TRUNC_EXPR:
1061 do_fix_trunc (r, op0);
1071 real_value_negate (const REAL_VALUE_TYPE *op0)
1074 real_arithmetic (&r, NEGATE_EXPR, op0, NULL);
1079 real_value_abs (const REAL_VALUE_TYPE *op0)
1082 real_arithmetic (&r, ABS_EXPR, op0, NULL);
1087 real_compare (int icode, const REAL_VALUE_TYPE *op0,
1088 const REAL_VALUE_TYPE *op1)
1090 enum tree_code code = (enum tree_code) icode;
1095 return do_compare (op0, op1, 1) < 0;
1097 return do_compare (op0, op1, 1) <= 0;
1099 return do_compare (op0, op1, -1) > 0;
1101 return do_compare (op0, op1, -1) >= 0;
1103 return do_compare (op0, op1, -1) == 0;
1105 return do_compare (op0, op1, -1) != 0;
1106 case UNORDERED_EXPR:
1107 return op0->cl == rvc_nan || op1->cl == rvc_nan;
1109 return op0->cl != rvc_nan && op1->cl != rvc_nan;
1111 return do_compare (op0, op1, -1) < 0;
1113 return do_compare (op0, op1, -1) <= 0;
1115 return do_compare (op0, op1, 1) > 0;
1117 return do_compare (op0, op1, 1) >= 0;
1119 return do_compare (op0, op1, 0) == 0;
1121 return do_compare (op0, op1, 0) != 0;
1128 /* Return floor log2(R). */
1131 real_exponent (const REAL_VALUE_TYPE *r)
1139 return (unsigned int)-1 >> 1;
1141 return REAL_EXP (r);
1147 /* R = OP0 * 2**EXP. */
1150 real_ldexp (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *op0, int exp)
1161 exp += REAL_EXP (op0);
1163 get_inf (r, r->sign);
1164 else if (exp < -MAX_EXP)
1165 get_zero (r, r->sign);
1167 SET_REAL_EXP (r, exp);
1175 /* Determine whether a floating-point value X is infinite. */
1178 real_isinf (const REAL_VALUE_TYPE *r)
1180 return (r->cl == rvc_inf);
1183 /* Determine whether a floating-point value X is a NaN. */
1186 real_isnan (const REAL_VALUE_TYPE *r)
1188 return (r->cl == rvc_nan);
1191 /* Determine whether a floating-point value X is finite. */
1194 real_isfinite (const REAL_VALUE_TYPE *r)
1196 return (r->cl != rvc_nan) && (r->cl != rvc_inf);
1199 /* Determine whether a floating-point value X is negative. */
1202 real_isneg (const REAL_VALUE_TYPE *r)
1207 /* Determine whether a floating-point value X is minus zero. */
1210 real_isnegzero (const REAL_VALUE_TYPE *r)
1212 return r->sign && r->cl == rvc_zero;
1215 /* Compare two floating-point objects for bitwise identity. */
1218 real_identical (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
1224 if (a->sign != b->sign)
1234 if (a->decimal != b->decimal)
1236 if (REAL_EXP (a) != REAL_EXP (b))
1241 if (a->signalling != b->signalling)
1243 /* The significand is ignored for canonical NaNs. */
1244 if (a->canonical || b->canonical)
1245 return a->canonical == b->canonical;
1252 for (i = 0; i < SIGSZ; ++i)
1253 if (a->sig[i] != b->sig[i])
1259 /* Try to change R into its exact multiplicative inverse in machine
1260 mode MODE. Return true if successful. */
1263 exact_real_inverse (enum machine_mode mode, REAL_VALUE_TYPE *r)
1265 const REAL_VALUE_TYPE *one = real_digit (1);
1269 if (r->cl != rvc_normal)
1272 /* Check for a power of two: all significand bits zero except the MSB. */
1273 for (i = 0; i < SIGSZ-1; ++i)
1276 if (r->sig[SIGSZ-1] != SIG_MSB)
1279 /* Find the inverse and truncate to the required mode. */
1280 do_divide (&u, one, r);
1281 real_convert (&u, mode, &u);
1283 /* The rounding may have overflowed. */
1284 if (u.cl != rvc_normal)
1286 for (i = 0; i < SIGSZ-1; ++i)
1289 if (u.sig[SIGSZ-1] != SIG_MSB)
1296 /* Return true if arithmetic on values in IMODE that were promoted
1297 from values in TMODE is equivalent to direct arithmetic on values
1301 real_can_shorten_arithmetic (enum machine_mode imode, enum machine_mode tmode)
1303 const struct real_format *tfmt, *ifmt;
1304 tfmt = REAL_MODE_FORMAT (tmode);
1305 ifmt = REAL_MODE_FORMAT (imode);
1306 /* These conditions are conservative rather than trying to catch the
1307 exact boundary conditions; the main case to allow is IEEE float
1309 return (ifmt->b == tfmt->b
1310 && ifmt->p > 2 * tfmt->p
1311 && ifmt->emin < 2 * tfmt->emin - tfmt->p - 2
1312 && ifmt->emin < tfmt->emin - tfmt->emax - tfmt->p - 2
1313 && ifmt->emax > 2 * tfmt->emax + 2
1314 && ifmt->emax > tfmt->emax - tfmt->emin + tfmt->p + 2
1315 && ifmt->round_towards_zero == tfmt->round_towards_zero
1316 && (ifmt->has_sign_dependent_rounding
1317 == tfmt->has_sign_dependent_rounding)
1318 && ifmt->has_nans >= tfmt->has_nans
1319 && ifmt->has_inf >= tfmt->has_inf
1320 && ifmt->has_signed_zero >= tfmt->has_signed_zero
1321 && !MODE_COMPOSITE_P (tmode)
1322 && !MODE_COMPOSITE_P (imode));
1325 /* Render R as an integer. */
1328 real_to_integer (const REAL_VALUE_TYPE *r)
1330 unsigned HOST_WIDE_INT i;
1341 i = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1348 return decimal_real_to_integer (r);
1350 if (REAL_EXP (r) <= 0)
1352 /* Only force overflow for unsigned overflow. Signed overflow is
1353 undefined, so it doesn't matter what we return, and some callers
1354 expect to be able to use this routine for both signed and
1355 unsigned conversions. */
1356 if (REAL_EXP (r) > HOST_BITS_PER_WIDE_INT)
1359 if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1360 i = r->sig[SIGSZ-1];
1363 gcc_assert (HOST_BITS_PER_WIDE_INT == 2 * HOST_BITS_PER_LONG);
1364 i = r->sig[SIGSZ-1];
1365 i = i << (HOST_BITS_PER_LONG - 1) << 1;
1366 i |= r->sig[SIGSZ-2];
1369 i >>= HOST_BITS_PER_WIDE_INT - REAL_EXP (r);
1380 /* Likewise, but to an integer pair, HI+LOW. */
1383 real_to_integer2 (HOST_WIDE_INT *plow, HOST_WIDE_INT *phigh,
1384 const REAL_VALUE_TYPE *r)
1387 HOST_WIDE_INT low, high;
1400 high = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1413 decimal_real_to_integer2 (plow, phigh, r);
1420 /* Only force overflow for unsigned overflow. Signed overflow is
1421 undefined, so it doesn't matter what we return, and some callers
1422 expect to be able to use this routine for both signed and
1423 unsigned conversions. */
1424 if (exp > HOST_BITS_PER_DOUBLE_INT)
1427 rshift_significand (&t, r, HOST_BITS_PER_DOUBLE_INT - exp);
1428 if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1430 high = t.sig[SIGSZ-1];
1431 low = t.sig[SIGSZ-2];
1435 gcc_assert (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG);
1436 high = t.sig[SIGSZ-1];
1437 high = high << (HOST_BITS_PER_LONG - 1) << 1;
1438 high |= t.sig[SIGSZ-2];
1440 low = t.sig[SIGSZ-3];
1441 low = low << (HOST_BITS_PER_LONG - 1) << 1;
1442 low |= t.sig[SIGSZ-4];
1450 low = -low, high = ~high;
1462 /* A subroutine of real_to_decimal. Compute the quotient and remainder
1463 of NUM / DEN. Return the quotient and place the remainder in NUM.
1464 It is expected that NUM / DEN are close enough that the quotient is
1467 static unsigned long
1468 rtd_divmod (REAL_VALUE_TYPE *num, REAL_VALUE_TYPE *den)
1470 unsigned long q, msb;
1471 int expn = REAL_EXP (num), expd = REAL_EXP (den);
1480 msb = num->sig[SIGSZ-1] & SIG_MSB;
1482 lshift_significand_1 (num, num);
1484 if (msb || cmp_significands (num, den) >= 0)
1486 sub_significands (num, num, den, 0);
1490 while (--expn >= expd);
1492 SET_REAL_EXP (num, expd);
1498 /* Render R as a decimal floating point constant. Emit DIGITS significant
1499 digits in the result, bounded by BUF_SIZE. If DIGITS is 0, choose the
1500 maximum for the representation. If CROP_TRAILING_ZEROS, strip trailing
1501 zeros. If MODE is VOIDmode, round to nearest value. Otherwise, round
1502 to a string that, when parsed back in mode MODE, yields the same value. */
1504 #define M_LOG10_2 0.30102999566398119521
1507 real_to_decimal_for_mode (char *str, const REAL_VALUE_TYPE *r_orig,
1508 size_t buf_size, size_t digits,
1509 int crop_trailing_zeros, enum machine_mode mode)
1511 const struct real_format *fmt = NULL;
1512 const REAL_VALUE_TYPE *one, *ten;
1513 REAL_VALUE_TYPE r, pten, u, v;
1514 int dec_exp, cmp_one, digit;
1516 char *p, *first, *last;
1520 if (mode != VOIDmode)
1522 fmt = REAL_MODE_FORMAT (mode);
1530 strcpy (str, (r.sign ? "-0.0" : "0.0"));
1535 strcpy (str, (r.sign ? "-Inf" : "+Inf"));
1538 /* ??? Print the significand as well, if not canonical? */
1539 sprintf (str, "%c%cNaN", (r_orig->sign ? '-' : '+'),
1540 (r_orig->signalling ? 'S' : 'Q'));
1548 decimal_real_to_decimal (str, &r, buf_size, digits, crop_trailing_zeros);
1552 /* Bound the number of digits printed by the size of the representation. */
1553 max_digits = SIGNIFICAND_BITS * M_LOG10_2;
1554 if (digits == 0 || digits > max_digits)
1555 digits = max_digits;
1557 /* Estimate the decimal exponent, and compute the length of the string it
1558 will print as. Be conservative and add one to account for possible
1559 overflow or rounding error. */
1560 dec_exp = REAL_EXP (&r) * M_LOG10_2;
1561 for (max_digits = 1; dec_exp ; max_digits++)
1564 /* Bound the number of digits printed by the size of the output buffer. */
1565 max_digits = buf_size - 1 - 1 - 2 - max_digits - 1;
1566 gcc_assert (max_digits <= buf_size);
1567 if (digits > max_digits)
1568 digits = max_digits;
1570 one = real_digit (1);
1571 ten = ten_to_ptwo (0);
1579 cmp_one = do_compare (&r, one, 0);
1584 /* Number is greater than one. Convert significand to an integer
1585 and strip trailing decimal zeros. */
1588 SET_REAL_EXP (&u, SIGNIFICAND_BITS - 1);
1590 /* Largest M, such that 10**2**M fits within SIGNIFICAND_BITS. */
1591 m = floor_log2 (max_digits);
1593 /* Iterate over the bits of the possible powers of 10 that might
1594 be present in U and eliminate them. That is, if we find that
1595 10**2**M divides U evenly, keep the division and increase
1601 do_divide (&t, &u, ten_to_ptwo (m));
1602 do_fix_trunc (&v, &t);
1603 if (cmp_significands (&v, &t) == 0)
1611 /* Revert the scaling to integer that we performed earlier. */
1612 SET_REAL_EXP (&u, REAL_EXP (&u) + REAL_EXP (&r)
1613 - (SIGNIFICAND_BITS - 1));
1616 /* Find power of 10. Do this by dividing out 10**2**M when
1617 this is larger than the current remainder. Fill PTEN with
1618 the power of 10 that we compute. */
1619 if (REAL_EXP (&r) > 0)
1621 m = floor_log2 ((int)(REAL_EXP (&r) * M_LOG10_2)) + 1;
1624 const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1625 if (do_compare (&u, ptentwo, 0) >= 0)
1627 do_divide (&u, &u, ptentwo);
1628 do_multiply (&pten, &pten, ptentwo);
1635 /* We managed to divide off enough tens in the above reduction
1636 loop that we've now got a negative exponent. Fall into the
1637 less-than-one code to compute the proper value for PTEN. */
1644 /* Number is less than one. Pad significand with leading
1650 /* Stop if we'd shift bits off the bottom. */
1654 do_multiply (&u, &v, ten);
1656 /* Stop if we're now >= 1. */
1657 if (REAL_EXP (&u) > 0)
1665 /* Find power of 10. Do this by multiplying in P=10**2**M when
1666 the current remainder is smaller than 1/P. Fill PTEN with the
1667 power of 10 that we compute. */
1668 m = floor_log2 ((int)(-REAL_EXP (&r) * M_LOG10_2)) + 1;
1671 const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1672 const REAL_VALUE_TYPE *ptenmtwo = ten_to_mptwo (m);
1674 if (do_compare (&v, ptenmtwo, 0) <= 0)
1676 do_multiply (&v, &v, ptentwo);
1677 do_multiply (&pten, &pten, ptentwo);
1683 /* Invert the positive power of 10 that we've collected so far. */
1684 do_divide (&pten, one, &pten);
1692 /* At this point, PTEN should contain the nearest power of 10 smaller
1693 than R, such that this division produces the first digit.
1695 Using a divide-step primitive that returns the complete integral
1696 remainder avoids the rounding error that would be produced if
1697 we were to use do_divide here and then simply multiply by 10 for
1698 each subsequent digit. */
1700 digit = rtd_divmod (&r, &pten);
1702 /* Be prepared for error in that division via underflow ... */
1703 if (digit == 0 && cmp_significand_0 (&r))
1705 /* Multiply by 10 and try again. */
1706 do_multiply (&r, &r, ten);
1707 digit = rtd_divmod (&r, &pten);
1709 gcc_assert (digit != 0);
1712 /* ... or overflow. */
1722 gcc_assert (digit <= 10);
1726 /* Generate subsequent digits. */
1727 while (--digits > 0)
1729 do_multiply (&r, &r, ten);
1730 digit = rtd_divmod (&r, &pten);
1735 /* Generate one more digit with which to do rounding. */
1736 do_multiply (&r, &r, ten);
1737 digit = rtd_divmod (&r, &pten);
1739 /* Round the result. */
1740 if (fmt && fmt->round_towards_zero)
1742 /* If the format uses round towards zero when parsing the string
1743 back in, we need to always round away from zero here. */
1744 if (cmp_significand_0 (&r))
1746 round_up = digit > 0;
1752 /* Round to nearest. If R is nonzero there are additional
1753 nonzero digits to be extracted. */
1754 if (cmp_significand_0 (&r))
1756 /* Round to even. */
1757 else if ((p[-1] - '0') & 1)
1761 round_up = digit > 5;
1778 /* Carry out of the first digit. This means we had all 9's and
1779 now have all 0's. "Prepend" a 1 by overwriting the first 0. */
1787 /* Insert the decimal point. */
1788 first[0] = first[1];
1791 /* If requested, drop trailing zeros. Never crop past "1.0". */
1792 if (crop_trailing_zeros)
1793 while (last > first + 3 && last[-1] == '0')
1796 /* Append the exponent. */
1797 sprintf (last, "e%+d", dec_exp);
1799 #ifdef ENABLE_CHECKING
1800 /* Verify that we can read the original value back in. */
1801 if (mode != VOIDmode)
1803 real_from_string (&r, str);
1804 real_convert (&r, mode, &r);
1805 gcc_assert (real_identical (&r, r_orig));
1810 /* Likewise, except always uses round-to-nearest. */
1813 real_to_decimal (char *str, const REAL_VALUE_TYPE *r_orig, size_t buf_size,
1814 size_t digits, int crop_trailing_zeros)
1816 real_to_decimal_for_mode (str, r_orig, buf_size,
1817 digits, crop_trailing_zeros, VOIDmode);
1820 /* Render R as a hexadecimal floating point constant. Emit DIGITS
1821 significant digits in the result, bounded by BUF_SIZE. If DIGITS is 0,
1822 choose the maximum for the representation. If CROP_TRAILING_ZEROS,
1823 strip trailing zeros. */
1826 real_to_hexadecimal (char *str, const REAL_VALUE_TYPE *r, size_t buf_size,
1827 size_t digits, int crop_trailing_zeros)
1829 int i, j, exp = REAL_EXP (r);
1842 strcpy (str, (r->sign ? "-Inf" : "+Inf"));
1845 /* ??? Print the significand as well, if not canonical? */
1846 sprintf (str, "%c%cNaN", (r->sign ? '-' : '+'),
1847 (r->signalling ? 'S' : 'Q'));
1855 /* Hexadecimal format for decimal floats is not interesting. */
1856 strcpy (str, "N/A");
1861 digits = SIGNIFICAND_BITS / 4;
1863 /* Bound the number of digits printed by the size of the output buffer. */
1865 sprintf (exp_buf, "p%+d", exp);
1866 max_digits = buf_size - strlen (exp_buf) - r->sign - 4 - 1;
1867 gcc_assert (max_digits <= buf_size);
1868 if (digits > max_digits)
1869 digits = max_digits;
1880 for (i = SIGSZ - 1; i >= 0; --i)
1881 for (j = HOST_BITS_PER_LONG - 4; j >= 0; j -= 4)
1883 *p++ = "0123456789abcdef"[(r->sig[i] >> j) & 15];
1889 if (crop_trailing_zeros)
1890 while (p > first + 1 && p[-1] == '0')
1893 sprintf (p, "p%+d", exp);
1896 /* Initialize R from a decimal or hexadecimal string. The string is
1897 assumed to have been syntax checked already. Return -1 if the
1898 value underflows, +1 if overflows, and 0 otherwise. */
1901 real_from_string (REAL_VALUE_TYPE *r, const char *str)
1913 else if (*str == '+')
1916 if (!strncmp (str, "QNaN", 4))
1918 get_canonical_qnan (r, sign);
1921 else if (!strncmp (str, "SNaN", 4))
1923 get_canonical_snan (r, sign);
1926 else if (!strncmp (str, "Inf", 3))
1932 if (str[0] == '0' && (str[1] == 'x' || str[1] == 'X'))
1934 /* Hexadecimal floating point. */
1935 int pos = SIGNIFICAND_BITS - 4, d;
1943 d = hex_value (*str);
1948 r->sig[pos / HOST_BITS_PER_LONG]
1949 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1953 /* Ensure correct rounding by setting last bit if there is
1954 a subsequent nonzero digit. */
1962 if (pos == SIGNIFICAND_BITS - 4)
1969 d = hex_value (*str);
1974 r->sig[pos / HOST_BITS_PER_LONG]
1975 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1979 /* Ensure correct rounding by setting last bit if there is
1980 a subsequent nonzero digit. */
1986 /* If the mantissa is zero, ignore the exponent. */
1987 if (!cmp_significand_0 (r))
1990 if (*str == 'p' || *str == 'P')
1992 bool exp_neg = false;
2000 else if (*str == '+')
2004 while (ISDIGIT (*str))
2010 /* Overflowed the exponent. */
2025 SET_REAL_EXP (r, exp);
2031 /* Decimal floating point. */
2032 const REAL_VALUE_TYPE *ten = ten_to_ptwo (0);
2037 while (ISDIGIT (*str))
2040 do_multiply (r, r, ten);
2042 do_add (r, r, real_digit (d), 0);
2047 if (r->cl == rvc_zero)
2052 while (ISDIGIT (*str))
2055 do_multiply (r, r, ten);
2057 do_add (r, r, real_digit (d), 0);
2062 /* If the mantissa is zero, ignore the exponent. */
2063 if (r->cl == rvc_zero)
2066 if (*str == 'e' || *str == 'E')
2068 bool exp_neg = false;
2076 else if (*str == '+')
2080 while (ISDIGIT (*str))
2086 /* Overflowed the exponent. */
2100 times_pten (r, exp);
2119 /* Legacy. Similar, but return the result directly. */
2122 real_from_string2 (const char *s, enum machine_mode mode)
2126 real_from_string (&r, s);
2127 if (mode != VOIDmode)
2128 real_convert (&r, mode, &r);
2133 /* Initialize R from string S and desired MODE. */
2136 real_from_string3 (REAL_VALUE_TYPE *r, const char *s, enum machine_mode mode)
2138 if (DECIMAL_FLOAT_MODE_P (mode))
2139 decimal_real_from_string (r, s);
2141 real_from_string (r, s);
2143 if (mode != VOIDmode)
2144 real_convert (r, mode, r);
2147 /* Initialize R from the integer pair HIGH+LOW. */
2150 real_from_integer (REAL_VALUE_TYPE *r, enum machine_mode mode,
2151 unsigned HOST_WIDE_INT low, HOST_WIDE_INT high,
2154 if (low == 0 && high == 0)
2158 memset (r, 0, sizeof (*r));
2160 r->sign = high < 0 && !unsigned_p;
2161 SET_REAL_EXP (r, HOST_BITS_PER_DOUBLE_INT);
2172 if (HOST_BITS_PER_LONG == HOST_BITS_PER_WIDE_INT)
2174 r->sig[SIGSZ-1] = high;
2175 r->sig[SIGSZ-2] = low;
2179 gcc_assert (HOST_BITS_PER_LONG*2 == HOST_BITS_PER_WIDE_INT);
2180 r->sig[SIGSZ-1] = high >> (HOST_BITS_PER_LONG - 1) >> 1;
2181 r->sig[SIGSZ-2] = high;
2182 r->sig[SIGSZ-3] = low >> (HOST_BITS_PER_LONG - 1) >> 1;
2183 r->sig[SIGSZ-4] = low;
2189 if (DECIMAL_FLOAT_MODE_P (mode))
2190 decimal_from_integer (r);
2191 else if (mode != VOIDmode)
2192 real_convert (r, mode, r);
2195 /* Render R, an integral value, as a floating point constant with no
2196 specified exponent. */
2199 decimal_integer_string (char *str, const REAL_VALUE_TYPE *r_orig,
2202 int dec_exp, digit, digits;
2203 REAL_VALUE_TYPE r, pten;
2209 if (r.cl == rvc_zero)
2218 dec_exp = REAL_EXP (&r) * M_LOG10_2;
2219 digits = dec_exp + 1;
2220 gcc_assert ((digits + 2) < (int)buf_size);
2222 pten = *real_digit (1);
2223 times_pten (&pten, dec_exp);
2229 digit = rtd_divmod (&r, &pten);
2230 gcc_assert (digit >= 0 && digit <= 9);
2232 while (--digits > 0)
2235 digit = rtd_divmod (&r, &pten);
2242 /* Convert a real with an integral value to decimal float. */
2245 decimal_from_integer (REAL_VALUE_TYPE *r)
2249 decimal_integer_string (str, r, sizeof (str) - 1);
2250 decimal_real_from_string (r, str);
2253 /* Returns 10**2**N. */
2255 static const REAL_VALUE_TYPE *
2258 static REAL_VALUE_TYPE tens[EXP_BITS];
2260 gcc_assert (n >= 0);
2261 gcc_assert (n < EXP_BITS);
2263 if (tens[n].cl == rvc_zero)
2265 if (n < (HOST_BITS_PER_WIDE_INT == 64 ? 5 : 4))
2267 HOST_WIDE_INT t = 10;
2270 for (i = 0; i < n; ++i)
2273 real_from_integer (&tens[n], VOIDmode, t, 0, 1);
2277 const REAL_VALUE_TYPE *t = ten_to_ptwo (n - 1);
2278 do_multiply (&tens[n], t, t);
2285 /* Returns 10**(-2**N). */
2287 static const REAL_VALUE_TYPE *
2288 ten_to_mptwo (int n)
2290 static REAL_VALUE_TYPE tens[EXP_BITS];
2292 gcc_assert (n >= 0);
2293 gcc_assert (n < EXP_BITS);
2295 if (tens[n].cl == rvc_zero)
2296 do_divide (&tens[n], real_digit (1), ten_to_ptwo (n));
2303 static const REAL_VALUE_TYPE *
2306 static REAL_VALUE_TYPE num[10];
2308 gcc_assert (n >= 0);
2309 gcc_assert (n <= 9);
2311 if (n > 0 && num[n].cl == rvc_zero)
2312 real_from_integer (&num[n], VOIDmode, n, 0, 1);
2317 /* Multiply R by 10**EXP. */
2320 times_pten (REAL_VALUE_TYPE *r, int exp)
2322 REAL_VALUE_TYPE pten, *rr;
2323 bool negative = (exp < 0);
2329 pten = *real_digit (1);
2335 for (i = 0; exp > 0; ++i, exp >>= 1)
2337 do_multiply (rr, rr, ten_to_ptwo (i));
2340 do_divide (r, r, &pten);
2343 /* Returns the special REAL_VALUE_TYPE corresponding to 'e'. */
2345 const REAL_VALUE_TYPE *
2348 static REAL_VALUE_TYPE value;
2350 /* Initialize mathematical constants for constant folding builtins.
2351 These constants need to be given to at least 160 bits precision. */
2352 if (value.cl == rvc_zero)
2355 mpfr_init2 (m, SIGNIFICAND_BITS);
2356 mpfr_set_ui (m, 1, GMP_RNDN);
2357 mpfr_exp (m, m, GMP_RNDN);
2358 real_from_mpfr (&value, m, NULL_TREE, GMP_RNDN);
2365 /* Returns the special REAL_VALUE_TYPE corresponding to 1/3. */
2367 const REAL_VALUE_TYPE *
2368 dconst_third_ptr (void)
2370 static REAL_VALUE_TYPE value;
2372 /* Initialize mathematical constants for constant folding builtins.
2373 These constants need to be given to at least 160 bits precision. */
2374 if (value.cl == rvc_zero)
2376 real_arithmetic (&value, RDIV_EXPR, &dconst1, real_digit (3));
2381 /* Returns the special REAL_VALUE_TYPE corresponding to sqrt(2). */
2383 const REAL_VALUE_TYPE *
2384 dconst_sqrt2_ptr (void)
2386 static REAL_VALUE_TYPE value;
2388 /* Initialize mathematical constants for constant folding builtins.
2389 These constants need to be given to at least 160 bits precision. */
2390 if (value.cl == rvc_zero)
2393 mpfr_init2 (m, SIGNIFICAND_BITS);
2394 mpfr_sqrt_ui (m, 2, GMP_RNDN);
2395 real_from_mpfr (&value, m, NULL_TREE, GMP_RNDN);
2401 /* Fills R with +Inf. */
2404 real_inf (REAL_VALUE_TYPE *r)
2409 /* Fills R with a NaN whose significand is described by STR. If QUIET,
2410 we force a QNaN, else we force an SNaN. The string, if not empty,
2411 is parsed as a number and placed in the significand. Return true
2412 if the string was successfully parsed. */
2415 real_nan (REAL_VALUE_TYPE *r, const char *str, int quiet,
2416 enum machine_mode mode)
2418 const struct real_format *fmt;
2420 fmt = REAL_MODE_FORMAT (mode);
2426 get_canonical_qnan (r, 0);
2428 get_canonical_snan (r, 0);
2434 memset (r, 0, sizeof (*r));
2437 /* Parse akin to strtol into the significand of R. */
2439 while (ISSPACE (*str))
2443 else if (*str == '+')
2448 if (*str == 'x' || *str == 'X')
2457 while ((d = hex_value (*str)) < base)
2464 lshift_significand (r, r, 3);
2467 lshift_significand (r, r, 4);
2470 lshift_significand_1 (&u, r);
2471 lshift_significand (r, r, 3);
2472 add_significands (r, r, &u);
2480 add_significands (r, r, &u);
2485 /* Must have consumed the entire string for success. */
2489 /* Shift the significand into place such that the bits
2490 are in the most significant bits for the format. */
2491 lshift_significand (r, r, SIGNIFICAND_BITS - fmt->pnan);
2493 /* Our MSB is always unset for NaNs. */
2494 r->sig[SIGSZ-1] &= ~SIG_MSB;
2496 /* Force quiet or signalling NaN. */
2497 r->signalling = !quiet;
2503 /* Fills R with the largest finite value representable in mode MODE.
2504 If SIGN is nonzero, R is set to the most negative finite value. */
2507 real_maxval (REAL_VALUE_TYPE *r, int sign, enum machine_mode mode)
2509 const struct real_format *fmt;
2512 fmt = REAL_MODE_FORMAT (mode);
2514 memset (r, 0, sizeof (*r));
2517 decimal_real_maxval (r, sign, mode);
2522 SET_REAL_EXP (r, fmt->emax);
2524 np2 = SIGNIFICAND_BITS - fmt->p;
2525 memset (r->sig, -1, SIGSZ * sizeof (unsigned long));
2526 clear_significand_below (r, np2);
2528 if (fmt->pnan < fmt->p)
2529 /* This is an IBM extended double format made up of two IEEE
2530 doubles. The value of the long double is the sum of the
2531 values of the two parts. The most significant part is
2532 required to be the value of the long double rounded to the
2533 nearest double. Rounding means we need a slightly smaller
2534 value for LDBL_MAX. */
2535 clear_significand_bit (r, SIGNIFICAND_BITS - fmt->pnan - 1);
2539 /* Fills R with 2**N. */
2542 real_2expN (REAL_VALUE_TYPE *r, int n, enum machine_mode fmode)
2544 memset (r, 0, sizeof (*r));
2549 else if (n < -MAX_EXP)
2554 SET_REAL_EXP (r, n);
2555 r->sig[SIGSZ-1] = SIG_MSB;
2557 if (DECIMAL_FLOAT_MODE_P (fmode))
2558 decimal_real_convert (r, fmode, r);
2563 round_for_format (const struct real_format *fmt, REAL_VALUE_TYPE *r)
2567 bool round_up = false;
2573 decimal_round_for_format (fmt, r);
2576 /* FIXME. We can come here via fp_easy_constant
2577 (e.g. -O0 on '_Decimal32 x = 1.0 + 2.0dd'), but have not
2578 investigated whether this convert needs to be here, or
2579 something else is missing. */
2580 decimal_real_convert (r, DFmode, r);
2584 emin2m1 = fmt->emin - 1;
2587 np2 = SIGNIFICAND_BITS - p2;
2591 get_zero (r, r->sign);
2593 if (!fmt->has_signed_zero)
2598 get_inf (r, r->sign);
2603 clear_significand_below (r, np2);
2613 /* Check the range of the exponent. If we're out of range,
2614 either underflow or overflow. */
2615 if (REAL_EXP (r) > emax2)
2617 else if (REAL_EXP (r) <= emin2m1)
2621 if (!fmt->has_denorm)
2623 /* Don't underflow completely until we've had a chance to round. */
2624 if (REAL_EXP (r) < emin2m1)
2629 diff = emin2m1 - REAL_EXP (r) + 1;
2633 /* De-normalize the significand. */
2634 r->sig[0] |= sticky_rshift_significand (r, r, diff);
2635 SET_REAL_EXP (r, REAL_EXP (r) + diff);
2639 if (!fmt->round_towards_zero)
2641 /* There are P2 true significand bits, followed by one guard bit,
2642 followed by one sticky bit, followed by stuff. Fold nonzero
2643 stuff into the sticky bit. */
2644 unsigned long sticky;
2648 for (i = 0, w = (np2 - 1) / HOST_BITS_PER_LONG; i < w; ++i)
2649 sticky |= r->sig[i];
2651 & (((unsigned long)1 << ((np2 - 1) % HOST_BITS_PER_LONG)) - 1);
2653 guard = test_significand_bit (r, np2 - 1);
2654 lsb = test_significand_bit (r, np2);
2656 /* Round to even. */
2657 round_up = guard && (sticky || lsb);
2664 set_significand_bit (&u, np2);
2666 if (add_significands (r, r, &u))
2668 /* Overflow. Means the significand had been all ones, and
2669 is now all zeros. Need to increase the exponent, and
2670 possibly re-normalize it. */
2671 SET_REAL_EXP (r, REAL_EXP (r) + 1);
2672 if (REAL_EXP (r) > emax2)
2674 r->sig[SIGSZ-1] = SIG_MSB;
2678 /* Catch underflow that we deferred until after rounding. */
2679 if (REAL_EXP (r) <= emin2m1)
2682 /* Clear out trailing garbage. */
2683 clear_significand_below (r, np2);
2686 /* Extend or truncate to a new mode. */
2689 real_convert (REAL_VALUE_TYPE *r, enum machine_mode mode,
2690 const REAL_VALUE_TYPE *a)
2692 const struct real_format *fmt;
2694 fmt = REAL_MODE_FORMAT (mode);
2699 if (a->decimal || fmt->b == 10)
2700 decimal_real_convert (r, mode, a);
2702 round_for_format (fmt, r);
2704 /* round_for_format de-normalizes denormals. Undo just that part. */
2705 if (r->cl == rvc_normal)
2709 /* Legacy. Likewise, except return the struct directly. */
2712 real_value_truncate (enum machine_mode mode, REAL_VALUE_TYPE a)
2715 real_convert (&r, mode, &a);
2719 /* Return true if truncating to MODE is exact. */
2722 exact_real_truncate (enum machine_mode mode, const REAL_VALUE_TYPE *a)
2724 const struct real_format *fmt;
2728 fmt = REAL_MODE_FORMAT (mode);
2731 /* Don't allow conversion to denormals. */
2732 emin2m1 = fmt->emin - 1;
2733 if (REAL_EXP (a) <= emin2m1)
2736 /* After conversion to the new mode, the value must be identical. */
2737 real_convert (&t, mode, a);
2738 return real_identical (&t, a);
2741 /* Write R to the given target format. Place the words of the result
2742 in target word order in BUF. There are always 32 bits in each
2743 long, no matter the size of the host long.
2745 Legacy: return word 0 for implementing REAL_VALUE_TO_TARGET_SINGLE. */
2748 real_to_target_fmt (long *buf, const REAL_VALUE_TYPE *r_orig,
2749 const struct real_format *fmt)
2755 round_for_format (fmt, &r);
2759 (*fmt->encode) (fmt, buf, &r);
2764 /* Similar, but look up the format from MODE. */
2767 real_to_target (long *buf, const REAL_VALUE_TYPE *r, enum machine_mode mode)
2769 const struct real_format *fmt;
2771 fmt = REAL_MODE_FORMAT (mode);
2774 return real_to_target_fmt (buf, r, fmt);
2777 /* Read R from the given target format. Read the words of the result
2778 in target word order in BUF. There are always 32 bits in each
2779 long, no matter the size of the host long. */
2782 real_from_target_fmt (REAL_VALUE_TYPE *r, const long *buf,
2783 const struct real_format *fmt)
2785 (*fmt->decode) (fmt, r, buf);
2788 /* Similar, but look up the format from MODE. */
2791 real_from_target (REAL_VALUE_TYPE *r, const long *buf, enum machine_mode mode)
2793 const struct real_format *fmt;
2795 fmt = REAL_MODE_FORMAT (mode);
2798 (*fmt->decode) (fmt, r, buf);
2801 /* Return the number of bits of the largest binary value that the
2802 significand of MODE will hold. */
2803 /* ??? Legacy. Should get access to real_format directly. */
2806 significand_size (enum machine_mode mode)
2808 const struct real_format *fmt;
2810 fmt = REAL_MODE_FORMAT (mode);
2816 /* Return the size in bits of the largest binary value that can be
2817 held by the decimal coefficient for this mode. This is one more
2818 than the number of bits required to hold the largest coefficient
2820 double log2_10 = 3.3219281;
2821 return fmt->p * log2_10;
2826 /* Return a hash value for the given real value. */
2827 /* ??? The "unsigned int" return value is intended to be hashval_t,
2828 but I didn't want to pull hashtab.h into real.h. */
2831 real_hash (const REAL_VALUE_TYPE *r)
2836 h = r->cl | (r->sign << 2);
2844 h |= REAL_EXP (r) << 3;
2849 h ^= (unsigned int)-1;
2858 if (sizeof(unsigned long) > sizeof(unsigned int))
2859 for (i = 0; i < SIGSZ; ++i)
2861 unsigned long s = r->sig[i];
2862 h ^= s ^ (s >> (HOST_BITS_PER_LONG / 2));
2865 for (i = 0; i < SIGSZ; ++i)
2871 /* IEEE single-precision format. */
2873 static void encode_ieee_single (const struct real_format *fmt,
2874 long *, const REAL_VALUE_TYPE *);
2875 static void decode_ieee_single (const struct real_format *,
2876 REAL_VALUE_TYPE *, const long *);
2879 encode_ieee_single (const struct real_format *fmt, long *buf,
2880 const REAL_VALUE_TYPE *r)
2882 unsigned long image, sig, exp;
2883 unsigned long sign = r->sign;
2884 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2887 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
2898 image |= 0x7fffffff;
2905 sig = (fmt->canonical_nan_lsbs_set ? (1 << 22) - 1 : 0);
2906 if (r->signalling == fmt->qnan_msb_set)
2917 image |= 0x7fffffff;
2921 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2922 whereas the intermediate representation is 0.F x 2**exp.
2923 Which means we're off by one. */
2927 exp = REAL_EXP (r) + 127 - 1;
2940 decode_ieee_single (const struct real_format *fmt, REAL_VALUE_TYPE *r,
2943 unsigned long image = buf[0] & 0xffffffff;
2944 bool sign = (image >> 31) & 1;
2945 int exp = (image >> 23) & 0xff;
2947 memset (r, 0, sizeof (*r));
2948 image <<= HOST_BITS_PER_LONG - 24;
2953 if (image && fmt->has_denorm)
2957 SET_REAL_EXP (r, -126);
2958 r->sig[SIGSZ-1] = image << 1;
2961 else if (fmt->has_signed_zero)
2964 else if (exp == 255 && (fmt->has_nans || fmt->has_inf))
2970 r->signalling = (((image >> (HOST_BITS_PER_LONG - 2)) & 1)
2971 ^ fmt->qnan_msb_set);
2972 r->sig[SIGSZ-1] = image;
2984 SET_REAL_EXP (r, exp - 127 + 1);
2985 r->sig[SIGSZ-1] = image | SIG_MSB;
2989 const struct real_format ieee_single_format =
3010 const struct real_format mips_single_format =
3031 const struct real_format motorola_single_format =
3052 /* SPU Single Precision (Extended-Range Mode) format is the same as IEEE
3053 single precision with the following differences:
3054 - Infinities are not supported. Instead MAX_FLOAT or MIN_FLOAT
3056 - NaNs are not supported.
3057 - The range of non-zero numbers in binary is
3058 (001)[1.]000...000 to (255)[1.]111...111.
3059 - Denormals can be represented, but are treated as +0.0 when
3060 used as an operand and are never generated as a result.
3061 - -0.0 can be represented, but a zero result is always +0.0.
3062 - the only supported rounding mode is trunction (towards zero). */
3063 const struct real_format spu_single_format =
3084 /* IEEE double-precision format. */
3086 static void encode_ieee_double (const struct real_format *fmt,
3087 long *, const REAL_VALUE_TYPE *);
3088 static void decode_ieee_double (const struct real_format *,
3089 REAL_VALUE_TYPE *, const long *);
3092 encode_ieee_double (const struct real_format *fmt, long *buf,
3093 const REAL_VALUE_TYPE *r)
3095 unsigned long image_lo, image_hi, sig_lo, sig_hi, exp;
3096 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3098 image_hi = r->sign << 31;
3101 if (HOST_BITS_PER_LONG == 64)
3103 sig_hi = r->sig[SIGSZ-1];
3104 sig_lo = (sig_hi >> (64 - 53)) & 0xffffffff;
3105 sig_hi = (sig_hi >> (64 - 53 + 1) >> 31) & 0xfffff;
3109 sig_hi = r->sig[SIGSZ-1];
3110 sig_lo = r->sig[SIGSZ-2];
3111 sig_lo = (sig_hi << 21) | (sig_lo >> 11);
3112 sig_hi = (sig_hi >> 11) & 0xfffff;
3122 image_hi |= 2047 << 20;
3125 image_hi |= 0x7fffffff;
3126 image_lo = 0xffffffff;
3135 if (fmt->canonical_nan_lsbs_set)
3137 sig_hi = (1 << 19) - 1;
3138 sig_lo = 0xffffffff;
3146 if (r->signalling == fmt->qnan_msb_set)
3147 sig_hi &= ~(1 << 19);
3150 if (sig_hi == 0 && sig_lo == 0)
3153 image_hi |= 2047 << 20;
3159 image_hi |= 0x7fffffff;
3160 image_lo = 0xffffffff;
3165 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3166 whereas the intermediate representation is 0.F x 2**exp.
3167 Which means we're off by one. */
3171 exp = REAL_EXP (r) + 1023 - 1;
3172 image_hi |= exp << 20;
3181 if (FLOAT_WORDS_BIG_ENDIAN)
3182 buf[0] = image_hi, buf[1] = image_lo;
3184 buf[0] = image_lo, buf[1] = image_hi;
3188 decode_ieee_double (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3191 unsigned long image_hi, image_lo;
3195 if (FLOAT_WORDS_BIG_ENDIAN)
3196 image_hi = buf[0], image_lo = buf[1];
3198 image_lo = buf[0], image_hi = buf[1];
3199 image_lo &= 0xffffffff;
3200 image_hi &= 0xffffffff;
3202 sign = (image_hi >> 31) & 1;
3203 exp = (image_hi >> 20) & 0x7ff;
3205 memset (r, 0, sizeof (*r));
3207 image_hi <<= 32 - 21;
3208 image_hi |= image_lo >> 21;
3209 image_hi &= 0x7fffffff;
3210 image_lo <<= 32 - 21;
3214 if ((image_hi || image_lo) && fmt->has_denorm)
3218 SET_REAL_EXP (r, -1022);
3219 if (HOST_BITS_PER_LONG == 32)
3221 image_hi = (image_hi << 1) | (image_lo >> 31);
3223 r->sig[SIGSZ-1] = image_hi;
3224 r->sig[SIGSZ-2] = image_lo;
3228 image_hi = (image_hi << 31 << 2) | (image_lo << 1);
3229 r->sig[SIGSZ-1] = image_hi;
3233 else if (fmt->has_signed_zero)
3236 else if (exp == 2047 && (fmt->has_nans || fmt->has_inf))
3238 if (image_hi || image_lo)
3242 r->signalling = ((image_hi >> 30) & 1) ^ fmt->qnan_msb_set;
3243 if (HOST_BITS_PER_LONG == 32)
3245 r->sig[SIGSZ-1] = image_hi;
3246 r->sig[SIGSZ-2] = image_lo;
3249 r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo;
3261 SET_REAL_EXP (r, exp - 1023 + 1);
3262 if (HOST_BITS_PER_LONG == 32)
3264 r->sig[SIGSZ-1] = image_hi | SIG_MSB;
3265 r->sig[SIGSZ-2] = image_lo;
3268 r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo | SIG_MSB;
3272 const struct real_format ieee_double_format =
3293 const struct real_format mips_double_format =
3314 const struct real_format motorola_double_format =
3335 /* IEEE extended real format. This comes in three flavors: Intel's as
3336 a 12 byte image, Intel's as a 16 byte image, and Motorola's. Intel
3337 12- and 16-byte images may be big- or little endian; Motorola's is
3338 always big endian. */
3340 /* Helper subroutine which converts from the internal format to the
3341 12-byte little-endian Intel format. Functions below adjust this
3342 for the other possible formats. */
3344 encode_ieee_extended (const struct real_format *fmt, long *buf,
3345 const REAL_VALUE_TYPE *r)
3347 unsigned long image_hi, sig_hi, sig_lo;
3348 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3350 image_hi = r->sign << 15;
3351 sig_hi = sig_lo = 0;
3363 /* Intel requires the explicit integer bit to be set, otherwise
3364 it considers the value a "pseudo-infinity". Motorola docs
3365 say it doesn't care. */
3366 sig_hi = 0x80000000;
3371 sig_lo = sig_hi = 0xffffffff;
3381 if (fmt->canonical_nan_lsbs_set)
3383 sig_hi = (1 << 30) - 1;
3384 sig_lo = 0xffffffff;
3387 else if (HOST_BITS_PER_LONG == 32)
3389 sig_hi = r->sig[SIGSZ-1];
3390 sig_lo = r->sig[SIGSZ-2];
3394 sig_lo = r->sig[SIGSZ-1];
3395 sig_hi = sig_lo >> 31 >> 1;
3396 sig_lo &= 0xffffffff;
3398 if (r->signalling == fmt->qnan_msb_set)
3399 sig_hi &= ~(1 << 30);
3402 if ((sig_hi & 0x7fffffff) == 0 && sig_lo == 0)
3405 /* Intel requires the explicit integer bit to be set, otherwise
3406 it considers the value a "pseudo-nan". Motorola docs say it
3408 sig_hi |= 0x80000000;
3413 sig_lo = sig_hi = 0xffffffff;
3419 int exp = REAL_EXP (r);
3421 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3422 whereas the intermediate representation is 0.F x 2**exp.
3423 Which means we're off by one.
3425 Except for Motorola, which consider exp=0 and explicit
3426 integer bit set to continue to be normalized. In theory
3427 this discrepancy has been taken care of by the difference
3428 in fmt->emin in round_for_format. */
3435 gcc_assert (exp >= 0);
3439 if (HOST_BITS_PER_LONG == 32)
3441 sig_hi = r->sig[SIGSZ-1];
3442 sig_lo = r->sig[SIGSZ-2];
3446 sig_lo = r->sig[SIGSZ-1];
3447 sig_hi = sig_lo >> 31 >> 1;
3448 sig_lo &= 0xffffffff;
3457 buf[0] = sig_lo, buf[1] = sig_hi, buf[2] = image_hi;
3460 /* Convert from the internal format to the 12-byte Motorola format
3461 for an IEEE extended real. */
3463 encode_ieee_extended_motorola (const struct real_format *fmt, long *buf,
3464 const REAL_VALUE_TYPE *r)
3467 encode_ieee_extended (fmt, intermed, r);
3469 /* Motorola chips are assumed always to be big-endian. Also, the
3470 padding in a Motorola extended real goes between the exponent and
3471 the mantissa. At this point the mantissa is entirely within
3472 elements 0 and 1 of intermed, and the exponent entirely within
3473 element 2, so all we have to do is swap the order around, and
3474 shift element 2 left 16 bits. */
3475 buf[0] = intermed[2] << 16;
3476 buf[1] = intermed[1];
3477 buf[2] = intermed[0];
3480 /* Convert from the internal format to the 12-byte Intel format for
3481 an IEEE extended real. */
3483 encode_ieee_extended_intel_96 (const struct real_format *fmt, long *buf,
3484 const REAL_VALUE_TYPE *r)
3486 if (FLOAT_WORDS_BIG_ENDIAN)
3488 /* All the padding in an Intel-format extended real goes at the high
3489 end, which in this case is after the mantissa, not the exponent.
3490 Therefore we must shift everything down 16 bits. */
3492 encode_ieee_extended (fmt, intermed, r);
3493 buf[0] = ((intermed[2] << 16) | ((unsigned long)(intermed[1] & 0xFFFF0000) >> 16));
3494 buf[1] = ((intermed[1] << 16) | ((unsigned long)(intermed[0] & 0xFFFF0000) >> 16));
3495 buf[2] = (intermed[0] << 16);
3498 /* encode_ieee_extended produces what we want directly. */
3499 encode_ieee_extended (fmt, buf, r);
3502 /* Convert from the internal format to the 16-byte Intel format for
3503 an IEEE extended real. */
3505 encode_ieee_extended_intel_128 (const struct real_format *fmt, long *buf,
3506 const REAL_VALUE_TYPE *r)
3508 /* All the padding in an Intel-format extended real goes at the high end. */
3509 encode_ieee_extended_intel_96 (fmt, buf, r);
3513 /* As above, we have a helper function which converts from 12-byte
3514 little-endian Intel format to internal format. Functions below
3515 adjust for the other possible formats. */
3517 decode_ieee_extended (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3520 unsigned long image_hi, sig_hi, sig_lo;
3524 sig_lo = buf[0], sig_hi = buf[1], image_hi = buf[2];
3525 sig_lo &= 0xffffffff;
3526 sig_hi &= 0xffffffff;
3527 image_hi &= 0xffffffff;
3529 sign = (image_hi >> 15) & 1;
3530 exp = image_hi & 0x7fff;
3532 memset (r, 0, sizeof (*r));
3536 if ((sig_hi || sig_lo) && fmt->has_denorm)
3541 /* When the IEEE format contains a hidden bit, we know that
3542 it's zero at this point, and so shift up the significand
3543 and decrease the exponent to match. In this case, Motorola
3544 defines the explicit integer bit to be valid, so we don't
3545 know whether the msb is set or not. */
3546 SET_REAL_EXP (r, fmt->emin);
3547 if (HOST_BITS_PER_LONG == 32)
3549 r->sig[SIGSZ-1] = sig_hi;
3550 r->sig[SIGSZ-2] = sig_lo;
3553 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3557 else if (fmt->has_signed_zero)
3560 else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3562 /* See above re "pseudo-infinities" and "pseudo-nans".
3563 Short summary is that the MSB will likely always be
3564 set, and that we don't care about it. */
3565 sig_hi &= 0x7fffffff;
3567 if (sig_hi || sig_lo)
3571 r->signalling = ((sig_hi >> 30) & 1) ^ fmt->qnan_msb_set;
3572 if (HOST_BITS_PER_LONG == 32)
3574 r->sig[SIGSZ-1] = sig_hi;
3575 r->sig[SIGSZ-2] = sig_lo;
3578 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3590 SET_REAL_EXP (r, exp - 16383 + 1);
3591 if (HOST_BITS_PER_LONG == 32)
3593 r->sig[SIGSZ-1] = sig_hi;
3594 r->sig[SIGSZ-2] = sig_lo;
3597 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3601 /* Convert from the internal format to the 12-byte Motorola format
3602 for an IEEE extended real. */
3604 decode_ieee_extended_motorola (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3609 /* Motorola chips are assumed always to be big-endian. Also, the
3610 padding in a Motorola extended real goes between the exponent and
3611 the mantissa; remove it. */
3612 intermed[0] = buf[2];
3613 intermed[1] = buf[1];
3614 intermed[2] = (unsigned long)buf[0] >> 16;
3616 decode_ieee_extended (fmt, r, intermed);
3619 /* Convert from the internal format to the 12-byte Intel format for
3620 an IEEE extended real. */
3622 decode_ieee_extended_intel_96 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3625 if (FLOAT_WORDS_BIG_ENDIAN)
3627 /* All the padding in an Intel-format extended real goes at the high
3628 end, which in this case is after the mantissa, not the exponent.
3629 Therefore we must shift everything up 16 bits. */
3632 intermed[0] = (((unsigned long)buf[2] >> 16) | (buf[1] << 16));
3633 intermed[1] = (((unsigned long)buf[1] >> 16) | (buf[0] << 16));
3634 intermed[2] = ((unsigned long)buf[0] >> 16);
3636 decode_ieee_extended (fmt, r, intermed);
3639 /* decode_ieee_extended produces what we want directly. */
3640 decode_ieee_extended (fmt, r, buf);
3643 /* Convert from the internal format to the 16-byte Intel format for
3644 an IEEE extended real. */
3646 decode_ieee_extended_intel_128 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3649 /* All the padding in an Intel-format extended real goes at the high end. */
3650 decode_ieee_extended_intel_96 (fmt, r, buf);
3653 const struct real_format ieee_extended_motorola_format =
3655 encode_ieee_extended_motorola,
3656 decode_ieee_extended_motorola,
3674 const struct real_format ieee_extended_intel_96_format =
3676 encode_ieee_extended_intel_96,
3677 decode_ieee_extended_intel_96,
3695 const struct real_format ieee_extended_intel_128_format =
3697 encode_ieee_extended_intel_128,
3698 decode_ieee_extended_intel_128,
3716 /* The following caters to i386 systems that set the rounding precision
3717 to 53 bits instead of 64, e.g. FreeBSD. */
3718 const struct real_format ieee_extended_intel_96_round_53_format =
3720 encode_ieee_extended_intel_96,
3721 decode_ieee_extended_intel_96,
3739 /* IBM 128-bit extended precision format: a pair of IEEE double precision
3740 numbers whose sum is equal to the extended precision value. The number
3741 with greater magnitude is first. This format has the same magnitude
3742 range as an IEEE double precision value, but effectively 106 bits of
3743 significand precision. Infinity and NaN are represented by their IEEE
3744 double precision value stored in the first number, the second number is
3745 +0.0 or -0.0 for Infinity and don't-care for NaN. */
3747 static void encode_ibm_extended (const struct real_format *fmt,
3748 long *, const REAL_VALUE_TYPE *);
3749 static void decode_ibm_extended (const struct real_format *,
3750 REAL_VALUE_TYPE *, const long *);
3753 encode_ibm_extended (const struct real_format *fmt, long *buf,
3754 const REAL_VALUE_TYPE *r)
3756 REAL_VALUE_TYPE u, normr, v;
3757 const struct real_format *base_fmt;
3759 base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3761 /* Renormalize R before doing any arithmetic on it. */
3763 if (normr.cl == rvc_normal)
3766 /* u = IEEE double precision portion of significand. */
3768 round_for_format (base_fmt, &u);
3769 encode_ieee_double (base_fmt, &buf[0], &u);
3771 if (u.cl == rvc_normal)
3773 do_add (&v, &normr, &u, 1);
3774 /* Call round_for_format since we might need to denormalize. */
3775 round_for_format (base_fmt, &v);
3776 encode_ieee_double (base_fmt, &buf[2], &v);
3780 /* Inf, NaN, 0 are all representable as doubles, so the
3781 least-significant part can be 0.0. */
3788 decode_ibm_extended (const struct real_format *fmt ATTRIBUTE_UNUSED, REAL_VALUE_TYPE *r,
3791 REAL_VALUE_TYPE u, v;
3792 const struct real_format *base_fmt;
3794 base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3795 decode_ieee_double (base_fmt, &u, &buf[0]);
3797 if (u.cl != rvc_zero && u.cl != rvc_inf && u.cl != rvc_nan)
3799 decode_ieee_double (base_fmt, &v, &buf[2]);
3800 do_add (r, &u, &v, 0);
3806 const struct real_format ibm_extended_format =
3808 encode_ibm_extended,
3809 decode_ibm_extended,
3827 const struct real_format mips_extended_format =
3829 encode_ibm_extended,
3830 decode_ibm_extended,
3849 /* IEEE quad precision format. */
3851 static void encode_ieee_quad (const struct real_format *fmt,
3852 long *, const REAL_VALUE_TYPE *);
3853 static void decode_ieee_quad (const struct real_format *,
3854 REAL_VALUE_TYPE *, const long *);
3857 encode_ieee_quad (const struct real_format *fmt, long *buf,
3858 const REAL_VALUE_TYPE *r)
3860 unsigned long image3, image2, image1, image0, exp;
3861 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3864 image3 = r->sign << 31;
3869 rshift_significand (&u, r, SIGNIFICAND_BITS - 113);
3878 image3 |= 32767 << 16;
3881 image3 |= 0x7fffffff;
3882 image2 = 0xffffffff;
3883 image1 = 0xffffffff;
3884 image0 = 0xffffffff;
3891 image3 |= 32767 << 16;
3895 if (fmt->canonical_nan_lsbs_set)
3898 image2 = image1 = image0 = 0xffffffff;
3901 else if (HOST_BITS_PER_LONG == 32)
3906 image3 |= u.sig[3] & 0xffff;
3911 image1 = image0 >> 31 >> 1;
3913 image3 |= (image2 >> 31 >> 1) & 0xffff;
3914 image0 &= 0xffffffff;
3915 image2 &= 0xffffffff;
3917 if (r->signalling == fmt->qnan_msb_set)
3921 if (((image3 & 0xffff) | image2 | image1 | image0) == 0)
3926 image3 |= 0x7fffffff;
3927 image2 = 0xffffffff;
3928 image1 = 0xffffffff;
3929 image0 = 0xffffffff;
3934 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3935 whereas the intermediate representation is 0.F x 2**exp.
3936 Which means we're off by one. */
3940 exp = REAL_EXP (r) + 16383 - 1;
3941 image3 |= exp << 16;
3943 if (HOST_BITS_PER_LONG == 32)
3948 image3 |= u.sig[3] & 0xffff;
3953 image1 = image0 >> 31 >> 1;
3955 image3 |= (image2 >> 31 >> 1) & 0xffff;
3956 image0 &= 0xffffffff;
3957 image2 &= 0xffffffff;
3965 if (FLOAT_WORDS_BIG_ENDIAN)
3982 decode_ieee_quad (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3985 unsigned long image3, image2, image1, image0;
3989 if (FLOAT_WORDS_BIG_ENDIAN)
4003 image0 &= 0xffffffff;
4004 image1 &= 0xffffffff;
4005 image2 &= 0xffffffff;
4007 sign = (image3 >> 31) & 1;
4008 exp = (image3 >> 16) & 0x7fff;
4011 memset (r, 0, sizeof (*r));
4015 if ((image3 | image2 | image1 | image0) && fmt->has_denorm)
4020 SET_REAL_EXP (r, -16382 + (SIGNIFICAND_BITS - 112));
4021 if (HOST_BITS_PER_LONG == 32)
4030 r->sig[0] = (image1 << 31 << 1) | image0;
4031 r->sig[1] = (image3 << 31 << 1) | image2;
4036 else if (fmt->has_signed_zero)
4039 else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
4041 if (image3 | image2 | image1 | image0)
4045 r->signalling = ((image3 >> 15) & 1) ^ fmt->qnan_msb_set;
4047 if (HOST_BITS_PER_LONG == 32)
4056 r->sig[0] = (image1 << 31 << 1) | image0;
4057 r->sig[1] = (image3 << 31 << 1) | image2;
4059 lshift_significand (r, r, SIGNIFICAND_BITS - 113);
4071 SET_REAL_EXP (r, exp - 16383 + 1);
4073 if (HOST_BITS_PER_LONG == 32)
4082 r->sig[0] = (image1 << 31 << 1) | image0;
4083 r->sig[1] = (image3 << 31 << 1) | image2;
4085 lshift_significand (r, r, SIGNIFICAND_BITS - 113);
4086 r->sig[SIGSZ-1] |= SIG_MSB;
4090 const struct real_format ieee_quad_format =
4111 const struct real_format mips_quad_format =
4132 /* Descriptions of VAX floating point formats can be found beginning at
4134 http://h71000.www7.hp.com/doc/73FINAL/4515/4515pro_013.html#f_floating_point_format
4136 The thing to remember is that they're almost IEEE, except for word
4137 order, exponent bias, and the lack of infinities, nans, and denormals.
4139 We don't implement the H_floating format here, simply because neither
4140 the VAX or Alpha ports use it. */
4142 static void encode_vax_f (const struct real_format *fmt,
4143 long *, const REAL_VALUE_TYPE *);
4144 static void decode_vax_f (const struct real_format *,
4145 REAL_VALUE_TYPE *, const long *);
4146 static void encode_vax_d (const struct real_format *fmt,
4147 long *, const REAL_VALUE_TYPE *);
4148 static void decode_vax_d (const struct real_format *,
4149 REAL_VALUE_TYPE *, const long *);
4150 static void encode_vax_g (const struct real_format *fmt,
4151 long *, const REAL_VALUE_TYPE *);
4152 static void decode_vax_g (const struct real_format *,
4153 REAL_VALUE_TYPE *, const long *);
4156 encode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4157 const REAL_VALUE_TYPE *r)
4159 unsigned long sign, exp, sig, image;
4161 sign = r->sign << 15;
4171 image = 0xffff7fff | sign;
4175 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
4176 exp = REAL_EXP (r) + 128;
4178 image = (sig << 16) & 0xffff0000;
4192 decode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED,
4193 REAL_VALUE_TYPE *r, const long *buf)
4195 unsigned long image = buf[0] & 0xffffffff;
4196 int exp = (image >> 7) & 0xff;
4198 memset (r, 0, sizeof (*r));
4203 r->sign = (image >> 15) & 1;
4204 SET_REAL_EXP (r, exp - 128);
4206 image = ((image & 0x7f) << 16) | ((image >> 16) & 0xffff);
4207 r->sig[SIGSZ-1] = (image << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
4212 encode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4213 const REAL_VALUE_TYPE *r)
4215 unsigned long image0, image1, sign = r->sign << 15;
4220 image0 = image1 = 0;
4225 image0 = 0xffff7fff | sign;
4226 image1 = 0xffffffff;
4230 /* Extract the significand into straight hi:lo. */
4231 if (HOST_BITS_PER_LONG == 64)
4233 image0 = r->sig[SIGSZ-1];
4234 image1 = (image0 >> (64 - 56)) & 0xffffffff;
4235 image0 = (image0 >> (64 - 56 + 1) >> 31) & 0x7fffff;
4239 image0 = r->sig[SIGSZ-1];
4240 image1 = r->sig[SIGSZ-2];
4241 image1 = (image0 << 24) | (image1 >> 8);
4242 image0 = (image0 >> 8) & 0xffffff;
4245 /* Rearrange the half-words of the significand to match the
4247 image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff007f;
4248 image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
4250 /* Add the sign and exponent. */
4252 image0 |= (REAL_EXP (r) + 128) << 7;
4259 if (FLOAT_WORDS_BIG_ENDIAN)
4260 buf[0] = image1, buf[1] = image0;
4262 buf[0] = image0, buf[1] = image1;
4266 decode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED,
4267 REAL_VALUE_TYPE *r, const long *buf)
4269 unsigned long image0, image1;
4272 if (FLOAT_WORDS_BIG_ENDIAN)
4273 image1 = buf[0], image0 = buf[1];
4275 image0 = buf[0], image1 = buf[1];
4276 image0 &= 0xffffffff;
4277 image1 &= 0xffffffff;
4279 exp = (image0 >> 7) & 0xff;
4281 memset (r, 0, sizeof (*r));
4286 r->sign = (image0 >> 15) & 1;
4287 SET_REAL_EXP (r, exp - 128);
4289 /* Rearrange the half-words of the external format into
4290 proper ascending order. */
4291 image0 = ((image0 & 0x7f) << 16) | ((image0 >> 16) & 0xffff);
4292 image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
4294 if (HOST_BITS_PER_LONG == 64)
4296 image0 = (image0 << 31 << 1) | image1;
4299 r->sig[SIGSZ-1] = image0;
4303 r->sig[SIGSZ-1] = image0;
4304 r->sig[SIGSZ-2] = image1;
4305 lshift_significand (r, r, 2*HOST_BITS_PER_LONG - 56);
4306 r->sig[SIGSZ-1] |= SIG_MSB;
4312 encode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4313 const REAL_VALUE_TYPE *r)
4315 unsigned long image0, image1, sign = r->sign << 15;
4320 image0 = image1 = 0;
4325 image0 = 0xffff7fff | sign;
4326 image1 = 0xffffffff;
4330 /* Extract the significand into straight hi:lo. */
4331 if (HOST_BITS_PER_LONG == 64)
4333 image0 = r->sig[SIGSZ-1];
4334 image1 = (image0 >> (64 - 53)) & 0xffffffff;
4335 image0 = (image0 >> (64 - 53 + 1) >> 31) & 0xfffff;
4339 image0 = r->sig[SIGSZ-1];
4340 image1 = r->sig[SIGSZ-2];
4341 image1 = (image0 << 21) | (image1 >> 11);
4342 image0 = (image0 >> 11) & 0xfffff;
4345 /* Rearrange the half-words of the significand to match the
4347 image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff000f;
4348 image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
4350 /* Add the sign and exponent. */
4352 image0 |= (REAL_EXP (r) + 1024) << 4;
4359 if (FLOAT_WORDS_BIG_ENDIAN)
4360 buf[0] = image1, buf[1] = image0;
4362 buf[0] = image0, buf[1] = image1;
4366 decode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED,
4367 REAL_VALUE_TYPE *r, const long *buf)
4369 unsigned long image0, image1;
4372 if (FLOAT_WORDS_BIG_ENDIAN)
4373 image1 = buf[0], image0 = buf[1];
4375 image0 = buf[0], image1 = buf[1];
4376 image0 &= 0xffffffff;
4377 image1 &= 0xffffffff;
4379 exp = (image0 >> 4) & 0x7ff;
4381 memset (r, 0, sizeof (*r));
4386 r->sign = (image0 >> 15) & 1;
4387 SET_REAL_EXP (r, exp - 1024);
4389 /* Rearrange the half-words of the external format into
4390 proper ascending order. */
4391 image0 = ((image0 & 0xf) << 16) | ((image0 >> 16) & 0xffff);
4392 image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
4394 if (HOST_BITS_PER_LONG == 64)
4396 image0 = (image0 << 31 << 1) | image1;
4399 r->sig[SIGSZ-1] = image0;
4403 r->sig[SIGSZ-1] = image0;
4404 r->sig[SIGSZ-2] = image1;
4405 lshift_significand (r, r, 64 - 53);
4406 r->sig[SIGSZ-1] |= SIG_MSB;
4411 const struct real_format vax_f_format =
4432 const struct real_format vax_d_format =
4453 const struct real_format vax_g_format =
4474 /* Encode real R into a single precision DFP value in BUF. */
4476 encode_decimal_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4477 long *buf ATTRIBUTE_UNUSED,
4478 const REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED)
4480 encode_decimal32 (fmt, buf, r);
4483 /* Decode a single precision DFP value in BUF into a real R. */
4485 decode_decimal_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4486 REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED,
4487 const long *buf ATTRIBUTE_UNUSED)
4489 decode_decimal32 (fmt, r, buf);
4492 /* Encode real R into a double precision DFP value in BUF. */
4494 encode_decimal_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4495 long *buf ATTRIBUTE_UNUSED,
4496 const REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED)
4498 encode_decimal64 (fmt, buf, r);
4501 /* Decode a double precision DFP value in BUF into a real R. */
4503 decode_decimal_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4504 REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED,
4505 const long *buf ATTRIBUTE_UNUSED)
4507 decode_decimal64 (fmt, r, buf);
4510 /* Encode real R into a quad precision DFP value in BUF. */
4512 encode_decimal_quad (const struct real_format *fmt ATTRIBUTE_UNUSED,
4513 long *buf ATTRIBUTE_UNUSED,
4514 const REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED)
4516 encode_decimal128 (fmt, buf, r);
4519 /* Decode a quad precision DFP value in BUF into a real R. */
4521 decode_decimal_quad (const struct real_format *fmt ATTRIBUTE_UNUSED,
4522 REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED,
4523 const long *buf ATTRIBUTE_UNUSED)
4525 decode_decimal128 (fmt, r, buf);
4528 /* Single precision decimal floating point (IEEE 754). */
4529 const struct real_format decimal_single_format =
4531 encode_decimal_single,
4532 decode_decimal_single,
4550 /* Double precision decimal floating point (IEEE 754). */
4551 const struct real_format decimal_double_format =
4553 encode_decimal_double,
4554 decode_decimal_double,
4572 /* Quad precision decimal floating point (IEEE 754). */
4573 const struct real_format decimal_quad_format =
4575 encode_decimal_quad,
4576 decode_decimal_quad,
4594 /* Encode half-precision floats. This routine is used both for the IEEE
4595 ARM alternative encodings. */
4597 encode_ieee_half (const struct real_format *fmt, long *buf,
4598 const REAL_VALUE_TYPE *r)
4600 unsigned long image, sig, exp;
4601 unsigned long sign = r->sign;
4602 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
4605 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 11)) & 0x3ff;
4623 sig = (fmt->canonical_nan_lsbs_set ? (1 << 9) - 1 : 0);
4624 if (r->signalling == fmt->qnan_msb_set)
4639 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
4640 whereas the intermediate representation is 0.F x 2**exp.
4641 Which means we're off by one. */
4645 exp = REAL_EXP (r) + 15 - 1;
4657 /* Decode half-precision floats. This routine is used both for the IEEE
4658 ARM alternative encodings. */
4660 decode_ieee_half (const struct real_format *fmt, REAL_VALUE_TYPE *r,
4663 unsigned long image = buf[0] & 0xffff;
4664 bool sign = (image >> 15) & 1;
4665 int exp = (image >> 10) & 0x1f;
4667 memset (r, 0, sizeof (*r));
4668 image <<= HOST_BITS_PER_LONG - 11;
4673 if (image && fmt->has_denorm)
4677 SET_REAL_EXP (r, -14);
4678 r->sig[SIGSZ-1] = image << 1;
4681 else if (fmt->has_signed_zero)
4684 else if (exp == 31 && (fmt->has_nans || fmt->has_inf))
4690 r->signalling = (((image >> (HOST_BITS_PER_LONG - 2)) & 1)
4691 ^ fmt->qnan_msb_set);
4692 r->sig[SIGSZ-1] = image;
4704 SET_REAL_EXP (r, exp - 15 + 1);
4705 r->sig[SIGSZ-1] = image | SIG_MSB;
4709 /* Half-precision format, as specified in IEEE 754R. */
4710 const struct real_format ieee_half_format =
4731 /* ARM's alternative half-precision format, similar to IEEE but with
4732 no reserved exponent value for NaNs and infinities; rather, it just
4733 extends the range of exponents by one. */
4734 const struct real_format arm_half_format =
4755 /* A synthetic "format" for internal arithmetic. It's the size of the
4756 internal significand minus the two bits needed for proper rounding.
4757 The encode and decode routines exist only to satisfy our paranoia
4760 static void encode_internal (const struct real_format *fmt,
4761 long *, const REAL_VALUE_TYPE *);
4762 static void decode_internal (const struct real_format *,
4763 REAL_VALUE_TYPE *, const long *);
4766 encode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4767 const REAL_VALUE_TYPE *r)
4769 memcpy (buf, r, sizeof (*r));
4773 decode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED,
4774 REAL_VALUE_TYPE *r, const long *buf)
4776 memcpy (r, buf, sizeof (*r));
4779 const struct real_format real_internal_format =
4784 SIGNIFICAND_BITS - 2,
4785 SIGNIFICAND_BITS - 2,
4800 /* Calculate the square root of X in mode MODE, and store the result
4801 in R. Return TRUE if the operation does not raise an exception.
4802 For details see "High Precision Division and Square Root",
4803 Alan H. Karp and Peter Markstein, HP Lab Report 93-93-42, June
4804 1993. http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf. */
4807 real_sqrt (REAL_VALUE_TYPE *r, enum machine_mode mode,
4808 const REAL_VALUE_TYPE *x)
4810 static REAL_VALUE_TYPE halfthree;
4811 static bool init = false;
4812 REAL_VALUE_TYPE h, t, i;
4815 /* sqrt(-0.0) is -0.0. */
4816 if (real_isnegzero (x))
4822 /* Negative arguments return NaN. */
4825 get_canonical_qnan (r, 0);
4829 /* Infinity and NaN return themselves. */
4830 if (!real_isfinite (x))
4838 do_add (&halfthree, &dconst1, &dconsthalf, 0);
4842 /* Initial guess for reciprocal sqrt, i. */
4843 exp = real_exponent (x);
4844 real_ldexp (&i, &dconst1, -exp/2);
4846 /* Newton's iteration for reciprocal sqrt, i. */
4847 for (iter = 0; iter < 16; iter++)
4849 /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x). */
4850 do_multiply (&t, x, &i);
4851 do_multiply (&h, &t, &i);
4852 do_multiply (&t, &h, &dconsthalf);
4853 do_add (&h, &halfthree, &t, 1);
4854 do_multiply (&t, &i, &h);
4856 /* Check for early convergence. */
4857 if (iter >= 6 && real_identical (&i, &t))
4860 /* ??? Unroll loop to avoid copying. */
4864 /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)). */
4865 do_multiply (&t, x, &i);
4866 do_multiply (&h, &t, &i);
4867 do_add (&i, &dconst1, &h, 1);
4868 do_multiply (&h, &t, &i);
4869 do_multiply (&i, &dconsthalf, &h);
4870 do_add (&h, &t, &i, 0);
4872 /* ??? We need a Tuckerman test to get the last bit. */
4874 real_convert (r, mode, &h);
4878 /* Calculate X raised to the integer exponent N in mode MODE and store
4879 the result in R. Return true if the result may be inexact due to
4880 loss of precision. The algorithm is the classic "left-to-right binary
4881 method" described in section 4.6.3 of Donald Knuth's "Seminumerical
4882 Algorithms", "The Art of Computer Programming", Volume 2. */
4885 real_powi (REAL_VALUE_TYPE *r, enum machine_mode mode,
4886 const REAL_VALUE_TYPE *x, HOST_WIDE_INT n)
4888 unsigned HOST_WIDE_INT bit;
4890 bool inexact = false;
4902 /* Don't worry about overflow, from now on n is unsigned. */
4910 bit = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
4911 for (i = 0; i < HOST_BITS_PER_WIDE_INT; i++)
4915 inexact |= do_multiply (&t, &t, &t);
4917 inexact |= do_multiply (&t, &t, x);
4925 inexact |= do_divide (&t, &dconst1, &t);
4927 real_convert (r, mode, &t);
4931 /* Round X to the nearest integer not larger in absolute value, i.e.
4932 towards zero, placing the result in R in mode MODE. */
4935 real_trunc (REAL_VALUE_TYPE *r, enum machine_mode mode,
4936 const REAL_VALUE_TYPE *x)
4938 do_fix_trunc (r, x);
4939 if (mode != VOIDmode)
4940 real_convert (r, mode, r);
4943 /* Round X to the largest integer not greater in value, i.e. round
4944 down, placing the result in R in mode MODE. */
4947 real_floor (REAL_VALUE_TYPE *r, enum machine_mode mode,
4948 const REAL_VALUE_TYPE *x)
4952 do_fix_trunc (&t, x);
4953 if (! real_identical (&t, x) && x->sign)
4954 do_add (&t, &t, &dconstm1, 0);
4955 if (mode != VOIDmode)
4956 real_convert (r, mode, &t);
4961 /* Round X to the smallest integer not less then argument, i.e. round
4962 up, placing the result in R in mode MODE. */
4965 real_ceil (REAL_VALUE_TYPE *r, enum machine_mode mode,
4966 const REAL_VALUE_TYPE *x)
4970 do_fix_trunc (&t, x);
4971 if (! real_identical (&t, x) && ! x->sign)
4972 do_add (&t, &t, &dconst1, 0);
4973 if (mode != VOIDmode)
4974 real_convert (r, mode, &t);
4979 /* Round X to the nearest integer, but round halfway cases away from
4983 real_round (REAL_VALUE_TYPE *r, enum machine_mode mode,
4984 const REAL_VALUE_TYPE *x)
4986 do_add (r, x, &dconsthalf, x->sign);
4987 do_fix_trunc (r, r);
4988 if (mode != VOIDmode)
4989 real_convert (r, mode, r);
4992 /* Set the sign of R to the sign of X. */
4995 real_copysign (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *x)
5000 /* Check whether the real constant value given is an integer. */
5003 real_isinteger (const REAL_VALUE_TYPE *c, enum machine_mode mode)
5005 REAL_VALUE_TYPE cint;
5007 real_trunc (&cint, mode, c);
5008 return real_identical (c, &cint);
5011 /* Write into BUF the maximum representable finite floating-point
5012 number, (1 - b**-p) * b**emax for a given FP format FMT as a hex
5013 float string. LEN is the size of BUF, and the buffer must be large
5014 enough to contain the resulting string. */
5017 get_max_float (const struct real_format *fmt, char *buf, size_t len)
5022 strcpy (buf, "0x0.");
5024 for (i = 0, p = buf + 4; i + 3 < n; i += 4)
5027 *p++ = "08ce"[n - i];
5028 sprintf (p, "p%d", fmt->emax);
5029 if (fmt->pnan < fmt->p)
5031 /* This is an IBM extended double format made up of two IEEE
5032 doubles. The value of the long double is the sum of the
5033 values of the two parts. The most significant part is
5034 required to be the value of the long double rounded to the
5035 nearest double. Rounding means we need a slightly smaller
5036 value for LDBL_MAX. */
5037 buf[4 + fmt->pnan / 4] = "7bde"[fmt->pnan % 4];
5040 gcc_assert (strlen (buf) < len);