1 /* real.c - software floating point emulation.
2 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 1999,
3 2000, 2002, 2003, 2004, 2005 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, 51 Franklin Street, Fifth Floor, Boston, MA
26 #include "coretypes.h"
34 /* The floating point model used internally is not exactly IEEE 754
35 compliant, and close to the description in the ISO C99 standard,
36 section 5.2.4.2.2 Characteristics of floating types.
40 x = s * b^e * \sum_{k=1}^p f_k * b^{-k}
44 b = base or radix, here always 2
46 p = precision (the number of base-b digits in the significand)
47 f_k = the digits of the significand.
49 We differ from typical IEEE 754 encodings in that the entire
50 significand is fractional. Normalized significands are in the
53 A requirement of the model is that P be larger than the largest
54 supported target floating-point type by at least 2 bits. This gives
55 us proper rounding when we truncate to the target type. In addition,
56 E must be large enough to hold the smallest supported denormal number
59 Both of these requirements are easily satisfied. The largest target
60 significand is 113 bits; we store at least 160. The smallest
61 denormal number fits in 17 exponent bits; we store 27.
63 Note that the decimal string conversion routines are sensitive to
64 rounding errors. Since the raw arithmetic routines do not themselves
65 have guard digits or rounding, the computation of 10**exp can
66 accumulate more than a few digits of error. The previous incarnation
67 of real.c successfully used a 144-bit fraction; given the current
68 layout of REAL_VALUE_TYPE we're forced to expand to at least 160 bits.
70 Target floating point models that use base 16 instead of base 2
71 (i.e. IBM 370), are handled during round_for_format, in which we
72 canonicalize the exponent to be a multiple of 4 (log2(16)), and
73 adjust the significand to match. */
76 /* Used to classify two numbers simultaneously. */
77 #define CLASS2(A, B) ((A) << 2 | (B))
79 #if HOST_BITS_PER_LONG != 64 && HOST_BITS_PER_LONG != 32
80 #error "Some constant folding done by hand to avoid shift count warnings"
83 static void get_zero (REAL_VALUE_TYPE *, int);
84 static void get_canonical_qnan (REAL_VALUE_TYPE *, int);
85 static void get_canonical_snan (REAL_VALUE_TYPE *, int);
86 static void get_inf (REAL_VALUE_TYPE *, int);
87 static bool sticky_rshift_significand (REAL_VALUE_TYPE *,
88 const REAL_VALUE_TYPE *, unsigned int);
89 static void rshift_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
91 static void lshift_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
93 static void lshift_significand_1 (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
94 static bool add_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *,
95 const REAL_VALUE_TYPE *);
96 static bool sub_significands (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
97 const REAL_VALUE_TYPE *, int);
98 static void neg_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
99 static int cmp_significands (const REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
100 static int cmp_significand_0 (const REAL_VALUE_TYPE *);
101 static void set_significand_bit (REAL_VALUE_TYPE *, unsigned int);
102 static void clear_significand_bit (REAL_VALUE_TYPE *, unsigned int);
103 static bool test_significand_bit (REAL_VALUE_TYPE *, unsigned int);
104 static void clear_significand_below (REAL_VALUE_TYPE *, unsigned int);
105 static bool div_significands (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
106 const REAL_VALUE_TYPE *);
107 static void normalize (REAL_VALUE_TYPE *);
109 static bool do_add (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
110 const REAL_VALUE_TYPE *, int);
111 static bool do_multiply (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
112 const REAL_VALUE_TYPE *);
113 static bool do_divide (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
114 const REAL_VALUE_TYPE *);
115 static int do_compare (const REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *, int);
116 static void do_fix_trunc (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
118 static unsigned long rtd_divmod (REAL_VALUE_TYPE *, REAL_VALUE_TYPE *);
120 static const REAL_VALUE_TYPE * ten_to_ptwo (int);
121 static const REAL_VALUE_TYPE * ten_to_mptwo (int);
122 static const REAL_VALUE_TYPE * real_digit (int);
123 static void times_pten (REAL_VALUE_TYPE *, int);
125 static void round_for_format (const struct real_format *, REAL_VALUE_TYPE *);
127 /* Initialize R with a positive zero. */
130 get_zero (REAL_VALUE_TYPE *r, int sign)
132 memset (r, 0, sizeof (*r));
136 /* Initialize R with the canonical quiet NaN. */
139 get_canonical_qnan (REAL_VALUE_TYPE *r, int sign)
141 memset (r, 0, sizeof (*r));
148 get_canonical_snan (REAL_VALUE_TYPE *r, int sign)
150 memset (r, 0, sizeof (*r));
158 get_inf (REAL_VALUE_TYPE *r, int sign)
160 memset (r, 0, sizeof (*r));
166 /* Right-shift the significand of A by N bits; put the result in the
167 significand of R. If any one bits are shifted out, return true. */
170 sticky_rshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
173 unsigned long sticky = 0;
174 unsigned int i, ofs = 0;
176 if (n >= HOST_BITS_PER_LONG)
178 for (i = 0, ofs = n / HOST_BITS_PER_LONG; i < ofs; ++i)
180 n &= HOST_BITS_PER_LONG - 1;
185 sticky |= a->sig[ofs] & (((unsigned long)1 << n) - 1);
186 for (i = 0; i < SIGSZ; ++i)
189 = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
190 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
191 << (HOST_BITS_PER_LONG - n)));
196 for (i = 0; ofs + i < SIGSZ; ++i)
197 r->sig[i] = a->sig[ofs + i];
198 for (; i < SIGSZ; ++i)
205 /* Right-shift the significand of A by N bits; put the result in the
209 rshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
212 unsigned int i, ofs = n / HOST_BITS_PER_LONG;
214 n &= HOST_BITS_PER_LONG - 1;
217 for (i = 0; i < SIGSZ; ++i)
220 = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
221 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
222 << (HOST_BITS_PER_LONG - n)));
227 for (i = 0; ofs + i < SIGSZ; ++i)
228 r->sig[i] = a->sig[ofs + i];
229 for (; i < SIGSZ; ++i)
234 /* Left-shift the significand of A by N bits; put the result in the
238 lshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
241 unsigned int i, ofs = n / HOST_BITS_PER_LONG;
243 n &= HOST_BITS_PER_LONG - 1;
246 for (i = 0; ofs + i < SIGSZ; ++i)
247 r->sig[SIGSZ-1-i] = a->sig[SIGSZ-1-i-ofs];
248 for (; i < SIGSZ; ++i)
249 r->sig[SIGSZ-1-i] = 0;
252 for (i = 0; i < SIGSZ; ++i)
255 = (((ofs + i >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs]) << n)
256 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs-1])
257 >> (HOST_BITS_PER_LONG - n)));
261 /* Likewise, but N is specialized to 1. */
264 lshift_significand_1 (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
268 for (i = SIGSZ - 1; i > 0; --i)
269 r->sig[i] = (a->sig[i] << 1) | (a->sig[i-1] >> (HOST_BITS_PER_LONG - 1));
270 r->sig[0] = a->sig[0] << 1;
273 /* Add the significands of A and B, placing the result in R. Return
274 true if there was carry out of the most significant word. */
277 add_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
278 const REAL_VALUE_TYPE *b)
283 for (i = 0; i < SIGSZ; ++i)
285 unsigned long ai = a->sig[i];
286 unsigned long ri = ai + b->sig[i];
302 /* Subtract the significands of A and B, placing the result in R. CARRY is
303 true if there's a borrow incoming to the least significant word.
304 Return true if there was borrow out of the most significant word. */
307 sub_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
308 const REAL_VALUE_TYPE *b, int carry)
312 for (i = 0; i < SIGSZ; ++i)
314 unsigned long ai = a->sig[i];
315 unsigned long ri = ai - b->sig[i];
331 /* Negate the significand A, placing the result in R. */
334 neg_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
339 for (i = 0; i < SIGSZ; ++i)
341 unsigned long ri, ai = a->sig[i];
360 /* Compare significands. Return tri-state vs zero. */
363 cmp_significands (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
367 for (i = SIGSZ - 1; i >= 0; --i)
369 unsigned long ai = a->sig[i];
370 unsigned long bi = b->sig[i];
381 /* Return true if A is nonzero. */
384 cmp_significand_0 (const REAL_VALUE_TYPE *a)
388 for (i = SIGSZ - 1; i >= 0; --i)
395 /* Set bit N of the significand of R. */
398 set_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
400 r->sig[n / HOST_BITS_PER_LONG]
401 |= (unsigned long)1 << (n % HOST_BITS_PER_LONG);
404 /* Clear bit N of the significand of R. */
407 clear_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
409 r->sig[n / HOST_BITS_PER_LONG]
410 &= ~((unsigned long)1 << (n % HOST_BITS_PER_LONG));
413 /* Test bit N of the significand of R. */
416 test_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
418 /* ??? Compiler bug here if we return this expression directly.
419 The conversion to bool strips the "&1" and we wind up testing
420 e.g. 2 != 0 -> true. Seen in gcc version 3.2 20020520. */
421 int t = (r->sig[n / HOST_BITS_PER_LONG] >> (n % HOST_BITS_PER_LONG)) & 1;
425 /* Clear bits 0..N-1 of the significand of R. */
428 clear_significand_below (REAL_VALUE_TYPE *r, unsigned int n)
430 int i, w = n / HOST_BITS_PER_LONG;
432 for (i = 0; i < w; ++i)
435 r->sig[w] &= ~(((unsigned long)1 << (n % HOST_BITS_PER_LONG)) - 1);
438 /* Divide the significands of A and B, placing the result in R. Return
439 true if the division was inexact. */
442 div_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
443 const REAL_VALUE_TYPE *b)
446 int i, bit = SIGNIFICAND_BITS - 1;
447 unsigned long msb, inexact;
450 memset (r->sig, 0, sizeof (r->sig));
456 msb = u.sig[SIGSZ-1] & SIG_MSB;
457 lshift_significand_1 (&u, &u);
459 if (msb || cmp_significands (&u, b) >= 0)
461 sub_significands (&u, &u, b, 0);
462 set_significand_bit (r, bit);
467 for (i = 0, inexact = 0; i < SIGSZ; i++)
473 /* Adjust the exponent and significand of R such that the most
474 significant bit is set. We underflow to zero and overflow to
475 infinity here, without denormals. (The intermediate representation
476 exponent is large enough to handle target denormals normalized.) */
479 normalize (REAL_VALUE_TYPE *r)
487 /* Find the first word that is nonzero. */
488 for (i = SIGSZ - 1; i >= 0; i--)
490 shift += HOST_BITS_PER_LONG;
494 /* Zero significand flushes to zero. */
502 /* Find the first bit that is nonzero. */
504 if (r->sig[i] & ((unsigned long)1 << (HOST_BITS_PER_LONG - 1 - j)))
510 exp = REAL_EXP (r) - shift;
512 get_inf (r, r->sign);
513 else if (exp < -MAX_EXP)
514 get_zero (r, r->sign);
517 SET_REAL_EXP (r, exp);
518 lshift_significand (r, r, shift);
523 /* Calculate R = A + (SUBTRACT_P ? -B : B). Return true if the
524 result may be inexact due to a loss of precision. */
527 do_add (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
528 const REAL_VALUE_TYPE *b, int subtract_p)
532 bool inexact = false;
534 /* Determine if we need to add or subtract. */
536 subtract_p = (sign ^ b->sign) ^ subtract_p;
538 switch (CLASS2 (a->cl, b->cl))
540 case CLASS2 (rvc_zero, rvc_zero):
541 /* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0. */
542 get_zero (r, sign & !subtract_p);
545 case CLASS2 (rvc_zero, rvc_normal):
546 case CLASS2 (rvc_zero, rvc_inf):
547 case CLASS2 (rvc_zero, rvc_nan):
549 case CLASS2 (rvc_normal, rvc_nan):
550 case CLASS2 (rvc_inf, rvc_nan):
551 case CLASS2 (rvc_nan, rvc_nan):
552 /* ANY + NaN = NaN. */
553 case CLASS2 (rvc_normal, rvc_inf):
556 r->sign = sign ^ subtract_p;
559 case CLASS2 (rvc_normal, rvc_zero):
560 case CLASS2 (rvc_inf, rvc_zero):
561 case CLASS2 (rvc_nan, rvc_zero):
563 case CLASS2 (rvc_nan, rvc_normal):
564 case CLASS2 (rvc_nan, rvc_inf):
565 /* NaN + ANY = NaN. */
566 case CLASS2 (rvc_inf, rvc_normal):
571 case CLASS2 (rvc_inf, rvc_inf):
573 /* Inf - Inf = NaN. */
574 get_canonical_qnan (r, 0);
576 /* Inf + Inf = Inf. */
580 case CLASS2 (rvc_normal, rvc_normal):
587 /* Swap the arguments such that A has the larger exponent. */
588 dexp = REAL_EXP (a) - REAL_EXP (b);
591 const REAL_VALUE_TYPE *t;
598 /* If the exponents are not identical, we need to shift the
599 significand of B down. */
602 /* If the exponents are too far apart, the significands
603 do not overlap, which makes the subtraction a noop. */
604 if (dexp >= SIGNIFICAND_BITS)
611 inexact |= sticky_rshift_significand (&t, b, dexp);
617 if (sub_significands (r, a, b, inexact))
619 /* We got a borrow out of the subtraction. That means that
620 A and B had the same exponent, and B had the larger
621 significand. We need to swap the sign and negate the
624 neg_significand (r, r);
629 if (add_significands (r, a, b))
631 /* We got carry out of the addition. This means we need to
632 shift the significand back down one bit and increase the
634 inexact |= sticky_rshift_significand (r, r, 1);
635 r->sig[SIGSZ-1] |= SIG_MSB;
646 SET_REAL_EXP (r, exp);
647 /* Zero out the remaining fields. */
652 /* Re-normalize the result. */
655 /* Special case: if the subtraction results in zero, the result
657 if (r->cl == rvc_zero)
660 r->sig[0] |= inexact;
665 /* Calculate R = A * B. Return true if the result may be inexact. */
668 do_multiply (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
669 const REAL_VALUE_TYPE *b)
671 REAL_VALUE_TYPE u, t, *rr;
672 unsigned int i, j, k;
673 int sign = a->sign ^ b->sign;
674 bool inexact = false;
676 switch (CLASS2 (a->cl, b->cl))
678 case CLASS2 (rvc_zero, rvc_zero):
679 case CLASS2 (rvc_zero, rvc_normal):
680 case CLASS2 (rvc_normal, rvc_zero):
681 /* +-0 * ANY = 0 with appropriate sign. */
685 case CLASS2 (rvc_zero, rvc_nan):
686 case CLASS2 (rvc_normal, rvc_nan):
687 case CLASS2 (rvc_inf, rvc_nan):
688 case CLASS2 (rvc_nan, rvc_nan):
689 /* ANY * NaN = NaN. */
694 case CLASS2 (rvc_nan, rvc_zero):
695 case CLASS2 (rvc_nan, rvc_normal):
696 case CLASS2 (rvc_nan, rvc_inf):
697 /* NaN * ANY = NaN. */
702 case CLASS2 (rvc_zero, rvc_inf):
703 case CLASS2 (rvc_inf, rvc_zero):
705 get_canonical_qnan (r, sign);
708 case CLASS2 (rvc_inf, rvc_inf):
709 case CLASS2 (rvc_normal, rvc_inf):
710 case CLASS2 (rvc_inf, rvc_normal):
711 /* Inf * Inf = Inf, R * Inf = Inf */
715 case CLASS2 (rvc_normal, rvc_normal):
722 if (r == a || r == b)
728 /* Collect all the partial products. Since we don't have sure access
729 to a widening multiply, we split each long into two half-words.
731 Consider the long-hand form of a four half-word multiplication:
741 We construct partial products of the widened half-word products
742 that are known to not overlap, e.g. DF+DH. Each such partial
743 product is given its proper exponent, which allows us to sum them
744 and obtain the finished product. */
746 for (i = 0; i < SIGSZ * 2; ++i)
748 unsigned long ai = a->sig[i / 2];
750 ai >>= HOST_BITS_PER_LONG / 2;
752 ai &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
757 for (j = 0; j < 2; ++j)
759 int exp = (REAL_EXP (a) - (2*SIGSZ-1-i)*(HOST_BITS_PER_LONG/2)
760 + (REAL_EXP (b) - (1-j)*(HOST_BITS_PER_LONG/2)));
769 /* Would underflow to zero, which we shouldn't bother adding. */
774 memset (&u, 0, sizeof (u));
776 SET_REAL_EXP (&u, exp);
778 for (k = j; k < SIGSZ * 2; k += 2)
780 unsigned long bi = b->sig[k / 2];
782 bi >>= HOST_BITS_PER_LONG / 2;
784 bi &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
786 u.sig[k / 2] = ai * bi;
790 inexact |= do_add (rr, rr, &u, 0);
801 /* Calculate R = A / B. Return true if the result may be inexact. */
804 do_divide (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
805 const REAL_VALUE_TYPE *b)
807 int exp, sign = a->sign ^ b->sign;
808 REAL_VALUE_TYPE t, *rr;
811 switch (CLASS2 (a->cl, b->cl))
813 case CLASS2 (rvc_zero, rvc_zero):
815 case CLASS2 (rvc_inf, rvc_inf):
816 /* Inf / Inf = NaN. */
817 get_canonical_qnan (r, sign);
820 case CLASS2 (rvc_zero, rvc_normal):
821 case CLASS2 (rvc_zero, rvc_inf):
823 case CLASS2 (rvc_normal, rvc_inf):
828 case CLASS2 (rvc_normal, rvc_zero):
830 case CLASS2 (rvc_inf, rvc_zero):
835 case CLASS2 (rvc_zero, rvc_nan):
836 case CLASS2 (rvc_normal, rvc_nan):
837 case CLASS2 (rvc_inf, rvc_nan):
838 case CLASS2 (rvc_nan, rvc_nan):
839 /* ANY / NaN = NaN. */
844 case CLASS2 (rvc_nan, rvc_zero):
845 case CLASS2 (rvc_nan, rvc_normal):
846 case CLASS2 (rvc_nan, rvc_inf):
847 /* NaN / ANY = NaN. */
852 case CLASS2 (rvc_inf, rvc_normal):
857 case CLASS2 (rvc_normal, rvc_normal):
864 if (r == a || r == b)
869 /* Make sure all fields in the result are initialized. */
874 exp = REAL_EXP (a) - REAL_EXP (b) + 1;
885 SET_REAL_EXP (rr, exp);
887 inexact = div_significands (rr, a, b);
889 /* Re-normalize the result. */
891 rr->sig[0] |= inexact;
899 /* Return a tri-state comparison of A vs B. Return NAN_RESULT if
900 one of the two operands is a NaN. */
903 do_compare (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b,
908 switch (CLASS2 (a->cl, b->cl))
910 case CLASS2 (rvc_zero, rvc_zero):
911 /* Sign of zero doesn't matter for compares. */
914 case CLASS2 (rvc_inf, rvc_zero):
915 case CLASS2 (rvc_inf, rvc_normal):
916 case CLASS2 (rvc_normal, rvc_zero):
917 return (a->sign ? -1 : 1);
919 case CLASS2 (rvc_inf, rvc_inf):
920 return -a->sign - -b->sign;
922 case CLASS2 (rvc_zero, rvc_normal):
923 case CLASS2 (rvc_zero, rvc_inf):
924 case CLASS2 (rvc_normal, rvc_inf):
925 return (b->sign ? 1 : -1);
927 case CLASS2 (rvc_zero, rvc_nan):
928 case CLASS2 (rvc_normal, rvc_nan):
929 case CLASS2 (rvc_inf, rvc_nan):
930 case CLASS2 (rvc_nan, rvc_nan):
931 case CLASS2 (rvc_nan, rvc_zero):
932 case CLASS2 (rvc_nan, rvc_normal):
933 case CLASS2 (rvc_nan, rvc_inf):
936 case CLASS2 (rvc_normal, rvc_normal):
943 if (a->sign != b->sign)
944 return -a->sign - -b->sign;
946 if (a->decimal || b->decimal)
947 return decimal_do_compare (a, b, nan_result);
949 if (REAL_EXP (a) > REAL_EXP (b))
951 else if (REAL_EXP (a) < REAL_EXP (b))
954 ret = cmp_significands (a, b);
956 return (a->sign ? -ret : ret);
959 /* Return A truncated to an integral value toward zero. */
962 do_fix_trunc (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
976 decimal_do_fix_trunc (r, a);
979 if (REAL_EXP (r) <= 0)
980 get_zero (r, r->sign);
981 else if (REAL_EXP (r) < SIGNIFICAND_BITS)
982 clear_significand_below (r, SIGNIFICAND_BITS - REAL_EXP (r));
990 /* Perform the binary or unary operation described by CODE.
991 For a unary operation, leave OP1 NULL. This function returns
992 true if the result may be inexact due to loss of precision. */
995 real_arithmetic (REAL_VALUE_TYPE *r, int icode, const REAL_VALUE_TYPE *op0,
996 const REAL_VALUE_TYPE *op1)
998 enum tree_code code = icode;
1000 if (op0->decimal || (op1 && op1->decimal))
1001 return decimal_real_arithmetic (r, icode, op0, op1);
1006 return do_add (r, op0, op1, 0);
1009 return do_add (r, op0, op1, 1);
1012 return do_multiply (r, op0, op1);
1015 return do_divide (r, op0, op1);
1018 if (op1->cl == rvc_nan)
1020 else if (do_compare (op0, op1, -1) < 0)
1027 if (op1->cl == rvc_nan)
1029 else if (do_compare (op0, op1, 1) < 0)
1045 case FIX_TRUNC_EXPR:
1046 do_fix_trunc (r, op0);
1055 /* Legacy. Similar, but return the result directly. */
1058 real_arithmetic2 (int icode, const REAL_VALUE_TYPE *op0,
1059 const REAL_VALUE_TYPE *op1)
1062 real_arithmetic (&r, icode, op0, op1);
1067 real_compare (int icode, const REAL_VALUE_TYPE *op0,
1068 const REAL_VALUE_TYPE *op1)
1070 enum tree_code code = icode;
1075 return do_compare (op0, op1, 1) < 0;
1077 return do_compare (op0, op1, 1) <= 0;
1079 return do_compare (op0, op1, -1) > 0;
1081 return do_compare (op0, op1, -1) >= 0;
1083 return do_compare (op0, op1, -1) == 0;
1085 return do_compare (op0, op1, -1) != 0;
1086 case UNORDERED_EXPR:
1087 return op0->cl == rvc_nan || op1->cl == rvc_nan;
1089 return op0->cl != rvc_nan && op1->cl != rvc_nan;
1091 return do_compare (op0, op1, -1) < 0;
1093 return do_compare (op0, op1, -1) <= 0;
1095 return do_compare (op0, op1, 1) > 0;
1097 return do_compare (op0, op1, 1) >= 0;
1099 return do_compare (op0, op1, 0) == 0;
1101 return do_compare (op0, op1, 0) != 0;
1108 /* Return floor log2(R). */
1111 real_exponent (const REAL_VALUE_TYPE *r)
1119 return (unsigned int)-1 >> 1;
1121 return REAL_EXP (r);
1127 /* R = OP0 * 2**EXP. */
1130 real_ldexp (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *op0, int exp)
1141 exp += REAL_EXP (op0);
1143 get_inf (r, r->sign);
1144 else if (exp < -MAX_EXP)
1145 get_zero (r, r->sign);
1147 SET_REAL_EXP (r, exp);
1155 /* Determine whether a floating-point value X is infinite. */
1158 real_isinf (const REAL_VALUE_TYPE *r)
1160 return (r->cl == rvc_inf);
1163 /* Determine whether a floating-point value X is a NaN. */
1166 real_isnan (const REAL_VALUE_TYPE *r)
1168 return (r->cl == rvc_nan);
1171 /* Determine whether a floating-point value X is negative. */
1174 real_isneg (const REAL_VALUE_TYPE *r)
1179 /* Determine whether a floating-point value X is minus zero. */
1182 real_isnegzero (const REAL_VALUE_TYPE *r)
1184 return r->sign && r->cl == rvc_zero;
1187 /* Compare two floating-point objects for bitwise identity. */
1190 real_identical (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
1196 if (a->sign != b->sign)
1206 if (a->decimal != b->decimal)
1208 if (REAL_EXP (a) != REAL_EXP (b))
1213 if (a->signalling != b->signalling)
1215 /* The significand is ignored for canonical NaNs. */
1216 if (a->canonical || b->canonical)
1217 return a->canonical == b->canonical;
1224 for (i = 0; i < SIGSZ; ++i)
1225 if (a->sig[i] != b->sig[i])
1231 /* Try to change R into its exact multiplicative inverse in machine
1232 mode MODE. Return true if successful. */
1235 exact_real_inverse (enum machine_mode mode, REAL_VALUE_TYPE *r)
1237 const REAL_VALUE_TYPE *one = real_digit (1);
1241 if (r->cl != rvc_normal)
1244 /* Check for a power of two: all significand bits zero except the MSB. */
1245 for (i = 0; i < SIGSZ-1; ++i)
1248 if (r->sig[SIGSZ-1] != SIG_MSB)
1251 /* Find the inverse and truncate to the required mode. */
1252 do_divide (&u, one, r);
1253 real_convert (&u, mode, &u);
1255 /* The rounding may have overflowed. */
1256 if (u.cl != rvc_normal)
1258 for (i = 0; i < SIGSZ-1; ++i)
1261 if (u.sig[SIGSZ-1] != SIG_MSB)
1268 /* Render R as an integer. */
1271 real_to_integer (const REAL_VALUE_TYPE *r)
1273 unsigned HOST_WIDE_INT i;
1284 i = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1291 return decimal_real_to_integer (r);
1293 if (REAL_EXP (r) <= 0)
1295 /* Only force overflow for unsigned overflow. Signed overflow is
1296 undefined, so it doesn't matter what we return, and some callers
1297 expect to be able to use this routine for both signed and
1298 unsigned conversions. */
1299 if (REAL_EXP (r) > HOST_BITS_PER_WIDE_INT)
1302 if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1303 i = r->sig[SIGSZ-1];
1306 gcc_assert (HOST_BITS_PER_WIDE_INT == 2 * HOST_BITS_PER_LONG);
1307 i = r->sig[SIGSZ-1];
1308 i = i << (HOST_BITS_PER_LONG - 1) << 1;
1309 i |= r->sig[SIGSZ-2];
1312 i >>= HOST_BITS_PER_WIDE_INT - REAL_EXP (r);
1323 /* Likewise, but to an integer pair, HI+LOW. */
1326 real_to_integer2 (HOST_WIDE_INT *plow, HOST_WIDE_INT *phigh,
1327 const REAL_VALUE_TYPE *r)
1330 HOST_WIDE_INT low, high;
1343 high = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1356 decimal_real_to_integer2 (plow, phigh, r);
1363 /* Only force overflow for unsigned overflow. Signed overflow is
1364 undefined, so it doesn't matter what we return, and some callers
1365 expect to be able to use this routine for both signed and
1366 unsigned conversions. */
1367 if (exp > 2*HOST_BITS_PER_WIDE_INT)
1370 rshift_significand (&t, r, 2*HOST_BITS_PER_WIDE_INT - exp);
1371 if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1373 high = t.sig[SIGSZ-1];
1374 low = t.sig[SIGSZ-2];
1378 gcc_assert (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG);
1379 high = t.sig[SIGSZ-1];
1380 high = high << (HOST_BITS_PER_LONG - 1) << 1;
1381 high |= t.sig[SIGSZ-2];
1383 low = t.sig[SIGSZ-3];
1384 low = low << (HOST_BITS_PER_LONG - 1) << 1;
1385 low |= t.sig[SIGSZ-4];
1393 low = -low, high = ~high;
1405 /* A subroutine of real_to_decimal. Compute the quotient and remainder
1406 of NUM / DEN. Return the quotient and place the remainder in NUM.
1407 It is expected that NUM / DEN are close enough that the quotient is
1410 static unsigned long
1411 rtd_divmod (REAL_VALUE_TYPE *num, REAL_VALUE_TYPE *den)
1413 unsigned long q, msb;
1414 int expn = REAL_EXP (num), expd = REAL_EXP (den);
1423 msb = num->sig[SIGSZ-1] & SIG_MSB;
1425 lshift_significand_1 (num, num);
1427 if (msb || cmp_significands (num, den) >= 0)
1429 sub_significands (num, num, den, 0);
1433 while (--expn >= expd);
1435 SET_REAL_EXP (num, expd);
1441 /* Render R as a decimal floating point constant. Emit DIGITS significant
1442 digits in the result, bounded by BUF_SIZE. If DIGITS is 0, choose the
1443 maximum for the representation. If CROP_TRAILING_ZEROS, strip trailing
1446 #define M_LOG10_2 0.30102999566398119521
1449 real_to_decimal (char *str, const REAL_VALUE_TYPE *r_orig, size_t buf_size,
1450 size_t digits, int crop_trailing_zeros)
1452 const REAL_VALUE_TYPE *one, *ten;
1453 REAL_VALUE_TYPE r, pten, u, v;
1454 int dec_exp, cmp_one, digit;
1456 char *p, *first, *last;
1463 strcpy (str, (r.sign ? "-0.0" : "0.0"));
1468 strcpy (str, (r.sign ? "-Inf" : "+Inf"));
1471 /* ??? Print the significand as well, if not canonical? */
1472 strcpy (str, (r.sign ? "-NaN" : "+NaN"));
1480 decimal_real_to_decimal (str, &r, buf_size, digits, crop_trailing_zeros);
1484 /* Bound the number of digits printed by the size of the representation. */
1485 max_digits = SIGNIFICAND_BITS * M_LOG10_2;
1486 if (digits == 0 || digits > max_digits)
1487 digits = max_digits;
1489 /* Estimate the decimal exponent, and compute the length of the string it
1490 will print as. Be conservative and add one to account for possible
1491 overflow or rounding error. */
1492 dec_exp = REAL_EXP (&r) * M_LOG10_2;
1493 for (max_digits = 1; dec_exp ; max_digits++)
1496 /* Bound the number of digits printed by the size of the output buffer. */
1497 max_digits = buf_size - 1 - 1 - 2 - max_digits - 1;
1498 gcc_assert (max_digits <= buf_size);
1499 if (digits > max_digits)
1500 digits = max_digits;
1502 one = real_digit (1);
1503 ten = ten_to_ptwo (0);
1511 cmp_one = do_compare (&r, one, 0);
1516 /* Number is greater than one. Convert significand to an integer
1517 and strip trailing decimal zeros. */
1520 SET_REAL_EXP (&u, SIGNIFICAND_BITS - 1);
1522 /* Largest M, such that 10**2**M fits within SIGNIFICAND_BITS. */
1523 m = floor_log2 (max_digits);
1525 /* Iterate over the bits of the possible powers of 10 that might
1526 be present in U and eliminate them. That is, if we find that
1527 10**2**M divides U evenly, keep the division and increase
1533 do_divide (&t, &u, ten_to_ptwo (m));
1534 do_fix_trunc (&v, &t);
1535 if (cmp_significands (&v, &t) == 0)
1543 /* Revert the scaling to integer that we performed earlier. */
1544 SET_REAL_EXP (&u, REAL_EXP (&u) + REAL_EXP (&r)
1545 - (SIGNIFICAND_BITS - 1));
1548 /* Find power of 10. Do this by dividing out 10**2**M when
1549 this is larger than the current remainder. Fill PTEN with
1550 the power of 10 that we compute. */
1551 if (REAL_EXP (&r) > 0)
1553 m = floor_log2 ((int)(REAL_EXP (&r) * M_LOG10_2)) + 1;
1556 const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1557 if (do_compare (&u, ptentwo, 0) >= 0)
1559 do_divide (&u, &u, ptentwo);
1560 do_multiply (&pten, &pten, ptentwo);
1567 /* We managed to divide off enough tens in the above reduction
1568 loop that we've now got a negative exponent. Fall into the
1569 less-than-one code to compute the proper value for PTEN. */
1576 /* Number is less than one. Pad significand with leading
1582 /* Stop if we'd shift bits off the bottom. */
1586 do_multiply (&u, &v, ten);
1588 /* Stop if we're now >= 1. */
1589 if (REAL_EXP (&u) > 0)
1597 /* Find power of 10. Do this by multiplying in P=10**2**M when
1598 the current remainder is smaller than 1/P. Fill PTEN with the
1599 power of 10 that we compute. */
1600 m = floor_log2 ((int)(-REAL_EXP (&r) * M_LOG10_2)) + 1;
1603 const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1604 const REAL_VALUE_TYPE *ptenmtwo = ten_to_mptwo (m);
1606 if (do_compare (&v, ptenmtwo, 0) <= 0)
1608 do_multiply (&v, &v, ptentwo);
1609 do_multiply (&pten, &pten, ptentwo);
1615 /* Invert the positive power of 10 that we've collected so far. */
1616 do_divide (&pten, one, &pten);
1624 /* At this point, PTEN should contain the nearest power of 10 smaller
1625 than R, such that this division produces the first digit.
1627 Using a divide-step primitive that returns the complete integral
1628 remainder avoids the rounding error that would be produced if
1629 we were to use do_divide here and then simply multiply by 10 for
1630 each subsequent digit. */
1632 digit = rtd_divmod (&r, &pten);
1634 /* Be prepared for error in that division via underflow ... */
1635 if (digit == 0 && cmp_significand_0 (&r))
1637 /* Multiply by 10 and try again. */
1638 do_multiply (&r, &r, ten);
1639 digit = rtd_divmod (&r, &pten);
1641 gcc_assert (digit != 0);
1644 /* ... or overflow. */
1654 gcc_assert (digit <= 10);
1658 /* Generate subsequent digits. */
1659 while (--digits > 0)
1661 do_multiply (&r, &r, ten);
1662 digit = rtd_divmod (&r, &pten);
1667 /* Generate one more digit with which to do rounding. */
1668 do_multiply (&r, &r, ten);
1669 digit = rtd_divmod (&r, &pten);
1671 /* Round the result. */
1674 /* Round to nearest. If R is nonzero there are additional
1675 nonzero digits to be extracted. */
1676 if (cmp_significand_0 (&r))
1678 /* Round to even. */
1679 else if ((p[-1] - '0') & 1)
1696 /* Carry out of the first digit. This means we had all 9's and
1697 now have all 0's. "Prepend" a 1 by overwriting the first 0. */
1705 /* Insert the decimal point. */
1706 first[0] = first[1];
1709 /* If requested, drop trailing zeros. Never crop past "1.0". */
1710 if (crop_trailing_zeros)
1711 while (last > first + 3 && last[-1] == '0')
1714 /* Append the exponent. */
1715 sprintf (last, "e%+d", dec_exp);
1718 /* Render R as a hexadecimal floating point constant. Emit DIGITS
1719 significant digits in the result, bounded by BUF_SIZE. If DIGITS is 0,
1720 choose the maximum for the representation. If CROP_TRAILING_ZEROS,
1721 strip trailing zeros. */
1724 real_to_hexadecimal (char *str, const REAL_VALUE_TYPE *r, size_t buf_size,
1725 size_t digits, int crop_trailing_zeros)
1727 int i, j, exp = REAL_EXP (r);
1740 strcpy (str, (r->sign ? "-Inf" : "+Inf"));
1743 /* ??? Print the significand as well, if not canonical? */
1744 strcpy (str, (r->sign ? "-NaN" : "+NaN"));
1752 /* Hexadecimal format for decimal floats is not interesting. */
1753 strcpy (str, "N/A");
1758 digits = SIGNIFICAND_BITS / 4;
1760 /* Bound the number of digits printed by the size of the output buffer. */
1762 sprintf (exp_buf, "p%+d", exp);
1763 max_digits = buf_size - strlen (exp_buf) - r->sign - 4 - 1;
1764 gcc_assert (max_digits <= buf_size);
1765 if (digits > max_digits)
1766 digits = max_digits;
1777 for (i = SIGSZ - 1; i >= 0; --i)
1778 for (j = HOST_BITS_PER_LONG - 4; j >= 0; j -= 4)
1780 *p++ = "0123456789abcdef"[(r->sig[i] >> j) & 15];
1786 if (crop_trailing_zeros)
1787 while (p > first + 1 && p[-1] == '0')
1790 sprintf (p, "p%+d", exp);
1793 /* Initialize R from a decimal or hexadecimal string. The string is
1794 assumed to have been syntax checked already. */
1797 real_from_string (REAL_VALUE_TYPE *r, const char *str)
1809 else if (*str == '+')
1812 if (str[0] == '0' && (str[1] == 'x' || str[1] == 'X'))
1814 /* Hexadecimal floating point. */
1815 int pos = SIGNIFICAND_BITS - 4, d;
1823 d = hex_value (*str);
1828 r->sig[pos / HOST_BITS_PER_LONG]
1829 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1833 /* Ensure correct rounding by setting last bit if there is
1834 a subsequent nonzero digit. */
1842 if (pos == SIGNIFICAND_BITS - 4)
1849 d = hex_value (*str);
1854 r->sig[pos / HOST_BITS_PER_LONG]
1855 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1859 /* Ensure correct rounding by setting last bit if there is
1860 a subsequent nonzero digit. */
1865 if (*str == 'p' || *str == 'P')
1867 bool exp_neg = false;
1875 else if (*str == '+')
1879 while (ISDIGIT (*str))
1885 /* Overflowed the exponent. */
1900 SET_REAL_EXP (r, exp);
1906 /* Decimal floating point. */
1907 const REAL_VALUE_TYPE *ten = ten_to_ptwo (0);
1912 while (ISDIGIT (*str))
1915 do_multiply (r, r, ten);
1917 do_add (r, r, real_digit (d), 0);
1922 if (r->cl == rvc_zero)
1927 while (ISDIGIT (*str))
1930 do_multiply (r, r, ten);
1932 do_add (r, r, real_digit (d), 0);
1937 if (*str == 'e' || *str == 'E')
1939 bool exp_neg = false;
1947 else if (*str == '+')
1951 while (ISDIGIT (*str))
1957 /* Overflowed the exponent. */
1971 times_pten (r, exp);
1986 /* Legacy. Similar, but return the result directly. */
1989 real_from_string2 (const char *s, enum machine_mode mode)
1993 real_from_string (&r, s);
1994 if (mode != VOIDmode)
1995 real_convert (&r, mode, &r);
2000 /* Initialize R from string S and desired MODE. */
2003 real_from_string3 (REAL_VALUE_TYPE *r, const char *s, enum machine_mode mode)
2005 if (DECIMAL_FLOAT_MODE_P (mode))
2006 decimal_real_from_string (r, s);
2008 real_from_string (r, s);
2010 if (mode != VOIDmode)
2011 real_convert (r, mode, r);
2014 /* Initialize R from the integer pair HIGH+LOW. */
2017 real_from_integer (REAL_VALUE_TYPE *r, enum machine_mode mode,
2018 unsigned HOST_WIDE_INT low, HOST_WIDE_INT high,
2021 if (low == 0 && high == 0)
2025 memset (r, 0, sizeof (*r));
2027 r->sign = high < 0 && !unsigned_p;
2028 SET_REAL_EXP (r, 2 * HOST_BITS_PER_WIDE_INT);
2039 if (HOST_BITS_PER_LONG == HOST_BITS_PER_WIDE_INT)
2041 r->sig[SIGSZ-1] = high;
2042 r->sig[SIGSZ-2] = low;
2046 gcc_assert (HOST_BITS_PER_LONG*2 == HOST_BITS_PER_WIDE_INT);
2047 r->sig[SIGSZ-1] = high >> (HOST_BITS_PER_LONG - 1) >> 1;
2048 r->sig[SIGSZ-2] = high;
2049 r->sig[SIGSZ-3] = low >> (HOST_BITS_PER_LONG - 1) >> 1;
2050 r->sig[SIGSZ-4] = low;
2056 if (mode != VOIDmode)
2057 real_convert (r, mode, r);
2060 /* Returns 10**2**N. */
2062 static const REAL_VALUE_TYPE *
2065 static REAL_VALUE_TYPE tens[EXP_BITS];
2067 gcc_assert (n >= 0);
2068 gcc_assert (n < EXP_BITS);
2070 if (tens[n].cl == rvc_zero)
2072 if (n < (HOST_BITS_PER_WIDE_INT == 64 ? 5 : 4))
2074 HOST_WIDE_INT t = 10;
2077 for (i = 0; i < n; ++i)
2080 real_from_integer (&tens[n], VOIDmode, t, 0, 1);
2084 const REAL_VALUE_TYPE *t = ten_to_ptwo (n - 1);
2085 do_multiply (&tens[n], t, t);
2092 /* Returns 10**(-2**N). */
2094 static const REAL_VALUE_TYPE *
2095 ten_to_mptwo (int n)
2097 static REAL_VALUE_TYPE tens[EXP_BITS];
2099 gcc_assert (n >= 0);
2100 gcc_assert (n < EXP_BITS);
2102 if (tens[n].cl == rvc_zero)
2103 do_divide (&tens[n], real_digit (1), ten_to_ptwo (n));
2110 static const REAL_VALUE_TYPE *
2113 static REAL_VALUE_TYPE num[10];
2115 gcc_assert (n >= 0);
2116 gcc_assert (n <= 9);
2118 if (n > 0 && num[n].cl == rvc_zero)
2119 real_from_integer (&num[n], VOIDmode, n, 0, 1);
2124 /* Multiply R by 10**EXP. */
2127 times_pten (REAL_VALUE_TYPE *r, int exp)
2129 REAL_VALUE_TYPE pten, *rr;
2130 bool negative = (exp < 0);
2136 pten = *real_digit (1);
2142 for (i = 0; exp > 0; ++i, exp >>= 1)
2144 do_multiply (rr, rr, ten_to_ptwo (i));
2147 do_divide (r, r, &pten);
2150 /* Fills R with +Inf. */
2153 real_inf (REAL_VALUE_TYPE *r)
2158 /* Fills R with a NaN whose significand is described by STR. If QUIET,
2159 we force a QNaN, else we force an SNaN. The string, if not empty,
2160 is parsed as a number and placed in the significand. Return true
2161 if the string was successfully parsed. */
2164 real_nan (REAL_VALUE_TYPE *r, const char *str, int quiet,
2165 enum machine_mode mode)
2167 const struct real_format *fmt;
2169 fmt = REAL_MODE_FORMAT (mode);
2175 get_canonical_qnan (r, 0);
2177 get_canonical_snan (r, 0);
2183 memset (r, 0, sizeof (*r));
2186 /* Parse akin to strtol into the significand of R. */
2188 while (ISSPACE (*str))
2192 else if (*str == '+')
2197 if (*str == 'x' || *str == 'X')
2206 while ((d = hex_value (*str)) < base)
2213 lshift_significand (r, r, 3);
2216 lshift_significand (r, r, 4);
2219 lshift_significand_1 (&u, r);
2220 lshift_significand (r, r, 3);
2221 add_significands (r, r, &u);
2229 add_significands (r, r, &u);
2234 /* Must have consumed the entire string for success. */
2238 /* Shift the significand into place such that the bits
2239 are in the most significant bits for the format. */
2240 lshift_significand (r, r, SIGNIFICAND_BITS - fmt->pnan);
2242 /* Our MSB is always unset for NaNs. */
2243 r->sig[SIGSZ-1] &= ~SIG_MSB;
2245 /* Force quiet or signalling NaN. */
2246 r->signalling = !quiet;
2252 /* Fills R with the largest finite value representable in mode MODE.
2253 If SIGN is nonzero, R is set to the most negative finite value. */
2256 real_maxval (REAL_VALUE_TYPE *r, int sign, enum machine_mode mode)
2258 const struct real_format *fmt;
2261 fmt = REAL_MODE_FORMAT (mode);
2263 memset (r, 0, sizeof (*r));
2266 decimal_real_maxval (r, sign, mode);
2271 SET_REAL_EXP (r, fmt->emax * fmt->log2_b);
2273 np2 = SIGNIFICAND_BITS - fmt->p * fmt->log2_b;
2274 memset (r->sig, -1, SIGSZ * sizeof (unsigned long));
2275 clear_significand_below (r, np2);
2279 /* Fills R with 2**N. */
2282 real_2expN (REAL_VALUE_TYPE *r, int n)
2284 memset (r, 0, sizeof (*r));
2289 else if (n < -MAX_EXP)
2294 SET_REAL_EXP (r, n);
2295 r->sig[SIGSZ-1] = SIG_MSB;
2301 round_for_format (const struct real_format *fmt, REAL_VALUE_TYPE *r)
2304 unsigned long sticky;
2312 decimal_round_for_format (fmt, r);
2315 /* FIXME. We can come here via fp_easy_constant
2316 (e.g. -O0 on '_Decimal32 x = 1.0 + 2.0dd'), but have not
2317 investigated whether this convert needs to be here, or
2318 something else is missing. */
2319 decimal_real_convert (r, DFmode, r);
2322 p2 = fmt->p * fmt->log2_b;
2323 emin2m1 = (fmt->emin - 1) * fmt->log2_b;
2324 emax2 = fmt->emax * fmt->log2_b;
2326 np2 = SIGNIFICAND_BITS - p2;
2330 get_zero (r, r->sign);
2332 if (!fmt->has_signed_zero)
2337 get_inf (r, r->sign);
2342 clear_significand_below (r, np2);
2352 /* If we're not base2, normalize the exponent to a multiple of
2354 if (fmt->log2_b != 1)
2358 gcc_assert (fmt->b != 10);
2359 shift = REAL_EXP (r) & (fmt->log2_b - 1);
2362 shift = fmt->log2_b - shift;
2363 r->sig[0] |= sticky_rshift_significand (r, r, shift);
2364 SET_REAL_EXP (r, REAL_EXP (r) + shift);
2368 /* Check the range of the exponent. If we're out of range,
2369 either underflow or overflow. */
2370 if (REAL_EXP (r) > emax2)
2372 else if (REAL_EXP (r) <= emin2m1)
2376 if (!fmt->has_denorm)
2378 /* Don't underflow completely until we've had a chance to round. */
2379 if (REAL_EXP (r) < emin2m1)
2384 diff = emin2m1 - REAL_EXP (r) + 1;
2388 /* De-normalize the significand. */
2389 r->sig[0] |= sticky_rshift_significand (r, r, diff);
2390 SET_REAL_EXP (r, REAL_EXP (r) + diff);
2394 /* There are P2 true significand bits, followed by one guard bit,
2395 followed by one sticky bit, followed by stuff. Fold nonzero
2396 stuff into the sticky bit. */
2399 for (i = 0, w = (np2 - 1) / HOST_BITS_PER_LONG; i < w; ++i)
2400 sticky |= r->sig[i];
2402 r->sig[w] & (((unsigned long)1 << ((np2 - 1) % HOST_BITS_PER_LONG)) - 1);
2404 guard = test_significand_bit (r, np2 - 1);
2405 lsb = test_significand_bit (r, np2);
2407 /* Round to even. */
2408 if (guard && (sticky || lsb))
2412 set_significand_bit (&u, np2);
2414 if (add_significands (r, r, &u))
2416 /* Overflow. Means the significand had been all ones, and
2417 is now all zeros. Need to increase the exponent, and
2418 possibly re-normalize it. */
2419 SET_REAL_EXP (r, REAL_EXP (r) + 1);
2420 if (REAL_EXP (r) > emax2)
2422 r->sig[SIGSZ-1] = SIG_MSB;
2424 if (fmt->log2_b != 1)
2426 int shift = REAL_EXP (r) & (fmt->log2_b - 1);
2429 shift = fmt->log2_b - shift;
2430 rshift_significand (r, r, shift);
2431 SET_REAL_EXP (r, REAL_EXP (r) + shift);
2432 if (REAL_EXP (r) > emax2)
2439 /* Catch underflow that we deferred until after rounding. */
2440 if (REAL_EXP (r) <= emin2m1)
2443 /* Clear out trailing garbage. */
2444 clear_significand_below (r, np2);
2447 /* Extend or truncate to a new mode. */
2450 real_convert (REAL_VALUE_TYPE *r, enum machine_mode mode,
2451 const REAL_VALUE_TYPE *a)
2453 const struct real_format *fmt;
2455 fmt = REAL_MODE_FORMAT (mode);
2460 if (a->decimal || fmt->b == 10)
2461 decimal_real_convert (r, mode, a);
2463 round_for_format (fmt, r);
2465 /* round_for_format de-normalizes denormals. Undo just that part. */
2466 if (r->cl == rvc_normal)
2470 /* Legacy. Likewise, except return the struct directly. */
2473 real_value_truncate (enum machine_mode mode, REAL_VALUE_TYPE a)
2476 real_convert (&r, mode, &a);
2480 /* Return true if truncating to MODE is exact. */
2483 exact_real_truncate (enum machine_mode mode, const REAL_VALUE_TYPE *a)
2485 const struct real_format *fmt;
2489 fmt = REAL_MODE_FORMAT (mode);
2492 /* Don't allow conversion to denormals. */
2493 emin2m1 = (fmt->emin - 1) * fmt->log2_b;
2494 if (REAL_EXP (a) <= emin2m1)
2497 /* After conversion to the new mode, the value must be identical. */
2498 real_convert (&t, mode, a);
2499 return real_identical (&t, a);
2502 /* Write R to the given target format. Place the words of the result
2503 in target word order in BUF. There are always 32 bits in each
2504 long, no matter the size of the host long.
2506 Legacy: return word 0 for implementing REAL_VALUE_TO_TARGET_SINGLE. */
2509 real_to_target_fmt (long *buf, const REAL_VALUE_TYPE *r_orig,
2510 const struct real_format *fmt)
2516 round_for_format (fmt, &r);
2520 (*fmt->encode) (fmt, buf, &r);
2525 /* Similar, but look up the format from MODE. */
2528 real_to_target (long *buf, const REAL_VALUE_TYPE *r, enum machine_mode mode)
2530 const struct real_format *fmt;
2532 fmt = REAL_MODE_FORMAT (mode);
2535 return real_to_target_fmt (buf, r, fmt);
2538 /* Read R from the given target format. Read the words of the result
2539 in target word order in BUF. There are always 32 bits in each
2540 long, no matter the size of the host long. */
2543 real_from_target_fmt (REAL_VALUE_TYPE *r, const long *buf,
2544 const struct real_format *fmt)
2546 (*fmt->decode) (fmt, r, buf);
2549 /* Similar, but look up the format from MODE. */
2552 real_from_target (REAL_VALUE_TYPE *r, const long *buf, enum machine_mode mode)
2554 const struct real_format *fmt;
2556 fmt = REAL_MODE_FORMAT (mode);
2559 (*fmt->decode) (fmt, r, buf);
2562 /* Return the number of bits of the largest binary value that the
2563 significand of MODE will hold. */
2564 /* ??? Legacy. Should get access to real_format directly. */
2567 significand_size (enum machine_mode mode)
2569 const struct real_format *fmt;
2571 fmt = REAL_MODE_FORMAT (mode);
2577 /* Return the size in bits of the largest binary value that can be
2578 held by the decimal coefficient for this mode. This is one more
2579 than the number of bits required to hold the largest coefficient
2581 double log2_10 = 3.3219281;
2582 return fmt->p * log2_10;
2584 return fmt->p * fmt->log2_b;
2587 /* Return a hash value for the given real value. */
2588 /* ??? The "unsigned int" return value is intended to be hashval_t,
2589 but I didn't want to pull hashtab.h into real.h. */
2592 real_hash (const REAL_VALUE_TYPE *r)
2597 h = r->cl | (r->sign << 2);
2605 h |= REAL_EXP (r) << 3;
2610 h ^= (unsigned int)-1;
2619 if (sizeof(unsigned long) > sizeof(unsigned int))
2620 for (i = 0; i < SIGSZ; ++i)
2622 unsigned long s = r->sig[i];
2623 h ^= s ^ (s >> (HOST_BITS_PER_LONG / 2));
2626 for (i = 0; i < SIGSZ; ++i)
2632 /* IEEE single-precision format. */
2634 static void encode_ieee_single (const struct real_format *fmt,
2635 long *, const REAL_VALUE_TYPE *);
2636 static void decode_ieee_single (const struct real_format *,
2637 REAL_VALUE_TYPE *, const long *);
2640 encode_ieee_single (const struct real_format *fmt, long *buf,
2641 const REAL_VALUE_TYPE *r)
2643 unsigned long image, sig, exp;
2644 unsigned long sign = r->sign;
2645 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2648 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
2659 image |= 0x7fffffff;
2667 if (r->signalling == fmt->qnan_msb_set)
2671 /* We overload qnan_msb_set here: it's only clear for
2672 mips_ieee_single, which wants all mantissa bits but the
2673 quiet/signalling one set in canonical NaNs (at least
2675 if (r->canonical && !fmt->qnan_msb_set)
2676 sig |= (1 << 22) - 1;
2684 image |= 0x7fffffff;
2688 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2689 whereas the intermediate representation is 0.F x 2**exp.
2690 Which means we're off by one. */
2694 exp = REAL_EXP (r) + 127 - 1;
2707 decode_ieee_single (const struct real_format *fmt, REAL_VALUE_TYPE *r,
2710 unsigned long image = buf[0] & 0xffffffff;
2711 bool sign = (image >> 31) & 1;
2712 int exp = (image >> 23) & 0xff;
2714 memset (r, 0, sizeof (*r));
2715 image <<= HOST_BITS_PER_LONG - 24;
2720 if (image && fmt->has_denorm)
2724 SET_REAL_EXP (r, -126);
2725 r->sig[SIGSZ-1] = image << 1;
2728 else if (fmt->has_signed_zero)
2731 else if (exp == 255 && (fmt->has_nans || fmt->has_inf))
2737 r->signalling = (((image >> (HOST_BITS_PER_LONG - 2)) & 1)
2738 ^ fmt->qnan_msb_set);
2739 r->sig[SIGSZ-1] = image;
2751 SET_REAL_EXP (r, exp - 127 + 1);
2752 r->sig[SIGSZ-1] = image | SIG_MSB;
2756 const struct real_format ieee_single_format =
2775 const struct real_format mips_single_format =
2795 /* IEEE double-precision format. */
2797 static void encode_ieee_double (const struct real_format *fmt,
2798 long *, const REAL_VALUE_TYPE *);
2799 static void decode_ieee_double (const struct real_format *,
2800 REAL_VALUE_TYPE *, const long *);
2803 encode_ieee_double (const struct real_format *fmt, long *buf,
2804 const REAL_VALUE_TYPE *r)
2806 unsigned long image_lo, image_hi, sig_lo, sig_hi, exp;
2807 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2809 image_hi = r->sign << 31;
2812 if (HOST_BITS_PER_LONG == 64)
2814 sig_hi = r->sig[SIGSZ-1];
2815 sig_lo = (sig_hi >> (64 - 53)) & 0xffffffff;
2816 sig_hi = (sig_hi >> (64 - 53 + 1) >> 31) & 0xfffff;
2820 sig_hi = r->sig[SIGSZ-1];
2821 sig_lo = r->sig[SIGSZ-2];
2822 sig_lo = (sig_hi << 21) | (sig_lo >> 11);
2823 sig_hi = (sig_hi >> 11) & 0xfffff;
2833 image_hi |= 2047 << 20;
2836 image_hi |= 0x7fffffff;
2837 image_lo = 0xffffffff;
2845 sig_hi = sig_lo = 0;
2846 if (r->signalling == fmt->qnan_msb_set)
2847 sig_hi &= ~(1 << 19);
2850 /* We overload qnan_msb_set here: it's only clear for
2851 mips_ieee_single, which wants all mantissa bits but the
2852 quiet/signalling one set in canonical NaNs (at least
2854 if (r->canonical && !fmt->qnan_msb_set)
2856 sig_hi |= (1 << 19) - 1;
2857 sig_lo = 0xffffffff;
2859 else if (sig_hi == 0 && sig_lo == 0)
2862 image_hi |= 2047 << 20;
2868 image_hi |= 0x7fffffff;
2869 image_lo = 0xffffffff;
2874 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2875 whereas the intermediate representation is 0.F x 2**exp.
2876 Which means we're off by one. */
2880 exp = REAL_EXP (r) + 1023 - 1;
2881 image_hi |= exp << 20;
2890 if (FLOAT_WORDS_BIG_ENDIAN)
2891 buf[0] = image_hi, buf[1] = image_lo;
2893 buf[0] = image_lo, buf[1] = image_hi;
2897 decode_ieee_double (const struct real_format *fmt, REAL_VALUE_TYPE *r,
2900 unsigned long image_hi, image_lo;
2904 if (FLOAT_WORDS_BIG_ENDIAN)
2905 image_hi = buf[0], image_lo = buf[1];
2907 image_lo = buf[0], image_hi = buf[1];
2908 image_lo &= 0xffffffff;
2909 image_hi &= 0xffffffff;
2911 sign = (image_hi >> 31) & 1;
2912 exp = (image_hi >> 20) & 0x7ff;
2914 memset (r, 0, sizeof (*r));
2916 image_hi <<= 32 - 21;
2917 image_hi |= image_lo >> 21;
2918 image_hi &= 0x7fffffff;
2919 image_lo <<= 32 - 21;
2923 if ((image_hi || image_lo) && fmt->has_denorm)
2927 SET_REAL_EXP (r, -1022);
2928 if (HOST_BITS_PER_LONG == 32)
2930 image_hi = (image_hi << 1) | (image_lo >> 31);
2932 r->sig[SIGSZ-1] = image_hi;
2933 r->sig[SIGSZ-2] = image_lo;
2937 image_hi = (image_hi << 31 << 2) | (image_lo << 1);
2938 r->sig[SIGSZ-1] = image_hi;
2942 else if (fmt->has_signed_zero)
2945 else if (exp == 2047 && (fmt->has_nans || fmt->has_inf))
2947 if (image_hi || image_lo)
2951 r->signalling = ((image_hi >> 30) & 1) ^ fmt->qnan_msb_set;
2952 if (HOST_BITS_PER_LONG == 32)
2954 r->sig[SIGSZ-1] = image_hi;
2955 r->sig[SIGSZ-2] = image_lo;
2958 r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo;
2970 SET_REAL_EXP (r, exp - 1023 + 1);
2971 if (HOST_BITS_PER_LONG == 32)
2973 r->sig[SIGSZ-1] = image_hi | SIG_MSB;
2974 r->sig[SIGSZ-2] = image_lo;
2977 r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo | SIG_MSB;
2981 const struct real_format ieee_double_format =
3000 const struct real_format mips_double_format =
3020 /* IEEE extended real format. This comes in three flavors: Intel's as
3021 a 12 byte image, Intel's as a 16 byte image, and Motorola's. Intel
3022 12- and 16-byte images may be big- or little endian; Motorola's is
3023 always big endian. */
3025 /* Helper subroutine which converts from the internal format to the
3026 12-byte little-endian Intel format. Functions below adjust this
3027 for the other possible formats. */
3029 encode_ieee_extended (const struct real_format *fmt, long *buf,
3030 const REAL_VALUE_TYPE *r)
3032 unsigned long image_hi, sig_hi, sig_lo;
3033 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3035 image_hi = r->sign << 15;
3036 sig_hi = sig_lo = 0;
3048 /* Intel requires the explicit integer bit to be set, otherwise
3049 it considers the value a "pseudo-infinity". Motorola docs
3050 say it doesn't care. */
3051 sig_hi = 0x80000000;
3056 sig_lo = sig_hi = 0xffffffff;
3064 if (HOST_BITS_PER_LONG == 32)
3066 sig_hi = r->sig[SIGSZ-1];
3067 sig_lo = r->sig[SIGSZ-2];
3071 sig_lo = r->sig[SIGSZ-1];
3072 sig_hi = sig_lo >> 31 >> 1;
3073 sig_lo &= 0xffffffff;
3075 if (r->signalling == fmt->qnan_msb_set)
3076 sig_hi &= ~(1 << 30);
3079 if ((sig_hi & 0x7fffffff) == 0 && sig_lo == 0)
3082 /* Intel requires the explicit integer bit to be set, otherwise
3083 it considers the value a "pseudo-nan". Motorola docs say it
3085 sig_hi |= 0x80000000;
3090 sig_lo = sig_hi = 0xffffffff;
3096 int exp = REAL_EXP (r);
3098 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3099 whereas the intermediate representation is 0.F x 2**exp.
3100 Which means we're off by one.
3102 Except for Motorola, which consider exp=0 and explicit
3103 integer bit set to continue to be normalized. In theory
3104 this discrepancy has been taken care of by the difference
3105 in fmt->emin in round_for_format. */
3112 gcc_assert (exp >= 0);
3116 if (HOST_BITS_PER_LONG == 32)
3118 sig_hi = r->sig[SIGSZ-1];
3119 sig_lo = r->sig[SIGSZ-2];
3123 sig_lo = r->sig[SIGSZ-1];
3124 sig_hi = sig_lo >> 31 >> 1;
3125 sig_lo &= 0xffffffff;
3134 buf[0] = sig_lo, buf[1] = sig_hi, buf[2] = image_hi;
3137 /* Convert from the internal format to the 12-byte Motorola format
3138 for an IEEE extended real. */
3140 encode_ieee_extended_motorola (const struct real_format *fmt, long *buf,
3141 const REAL_VALUE_TYPE *r)
3144 encode_ieee_extended (fmt, intermed, r);
3146 /* Motorola chips are assumed always to be big-endian. Also, the
3147 padding in a Motorola extended real goes between the exponent and
3148 the mantissa. At this point the mantissa is entirely within
3149 elements 0 and 1 of intermed, and the exponent entirely within
3150 element 2, so all we have to do is swap the order around, and
3151 shift element 2 left 16 bits. */
3152 buf[0] = intermed[2] << 16;
3153 buf[1] = intermed[1];
3154 buf[2] = intermed[0];
3157 /* Convert from the internal format to the 12-byte Intel format for
3158 an IEEE extended real. */
3160 encode_ieee_extended_intel_96 (const struct real_format *fmt, long *buf,
3161 const REAL_VALUE_TYPE *r)
3163 if (FLOAT_WORDS_BIG_ENDIAN)
3165 /* All the padding in an Intel-format extended real goes at the high
3166 end, which in this case is after the mantissa, not the exponent.
3167 Therefore we must shift everything down 16 bits. */
3169 encode_ieee_extended (fmt, intermed, r);
3170 buf[0] = ((intermed[2] << 16) | ((unsigned long)(intermed[1] & 0xFFFF0000) >> 16));
3171 buf[1] = ((intermed[1] << 16) | ((unsigned long)(intermed[0] & 0xFFFF0000) >> 16));
3172 buf[2] = (intermed[0] << 16);
3175 /* encode_ieee_extended produces what we want directly. */
3176 encode_ieee_extended (fmt, buf, r);
3179 /* Convert from the internal format to the 16-byte Intel format for
3180 an IEEE extended real. */
3182 encode_ieee_extended_intel_128 (const struct real_format *fmt, long *buf,
3183 const REAL_VALUE_TYPE *r)
3185 /* All the padding in an Intel-format extended real goes at the high end. */
3186 encode_ieee_extended_intel_96 (fmt, buf, r);
3190 /* As above, we have a helper function which converts from 12-byte
3191 little-endian Intel format to internal format. Functions below
3192 adjust for the other possible formats. */
3194 decode_ieee_extended (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3197 unsigned long image_hi, sig_hi, sig_lo;
3201 sig_lo = buf[0], sig_hi = buf[1], image_hi = buf[2];
3202 sig_lo &= 0xffffffff;
3203 sig_hi &= 0xffffffff;
3204 image_hi &= 0xffffffff;
3206 sign = (image_hi >> 15) & 1;
3207 exp = image_hi & 0x7fff;
3209 memset (r, 0, sizeof (*r));
3213 if ((sig_hi || sig_lo) && fmt->has_denorm)
3218 /* When the IEEE format contains a hidden bit, we know that
3219 it's zero at this point, and so shift up the significand
3220 and decrease the exponent to match. In this case, Motorola
3221 defines the explicit integer bit to be valid, so we don't
3222 know whether the msb is set or not. */
3223 SET_REAL_EXP (r, fmt->emin);
3224 if (HOST_BITS_PER_LONG == 32)
3226 r->sig[SIGSZ-1] = sig_hi;
3227 r->sig[SIGSZ-2] = sig_lo;
3230 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3234 else if (fmt->has_signed_zero)
3237 else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3239 /* See above re "pseudo-infinities" and "pseudo-nans".
3240 Short summary is that the MSB will likely always be
3241 set, and that we don't care about it. */
3242 sig_hi &= 0x7fffffff;
3244 if (sig_hi || sig_lo)
3248 r->signalling = ((sig_hi >> 30) & 1) ^ fmt->qnan_msb_set;
3249 if (HOST_BITS_PER_LONG == 32)
3251 r->sig[SIGSZ-1] = sig_hi;
3252 r->sig[SIGSZ-2] = sig_lo;
3255 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3267 SET_REAL_EXP (r, exp - 16383 + 1);
3268 if (HOST_BITS_PER_LONG == 32)
3270 r->sig[SIGSZ-1] = sig_hi;
3271 r->sig[SIGSZ-2] = sig_lo;
3274 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3278 /* Convert from the internal format to the 12-byte Motorola format
3279 for an IEEE extended real. */
3281 decode_ieee_extended_motorola (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3286 /* Motorola chips are assumed always to be big-endian. Also, the
3287 padding in a Motorola extended real goes between the exponent and
3288 the mantissa; remove it. */
3289 intermed[0] = buf[2];
3290 intermed[1] = buf[1];
3291 intermed[2] = (unsigned long)buf[0] >> 16;
3293 decode_ieee_extended (fmt, r, intermed);
3296 /* Convert from the internal format to the 12-byte Intel format for
3297 an IEEE extended real. */
3299 decode_ieee_extended_intel_96 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3302 if (FLOAT_WORDS_BIG_ENDIAN)
3304 /* All the padding in an Intel-format extended real goes at the high
3305 end, which in this case is after the mantissa, not the exponent.
3306 Therefore we must shift everything up 16 bits. */
3309 intermed[0] = (((unsigned long)buf[2] >> 16) | (buf[1] << 16));
3310 intermed[1] = (((unsigned long)buf[1] >> 16) | (buf[0] << 16));
3311 intermed[2] = ((unsigned long)buf[0] >> 16);
3313 decode_ieee_extended (fmt, r, intermed);
3316 /* decode_ieee_extended produces what we want directly. */
3317 decode_ieee_extended (fmt, r, buf);
3320 /* Convert from the internal format to the 16-byte Intel format for
3321 an IEEE extended real. */
3323 decode_ieee_extended_intel_128 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3326 /* All the padding in an Intel-format extended real goes at the high end. */
3327 decode_ieee_extended_intel_96 (fmt, r, buf);
3330 const struct real_format ieee_extended_motorola_format =
3332 encode_ieee_extended_motorola,
3333 decode_ieee_extended_motorola,
3349 const struct real_format ieee_extended_intel_96_format =
3351 encode_ieee_extended_intel_96,
3352 decode_ieee_extended_intel_96,
3368 const struct real_format ieee_extended_intel_128_format =
3370 encode_ieee_extended_intel_128,
3371 decode_ieee_extended_intel_128,
3387 /* The following caters to i386 systems that set the rounding precision
3388 to 53 bits instead of 64, e.g. FreeBSD. */
3389 const struct real_format ieee_extended_intel_96_round_53_format =
3391 encode_ieee_extended_intel_96,
3392 decode_ieee_extended_intel_96,
3408 /* IBM 128-bit extended precision format: a pair of IEEE double precision
3409 numbers whose sum is equal to the extended precision value. The number
3410 with greater magnitude is first. This format has the same magnitude
3411 range as an IEEE double precision value, but effectively 106 bits of
3412 significand precision. Infinity and NaN are represented by their IEEE
3413 double precision value stored in the first number, the second number is
3414 +0.0 or -0.0 for Infinity and don't-care for NaN. */
3416 static void encode_ibm_extended (const struct real_format *fmt,
3417 long *, const REAL_VALUE_TYPE *);
3418 static void decode_ibm_extended (const struct real_format *,
3419 REAL_VALUE_TYPE *, const long *);
3422 encode_ibm_extended (const struct real_format *fmt, long *buf,
3423 const REAL_VALUE_TYPE *r)
3425 REAL_VALUE_TYPE u, normr, v;
3426 const struct real_format *base_fmt;
3428 base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3430 /* Renormlize R before doing any arithmetic on it. */
3432 if (normr.cl == rvc_normal)
3435 /* u = IEEE double precision portion of significand. */
3437 round_for_format (base_fmt, &u);
3438 encode_ieee_double (base_fmt, &buf[0], &u);
3440 if (u.cl == rvc_normal)
3442 do_add (&v, &normr, &u, 1);
3443 /* Call round_for_format since we might need to denormalize. */
3444 round_for_format (base_fmt, &v);
3445 encode_ieee_double (base_fmt, &buf[2], &v);
3449 /* Inf, NaN, 0 are all representable as doubles, so the
3450 least-significant part can be 0.0. */
3457 decode_ibm_extended (const struct real_format *fmt ATTRIBUTE_UNUSED, REAL_VALUE_TYPE *r,
3460 REAL_VALUE_TYPE u, v;
3461 const struct real_format *base_fmt;
3463 base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3464 decode_ieee_double (base_fmt, &u, &buf[0]);
3466 if (u.cl != rvc_zero && u.cl != rvc_inf && u.cl != rvc_nan)
3468 decode_ieee_double (base_fmt, &v, &buf[2]);
3469 do_add (r, &u, &v, 0);
3475 const struct real_format ibm_extended_format =
3477 encode_ibm_extended,
3478 decode_ibm_extended,
3494 const struct real_format mips_extended_format =
3496 encode_ibm_extended,
3497 decode_ibm_extended,
3514 /* IEEE quad precision format. */
3516 static void encode_ieee_quad (const struct real_format *fmt,
3517 long *, const REAL_VALUE_TYPE *);
3518 static void decode_ieee_quad (const struct real_format *,
3519 REAL_VALUE_TYPE *, const long *);
3522 encode_ieee_quad (const struct real_format *fmt, long *buf,
3523 const REAL_VALUE_TYPE *r)
3525 unsigned long image3, image2, image1, image0, exp;
3526 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3529 image3 = r->sign << 31;
3534 rshift_significand (&u, r, SIGNIFICAND_BITS - 113);
3543 image3 |= 32767 << 16;
3546 image3 |= 0x7fffffff;
3547 image2 = 0xffffffff;
3548 image1 = 0xffffffff;
3549 image0 = 0xffffffff;
3556 image3 |= 32767 << 16;
3560 /* Don't use bits from the significand. The
3561 initialization above is right. */
3563 else if (HOST_BITS_PER_LONG == 32)
3568 image3 |= u.sig[3] & 0xffff;
3573 image1 = image0 >> 31 >> 1;
3575 image3 |= (image2 >> 31 >> 1) & 0xffff;
3576 image0 &= 0xffffffff;
3577 image2 &= 0xffffffff;
3579 if (r->signalling == fmt->qnan_msb_set)
3583 /* We overload qnan_msb_set here: it's only clear for
3584 mips_ieee_single, which wants all mantissa bits but the
3585 quiet/signalling one set in canonical NaNs (at least
3587 if (r->canonical && !fmt->qnan_msb_set)
3590 image2 = image1 = image0 = 0xffffffff;
3592 else if (((image3 & 0xffff) | image2 | image1 | image0) == 0)
3597 image3 |= 0x7fffffff;
3598 image2 = 0xffffffff;
3599 image1 = 0xffffffff;
3600 image0 = 0xffffffff;
3605 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3606 whereas the intermediate representation is 0.F x 2**exp.
3607 Which means we're off by one. */
3611 exp = REAL_EXP (r) + 16383 - 1;
3612 image3 |= exp << 16;
3614 if (HOST_BITS_PER_LONG == 32)
3619 image3 |= u.sig[3] & 0xffff;
3624 image1 = image0 >> 31 >> 1;
3626 image3 |= (image2 >> 31 >> 1) & 0xffff;
3627 image0 &= 0xffffffff;
3628 image2 &= 0xffffffff;
3636 if (FLOAT_WORDS_BIG_ENDIAN)
3653 decode_ieee_quad (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3656 unsigned long image3, image2, image1, image0;
3660 if (FLOAT_WORDS_BIG_ENDIAN)
3674 image0 &= 0xffffffff;
3675 image1 &= 0xffffffff;
3676 image2 &= 0xffffffff;
3678 sign = (image3 >> 31) & 1;
3679 exp = (image3 >> 16) & 0x7fff;
3682 memset (r, 0, sizeof (*r));
3686 if ((image3 | image2 | image1 | image0) && fmt->has_denorm)
3691 SET_REAL_EXP (r, -16382 + (SIGNIFICAND_BITS - 112));
3692 if (HOST_BITS_PER_LONG == 32)
3701 r->sig[0] = (image1 << 31 << 1) | image0;
3702 r->sig[1] = (image3 << 31 << 1) | image2;
3707 else if (fmt->has_signed_zero)
3710 else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3712 if (image3 | image2 | image1 | image0)
3716 r->signalling = ((image3 >> 15) & 1) ^ fmt->qnan_msb_set;
3718 if (HOST_BITS_PER_LONG == 32)
3727 r->sig[0] = (image1 << 31 << 1) | image0;
3728 r->sig[1] = (image3 << 31 << 1) | image2;
3730 lshift_significand (r, r, SIGNIFICAND_BITS - 113);
3742 SET_REAL_EXP (r, exp - 16383 + 1);
3744 if (HOST_BITS_PER_LONG == 32)
3753 r->sig[0] = (image1 << 31 << 1) | image0;
3754 r->sig[1] = (image3 << 31 << 1) | image2;
3756 lshift_significand (r, r, SIGNIFICAND_BITS - 113);
3757 r->sig[SIGSZ-1] |= SIG_MSB;
3761 const struct real_format ieee_quad_format =
3780 const struct real_format mips_quad_format =
3799 /* Descriptions of VAX floating point formats can be found beginning at
3801 http://h71000.www7.hp.com/doc/73FINAL/4515/4515pro_013.html#f_floating_point_format
3803 The thing to remember is that they're almost IEEE, except for word
3804 order, exponent bias, and the lack of infinities, nans, and denormals.
3806 We don't implement the H_floating format here, simply because neither
3807 the VAX or Alpha ports use it. */
3809 static void encode_vax_f (const struct real_format *fmt,
3810 long *, const REAL_VALUE_TYPE *);
3811 static void decode_vax_f (const struct real_format *,
3812 REAL_VALUE_TYPE *, const long *);
3813 static void encode_vax_d (const struct real_format *fmt,
3814 long *, const REAL_VALUE_TYPE *);
3815 static void decode_vax_d (const struct real_format *,
3816 REAL_VALUE_TYPE *, const long *);
3817 static void encode_vax_g (const struct real_format *fmt,
3818 long *, const REAL_VALUE_TYPE *);
3819 static void decode_vax_g (const struct real_format *,
3820 REAL_VALUE_TYPE *, const long *);
3823 encode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3824 const REAL_VALUE_TYPE *r)
3826 unsigned long sign, exp, sig, image;
3828 sign = r->sign << 15;
3838 image = 0xffff7fff | sign;
3842 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
3843 exp = REAL_EXP (r) + 128;
3845 image = (sig << 16) & 0xffff0000;
3859 decode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED,
3860 REAL_VALUE_TYPE *r, const long *buf)
3862 unsigned long image = buf[0] & 0xffffffff;
3863 int exp = (image >> 7) & 0xff;
3865 memset (r, 0, sizeof (*r));
3870 r->sign = (image >> 15) & 1;
3871 SET_REAL_EXP (r, exp - 128);
3873 image = ((image & 0x7f) << 16) | ((image >> 16) & 0xffff);
3874 r->sig[SIGSZ-1] = (image << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
3879 encode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3880 const REAL_VALUE_TYPE *r)
3882 unsigned long image0, image1, sign = r->sign << 15;
3887 image0 = image1 = 0;
3892 image0 = 0xffff7fff | sign;
3893 image1 = 0xffffffff;
3897 /* Extract the significand into straight hi:lo. */
3898 if (HOST_BITS_PER_LONG == 64)
3900 image0 = r->sig[SIGSZ-1];
3901 image1 = (image0 >> (64 - 56)) & 0xffffffff;
3902 image0 = (image0 >> (64 - 56 + 1) >> 31) & 0x7fffff;
3906 image0 = r->sig[SIGSZ-1];
3907 image1 = r->sig[SIGSZ-2];
3908 image1 = (image0 << 24) | (image1 >> 8);
3909 image0 = (image0 >> 8) & 0xffffff;
3912 /* Rearrange the half-words of the significand to match the
3914 image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff007f;
3915 image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
3917 /* Add the sign and exponent. */
3919 image0 |= (REAL_EXP (r) + 128) << 7;
3926 if (FLOAT_WORDS_BIG_ENDIAN)
3927 buf[0] = image1, buf[1] = image0;
3929 buf[0] = image0, buf[1] = image1;
3933 decode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED,
3934 REAL_VALUE_TYPE *r, const long *buf)
3936 unsigned long image0, image1;
3939 if (FLOAT_WORDS_BIG_ENDIAN)
3940 image1 = buf[0], image0 = buf[1];
3942 image0 = buf[0], image1 = buf[1];
3943 image0 &= 0xffffffff;
3944 image1 &= 0xffffffff;
3946 exp = (image0 >> 7) & 0xff;
3948 memset (r, 0, sizeof (*r));
3953 r->sign = (image0 >> 15) & 1;
3954 SET_REAL_EXP (r, exp - 128);
3956 /* Rearrange the half-words of the external format into
3957 proper ascending order. */
3958 image0 = ((image0 & 0x7f) << 16) | ((image0 >> 16) & 0xffff);
3959 image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
3961 if (HOST_BITS_PER_LONG == 64)
3963 image0 = (image0 << 31 << 1) | image1;
3966 r->sig[SIGSZ-1] = image0;
3970 r->sig[SIGSZ-1] = image0;
3971 r->sig[SIGSZ-2] = image1;
3972 lshift_significand (r, r, 2*HOST_BITS_PER_LONG - 56);
3973 r->sig[SIGSZ-1] |= SIG_MSB;
3979 encode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3980 const REAL_VALUE_TYPE *r)
3982 unsigned long image0, image1, sign = r->sign << 15;
3987 image0 = image1 = 0;
3992 image0 = 0xffff7fff | sign;
3993 image1 = 0xffffffff;
3997 /* Extract the significand into straight hi:lo. */
3998 if (HOST_BITS_PER_LONG == 64)
4000 image0 = r->sig[SIGSZ-1];
4001 image1 = (image0 >> (64 - 53)) & 0xffffffff;
4002 image0 = (image0 >> (64 - 53 + 1) >> 31) & 0xfffff;
4006 image0 = r->sig[SIGSZ-1];
4007 image1 = r->sig[SIGSZ-2];
4008 image1 = (image0 << 21) | (image1 >> 11);
4009 image0 = (image0 >> 11) & 0xfffff;
4012 /* Rearrange the half-words of the significand to match the
4014 image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff000f;
4015 image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
4017 /* Add the sign and exponent. */
4019 image0 |= (REAL_EXP (r) + 1024) << 4;
4026 if (FLOAT_WORDS_BIG_ENDIAN)
4027 buf[0] = image1, buf[1] = image0;
4029 buf[0] = image0, buf[1] = image1;
4033 decode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED,
4034 REAL_VALUE_TYPE *r, const long *buf)
4036 unsigned long image0, image1;
4039 if (FLOAT_WORDS_BIG_ENDIAN)
4040 image1 = buf[0], image0 = buf[1];
4042 image0 = buf[0], image1 = buf[1];
4043 image0 &= 0xffffffff;
4044 image1 &= 0xffffffff;
4046 exp = (image0 >> 4) & 0x7ff;
4048 memset (r, 0, sizeof (*r));
4053 r->sign = (image0 >> 15) & 1;
4054 SET_REAL_EXP (r, exp - 1024);
4056 /* Rearrange the half-words of the external format into
4057 proper ascending order. */
4058 image0 = ((image0 & 0xf) << 16) | ((image0 >> 16) & 0xffff);
4059 image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
4061 if (HOST_BITS_PER_LONG == 64)
4063 image0 = (image0 << 31 << 1) | image1;
4066 r->sig[SIGSZ-1] = image0;
4070 r->sig[SIGSZ-1] = image0;
4071 r->sig[SIGSZ-2] = image1;
4072 lshift_significand (r, r, 64 - 53);
4073 r->sig[SIGSZ-1] |= SIG_MSB;
4078 const struct real_format vax_f_format =
4097 const struct real_format vax_d_format =
4116 const struct real_format vax_g_format =
4135 /* A good reference for these can be found in chapter 9 of
4136 "ESA/390 Principles of Operation", IBM document number SA22-7201-01.
4137 An on-line version can be found here:
4139 http://publibz.boulder.ibm.com/cgi-bin/bookmgr_OS390/BOOKS/DZ9AR001/9.1?DT=19930923083613
4142 static void encode_i370_single (const struct real_format *fmt,
4143 long *, const REAL_VALUE_TYPE *);
4144 static void decode_i370_single (const struct real_format *,
4145 REAL_VALUE_TYPE *, const long *);
4146 static void encode_i370_double (const struct real_format *fmt,
4147 long *, const REAL_VALUE_TYPE *);
4148 static void decode_i370_double (const struct real_format *,
4149 REAL_VALUE_TYPE *, const long *);
4152 encode_i370_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4153 long *buf, const REAL_VALUE_TYPE *r)
4155 unsigned long sign, exp, sig, image;
4157 sign = r->sign << 31;
4167 image = 0x7fffffff | sign;
4171 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0xffffff;
4172 exp = ((REAL_EXP (r) / 4) + 64) << 24;
4173 image = sign | exp | sig;
4184 decode_i370_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4185 REAL_VALUE_TYPE *r, const long *buf)
4187 unsigned long sign, sig, image = buf[0];
4190 sign = (image >> 31) & 1;
4191 exp = (image >> 24) & 0x7f;
4192 sig = image & 0xffffff;
4194 memset (r, 0, sizeof (*r));
4200 SET_REAL_EXP (r, (exp - 64) * 4);
4201 r->sig[SIGSZ-1] = sig << (HOST_BITS_PER_LONG - 24);
4207 encode_i370_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4208 long *buf, const REAL_VALUE_TYPE *r)
4210 unsigned long sign, exp, image_hi, image_lo;
4212 sign = r->sign << 31;
4217 image_hi = image_lo = 0;
4222 image_hi = 0x7fffffff | sign;
4223 image_lo = 0xffffffff;
4227 if (HOST_BITS_PER_LONG == 64)
4229 image_hi = r->sig[SIGSZ-1];
4230 image_lo = (image_hi >> (64 - 56)) & 0xffffffff;
4231 image_hi = (image_hi >> (64 - 56 + 1) >> 31) & 0xffffff;
4235 image_hi = r->sig[SIGSZ-1];
4236 image_lo = r->sig[SIGSZ-2];
4237 image_lo = (image_lo >> 8) | (image_hi << 24);
4241 exp = ((REAL_EXP (r) / 4) + 64) << 24;
4242 image_hi |= sign | exp;
4249 if (FLOAT_WORDS_BIG_ENDIAN)
4250 buf[0] = image_hi, buf[1] = image_lo;
4252 buf[0] = image_lo, buf[1] = image_hi;
4256 decode_i370_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4257 REAL_VALUE_TYPE *r, const long *buf)
4259 unsigned long sign, image_hi, image_lo;
4262 if (FLOAT_WORDS_BIG_ENDIAN)
4263 image_hi = buf[0], image_lo = buf[1];
4265 image_lo = buf[0], image_hi = buf[1];
4267 sign = (image_hi >> 31) & 1;
4268 exp = (image_hi >> 24) & 0x7f;
4269 image_hi &= 0xffffff;
4270 image_lo &= 0xffffffff;
4272 memset (r, 0, sizeof (*r));
4274 if (exp || image_hi || image_lo)
4278 SET_REAL_EXP (r, (exp - 64) * 4 + (SIGNIFICAND_BITS - 56));
4280 if (HOST_BITS_PER_LONG == 32)
4282 r->sig[0] = image_lo;
4283 r->sig[1] = image_hi;
4286 r->sig[0] = image_lo | (image_hi << 31 << 1);
4292 const struct real_format i370_single_format =
4306 false, /* ??? The encoding does allow for "unnormals". */
4307 false, /* ??? The encoding does allow for "unnormals". */
4311 const struct real_format i370_double_format =
4325 false, /* ??? The encoding does allow for "unnormals". */
4326 false, /* ??? The encoding does allow for "unnormals". */
4330 /* Encode real R into a single precision DFP value in BUF. */
4332 encode_decimal_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4333 long *buf ATTRIBUTE_UNUSED,
4334 const REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED)
4336 encode_decimal32 (fmt, buf, r);
4339 /* Decode a single precision DFP value in BUF into a real R. */
4341 decode_decimal_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4342 REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED,
4343 const long *buf ATTRIBUTE_UNUSED)
4345 decode_decimal32 (fmt, r, buf);
4348 /* Encode real R into a double precision DFP value in BUF. */
4350 encode_decimal_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4351 long *buf ATTRIBUTE_UNUSED,
4352 const REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED)
4354 encode_decimal64 (fmt, buf, r);
4357 /* Decode a double precision DFP value in BUF into a real R. */
4359 decode_decimal_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4360 REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED,
4361 const long *buf ATTRIBUTE_UNUSED)
4363 decode_decimal64 (fmt, r, buf);
4366 /* Encode real R into a quad precision DFP value in BUF. */
4368 encode_decimal_quad (const struct real_format *fmt ATTRIBUTE_UNUSED,
4369 long *buf ATTRIBUTE_UNUSED,
4370 const REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED)
4372 encode_decimal128 (fmt, buf, r);
4375 /* Decode a quad precision DFP value in BUF into a real R. */
4377 decode_decimal_quad (const struct real_format *fmt ATTRIBUTE_UNUSED,
4378 REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED,
4379 const long *buf ATTRIBUTE_UNUSED)
4381 decode_decimal128 (fmt, r, buf);
4384 /* Single precision decimal floating point (IEEE 754R). */
4385 const struct real_format decimal_single_format =
4387 encode_decimal_single,
4388 decode_decimal_single,
4404 /* Double precision decimal floating point (IEEE 754R). */
4405 const struct real_format decimal_double_format =
4407 encode_decimal_double,
4408 decode_decimal_double,
4424 /* Quad precision decimal floating point (IEEE 754R). */
4425 const struct real_format decimal_quad_format =
4427 encode_decimal_quad,
4428 decode_decimal_quad,
4444 /* The "twos-complement" c4x format is officially defined as
4448 This is rather misleading. One must remember that F is signed.
4449 A better description would be
4451 x = -1**s * ((s + 1 + .f) * 2**e
4453 So if we have a (4 bit) fraction of .1000 with a sign bit of 1,
4454 that's -1 * (1+1+(-.5)) == -1.5. I think.
4456 The constructions here are taken from Tables 5-1 and 5-2 of the
4457 TMS320C4x User's Guide wherein step-by-step instructions for
4458 conversion from IEEE are presented. That's close enough to our
4459 internal representation so as to make things easy.
4461 See http://www-s.ti.com/sc/psheets/spru063c/spru063c.pdf */
4463 static void encode_c4x_single (const struct real_format *fmt,
4464 long *, const REAL_VALUE_TYPE *);
4465 static void decode_c4x_single (const struct real_format *,
4466 REAL_VALUE_TYPE *, const long *);
4467 static void encode_c4x_extended (const struct real_format *fmt,
4468 long *, const REAL_VALUE_TYPE *);
4469 static void decode_c4x_extended (const struct real_format *,
4470 REAL_VALUE_TYPE *, const long *);
4473 encode_c4x_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4474 long *buf, const REAL_VALUE_TYPE *r)
4476 unsigned long image, exp, sig;
4488 sig = 0x800000 - r->sign;
4492 exp = REAL_EXP (r) - 1;
4493 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
4508 image = ((exp & 0xff) << 24) | (sig & 0xffffff);
4513 decode_c4x_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4514 REAL_VALUE_TYPE *r, const long *buf)
4516 unsigned long image = buf[0];
4520 exp = (((image >> 24) & 0xff) ^ 0x80) - 0x80;
4521 sf = ((image & 0xffffff) ^ 0x800000) - 0x800000;
4523 memset (r, 0, sizeof (*r));
4529 sig = sf & 0x7fffff;
4538 sig = (sig << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
4540 SET_REAL_EXP (r, exp + 1);
4541 r->sig[SIGSZ-1] = sig;
4546 encode_c4x_extended (const struct real_format *fmt ATTRIBUTE_UNUSED,
4547 long *buf, const REAL_VALUE_TYPE *r)
4549 unsigned long exp, sig;
4561 sig = 0x80000000 - r->sign;
4565 exp = REAL_EXP (r) - 1;
4567 sig = r->sig[SIGSZ-1];
4568 if (HOST_BITS_PER_LONG == 64)
4569 sig = sig >> 1 >> 31;
4586 exp = (exp & 0xff) << 24;
4589 if (FLOAT_WORDS_BIG_ENDIAN)
4590 buf[0] = exp, buf[1] = sig;
4592 buf[0] = sig, buf[0] = exp;
4596 decode_c4x_extended (const struct real_format *fmt ATTRIBUTE_UNUSED,
4597 REAL_VALUE_TYPE *r, const long *buf)
4602 if (FLOAT_WORDS_BIG_ENDIAN)
4603 exp = buf[0], sf = buf[1];
4605 sf = buf[0], exp = buf[1];
4607 exp = (((exp >> 24) & 0xff) & 0x80) - 0x80;
4608 sf = ((sf & 0xffffffff) ^ 0x80000000) - 0x80000000;
4610 memset (r, 0, sizeof (*r));
4616 sig = sf & 0x7fffffff;
4625 if (HOST_BITS_PER_LONG == 64)
4626 sig = sig << 1 << 31;
4629 SET_REAL_EXP (r, exp + 1);
4630 r->sig[SIGSZ-1] = sig;
4634 const struct real_format c4x_single_format =
4653 const struct real_format c4x_extended_format =
4655 encode_c4x_extended,
4656 decode_c4x_extended,
4673 /* A synthetic "format" for internal arithmetic. It's the size of the
4674 internal significand minus the two bits needed for proper rounding.
4675 The encode and decode routines exist only to satisfy our paranoia
4678 static void encode_internal (const struct real_format *fmt,
4679 long *, const REAL_VALUE_TYPE *);
4680 static void decode_internal (const struct real_format *,
4681 REAL_VALUE_TYPE *, const long *);
4684 encode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4685 const REAL_VALUE_TYPE *r)
4687 memcpy (buf, r, sizeof (*r));
4691 decode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED,
4692 REAL_VALUE_TYPE *r, const long *buf)
4694 memcpy (r, buf, sizeof (*r));
4697 const struct real_format real_internal_format =
4703 SIGNIFICAND_BITS - 2,
4704 SIGNIFICAND_BITS - 2,
4716 /* Calculate the square root of X in mode MODE, and store the result
4717 in R. Return TRUE if the operation does not raise an exception.
4718 For details see "High Precision Division and Square Root",
4719 Alan H. Karp and Peter Markstein, HP Lab Report 93-93-42, June
4720 1993. http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf. */
4723 real_sqrt (REAL_VALUE_TYPE *r, enum machine_mode mode,
4724 const REAL_VALUE_TYPE *x)
4726 static REAL_VALUE_TYPE halfthree;
4727 static bool init = false;
4728 REAL_VALUE_TYPE h, t, i;
4731 /* sqrt(-0.0) is -0.0. */
4732 if (real_isnegzero (x))
4738 /* Negative arguments return NaN. */
4741 get_canonical_qnan (r, 0);
4745 /* Infinity and NaN return themselves. */
4746 if (real_isinf (x) || real_isnan (x))
4754 do_add (&halfthree, &dconst1, &dconsthalf, 0);
4758 /* Initial guess for reciprocal sqrt, i. */
4759 exp = real_exponent (x);
4760 real_ldexp (&i, &dconst1, -exp/2);
4762 /* Newton's iteration for reciprocal sqrt, i. */
4763 for (iter = 0; iter < 16; iter++)
4765 /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x). */
4766 do_multiply (&t, x, &i);
4767 do_multiply (&h, &t, &i);
4768 do_multiply (&t, &h, &dconsthalf);
4769 do_add (&h, &halfthree, &t, 1);
4770 do_multiply (&t, &i, &h);
4772 /* Check for early convergence. */
4773 if (iter >= 6 && real_identical (&i, &t))
4776 /* ??? Unroll loop to avoid copying. */
4780 /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)). */
4781 do_multiply (&t, x, &i);
4782 do_multiply (&h, &t, &i);
4783 do_add (&i, &dconst1, &h, 1);
4784 do_multiply (&h, &t, &i);
4785 do_multiply (&i, &dconsthalf, &h);
4786 do_add (&h, &t, &i, 0);
4788 /* ??? We need a Tuckerman test to get the last bit. */
4790 real_convert (r, mode, &h);
4794 /* Calculate X raised to the integer exponent N in mode MODE and store
4795 the result in R. Return true if the result may be inexact due to
4796 loss of precision. The algorithm is the classic "left-to-right binary
4797 method" described in section 4.6.3 of Donald Knuth's "Seminumerical
4798 Algorithms", "The Art of Computer Programming", Volume 2. */
4801 real_powi (REAL_VALUE_TYPE *r, enum machine_mode mode,
4802 const REAL_VALUE_TYPE *x, HOST_WIDE_INT n)
4804 unsigned HOST_WIDE_INT bit;
4806 bool inexact = false;
4818 /* Don't worry about overflow, from now on n is unsigned. */
4826 bit = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
4827 for (i = 0; i < HOST_BITS_PER_WIDE_INT; i++)
4831 inexact |= do_multiply (&t, &t, &t);
4833 inexact |= do_multiply (&t, &t, x);
4841 inexact |= do_divide (&t, &dconst1, &t);
4843 real_convert (r, mode, &t);
4847 /* Round X to the nearest integer not larger in absolute value, i.e.
4848 towards zero, placing the result in R in mode MODE. */
4851 real_trunc (REAL_VALUE_TYPE *r, enum machine_mode mode,
4852 const REAL_VALUE_TYPE *x)
4854 do_fix_trunc (r, x);
4855 if (mode != VOIDmode)
4856 real_convert (r, mode, r);
4859 /* Round X to the largest integer not greater in value, i.e. round
4860 down, placing the result in R in mode MODE. */
4863 real_floor (REAL_VALUE_TYPE *r, enum machine_mode mode,
4864 const REAL_VALUE_TYPE *x)
4868 do_fix_trunc (&t, x);
4869 if (! real_identical (&t, x) && x->sign)
4870 do_add (&t, &t, &dconstm1, 0);
4871 if (mode != VOIDmode)
4872 real_convert (r, mode, &t);
4877 /* Round X to the smallest integer not less then argument, i.e. round
4878 up, placing the result in R in mode MODE. */
4881 real_ceil (REAL_VALUE_TYPE *r, enum machine_mode mode,
4882 const REAL_VALUE_TYPE *x)
4886 do_fix_trunc (&t, x);
4887 if (! real_identical (&t, x) && ! x->sign)
4888 do_add (&t, &t, &dconst1, 0);
4889 if (mode != VOIDmode)
4890 real_convert (r, mode, &t);
4895 /* Round X to the nearest integer, but round halfway cases away from
4899 real_round (REAL_VALUE_TYPE *r, enum machine_mode mode,
4900 const REAL_VALUE_TYPE *x)
4902 do_add (r, x, &dconsthalf, x->sign);
4903 do_fix_trunc (r, r);
4904 if (mode != VOIDmode)
4905 real_convert (r, mode, r);
4908 /* Set the sign of R to the sign of X. */
4911 real_copysign (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *x)