* arm.md (addsf3, adddf3, subsf3, subdf3, mulsf3, muldf3, negsf2)
[official-gcc.git] / gcc / real.c
blob1e70bf7d4d95d94548a3c431b40dba3f6e83f19e
1 /* real.c - software floating point emulation.
2 Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 1999,
3 2000, 2002, 2003 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
12 version.
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
17 for more details.
19 You should have received a copy of the GNU General Public License
20 along with GCC; see the file COPYING. If not, write to the Free
21 Software Foundation, 59 Temple Place - Suite 330, Boston, MA
22 02111-1307, USA. */
24 #include "config.h"
25 #include "system.h"
26 #include "coretypes.h"
27 #include "tm.h"
28 #include "tree.h"
29 #include "toplev.h"
30 #include "real.h"
31 #include "tm_p.h"
33 /* The floating point model used internally is not exactly IEEE 754
34 compliant, and close to the description in the ISO C standard,
35 section 5.2.4.2.2 Characteristics of floating types.
37 Specifically
39 x = s * b^e * \sum_{k=1}^p f_k * b^{-k}
41 where
42 s = sign (+- 1)
43 b = base or radix, here always 2
44 e = exponent
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
50 range [0.5, 1.0).
52 A requirement of the model is that P be larger than than the
53 largest supported target floating-point type by at least 2 bits.
54 This gives us proper rounding when we truncate to the target type.
55 In addition, E must be large enough to hold the smallest supported
56 denormal number in a normalized form.
58 Both of these requirements are easily satisfied. The largest target
59 significand is 113 bits; we store at least 160. The smallest
60 denormal number fits in 17 exponent bits; we store 29.
62 Note that the decimal string conversion routines are sensitive to
63 rounding error. Since the raw arithmetic routines do not themselves
64 have guard digits or rounding, the computation of 10**exp can
65 accumulate more than a few digits of error. The previous incarnation
66 of real.c successfully used a 144 bit fraction; given the current
67 layout of REAL_VALUE_TYPE we're forced to expand to at least 160 bits.
69 Target floating point models that use base 16 instead of base 2
70 (i.e. IBM 370), are handled during round_for_format, in which we
71 canonicalize the exponent to be a multiple of 4 (log2(16)), and
72 adjust the significand to match. */
75 /* Used to classify two numbers simultaneously. */
76 #define CLASS2(A, B) ((A) << 2 | (B))
78 #if HOST_BITS_PER_LONG != 64 && HOST_BITS_PER_LONG != 32
79 #error "Some constant folding done by hand to avoid shift count warnings"
80 #endif
82 static void get_zero PARAMS ((REAL_VALUE_TYPE *, int));
83 static void get_canonical_qnan PARAMS ((REAL_VALUE_TYPE *, int));
84 static void get_canonical_snan PARAMS ((REAL_VALUE_TYPE *, int));
85 static void get_inf PARAMS ((REAL_VALUE_TYPE *, int));
86 static bool sticky_rshift_significand PARAMS ((REAL_VALUE_TYPE *,
87 const REAL_VALUE_TYPE *,
88 unsigned int));
89 static void rshift_significand PARAMS ((REAL_VALUE_TYPE *,
90 const REAL_VALUE_TYPE *,
91 unsigned int));
92 static void lshift_significand PARAMS ((REAL_VALUE_TYPE *,
93 const REAL_VALUE_TYPE *,
94 unsigned int));
95 static void lshift_significand_1 PARAMS ((REAL_VALUE_TYPE *,
96 const REAL_VALUE_TYPE *));
97 static bool add_significands PARAMS ((REAL_VALUE_TYPE *r,
98 const REAL_VALUE_TYPE *,
99 const REAL_VALUE_TYPE *));
100 static bool sub_significands PARAMS ((REAL_VALUE_TYPE *,
101 const REAL_VALUE_TYPE *,
102 const REAL_VALUE_TYPE *, int));
103 static void neg_significand PARAMS ((REAL_VALUE_TYPE *,
104 const REAL_VALUE_TYPE *));
105 static int cmp_significands PARAMS ((const REAL_VALUE_TYPE *,
106 const REAL_VALUE_TYPE *));
107 static int cmp_significand_0 PARAMS ((const REAL_VALUE_TYPE *));
108 static void set_significand_bit PARAMS ((REAL_VALUE_TYPE *, unsigned int));
109 static void clear_significand_bit PARAMS ((REAL_VALUE_TYPE *, unsigned int));
110 static bool test_significand_bit PARAMS ((REAL_VALUE_TYPE *, unsigned int));
111 static void clear_significand_below PARAMS ((REAL_VALUE_TYPE *,
112 unsigned int));
113 static bool div_significands PARAMS ((REAL_VALUE_TYPE *,
114 const REAL_VALUE_TYPE *,
115 const REAL_VALUE_TYPE *));
116 static void normalize PARAMS ((REAL_VALUE_TYPE *));
118 static void do_add PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
119 const REAL_VALUE_TYPE *, int));
120 static void do_multiply PARAMS ((REAL_VALUE_TYPE *,
121 const REAL_VALUE_TYPE *,
122 const REAL_VALUE_TYPE *));
123 static void do_divide PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
124 const REAL_VALUE_TYPE *));
125 static int do_compare PARAMS ((const REAL_VALUE_TYPE *,
126 const REAL_VALUE_TYPE *, int));
127 static void do_fix_trunc PARAMS ((REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *));
129 static unsigned long rtd_divmod PARAMS ((REAL_VALUE_TYPE *,
130 REAL_VALUE_TYPE *));
132 static const REAL_VALUE_TYPE * ten_to_ptwo PARAMS ((int));
133 static const REAL_VALUE_TYPE * ten_to_mptwo PARAMS ((int));
134 static const REAL_VALUE_TYPE * real_digit PARAMS ((int));
135 static void times_pten PARAMS ((REAL_VALUE_TYPE *, int));
137 static void round_for_format PARAMS ((const struct real_format *,
138 REAL_VALUE_TYPE *));
140 /* Initialize R with a positive zero. */
142 static inline void
143 get_zero (r, sign)
144 REAL_VALUE_TYPE *r;
145 int sign;
147 memset (r, 0, sizeof (*r));
148 r->sign = sign;
151 /* Initialize R with the canonical quiet NaN. */
153 static inline void
154 get_canonical_qnan (r, sign)
155 REAL_VALUE_TYPE *r;
156 int sign;
158 memset (r, 0, sizeof (*r));
159 r->class = rvc_nan;
160 r->sign = sign;
161 r->sig[SIGSZ-1] = SIG_MSB >> 1;
164 static inline void
165 get_canonical_snan (r, sign)
166 REAL_VALUE_TYPE *r;
167 int sign;
169 memset (r, 0, sizeof (*r));
170 r->class = rvc_nan;
171 r->sign = sign;
172 r->sig[SIGSZ-1] = SIG_MSB >> 2;
175 static inline void
176 get_inf (r, sign)
177 REAL_VALUE_TYPE *r;
178 int sign;
180 memset (r, 0, sizeof (*r));
181 r->class = rvc_inf;
182 r->sign = sign;
186 /* Right-shift the significand of A by N bits; put the result in the
187 significand of R. If any one bits are shifted out, return true. */
189 static bool
190 sticky_rshift_significand (r, a, n)
191 REAL_VALUE_TYPE *r;
192 const REAL_VALUE_TYPE *a;
193 unsigned int n;
195 unsigned long sticky = 0;
196 unsigned int i, ofs = 0;
198 if (n >= HOST_BITS_PER_LONG)
200 for (i = 0, ofs = n / HOST_BITS_PER_LONG; i < ofs; ++i)
201 sticky |= a->sig[i];
202 n &= HOST_BITS_PER_LONG - 1;
205 if (n != 0)
207 sticky |= a->sig[ofs] & (((unsigned long)1 << n) - 1);
208 for (i = 0; i < SIGSZ; ++i)
210 r->sig[i]
211 = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
212 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
213 << (HOST_BITS_PER_LONG - n)));
216 else
218 for (i = 0; ofs + i < SIGSZ; ++i)
219 r->sig[i] = a->sig[ofs + i];
220 for (; i < SIGSZ; ++i)
221 r->sig[i] = 0;
224 return sticky != 0;
227 /* Right-shift the significand of A by N bits; put the result in the
228 significand of R. */
230 static void
231 rshift_significand (r, a, n)
232 REAL_VALUE_TYPE *r;
233 const REAL_VALUE_TYPE *a;
234 unsigned int n;
236 unsigned int i, ofs = n / HOST_BITS_PER_LONG;
238 n &= HOST_BITS_PER_LONG - 1;
239 if (n != 0)
241 for (i = 0; i < SIGSZ; ++i)
243 r->sig[i]
244 = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
245 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
246 << (HOST_BITS_PER_LONG - n)));
249 else
251 for (i = 0; ofs + i < SIGSZ; ++i)
252 r->sig[i] = a->sig[ofs + i];
253 for (; i < SIGSZ; ++i)
254 r->sig[i] = 0;
258 /* Left-shift the significand of A by N bits; put the result in the
259 significand of R. */
261 static void
262 lshift_significand (r, a, n)
263 REAL_VALUE_TYPE *r;
264 const REAL_VALUE_TYPE *a;
265 unsigned int n;
267 unsigned int i, ofs = n / HOST_BITS_PER_LONG;
269 n &= HOST_BITS_PER_LONG - 1;
270 if (n == 0)
272 for (i = 0; ofs + i < SIGSZ; ++i)
273 r->sig[SIGSZ-1-i] = a->sig[SIGSZ-1-i-ofs];
274 for (; i < SIGSZ; ++i)
275 r->sig[SIGSZ-1-i] = 0;
277 else
278 for (i = 0; i < SIGSZ; ++i)
280 r->sig[SIGSZ-1-i]
281 = (((ofs + i >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs]) << n)
282 | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs-1])
283 >> (HOST_BITS_PER_LONG - n)));
287 /* Likewise, but N is specialized to 1. */
289 static inline void
290 lshift_significand_1 (r, a)
291 REAL_VALUE_TYPE *r;
292 const REAL_VALUE_TYPE *a;
294 unsigned int i;
296 for (i = SIGSZ - 1; i > 0; --i)
297 r->sig[i] = (a->sig[i] << 1) | (a->sig[i-1] >> (HOST_BITS_PER_LONG - 1));
298 r->sig[0] = a->sig[0] << 1;
301 /* Add the significands of A and B, placing the result in R. Return
302 true if there was carry out of the most significant word. */
304 static inline bool
305 add_significands (r, a, b)
306 REAL_VALUE_TYPE *r;
307 const REAL_VALUE_TYPE *a, *b;
309 bool carry = false;
310 int i;
312 for (i = 0; i < SIGSZ; ++i)
314 unsigned long ai = a->sig[i];
315 unsigned long ri = ai + b->sig[i];
317 if (carry)
319 carry = ri < ai;
320 carry |= ++ri == 0;
322 else
323 carry = ri < ai;
325 r->sig[i] = ri;
328 return carry;
331 /* Subtract the significands of A and B, placing the result in R. CARRY is
332 true if there's a borrow incoming to the least significant word.
333 Return true if there was borrow out of the most significant word. */
335 static inline bool
336 sub_significands (r, a, b, carry)
337 REAL_VALUE_TYPE *r;
338 const REAL_VALUE_TYPE *a, *b;
339 int carry;
341 int i;
343 for (i = 0; i < SIGSZ; ++i)
345 unsigned long ai = a->sig[i];
346 unsigned long ri = ai - b->sig[i];
348 if (carry)
350 carry = ri > ai;
351 carry |= ~--ri == 0;
353 else
354 carry = ri > ai;
356 r->sig[i] = ri;
359 return carry;
362 /* Negate the significand A, placing the result in R. */
364 static inline void
365 neg_significand (r, a)
366 REAL_VALUE_TYPE *r;
367 const REAL_VALUE_TYPE *a;
369 bool carry = true;
370 int i;
372 for (i = 0; i < SIGSZ; ++i)
374 unsigned long ri, ai = a->sig[i];
376 if (carry)
378 if (ai)
380 ri = -ai;
381 carry = false;
383 else
384 ri = ai;
386 else
387 ri = ~ai;
389 r->sig[i] = ri;
393 /* Compare significands. Return tri-state vs zero. */
395 static inline int
396 cmp_significands (a, b)
397 const REAL_VALUE_TYPE *a, *b;
399 int i;
401 for (i = SIGSZ - 1; i >= 0; --i)
403 unsigned long ai = a->sig[i];
404 unsigned long bi = b->sig[i];
406 if (ai > bi)
407 return 1;
408 if (ai < bi)
409 return -1;
412 return 0;
415 /* Return true if A is nonzero. */
417 static inline int
418 cmp_significand_0 (a)
419 const REAL_VALUE_TYPE *a;
421 int i;
423 for (i = SIGSZ - 1; i >= 0; --i)
424 if (a->sig[i])
425 return 1;
427 return 0;
430 /* Set bit N of the significand of R. */
432 static inline void
433 set_significand_bit (r, n)
434 REAL_VALUE_TYPE *r;
435 unsigned int n;
437 r->sig[n / HOST_BITS_PER_LONG]
438 |= (unsigned long)1 << (n % HOST_BITS_PER_LONG);
441 /* Clear bit N of the significand of R. */
443 static inline void
444 clear_significand_bit (r, n)
445 REAL_VALUE_TYPE *r;
446 unsigned int n;
448 r->sig[n / HOST_BITS_PER_LONG]
449 &= ~((unsigned long)1 << (n % HOST_BITS_PER_LONG));
452 /* Test bit N of the significand of R. */
454 static inline bool
455 test_significand_bit (r, n)
456 REAL_VALUE_TYPE *r;
457 unsigned int n;
459 /* ??? Compiler bug here if we return this expression directly.
460 The conversion to bool strips the "&1" and we wind up testing
461 e.g. 2 != 0 -> true. Seen in gcc version 3.2 20020520. */
462 int t = (r->sig[n / HOST_BITS_PER_LONG] >> (n % HOST_BITS_PER_LONG)) & 1;
463 return t;
466 /* Clear bits 0..N-1 of the significand of R. */
468 static void
469 clear_significand_below (r, n)
470 REAL_VALUE_TYPE *r;
471 unsigned int n;
473 int i, w = n / HOST_BITS_PER_LONG;
475 for (i = 0; i < w; ++i)
476 r->sig[i] = 0;
478 r->sig[w] &= ~(((unsigned long)1 << (n % HOST_BITS_PER_LONG)) - 1);
481 /* Divide the significands of A and B, placing the result in R. Return
482 true if the division was inexact. */
484 static inline bool
485 div_significands (r, a, b)
486 REAL_VALUE_TYPE *r;
487 const REAL_VALUE_TYPE *a, *b;
489 REAL_VALUE_TYPE u;
490 int i, bit = SIGNIFICAND_BITS - 1;
491 unsigned long msb, inexact;
493 u = *a;
494 memset (r->sig, 0, sizeof (r->sig));
496 msb = 0;
497 goto start;
500 msb = u.sig[SIGSZ-1] & SIG_MSB;
501 lshift_significand_1 (&u, &u);
502 start:
503 if (msb || cmp_significands (&u, b) >= 0)
505 sub_significands (&u, &u, b, 0);
506 set_significand_bit (r, bit);
509 while (--bit >= 0);
511 for (i = 0, inexact = 0; i < SIGSZ; i++)
512 inexact |= u.sig[i];
514 return inexact != 0;
517 /* Adjust the exponent and significand of R such that the most
518 significant bit is set. We underflow to zero and overflow to
519 infinity here, without denormals. (The intermediate representation
520 exponent is large enough to handle target denormals normalized.) */
522 static void
523 normalize (r)
524 REAL_VALUE_TYPE *r;
526 int shift = 0, exp;
527 int i, j;
529 /* Find the first word that is nonzero. */
530 for (i = SIGSZ - 1; i >= 0; i--)
531 if (r->sig[i] == 0)
532 shift += HOST_BITS_PER_LONG;
533 else
534 break;
536 /* Zero significand flushes to zero. */
537 if (i < 0)
539 r->class = rvc_zero;
540 r->exp = 0;
541 return;
544 /* Find the first bit that is nonzero. */
545 for (j = 0; ; j++)
546 if (r->sig[i] & ((unsigned long)1 << (HOST_BITS_PER_LONG - 1 - j)))
547 break;
548 shift += j;
550 if (shift > 0)
552 exp = r->exp - shift;
553 if (exp > MAX_EXP)
554 get_inf (r, r->sign);
555 else if (exp < -MAX_EXP)
556 get_zero (r, r->sign);
557 else
559 r->exp = exp;
560 lshift_significand (r, r, shift);
565 /* Return R = A + (SUBTRACT_P ? -B : B). */
567 static void
568 do_add (r, a, b, subtract_p)
569 REAL_VALUE_TYPE *r;
570 const REAL_VALUE_TYPE *a, *b;
571 int subtract_p;
573 int dexp, sign, exp;
574 REAL_VALUE_TYPE t;
575 bool inexact = false;
577 /* Determine if we need to add or subtract. */
578 sign = a->sign;
579 subtract_p = (sign ^ b->sign) ^ subtract_p;
581 switch (CLASS2 (a->class, b->class))
583 case CLASS2 (rvc_zero, rvc_zero):
584 /* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0. */
585 get_zero (r, sign & !subtract_p);
586 return;
588 case CLASS2 (rvc_zero, rvc_normal):
589 case CLASS2 (rvc_zero, rvc_inf):
590 case CLASS2 (rvc_zero, rvc_nan):
591 /* 0 + ANY = ANY. */
592 case CLASS2 (rvc_normal, rvc_nan):
593 case CLASS2 (rvc_inf, rvc_nan):
594 case CLASS2 (rvc_nan, rvc_nan):
595 /* ANY + NaN = NaN. */
596 case CLASS2 (rvc_normal, rvc_inf):
597 /* R + Inf = Inf. */
598 *r = *b;
599 r->sign = sign ^ subtract_p;
600 return;
602 case CLASS2 (rvc_normal, rvc_zero):
603 case CLASS2 (rvc_inf, rvc_zero):
604 case CLASS2 (rvc_nan, rvc_zero):
605 /* ANY + 0 = ANY. */
606 case CLASS2 (rvc_nan, rvc_normal):
607 case CLASS2 (rvc_nan, rvc_inf):
608 /* NaN + ANY = NaN. */
609 case CLASS2 (rvc_inf, rvc_normal):
610 /* Inf + R = Inf. */
611 *r = *a;
612 return;
614 case CLASS2 (rvc_inf, rvc_inf):
615 if (subtract_p)
616 /* Inf - Inf = NaN. */
617 get_canonical_qnan (r, 0);
618 else
619 /* Inf + Inf = Inf. */
620 *r = *a;
621 return;
623 case CLASS2 (rvc_normal, rvc_normal):
624 break;
626 default:
627 abort ();
630 /* Swap the arguments such that A has the larger exponent. */
631 dexp = a->exp - b->exp;
632 if (dexp < 0)
634 const REAL_VALUE_TYPE *t;
635 t = a, a = b, b = t;
636 dexp = -dexp;
637 sign ^= subtract_p;
639 exp = a->exp;
641 /* If the exponents are not identical, we need to shift the
642 significand of B down. */
643 if (dexp > 0)
645 /* If the exponents are too far apart, the significands
646 do not overlap, which makes the subtraction a noop. */
647 if (dexp >= SIGNIFICAND_BITS)
649 *r = *a;
650 r->sign = sign;
651 return;
654 inexact |= sticky_rshift_significand (&t, b, dexp);
655 b = &t;
658 if (subtract_p)
660 if (sub_significands (r, a, b, inexact))
662 /* We got a borrow out of the subtraction. That means that
663 A and B had the same exponent, and B had the larger
664 significand. We need to swap the sign and negate the
665 significand. */
666 sign ^= 1;
667 neg_significand (r, r);
670 else
672 if (add_significands (r, a, b))
674 /* We got carry out of the addition. This means we need to
675 shift the significand back down one bit and increase the
676 exponent. */
677 inexact |= sticky_rshift_significand (r, r, 1);
678 r->sig[SIGSZ-1] |= SIG_MSB;
679 if (++exp > MAX_EXP)
681 get_inf (r, sign);
682 return;
687 r->class = rvc_normal;
688 r->sign = sign;
689 r->exp = exp;
691 /* Re-normalize the result. */
692 normalize (r);
694 /* Special case: if the subtraction results in zero, the result
695 is positive. */
696 if (r->class == rvc_zero)
697 r->sign = 0;
698 else
699 r->sig[0] |= inexact;
702 /* Return R = A * B. */
704 static void
705 do_multiply (r, a, b)
706 REAL_VALUE_TYPE *r;
707 const REAL_VALUE_TYPE *a, *b;
709 REAL_VALUE_TYPE u, t, *rr;
710 unsigned int i, j, k;
711 int sign = a->sign ^ b->sign;
713 switch (CLASS2 (a->class, b->class))
715 case CLASS2 (rvc_zero, rvc_zero):
716 case CLASS2 (rvc_zero, rvc_normal):
717 case CLASS2 (rvc_normal, rvc_zero):
718 /* +-0 * ANY = 0 with appropriate sign. */
719 get_zero (r, sign);
720 return;
722 case CLASS2 (rvc_zero, rvc_nan):
723 case CLASS2 (rvc_normal, rvc_nan):
724 case CLASS2 (rvc_inf, rvc_nan):
725 case CLASS2 (rvc_nan, rvc_nan):
726 /* ANY * NaN = NaN. */
727 *r = *b;
728 r->sign = sign;
729 return;
731 case CLASS2 (rvc_nan, rvc_zero):
732 case CLASS2 (rvc_nan, rvc_normal):
733 case CLASS2 (rvc_nan, rvc_inf):
734 /* NaN * ANY = NaN. */
735 *r = *a;
736 r->sign = sign;
737 return;
739 case CLASS2 (rvc_zero, rvc_inf):
740 case CLASS2 (rvc_inf, rvc_zero):
741 /* 0 * Inf = NaN */
742 get_canonical_qnan (r, sign);
743 return;
745 case CLASS2 (rvc_inf, rvc_inf):
746 case CLASS2 (rvc_normal, rvc_inf):
747 case CLASS2 (rvc_inf, rvc_normal):
748 /* Inf * Inf = Inf, R * Inf = Inf */
749 overflow:
750 get_inf (r, sign);
751 return;
753 case CLASS2 (rvc_normal, rvc_normal):
754 break;
756 default:
757 abort ();
760 if (r == a || r == b)
761 rr = &t;
762 else
763 rr = r;
764 get_zero (rr, 0);
766 /* Collect all the partial products. Since we don't have sure access
767 to a widening multiply, we split each long into two half-words.
769 Consider the long-hand form of a four half-word multiplication:
771 A B C D
772 * E F G H
773 --------------
774 DE DF DG DH
775 CE CF CG CH
776 BE BF BG BH
777 AE AF AG AH
779 We construct partial products of the widened half-word products
780 that are known to not overlap, e.g. DF+DH. Each such partial
781 product is given its proper exponent, which allows us to sum them
782 and obtain the finished product. */
784 for (i = 0; i < SIGSZ * 2; ++i)
786 unsigned long ai = a->sig[i / 2];
787 if (i & 1)
788 ai >>= HOST_BITS_PER_LONG / 2;
789 else
790 ai &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
792 if (ai == 0)
793 continue;
795 for (j = 0; j < 2; ++j)
797 int exp = (a->exp - (2*SIGSZ-1-i)*(HOST_BITS_PER_LONG/2)
798 + (b->exp - (1-j)*(HOST_BITS_PER_LONG/2)));
800 if (exp > MAX_EXP)
801 goto overflow;
802 if (exp < -MAX_EXP)
803 /* Would underflow to zero, which we shouldn't bother adding. */
804 continue;
806 u.class = rvc_normal;
807 u.sign = 0;
808 u.exp = exp;
810 for (k = j; k < SIGSZ * 2; k += 2)
812 unsigned long bi = b->sig[k / 2];
813 if (k & 1)
814 bi >>= HOST_BITS_PER_LONG / 2;
815 else
816 bi &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
818 u.sig[k / 2] = ai * bi;
821 normalize (&u);
822 do_add (rr, rr, &u, 0);
826 rr->sign = sign;
827 if (rr != r)
828 *r = t;
831 /* Return R = A / B. */
833 static void
834 do_divide (r, a, b)
835 REAL_VALUE_TYPE *r;
836 const REAL_VALUE_TYPE *a, *b;
838 int exp, sign = a->sign ^ b->sign;
839 REAL_VALUE_TYPE t, *rr;
840 bool inexact;
842 switch (CLASS2 (a->class, b->class))
844 case CLASS2 (rvc_zero, rvc_zero):
845 /* 0 / 0 = NaN. */
846 case CLASS2 (rvc_inf, rvc_inf):
847 /* Inf / Inf = NaN. */
848 get_canonical_qnan (r, sign);
849 return;
851 case CLASS2 (rvc_zero, rvc_normal):
852 case CLASS2 (rvc_zero, rvc_inf):
853 /* 0 / ANY = 0. */
854 case CLASS2 (rvc_normal, rvc_inf):
855 /* R / Inf = 0. */
856 underflow:
857 get_zero (r, sign);
858 return;
860 case CLASS2 (rvc_normal, rvc_zero):
861 /* R / 0 = Inf. */
862 case CLASS2 (rvc_inf, rvc_zero):
863 /* Inf / 0 = Inf. */
864 get_inf (r, sign);
865 return;
867 case CLASS2 (rvc_zero, rvc_nan):
868 case CLASS2 (rvc_normal, rvc_nan):
869 case CLASS2 (rvc_inf, rvc_nan):
870 case CLASS2 (rvc_nan, rvc_nan):
871 /* ANY / NaN = NaN. */
872 *r = *b;
873 r->sign = sign;
874 return;
876 case CLASS2 (rvc_nan, rvc_zero):
877 case CLASS2 (rvc_nan, rvc_normal):
878 case CLASS2 (rvc_nan, rvc_inf):
879 /* NaN / ANY = NaN. */
880 *r = *a;
881 r->sign = sign;
882 return;
884 case CLASS2 (rvc_inf, rvc_normal):
885 /* Inf / R = Inf. */
886 overflow:
887 get_inf (r, sign);
888 return;
890 case CLASS2 (rvc_normal, rvc_normal):
891 break;
893 default:
894 abort ();
897 if (r == a || r == b)
898 rr = &t;
899 else
900 rr = r;
902 rr->class = rvc_normal;
903 rr->sign = sign;
905 exp = a->exp - b->exp + 1;
906 if (exp > MAX_EXP)
907 goto overflow;
908 if (exp < -MAX_EXP)
909 goto underflow;
910 rr->exp = exp;
912 inexact = div_significands (rr, a, b);
914 /* Re-normalize the result. */
915 normalize (rr);
916 rr->sig[0] |= inexact;
918 if (rr != r)
919 *r = t;
922 /* Return a tri-state comparison of A vs B. Return NAN_RESULT if
923 one of the two operands is a NaN. */
925 static int
926 do_compare (a, b, nan_result)
927 const REAL_VALUE_TYPE *a, *b;
928 int nan_result;
930 int ret;
932 switch (CLASS2 (a->class, b->class))
934 case CLASS2 (rvc_zero, rvc_zero):
935 /* Sign of zero doesn't matter for compares. */
936 return 0;
938 case CLASS2 (rvc_inf, rvc_zero):
939 case CLASS2 (rvc_inf, rvc_normal):
940 case CLASS2 (rvc_normal, rvc_zero):
941 return (a->sign ? -1 : 1);
943 case CLASS2 (rvc_inf, rvc_inf):
944 return -a->sign - -b->sign;
946 case CLASS2 (rvc_zero, rvc_normal):
947 case CLASS2 (rvc_zero, rvc_inf):
948 case CLASS2 (rvc_normal, rvc_inf):
949 return (b->sign ? 1 : -1);
951 case CLASS2 (rvc_zero, rvc_nan):
952 case CLASS2 (rvc_normal, rvc_nan):
953 case CLASS2 (rvc_inf, rvc_nan):
954 case CLASS2 (rvc_nan, rvc_nan):
955 case CLASS2 (rvc_nan, rvc_zero):
956 case CLASS2 (rvc_nan, rvc_normal):
957 case CLASS2 (rvc_nan, rvc_inf):
958 return nan_result;
960 case CLASS2 (rvc_normal, rvc_normal):
961 break;
963 default:
964 abort ();
967 if (a->sign != b->sign)
968 return -a->sign - -b->sign;
970 if (a->exp > b->exp)
971 ret = 1;
972 else if (a->exp < b->exp)
973 ret = -1;
974 else
975 ret = cmp_significands (a, b);
977 return (a->sign ? -ret : ret);
980 /* Return A truncated to an integral value toward zero. */
982 static void
983 do_fix_trunc (r, a)
984 REAL_VALUE_TYPE *r;
985 const REAL_VALUE_TYPE *a;
987 *r = *a;
989 switch (r->class)
991 case rvc_zero:
992 case rvc_inf:
993 case rvc_nan:
994 break;
996 case rvc_normal:
997 if (r->exp <= 0)
998 get_zero (r, r->sign);
999 else if (r->exp < SIGNIFICAND_BITS)
1000 clear_significand_below (r, SIGNIFICAND_BITS - r->exp);
1001 break;
1003 default:
1004 abort ();
1008 /* Perform the binary or unary operation described by CODE.
1009 For a unary operation, leave OP1 NULL. */
1011 void
1012 real_arithmetic (r, icode, op0, op1)
1013 REAL_VALUE_TYPE *r;
1014 int icode;
1015 const REAL_VALUE_TYPE *op0, *op1;
1017 enum tree_code code = icode;
1019 switch (code)
1021 case PLUS_EXPR:
1022 do_add (r, op0, op1, 0);
1023 break;
1025 case MINUS_EXPR:
1026 do_add (r, op0, op1, 1);
1027 break;
1029 case MULT_EXPR:
1030 do_multiply (r, op0, op1);
1031 break;
1033 case RDIV_EXPR:
1034 do_divide (r, op0, op1);
1035 break;
1037 case MIN_EXPR:
1038 if (op1->class == rvc_nan)
1039 *r = *op1;
1040 else if (do_compare (op0, op1, -1) < 0)
1041 *r = *op0;
1042 else
1043 *r = *op1;
1044 break;
1046 case MAX_EXPR:
1047 if (op1->class == rvc_nan)
1048 *r = *op1;
1049 else if (do_compare (op0, op1, 1) < 0)
1050 *r = *op1;
1051 else
1052 *r = *op0;
1053 break;
1055 case NEGATE_EXPR:
1056 *r = *op0;
1057 r->sign ^= 1;
1058 break;
1060 case ABS_EXPR:
1061 *r = *op0;
1062 r->sign = 0;
1063 break;
1065 case FIX_TRUNC_EXPR:
1066 do_fix_trunc (r, op0);
1067 break;
1069 default:
1070 abort ();
1074 /* Legacy. Similar, but return the result directly. */
1076 REAL_VALUE_TYPE
1077 real_arithmetic2 (icode, op0, op1)
1078 int icode;
1079 const REAL_VALUE_TYPE *op0, *op1;
1081 REAL_VALUE_TYPE r;
1082 real_arithmetic (&r, icode, op0, op1);
1083 return r;
1086 bool
1087 real_compare (icode, op0, op1)
1088 int icode;
1089 const REAL_VALUE_TYPE *op0, *op1;
1091 enum tree_code code = icode;
1093 switch (code)
1095 case LT_EXPR:
1096 return do_compare (op0, op1, 1) < 0;
1097 case LE_EXPR:
1098 return do_compare (op0, op1, 1) <= 0;
1099 case GT_EXPR:
1100 return do_compare (op0, op1, -1) > 0;
1101 case GE_EXPR:
1102 return do_compare (op0, op1, -1) >= 0;
1103 case EQ_EXPR:
1104 return do_compare (op0, op1, -1) == 0;
1105 case NE_EXPR:
1106 return do_compare (op0, op1, -1) != 0;
1107 case UNORDERED_EXPR:
1108 return op0->class == rvc_nan || op1->class == rvc_nan;
1109 case ORDERED_EXPR:
1110 return op0->class != rvc_nan && op1->class != rvc_nan;
1111 case UNLT_EXPR:
1112 return do_compare (op0, op1, -1) < 0;
1113 case UNLE_EXPR:
1114 return do_compare (op0, op1, -1) <= 0;
1115 case UNGT_EXPR:
1116 return do_compare (op0, op1, 1) > 0;
1117 case UNGE_EXPR:
1118 return do_compare (op0, op1, 1) >= 0;
1119 case UNEQ_EXPR:
1120 return do_compare (op0, op1, 0) == 0;
1122 default:
1123 abort ();
1127 /* Return floor log2(R). */
1130 real_exponent (r)
1131 const REAL_VALUE_TYPE *r;
1133 switch (r->class)
1135 case rvc_zero:
1136 return 0;
1137 case rvc_inf:
1138 case rvc_nan:
1139 return (unsigned int)-1 >> 1;
1140 case rvc_normal:
1141 return r->exp;
1142 default:
1143 abort ();
1147 /* R = OP0 * 2**EXP. */
1149 void
1150 real_ldexp (r, op0, exp)
1151 REAL_VALUE_TYPE *r;
1152 const REAL_VALUE_TYPE *op0;
1153 int exp;
1155 *r = *op0;
1156 switch (r->class)
1158 case rvc_zero:
1159 case rvc_inf:
1160 case rvc_nan:
1161 break;
1163 case rvc_normal:
1164 exp += op0->exp;
1165 if (exp > MAX_EXP)
1166 get_inf (r, r->sign);
1167 else if (exp < -MAX_EXP)
1168 get_zero (r, r->sign);
1169 else
1170 r->exp = exp;
1171 break;
1173 default:
1174 abort ();
1178 /* Determine whether a floating-point value X is infinite. */
1180 bool
1181 real_isinf (r)
1182 const REAL_VALUE_TYPE *r;
1184 return (r->class == rvc_inf);
1187 /* Determine whether a floating-point value X is a NaN. */
1189 bool
1190 real_isnan (r)
1191 const REAL_VALUE_TYPE *r;
1193 return (r->class == rvc_nan);
1196 /* Determine whether a floating-point value X is negative. */
1198 bool
1199 real_isneg (r)
1200 const REAL_VALUE_TYPE *r;
1202 return r->sign;
1205 /* Determine whether a floating-point value X is minus zero. */
1207 bool
1208 real_isnegzero (r)
1209 const REAL_VALUE_TYPE *r;
1211 return r->sign && r->class == rvc_zero;
1214 /* Compare two floating-point objects for bitwise identity. */
1216 extern bool
1217 real_identical (a, b)
1218 const REAL_VALUE_TYPE *a, *b;
1220 int i;
1222 if (a->class != b->class)
1223 return false;
1224 if (a->sign != b->sign)
1225 return false;
1227 switch (a->class)
1229 case rvc_zero:
1230 case rvc_inf:
1231 break;
1233 case rvc_normal:
1234 if (a->exp != b->exp)
1235 return false;
1236 /* FALLTHRU */
1237 case rvc_nan:
1238 for (i = 0; i < SIGSZ; ++i)
1239 if (a->sig[i] != b->sig[i])
1240 return false;
1241 break;
1243 default:
1244 abort ();
1247 return true;
1250 /* Try to change R into its exact multiplicative inverse in machine
1251 mode MODE. Return true if successful. */
1253 bool
1254 exact_real_inverse (mode, r)
1255 enum machine_mode mode;
1256 REAL_VALUE_TYPE *r;
1258 const REAL_VALUE_TYPE *one = real_digit (1);
1259 REAL_VALUE_TYPE u;
1260 int i;
1262 if (r->class != rvc_normal)
1263 return false;
1265 /* Check for a power of two: all significand bits zero except the MSB. */
1266 for (i = 0; i < SIGSZ-1; ++i)
1267 if (r->sig[i] != 0)
1268 return false;
1269 if (r->sig[SIGSZ-1] != SIG_MSB)
1270 return false;
1272 /* Find the inverse and truncate to the required mode. */
1273 do_divide (&u, one, r);
1274 real_convert (&u, mode, &u);
1276 /* The rounding may have overflowed. */
1277 if (u.class != rvc_normal)
1278 return false;
1279 for (i = 0; i < SIGSZ-1; ++i)
1280 if (u.sig[i] != 0)
1281 return false;
1282 if (u.sig[SIGSZ-1] != SIG_MSB)
1283 return false;
1285 *r = u;
1286 return true;
1289 /* Render R as an integer. */
1291 HOST_WIDE_INT
1292 real_to_integer (r)
1293 const REAL_VALUE_TYPE *r;
1295 unsigned HOST_WIDE_INT i;
1297 switch (r->class)
1299 case rvc_zero:
1300 underflow:
1301 return 0;
1303 case rvc_inf:
1304 case rvc_nan:
1305 overflow:
1306 i = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1307 if (!r->sign)
1308 i--;
1309 return i;
1311 case rvc_normal:
1312 if (r->exp <= 0)
1313 goto underflow;
1314 /* Only force overflow for unsigned overflow. Signed overflow is
1315 undefined, so it doesn't matter what we return, and some callers
1316 expect to be able to use this routine for both signed and
1317 unsigned conversions. */
1318 if (r->exp > HOST_BITS_PER_WIDE_INT)
1319 goto overflow;
1321 if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1322 i = r->sig[SIGSZ-1];
1323 else if (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG)
1325 i = r->sig[SIGSZ-1];
1326 i = i << (HOST_BITS_PER_LONG - 1) << 1;
1327 i |= r->sig[SIGSZ-2];
1329 else
1330 abort ();
1332 i >>= HOST_BITS_PER_WIDE_INT - r->exp;
1334 if (r->sign)
1335 i = -i;
1336 return i;
1338 default:
1339 abort ();
1343 /* Likewise, but to an integer pair, HI+LOW. */
1345 void
1346 real_to_integer2 (plow, phigh, r)
1347 HOST_WIDE_INT *plow, *phigh;
1348 const REAL_VALUE_TYPE *r;
1350 REAL_VALUE_TYPE t;
1351 HOST_WIDE_INT low, high;
1352 int exp;
1354 switch (r->class)
1356 case rvc_zero:
1357 underflow:
1358 low = high = 0;
1359 break;
1361 case rvc_inf:
1362 case rvc_nan:
1363 overflow:
1364 high = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1365 if (r->sign)
1366 low = 0;
1367 else
1369 high--;
1370 low = -1;
1372 break;
1374 case rvc_normal:
1375 exp = r->exp;
1376 if (exp <= 0)
1377 goto underflow;
1378 /* Only force overflow for unsigned overflow. Signed overflow is
1379 undefined, so it doesn't matter what we return, and some callers
1380 expect to be able to use this routine for both signed and
1381 unsigned conversions. */
1382 if (exp > 2*HOST_BITS_PER_WIDE_INT)
1383 goto overflow;
1385 rshift_significand (&t, r, 2*HOST_BITS_PER_WIDE_INT - exp);
1386 if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1388 high = t.sig[SIGSZ-1];
1389 low = t.sig[SIGSZ-2];
1391 else if (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG)
1393 high = t.sig[SIGSZ-1];
1394 high = high << (HOST_BITS_PER_LONG - 1) << 1;
1395 high |= t.sig[SIGSZ-2];
1397 low = t.sig[SIGSZ-3];
1398 low = low << (HOST_BITS_PER_LONG - 1) << 1;
1399 low |= t.sig[SIGSZ-4];
1401 else
1402 abort ();
1404 if (r->sign)
1406 if (low == 0)
1407 high = -high;
1408 else
1409 low = -low, high = ~high;
1411 break;
1413 default:
1414 abort ();
1417 *plow = low;
1418 *phigh = high;
1421 /* A subroutine of real_to_decimal. Compute the quotient and remainder
1422 of NUM / DEN. Return the quotient and place the remainder in NUM.
1423 It is expected that NUM / DEN are close enough that the quotient is
1424 small. */
1426 static unsigned long
1427 rtd_divmod (num, den)
1428 REAL_VALUE_TYPE *num, *den;
1430 unsigned long q, msb;
1431 int expn = num->exp, expd = den->exp;
1433 if (expn < expd)
1434 return 0;
1436 q = msb = 0;
1437 goto start;
1440 msb = num->sig[SIGSZ-1] & SIG_MSB;
1441 q <<= 1;
1442 lshift_significand_1 (num, num);
1443 start:
1444 if (msb || cmp_significands (num, den) >= 0)
1446 sub_significands (num, num, den, 0);
1447 q |= 1;
1450 while (--expn >= expd);
1452 num->exp = expd;
1453 normalize (num);
1455 return q;
1458 /* Render R as a decimal floating point constant. Emit DIGITS significant
1459 digits in the result, bounded by BUF_SIZE. If DIGITS is 0, choose the
1460 maximum for the representation. If CROP_TRAILING_ZEROS, strip trailing
1461 zeros. */
1463 #define M_LOG10_2 0.30102999566398119521
1465 void
1466 real_to_decimal (str, r_orig, buf_size, digits, crop_trailing_zeros)
1467 char *str;
1468 const REAL_VALUE_TYPE *r_orig;
1469 size_t buf_size, digits;
1470 int crop_trailing_zeros;
1472 const REAL_VALUE_TYPE *one, *ten;
1473 REAL_VALUE_TYPE r, pten, u, v;
1474 int dec_exp, cmp_one, digit;
1475 size_t max_digits;
1476 char *p, *first, *last;
1477 bool sign;
1479 r = *r_orig;
1480 switch (r.class)
1482 case rvc_zero:
1483 strcpy (str, (r.sign ? "-0.0" : "0.0"));
1484 return;
1485 case rvc_normal:
1486 break;
1487 case rvc_inf:
1488 strcpy (str, (r.sign ? "-Inf" : "+Inf"));
1489 return;
1490 case rvc_nan:
1491 /* ??? Print the significand as well, if not canonical? */
1492 strcpy (str, (r.sign ? "-NaN" : "+NaN"));
1493 return;
1494 default:
1495 abort ();
1498 /* Bound the number of digits printed by the size of the representation. */
1499 max_digits = SIGNIFICAND_BITS * M_LOG10_2;
1500 if (digits == 0 || digits > max_digits)
1501 digits = max_digits;
1503 /* Estimate the decimal exponent, and compute the length of the string it
1504 will print as. Be conservative and add one to account for possible
1505 overflow or rounding error. */
1506 dec_exp = r.exp * M_LOG10_2;
1507 for (max_digits = 1; dec_exp ; max_digits++)
1508 dec_exp /= 10;
1510 /* Bound the number of digits printed by the size of the output buffer. */
1511 max_digits = buf_size - 1 - 1 - 2 - max_digits - 1;
1512 if (max_digits > buf_size)
1513 abort ();
1514 if (digits > max_digits)
1515 digits = max_digits;
1517 one = real_digit (1);
1518 ten = ten_to_ptwo (0);
1520 sign = r.sign;
1521 r.sign = 0;
1523 dec_exp = 0;
1524 pten = *one;
1526 cmp_one = do_compare (&r, one, 0);
1527 if (cmp_one > 0)
1529 int m;
1531 /* Number is greater than one. Convert significand to an integer
1532 and strip trailing decimal zeros. */
1534 u = r;
1535 u.exp = SIGNIFICAND_BITS - 1;
1537 /* Largest M, such that 10**2**M fits within SIGNIFICAND_BITS. */
1538 m = floor_log2 (max_digits);
1540 /* Iterate over the bits of the possible powers of 10 that might
1541 be present in U and eliminate them. That is, if we find that
1542 10**2**M divides U evenly, keep the division and increase
1543 DEC_EXP by 2**M. */
1546 REAL_VALUE_TYPE t;
1548 do_divide (&t, &u, ten_to_ptwo (m));
1549 do_fix_trunc (&v, &t);
1550 if (cmp_significands (&v, &t) == 0)
1552 u = t;
1553 dec_exp += 1 << m;
1556 while (--m >= 0);
1558 /* Revert the scaling to integer that we performed earlier. */
1559 u.exp += r.exp - (SIGNIFICAND_BITS - 1);
1560 r = u;
1562 /* Find power of 10. Do this by dividing out 10**2**M when
1563 this is larger than the current remainder. Fill PTEN with
1564 the power of 10 that we compute. */
1565 if (r.exp > 0)
1567 m = floor_log2 ((int)(r.exp * M_LOG10_2)) + 1;
1570 const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1571 if (do_compare (&u, ptentwo, 0) >= 0)
1573 do_divide (&u, &u, ptentwo);
1574 do_multiply (&pten, &pten, ptentwo);
1575 dec_exp += 1 << m;
1578 while (--m >= 0);
1580 else
1581 /* We managed to divide off enough tens in the above reduction
1582 loop that we've now got a negative exponent. Fall into the
1583 less-than-one code to compute the proper value for PTEN. */
1584 cmp_one = -1;
1586 if (cmp_one < 0)
1588 int m;
1590 /* Number is less than one. Pad significand with leading
1591 decimal zeros. */
1593 v = r;
1594 while (1)
1596 /* Stop if we'd shift bits off the bottom. */
1597 if (v.sig[0] & 7)
1598 break;
1600 do_multiply (&u, &v, ten);
1602 /* Stop if we're now >= 1. */
1603 if (u.exp > 0)
1604 break;
1606 v = u;
1607 dec_exp -= 1;
1609 r = v;
1611 /* Find power of 10. Do this by multiplying in P=10**2**M when
1612 the current remainder is smaller than 1/P. Fill PTEN with the
1613 power of 10 that we compute. */
1614 m = floor_log2 ((int)(-r.exp * M_LOG10_2)) + 1;
1617 const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1618 const REAL_VALUE_TYPE *ptenmtwo = ten_to_mptwo (m);
1620 if (do_compare (&v, ptenmtwo, 0) <= 0)
1622 do_multiply (&v, &v, ptentwo);
1623 do_multiply (&pten, &pten, ptentwo);
1624 dec_exp -= 1 << m;
1627 while (--m >= 0);
1629 /* Invert the positive power of 10 that we've collected so far. */
1630 do_divide (&pten, one, &pten);
1633 p = str;
1634 if (sign)
1635 *p++ = '-';
1636 first = p++;
1638 /* At this point, PTEN should contain the nearest power of 10 smaller
1639 than R, such that this division produces the first digit.
1641 Using a divide-step primitive that returns the complete integral
1642 remainder avoids the rounding error that would be produced if
1643 we were to use do_divide here and then simply multiply by 10 for
1644 each subsequent digit. */
1646 digit = rtd_divmod (&r, &pten);
1648 /* Be prepared for error in that division via underflow ... */
1649 if (digit == 0 && cmp_significand_0 (&r))
1651 /* Multiply by 10 and try again. */
1652 do_multiply (&r, &r, ten);
1653 digit = rtd_divmod (&r, &pten);
1654 dec_exp -= 1;
1655 if (digit == 0)
1656 abort ();
1659 /* ... or overflow. */
1660 if (digit == 10)
1662 *p++ = '1';
1663 if (--digits > 0)
1664 *p++ = '0';
1665 dec_exp += 1;
1667 else if (digit > 10)
1668 abort ();
1669 else
1670 *p++ = digit + '0';
1672 /* Generate subsequent digits. */
1673 while (--digits > 0)
1675 do_multiply (&r, &r, ten);
1676 digit = rtd_divmod (&r, &pten);
1677 *p++ = digit + '0';
1679 last = p;
1681 /* Generate one more digit with which to do rounding. */
1682 do_multiply (&r, &r, ten);
1683 digit = rtd_divmod (&r, &pten);
1685 /* Round the result. */
1686 if (digit == 5)
1688 /* Round to nearest. If R is nonzero there are additional
1689 nonzero digits to be extracted. */
1690 if (cmp_significand_0 (&r))
1691 digit++;
1692 /* Round to even. */
1693 else if ((p[-1] - '0') & 1)
1694 digit++;
1696 if (digit > 5)
1698 while (p > first)
1700 digit = *--p;
1701 if (digit == '9')
1702 *p = '0';
1703 else
1705 *p = digit + 1;
1706 break;
1710 /* Carry out of the first digit. This means we had all 9's and
1711 now have all 0's. "Prepend" a 1 by overwriting the first 0. */
1712 if (p == first)
1714 first[1] = '1';
1715 dec_exp++;
1719 /* Insert the decimal point. */
1720 first[0] = first[1];
1721 first[1] = '.';
1723 /* If requested, drop trailing zeros. Never crop past "1.0". */
1724 if (crop_trailing_zeros)
1725 while (last > first + 3 && last[-1] == '0')
1726 last--;
1728 /* Append the exponent. */
1729 sprintf (last, "e%+d", dec_exp);
1732 /* Render R as a hexadecimal floating point constant. Emit DIGITS
1733 significant digits in the result, bounded by BUF_SIZE. If DIGITS is 0,
1734 choose the maximum for the representation. If CROP_TRAILING_ZEROS,
1735 strip trailing zeros. */
1737 void
1738 real_to_hexadecimal (str, r, buf_size, digits, crop_trailing_zeros)
1739 char *str;
1740 const REAL_VALUE_TYPE *r;
1741 size_t buf_size, digits;
1742 int crop_trailing_zeros;
1744 int i, j, exp = r->exp;
1745 char *p, *first;
1746 char exp_buf[16];
1747 size_t max_digits;
1749 switch (r->class)
1751 case rvc_zero:
1752 exp = 0;
1753 break;
1754 case rvc_normal:
1755 break;
1756 case rvc_inf:
1757 strcpy (str, (r->sign ? "-Inf" : "+Inf"));
1758 return;
1759 case rvc_nan:
1760 /* ??? Print the significand as well, if not canonical? */
1761 strcpy (str, (r->sign ? "-NaN" : "+NaN"));
1762 return;
1763 default:
1764 abort ();
1767 if (digits == 0)
1768 digits = SIGNIFICAND_BITS / 4;
1770 /* Bound the number of digits printed by the size of the output buffer. */
1772 sprintf (exp_buf, "p%+d", exp);
1773 max_digits = buf_size - strlen (exp_buf) - r->sign - 4 - 1;
1774 if (max_digits > buf_size)
1775 abort ();
1776 if (digits > max_digits)
1777 digits = max_digits;
1779 p = str;
1780 if (r->sign)
1781 *p++ = '-';
1782 *p++ = '0';
1783 *p++ = 'x';
1784 *p++ = '0';
1785 *p++ = '.';
1786 first = p;
1788 for (i = SIGSZ - 1; i >= 0; --i)
1789 for (j = HOST_BITS_PER_LONG - 4; j >= 0; j -= 4)
1791 *p++ = "0123456789abcdef"[(r->sig[i] >> j) & 15];
1792 if (--digits == 0)
1793 goto out;
1796 out:
1797 if (crop_trailing_zeros)
1798 while (p > first + 1 && p[-1] == '0')
1799 p--;
1801 sprintf (p, "p%+d", exp);
1804 /* Initialize R from a decimal or hexadecimal string. The string is
1805 assumed to have been syntax checked already. */
1807 void
1808 real_from_string (r, str)
1809 REAL_VALUE_TYPE *r;
1810 const char *str;
1812 int exp = 0;
1813 bool sign = false;
1815 get_zero (r, 0);
1817 if (*str == '-')
1819 sign = true;
1820 str++;
1822 else if (*str == '+')
1823 str++;
1825 if (str[0] == '0' && str[1] == 'x')
1827 /* Hexadecimal floating point. */
1828 int pos = SIGNIFICAND_BITS - 4, d;
1830 str += 2;
1832 while (*str == '0')
1833 str++;
1834 while (1)
1836 d = hex_value (*str);
1837 if (d == _hex_bad)
1838 break;
1839 if (pos >= 0)
1841 r->sig[pos / HOST_BITS_PER_LONG]
1842 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1843 pos -= 4;
1845 exp += 4;
1846 str++;
1848 if (*str == '.')
1850 str++;
1851 if (pos == SIGNIFICAND_BITS - 4)
1853 while (*str == '0')
1854 str++, exp -= 4;
1856 while (1)
1858 d = hex_value (*str);
1859 if (d == _hex_bad)
1860 break;
1861 if (pos >= 0)
1863 r->sig[pos / HOST_BITS_PER_LONG]
1864 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1865 pos -= 4;
1867 str++;
1870 if (*str == 'p' || *str == 'P')
1872 bool exp_neg = false;
1874 str++;
1875 if (*str == '-')
1877 exp_neg = true;
1878 str++;
1880 else if (*str == '+')
1881 str++;
1883 d = 0;
1884 while (ISDIGIT (*str))
1886 d *= 10;
1887 d += *str - '0';
1888 if (d > MAX_EXP)
1890 /* Overflowed the exponent. */
1891 if (exp_neg)
1892 goto underflow;
1893 else
1894 goto overflow;
1896 str++;
1898 if (exp_neg)
1899 d = -d;
1901 exp += d;
1904 r->class = rvc_normal;
1905 r->exp = exp;
1907 normalize (r);
1909 else
1911 /* Decimal floating point. */
1912 const REAL_VALUE_TYPE *ten = ten_to_ptwo (0);
1913 int d;
1915 while (*str == '0')
1916 str++;
1917 while (ISDIGIT (*str))
1919 d = *str++ - '0';
1920 do_multiply (r, r, ten);
1921 if (d)
1922 do_add (r, r, real_digit (d), 0);
1924 if (*str == '.')
1926 str++;
1927 if (r->class == rvc_zero)
1929 while (*str == '0')
1930 str++, exp--;
1932 while (ISDIGIT (*str))
1934 d = *str++ - '0';
1935 do_multiply (r, r, ten);
1936 if (d)
1937 do_add (r, r, real_digit (d), 0);
1938 exp--;
1942 if (*str == 'e' || *str == 'E')
1944 bool exp_neg = false;
1946 str++;
1947 if (*str == '-')
1949 exp_neg = true;
1950 str++;
1952 else if (*str == '+')
1953 str++;
1955 d = 0;
1956 while (ISDIGIT (*str))
1958 d *= 10;
1959 d += *str - '0';
1960 if (d > MAX_EXP)
1962 /* Overflowed the exponent. */
1963 if (exp_neg)
1964 goto underflow;
1965 else
1966 goto overflow;
1968 str++;
1970 if (exp_neg)
1971 d = -d;
1972 exp += d;
1975 if (exp)
1976 times_pten (r, exp);
1979 r->sign = sign;
1980 return;
1982 underflow:
1983 get_zero (r, sign);
1984 return;
1986 overflow:
1987 get_inf (r, sign);
1988 return;
1991 /* Legacy. Similar, but return the result directly. */
1993 REAL_VALUE_TYPE
1994 real_from_string2 (s, mode)
1995 const char *s;
1996 enum machine_mode mode;
1998 REAL_VALUE_TYPE r;
2000 real_from_string (&r, s);
2001 if (mode != VOIDmode)
2002 real_convert (&r, mode, &r);
2004 return r;
2007 /* Initialize R from the integer pair HIGH+LOW. */
2009 void
2010 real_from_integer (r, mode, low, high, unsigned_p)
2011 REAL_VALUE_TYPE *r;
2012 enum machine_mode mode;
2013 unsigned HOST_WIDE_INT low;
2014 HOST_WIDE_INT high;
2015 int unsigned_p;
2017 if (low == 0 && high == 0)
2018 get_zero (r, 0);
2019 else
2021 r->class = rvc_normal;
2022 r->sign = high < 0 && !unsigned_p;
2023 r->exp = 2 * HOST_BITS_PER_WIDE_INT;
2025 if (r->sign)
2027 high = ~high;
2028 if (low == 0)
2029 high += 1;
2030 else
2031 low = -low;
2034 if (HOST_BITS_PER_LONG == HOST_BITS_PER_WIDE_INT)
2036 r->sig[SIGSZ-1] = high;
2037 r->sig[SIGSZ-2] = low;
2038 memset (r->sig, 0, sizeof(long)*(SIGSZ-2));
2040 else if (HOST_BITS_PER_LONG*2 == HOST_BITS_PER_WIDE_INT)
2042 r->sig[SIGSZ-1] = high >> (HOST_BITS_PER_LONG - 1) >> 1;
2043 r->sig[SIGSZ-2] = high;
2044 r->sig[SIGSZ-3] = low >> (HOST_BITS_PER_LONG - 1) >> 1;
2045 r->sig[SIGSZ-4] = low;
2046 if (SIGSZ > 4)
2047 memset (r->sig, 0, sizeof(long)*(SIGSZ-4));
2049 else
2050 abort ();
2052 normalize (r);
2055 if (mode != VOIDmode)
2056 real_convert (r, mode, r);
2059 /* Returns 10**2**N. */
2061 static const REAL_VALUE_TYPE *
2062 ten_to_ptwo (n)
2063 int n;
2065 static REAL_VALUE_TYPE tens[EXP_BITS];
2067 if (n < 0 || n >= EXP_BITS)
2068 abort ();
2070 if (tens[n].class == rvc_zero)
2072 if (n < (HOST_BITS_PER_WIDE_INT == 64 ? 5 : 4))
2074 HOST_WIDE_INT t = 10;
2075 int i;
2077 for (i = 0; i < n; ++i)
2078 t *= t;
2080 real_from_integer (&tens[n], VOIDmode, t, 0, 1);
2082 else
2084 const REAL_VALUE_TYPE *t = ten_to_ptwo (n - 1);
2085 do_multiply (&tens[n], t, t);
2089 return &tens[n];
2092 /* Returns 10**(-2**N). */
2094 static const REAL_VALUE_TYPE *
2095 ten_to_mptwo (n)
2096 int n;
2098 static REAL_VALUE_TYPE tens[EXP_BITS];
2100 if (n < 0 || n >= EXP_BITS)
2101 abort ();
2103 if (tens[n].class == rvc_zero)
2104 do_divide (&tens[n], real_digit (1), ten_to_ptwo (n));
2106 return &tens[n];
2109 /* Returns N. */
2111 static const REAL_VALUE_TYPE *
2112 real_digit (n)
2113 int n;
2115 static REAL_VALUE_TYPE num[10];
2117 if (n < 0 || n > 9)
2118 abort ();
2120 if (n > 0 && num[n].class == rvc_zero)
2121 real_from_integer (&num[n], VOIDmode, n, 0, 1);
2123 return &num[n];
2126 /* Multiply R by 10**EXP. */
2128 static void
2129 times_pten (r, exp)
2130 REAL_VALUE_TYPE *r;
2131 int exp;
2133 REAL_VALUE_TYPE pten, *rr;
2134 bool negative = (exp < 0);
2135 int i;
2137 if (negative)
2139 exp = -exp;
2140 pten = *real_digit (1);
2141 rr = &pten;
2143 else
2144 rr = r;
2146 for (i = 0; exp > 0; ++i, exp >>= 1)
2147 if (exp & 1)
2148 do_multiply (rr, rr, ten_to_ptwo (i));
2150 if (negative)
2151 do_divide (r, r, &pten);
2154 /* Fills R with +Inf. */
2156 void
2157 real_inf (r)
2158 REAL_VALUE_TYPE *r;
2160 get_inf (r, 0);
2163 /* Fills R with a NaN whose significand is described by STR. If QUIET,
2164 we force a QNaN, else we force an SNaN. The string, if not empty,
2165 is parsed as a number and placed in the significand. Return true
2166 if the string was successfully parsed. */
2168 bool
2169 real_nan (r, str, quiet, mode)
2170 REAL_VALUE_TYPE *r;
2171 const char *str;
2172 int quiet;
2173 enum machine_mode mode;
2175 const struct real_format *fmt;
2177 fmt = real_format_for_mode[mode - QFmode];
2178 if (fmt == NULL)
2179 abort ();
2181 if (*str == 0)
2183 if (quiet)
2184 get_canonical_qnan (r, 0);
2185 else
2186 get_canonical_snan (r, 0);
2188 else
2190 int base = 10, d;
2191 bool neg = false;
2193 memset (r, 0, sizeof (*r));
2194 r->class = rvc_nan;
2196 /* Parse akin to strtol into the significand of R. */
2198 while (ISSPACE (*str))
2199 str++;
2200 if (*str == '-')
2201 str++, neg = true;
2202 else if (*str == '+')
2203 str++;
2204 if (*str == '0')
2206 if (*++str == 'x')
2207 str++, base = 16;
2208 else
2209 base = 8;
2212 while ((d = hex_value (*str)) < base)
2214 REAL_VALUE_TYPE u;
2216 switch (base)
2218 case 8:
2219 lshift_significand (r, r, 3);
2220 break;
2221 case 16:
2222 lshift_significand (r, r, 4);
2223 break;
2224 case 10:
2225 lshift_significand_1 (&u, r);
2226 lshift_significand (r, r, 3);
2227 add_significands (r, r, &u);
2228 break;
2229 default:
2230 abort ();
2233 get_zero (&u, 0);
2234 u.sig[0] = d;
2235 add_significands (r, r, &u);
2237 str++;
2240 /* Must have consumed the entire string for success. */
2241 if (*str != 0)
2242 return false;
2244 /* Shift the significand into place such that the bits
2245 are in the most significant bits for the format. */
2246 lshift_significand (r, r, SIGNIFICAND_BITS - fmt->p);
2248 /* Our MSB is always unset for NaNs. */
2249 r->sig[SIGSZ-1] &= ~SIG_MSB;
2251 /* Force quiet or signalling NaN. */
2252 if (quiet)
2253 r->sig[SIGSZ-1] |= SIG_MSB >> 1;
2254 else
2255 r->sig[SIGSZ-1] &= ~(SIG_MSB >> 1);
2257 /* Force at least one bit of the significand set. */
2258 for (d = 0; d < SIGSZ; ++d)
2259 if (r->sig[d])
2260 break;
2261 if (d == SIGSZ)
2262 r->sig[SIGSZ-1] |= SIG_MSB >> 2;
2264 /* Our intermediate format forces QNaNs to have MSB-1 set.
2265 If the target format has QNaNs with the top bit unset,
2266 mirror the output routines and invert the top two bits. */
2267 if (!fmt->qnan_msb_set)
2268 r->sig[SIGSZ-1] ^= (SIG_MSB >> 1) | (SIG_MSB >> 2);
2271 return true;
2274 /* Fills R with 2**N. */
2276 void
2277 real_2expN (r, n)
2278 REAL_VALUE_TYPE *r;
2279 int n;
2281 memset (r, 0, sizeof (*r));
2283 n++;
2284 if (n > MAX_EXP)
2285 r->class = rvc_inf;
2286 else if (n < -MAX_EXP)
2288 else
2290 r->class = rvc_normal;
2291 r->exp = n;
2292 r->sig[SIGSZ-1] = SIG_MSB;
2297 static void
2298 round_for_format (fmt, r)
2299 const struct real_format *fmt;
2300 REAL_VALUE_TYPE *r;
2302 int p2, np2, i, w;
2303 unsigned long sticky;
2304 bool guard, lsb;
2305 int emin2m1, emax2;
2307 p2 = fmt->p * fmt->log2_b;
2308 emin2m1 = (fmt->emin - 1) * fmt->log2_b;
2309 emax2 = fmt->emax * fmt->log2_b;
2311 np2 = SIGNIFICAND_BITS - p2;
2312 switch (r->class)
2314 underflow:
2315 get_zero (r, r->sign);
2316 case rvc_zero:
2317 if (!fmt->has_signed_zero)
2318 r->sign = 0;
2319 return;
2321 overflow:
2322 get_inf (r, r->sign);
2323 case rvc_inf:
2324 return;
2326 case rvc_nan:
2327 clear_significand_below (r, np2);
2329 /* If we've cleared the entire significand, we need one bit
2330 set for this to continue to be a NaN. */
2331 for (i = 0; i < SIGSZ; ++i)
2332 if (r->sig[i])
2333 break;
2334 if (i == SIGSZ)
2335 r->sig[SIGSZ-1] = SIG_MSB >> 2;
2336 return;
2338 case rvc_normal:
2339 break;
2341 default:
2342 abort ();
2345 /* If we're not base2, normalize the exponent to a multiple of
2346 the true base. */
2347 if (fmt->log2_b != 1)
2349 int shift = r->exp & (fmt->log2_b - 1);
2350 if (shift)
2352 shift = fmt->log2_b - shift;
2353 r->sig[0] |= sticky_rshift_significand (r, r, shift);
2354 r->exp += shift;
2358 /* Check the range of the exponent. If we're out of range,
2359 either underflow or overflow. */
2360 if (r->exp > emax2)
2361 goto overflow;
2362 else if (r->exp <= emin2m1)
2364 int diff;
2366 if (!fmt->has_denorm)
2368 /* Don't underflow completely until we've had a chance to round. */
2369 if (r->exp < emin2m1)
2370 goto underflow;
2372 else
2374 diff = emin2m1 - r->exp + 1;
2375 if (diff > p2)
2376 goto underflow;
2378 /* De-normalize the significand. */
2379 r->sig[0] |= sticky_rshift_significand (r, r, diff);
2380 r->exp += diff;
2384 /* There are P2 true significand bits, followed by one guard bit,
2385 followed by one sticky bit, followed by stuff. Fold nonzero
2386 stuff into the sticky bit. */
2388 sticky = 0;
2389 for (i = 0, w = (np2 - 1) / HOST_BITS_PER_LONG; i < w; ++i)
2390 sticky |= r->sig[i];
2391 sticky |=
2392 r->sig[w] & (((unsigned long)1 << ((np2 - 1) % HOST_BITS_PER_LONG)) - 1);
2394 guard = test_significand_bit (r, np2 - 1);
2395 lsb = test_significand_bit (r, np2);
2397 /* Round to even. */
2398 if (guard && (sticky || lsb))
2400 REAL_VALUE_TYPE u;
2401 get_zero (&u, 0);
2402 set_significand_bit (&u, np2);
2404 if (add_significands (r, r, &u))
2406 /* Overflow. Means the significand had been all ones, and
2407 is now all zeros. Need to increase the exponent, and
2408 possibly re-normalize it. */
2409 if (++r->exp > emax2)
2410 goto overflow;
2411 r->sig[SIGSZ-1] = SIG_MSB;
2413 if (fmt->log2_b != 1)
2415 int shift = r->exp & (fmt->log2_b - 1);
2416 if (shift)
2418 shift = fmt->log2_b - shift;
2419 rshift_significand (r, r, shift);
2420 r->exp += shift;
2421 if (r->exp > emax2)
2422 goto overflow;
2428 /* Catch underflow that we deferred until after rounding. */
2429 if (r->exp <= emin2m1)
2430 goto underflow;
2432 /* Clear out trailing garbage. */
2433 clear_significand_below (r, np2);
2436 /* Extend or truncate to a new mode. */
2438 void
2439 real_convert (r, mode, a)
2440 REAL_VALUE_TYPE *r;
2441 enum machine_mode mode;
2442 const REAL_VALUE_TYPE *a;
2444 const struct real_format *fmt;
2446 fmt = real_format_for_mode[mode - QFmode];
2447 if (fmt == NULL)
2448 abort ();
2450 *r = *a;
2451 round_for_format (fmt, r);
2453 /* round_for_format de-normalizes denormals. Undo just that part. */
2454 if (r->class == rvc_normal)
2455 normalize (r);
2458 /* Legacy. Likewise, except return the struct directly. */
2460 REAL_VALUE_TYPE
2461 real_value_truncate (mode, a)
2462 enum machine_mode mode;
2463 REAL_VALUE_TYPE a;
2465 REAL_VALUE_TYPE r;
2466 real_convert (&r, mode, &a);
2467 return r;
2470 /* Return true if truncating to MODE is exact. */
2472 bool
2473 exact_real_truncate (mode, a)
2474 enum machine_mode mode;
2475 const REAL_VALUE_TYPE *a;
2477 REAL_VALUE_TYPE t;
2478 real_convert (&t, mode, a);
2479 return real_identical (&t, a);
2482 /* Write R to the given target format. Place the words of the result
2483 in target word order in BUF. There are always 32 bits in each
2484 long, no matter the size of the host long.
2486 Legacy: return word 0 for implementing REAL_VALUE_TO_TARGET_SINGLE. */
2488 long
2489 real_to_target_fmt (buf, r_orig, fmt)
2490 long *buf;
2491 const REAL_VALUE_TYPE *r_orig;
2492 const struct real_format *fmt;
2494 REAL_VALUE_TYPE r;
2495 long buf1;
2497 r = *r_orig;
2498 round_for_format (fmt, &r);
2500 if (!buf)
2501 buf = &buf1;
2502 (*fmt->encode) (fmt, buf, &r);
2504 return *buf;
2507 /* Similar, but look up the format from MODE. */
2509 long
2510 real_to_target (buf, r, mode)
2511 long *buf;
2512 const REAL_VALUE_TYPE *r;
2513 enum machine_mode mode;
2515 const struct real_format *fmt;
2517 fmt = real_format_for_mode[mode - QFmode];
2518 if (fmt == NULL)
2519 abort ();
2521 return real_to_target_fmt (buf, r, fmt);
2524 /* Read R from the given target format. Read the words of the result
2525 in target word order in BUF. There are always 32 bits in each
2526 long, no matter the size of the host long. */
2528 void
2529 real_from_target_fmt (r, buf, fmt)
2530 REAL_VALUE_TYPE *r;
2531 const long *buf;
2532 const struct real_format *fmt;
2534 (*fmt->decode) (fmt, r, buf);
2537 /* Similar, but look up the format from MODE. */
2539 void
2540 real_from_target (r, buf, mode)
2541 REAL_VALUE_TYPE *r;
2542 const long *buf;
2543 enum machine_mode mode;
2545 const struct real_format *fmt;
2547 fmt = real_format_for_mode[mode - QFmode];
2548 if (fmt == NULL)
2549 abort ();
2551 (*fmt->decode) (fmt, r, buf);
2554 /* Return the number of bits in the significand for MODE. */
2555 /* ??? Legacy. Should get access to real_format directly. */
2558 significand_size (mode)
2559 enum machine_mode mode;
2561 const struct real_format *fmt;
2563 fmt = real_format_for_mode[mode - QFmode];
2564 if (fmt == NULL)
2565 return 0;
2567 return fmt->p * fmt->log2_b;
2570 /* Return a hash value for the given real value. */
2571 /* ??? The "unsigned int" return value is intended to be hashval_t,
2572 but I didn't want to pull hashtab.h into real.h. */
2574 unsigned int
2575 real_hash (r)
2576 const REAL_VALUE_TYPE *r;
2578 unsigned int h;
2579 size_t i;
2581 h = r->class | (r->sign << 2);
2582 switch (r->class)
2584 case rvc_zero:
2585 case rvc_inf:
2586 break;
2588 case rvc_normal:
2589 h |= r->exp << 3;
2590 /* FALLTHRU */
2592 case rvc_nan:
2593 if (sizeof(unsigned long) > sizeof(unsigned int))
2594 for (i = 0; i < SIGSZ; ++i)
2596 unsigned long s = r->sig[i];
2597 h ^= s ^ (s >> (HOST_BITS_PER_LONG / 2));
2599 else
2600 for (i = 0; i < SIGSZ; ++i)
2601 h ^= r->sig[i];
2602 break;
2604 default:
2605 abort ();
2608 return h;
2611 /* IEEE single-precision format. */
2613 static void encode_ieee_single PARAMS ((const struct real_format *fmt,
2614 long *, const REAL_VALUE_TYPE *));
2615 static void decode_ieee_single PARAMS ((const struct real_format *,
2616 REAL_VALUE_TYPE *, const long *));
2618 static void
2619 encode_ieee_single (fmt, buf, r)
2620 const struct real_format *fmt;
2621 long *buf;
2622 const REAL_VALUE_TYPE *r;
2624 unsigned long image, sig, exp;
2625 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2627 image = r->sign << 31;
2628 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
2630 switch (r->class)
2632 case rvc_zero:
2633 break;
2635 case rvc_inf:
2636 if (fmt->has_inf)
2637 image |= 255 << 23;
2638 else
2639 image |= 0x7fffffff;
2640 break;
2642 case rvc_nan:
2643 if (fmt->has_nans)
2645 image |= 255 << 23;
2646 image |= sig;
2647 if (!fmt->qnan_msb_set)
2648 image ^= 1 << 23 | 1 << 22;
2650 else
2651 image |= 0x7fffffff;
2652 break;
2654 case rvc_normal:
2655 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2656 whereas the intermediate representation is 0.F x 2**exp.
2657 Which means we're off by one. */
2658 if (denormal)
2659 exp = 0;
2660 else
2661 exp = r->exp + 127 - 1;
2662 image |= exp << 23;
2663 image |= sig;
2664 break;
2666 default:
2667 abort ();
2670 buf[0] = image;
2673 static void
2674 decode_ieee_single (fmt, r, buf)
2675 const struct real_format *fmt;
2676 REAL_VALUE_TYPE *r;
2677 const long *buf;
2679 unsigned long image = buf[0] & 0xffffffff;
2680 bool sign = (image >> 31) & 1;
2681 int exp = (image >> 23) & 0xff;
2683 memset (r, 0, sizeof (*r));
2684 image <<= HOST_BITS_PER_LONG - 24;
2685 image &= ~SIG_MSB;
2687 if (exp == 0)
2689 if (image && fmt->has_denorm)
2691 r->class = rvc_normal;
2692 r->sign = sign;
2693 r->exp = -126;
2694 r->sig[SIGSZ-1] = image << 1;
2695 normalize (r);
2697 else if (fmt->has_signed_zero)
2698 r->sign = sign;
2700 else if (exp == 255 && (fmt->has_nans || fmt->has_inf))
2702 if (image)
2704 r->class = rvc_nan;
2705 r->sign = sign;
2706 if (!fmt->qnan_msb_set)
2707 image ^= (SIG_MSB >> 1 | SIG_MSB >> 2);
2708 r->sig[SIGSZ-1] = image;
2710 else
2712 r->class = rvc_inf;
2713 r->sign = sign;
2716 else
2718 r->class = rvc_normal;
2719 r->sign = sign;
2720 r->exp = exp - 127 + 1;
2721 r->sig[SIGSZ-1] = image | SIG_MSB;
2725 const struct real_format ieee_single_format =
2727 encode_ieee_single,
2728 decode_ieee_single,
2732 -125,
2733 128,
2735 true,
2736 true,
2737 true,
2738 true,
2739 true
2743 /* IEEE double-precision format. */
2745 static void encode_ieee_double PARAMS ((const struct real_format *fmt,
2746 long *, const REAL_VALUE_TYPE *));
2747 static void decode_ieee_double PARAMS ((const struct real_format *,
2748 REAL_VALUE_TYPE *, const long *));
2750 static void
2751 encode_ieee_double (fmt, buf, r)
2752 const struct real_format *fmt;
2753 long *buf;
2754 const REAL_VALUE_TYPE *r;
2756 unsigned long image_lo, image_hi, sig_lo, sig_hi, exp;
2757 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2759 image_hi = r->sign << 31;
2760 image_lo = 0;
2762 if (HOST_BITS_PER_LONG == 64)
2764 sig_hi = r->sig[SIGSZ-1];
2765 sig_lo = (sig_hi >> (64 - 53)) & 0xffffffff;
2766 sig_hi = (sig_hi >> (64 - 53 + 1) >> 31) & 0xfffff;
2768 else
2770 sig_hi = r->sig[SIGSZ-1];
2771 sig_lo = r->sig[SIGSZ-2];
2772 sig_lo = (sig_hi << 21) | (sig_lo >> 11);
2773 sig_hi = (sig_hi >> 11) & 0xfffff;
2776 switch (r->class)
2778 case rvc_zero:
2779 break;
2781 case rvc_inf:
2782 if (fmt->has_inf)
2783 image_hi |= 2047 << 20;
2784 else
2786 image_hi |= 0x7fffffff;
2787 image_lo = 0xffffffff;
2789 break;
2791 case rvc_nan:
2792 if (fmt->has_nans)
2794 image_hi |= 2047 << 20;
2795 image_hi |= sig_hi;
2796 if (!fmt->qnan_msb_set)
2797 image_hi ^= 1 << 19 | 1 << 18;
2798 image_lo = sig_lo;
2800 else
2802 image_hi |= 0x7fffffff;
2803 image_lo = 0xffffffff;
2805 break;
2807 case rvc_normal:
2808 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2809 whereas the intermediate representation is 0.F x 2**exp.
2810 Which means we're off by one. */
2811 if (denormal)
2812 exp = 0;
2813 else
2814 exp = r->exp + 1023 - 1;
2815 image_hi |= exp << 20;
2816 image_hi |= sig_hi;
2817 image_lo = sig_lo;
2818 break;
2820 default:
2821 abort ();
2824 if (FLOAT_WORDS_BIG_ENDIAN)
2825 buf[0] = image_hi, buf[1] = image_lo;
2826 else
2827 buf[0] = image_lo, buf[1] = image_hi;
2830 static void
2831 decode_ieee_double (fmt, r, buf)
2832 const struct real_format *fmt;
2833 REAL_VALUE_TYPE *r;
2834 const long *buf;
2836 unsigned long image_hi, image_lo;
2837 bool sign;
2838 int exp;
2840 if (FLOAT_WORDS_BIG_ENDIAN)
2841 image_hi = buf[0], image_lo = buf[1];
2842 else
2843 image_lo = buf[0], image_hi = buf[1];
2844 image_lo &= 0xffffffff;
2845 image_hi &= 0xffffffff;
2847 sign = (image_hi >> 31) & 1;
2848 exp = (image_hi >> 20) & 0x7ff;
2850 memset (r, 0, sizeof (*r));
2852 image_hi <<= 32 - 21;
2853 image_hi |= image_lo >> 21;
2854 image_hi &= 0x7fffffff;
2855 image_lo <<= 32 - 21;
2857 if (exp == 0)
2859 if ((image_hi || image_lo) && fmt->has_denorm)
2861 r->class = rvc_normal;
2862 r->sign = sign;
2863 r->exp = -1022;
2864 if (HOST_BITS_PER_LONG == 32)
2866 image_hi = (image_hi << 1) | (image_lo >> 31);
2867 image_lo <<= 1;
2868 r->sig[SIGSZ-1] = image_hi;
2869 r->sig[SIGSZ-2] = image_lo;
2871 else
2873 image_hi = (image_hi << 31 << 2) | (image_lo << 1);
2874 r->sig[SIGSZ-1] = image_hi;
2876 normalize (r);
2878 else if (fmt->has_signed_zero)
2879 r->sign = sign;
2881 else if (exp == 2047 && (fmt->has_nans || fmt->has_inf))
2883 if (image_hi || image_lo)
2885 r->class = rvc_nan;
2886 r->sign = sign;
2887 if (HOST_BITS_PER_LONG == 32)
2889 r->sig[SIGSZ-1] = image_hi;
2890 r->sig[SIGSZ-2] = image_lo;
2892 else
2893 r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo;
2895 if (!fmt->qnan_msb_set)
2896 r->sig[SIGSZ-1] ^= (SIG_MSB >> 1 | SIG_MSB >> 2);
2898 else
2900 r->class = rvc_inf;
2901 r->sign = sign;
2904 else
2906 r->class = rvc_normal;
2907 r->sign = sign;
2908 r->exp = exp - 1023 + 1;
2909 if (HOST_BITS_PER_LONG == 32)
2911 r->sig[SIGSZ-1] = image_hi | SIG_MSB;
2912 r->sig[SIGSZ-2] = image_lo;
2914 else
2915 r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo | SIG_MSB;
2919 const struct real_format ieee_double_format =
2921 encode_ieee_double,
2922 decode_ieee_double,
2926 -1021,
2927 1024,
2929 true,
2930 true,
2931 true,
2932 true,
2933 true
2937 /* IEEE extended double precision format. This comes in three
2938 flavours: Intel's as a 12 byte image, Intel's as a 16 byte image,
2939 and Motorola's. */
2941 static void encode_ieee_extended PARAMS ((const struct real_format *fmt,
2942 long *, const REAL_VALUE_TYPE *));
2943 static void decode_ieee_extended PARAMS ((const struct real_format *,
2944 REAL_VALUE_TYPE *, const long *));
2946 static void encode_ieee_extended_128 PARAMS ((const struct real_format *fmt,
2947 long *,
2948 const REAL_VALUE_TYPE *));
2949 static void decode_ieee_extended_128 PARAMS ((const struct real_format *,
2950 REAL_VALUE_TYPE *,
2951 const long *));
2953 static void
2954 encode_ieee_extended (fmt, buf, r)
2955 const struct real_format *fmt;
2956 long *buf;
2957 const REAL_VALUE_TYPE *r;
2959 unsigned long image_hi, sig_hi, sig_lo;
2960 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2962 image_hi = r->sign << 15;
2963 sig_hi = sig_lo = 0;
2965 switch (r->class)
2967 case rvc_zero:
2968 break;
2970 case rvc_inf:
2971 if (fmt->has_inf)
2973 image_hi |= 32767;
2975 /* Intel requires the explicit integer bit to be set, otherwise
2976 it considers the value a "pseudo-infinity". Motorola docs
2977 say it doesn't care. */
2978 sig_hi = 0x80000000;
2980 else
2982 image_hi |= 32767;
2983 sig_lo = sig_hi = 0xffffffff;
2985 break;
2987 case rvc_nan:
2988 if (fmt->has_nans)
2990 image_hi |= 32767;
2991 if (HOST_BITS_PER_LONG == 32)
2993 sig_hi = r->sig[SIGSZ-1];
2994 sig_lo = r->sig[SIGSZ-2];
2996 else
2998 sig_lo = r->sig[SIGSZ-1];
2999 sig_hi = sig_lo >> 31 >> 1;
3000 sig_lo &= 0xffffffff;
3002 if (!fmt->qnan_msb_set)
3003 sig_hi ^= 1 << 30 | 1 << 29;
3005 /* Intel requires the explicit integer bit to be set, otherwise
3006 it considers the value a "pseudo-nan". Motorola docs say it
3007 doesn't care. */
3008 sig_hi |= 0x80000000;
3010 else
3012 image_hi |= 32767;
3013 sig_lo = sig_hi = 0xffffffff;
3015 break;
3017 case rvc_normal:
3019 int exp = r->exp;
3021 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3022 whereas the intermediate representation is 0.F x 2**exp.
3023 Which means we're off by one.
3025 Except for Motorola, which consider exp=0 and explicit
3026 integer bit set to continue to be normalized. In theory
3027 this discrepancy has been taken care of by the difference
3028 in fmt->emin in round_for_format. */
3030 if (denormal)
3031 exp = 0;
3032 else
3034 exp += 16383 - 1;
3035 if (exp < 0)
3036 abort ();
3038 image_hi |= exp;
3040 if (HOST_BITS_PER_LONG == 32)
3042 sig_hi = r->sig[SIGSZ-1];
3043 sig_lo = r->sig[SIGSZ-2];
3045 else
3047 sig_lo = r->sig[SIGSZ-1];
3048 sig_hi = sig_lo >> 31 >> 1;
3049 sig_lo &= 0xffffffff;
3052 break;
3054 default:
3055 abort ();
3058 if (FLOAT_WORDS_BIG_ENDIAN)
3059 buf[0] = image_hi << 16, buf[1] = sig_hi, buf[2] = sig_lo;
3060 else
3061 buf[0] = sig_lo, buf[1] = sig_hi, buf[2] = image_hi;
3064 static void
3065 encode_ieee_extended_128 (fmt, buf, r)
3066 const struct real_format *fmt;
3067 long *buf;
3068 const REAL_VALUE_TYPE *r;
3070 buf[3 * !FLOAT_WORDS_BIG_ENDIAN] = 0;
3071 encode_ieee_extended (fmt, buf+!!FLOAT_WORDS_BIG_ENDIAN, r);
3074 static void
3075 decode_ieee_extended (fmt, r, buf)
3076 const struct real_format *fmt;
3077 REAL_VALUE_TYPE *r;
3078 const long *buf;
3080 unsigned long image_hi, sig_hi, sig_lo;
3081 bool sign;
3082 int exp;
3084 if (FLOAT_WORDS_BIG_ENDIAN)
3085 image_hi = buf[0] >> 16, sig_hi = buf[1], sig_lo = buf[2];
3086 else
3087 sig_lo = buf[0], sig_hi = buf[1], image_hi = buf[2];
3088 sig_lo &= 0xffffffff;
3089 sig_hi &= 0xffffffff;
3090 image_hi &= 0xffffffff;
3092 sign = (image_hi >> 15) & 1;
3093 exp = image_hi & 0x7fff;
3095 memset (r, 0, sizeof (*r));
3097 if (exp == 0)
3099 if ((sig_hi || sig_lo) && fmt->has_denorm)
3101 r->class = rvc_normal;
3102 r->sign = sign;
3104 /* When the IEEE format contains a hidden bit, we know that
3105 it's zero at this point, and so shift up the significand
3106 and decrease the exponent to match. In this case, Motorola
3107 defines the explicit integer bit to be valid, so we don't
3108 know whether the msb is set or not. */
3109 r->exp = fmt->emin;
3110 if (HOST_BITS_PER_LONG == 32)
3112 r->sig[SIGSZ-1] = sig_hi;
3113 r->sig[SIGSZ-2] = sig_lo;
3115 else
3116 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3118 normalize (r);
3120 else if (fmt->has_signed_zero)
3121 r->sign = sign;
3123 else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3125 /* See above re "pseudo-infinities" and "pseudo-nans".
3126 Short summary is that the MSB will likely always be
3127 set, and that we don't care about it. */
3128 sig_hi &= 0x7fffffff;
3130 if (sig_hi || sig_lo)
3132 r->class = rvc_nan;
3133 r->sign = sign;
3134 if (HOST_BITS_PER_LONG == 32)
3136 r->sig[SIGSZ-1] = sig_hi;
3137 r->sig[SIGSZ-2] = sig_lo;
3139 else
3140 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3142 if (!fmt->qnan_msb_set)
3143 r->sig[SIGSZ-1] ^= (SIG_MSB >> 1 | SIG_MSB >> 2);
3145 else
3147 r->class = rvc_inf;
3148 r->sign = sign;
3151 else
3153 r->class = rvc_normal;
3154 r->sign = sign;
3155 r->exp = exp - 16383 + 1;
3156 if (HOST_BITS_PER_LONG == 32)
3158 r->sig[SIGSZ-1] = sig_hi;
3159 r->sig[SIGSZ-2] = sig_lo;
3161 else
3162 r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3166 static void
3167 decode_ieee_extended_128 (fmt, r, buf)
3168 const struct real_format *fmt;
3169 REAL_VALUE_TYPE *r;
3170 const long *buf;
3172 decode_ieee_extended (fmt, r, buf+!!FLOAT_WORDS_BIG_ENDIAN);
3175 const struct real_format ieee_extended_motorola_format =
3177 encode_ieee_extended,
3178 decode_ieee_extended,
3182 -16382,
3183 16384,
3185 true,
3186 true,
3187 true,
3188 true,
3189 true
3192 const struct real_format ieee_extended_intel_96_format =
3194 encode_ieee_extended,
3195 decode_ieee_extended,
3199 -16381,
3200 16384,
3202 true,
3203 true,
3204 true,
3205 true,
3206 true
3209 const struct real_format ieee_extended_intel_128_format =
3211 encode_ieee_extended_128,
3212 decode_ieee_extended_128,
3216 -16381,
3217 16384,
3219 true,
3220 true,
3221 true,
3222 true,
3223 true
3227 /* IBM 128-bit extended precision format: a pair of IEEE double precision
3228 numbers whose sum is equal to the extended precision value. The number
3229 with greater magnitude is first. This format has the same magnitude
3230 range as an IEEE double precision value, but effectively 106 bits of
3231 significand precision. Infinity and NaN are represented by their IEEE
3232 double precision value stored in the first number, the second number is
3233 ignored. Zeroes, Infinities, and NaNs are set in both doubles
3234 due to precedent. */
3236 static void encode_ibm_extended PARAMS ((const struct real_format *fmt,
3237 long *, const REAL_VALUE_TYPE *));
3238 static void decode_ibm_extended PARAMS ((const struct real_format *,
3239 REAL_VALUE_TYPE *, const long *));
3241 static void
3242 encode_ibm_extended (fmt, buf, r)
3243 const struct real_format *fmt ATTRIBUTE_UNUSED;
3244 long *buf;
3245 const REAL_VALUE_TYPE *r;
3247 REAL_VALUE_TYPE u, v;
3249 switch (r->class)
3251 case rvc_zero:
3252 /* Both doubles have sign bit set. */
3253 buf[0] = FLOAT_WORDS_BIG_ENDIAN ? r->sign << 31 : 0;
3254 buf[1] = FLOAT_WORDS_BIG_ENDIAN ? 0 : r->sign << 31;
3255 buf[2] = buf[0];
3256 buf[3] = buf[1];
3257 break;
3259 case rvc_inf:
3260 case rvc_nan:
3261 /* Both doubles set to Inf / NaN. */
3262 encode_ieee_double (&ieee_double_format, &buf[0], r);
3263 buf[2] = buf[0];
3264 buf[3] = buf[1];
3265 return;
3267 case rvc_normal:
3268 /* u = IEEE double precision portion of significand. */
3269 u = *r;
3270 clear_significand_below (&u, SIGNIFICAND_BITS - 53);
3272 normalize (&u);
3273 /* If the upper double is zero, we have a denormal double, so
3274 move it to the first double and leave the second as zero. */
3275 if (u.class == rvc_zero)
3277 v = u;
3278 u = *r;
3279 normalize (&u);
3281 else
3283 /* v = remainder containing additional 53 bits of significand. */
3284 do_add (&v, r, &u, 1);
3285 round_for_format (&ieee_double_format, &v);
3288 round_for_format (&ieee_double_format, &u);
3290 encode_ieee_double (&ieee_double_format, &buf[0], &u);
3291 encode_ieee_double (&ieee_double_format, &buf[2], &v);
3292 break;
3294 default:
3295 abort ();
3299 static void
3300 decode_ibm_extended (fmt, r, buf)
3301 const struct real_format *fmt ATTRIBUTE_UNUSED;
3302 REAL_VALUE_TYPE *r;
3303 const long *buf;
3305 REAL_VALUE_TYPE u, v;
3307 decode_ieee_double (&ieee_double_format, &u, &buf[0]);
3309 if (u.class != rvc_zero && u.class != rvc_inf && u.class != rvc_nan)
3311 decode_ieee_double (&ieee_double_format, &v, &buf[2]);
3312 do_add (r, &u, &v, 0);
3314 else
3315 *r = u;
3318 const struct real_format ibm_extended_format =
3320 encode_ibm_extended,
3321 decode_ibm_extended,
3324 53 + 53,
3325 -1021 + 53,
3326 1024,
3328 true,
3329 true,
3330 true,
3331 true,
3332 true
3336 /* IEEE quad precision format. */
3338 static void encode_ieee_quad PARAMS ((const struct real_format *fmt,
3339 long *, const REAL_VALUE_TYPE *));
3340 static void decode_ieee_quad PARAMS ((const struct real_format *,
3341 REAL_VALUE_TYPE *, const long *));
3343 static void
3344 encode_ieee_quad (fmt, buf, r)
3345 const struct real_format *fmt;
3346 long *buf;
3347 const REAL_VALUE_TYPE *r;
3349 unsigned long image3, image2, image1, image0, exp;
3350 bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3351 REAL_VALUE_TYPE u;
3353 image3 = r->sign << 31;
3354 image2 = 0;
3355 image1 = 0;
3356 image0 = 0;
3358 rshift_significand (&u, r, SIGNIFICAND_BITS - 113);
3360 switch (r->class)
3362 case rvc_zero:
3363 break;
3365 case rvc_inf:
3366 if (fmt->has_inf)
3367 image3 |= 32767 << 16;
3368 else
3370 image3 |= 0x7fffffff;
3371 image2 = 0xffffffff;
3372 image1 = 0xffffffff;
3373 image0 = 0xffffffff;
3375 break;
3377 case rvc_nan:
3378 if (fmt->has_nans)
3380 image3 |= 32767 << 16;
3382 if (HOST_BITS_PER_LONG == 32)
3384 image0 = u.sig[0];
3385 image1 = u.sig[1];
3386 image2 = u.sig[2];
3387 image3 |= u.sig[3] & 0xffff;
3389 else
3391 image0 = u.sig[0];
3392 image1 = image0 >> 31 >> 1;
3393 image2 = u.sig[1];
3394 image3 |= (image2 >> 31 >> 1) & 0xffff;
3395 image0 &= 0xffffffff;
3396 image2 &= 0xffffffff;
3399 if (!fmt->qnan_msb_set)
3400 image3 ^= 1 << 15 | 1 << 14;
3402 else
3404 image3 |= 0x7fffffff;
3405 image2 = 0xffffffff;
3406 image1 = 0xffffffff;
3407 image0 = 0xffffffff;
3409 break;
3411 case rvc_normal:
3412 /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3413 whereas the intermediate representation is 0.F x 2**exp.
3414 Which means we're off by one. */
3415 if (denormal)
3416 exp = 0;
3417 else
3418 exp = r->exp + 16383 - 1;
3419 image3 |= exp << 16;
3421 if (HOST_BITS_PER_LONG == 32)
3423 image0 = u.sig[0];
3424 image1 = u.sig[1];
3425 image2 = u.sig[2];
3426 image3 |= u.sig[3] & 0xffff;
3428 else
3430 image0 = u.sig[0];
3431 image1 = image0 >> 31 >> 1;
3432 image2 = u.sig[1];
3433 image3 |= (image2 >> 31 >> 1) & 0xffff;
3434 image0 &= 0xffffffff;
3435 image2 &= 0xffffffff;
3437 break;
3439 default:
3440 abort ();
3443 if (FLOAT_WORDS_BIG_ENDIAN)
3445 buf[0] = image3;
3446 buf[1] = image2;
3447 buf[2] = image1;
3448 buf[3] = image0;
3450 else
3452 buf[0] = image0;
3453 buf[1] = image1;
3454 buf[2] = image2;
3455 buf[3] = image3;
3459 static void
3460 decode_ieee_quad (fmt, r, buf)
3461 const struct real_format *fmt;
3462 REAL_VALUE_TYPE *r;
3463 const long *buf;
3465 unsigned long image3, image2, image1, image0;
3466 bool sign;
3467 int exp;
3469 if (FLOAT_WORDS_BIG_ENDIAN)
3471 image3 = buf[0];
3472 image2 = buf[1];
3473 image1 = buf[2];
3474 image0 = buf[3];
3476 else
3478 image0 = buf[0];
3479 image1 = buf[1];
3480 image2 = buf[2];
3481 image3 = buf[3];
3483 image0 &= 0xffffffff;
3484 image1 &= 0xffffffff;
3485 image2 &= 0xffffffff;
3487 sign = (image3 >> 31) & 1;
3488 exp = (image3 >> 16) & 0x7fff;
3489 image3 &= 0xffff;
3491 memset (r, 0, sizeof (*r));
3493 if (exp == 0)
3495 if ((image3 | image2 | image1 | image0) && fmt->has_denorm)
3497 r->class = rvc_normal;
3498 r->sign = sign;
3500 r->exp = -16382 + (SIGNIFICAND_BITS - 112);
3501 if (HOST_BITS_PER_LONG == 32)
3503 r->sig[0] = image0;
3504 r->sig[1] = image1;
3505 r->sig[2] = image2;
3506 r->sig[3] = image3;
3508 else
3510 r->sig[0] = (image1 << 31 << 1) | image0;
3511 r->sig[1] = (image3 << 31 << 1) | image2;
3514 normalize (r);
3516 else if (fmt->has_signed_zero)
3517 r->sign = sign;
3519 else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3521 if (image3 | image2 | image1 | image0)
3523 r->class = rvc_nan;
3524 r->sign = sign;
3526 if (HOST_BITS_PER_LONG == 32)
3528 r->sig[0] = image0;
3529 r->sig[1] = image1;
3530 r->sig[2] = image2;
3531 r->sig[3] = image3;
3533 else
3535 r->sig[0] = (image1 << 31 << 1) | image0;
3536 r->sig[1] = (image3 << 31 << 1) | image2;
3538 lshift_significand (r, r, SIGNIFICAND_BITS - 113);
3540 if (!fmt->qnan_msb_set)
3541 r->sig[SIGSZ-1] ^= (SIG_MSB >> 1 | SIG_MSB >> 2);
3543 else
3545 r->class = rvc_inf;
3546 r->sign = sign;
3549 else
3551 r->class = rvc_normal;
3552 r->sign = sign;
3553 r->exp = exp - 16383 + 1;
3555 if (HOST_BITS_PER_LONG == 32)
3557 r->sig[0] = image0;
3558 r->sig[1] = image1;
3559 r->sig[2] = image2;
3560 r->sig[3] = image3;
3562 else
3564 r->sig[0] = (image1 << 31 << 1) | image0;
3565 r->sig[1] = (image3 << 31 << 1) | image2;
3567 lshift_significand (r, r, SIGNIFICAND_BITS - 113);
3568 r->sig[SIGSZ-1] |= SIG_MSB;
3572 const struct real_format ieee_quad_format =
3574 encode_ieee_quad,
3575 decode_ieee_quad,
3578 113,
3579 -16381,
3580 16384,
3581 127,
3582 true,
3583 true,
3584 true,
3585 true,
3586 true
3589 /* Descriptions of VAX floating point formats can be found beginning at
3591 http://www.openvms.compaq.com:8000/73final/4515/4515pro_013.html#f_floating_point_format
3593 The thing to remember is that they're almost IEEE, except for word
3594 order, exponent bias, and the lack of infinities, nans, and denormals.
3596 We don't implement the H_floating format here, simply because neither
3597 the VAX or Alpha ports use it. */
3599 static void encode_vax_f PARAMS ((const struct real_format *fmt,
3600 long *, const REAL_VALUE_TYPE *));
3601 static void decode_vax_f PARAMS ((const struct real_format *,
3602 REAL_VALUE_TYPE *, const long *));
3603 static void encode_vax_d PARAMS ((const struct real_format *fmt,
3604 long *, const REAL_VALUE_TYPE *));
3605 static void decode_vax_d PARAMS ((const struct real_format *,
3606 REAL_VALUE_TYPE *, const long *));
3607 static void encode_vax_g PARAMS ((const struct real_format *fmt,
3608 long *, const REAL_VALUE_TYPE *));
3609 static void decode_vax_g PARAMS ((const struct real_format *,
3610 REAL_VALUE_TYPE *, const long *));
3612 static void
3613 encode_vax_f (fmt, buf, r)
3614 const struct real_format *fmt ATTRIBUTE_UNUSED;
3615 long *buf;
3616 const REAL_VALUE_TYPE *r;
3618 unsigned long sign, exp, sig, image;
3620 sign = r->sign << 15;
3622 switch (r->class)
3624 case rvc_zero:
3625 image = 0;
3626 break;
3628 case rvc_inf:
3629 case rvc_nan:
3630 image = 0xffff7fff | sign;
3631 break;
3633 case rvc_normal:
3634 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
3635 exp = r->exp + 128;
3637 image = (sig << 16) & 0xffff0000;
3638 image |= sign;
3639 image |= exp << 7;
3640 image |= sig >> 16;
3641 break;
3643 default:
3644 abort ();
3647 buf[0] = image;
3650 static void
3651 decode_vax_f (fmt, r, buf)
3652 const struct real_format *fmt ATTRIBUTE_UNUSED;
3653 REAL_VALUE_TYPE *r;
3654 const long *buf;
3656 unsigned long image = buf[0] & 0xffffffff;
3657 int exp = (image >> 7) & 0xff;
3659 memset (r, 0, sizeof (*r));
3661 if (exp != 0)
3663 r->class = rvc_normal;
3664 r->sign = (image >> 15) & 1;
3665 r->exp = exp - 128;
3667 image = ((image & 0x7f) << 16) | ((image >> 16) & 0xffff);
3668 r->sig[SIGSZ-1] = (image << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
3672 static void
3673 encode_vax_d (fmt, buf, r)
3674 const struct real_format *fmt ATTRIBUTE_UNUSED;
3675 long *buf;
3676 const REAL_VALUE_TYPE *r;
3678 unsigned long image0, image1, sign = r->sign << 15;
3680 switch (r->class)
3682 case rvc_zero:
3683 image0 = image1 = 0;
3684 break;
3686 case rvc_inf:
3687 case rvc_nan:
3688 image0 = 0xffff7fff | sign;
3689 image1 = 0xffffffff;
3690 break;
3692 case rvc_normal:
3693 /* Extract the significand into straight hi:lo. */
3694 if (HOST_BITS_PER_LONG == 64)
3696 image0 = r->sig[SIGSZ-1];
3697 image1 = (image0 >> (64 - 56)) & 0xffffffff;
3698 image0 = (image0 >> (64 - 56 + 1) >> 31) & 0x7fffff;
3700 else
3702 image0 = r->sig[SIGSZ-1];
3703 image1 = r->sig[SIGSZ-2];
3704 image1 = (image0 << 24) | (image1 >> 8);
3705 image0 = (image0 >> 8) & 0xffffff;
3708 /* Rearrange the half-words of the significand to match the
3709 external format. */
3710 image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff007f;
3711 image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
3713 /* Add the sign and exponent. */
3714 image0 |= sign;
3715 image0 |= (r->exp + 128) << 7;
3716 break;
3718 default:
3719 abort ();
3722 if (FLOAT_WORDS_BIG_ENDIAN)
3723 buf[0] = image1, buf[1] = image0;
3724 else
3725 buf[0] = image0, buf[1] = image1;
3728 static void
3729 decode_vax_d (fmt, r, buf)
3730 const struct real_format *fmt ATTRIBUTE_UNUSED;
3731 REAL_VALUE_TYPE *r;
3732 const long *buf;
3734 unsigned long image0, image1;
3735 int exp;
3737 if (FLOAT_WORDS_BIG_ENDIAN)
3738 image1 = buf[0], image0 = buf[1];
3739 else
3740 image0 = buf[0], image1 = buf[1];
3741 image0 &= 0xffffffff;
3742 image1 &= 0xffffffff;
3744 exp = (image0 >> 7) & 0x7f;
3746 memset (r, 0, sizeof (*r));
3748 if (exp != 0)
3750 r->class = rvc_normal;
3751 r->sign = (image0 >> 15) & 1;
3752 r->exp = exp - 128;
3754 /* Rearrange the half-words of the external format into
3755 proper ascending order. */
3756 image0 = ((image0 & 0x7f) << 16) | ((image0 >> 16) & 0xffff);
3757 image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
3759 if (HOST_BITS_PER_LONG == 64)
3761 image0 = (image0 << 31 << 1) | image1;
3762 image0 <<= 64 - 56;
3763 image0 |= SIG_MSB;
3764 r->sig[SIGSZ-1] = image0;
3766 else
3768 r->sig[SIGSZ-1] = image0;
3769 r->sig[SIGSZ-2] = image1;
3770 lshift_significand (r, r, 2*HOST_BITS_PER_LONG - 56);
3771 r->sig[SIGSZ-1] |= SIG_MSB;
3776 static void
3777 encode_vax_g (fmt, buf, r)
3778 const struct real_format *fmt ATTRIBUTE_UNUSED;
3779 long *buf;
3780 const REAL_VALUE_TYPE *r;
3782 unsigned long image0, image1, sign = r->sign << 15;
3784 switch (r->class)
3786 case rvc_zero:
3787 image0 = image1 = 0;
3788 break;
3790 case rvc_inf:
3791 case rvc_nan:
3792 image0 = 0xffff7fff | sign;
3793 image1 = 0xffffffff;
3794 break;
3796 case rvc_normal:
3797 /* Extract the significand into straight hi:lo. */
3798 if (HOST_BITS_PER_LONG == 64)
3800 image0 = r->sig[SIGSZ-1];
3801 image1 = (image0 >> (64 - 53)) & 0xffffffff;
3802 image0 = (image0 >> (64 - 53 + 1) >> 31) & 0xfffff;
3804 else
3806 image0 = r->sig[SIGSZ-1];
3807 image1 = r->sig[SIGSZ-2];
3808 image1 = (image0 << 21) | (image1 >> 11);
3809 image0 = (image0 >> 11) & 0xfffff;
3812 /* Rearrange the half-words of the significand to match the
3813 external format. */
3814 image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff000f;
3815 image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
3817 /* Add the sign and exponent. */
3818 image0 |= sign;
3819 image0 |= (r->exp + 1024) << 4;
3820 break;
3822 default:
3823 abort ();
3826 if (FLOAT_WORDS_BIG_ENDIAN)
3827 buf[0] = image1, buf[1] = image0;
3828 else
3829 buf[0] = image0, buf[1] = image1;
3832 static void
3833 decode_vax_g (fmt, r, buf)
3834 const struct real_format *fmt ATTRIBUTE_UNUSED;
3835 REAL_VALUE_TYPE *r;
3836 const long *buf;
3838 unsigned long image0, image1;
3839 int exp;
3841 if (FLOAT_WORDS_BIG_ENDIAN)
3842 image1 = buf[0], image0 = buf[1];
3843 else
3844 image0 = buf[0], image1 = buf[1];
3845 image0 &= 0xffffffff;
3846 image1 &= 0xffffffff;
3848 exp = (image0 >> 4) & 0x7ff;
3850 memset (r, 0, sizeof (*r));
3852 if (exp != 0)
3854 r->class = rvc_normal;
3855 r->sign = (image0 >> 15) & 1;
3856 r->exp = exp - 1024;
3858 /* Rearrange the half-words of the external format into
3859 proper ascending order. */
3860 image0 = ((image0 & 0xf) << 16) | ((image0 >> 16) & 0xffff);
3861 image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
3863 if (HOST_BITS_PER_LONG == 64)
3865 image0 = (image0 << 31 << 1) | image1;
3866 image0 <<= 64 - 53;
3867 image0 |= SIG_MSB;
3868 r->sig[SIGSZ-1] = image0;
3870 else
3872 r->sig[SIGSZ-1] = image0;
3873 r->sig[SIGSZ-2] = image1;
3874 lshift_significand (r, r, 64 - 53);
3875 r->sig[SIGSZ-1] |= SIG_MSB;
3880 const struct real_format vax_f_format =
3882 encode_vax_f,
3883 decode_vax_f,
3887 -127,
3888 127,
3890 false,
3891 false,
3892 false,
3893 false,
3894 false
3897 const struct real_format vax_d_format =
3899 encode_vax_d,
3900 decode_vax_d,
3904 -127,
3905 127,
3907 false,
3908 false,
3909 false,
3910 false,
3911 false
3914 const struct real_format vax_g_format =
3916 encode_vax_g,
3917 decode_vax_g,
3921 -1023,
3922 1023,
3924 false,
3925 false,
3926 false,
3927 false,
3928 false
3931 /* A good reference for these can be found in chapter 9 of
3932 "ESA/390 Principles of Operation", IBM document number SA22-7201-01.
3933 An on-line version can be found here:
3935 http://publibz.boulder.ibm.com/cgi-bin/bookmgr_OS390/BOOKS/DZ9AR001/9.1?DT=19930923083613
3938 static void encode_i370_single PARAMS ((const struct real_format *fmt,
3939 long *, const REAL_VALUE_TYPE *));
3940 static void decode_i370_single PARAMS ((const struct real_format *,
3941 REAL_VALUE_TYPE *, const long *));
3942 static void encode_i370_double PARAMS ((const struct real_format *fmt,
3943 long *, const REAL_VALUE_TYPE *));
3944 static void decode_i370_double PARAMS ((const struct real_format *,
3945 REAL_VALUE_TYPE *, const long *));
3947 static void
3948 encode_i370_single (fmt, buf, r)
3949 const struct real_format *fmt ATTRIBUTE_UNUSED;
3950 long *buf;
3951 const REAL_VALUE_TYPE *r;
3953 unsigned long sign, exp, sig, image;
3955 sign = r->sign << 31;
3957 switch (r->class)
3959 case rvc_zero:
3960 image = 0;
3961 break;
3963 case rvc_inf:
3964 case rvc_nan:
3965 image = 0x7fffffff | sign;
3966 break;
3968 case rvc_normal:
3969 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0xffffff;
3970 exp = ((r->exp / 4) + 64) << 24;
3971 image = sign | exp | sig;
3972 break;
3974 default:
3975 abort ();
3978 buf[0] = image;
3981 static void
3982 decode_i370_single (fmt, r, buf)
3983 const struct real_format *fmt ATTRIBUTE_UNUSED;
3984 REAL_VALUE_TYPE *r;
3985 const long *buf;
3987 unsigned long sign, sig, image = buf[0];
3988 int exp;
3990 sign = (image >> 31) & 1;
3991 exp = (image >> 24) & 0x7f;
3992 sig = image & 0xffffff;
3994 memset (r, 0, sizeof (*r));
3996 if (exp || sig)
3998 r->class = rvc_normal;
3999 r->sign = sign;
4000 r->exp = (exp - 64) * 4;
4001 r->sig[SIGSZ-1] = sig << (HOST_BITS_PER_LONG - 24);
4002 normalize (r);
4006 static void
4007 encode_i370_double (fmt, buf, r)
4008 const struct real_format *fmt ATTRIBUTE_UNUSED;
4009 long *buf;
4010 const REAL_VALUE_TYPE *r;
4012 unsigned long sign, exp, image_hi, image_lo;
4014 sign = r->sign << 31;
4016 switch (r->class)
4018 case rvc_zero:
4019 image_hi = image_lo = 0;
4020 break;
4022 case rvc_inf:
4023 case rvc_nan:
4024 image_hi = 0x7fffffff | sign;
4025 image_lo = 0xffffffff;
4026 break;
4028 case rvc_normal:
4029 if (HOST_BITS_PER_LONG == 64)
4031 image_hi = r->sig[SIGSZ-1];
4032 image_lo = (image_hi >> (64 - 56)) & 0xffffffff;
4033 image_hi = (image_hi >> (64 - 56 + 1) >> 31) & 0xffffff;
4035 else
4037 image_hi = r->sig[SIGSZ-1];
4038 image_lo = r->sig[SIGSZ-2];
4039 image_lo = (image_lo >> 8) | (image_hi << 24);
4040 image_hi >>= 8;
4043 exp = ((r->exp / 4) + 64) << 24;
4044 image_hi |= sign | exp;
4045 break;
4047 default:
4048 abort ();
4051 if (FLOAT_WORDS_BIG_ENDIAN)
4052 buf[0] = image_hi, buf[1] = image_lo;
4053 else
4054 buf[0] = image_lo, buf[1] = image_hi;
4057 static void
4058 decode_i370_double (fmt, r, buf)
4059 const struct real_format *fmt ATTRIBUTE_UNUSED;
4060 REAL_VALUE_TYPE *r;
4061 const long *buf;
4063 unsigned long sign, image_hi, image_lo;
4064 int exp;
4066 if (FLOAT_WORDS_BIG_ENDIAN)
4067 image_hi = buf[0], image_lo = buf[1];
4068 else
4069 image_lo = buf[0], image_hi = buf[1];
4071 sign = (image_hi >> 31) & 1;
4072 exp = (image_hi >> 24) & 0x7f;
4073 image_hi &= 0xffffff;
4074 image_lo &= 0xffffffff;
4076 memset (r, 0, sizeof (*r));
4078 if (exp || image_hi || image_lo)
4080 r->class = rvc_normal;
4081 r->sign = sign;
4082 r->exp = (exp - 64) * 4 + (SIGNIFICAND_BITS - 56);
4084 if (HOST_BITS_PER_LONG == 32)
4086 r->sig[0] = image_lo;
4087 r->sig[1] = image_hi;
4089 else
4090 r->sig[0] = image_lo | (image_hi << 31 << 1);
4092 normalize (r);
4096 const struct real_format i370_single_format =
4098 encode_i370_single,
4099 decode_i370_single,
4103 -64,
4106 false,
4107 false,
4108 false, /* ??? The encoding does allow for "unnormals". */
4109 false, /* ??? The encoding does allow for "unnormals". */
4110 false
4113 const struct real_format i370_double_format =
4115 encode_i370_double,
4116 decode_i370_double,
4120 -64,
4123 false,
4124 false,
4125 false, /* ??? The encoding does allow for "unnormals". */
4126 false, /* ??? The encoding does allow for "unnormals". */
4127 false
4130 /* The "twos-complement" c4x format is officially defined as
4132 x = s(~s).f * 2**e
4134 This is rather misleading. One must remember that F is signed.
4135 A better description would be
4137 x = -1**s * ((s + 1 + .f) * 2**e
4139 So if we have a (4 bit) fraction of .1000 with a sign bit of 1,
4140 that's -1 * (1+1+(-.5)) == -1.5. I think.
4142 The constructions here are taken from Tables 5-1 and 5-2 of the
4143 TMS320C4x User's Guide wherein step-by-step instructions for
4144 conversion from IEEE are presented. That's close enough to our
4145 internal representation so as to make things easy.
4147 See http://www-s.ti.com/sc/psheets/spru063c/spru063c.pdf */
4149 static void encode_c4x_single PARAMS ((const struct real_format *fmt,
4150 long *, const REAL_VALUE_TYPE *));
4151 static void decode_c4x_single PARAMS ((const struct real_format *,
4152 REAL_VALUE_TYPE *, const long *));
4153 static void encode_c4x_extended PARAMS ((const struct real_format *fmt,
4154 long *, const REAL_VALUE_TYPE *));
4155 static void decode_c4x_extended PARAMS ((const struct real_format *,
4156 REAL_VALUE_TYPE *, const long *));
4158 static void
4159 encode_c4x_single (fmt, buf, r)
4160 const struct real_format *fmt ATTRIBUTE_UNUSED;
4161 long *buf;
4162 const REAL_VALUE_TYPE *r;
4164 unsigned long image, exp, sig;
4166 switch (r->class)
4168 case rvc_zero:
4169 exp = -128;
4170 sig = 0;
4171 break;
4173 case rvc_inf:
4174 case rvc_nan:
4175 exp = 127;
4176 sig = 0x800000 - r->sign;
4177 break;
4179 case rvc_normal:
4180 exp = r->exp - 1;
4181 sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
4182 if (r->sign)
4184 if (sig)
4185 sig = -sig;
4186 else
4187 exp--;
4188 sig |= 0x800000;
4190 break;
4192 default:
4193 abort ();
4196 image = ((exp & 0xff) << 24) | (sig & 0xffffff);
4197 buf[0] = image;
4200 static void
4201 decode_c4x_single (fmt, r, buf)
4202 const struct real_format *fmt ATTRIBUTE_UNUSED;
4203 REAL_VALUE_TYPE *r;
4204 const long *buf;
4206 unsigned long image = buf[0];
4207 unsigned long sig;
4208 int exp, sf;
4210 exp = (((image >> 24) & 0xff) ^ 0x80) - 0x80;
4211 sf = ((image & 0xffffff) ^ 0x800000) - 0x800000;
4213 memset (r, 0, sizeof (*r));
4215 if (exp != -128)
4217 r->class = rvc_normal;
4219 sig = sf & 0x7fffff;
4220 if (sf < 0)
4222 r->sign = 1;
4223 if (sig)
4224 sig = -sig;
4225 else
4226 exp++;
4228 sig = (sig << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
4230 r->exp = exp + 1;
4231 r->sig[SIGSZ-1] = sig;
4235 static void
4236 encode_c4x_extended (fmt, buf, r)
4237 const struct real_format *fmt ATTRIBUTE_UNUSED;
4238 long *buf;
4239 const REAL_VALUE_TYPE *r;
4241 unsigned long exp, sig;
4243 switch (r->class)
4245 case rvc_zero:
4246 exp = -128;
4247 sig = 0;
4248 break;
4250 case rvc_inf:
4251 case rvc_nan:
4252 exp = 127;
4253 sig = 0x80000000 - r->sign;
4254 break;
4256 case rvc_normal:
4257 exp = r->exp - 1;
4259 sig = r->sig[SIGSZ-1];
4260 if (HOST_BITS_PER_LONG == 64)
4261 sig = sig >> 1 >> 31;
4262 sig &= 0x7fffffff;
4264 if (r->sign)
4266 if (sig)
4267 sig = -sig;
4268 else
4269 exp--;
4270 sig |= 0x80000000;
4272 break;
4274 default:
4275 abort ();
4278 exp = (exp & 0xff) << 24;
4279 sig &= 0xffffffff;
4281 if (FLOAT_WORDS_BIG_ENDIAN)
4282 buf[0] = exp, buf[1] = sig;
4283 else
4284 buf[0] = sig, buf[0] = exp;
4287 static void
4288 decode_c4x_extended (fmt, r, buf)
4289 const struct real_format *fmt ATTRIBUTE_UNUSED;
4290 REAL_VALUE_TYPE *r;
4291 const long *buf;
4293 unsigned long sig;
4294 int exp, sf;
4296 if (FLOAT_WORDS_BIG_ENDIAN)
4297 exp = buf[0], sf = buf[1];
4298 else
4299 sf = buf[0], exp = buf[1];
4301 exp = (((exp >> 24) & 0xff) & 0x80) - 0x80;
4302 sf = ((sf & 0xffffffff) ^ 0x80000000) - 0x80000000;
4304 memset (r, 0, sizeof (*r));
4306 if (exp != -128)
4308 r->class = rvc_normal;
4310 sig = sf & 0x7fffffff;
4311 if (sf < 0)
4313 r->sign = 1;
4314 if (sig)
4315 sig = -sig;
4316 else
4317 exp++;
4319 if (HOST_BITS_PER_LONG == 64)
4320 sig = sig << 1 << 31;
4321 sig |= SIG_MSB;
4323 r->exp = exp + 1;
4324 r->sig[SIGSZ-1] = sig;
4328 const struct real_format c4x_single_format =
4330 encode_c4x_single,
4331 decode_c4x_single,
4335 -126,
4336 128,
4338 false,
4339 false,
4340 false,
4341 false,
4342 false
4345 const struct real_format c4x_extended_format =
4347 encode_c4x_extended,
4348 decode_c4x_extended,
4352 -126,
4353 128,
4355 false,
4356 false,
4357 false,
4358 false,
4359 false
4363 /* A synthetic "format" for internal arithmetic. It's the size of the
4364 internal significand minus the two bits needed for proper rounding.
4365 The encode and decode routines exist only to satisfy our paranoia
4366 harness. */
4368 static void encode_internal PARAMS ((const struct real_format *fmt,
4369 long *, const REAL_VALUE_TYPE *));
4370 static void decode_internal PARAMS ((const struct real_format *,
4371 REAL_VALUE_TYPE *, const long *));
4373 static void
4374 encode_internal (fmt, buf, r)
4375 const struct real_format *fmt ATTRIBUTE_UNUSED;
4376 long *buf;
4377 const REAL_VALUE_TYPE *r;
4379 memcpy (buf, r, sizeof (*r));
4382 static void
4383 decode_internal (fmt, r, buf)
4384 const struct real_format *fmt ATTRIBUTE_UNUSED;
4385 REAL_VALUE_TYPE *r;
4386 const long *buf;
4388 memcpy (r, buf, sizeof (*r));
4391 const struct real_format real_internal_format =
4393 encode_internal,
4394 decode_internal,
4397 SIGNIFICAND_BITS - 2,
4398 -MAX_EXP,
4399 MAX_EXP,
4401 true,
4402 true,
4403 false,
4404 true,
4405 true
4408 /* Set up default mode to format mapping for IEEE. Everyone else has
4409 to set these values in OVERRIDE_OPTIONS. */
4411 const struct real_format *real_format_for_mode[TFmode - QFmode + 1] =
4413 NULL, /* QFmode */
4414 NULL, /* HFmode */
4415 NULL, /* TQFmode */
4416 &ieee_single_format, /* SFmode */
4417 &ieee_double_format, /* DFmode */
4419 /* We explicitly don't handle XFmode. There are two formats,
4420 pretty much equally common. Choose one in OVERRIDE_OPTIONS. */
4421 NULL, /* XFmode */
4422 &ieee_quad_format /* TFmode */
4426 /* Calculate the square root of X in mode MODE, and store the result
4427 in R. Return TRUE if the operation does not raise an exception.
4428 For details see "High Precision Division and Square Root",
4429 Alan H. Karp and Peter Markstein, HP Lab Report 93-93-42, June
4430 1993. http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf. */
4432 bool
4433 real_sqrt (r, mode, x)
4434 REAL_VALUE_TYPE *r;
4435 enum machine_mode mode;
4436 const REAL_VALUE_TYPE *x;
4438 static REAL_VALUE_TYPE halfthree;
4439 static REAL_VALUE_TYPE half;
4440 static bool init = false;
4441 REAL_VALUE_TYPE h, t, i;
4442 int iter, exp;
4444 /* sqrt(-0.0) is -0.0. */
4445 if (real_isnegzero (x))
4447 *r = *x;
4448 return false;
4451 /* Negative arguments return NaN. */
4452 if (real_isneg (x))
4454 /* Mode is ignored for canonical NaN. */
4455 real_nan (r, "", 1, SFmode);
4456 return false;
4459 /* Infinity and NaN return themselves. */
4460 if (real_isinf (x) || real_isnan (x))
4462 *r = *x;
4463 return false;
4466 if (!init)
4468 real_arithmetic (&half, RDIV_EXPR, &dconst1, &dconst2);
4469 real_arithmetic (&halfthree, PLUS_EXPR, &dconst1, &half);
4470 init = true;
4473 /* Initial guess for reciprocal sqrt, i. */
4474 exp = real_exponent (x);
4475 real_ldexp (&i, &dconst1, -exp/2);
4477 /* Newton's iteration for reciprocal sqrt, i. */
4478 for (iter = 0; iter < 16; iter++)
4480 /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x). */
4481 real_arithmetic (&t, MULT_EXPR, x, &i);
4482 real_arithmetic (&h, MULT_EXPR, &t, &i);
4483 real_arithmetic (&t, MULT_EXPR, &h, &half);
4484 real_arithmetic (&h, MINUS_EXPR, &halfthree, &t);
4485 real_arithmetic (&t, MULT_EXPR, &i, &h);
4487 /* Check for early convergence. */
4488 if (iter >= 6 && real_identical (&i, &t))
4489 break;
4491 /* ??? Unroll loop to avoid copying. */
4492 i = t;
4495 /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)). */
4496 real_arithmetic (&t, MULT_EXPR, x, &i);
4497 real_arithmetic (&h, MULT_EXPR, &t, &i);
4498 real_arithmetic (&i, MINUS_EXPR, &dconst1, &h);
4499 real_arithmetic (&h, MULT_EXPR, &t, &i);
4500 real_arithmetic (&i, MULT_EXPR, &half, &h);
4501 real_arithmetic (&h, PLUS_EXPR, &t, &i);
4503 /* ??? We need a Tuckerman test to get the last bit. */
4505 real_convert (r, mode, &h);
4506 return true;