OSDN Git Service

* real.c (encode_ieee_single): Ensure proper promotion.
[pf3gnuchains/gcc-fork.git] / gcc / real.c
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>
6
7    This file is part of GCC.
8
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.
13
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.
18
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.  */
23
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"
32
33 /* The floating point model used internally is not exactly IEEE 754
34    compliant, and close to the description in the ISO C99 standard,
35    section 5.2.4.2.2 Characteristics of floating types.
36
37    Specifically
38
39         x = s * b^e * \sum_{k=1}^p f_k * b^{-k}
40
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.
47
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).
51
52    A requirement of the model is that P be larger than the largest
53    supported target floating-point type by at least 2 bits.  This gives
54    us proper rounding when we truncate to the target type.  In addition,
55    E must be large enough to hold the smallest supported denormal number
56    in a normalized form.
57
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.
61
62    Note that the decimal string conversion routines are sensitive to
63    rounding errors.  Since the raw arithmetic routines do not themselves
64    have guard digits or rounding, the computation of 10**exp can
65    accumulate more than a few digits of error.  The previous incarnation
66    of real.c successfully used a 144-bit fraction; given the current
67    layout of REAL_VALUE_TYPE we're forced to expand to at least 160 bits.
68
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.  */
73
74
75 /* Used to classify two numbers simultaneously.  */
76 #define CLASS2(A, B)  ((A) << 2 | (B))
77
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
81
82 static void get_zero (REAL_VALUE_TYPE *, int);
83 static void get_canonical_qnan (REAL_VALUE_TYPE *, int);
84 static void get_canonical_snan (REAL_VALUE_TYPE *, int);
85 static void get_inf (REAL_VALUE_TYPE *, int);
86 static bool sticky_rshift_significand (REAL_VALUE_TYPE *,
87                                        const REAL_VALUE_TYPE *, unsigned int);
88 static void rshift_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
89                                 unsigned int);
90 static void lshift_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
91                                 unsigned int);
92 static void lshift_significand_1 (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
93 static bool add_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *,
94                               const REAL_VALUE_TYPE *);
95 static bool sub_significands (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
96                               const REAL_VALUE_TYPE *, int);
97 static void neg_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
98 static int cmp_significands (const REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
99 static int cmp_significand_0 (const REAL_VALUE_TYPE *);
100 static void set_significand_bit (REAL_VALUE_TYPE *, unsigned int);
101 static void clear_significand_bit (REAL_VALUE_TYPE *, unsigned int);
102 static bool test_significand_bit (REAL_VALUE_TYPE *, unsigned int);
103 static void clear_significand_below (REAL_VALUE_TYPE *, unsigned int);
104 static bool div_significands (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
105                               const REAL_VALUE_TYPE *);
106 static void normalize (REAL_VALUE_TYPE *);
107
108 static bool do_add (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
109                     const REAL_VALUE_TYPE *, int);
110 static bool do_multiply (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
111                          const REAL_VALUE_TYPE *);
112 static bool do_divide (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
113                        const REAL_VALUE_TYPE *);
114 static int do_compare (const REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *, int);
115 static void do_fix_trunc (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
116
117 static unsigned long rtd_divmod (REAL_VALUE_TYPE *, REAL_VALUE_TYPE *);
118
119 static const REAL_VALUE_TYPE * ten_to_ptwo (int);
120 static const REAL_VALUE_TYPE * ten_to_mptwo (int);
121 static const REAL_VALUE_TYPE * real_digit (int);
122 static void times_pten (REAL_VALUE_TYPE *, int);
123
124 static void round_for_format (const struct real_format *, REAL_VALUE_TYPE *);
125 \f
126 /* Initialize R with a positive zero.  */
127
128 static inline void
129 get_zero (REAL_VALUE_TYPE *r, int sign)
130 {
131   memset (r, 0, sizeof (*r));
132   r->sign = sign;
133 }
134
135 /* Initialize R with the canonical quiet NaN.  */
136
137 static inline void
138 get_canonical_qnan (REAL_VALUE_TYPE *r, int sign)
139 {
140   memset (r, 0, sizeof (*r));
141   r->class = rvc_nan;
142   r->sign = sign;
143   r->canonical = 1;
144 }
145
146 static inline void
147 get_canonical_snan (REAL_VALUE_TYPE *r, int sign)
148 {
149   memset (r, 0, sizeof (*r));
150   r->class = rvc_nan;
151   r->sign = sign;
152   r->signalling = 1;
153   r->canonical = 1;
154 }
155
156 static inline void
157 get_inf (REAL_VALUE_TYPE *r, int sign)
158 {
159   memset (r, 0, sizeof (*r));
160   r->class = rvc_inf;
161   r->sign = sign;
162 }
163
164 \f
165 /* Right-shift the significand of A by N bits; put the result in the
166    significand of R.  If any one bits are shifted out, return true.  */
167
168 static bool
169 sticky_rshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
170                            unsigned int n)
171 {
172   unsigned long sticky = 0;
173   unsigned int i, ofs = 0;
174
175   if (n >= HOST_BITS_PER_LONG)
176     {
177       for (i = 0, ofs = n / HOST_BITS_PER_LONG; i < ofs; ++i)
178         sticky |= a->sig[i];
179       n &= HOST_BITS_PER_LONG - 1;
180     }
181
182   if (n != 0)
183     {
184       sticky |= a->sig[ofs] & (((unsigned long)1 << n) - 1);
185       for (i = 0; i < SIGSZ; ++i)
186         {
187           r->sig[i]
188             = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
189                | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
190                   << (HOST_BITS_PER_LONG - n)));
191         }
192     }
193   else
194     {
195       for (i = 0; ofs + i < SIGSZ; ++i)
196         r->sig[i] = a->sig[ofs + i];
197       for (; i < SIGSZ; ++i)
198         r->sig[i] = 0;
199     }
200
201   return sticky != 0;
202 }
203
204 /* Right-shift the significand of A by N bits; put the result in the
205    significand of R.  */
206
207 static void
208 rshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
209                     unsigned int n)
210 {
211   unsigned int i, ofs = n / HOST_BITS_PER_LONG;
212
213   n &= HOST_BITS_PER_LONG - 1;
214   if (n != 0)
215     {
216       for (i = 0; i < SIGSZ; ++i)
217         {
218           r->sig[i]
219             = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
220                | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
221                   << (HOST_BITS_PER_LONG - n)));
222         }
223     }
224   else
225     {
226       for (i = 0; ofs + i < SIGSZ; ++i)
227         r->sig[i] = a->sig[ofs + i];
228       for (; i < SIGSZ; ++i)
229         r->sig[i] = 0;
230     }
231 }
232
233 /* Left-shift the significand of A by N bits; put the result in the
234    significand of R.  */
235
236 static void
237 lshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
238                     unsigned int n)
239 {
240   unsigned int i, ofs = n / HOST_BITS_PER_LONG;
241
242   n &= HOST_BITS_PER_LONG - 1;
243   if (n == 0)
244     {
245       for (i = 0; ofs + i < SIGSZ; ++i)
246         r->sig[SIGSZ-1-i] = a->sig[SIGSZ-1-i-ofs];
247       for (; i < SIGSZ; ++i)
248         r->sig[SIGSZ-1-i] = 0;
249     }
250   else
251     for (i = 0; i < SIGSZ; ++i)
252       {
253         r->sig[SIGSZ-1-i]
254           = (((ofs + i >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs]) << n)
255              | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs-1])
256                 >> (HOST_BITS_PER_LONG - n)));
257       }
258 }
259
260 /* Likewise, but N is specialized to 1.  */
261
262 static inline void
263 lshift_significand_1 (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
264 {
265   unsigned int i;
266
267   for (i = SIGSZ - 1; i > 0; --i)
268     r->sig[i] = (a->sig[i] << 1) | (a->sig[i-1] >> (HOST_BITS_PER_LONG - 1));
269   r->sig[0] = a->sig[0] << 1;
270 }
271
272 /* Add the significands of A and B, placing the result in R.  Return
273    true if there was carry out of the most significant word.  */
274
275 static inline bool
276 add_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
277                   const REAL_VALUE_TYPE *b)
278 {
279   bool carry = false;
280   int i;
281
282   for (i = 0; i < SIGSZ; ++i)
283     {
284       unsigned long ai = a->sig[i];
285       unsigned long ri = ai + b->sig[i];
286
287       if (carry)
288         {
289           carry = ri < ai;
290           carry |= ++ri == 0;
291         }
292       else
293         carry = ri < ai;
294
295       r->sig[i] = ri;
296     }
297
298   return carry;
299 }
300
301 /* Subtract the significands of A and B, placing the result in R.  CARRY is
302    true if there's a borrow incoming to the least significant word.
303    Return true if there was borrow out of the most significant word.  */
304
305 static inline bool
306 sub_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
307                   const REAL_VALUE_TYPE *b, int carry)
308 {
309   int i;
310
311   for (i = 0; i < SIGSZ; ++i)
312     {
313       unsigned long ai = a->sig[i];
314       unsigned long ri = ai - b->sig[i];
315
316       if (carry)
317         {
318           carry = ri > ai;
319           carry |= ~--ri == 0;
320         }
321       else
322         carry = ri > ai;
323
324       r->sig[i] = ri;
325     }
326
327   return carry;
328 }
329
330 /* Negate the significand A, placing the result in R.  */
331
332 static inline void
333 neg_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
334 {
335   bool carry = true;
336   int i;
337
338   for (i = 0; i < SIGSZ; ++i)
339     {
340       unsigned long ri, ai = a->sig[i];
341
342       if (carry)
343         {
344           if (ai)
345             {
346               ri = -ai;
347               carry = false;
348             }
349           else
350             ri = ai;
351         }
352       else
353         ri = ~ai;
354
355       r->sig[i] = ri;
356     }
357 }
358
359 /* Compare significands.  Return tri-state vs zero.  */
360
361 static inline int
362 cmp_significands (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
363 {
364   int i;
365
366   for (i = SIGSZ - 1; i >= 0; --i)
367     {
368       unsigned long ai = a->sig[i];
369       unsigned long bi = b->sig[i];
370
371       if (ai > bi)
372         return 1;
373       if (ai < bi)
374         return -1;
375     }
376
377   return 0;
378 }
379
380 /* Return true if A is nonzero.  */
381
382 static inline int
383 cmp_significand_0 (const REAL_VALUE_TYPE *a)
384 {
385   int i;
386
387   for (i = SIGSZ - 1; i >= 0; --i)
388     if (a->sig[i])
389       return 1;
390
391   return 0;
392 }
393
394 /* Set bit N of the significand of R.  */
395
396 static inline void
397 set_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
398 {
399   r->sig[n / HOST_BITS_PER_LONG]
400     |= (unsigned long)1 << (n % HOST_BITS_PER_LONG);
401 }
402
403 /* Clear bit N of the significand of R.  */
404
405 static inline void
406 clear_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
407 {
408   r->sig[n / HOST_BITS_PER_LONG]
409     &= ~((unsigned long)1 << (n % HOST_BITS_PER_LONG));
410 }
411
412 /* Test bit N of the significand of R.  */
413
414 static inline bool
415 test_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
416 {
417   /* ??? Compiler bug here if we return this expression directly.
418      The conversion to bool strips the "&1" and we wind up testing
419      e.g. 2 != 0 -> true.  Seen in gcc version 3.2 20020520.  */
420   int t = (r->sig[n / HOST_BITS_PER_LONG] >> (n % HOST_BITS_PER_LONG)) & 1;
421   return t;
422 }
423
424 /* Clear bits 0..N-1 of the significand of R.  */
425
426 static void
427 clear_significand_below (REAL_VALUE_TYPE *r, unsigned int n)
428 {
429   int i, w = n / HOST_BITS_PER_LONG;
430
431   for (i = 0; i < w; ++i)
432     r->sig[i] = 0;
433
434   r->sig[w] &= ~(((unsigned long)1 << (n % HOST_BITS_PER_LONG)) - 1);
435 }
436
437 /* Divide the significands of A and B, placing the result in R.  Return
438    true if the division was inexact.  */
439
440 static inline bool
441 div_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
442                   const REAL_VALUE_TYPE *b)
443 {
444   REAL_VALUE_TYPE u;
445   int i, bit = SIGNIFICAND_BITS - 1;
446   unsigned long msb, inexact;
447
448   u = *a;
449   memset (r->sig, 0, sizeof (r->sig));
450
451   msb = 0;
452   goto start;
453   do
454     {
455       msb = u.sig[SIGSZ-1] & SIG_MSB;
456       lshift_significand_1 (&u, &u);
457     start:
458       if (msb || cmp_significands (&u, b) >= 0)
459         {
460           sub_significands (&u, &u, b, 0);
461           set_significand_bit (r, bit);
462         }
463     }
464   while (--bit >= 0);
465
466   for (i = 0, inexact = 0; i < SIGSZ; i++)
467     inexact |= u.sig[i];
468
469   return inexact != 0;
470 }
471
472 /* Adjust the exponent and significand of R such that the most
473    significant bit is set.  We underflow to zero and overflow to
474    infinity here, without denormals.  (The intermediate representation
475    exponent is large enough to handle target denormals normalized.)  */
476
477 static void
478 normalize (REAL_VALUE_TYPE *r)
479 {
480   int shift = 0, exp;
481   int i, j;
482
483   /* Find the first word that is nonzero.  */
484   for (i = SIGSZ - 1; i >= 0; i--)
485     if (r->sig[i] == 0)
486       shift += HOST_BITS_PER_LONG;
487     else
488       break;
489
490   /* Zero significand flushes to zero.  */
491   if (i < 0)
492     {
493       r->class = rvc_zero;
494       r->exp = 0;
495       return;
496     }
497
498   /* Find the first bit that is nonzero.  */
499   for (j = 0; ; j++)
500     if (r->sig[i] & ((unsigned long)1 << (HOST_BITS_PER_LONG - 1 - j)))
501       break;
502   shift += j;
503
504   if (shift > 0)
505     {
506       exp = r->exp - shift;
507       if (exp > MAX_EXP)
508         get_inf (r, r->sign);
509       else if (exp < -MAX_EXP)
510         get_zero (r, r->sign);
511       else
512         {
513           r->exp = exp;
514           lshift_significand (r, r, shift);
515         }
516     }
517 }
518 \f
519 /* Calculate R = A + (SUBTRACT_P ? -B : B).  Return true if the
520    result may be inexact due to a loss of precision.  */
521
522 static bool
523 do_add (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
524         const REAL_VALUE_TYPE *b, int subtract_p)
525 {
526   int dexp, sign, exp;
527   REAL_VALUE_TYPE t;
528   bool inexact = false;
529
530   /* Determine if we need to add or subtract.  */
531   sign = a->sign;
532   subtract_p = (sign ^ b->sign) ^ subtract_p;
533
534   switch (CLASS2 (a->class, b->class))
535     {
536     case CLASS2 (rvc_zero, rvc_zero):
537       /* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0.  */
538       get_zero (r, sign & !subtract_p);
539       return false;
540
541     case CLASS2 (rvc_zero, rvc_normal):
542     case CLASS2 (rvc_zero, rvc_inf):
543     case CLASS2 (rvc_zero, rvc_nan):
544       /* 0 + ANY = ANY.  */
545     case CLASS2 (rvc_normal, rvc_nan):
546     case CLASS2 (rvc_inf, rvc_nan):
547     case CLASS2 (rvc_nan, rvc_nan):
548       /* ANY + NaN = NaN.  */
549     case CLASS2 (rvc_normal, rvc_inf):
550       /* R + Inf = Inf.  */
551       *r = *b;
552       r->sign = sign ^ subtract_p;
553       return false;
554
555     case CLASS2 (rvc_normal, rvc_zero):
556     case CLASS2 (rvc_inf, rvc_zero):
557     case CLASS2 (rvc_nan, rvc_zero):
558       /* ANY + 0 = ANY.  */
559     case CLASS2 (rvc_nan, rvc_normal):
560     case CLASS2 (rvc_nan, rvc_inf):
561       /* NaN + ANY = NaN.  */
562     case CLASS2 (rvc_inf, rvc_normal):
563       /* Inf + R = Inf.  */
564       *r = *a;
565       return false;
566
567     case CLASS2 (rvc_inf, rvc_inf):
568       if (subtract_p)
569         /* Inf - Inf = NaN.  */
570         get_canonical_qnan (r, 0);
571       else
572         /* Inf + Inf = Inf.  */
573         *r = *a;
574       return false;
575
576     case CLASS2 (rvc_normal, rvc_normal):
577       break;
578
579     default:
580       abort ();
581     }
582
583   /* Swap the arguments such that A has the larger exponent.  */
584   dexp = a->exp - b->exp;
585   if (dexp < 0)
586     {
587       const REAL_VALUE_TYPE *t;
588       t = a, a = b, b = t;
589       dexp = -dexp;
590       sign ^= subtract_p;
591     }
592   exp = a->exp;
593
594   /* If the exponents are not identical, we need to shift the
595      significand of B down.  */
596   if (dexp > 0)
597     {
598       /* If the exponents are too far apart, the significands
599          do not overlap, which makes the subtraction a noop.  */
600       if (dexp >= SIGNIFICAND_BITS)
601         {
602           *r = *a;
603           r->sign = sign;
604           return true;
605         }
606
607       inexact |= sticky_rshift_significand (&t, b, dexp);
608       b = &t;
609     }
610
611   if (subtract_p)
612     {
613       if (sub_significands (r, a, b, inexact))
614         {
615           /* We got a borrow out of the subtraction.  That means that
616              A and B had the same exponent, and B had the larger
617              significand.  We need to swap the sign and negate the
618              significand.  */
619           sign ^= 1;
620           neg_significand (r, r);
621         }
622     }
623   else
624     {
625       if (add_significands (r, a, b))
626         {
627           /* We got carry out of the addition.  This means we need to
628              shift the significand back down one bit and increase the
629              exponent.  */
630           inexact |= sticky_rshift_significand (r, r, 1);
631           r->sig[SIGSZ-1] |= SIG_MSB;
632           if (++exp > MAX_EXP)
633             {
634               get_inf (r, sign);
635               return true;
636             }
637         }
638     }
639
640   r->class = rvc_normal;
641   r->sign = sign;
642   r->exp = exp;
643
644   /* Re-normalize the result.  */
645   normalize (r);
646
647   /* Special case: if the subtraction results in zero, the result
648      is positive.  */
649   if (r->class == rvc_zero)
650     r->sign = 0;
651   else
652     r->sig[0] |= inexact;
653
654   return inexact;
655 }
656
657 /* Calculate R = A * B.  Return true if the result may be inexact.  */
658
659 static bool
660 do_multiply (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
661              const REAL_VALUE_TYPE *b)
662 {
663   REAL_VALUE_TYPE u, t, *rr;
664   unsigned int i, j, k;
665   int sign = a->sign ^ b->sign;
666   bool inexact = false;
667
668   switch (CLASS2 (a->class, b->class))
669     {
670     case CLASS2 (rvc_zero, rvc_zero):
671     case CLASS2 (rvc_zero, rvc_normal):
672     case CLASS2 (rvc_normal, rvc_zero):
673       /* +-0 * ANY = 0 with appropriate sign.  */
674       get_zero (r, sign);
675       return false;
676
677     case CLASS2 (rvc_zero, rvc_nan):
678     case CLASS2 (rvc_normal, rvc_nan):
679     case CLASS2 (rvc_inf, rvc_nan):
680     case CLASS2 (rvc_nan, rvc_nan):
681       /* ANY * NaN = NaN.  */
682       *r = *b;
683       r->sign = sign;
684       return false;
685
686     case CLASS2 (rvc_nan, rvc_zero):
687     case CLASS2 (rvc_nan, rvc_normal):
688     case CLASS2 (rvc_nan, rvc_inf):
689       /* NaN * ANY = NaN.  */
690       *r = *a;
691       r->sign = sign;
692       return false;
693
694     case CLASS2 (rvc_zero, rvc_inf):
695     case CLASS2 (rvc_inf, rvc_zero):
696       /* 0 * Inf = NaN */
697       get_canonical_qnan (r, sign);
698       return false;
699
700     case CLASS2 (rvc_inf, rvc_inf):
701     case CLASS2 (rvc_normal, rvc_inf):
702     case CLASS2 (rvc_inf, rvc_normal):
703       /* Inf * Inf = Inf, R * Inf = Inf */
704       get_inf (r, sign);
705       return false;
706
707     case CLASS2 (rvc_normal, rvc_normal):
708       break;
709
710     default:
711       abort ();
712     }
713
714   if (r == a || r == b)
715     rr = &t;
716   else
717     rr = r;
718   get_zero (rr, 0);
719
720   /* Collect all the partial products.  Since we don't have sure access
721      to a widening multiply, we split each long into two half-words.
722
723      Consider the long-hand form of a four half-word multiplication:
724
725                  A  B  C  D
726               *  E  F  G  H
727              --------------
728                 DE DF DG DH
729              CE CF CG CH
730           BE BF BG BH
731        AE AF AG AH
732
733      We construct partial products of the widened half-word products
734      that are known to not overlap, e.g. DF+DH.  Each such partial
735      product is given its proper exponent, which allows us to sum them
736      and obtain the finished product.  */
737
738   for (i = 0; i < SIGSZ * 2; ++i)
739     {
740       unsigned long ai = a->sig[i / 2];
741       if (i & 1)
742         ai >>= HOST_BITS_PER_LONG / 2;
743       else
744         ai &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
745
746       if (ai == 0)
747         continue;
748
749       for (j = 0; j < 2; ++j)
750         {
751           int exp = (a->exp - (2*SIGSZ-1-i)*(HOST_BITS_PER_LONG/2)
752                      + (b->exp - (1-j)*(HOST_BITS_PER_LONG/2)));
753
754           if (exp > MAX_EXP)
755             {
756               get_inf (r, sign);
757               return true;
758             }
759           if (exp < -MAX_EXP)
760             {
761               /* Would underflow to zero, which we shouldn't bother adding.  */
762               inexact = true;
763               continue;
764             }
765
766           memset (&u, 0, sizeof (u));
767           u.class = rvc_normal;
768           u.exp = exp;
769
770           for (k = j; k < SIGSZ * 2; k += 2)
771             {
772               unsigned long bi = b->sig[k / 2];
773               if (k & 1)
774                 bi >>= HOST_BITS_PER_LONG / 2;
775               else
776                 bi &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
777
778               u.sig[k / 2] = ai * bi;
779             }
780
781           normalize (&u);
782           inexact |= do_add (rr, rr, &u, 0);
783         }
784     }
785
786   rr->sign = sign;
787   if (rr != r)
788     *r = t;
789
790   return inexact;
791 }
792
793 /* Calculate R = A / B.  Return true if the result may be inexact.  */
794
795 static bool
796 do_divide (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
797            const REAL_VALUE_TYPE *b)
798 {
799   int exp, sign = a->sign ^ b->sign;
800   REAL_VALUE_TYPE t, *rr;
801   bool inexact;
802
803   switch (CLASS2 (a->class, b->class))
804     {
805     case CLASS2 (rvc_zero, rvc_zero):
806       /* 0 / 0 = NaN.  */
807     case CLASS2 (rvc_inf, rvc_inf):
808       /* Inf / Inf = NaN.  */
809       get_canonical_qnan (r, sign);
810       return false;
811
812     case CLASS2 (rvc_zero, rvc_normal):
813     case CLASS2 (rvc_zero, rvc_inf):
814       /* 0 / ANY = 0.  */
815     case CLASS2 (rvc_normal, rvc_inf):
816       /* R / Inf = 0.  */
817       get_zero (r, sign);
818       return false;
819
820     case CLASS2 (rvc_normal, rvc_zero):
821       /* R / 0 = Inf.  */
822     case CLASS2 (rvc_inf, rvc_zero):
823       /* Inf / 0 = Inf.  */
824       get_inf (r, sign);
825       return false;
826
827     case CLASS2 (rvc_zero, rvc_nan):
828     case CLASS2 (rvc_normal, rvc_nan):
829     case CLASS2 (rvc_inf, rvc_nan):
830     case CLASS2 (rvc_nan, rvc_nan):
831       /* ANY / NaN = NaN.  */
832       *r = *b;
833       r->sign = sign;
834       return false;
835
836     case CLASS2 (rvc_nan, rvc_zero):
837     case CLASS2 (rvc_nan, rvc_normal):
838     case CLASS2 (rvc_nan, rvc_inf):
839       /* NaN / ANY = NaN.  */
840       *r = *a;
841       r->sign = sign;
842       return false;
843
844     case CLASS2 (rvc_inf, rvc_normal):
845       /* Inf / R = Inf.  */
846       get_inf (r, sign);
847       return false;
848
849     case CLASS2 (rvc_normal, rvc_normal):
850       break;
851
852     default:
853       abort ();
854     }
855
856   if (r == a || r == b)
857     rr = &t;
858   else
859     rr = r;
860
861   /* Make sure all fields in the result are initialized.  */
862   get_zero (rr, 0);
863   rr->class = rvc_normal;
864   rr->sign = sign;
865
866   exp = a->exp - b->exp + 1;
867   if (exp > MAX_EXP)
868     {
869       get_inf (r, sign);
870       return true;
871     }
872   if (exp < -MAX_EXP)
873     {
874       get_zero (r, sign);
875       return true;
876     }
877   rr->exp = exp;
878
879   inexact = div_significands (rr, a, b);
880
881   /* Re-normalize the result.  */
882   normalize (rr);
883   rr->sig[0] |= inexact;
884
885   if (rr != r)
886     *r = t;
887
888   return inexact;
889 }
890
891 /* Return a tri-state comparison of A vs B.  Return NAN_RESULT if
892    one of the two operands is a NaN.  */
893
894 static int
895 do_compare (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b,
896             int nan_result)
897 {
898   int ret;
899
900   switch (CLASS2 (a->class, b->class))
901     {
902     case CLASS2 (rvc_zero, rvc_zero):
903       /* Sign of zero doesn't matter for compares.  */
904       return 0;
905
906     case CLASS2 (rvc_inf, rvc_zero):
907     case CLASS2 (rvc_inf, rvc_normal):
908     case CLASS2 (rvc_normal, rvc_zero):
909       return (a->sign ? -1 : 1);
910
911     case CLASS2 (rvc_inf, rvc_inf):
912       return -a->sign - -b->sign;
913
914     case CLASS2 (rvc_zero, rvc_normal):
915     case CLASS2 (rvc_zero, rvc_inf):
916     case CLASS2 (rvc_normal, rvc_inf):
917       return (b->sign ? 1 : -1);
918
919     case CLASS2 (rvc_zero, rvc_nan):
920     case CLASS2 (rvc_normal, rvc_nan):
921     case CLASS2 (rvc_inf, rvc_nan):
922     case CLASS2 (rvc_nan, rvc_nan):
923     case CLASS2 (rvc_nan, rvc_zero):
924     case CLASS2 (rvc_nan, rvc_normal):
925     case CLASS2 (rvc_nan, rvc_inf):
926       return nan_result;
927
928     case CLASS2 (rvc_normal, rvc_normal):
929       break;
930
931     default:
932       abort ();
933     }
934
935   if (a->sign != b->sign)
936     return -a->sign - -b->sign;
937
938   if (a->exp > b->exp)
939     ret = 1;
940   else if (a->exp < b->exp)
941     ret = -1;
942   else
943     ret = cmp_significands (a, b);
944
945   return (a->sign ? -ret : ret);
946 }
947
948 /* Return A truncated to an integral value toward zero.  */
949
950 static void
951 do_fix_trunc (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
952 {
953   *r = *a;
954
955   switch (r->class)
956     {
957     case rvc_zero:
958     case rvc_inf:
959     case rvc_nan:
960       break;
961
962     case rvc_normal:
963       if (r->exp <= 0)
964         get_zero (r, r->sign);
965       else if (r->exp < SIGNIFICAND_BITS)
966         clear_significand_below (r, SIGNIFICAND_BITS - r->exp);
967       break;
968
969     default:
970       abort ();
971     }
972 }
973
974 /* Perform the binary or unary operation described by CODE.
975    For a unary operation, leave OP1 NULL.  */
976
977 void
978 real_arithmetic (REAL_VALUE_TYPE *r, int icode, const REAL_VALUE_TYPE *op0,
979                  const REAL_VALUE_TYPE *op1)
980 {
981   enum tree_code code = icode;
982
983   switch (code)
984     {
985     case PLUS_EXPR:
986       do_add (r, op0, op1, 0);
987       break;
988
989     case MINUS_EXPR:
990       do_add (r, op0, op1, 1);
991       break;
992
993     case MULT_EXPR:
994       do_multiply (r, op0, op1);
995       break;
996
997     case RDIV_EXPR:
998       do_divide (r, op0, op1);
999       break;
1000
1001     case MIN_EXPR:
1002       if (op1->class == rvc_nan)
1003         *r = *op1;
1004       else if (do_compare (op0, op1, -1) < 0)
1005         *r = *op0;
1006       else
1007         *r = *op1;
1008       break;
1009
1010     case MAX_EXPR:
1011       if (op1->class == rvc_nan)
1012         *r = *op1;
1013       else if (do_compare (op0, op1, 1) < 0)
1014         *r = *op1;
1015       else
1016         *r = *op0;
1017       break;
1018
1019     case NEGATE_EXPR:
1020       *r = *op0;
1021       r->sign ^= 1;
1022       break;
1023
1024     case ABS_EXPR:
1025       *r = *op0;
1026       r->sign = 0;
1027       break;
1028
1029     case FIX_TRUNC_EXPR:
1030       do_fix_trunc (r, op0);
1031       break;
1032
1033     default:
1034       abort ();
1035     }
1036 }
1037
1038 /* Legacy.  Similar, but return the result directly.  */
1039
1040 REAL_VALUE_TYPE
1041 real_arithmetic2 (int icode, const REAL_VALUE_TYPE *op0,
1042                   const REAL_VALUE_TYPE *op1)
1043 {
1044   REAL_VALUE_TYPE r;
1045   real_arithmetic (&r, icode, op0, op1);
1046   return r;
1047 }
1048
1049 bool
1050 real_compare (int icode, const REAL_VALUE_TYPE *op0,
1051               const REAL_VALUE_TYPE *op1)
1052 {
1053   enum tree_code code = icode;
1054
1055   switch (code)
1056     {
1057     case LT_EXPR:
1058       return do_compare (op0, op1, 1) < 0;
1059     case LE_EXPR:
1060       return do_compare (op0, op1, 1) <= 0;
1061     case GT_EXPR:
1062       return do_compare (op0, op1, -1) > 0;
1063     case GE_EXPR:
1064       return do_compare (op0, op1, -1) >= 0;
1065     case EQ_EXPR:
1066       return do_compare (op0, op1, -1) == 0;
1067     case NE_EXPR:
1068       return do_compare (op0, op1, -1) != 0;
1069     case UNORDERED_EXPR:
1070       return op0->class == rvc_nan || op1->class == rvc_nan;
1071     case ORDERED_EXPR:
1072       return op0->class != rvc_nan && op1->class != rvc_nan;
1073     case UNLT_EXPR:
1074       return do_compare (op0, op1, -1) < 0;
1075     case UNLE_EXPR:
1076       return do_compare (op0, op1, -1) <= 0;
1077     case UNGT_EXPR:
1078       return do_compare (op0, op1, 1) > 0;
1079     case UNGE_EXPR:
1080       return do_compare (op0, op1, 1) >= 0;
1081     case UNEQ_EXPR:
1082       return do_compare (op0, op1, 0) == 0;
1083
1084     default:
1085       abort ();
1086     }
1087 }
1088
1089 /* Return floor log2(R).  */
1090
1091 int
1092 real_exponent (const REAL_VALUE_TYPE *r)
1093 {
1094   switch (r->class)
1095     {
1096     case rvc_zero:
1097       return 0;
1098     case rvc_inf:
1099     case rvc_nan:
1100       return (unsigned int)-1 >> 1;
1101     case rvc_normal:
1102       return r->exp;
1103     default:
1104       abort ();
1105     }
1106 }
1107
1108 /* R = OP0 * 2**EXP.  */
1109
1110 void
1111 real_ldexp (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *op0, int exp)
1112 {
1113   *r = *op0;
1114   switch (r->class)
1115     {
1116     case rvc_zero:
1117     case rvc_inf:
1118     case rvc_nan:
1119       break;
1120
1121     case rvc_normal:
1122       exp += op0->exp;
1123       if (exp > MAX_EXP)
1124         get_inf (r, r->sign);
1125       else if (exp < -MAX_EXP)
1126         get_zero (r, r->sign);
1127       else
1128         r->exp = exp;
1129       break;
1130
1131     default:
1132       abort ();
1133     }
1134 }
1135
1136 /* Determine whether a floating-point value X is infinite.  */
1137
1138 bool
1139 real_isinf (const REAL_VALUE_TYPE *r)
1140 {
1141   return (r->class == rvc_inf);
1142 }
1143
1144 /* Determine whether a floating-point value X is a NaN.  */
1145
1146 bool
1147 real_isnan (const REAL_VALUE_TYPE *r)
1148 {
1149   return (r->class == rvc_nan);
1150 }
1151
1152 /* Determine whether a floating-point value X is negative.  */
1153
1154 bool
1155 real_isneg (const REAL_VALUE_TYPE *r)
1156 {
1157   return r->sign;
1158 }
1159
1160 /* Determine whether a floating-point value X is minus zero.  */
1161
1162 bool
1163 real_isnegzero (const REAL_VALUE_TYPE *r)
1164 {
1165   return r->sign && r->class == rvc_zero;
1166 }
1167
1168 /* Compare two floating-point objects for bitwise identity.  */
1169
1170 bool
1171 real_identical (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
1172 {
1173   int i;
1174
1175   if (a->class != b->class)
1176     return false;
1177   if (a->sign != b->sign)
1178     return false;
1179
1180   switch (a->class)
1181     {
1182     case rvc_zero:
1183     case rvc_inf:
1184       return true;
1185
1186     case rvc_normal:
1187       if (a->exp != b->exp)
1188         return false;
1189       break;
1190
1191     case rvc_nan:
1192       if (a->signalling != b->signalling)
1193         return false;
1194       /* The significand is ignored for canonical NaNs.  */
1195       if (a->canonical || b->canonical)
1196         return a->canonical == b->canonical;
1197       break;
1198
1199     default:
1200       abort ();
1201     }
1202
1203   for (i = 0; i < SIGSZ; ++i)
1204     if (a->sig[i] != b->sig[i])
1205       return false;
1206
1207   return true;
1208 }
1209
1210 /* Try to change R into its exact multiplicative inverse in machine
1211    mode MODE.  Return true if successful.  */
1212
1213 bool
1214 exact_real_inverse (enum machine_mode mode, REAL_VALUE_TYPE *r)
1215 {
1216   const REAL_VALUE_TYPE *one = real_digit (1);
1217   REAL_VALUE_TYPE u;
1218   int i;
1219
1220   if (r->class != rvc_normal)
1221     return false;
1222
1223   /* Check for a power of two: all significand bits zero except the MSB.  */
1224   for (i = 0; i < SIGSZ-1; ++i)
1225     if (r->sig[i] != 0)
1226       return false;
1227   if (r->sig[SIGSZ-1] != SIG_MSB)
1228     return false;
1229
1230   /* Find the inverse and truncate to the required mode.  */
1231   do_divide (&u, one, r);
1232   real_convert (&u, mode, &u);
1233
1234   /* The rounding may have overflowed.  */
1235   if (u.class != rvc_normal)
1236     return false;
1237   for (i = 0; i < SIGSZ-1; ++i)
1238     if (u.sig[i] != 0)
1239       return false;
1240   if (u.sig[SIGSZ-1] != SIG_MSB)
1241     return false;
1242
1243   *r = u;
1244   return true;
1245 }
1246 \f
1247 /* Render R as an integer.  */
1248
1249 HOST_WIDE_INT
1250 real_to_integer (const REAL_VALUE_TYPE *r)
1251 {
1252   unsigned HOST_WIDE_INT i;
1253
1254   switch (r->class)
1255     {
1256     case rvc_zero:
1257     underflow:
1258       return 0;
1259
1260     case rvc_inf:
1261     case rvc_nan:
1262     overflow:
1263       i = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1264       if (!r->sign)
1265         i--;
1266       return i;
1267
1268     case rvc_normal:
1269       if (r->exp <= 0)
1270         goto underflow;
1271       /* Only force overflow for unsigned overflow.  Signed overflow is
1272          undefined, so it doesn't matter what we return, and some callers
1273          expect to be able to use this routine for both signed and
1274          unsigned conversions.  */
1275       if (r->exp > HOST_BITS_PER_WIDE_INT)
1276         goto overflow;
1277
1278       if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1279         i = r->sig[SIGSZ-1];
1280       else if (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG)
1281         {
1282           i = r->sig[SIGSZ-1];
1283           i = i << (HOST_BITS_PER_LONG - 1) << 1;
1284           i |= r->sig[SIGSZ-2];
1285         }
1286       else
1287         abort ();
1288
1289       i >>= HOST_BITS_PER_WIDE_INT - r->exp;
1290
1291       if (r->sign)
1292         i = -i;
1293       return i;
1294
1295     default:
1296       abort ();
1297     }
1298 }
1299
1300 /* Likewise, but to an integer pair, HI+LOW.  */
1301
1302 void
1303 real_to_integer2 (HOST_WIDE_INT *plow, HOST_WIDE_INT *phigh,
1304                   const REAL_VALUE_TYPE *r)
1305 {
1306   REAL_VALUE_TYPE t;
1307   HOST_WIDE_INT low, high;
1308   int exp;
1309
1310   switch (r->class)
1311     {
1312     case rvc_zero:
1313     underflow:
1314       low = high = 0;
1315       break;
1316
1317     case rvc_inf:
1318     case rvc_nan:
1319     overflow:
1320       high = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1321       if (r->sign)
1322         low = 0;
1323       else
1324         {
1325           high--;
1326           low = -1;
1327         }
1328       break;
1329
1330     case rvc_normal:
1331       exp = r->exp;
1332       if (exp <= 0)
1333         goto underflow;
1334       /* Only force overflow for unsigned overflow.  Signed overflow is
1335          undefined, so it doesn't matter what we return, and some callers
1336          expect to be able to use this routine for both signed and
1337          unsigned conversions.  */
1338       if (exp > 2*HOST_BITS_PER_WIDE_INT)
1339         goto overflow;
1340
1341       rshift_significand (&t, r, 2*HOST_BITS_PER_WIDE_INT - exp);
1342       if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1343         {
1344           high = t.sig[SIGSZ-1];
1345           low = t.sig[SIGSZ-2];
1346         }
1347       else if (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG)
1348         {
1349           high = t.sig[SIGSZ-1];
1350           high = high << (HOST_BITS_PER_LONG - 1) << 1;
1351           high |= t.sig[SIGSZ-2];
1352
1353           low = t.sig[SIGSZ-3];
1354           low = low << (HOST_BITS_PER_LONG - 1) << 1;
1355           low |= t.sig[SIGSZ-4];
1356         }
1357       else
1358         abort ();
1359
1360       if (r->sign)
1361         {
1362           if (low == 0)
1363             high = -high;
1364           else
1365             low = -low, high = ~high;
1366         }
1367       break;
1368
1369     default:
1370       abort ();
1371     }
1372
1373   *plow = low;
1374   *phigh = high;
1375 }
1376
1377 /* A subroutine of real_to_decimal.  Compute the quotient and remainder
1378    of NUM / DEN.  Return the quotient and place the remainder in NUM.
1379    It is expected that NUM / DEN are close enough that the quotient is
1380    small.  */
1381
1382 static unsigned long
1383 rtd_divmod (REAL_VALUE_TYPE *num, REAL_VALUE_TYPE *den)
1384 {
1385   unsigned long q, msb;
1386   int expn = num->exp, expd = den->exp;
1387
1388   if (expn < expd)
1389     return 0;
1390
1391   q = msb = 0;
1392   goto start;
1393   do
1394     {
1395       msb = num->sig[SIGSZ-1] & SIG_MSB;
1396       q <<= 1;
1397       lshift_significand_1 (num, num);
1398     start:
1399       if (msb || cmp_significands (num, den) >= 0)
1400         {
1401           sub_significands (num, num, den, 0);
1402           q |= 1;
1403         }
1404     }
1405   while (--expn >= expd);
1406
1407   num->exp = expd;
1408   normalize (num);
1409
1410   return q;
1411 }
1412
1413 /* Render R as a decimal floating point constant.  Emit DIGITS significant
1414    digits in the result, bounded by BUF_SIZE.  If DIGITS is 0, choose the
1415    maximum for the representation.  If CROP_TRAILING_ZEROS, strip trailing
1416    zeros.  */
1417
1418 #define M_LOG10_2       0.30102999566398119521
1419
1420 void
1421 real_to_decimal (char *str, const REAL_VALUE_TYPE *r_orig, size_t buf_size,
1422                  size_t digits, int crop_trailing_zeros)
1423 {
1424   const REAL_VALUE_TYPE *one, *ten;
1425   REAL_VALUE_TYPE r, pten, u, v;
1426   int dec_exp, cmp_one, digit;
1427   size_t max_digits;
1428   char *p, *first, *last;
1429   bool sign;
1430
1431   r = *r_orig;
1432   switch (r.class)
1433     {
1434     case rvc_zero:
1435       strcpy (str, (r.sign ? "-0.0" : "0.0"));
1436       return;
1437     case rvc_normal:
1438       break;
1439     case rvc_inf:
1440       strcpy (str, (r.sign ? "-Inf" : "+Inf"));
1441       return;
1442     case rvc_nan:
1443       /* ??? Print the significand as well, if not canonical?  */
1444       strcpy (str, (r.sign ? "-NaN" : "+NaN"));
1445       return;
1446     default:
1447       abort ();
1448     }
1449
1450   /* Bound the number of digits printed by the size of the representation.  */
1451   max_digits = SIGNIFICAND_BITS * M_LOG10_2;
1452   if (digits == 0 || digits > max_digits)
1453     digits = max_digits;
1454
1455   /* Estimate the decimal exponent, and compute the length of the string it
1456      will print as.  Be conservative and add one to account for possible
1457      overflow or rounding error.  */
1458   dec_exp = r.exp * M_LOG10_2;
1459   for (max_digits = 1; dec_exp ; max_digits++)
1460     dec_exp /= 10;
1461
1462   /* Bound the number of digits printed by the size of the output buffer.  */
1463   max_digits = buf_size - 1 - 1 - 2 - max_digits - 1;
1464   if (max_digits > buf_size)
1465     abort ();
1466   if (digits > max_digits)
1467     digits = max_digits;
1468
1469   one = real_digit (1);
1470   ten = ten_to_ptwo (0);
1471
1472   sign = r.sign;
1473   r.sign = 0;
1474
1475   dec_exp = 0;
1476   pten = *one;
1477
1478   cmp_one = do_compare (&r, one, 0);
1479   if (cmp_one > 0)
1480     {
1481       int m;
1482
1483       /* Number is greater than one.  Convert significand to an integer
1484          and strip trailing decimal zeros.  */
1485
1486       u = r;
1487       u.exp = SIGNIFICAND_BITS - 1;
1488
1489       /* Largest M, such that 10**2**M fits within SIGNIFICAND_BITS.  */
1490       m = floor_log2 (max_digits);
1491
1492       /* Iterate over the bits of the possible powers of 10 that might
1493          be present in U and eliminate them.  That is, if we find that
1494          10**2**M divides U evenly, keep the division and increase
1495          DEC_EXP by 2**M.  */
1496       do
1497         {
1498           REAL_VALUE_TYPE t;
1499
1500           do_divide (&t, &u, ten_to_ptwo (m));
1501           do_fix_trunc (&v, &t);
1502           if (cmp_significands (&v, &t) == 0)
1503             {
1504               u = t;
1505               dec_exp += 1 << m;
1506             }
1507         }
1508       while (--m >= 0);
1509
1510       /* Revert the scaling to integer that we performed earlier.  */
1511       u.exp += r.exp - (SIGNIFICAND_BITS - 1);
1512       r = u;
1513
1514       /* Find power of 10.  Do this by dividing out 10**2**M when
1515          this is larger than the current remainder.  Fill PTEN with
1516          the power of 10 that we compute.  */
1517       if (r.exp > 0)
1518         {
1519           m = floor_log2 ((int)(r.exp * M_LOG10_2)) + 1;
1520           do
1521             {
1522               const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1523               if (do_compare (&u, ptentwo, 0) >= 0)
1524                 {
1525                   do_divide (&u, &u, ptentwo);
1526                   do_multiply (&pten, &pten, ptentwo);
1527                   dec_exp += 1 << m;
1528                 }
1529             }
1530           while (--m >= 0);
1531         }
1532       else
1533         /* We managed to divide off enough tens in the above reduction
1534            loop that we've now got a negative exponent.  Fall into the
1535            less-than-one code to compute the proper value for PTEN.  */
1536         cmp_one = -1;
1537     }
1538   if (cmp_one < 0)
1539     {
1540       int m;
1541
1542       /* Number is less than one.  Pad significand with leading
1543          decimal zeros.  */
1544
1545       v = r;
1546       while (1)
1547         {
1548           /* Stop if we'd shift bits off the bottom.  */
1549           if (v.sig[0] & 7)
1550             break;
1551
1552           do_multiply (&u, &v, ten);
1553
1554           /* Stop if we're now >= 1.  */
1555           if (u.exp > 0)
1556             break;
1557
1558           v = u;
1559           dec_exp -= 1;
1560         }
1561       r = v;
1562
1563       /* Find power of 10.  Do this by multiplying in P=10**2**M when
1564          the current remainder is smaller than 1/P.  Fill PTEN with the
1565          power of 10 that we compute.  */
1566       m = floor_log2 ((int)(-r.exp * M_LOG10_2)) + 1;
1567       do
1568         {
1569           const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1570           const REAL_VALUE_TYPE *ptenmtwo = ten_to_mptwo (m);
1571
1572           if (do_compare (&v, ptenmtwo, 0) <= 0)
1573             {
1574               do_multiply (&v, &v, ptentwo);
1575               do_multiply (&pten, &pten, ptentwo);
1576               dec_exp -= 1 << m;
1577             }
1578         }
1579       while (--m >= 0);
1580
1581       /* Invert the positive power of 10 that we've collected so far.  */
1582       do_divide (&pten, one, &pten);
1583     }
1584
1585   p = str;
1586   if (sign)
1587     *p++ = '-';
1588   first = p++;
1589
1590   /* At this point, PTEN should contain the nearest power of 10 smaller
1591      than R, such that this division produces the first digit.
1592
1593      Using a divide-step primitive that returns the complete integral
1594      remainder avoids the rounding error that would be produced if
1595      we were to use do_divide here and then simply multiply by 10 for
1596      each subsequent digit.  */
1597
1598   digit = rtd_divmod (&r, &pten);
1599
1600   /* Be prepared for error in that division via underflow ...  */
1601   if (digit == 0 && cmp_significand_0 (&r))
1602     {
1603       /* Multiply by 10 and try again.  */
1604       do_multiply (&r, &r, ten);
1605       digit = rtd_divmod (&r, &pten);
1606       dec_exp -= 1;
1607       if (digit == 0)
1608         abort ();
1609     }
1610
1611   /* ... or overflow.  */
1612   if (digit == 10)
1613     {
1614       *p++ = '1';
1615       if (--digits > 0)
1616         *p++ = '0';
1617       dec_exp += 1;
1618     }
1619   else if (digit > 10)
1620     abort ();
1621   else
1622     *p++ = digit + '0';
1623
1624   /* Generate subsequent digits.  */
1625   while (--digits > 0)
1626     {
1627       do_multiply (&r, &r, ten);
1628       digit = rtd_divmod (&r, &pten);
1629       *p++ = digit + '0';
1630     }
1631   last = p;
1632
1633   /* Generate one more digit with which to do rounding.  */
1634   do_multiply (&r, &r, ten);
1635   digit = rtd_divmod (&r, &pten);
1636
1637   /* Round the result.  */
1638   if (digit == 5)
1639     {
1640       /* Round to nearest.  If R is nonzero there are additional
1641          nonzero digits to be extracted.  */
1642       if (cmp_significand_0 (&r))
1643         digit++;
1644       /* Round to even.  */
1645       else if ((p[-1] - '0') & 1)
1646         digit++;
1647     }
1648   if (digit > 5)
1649     {
1650       while (p > first)
1651         {
1652           digit = *--p;
1653           if (digit == '9')
1654             *p = '0';
1655           else
1656             {
1657               *p = digit + 1;
1658               break;
1659             }
1660         }
1661
1662       /* Carry out of the first digit.  This means we had all 9's and
1663          now have all 0's.  "Prepend" a 1 by overwriting the first 0.  */
1664       if (p == first)
1665         {
1666           first[1] = '1';
1667           dec_exp++;
1668         }
1669     }
1670
1671   /* Insert the decimal point.  */
1672   first[0] = first[1];
1673   first[1] = '.';
1674
1675   /* If requested, drop trailing zeros.  Never crop past "1.0".  */
1676   if (crop_trailing_zeros)
1677     while (last > first + 3 && last[-1] == '0')
1678       last--;
1679
1680   /* Append the exponent.  */
1681   sprintf (last, "e%+d", dec_exp);
1682 }
1683
1684 /* Render R as a hexadecimal floating point constant.  Emit DIGITS
1685    significant digits in the result, bounded by BUF_SIZE.  If DIGITS is 0,
1686    choose the maximum for the representation.  If CROP_TRAILING_ZEROS,
1687    strip trailing zeros.  */
1688
1689 void
1690 real_to_hexadecimal (char *str, const REAL_VALUE_TYPE *r, size_t buf_size,
1691                      size_t digits, int crop_trailing_zeros)
1692 {
1693   int i, j, exp = r->exp;
1694   char *p, *first;
1695   char exp_buf[16];
1696   size_t max_digits;
1697
1698   switch (r->class)
1699     {
1700     case rvc_zero:
1701       exp = 0;
1702       break;
1703     case rvc_normal:
1704       break;
1705     case rvc_inf:
1706       strcpy (str, (r->sign ? "-Inf" : "+Inf"));
1707       return;
1708     case rvc_nan:
1709       /* ??? Print the significand as well, if not canonical?  */
1710       strcpy (str, (r->sign ? "-NaN" : "+NaN"));
1711       return;
1712     default:
1713       abort ();
1714     }
1715
1716   if (digits == 0)
1717     digits = SIGNIFICAND_BITS / 4;
1718
1719   /* Bound the number of digits printed by the size of the output buffer.  */
1720
1721   sprintf (exp_buf, "p%+d", exp);
1722   max_digits = buf_size - strlen (exp_buf) - r->sign - 4 - 1;
1723   if (max_digits > buf_size)
1724     abort ();
1725   if (digits > max_digits)
1726     digits = max_digits;
1727
1728   p = str;
1729   if (r->sign)
1730     *p++ = '-';
1731   *p++ = '0';
1732   *p++ = 'x';
1733   *p++ = '0';
1734   *p++ = '.';
1735   first = p;
1736
1737   for (i = SIGSZ - 1; i >= 0; --i)
1738     for (j = HOST_BITS_PER_LONG - 4; j >= 0; j -= 4)
1739       {
1740         *p++ = "0123456789abcdef"[(r->sig[i] >> j) & 15];
1741         if (--digits == 0)
1742           goto out;
1743       }
1744
1745  out:
1746   if (crop_trailing_zeros)
1747     while (p > first + 1 && p[-1] == '0')
1748       p--;
1749
1750   sprintf (p, "p%+d", exp);
1751 }
1752
1753 /* Initialize R from a decimal or hexadecimal string.  The string is
1754    assumed to have been syntax checked already.  */
1755
1756 void
1757 real_from_string (REAL_VALUE_TYPE *r, const char *str)
1758 {
1759   int exp = 0;
1760   bool sign = false;
1761
1762   get_zero (r, 0);
1763
1764   if (*str == '-')
1765     {
1766       sign = true;
1767       str++;
1768     }
1769   else if (*str == '+')
1770     str++;
1771
1772   if (str[0] == '0' && str[1] == 'x')
1773     {
1774       /* Hexadecimal floating point.  */
1775       int pos = SIGNIFICAND_BITS - 4, d;
1776
1777       str += 2;
1778
1779       while (*str == '0')
1780         str++;
1781       while (1)
1782         {
1783           d = hex_value (*str);
1784           if (d == _hex_bad)
1785             break;
1786           if (pos >= 0)
1787             {
1788               r->sig[pos / HOST_BITS_PER_LONG]
1789                 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1790               pos -= 4;
1791             }
1792           exp += 4;
1793           str++;
1794         }
1795       if (*str == '.')
1796         {
1797           str++;
1798           if (pos == SIGNIFICAND_BITS - 4)
1799             {
1800               while (*str == '0')
1801                 str++, exp -= 4;
1802             }
1803           while (1)
1804             {
1805               d = hex_value (*str);
1806               if (d == _hex_bad)
1807                 break;
1808               if (pos >= 0)
1809                 {
1810                   r->sig[pos / HOST_BITS_PER_LONG]
1811                     |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1812                   pos -= 4;
1813                 }
1814               str++;
1815             }
1816         }
1817       if (*str == 'p' || *str == 'P')
1818         {
1819           bool exp_neg = false;
1820
1821           str++;
1822           if (*str == '-')
1823             {
1824               exp_neg = true;
1825               str++;
1826             }
1827           else if (*str == '+')
1828             str++;
1829
1830           d = 0;
1831           while (ISDIGIT (*str))
1832             {
1833               d *= 10;
1834               d += *str - '0';
1835               if (d > MAX_EXP)
1836                 {
1837                   /* Overflowed the exponent.  */
1838                   if (exp_neg)
1839                     goto underflow;
1840                   else
1841                     goto overflow;
1842                 }
1843               str++;
1844             }
1845           if (exp_neg)
1846             d = -d;
1847
1848           exp += d;
1849         }
1850
1851       r->class = rvc_normal;
1852       r->exp = exp;
1853
1854       normalize (r);
1855     }
1856   else
1857     {
1858       /* Decimal floating point.  */
1859       const REAL_VALUE_TYPE *ten = ten_to_ptwo (0);
1860       int d;
1861
1862       while (*str == '0')
1863         str++;
1864       while (ISDIGIT (*str))
1865         {
1866           d = *str++ - '0';
1867           do_multiply (r, r, ten);
1868           if (d)
1869             do_add (r, r, real_digit (d), 0);
1870         }
1871       if (*str == '.')
1872         {
1873           str++;
1874           if (r->class == rvc_zero)
1875             {
1876               while (*str == '0')
1877                 str++, exp--;
1878             }
1879           while (ISDIGIT (*str))
1880             {
1881               d = *str++ - '0';
1882               do_multiply (r, r, ten);
1883               if (d)
1884                 do_add (r, r, real_digit (d), 0);
1885               exp--;
1886             }
1887         }
1888
1889       if (*str == 'e' || *str == 'E')
1890         {
1891           bool exp_neg = false;
1892
1893           str++;
1894           if (*str == '-')
1895             {
1896               exp_neg = true;
1897               str++;
1898             }
1899           else if (*str == '+')
1900             str++;
1901
1902           d = 0;
1903           while (ISDIGIT (*str))
1904             {
1905               d *= 10;
1906               d += *str - '0';
1907               if (d > MAX_EXP)
1908                 {
1909                   /* Overflowed the exponent.  */
1910                   if (exp_neg)
1911                     goto underflow;
1912                   else
1913                     goto overflow;
1914                 }
1915               str++;
1916             }
1917           if (exp_neg)
1918             d = -d;
1919           exp += d;
1920         }
1921
1922       if (exp)
1923         times_pten (r, exp);
1924     }
1925
1926   r->sign = sign;
1927   return;
1928
1929  underflow:
1930   get_zero (r, sign);
1931   return;
1932
1933  overflow:
1934   get_inf (r, sign);
1935   return;
1936 }
1937
1938 /* Legacy.  Similar, but return the result directly.  */
1939
1940 REAL_VALUE_TYPE
1941 real_from_string2 (const char *s, enum machine_mode mode)
1942 {
1943   REAL_VALUE_TYPE r;
1944
1945   real_from_string (&r, s);
1946   if (mode != VOIDmode)
1947     real_convert (&r, mode, &r);
1948
1949   return r;
1950 }
1951
1952 /* Initialize R from the integer pair HIGH+LOW.  */
1953
1954 void
1955 real_from_integer (REAL_VALUE_TYPE *r, enum machine_mode mode,
1956                    unsigned HOST_WIDE_INT low, HOST_WIDE_INT high,
1957                    int unsigned_p)
1958 {
1959   if (low == 0 && high == 0)
1960     get_zero (r, 0);
1961   else
1962     {
1963       r->class = rvc_normal;
1964       r->sign = high < 0 && !unsigned_p;
1965       r->exp = 2 * HOST_BITS_PER_WIDE_INT;
1966
1967       if (r->sign)
1968         {
1969           high = ~high;
1970           if (low == 0)
1971             high += 1;
1972           else
1973             low = -low;
1974         }
1975
1976       if (HOST_BITS_PER_LONG == HOST_BITS_PER_WIDE_INT)
1977         {
1978           r->sig[SIGSZ-1] = high;
1979           r->sig[SIGSZ-2] = low;
1980           memset (r->sig, 0, sizeof(long)*(SIGSZ-2));
1981         }
1982       else if (HOST_BITS_PER_LONG*2 == HOST_BITS_PER_WIDE_INT)
1983         {
1984           r->sig[SIGSZ-1] = high >> (HOST_BITS_PER_LONG - 1) >> 1;
1985           r->sig[SIGSZ-2] = high;
1986           r->sig[SIGSZ-3] = low >> (HOST_BITS_PER_LONG - 1) >> 1;
1987           r->sig[SIGSZ-4] = low;
1988           if (SIGSZ > 4)
1989             memset (r->sig, 0, sizeof(long)*(SIGSZ-4));
1990         }
1991       else
1992         abort ();
1993
1994       normalize (r);
1995     }
1996
1997   if (mode != VOIDmode)
1998     real_convert (r, mode, r);
1999 }
2000
2001 /* Returns 10**2**N.  */
2002
2003 static const REAL_VALUE_TYPE *
2004 ten_to_ptwo (int n)
2005 {
2006   static REAL_VALUE_TYPE tens[EXP_BITS];
2007
2008   if (n < 0 || n >= EXP_BITS)
2009     abort ();
2010
2011   if (tens[n].class == rvc_zero)
2012     {
2013       if (n < (HOST_BITS_PER_WIDE_INT == 64 ? 5 : 4))
2014         {
2015           HOST_WIDE_INT t = 10;
2016           int i;
2017
2018           for (i = 0; i < n; ++i)
2019             t *= t;
2020
2021           real_from_integer (&tens[n], VOIDmode, t, 0, 1);
2022         }
2023       else
2024         {
2025           const REAL_VALUE_TYPE *t = ten_to_ptwo (n - 1);
2026           do_multiply (&tens[n], t, t);
2027         }
2028     }
2029
2030   return &tens[n];
2031 }
2032
2033 /* Returns 10**(-2**N).  */
2034
2035 static const REAL_VALUE_TYPE *
2036 ten_to_mptwo (int n)
2037 {
2038   static REAL_VALUE_TYPE tens[EXP_BITS];
2039
2040   if (n < 0 || n >= EXP_BITS)
2041     abort ();
2042
2043   if (tens[n].class == rvc_zero)
2044     do_divide (&tens[n], real_digit (1), ten_to_ptwo (n));
2045
2046   return &tens[n];
2047 }
2048
2049 /* Returns N.  */
2050
2051 static const REAL_VALUE_TYPE *
2052 real_digit (int n)
2053 {
2054   static REAL_VALUE_TYPE num[10];
2055
2056   if (n < 0 || n > 9)
2057     abort ();
2058
2059   if (n > 0 && num[n].class == rvc_zero)
2060     real_from_integer (&num[n], VOIDmode, n, 0, 1);
2061
2062   return &num[n];
2063 }
2064
2065 /* Multiply R by 10**EXP.  */
2066
2067 static void
2068 times_pten (REAL_VALUE_TYPE *r, int exp)
2069 {
2070   REAL_VALUE_TYPE pten, *rr;
2071   bool negative = (exp < 0);
2072   int i;
2073
2074   if (negative)
2075     {
2076       exp = -exp;
2077       pten = *real_digit (1);
2078       rr = &pten;
2079     }
2080   else
2081     rr = r;
2082
2083   for (i = 0; exp > 0; ++i, exp >>= 1)
2084     if (exp & 1)
2085       do_multiply (rr, rr, ten_to_ptwo (i));
2086
2087   if (negative)
2088     do_divide (r, r, &pten);
2089 }
2090
2091 /* Fills R with +Inf.  */
2092
2093 void
2094 real_inf (REAL_VALUE_TYPE *r)
2095 {
2096   get_inf (r, 0);
2097 }
2098
2099 /* Fills R with a NaN whose significand is described by STR.  If QUIET,
2100    we force a QNaN, else we force an SNaN.  The string, if not empty,
2101    is parsed as a number and placed in the significand.  Return true
2102    if the string was successfully parsed.  */
2103
2104 bool
2105 real_nan (REAL_VALUE_TYPE *r, const char *str, int quiet,
2106           enum machine_mode mode)
2107 {
2108   const struct real_format *fmt;
2109
2110   fmt = REAL_MODE_FORMAT (mode);
2111   if (fmt == NULL)
2112     abort ();
2113
2114   if (*str == 0)
2115     {
2116       if (quiet)
2117         get_canonical_qnan (r, 0);
2118       else
2119         get_canonical_snan (r, 0);
2120     }
2121   else
2122     {
2123       int base = 10, d;
2124       bool neg = false;
2125
2126       memset (r, 0, sizeof (*r));
2127       r->class = rvc_nan;
2128
2129       /* Parse akin to strtol into the significand of R.  */
2130
2131       while (ISSPACE (*str))
2132         str++;
2133       if (*str == '-')
2134         str++, neg = true;
2135       else if (*str == '+')
2136         str++;
2137       if (*str == '0')
2138         {
2139           if (*++str == 'x')
2140             str++, base = 16;
2141           else
2142             base = 8;
2143         }
2144
2145       while ((d = hex_value (*str)) < base)
2146         {
2147           REAL_VALUE_TYPE u;
2148
2149           switch (base)
2150             {
2151             case 8:
2152               lshift_significand (r, r, 3);
2153               break;
2154             case 16:
2155               lshift_significand (r, r, 4);
2156               break;
2157             case 10:
2158               lshift_significand_1 (&u, r);
2159               lshift_significand (r, r, 3);
2160               add_significands (r, r, &u);
2161               break;
2162             default:
2163               abort ();
2164             }
2165
2166           get_zero (&u, 0);
2167           u.sig[0] = d;
2168           add_significands (r, r, &u);
2169
2170           str++;
2171         }
2172
2173       /* Must have consumed the entire string for success.  */
2174       if (*str != 0)
2175         return false;
2176
2177       /* Shift the significand into place such that the bits
2178          are in the most significant bits for the format.  */
2179       lshift_significand (r, r, SIGNIFICAND_BITS - fmt->pnan);
2180
2181       /* Our MSB is always unset for NaNs.  */
2182       r->sig[SIGSZ-1] &= ~SIG_MSB;
2183
2184       /* Force quiet or signalling NaN.  */
2185       r->signalling = !quiet;
2186     }
2187
2188   return true;
2189 }
2190
2191 /* Fills R with the largest finite value representable in mode MODE.
2192    If SIGN is nonzero, R is set to the most negative finite value.  */
2193
2194 void
2195 real_maxval (REAL_VALUE_TYPE *r, int sign, enum machine_mode mode)
2196 {
2197   const struct real_format *fmt;
2198   int np2;
2199
2200   fmt = REAL_MODE_FORMAT (mode);
2201   if (fmt == NULL)
2202     abort ();
2203
2204   r->class = rvc_normal;
2205   r->sign = sign;
2206   r->signalling = 0;
2207   r->canonical = 0;
2208   r->exp = fmt->emax * fmt->log2_b;
2209
2210   np2 = SIGNIFICAND_BITS - fmt->p * fmt->log2_b;
2211   memset (r->sig, -1, SIGSZ * sizeof (unsigned long));
2212   clear_significand_below (r, np2);
2213 }
2214
2215 /* Fills R with 2**N.  */
2216
2217 void
2218 real_2expN (REAL_VALUE_TYPE *r, int n)
2219 {
2220   memset (r, 0, sizeof (*r));
2221
2222   n++;
2223   if (n > MAX_EXP)
2224     r->class = rvc_inf;
2225   else if (n < -MAX_EXP)
2226     ;
2227   else
2228     {
2229       r->class = rvc_normal;
2230       r->exp = n;
2231       r->sig[SIGSZ-1] = SIG_MSB;
2232     }
2233 }
2234
2235 \f
2236 static void
2237 round_for_format (const struct real_format *fmt, REAL_VALUE_TYPE *r)
2238 {
2239   int p2, np2, i, w;
2240   unsigned long sticky;
2241   bool guard, lsb;
2242   int emin2m1, emax2;
2243
2244   p2 = fmt->p * fmt->log2_b;
2245   emin2m1 = (fmt->emin - 1) * fmt->log2_b;
2246   emax2 = fmt->emax * fmt->log2_b;
2247
2248   np2 = SIGNIFICAND_BITS - p2;
2249   switch (r->class)
2250     {
2251     underflow:
2252       get_zero (r, r->sign);
2253     case rvc_zero:
2254       if (!fmt->has_signed_zero)
2255         r->sign = 0;
2256       return;
2257
2258     overflow:
2259       get_inf (r, r->sign);
2260     case rvc_inf:
2261       return;
2262
2263     case rvc_nan:
2264       clear_significand_below (r, np2);
2265       return;
2266
2267     case rvc_normal:
2268       break;
2269
2270     default:
2271       abort ();
2272     }
2273
2274   /* If we're not base2, normalize the exponent to a multiple of
2275      the true base.  */
2276   if (fmt->log2_b != 1)
2277     {
2278       int shift = r->exp & (fmt->log2_b - 1);
2279       if (shift)
2280         {
2281           shift = fmt->log2_b - shift;
2282           r->sig[0] |= sticky_rshift_significand (r, r, shift);
2283           r->exp += shift;
2284         }
2285     }
2286
2287   /* Check the range of the exponent.  If we're out of range,
2288      either underflow or overflow.  */
2289   if (r->exp > emax2)
2290     goto overflow;
2291   else if (r->exp <= emin2m1)
2292     {
2293       int diff;
2294
2295       if (!fmt->has_denorm)
2296         {
2297           /* Don't underflow completely until we've had a chance to round.  */
2298           if (r->exp < emin2m1)
2299             goto underflow;
2300         }
2301       else
2302         {
2303           diff = emin2m1 - r->exp + 1;
2304           if (diff > p2)
2305             goto underflow;
2306
2307           /* De-normalize the significand.  */
2308           r->sig[0] |= sticky_rshift_significand (r, r, diff);
2309           r->exp += diff;
2310         }
2311     }
2312
2313   /* There are P2 true significand bits, followed by one guard bit,
2314      followed by one sticky bit, followed by stuff.  Fold nonzero
2315      stuff into the sticky bit.  */
2316
2317   sticky = 0;
2318   for (i = 0, w = (np2 - 1) / HOST_BITS_PER_LONG; i < w; ++i)
2319     sticky |= r->sig[i];
2320   sticky |=
2321     r->sig[w] & (((unsigned long)1 << ((np2 - 1) % HOST_BITS_PER_LONG)) - 1);
2322
2323   guard = test_significand_bit (r, np2 - 1);
2324   lsb = test_significand_bit (r, np2);
2325
2326   /* Round to even.  */
2327   if (guard && (sticky || lsb))
2328     {
2329       REAL_VALUE_TYPE u;
2330       get_zero (&u, 0);
2331       set_significand_bit (&u, np2);
2332
2333       if (add_significands (r, r, &u))
2334         {
2335           /* Overflow.  Means the significand had been all ones, and
2336              is now all zeros.  Need to increase the exponent, and
2337              possibly re-normalize it.  */
2338           if (++r->exp > emax2)
2339             goto overflow;
2340           r->sig[SIGSZ-1] = SIG_MSB;
2341
2342           if (fmt->log2_b != 1)
2343             {
2344               int shift = r->exp & (fmt->log2_b - 1);
2345               if (shift)
2346                 {
2347                   shift = fmt->log2_b - shift;
2348                   rshift_significand (r, r, shift);
2349                   r->exp += shift;
2350                   if (r->exp > emax2)
2351                     goto overflow;
2352                 }
2353             }
2354         }
2355     }
2356
2357   /* Catch underflow that we deferred until after rounding.  */
2358   if (r->exp <= emin2m1)
2359     goto underflow;
2360
2361   /* Clear out trailing garbage.  */
2362   clear_significand_below (r, np2);
2363 }
2364
2365 /* Extend or truncate to a new mode.  */
2366
2367 void
2368 real_convert (REAL_VALUE_TYPE *r, enum machine_mode mode,
2369               const REAL_VALUE_TYPE *a)
2370 {
2371   const struct real_format *fmt;
2372
2373   fmt = REAL_MODE_FORMAT (mode);
2374   if (fmt == NULL)
2375     abort ();
2376
2377   *r = *a;
2378   round_for_format (fmt, r);
2379
2380   /* round_for_format de-normalizes denormals.  Undo just that part.  */
2381   if (r->class == rvc_normal)
2382     normalize (r);
2383 }
2384
2385 /* Legacy.  Likewise, except return the struct directly.  */
2386
2387 REAL_VALUE_TYPE
2388 real_value_truncate (enum machine_mode mode, REAL_VALUE_TYPE a)
2389 {
2390   REAL_VALUE_TYPE r;
2391   real_convert (&r, mode, &a);
2392   return r;
2393 }
2394
2395 /* Return true if truncating to MODE is exact.  */
2396
2397 bool
2398 exact_real_truncate (enum machine_mode mode, const REAL_VALUE_TYPE *a)
2399 {
2400   REAL_VALUE_TYPE t;
2401   real_convert (&t, mode, a);
2402   return real_identical (&t, a);
2403 }
2404
2405 /* Write R to the given target format.  Place the words of the result
2406    in target word order in BUF.  There are always 32 bits in each
2407    long, no matter the size of the host long.
2408
2409    Legacy: return word 0 for implementing REAL_VALUE_TO_TARGET_SINGLE.  */
2410
2411 long
2412 real_to_target_fmt (long *buf, const REAL_VALUE_TYPE *r_orig,
2413                     const struct real_format *fmt)
2414 {
2415   REAL_VALUE_TYPE r;
2416   long buf1;
2417
2418   r = *r_orig;
2419   round_for_format (fmt, &r);
2420
2421   if (!buf)
2422     buf = &buf1;
2423   (*fmt->encode) (fmt, buf, &r);
2424
2425   return *buf;
2426 }
2427
2428 /* Similar, but look up the format from MODE.  */
2429
2430 long
2431 real_to_target (long *buf, const REAL_VALUE_TYPE *r, enum machine_mode mode)
2432 {
2433   const struct real_format *fmt;
2434
2435   fmt = REAL_MODE_FORMAT (mode);
2436   if (fmt == NULL)
2437     abort ();
2438
2439   return real_to_target_fmt (buf, r, fmt);
2440 }
2441
2442 /* Read R from the given target format.  Read the words of the result
2443    in target word order in BUF.  There are always 32 bits in each
2444    long, no matter the size of the host long.  */
2445
2446 void
2447 real_from_target_fmt (REAL_VALUE_TYPE *r, const long *buf,
2448                       const struct real_format *fmt)
2449 {
2450   (*fmt->decode) (fmt, r, buf);
2451 }
2452
2453 /* Similar, but look up the format from MODE.  */
2454
2455 void
2456 real_from_target (REAL_VALUE_TYPE *r, const long *buf, enum machine_mode mode)
2457 {
2458   const struct real_format *fmt;
2459
2460   fmt = REAL_MODE_FORMAT (mode);
2461   if (fmt == NULL)
2462     abort ();
2463
2464   (*fmt->decode) (fmt, r, buf);
2465 }
2466
2467 /* Return the number of bits in the significand for MODE.  */
2468 /* ??? Legacy.  Should get access to real_format directly.  */
2469
2470 int
2471 significand_size (enum machine_mode mode)
2472 {
2473   const struct real_format *fmt;
2474
2475   fmt = REAL_MODE_FORMAT (mode);
2476   if (fmt == NULL)
2477     return 0;
2478
2479   return fmt->p * fmt->log2_b;
2480 }
2481
2482 /* Return a hash value for the given real value.  */
2483 /* ??? The "unsigned int" return value is intended to be hashval_t,
2484    but I didn't want to pull hashtab.h into real.h.  */
2485
2486 unsigned int
2487 real_hash (const REAL_VALUE_TYPE *r)
2488 {
2489   unsigned int h;
2490   size_t i;
2491
2492   h = r->class | (r->sign << 2);
2493   switch (r->class)
2494     {
2495     case rvc_zero:
2496     case rvc_inf:
2497       return h;
2498
2499     case rvc_normal:
2500       h |= r->exp << 3;
2501       break;
2502
2503     case rvc_nan:
2504       if (r->signalling)
2505         h ^= (unsigned int)-1;
2506       if (r->canonical)
2507         return h;
2508       break;
2509
2510     default:
2511       abort ();
2512     }
2513
2514   if (sizeof(unsigned long) > sizeof(unsigned int))
2515     for (i = 0; i < SIGSZ; ++i)
2516       {
2517         unsigned long s = r->sig[i];
2518         h ^= s ^ (s >> (HOST_BITS_PER_LONG / 2));
2519       }
2520   else
2521     for (i = 0; i < SIGSZ; ++i)
2522       h ^= r->sig[i];
2523
2524   return h;
2525 }
2526 \f
2527 /* IEEE single-precision format.  */
2528
2529 static void encode_ieee_single (const struct real_format *fmt,
2530                                 long *, const REAL_VALUE_TYPE *);
2531 static void decode_ieee_single (const struct real_format *,
2532                                 REAL_VALUE_TYPE *, const long *);
2533
2534 static void
2535 encode_ieee_single (const struct real_format *fmt, long *buf,
2536                     const REAL_VALUE_TYPE *r)
2537 {
2538   unsigned long image, sig, exp;
2539   unsigned long sign = r->sign;
2540   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2541
2542   image = sign << 31;
2543   sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
2544
2545   switch (r->class)
2546     {
2547     case rvc_zero:
2548       break;
2549
2550     case rvc_inf:
2551       if (fmt->has_inf)
2552         image |= 255 << 23;
2553       else
2554         image |= 0x7fffffff;
2555       break;
2556
2557     case rvc_nan:
2558       if (fmt->has_nans)
2559         {
2560           if (r->canonical)
2561             sig = 0;
2562           if (r->signalling == fmt->qnan_msb_set)
2563             sig &= ~(1 << 22);
2564           else
2565             sig |= 1 << 22;
2566           /* We overload qnan_msb_set here: it's only clear for
2567              mips_ieee_single, which wants all mantissa bits but the
2568              quiet/signalling one set in canonical NaNs (at least
2569              Quiet ones).  */
2570           if (r->canonical && !fmt->qnan_msb_set)
2571             sig |= (1 << 22) - 1;
2572           else if (sig == 0)
2573             sig = 1 << 21;
2574
2575           image |= 255 << 23;
2576           image |= sig;
2577         }
2578       else
2579         image |= 0x7fffffff;
2580       break;
2581
2582     case rvc_normal:
2583       /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2584          whereas the intermediate representation is 0.F x 2**exp.
2585          Which means we're off by one.  */
2586       if (denormal)
2587         exp = 0;
2588       else
2589       exp = r->exp + 127 - 1;
2590       image |= exp << 23;
2591       image |= sig;
2592       break;
2593
2594     default:
2595       abort ();
2596     }
2597
2598   buf[0] = image;
2599 }
2600
2601 static void
2602 decode_ieee_single (const struct real_format *fmt, REAL_VALUE_TYPE *r,
2603                     const long *buf)
2604 {
2605   unsigned long image = buf[0] & 0xffffffff;
2606   bool sign = (image >> 31) & 1;
2607   int exp = (image >> 23) & 0xff;
2608
2609   memset (r, 0, sizeof (*r));
2610   image <<= HOST_BITS_PER_LONG - 24;
2611   image &= ~SIG_MSB;
2612
2613   if (exp == 0)
2614     {
2615       if (image && fmt->has_denorm)
2616         {
2617           r->class = rvc_normal;
2618           r->sign = sign;
2619           r->exp = -126;
2620           r->sig[SIGSZ-1] = image << 1;
2621           normalize (r);
2622         }
2623       else if (fmt->has_signed_zero)
2624         r->sign = sign;
2625     }
2626   else if (exp == 255 && (fmt->has_nans || fmt->has_inf))
2627     {
2628       if (image)
2629         {
2630           r->class = rvc_nan;
2631           r->sign = sign;
2632           r->signalling = (((image >> (HOST_BITS_PER_LONG - 2)) & 1)
2633                            ^ fmt->qnan_msb_set);
2634           r->sig[SIGSZ-1] = image;
2635         }
2636       else
2637         {
2638           r->class = rvc_inf;
2639           r->sign = sign;
2640         }
2641     }
2642   else
2643     {
2644       r->class = rvc_normal;
2645       r->sign = sign;
2646       r->exp = exp - 127 + 1;
2647       r->sig[SIGSZ-1] = image | SIG_MSB;
2648     }
2649 }
2650
2651 const struct real_format ieee_single_format =
2652   {
2653     encode_ieee_single,
2654     decode_ieee_single,
2655     2,
2656     1,
2657     24,
2658     24,
2659     -125,
2660     128,
2661     31,
2662     true,
2663     true,
2664     true,
2665     true,
2666     true
2667   };
2668
2669 const struct real_format mips_single_format =
2670   {
2671     encode_ieee_single,
2672     decode_ieee_single,
2673     2,
2674     1,
2675     24,
2676     24,
2677     -125,
2678     128,
2679     31,
2680     true,
2681     true,
2682     true,
2683     true,
2684     false
2685   };
2686
2687 \f
2688 /* IEEE double-precision format.  */
2689
2690 static void encode_ieee_double (const struct real_format *fmt,
2691                                 long *, const REAL_VALUE_TYPE *);
2692 static void decode_ieee_double (const struct real_format *,
2693                                 REAL_VALUE_TYPE *, const long *);
2694
2695 static void
2696 encode_ieee_double (const struct real_format *fmt, long *buf,
2697                     const REAL_VALUE_TYPE *r)
2698 {
2699   unsigned long image_lo, image_hi, sig_lo, sig_hi, exp;
2700   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2701
2702   image_hi = r->sign << 31;
2703   image_lo = 0;
2704
2705   if (HOST_BITS_PER_LONG == 64)
2706     {
2707       sig_hi = r->sig[SIGSZ-1];
2708       sig_lo = (sig_hi >> (64 - 53)) & 0xffffffff;
2709       sig_hi = (sig_hi >> (64 - 53 + 1) >> 31) & 0xfffff;
2710     }
2711   else
2712     {
2713       sig_hi = r->sig[SIGSZ-1];
2714       sig_lo = r->sig[SIGSZ-2];
2715       sig_lo = (sig_hi << 21) | (sig_lo >> 11);
2716       sig_hi = (sig_hi >> 11) & 0xfffff;
2717     }
2718
2719   switch (r->class)
2720     {
2721     case rvc_zero:
2722       break;
2723
2724     case rvc_inf:
2725       if (fmt->has_inf)
2726         image_hi |= 2047 << 20;
2727       else
2728         {
2729           image_hi |= 0x7fffffff;
2730           image_lo = 0xffffffff;
2731         }
2732       break;
2733
2734     case rvc_nan:
2735       if (fmt->has_nans)
2736         {
2737           if (r->canonical)
2738             sig_hi = sig_lo = 0;
2739           if (r->signalling == fmt->qnan_msb_set)
2740             sig_hi &= ~(1 << 19);
2741           else
2742             sig_hi |= 1 << 19;
2743           /* We overload qnan_msb_set here: it's only clear for
2744              mips_ieee_single, which wants all mantissa bits but the
2745              quiet/signalling one set in canonical NaNs (at least
2746              Quiet ones).  */
2747           if (r->canonical && !fmt->qnan_msb_set)
2748             {
2749               sig_hi |= (1 << 19) - 1;
2750               sig_lo = 0xffffffff;
2751             }
2752           else if (sig_hi == 0 && sig_lo == 0)
2753             sig_hi = 1 << 18;
2754
2755           image_hi |= 2047 << 20;
2756           image_hi |= sig_hi;
2757           image_lo = sig_lo;
2758         }
2759       else
2760         {
2761           image_hi |= 0x7fffffff;
2762           image_lo = 0xffffffff;
2763         }
2764       break;
2765
2766     case rvc_normal:
2767       /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2768          whereas the intermediate representation is 0.F x 2**exp.
2769          Which means we're off by one.  */
2770       if (denormal)
2771         exp = 0;
2772       else
2773         exp = r->exp + 1023 - 1;
2774       image_hi |= exp << 20;
2775       image_hi |= sig_hi;
2776       image_lo = sig_lo;
2777       break;
2778
2779     default:
2780       abort ();
2781     }
2782
2783   if (FLOAT_WORDS_BIG_ENDIAN)
2784     buf[0] = image_hi, buf[1] = image_lo;
2785   else
2786     buf[0] = image_lo, buf[1] = image_hi;
2787 }
2788
2789 static void
2790 decode_ieee_double (const struct real_format *fmt, REAL_VALUE_TYPE *r,
2791                     const long *buf)
2792 {
2793   unsigned long image_hi, image_lo;
2794   bool sign;
2795   int exp;
2796
2797   if (FLOAT_WORDS_BIG_ENDIAN)
2798     image_hi = buf[0], image_lo = buf[1];
2799   else
2800     image_lo = buf[0], image_hi = buf[1];
2801   image_lo &= 0xffffffff;
2802   image_hi &= 0xffffffff;
2803
2804   sign = (image_hi >> 31) & 1;
2805   exp = (image_hi >> 20) & 0x7ff;
2806
2807   memset (r, 0, sizeof (*r));
2808
2809   image_hi <<= 32 - 21;
2810   image_hi |= image_lo >> 21;
2811   image_hi &= 0x7fffffff;
2812   image_lo <<= 32 - 21;
2813
2814   if (exp == 0)
2815     {
2816       if ((image_hi || image_lo) && fmt->has_denorm)
2817         {
2818           r->class = rvc_normal;
2819           r->sign = sign;
2820           r->exp = -1022;
2821           if (HOST_BITS_PER_LONG == 32)
2822             {
2823               image_hi = (image_hi << 1) | (image_lo >> 31);
2824               image_lo <<= 1;
2825               r->sig[SIGSZ-1] = image_hi;
2826               r->sig[SIGSZ-2] = image_lo;
2827             }
2828           else
2829             {
2830               image_hi = (image_hi << 31 << 2) | (image_lo << 1);
2831               r->sig[SIGSZ-1] = image_hi;
2832             }
2833           normalize (r);
2834         }
2835       else if (fmt->has_signed_zero)
2836         r->sign = sign;
2837     }
2838   else if (exp == 2047 && (fmt->has_nans || fmt->has_inf))
2839     {
2840       if (image_hi || image_lo)
2841         {
2842           r->class = rvc_nan;
2843           r->sign = sign;
2844           r->signalling = ((image_hi >> 30) & 1) ^ fmt->qnan_msb_set;
2845           if (HOST_BITS_PER_LONG == 32)
2846             {
2847               r->sig[SIGSZ-1] = image_hi;
2848               r->sig[SIGSZ-2] = image_lo;
2849             }
2850           else
2851             r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo;
2852         }
2853       else
2854         {
2855           r->class = rvc_inf;
2856           r->sign = sign;
2857         }
2858     }
2859   else
2860     {
2861       r->class = rvc_normal;
2862       r->sign = sign;
2863       r->exp = exp - 1023 + 1;
2864       if (HOST_BITS_PER_LONG == 32)
2865         {
2866           r->sig[SIGSZ-1] = image_hi | SIG_MSB;
2867           r->sig[SIGSZ-2] = image_lo;
2868         }
2869       else
2870         r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo | SIG_MSB;
2871     }
2872 }
2873
2874 const struct real_format ieee_double_format =
2875   {
2876     encode_ieee_double,
2877     decode_ieee_double,
2878     2,
2879     1,
2880     53,
2881     53,
2882     -1021,
2883     1024,
2884     63,
2885     true,
2886     true,
2887     true,
2888     true,
2889     true
2890   };
2891
2892 const struct real_format mips_double_format =
2893   {
2894     encode_ieee_double,
2895     decode_ieee_double,
2896     2,
2897     1,
2898     53,
2899     53,
2900     -1021,
2901     1024,
2902     63,
2903     true,
2904     true,
2905     true,
2906     true,
2907     false
2908   };
2909
2910 \f
2911 /* IEEE extended double precision format.  This comes in three
2912    flavors: Intel's as a 12 byte image, Intel's as a 16 byte image,
2913    and Motorola's.  */
2914
2915 static void encode_ieee_extended (const struct real_format *fmt,
2916                                   long *, const REAL_VALUE_TYPE *);
2917 static void decode_ieee_extended (const struct real_format *,
2918                                   REAL_VALUE_TYPE *, const long *);
2919
2920 static void encode_ieee_extended_128 (const struct real_format *fmt,
2921                                       long *, const REAL_VALUE_TYPE *);
2922 static void decode_ieee_extended_128 (const struct real_format *,
2923                                       REAL_VALUE_TYPE *, const long *);
2924
2925 static void
2926 encode_ieee_extended (const struct real_format *fmt, long *buf,
2927                       const REAL_VALUE_TYPE *r)
2928 {
2929   unsigned long image_hi, sig_hi, sig_lo;
2930   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2931
2932   image_hi = r->sign << 15;
2933   sig_hi = sig_lo = 0;
2934
2935   switch (r->class)
2936     {
2937     case rvc_zero:
2938       break;
2939
2940     case rvc_inf:
2941       if (fmt->has_inf)
2942         {
2943           image_hi |= 32767;
2944
2945           /* Intel requires the explicit integer bit to be set, otherwise
2946              it considers the value a "pseudo-infinity".  Motorola docs
2947              say it doesn't care.  */
2948           sig_hi = 0x80000000;
2949         }
2950       else
2951         {
2952           image_hi |= 32767;
2953           sig_lo = sig_hi = 0xffffffff;
2954         }
2955       break;
2956
2957     case rvc_nan:
2958       if (fmt->has_nans)
2959         {
2960           image_hi |= 32767;
2961           if (HOST_BITS_PER_LONG == 32)
2962             {
2963               sig_hi = r->sig[SIGSZ-1];
2964               sig_lo = r->sig[SIGSZ-2];
2965             }
2966           else
2967             {
2968               sig_lo = r->sig[SIGSZ-1];
2969               sig_hi = sig_lo >> 31 >> 1;
2970               sig_lo &= 0xffffffff;
2971             }
2972           if (r->signalling == fmt->qnan_msb_set)
2973             sig_hi &= ~(1 << 30);
2974           else
2975             sig_hi |= 1 << 30;
2976           if ((sig_hi & 0x7fffffff) == 0 && sig_lo == 0)
2977             sig_hi = 1 << 29;
2978
2979           /* Intel requires the explicit integer bit to be set, otherwise
2980              it considers the value a "pseudo-nan".  Motorola docs say it
2981              doesn't care.  */
2982           sig_hi |= 0x80000000;
2983         }
2984       else
2985         {
2986           image_hi |= 32767;
2987           sig_lo = sig_hi = 0xffffffff;
2988         }
2989       break;
2990
2991     case rvc_normal:
2992       {
2993         int exp = r->exp;
2994
2995         /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2996            whereas the intermediate representation is 0.F x 2**exp.
2997            Which means we're off by one.
2998
2999            Except for Motorola, which consider exp=0 and explicit
3000            integer bit set to continue to be normalized.  In theory
3001            this discrepancy has been taken care of by the difference
3002            in fmt->emin in round_for_format.  */
3003
3004         if (denormal)
3005           exp = 0;
3006         else
3007           {
3008             exp += 16383 - 1;
3009             if (exp < 0)
3010               abort ();
3011           }
3012         image_hi |= exp;
3013
3014         if (HOST_BITS_PER_LONG == 32)
3015           {
3016             sig_hi = r->sig[SIGSZ-1];
3017             sig_lo = r->sig[SIGSZ-2];
3018           }
3019         else
3020           {
3021             sig_lo = r->sig[SIGSZ-1];
3022             sig_hi = sig_lo >> 31 >> 1;
3023             sig_lo &= 0xffffffff;
3024           }
3025       }
3026       break;
3027
3028     default:
3029       abort ();
3030     }
3031
3032   if (FLOAT_WORDS_BIG_ENDIAN)
3033     buf[0] = image_hi << 16, buf[1] = sig_hi, buf[2] = sig_lo;
3034   else
3035     buf[0] = sig_lo, buf[1] = sig_hi, buf[2] = image_hi;
3036 }
3037
3038 static void
3039 encode_ieee_extended_128 (const struct real_format *fmt, long *buf,
3040                           const REAL_VALUE_TYPE *r)
3041 {
3042   buf[3 * !FLOAT_WORDS_BIG_ENDIAN] = 0;
3043   encode_ieee_extended (fmt, buf+!!FLOAT_WORDS_BIG_ENDIAN, r);
3044 }
3045
3046 static void
3047 decode_ieee_extended (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3048                       const long *buf)
3049 {
3050   unsigned long image_hi, sig_hi, sig_lo;
3051   bool sign;
3052   int exp;
3053
3054   if (FLOAT_WORDS_BIG_ENDIAN)
3055     image_hi = buf[0] >> 16, sig_hi = buf[1], sig_lo = buf[2];
3056   else
3057     sig_lo = buf[0], sig_hi = buf[1], image_hi = buf[2];
3058   sig_lo &= 0xffffffff;
3059   sig_hi &= 0xffffffff;
3060   image_hi &= 0xffffffff;
3061
3062   sign = (image_hi >> 15) & 1;
3063   exp = image_hi & 0x7fff;
3064
3065   memset (r, 0, sizeof (*r));
3066
3067   if (exp == 0)
3068     {
3069       if ((sig_hi || sig_lo) && fmt->has_denorm)
3070         {
3071           r->class = rvc_normal;
3072           r->sign = sign;
3073
3074           /* When the IEEE format contains a hidden bit, we know that
3075              it's zero at this point, and so shift up the significand
3076              and decrease the exponent to match.  In this case, Motorola
3077              defines the explicit integer bit to be valid, so we don't
3078              know whether the msb is set or not.  */
3079           r->exp = fmt->emin;
3080           if (HOST_BITS_PER_LONG == 32)
3081             {
3082               r->sig[SIGSZ-1] = sig_hi;
3083               r->sig[SIGSZ-2] = sig_lo;
3084             }
3085           else
3086             r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3087
3088           normalize (r);
3089         }
3090       else if (fmt->has_signed_zero)
3091         r->sign = sign;
3092     }
3093   else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3094     {
3095       /* See above re "pseudo-infinities" and "pseudo-nans".
3096          Short summary is that the MSB will likely always be
3097          set, and that we don't care about it.  */
3098       sig_hi &= 0x7fffffff;
3099
3100       if (sig_hi || sig_lo)
3101         {
3102           r->class = rvc_nan;
3103           r->sign = sign;
3104           r->signalling = ((sig_hi >> 30) & 1) ^ fmt->qnan_msb_set;
3105           if (HOST_BITS_PER_LONG == 32)
3106             {
3107               r->sig[SIGSZ-1] = sig_hi;
3108               r->sig[SIGSZ-2] = sig_lo;
3109             }
3110           else
3111             r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3112         }
3113       else
3114         {
3115           r->class = rvc_inf;
3116           r->sign = sign;
3117         }
3118     }
3119   else
3120     {
3121       r->class = rvc_normal;
3122       r->sign = sign;
3123       r->exp = exp - 16383 + 1;
3124       if (HOST_BITS_PER_LONG == 32)
3125         {
3126           r->sig[SIGSZ-1] = sig_hi;
3127           r->sig[SIGSZ-2] = sig_lo;
3128         }
3129       else
3130         r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3131     }
3132 }
3133
3134 static void
3135 decode_ieee_extended_128 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3136                           const long *buf)
3137 {
3138   decode_ieee_extended (fmt, r, buf+!!FLOAT_WORDS_BIG_ENDIAN);
3139 }
3140
3141 const struct real_format ieee_extended_motorola_format =
3142   {
3143     encode_ieee_extended,
3144     decode_ieee_extended,
3145     2,
3146     1,
3147     64,
3148     64,
3149     -16382,
3150     16384,
3151     95,
3152     true,
3153     true,
3154     true,
3155     true,
3156     true
3157   };
3158
3159 const struct real_format ieee_extended_intel_96_format =
3160   {
3161     encode_ieee_extended,
3162     decode_ieee_extended,
3163     2,
3164     1,
3165     64,
3166     64,
3167     -16381,
3168     16384,
3169     79,
3170     true,
3171     true,
3172     true,
3173     true,
3174     true
3175   };
3176
3177 const struct real_format ieee_extended_intel_128_format =
3178   {
3179     encode_ieee_extended_128,
3180     decode_ieee_extended_128,
3181     2,
3182     1,
3183     64,
3184     64,
3185     -16381,
3186     16384,
3187     79,
3188     true,
3189     true,
3190     true,
3191     true,
3192     true
3193   };
3194
3195 /* The following caters to i386 systems that set the rounding precision
3196    to 53 bits instead of 64, e.g. FreeBSD.  */
3197 const struct real_format ieee_extended_intel_96_round_53_format =
3198   {
3199     encode_ieee_extended,
3200     decode_ieee_extended,
3201     2,
3202     1,
3203     53,
3204     53,
3205     -16381,
3206     16384,
3207     79,
3208     true,
3209     true,
3210     true,
3211     true,
3212     true
3213   };
3214 \f
3215 /* IBM 128-bit extended precision format: a pair of IEEE double precision
3216    numbers whose sum is equal to the extended precision value.  The number
3217    with greater magnitude is first.  This format has the same magnitude
3218    range as an IEEE double precision value, but effectively 106 bits of
3219    significand precision.  Infinity and NaN are represented by their IEEE
3220    double precision value stored in the first number, the second number is
3221    ignored.  Zeroes, Infinities, and NaNs are set in both doubles
3222    due to precedent.  */
3223
3224 static void encode_ibm_extended (const struct real_format *fmt,
3225                                  long *, const REAL_VALUE_TYPE *);
3226 static void decode_ibm_extended (const struct real_format *,
3227                                  REAL_VALUE_TYPE *, const long *);
3228
3229 static void
3230 encode_ibm_extended (const struct real_format *fmt, long *buf,
3231                      const REAL_VALUE_TYPE *r)
3232 {
3233   REAL_VALUE_TYPE u, v;
3234   const struct real_format *base_fmt;
3235
3236   base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3237
3238   switch (r->class)
3239     {
3240     case rvc_zero:
3241       /* Both doubles have sign bit set.  */
3242       buf[0] = FLOAT_WORDS_BIG_ENDIAN ? r->sign << 31 : 0;
3243       buf[1] = FLOAT_WORDS_BIG_ENDIAN ? 0 : r->sign << 31;
3244       buf[2] = buf[0];
3245       buf[3] = buf[1];
3246       break;
3247
3248     case rvc_inf:
3249     case rvc_nan:
3250       /* Both doubles set to Inf / NaN.  */
3251       encode_ieee_double (base_fmt, &buf[0], r);
3252       buf[2] = buf[0];
3253       buf[3] = buf[1];
3254       return;
3255
3256     case rvc_normal:
3257       /* u = IEEE double precision portion of significand.  */
3258       u = *r;
3259       clear_significand_below (&u, SIGNIFICAND_BITS - 53);
3260
3261       normalize (&u);
3262       /* If the upper double is zero, we have a denormal double, so
3263          move it to the first double and leave the second as zero.  */
3264       if (u.class == rvc_zero)
3265         {
3266           v = u;
3267           u = *r;
3268           normalize (&u);
3269         }
3270       else
3271         {
3272           /* v = remainder containing additional 53 bits of significand.  */
3273           do_add (&v, r, &u, 1);
3274           round_for_format (base_fmt, &v);
3275         }
3276
3277       round_for_format (base_fmt, &u);
3278
3279       encode_ieee_double (base_fmt, &buf[0], &u);
3280       encode_ieee_double (base_fmt, &buf[2], &v);
3281       break;
3282
3283     default:
3284       abort ();
3285     }
3286 }
3287
3288 static void
3289 decode_ibm_extended (const struct real_format *fmt ATTRIBUTE_UNUSED, REAL_VALUE_TYPE *r,
3290                      const long *buf)
3291 {
3292   REAL_VALUE_TYPE u, v;
3293   const struct real_format *base_fmt;
3294
3295   base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3296   decode_ieee_double (base_fmt, &u, &buf[0]);
3297
3298   if (u.class != rvc_zero && u.class != rvc_inf && u.class != rvc_nan)
3299     {
3300       decode_ieee_double (base_fmt, &v, &buf[2]);
3301       do_add (r, &u, &v, 0);
3302     }
3303   else
3304     *r = u;
3305 }
3306
3307 const struct real_format ibm_extended_format =
3308   {
3309     encode_ibm_extended,
3310     decode_ibm_extended,
3311     2,
3312     1,
3313     53 + 53,
3314     53,
3315     -1021 + 53,
3316     1024,
3317     -1,
3318     true,
3319     true,
3320     true,
3321     true,
3322     true
3323   };
3324
3325 const struct real_format mips_extended_format =
3326   {
3327     encode_ibm_extended,
3328     decode_ibm_extended,
3329     2,
3330     1,
3331     53 + 53,
3332     53,
3333     -1021 + 53,
3334     1024,
3335     -1,
3336     true,
3337     true,
3338     true,
3339     true,
3340     false
3341   };
3342
3343 \f
3344 /* IEEE quad precision format.  */
3345
3346 static void encode_ieee_quad (const struct real_format *fmt,
3347                               long *, const REAL_VALUE_TYPE *);
3348 static void decode_ieee_quad (const struct real_format *,
3349                               REAL_VALUE_TYPE *, const long *);
3350
3351 static void
3352 encode_ieee_quad (const struct real_format *fmt, long *buf,
3353                   const REAL_VALUE_TYPE *r)
3354 {
3355   unsigned long image3, image2, image1, image0, exp;
3356   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3357   REAL_VALUE_TYPE u;
3358
3359   image3 = r->sign << 31;
3360   image2 = 0;
3361   image1 = 0;
3362   image0 = 0;
3363
3364   rshift_significand (&u, r, SIGNIFICAND_BITS - 113);
3365
3366   switch (r->class)
3367     {
3368     case rvc_zero:
3369       break;
3370
3371     case rvc_inf:
3372       if (fmt->has_inf)
3373         image3 |= 32767 << 16;
3374       else
3375         {
3376           image3 |= 0x7fffffff;
3377           image2 = 0xffffffff;
3378           image1 = 0xffffffff;
3379           image0 = 0xffffffff;
3380         }
3381       break;
3382
3383     case rvc_nan:
3384       if (fmt->has_nans)
3385         {
3386           image3 |= 32767 << 16;
3387
3388           if (r->canonical)
3389             {
3390               /* Don't use bits from the significand.  The
3391                  initialization above is right.  */
3392             }
3393           else if (HOST_BITS_PER_LONG == 32)
3394             {
3395               image0 = u.sig[0];
3396               image1 = u.sig[1];
3397               image2 = u.sig[2];
3398               image3 |= u.sig[3] & 0xffff;
3399             }
3400           else
3401             {
3402               image0 = u.sig[0];
3403               image1 = image0 >> 31 >> 1;
3404               image2 = u.sig[1];
3405               image3 |= (image2 >> 31 >> 1) & 0xffff;
3406               image0 &= 0xffffffff;
3407               image2 &= 0xffffffff;
3408             }
3409           if (r->signalling == fmt->qnan_msb_set)
3410             image3 &= ~0x8000;
3411           else
3412             image3 |= 0x8000;
3413           /* We overload qnan_msb_set here: it's only clear for
3414              mips_ieee_single, which wants all mantissa bits but the
3415              quiet/signalling one set in canonical NaNs (at least
3416              Quiet ones).  */
3417           if (r->canonical && !fmt->qnan_msb_set)
3418             {
3419               image3 |= 0x7fff;
3420               image2 = image1 = image0 = 0xffffffff;
3421             }
3422           else if (((image3 & 0xffff) | image2 | image1 | image0) == 0)
3423             image3 |= 0x4000;
3424         }
3425       else
3426         {
3427           image3 |= 0x7fffffff;
3428           image2 = 0xffffffff;
3429           image1 = 0xffffffff;
3430           image0 = 0xffffffff;
3431         }
3432       break;
3433
3434     case rvc_normal:
3435       /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3436          whereas the intermediate representation is 0.F x 2**exp.
3437          Which means we're off by one.  */
3438       if (denormal)
3439         exp = 0;
3440       else
3441         exp = r->exp + 16383 - 1;
3442       image3 |= exp << 16;
3443
3444       if (HOST_BITS_PER_LONG == 32)
3445         {
3446           image0 = u.sig[0];
3447           image1 = u.sig[1];
3448           image2 = u.sig[2];
3449           image3 |= u.sig[3] & 0xffff;
3450         }
3451       else
3452         {
3453           image0 = u.sig[0];
3454           image1 = image0 >> 31 >> 1;
3455           image2 = u.sig[1];
3456           image3 |= (image2 >> 31 >> 1) & 0xffff;
3457           image0 &= 0xffffffff;
3458           image2 &= 0xffffffff;
3459         }
3460       break;
3461
3462     default:
3463       abort ();
3464     }
3465
3466   if (FLOAT_WORDS_BIG_ENDIAN)
3467     {
3468       buf[0] = image3;
3469       buf[1] = image2;
3470       buf[2] = image1;
3471       buf[3] = image0;
3472     }
3473   else
3474     {
3475       buf[0] = image0;
3476       buf[1] = image1;
3477       buf[2] = image2;
3478       buf[3] = image3;
3479     }
3480 }
3481
3482 static void
3483 decode_ieee_quad (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3484                   const long *buf)
3485 {
3486   unsigned long image3, image2, image1, image0;
3487   bool sign;
3488   int exp;
3489
3490   if (FLOAT_WORDS_BIG_ENDIAN)
3491     {
3492       image3 = buf[0];
3493       image2 = buf[1];
3494       image1 = buf[2];
3495       image0 = buf[3];
3496     }
3497   else
3498     {
3499       image0 = buf[0];
3500       image1 = buf[1];
3501       image2 = buf[2];
3502       image3 = buf[3];
3503     }
3504   image0 &= 0xffffffff;
3505   image1 &= 0xffffffff;
3506   image2 &= 0xffffffff;
3507
3508   sign = (image3 >> 31) & 1;
3509   exp = (image3 >> 16) & 0x7fff;
3510   image3 &= 0xffff;
3511
3512   memset (r, 0, sizeof (*r));
3513
3514   if (exp == 0)
3515     {
3516       if ((image3 | image2 | image1 | image0) && fmt->has_denorm)
3517         {
3518           r->class = rvc_normal;
3519           r->sign = sign;
3520
3521           r->exp = -16382 + (SIGNIFICAND_BITS - 112);
3522           if (HOST_BITS_PER_LONG == 32)
3523             {
3524               r->sig[0] = image0;
3525               r->sig[1] = image1;
3526               r->sig[2] = image2;
3527               r->sig[3] = image3;
3528             }
3529           else
3530             {
3531               r->sig[0] = (image1 << 31 << 1) | image0;
3532               r->sig[1] = (image3 << 31 << 1) | image2;
3533             }
3534
3535           normalize (r);
3536         }
3537       else if (fmt->has_signed_zero)
3538         r->sign = sign;
3539     }
3540   else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3541     {
3542       if (image3 | image2 | image1 | image0)
3543         {
3544           r->class = rvc_nan;
3545           r->sign = sign;
3546           r->signalling = ((image3 >> 15) & 1) ^ fmt->qnan_msb_set;
3547
3548           if (HOST_BITS_PER_LONG == 32)
3549             {
3550               r->sig[0] = image0;
3551               r->sig[1] = image1;
3552               r->sig[2] = image2;
3553               r->sig[3] = image3;
3554             }
3555           else
3556             {
3557               r->sig[0] = (image1 << 31 << 1) | image0;
3558               r->sig[1] = (image3 << 31 << 1) | image2;
3559             }
3560           lshift_significand (r, r, SIGNIFICAND_BITS - 113);
3561         }
3562       else
3563         {
3564           r->class = rvc_inf;
3565           r->sign = sign;
3566         }
3567     }
3568   else
3569     {
3570       r->class = rvc_normal;
3571       r->sign = sign;
3572       r->exp = exp - 16383 + 1;
3573
3574       if (HOST_BITS_PER_LONG == 32)
3575         {
3576           r->sig[0] = image0;
3577           r->sig[1] = image1;
3578           r->sig[2] = image2;
3579           r->sig[3] = image3;
3580         }
3581       else
3582         {
3583           r->sig[0] = (image1 << 31 << 1) | image0;
3584           r->sig[1] = (image3 << 31 << 1) | image2;
3585         }
3586       lshift_significand (r, r, SIGNIFICAND_BITS - 113);
3587       r->sig[SIGSZ-1] |= SIG_MSB;
3588     }
3589 }
3590
3591 const struct real_format ieee_quad_format =
3592   {
3593     encode_ieee_quad,
3594     decode_ieee_quad,
3595     2,
3596     1,
3597     113,
3598     113,
3599     -16381,
3600     16384,
3601     127,
3602     true,
3603     true,
3604     true,
3605     true,
3606     true
3607   };
3608
3609 const struct real_format mips_quad_format =
3610   {
3611     encode_ieee_quad,
3612     decode_ieee_quad,
3613     2,
3614     1,
3615     113,
3616     113,
3617     -16381,
3618     16384,
3619     127,
3620     true,
3621     true,
3622     true,
3623     true,
3624     false
3625   };
3626 \f
3627 /* Descriptions of VAX floating point formats can be found beginning at
3628
3629    http://h71000.www7.hp.com/doc/73FINAL/4515/4515pro_013.html#f_floating_point_format
3630
3631    The thing to remember is that they're almost IEEE, except for word
3632    order, exponent bias, and the lack of infinities, nans, and denormals.
3633
3634    We don't implement the H_floating format here, simply because neither
3635    the VAX or Alpha ports use it.  */
3636
3637 static void encode_vax_f (const struct real_format *fmt,
3638                           long *, const REAL_VALUE_TYPE *);
3639 static void decode_vax_f (const struct real_format *,
3640                           REAL_VALUE_TYPE *, const long *);
3641 static void encode_vax_d (const struct real_format *fmt,
3642                           long *, const REAL_VALUE_TYPE *);
3643 static void decode_vax_d (const struct real_format *,
3644                           REAL_VALUE_TYPE *, const long *);
3645 static void encode_vax_g (const struct real_format *fmt,
3646                           long *, const REAL_VALUE_TYPE *);
3647 static void decode_vax_g (const struct real_format *,
3648                           REAL_VALUE_TYPE *, const long *);
3649
3650 static void
3651 encode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3652               const REAL_VALUE_TYPE *r)
3653 {
3654   unsigned long sign, exp, sig, image;
3655
3656   sign = r->sign << 15;
3657
3658   switch (r->class)
3659     {
3660     case rvc_zero:
3661       image = 0;
3662       break;
3663
3664     case rvc_inf:
3665     case rvc_nan:
3666       image = 0xffff7fff | sign;
3667       break;
3668
3669     case rvc_normal:
3670       sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
3671       exp = r->exp + 128;
3672
3673       image = (sig << 16) & 0xffff0000;
3674       image |= sign;
3675       image |= exp << 7;
3676       image |= sig >> 16;
3677       break;
3678
3679     default:
3680       abort ();
3681     }
3682
3683   buf[0] = image;
3684 }
3685
3686 static void
3687 decode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED,
3688               REAL_VALUE_TYPE *r, const long *buf)
3689 {
3690   unsigned long image = buf[0] & 0xffffffff;
3691   int exp = (image >> 7) & 0xff;
3692
3693   memset (r, 0, sizeof (*r));
3694
3695   if (exp != 0)
3696     {
3697       r->class = rvc_normal;
3698       r->sign = (image >> 15) & 1;
3699       r->exp = exp - 128;
3700
3701       image = ((image & 0x7f) << 16) | ((image >> 16) & 0xffff);
3702       r->sig[SIGSZ-1] = (image << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
3703     }
3704 }
3705
3706 static void
3707 encode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3708               const REAL_VALUE_TYPE *r)
3709 {
3710   unsigned long image0, image1, sign = r->sign << 15;
3711
3712   switch (r->class)
3713     {
3714     case rvc_zero:
3715       image0 = image1 = 0;
3716       break;
3717
3718     case rvc_inf:
3719     case rvc_nan:
3720       image0 = 0xffff7fff | sign;
3721       image1 = 0xffffffff;
3722       break;
3723
3724     case rvc_normal:
3725       /* Extract the significand into straight hi:lo.  */
3726       if (HOST_BITS_PER_LONG == 64)
3727         {
3728           image0 = r->sig[SIGSZ-1];
3729           image1 = (image0 >> (64 - 56)) & 0xffffffff;
3730           image0 = (image0 >> (64 - 56 + 1) >> 31) & 0x7fffff;
3731         }
3732       else
3733         {
3734           image0 = r->sig[SIGSZ-1];
3735           image1 = r->sig[SIGSZ-2];
3736           image1 = (image0 << 24) | (image1 >> 8);
3737           image0 = (image0 >> 8) & 0xffffff;
3738         }
3739
3740       /* Rearrange the half-words of the significand to match the
3741          external format.  */
3742       image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff007f;
3743       image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
3744
3745       /* Add the sign and exponent.  */
3746       image0 |= sign;
3747       image0 |= (r->exp + 128) << 7;
3748       break;
3749
3750     default:
3751       abort ();
3752     }
3753
3754   if (FLOAT_WORDS_BIG_ENDIAN)
3755     buf[0] = image1, buf[1] = image0;
3756   else
3757     buf[0] = image0, buf[1] = image1;
3758 }
3759
3760 static void
3761 decode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED,
3762               REAL_VALUE_TYPE *r, const long *buf)
3763 {
3764   unsigned long image0, image1;
3765   int exp;
3766
3767   if (FLOAT_WORDS_BIG_ENDIAN)
3768     image1 = buf[0], image0 = buf[1];
3769   else
3770     image0 = buf[0], image1 = buf[1];
3771   image0 &= 0xffffffff;
3772   image1 &= 0xffffffff;
3773
3774   exp = (image0 >> 7) & 0xff;
3775
3776   memset (r, 0, sizeof (*r));
3777
3778   if (exp != 0)
3779     {
3780       r->class = rvc_normal;
3781       r->sign = (image0 >> 15) & 1;
3782       r->exp = exp - 128;
3783
3784       /* Rearrange the half-words of the external format into
3785          proper ascending order.  */
3786       image0 = ((image0 & 0x7f) << 16) | ((image0 >> 16) & 0xffff);
3787       image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);