Merge branch 'master' of ssh://crater.dragonflybsd.org/repository/git/dragonfly
[dragonfly.git] / contrib / gcc-3.4 / gcc / real.c
1 /* real.c - software floating point emulation.
2    Copyright (C) 1993, 1994, 1995, 1996, 1997, 1998, 1999,
3    2000, 2002, 2003, 2004 Free Software Foundation, Inc.
4    Contributed by Stephen L. Moshier (moshier@world.std.com).
5    Re-written by Richard Henderson <rth@redhat.com>
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 27.
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   /* Zero out the remaining fields.  */
644   r->signalling = 0;
645   r->canonical = 0;
646
647   /* Re-normalize the result.  */
648   normalize (r);
649
650   /* Special case: if the subtraction results in zero, the result
651      is positive.  */
652   if (r->class == rvc_zero)
653     r->sign = 0;
654   else
655     r->sig[0] |= inexact;
656
657   return inexact;
658 }
659
660 /* Calculate R = A * B.  Return true if the result may be inexact.  */
661
662 static bool
663 do_multiply (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
664              const REAL_VALUE_TYPE *b)
665 {
666   REAL_VALUE_TYPE u, t, *rr;
667   unsigned int i, j, k;
668   int sign = a->sign ^ b->sign;
669   bool inexact = false;
670
671   switch (CLASS2 (a->class, b->class))
672     {
673     case CLASS2 (rvc_zero, rvc_zero):
674     case CLASS2 (rvc_zero, rvc_normal):
675     case CLASS2 (rvc_normal, rvc_zero):
676       /* +-0 * ANY = 0 with appropriate sign.  */
677       get_zero (r, sign);
678       return false;
679
680     case CLASS2 (rvc_zero, rvc_nan):
681     case CLASS2 (rvc_normal, rvc_nan):
682     case CLASS2 (rvc_inf, rvc_nan):
683     case CLASS2 (rvc_nan, rvc_nan):
684       /* ANY * NaN = NaN.  */
685       *r = *b;
686       r->sign = sign;
687       return false;
688
689     case CLASS2 (rvc_nan, rvc_zero):
690     case CLASS2 (rvc_nan, rvc_normal):
691     case CLASS2 (rvc_nan, rvc_inf):
692       /* NaN * ANY = NaN.  */
693       *r = *a;
694       r->sign = sign;
695       return false;
696
697     case CLASS2 (rvc_zero, rvc_inf):
698     case CLASS2 (rvc_inf, rvc_zero):
699       /* 0 * Inf = NaN */
700       get_canonical_qnan (r, sign);
701       return false;
702
703     case CLASS2 (rvc_inf, rvc_inf):
704     case CLASS2 (rvc_normal, rvc_inf):
705     case CLASS2 (rvc_inf, rvc_normal):
706       /* Inf * Inf = Inf, R * Inf = Inf */
707       get_inf (r, sign);
708       return false;
709
710     case CLASS2 (rvc_normal, rvc_normal):
711       break;
712
713     default:
714       abort ();
715     }
716
717   if (r == a || r == b)
718     rr = &t;
719   else
720     rr = r;
721   get_zero (rr, 0);
722
723   /* Collect all the partial products.  Since we don't have sure access
724      to a widening multiply, we split each long into two half-words.
725
726      Consider the long-hand form of a four half-word multiplication:
727
728                  A  B  C  D
729               *  E  F  G  H
730              --------------
731                 DE DF DG DH
732              CE CF CG CH
733           BE BF BG BH
734        AE AF AG AH
735
736      We construct partial products of the widened half-word products
737      that are known to not overlap, e.g. DF+DH.  Each such partial
738      product is given its proper exponent, which allows us to sum them
739      and obtain the finished product.  */
740
741   for (i = 0; i < SIGSZ * 2; ++i)
742     {
743       unsigned long ai = a->sig[i / 2];
744       if (i & 1)
745         ai >>= HOST_BITS_PER_LONG / 2;
746       else
747         ai &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
748
749       if (ai == 0)
750         continue;
751
752       for (j = 0; j < 2; ++j)
753         {
754           int exp = (a->exp - (2*SIGSZ-1-i)*(HOST_BITS_PER_LONG/2)
755                      + (b->exp - (1-j)*(HOST_BITS_PER_LONG/2)));
756
757           if (exp > MAX_EXP)
758             {
759               get_inf (r, sign);
760               return true;
761             }
762           if (exp < -MAX_EXP)
763             {
764               /* Would underflow to zero, which we shouldn't bother adding.  */
765               inexact = true;
766               continue;
767             }
768
769           memset (&u, 0, sizeof (u));
770           u.class = rvc_normal;
771           u.exp = exp;
772
773           for (k = j; k < SIGSZ * 2; k += 2)
774             {
775               unsigned long bi = b->sig[k / 2];
776               if (k & 1)
777                 bi >>= HOST_BITS_PER_LONG / 2;
778               else
779                 bi &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
780
781               u.sig[k / 2] = ai * bi;
782             }
783
784           normalize (&u);
785           inexact |= do_add (rr, rr, &u, 0);
786         }
787     }
788
789   rr->sign = sign;
790   if (rr != r)
791     *r = t;
792
793   return inexact;
794 }
795
796 /* Calculate R = A / B.  Return true if the result may be inexact.  */
797
798 static bool
799 do_divide (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
800            const REAL_VALUE_TYPE *b)
801 {
802   int exp, sign = a->sign ^ b->sign;
803   REAL_VALUE_TYPE t, *rr;
804   bool inexact;
805
806   switch (CLASS2 (a->class, b->class))
807     {
808     case CLASS2 (rvc_zero, rvc_zero):
809       /* 0 / 0 = NaN.  */
810     case CLASS2 (rvc_inf, rvc_inf):
811       /* Inf / Inf = NaN.  */
812       get_canonical_qnan (r, sign);
813       return false;
814
815     case CLASS2 (rvc_zero, rvc_normal):
816     case CLASS2 (rvc_zero, rvc_inf):
817       /* 0 / ANY = 0.  */
818     case CLASS2 (rvc_normal, rvc_inf):
819       /* R / Inf = 0.  */
820       get_zero (r, sign);
821       return false;
822
823     case CLASS2 (rvc_normal, rvc_zero):
824       /* R / 0 = Inf.  */
825     case CLASS2 (rvc_inf, rvc_zero):
826       /* Inf / 0 = Inf.  */
827       get_inf (r, sign);
828       return false;
829
830     case CLASS2 (rvc_zero, rvc_nan):
831     case CLASS2 (rvc_normal, rvc_nan):
832     case CLASS2 (rvc_inf, rvc_nan):
833     case CLASS2 (rvc_nan, rvc_nan):
834       /* ANY / NaN = NaN.  */
835       *r = *b;
836       r->sign = sign;
837       return false;
838
839     case CLASS2 (rvc_nan, rvc_zero):
840     case CLASS2 (rvc_nan, rvc_normal):
841     case CLASS2 (rvc_nan, rvc_inf):
842       /* NaN / ANY = NaN.  */
843       *r = *a;
844       r->sign = sign;
845       return false;
846
847     case CLASS2 (rvc_inf, rvc_normal):
848       /* Inf / R = Inf.  */
849       get_inf (r, sign);
850       return false;
851
852     case CLASS2 (rvc_normal, rvc_normal):
853       break;
854
855     default:
856       abort ();
857     }
858
859   if (r == a || r == b)
860     rr = &t;
861   else
862     rr = r;
863
864   /* Make sure all fields in the result are initialized.  */
865   get_zero (rr, 0);
866   rr->class = rvc_normal;
867   rr->sign = sign;
868
869   exp = a->exp - b->exp + 1;
870   if (exp > MAX_EXP)
871     {
872       get_inf (r, sign);
873       return true;
874     }
875   if (exp < -MAX_EXP)
876     {
877       get_zero (r, sign);
878       return true;
879     }
880   rr->exp = exp;
881
882   inexact = div_significands (rr, a, b);
883
884   /* Re-normalize the result.  */
885   normalize (rr);
886   rr->sig[0] |= inexact;
887
888   if (rr != r)
889     *r = t;
890
891   return inexact;
892 }
893
894 /* Return a tri-state comparison of A vs B.  Return NAN_RESULT if
895    one of the two operands is a NaN.  */
896
897 static int
898 do_compare (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b,
899             int nan_result)
900 {
901   int ret;
902
903   switch (CLASS2 (a->class, b->class))
904     {
905     case CLASS2 (rvc_zero, rvc_zero):
906       /* Sign of zero doesn't matter for compares.  */
907       return 0;
908
909     case CLASS2 (rvc_inf, rvc_zero):
910     case CLASS2 (rvc_inf, rvc_normal):
911     case CLASS2 (rvc_normal, rvc_zero):
912       return (a->sign ? -1 : 1);
913
914     case CLASS2 (rvc_inf, rvc_inf):
915       return -a->sign - -b->sign;
916
917     case CLASS2 (rvc_zero, rvc_normal):
918     case CLASS2 (rvc_zero, rvc_inf):
919     case CLASS2 (rvc_normal, rvc_inf):
920       return (b->sign ? 1 : -1);
921
922     case CLASS2 (rvc_zero, rvc_nan):
923     case CLASS2 (rvc_normal, rvc_nan):
924     case CLASS2 (rvc_inf, rvc_nan):
925     case CLASS2 (rvc_nan, rvc_nan):
926     case CLASS2 (rvc_nan, rvc_zero):
927     case CLASS2 (rvc_nan, rvc_normal):
928     case CLASS2 (rvc_nan, rvc_inf):
929       return nan_result;
930
931     case CLASS2 (rvc_normal, rvc_normal):
932       break;
933
934     default:
935       abort ();
936     }
937
938   if (a->sign != b->sign)
939     return -a->sign - -b->sign;
940
941   if (a->exp > b->exp)
942     ret = 1;
943   else if (a->exp < b->exp)
944     ret = -1;
945   else
946     ret = cmp_significands (a, b);
947
948   return (a->sign ? -ret : ret);
949 }
950
951 /* Return A truncated to an integral value toward zero.  */
952
953 static void
954 do_fix_trunc (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
955 {
956   *r = *a;
957
958   switch (r->class)
959     {
960     case rvc_zero:
961     case rvc_inf:
962     case rvc_nan:
963       break;
964
965     case rvc_normal:
966       if (r->exp <= 0)
967         get_zero (r, r->sign);
968       else if (r->exp < SIGNIFICAND_BITS)
969         clear_significand_below (r, SIGNIFICAND_BITS - r->exp);
970       break;
971
972     default:
973       abort ();
974     }
975 }
976
977 /* Perform the binary or unary operation described by CODE.
978    For a unary operation, leave OP1 NULL.  */
979
980 void
981 real_arithmetic (REAL_VALUE_TYPE *r, int icode, const REAL_VALUE_TYPE *op0,
982                  const REAL_VALUE_TYPE *op1)
983 {
984   enum tree_code code = icode;
985
986   switch (code)
987     {
988     case PLUS_EXPR:
989       do_add (r, op0, op1, 0);
990       break;
991
992     case MINUS_EXPR:
993       do_add (r, op0, op1, 1);
994       break;
995
996     case MULT_EXPR:
997       do_multiply (r, op0, op1);
998       break;
999
1000     case RDIV_EXPR:
1001       do_divide (r, op0, op1);
1002       break;
1003
1004     case MIN_EXPR:
1005       if (op1->class == rvc_nan)
1006         *r = *op1;
1007       else if (do_compare (op0, op1, -1) < 0)
1008         *r = *op0;
1009       else
1010         *r = *op1;
1011       break;
1012
1013     case MAX_EXPR:
1014       if (op1->class == rvc_nan)
1015         *r = *op1;
1016       else if (do_compare (op0, op1, 1) < 0)
1017         *r = *op1;
1018       else
1019         *r = *op0;
1020       break;
1021
1022     case NEGATE_EXPR:
1023       *r = *op0;
1024       r->sign ^= 1;
1025       break;
1026
1027     case ABS_EXPR:
1028       *r = *op0;
1029       r->sign = 0;
1030       break;
1031
1032     case FIX_TRUNC_EXPR:
1033       do_fix_trunc (r, op0);
1034       break;
1035
1036     default:
1037       abort ();
1038     }
1039 }
1040
1041 /* Legacy.  Similar, but return the result directly.  */
1042
1043 REAL_VALUE_TYPE
1044 real_arithmetic2 (int icode, const REAL_VALUE_TYPE *op0,
1045                   const REAL_VALUE_TYPE *op1)
1046 {
1047   REAL_VALUE_TYPE r;
1048   real_arithmetic (&r, icode, op0, op1);
1049   return r;
1050 }
1051
1052 bool
1053 real_compare (int icode, const REAL_VALUE_TYPE *op0,
1054               const REAL_VALUE_TYPE *op1)
1055 {
1056   enum tree_code code = icode;
1057
1058   switch (code)
1059     {
1060     case LT_EXPR:
1061       return do_compare (op0, op1, 1) < 0;
1062     case LE_EXPR:
1063       return do_compare (op0, op1, 1) <= 0;
1064     case GT_EXPR:
1065       return do_compare (op0, op1, -1) > 0;
1066     case GE_EXPR:
1067       return do_compare (op0, op1, -1) >= 0;
1068     case EQ_EXPR:
1069       return do_compare (op0, op1, -1) == 0;
1070     case NE_EXPR:
1071       return do_compare (op0, op1, -1) != 0;
1072     case UNORDERED_EXPR:
1073       return op0->class == rvc_nan || op1->class == rvc_nan;
1074     case ORDERED_EXPR:
1075       return op0->class != rvc_nan && op1->class != rvc_nan;
1076     case UNLT_EXPR:
1077       return do_compare (op0, op1, -1) < 0;
1078     case UNLE_EXPR:
1079       return do_compare (op0, op1, -1) <= 0;
1080     case UNGT_EXPR:
1081       return do_compare (op0, op1, 1) > 0;
1082     case UNGE_EXPR:
1083       return do_compare (op0, op1, 1) >= 0;
1084     case UNEQ_EXPR:
1085       return do_compare (op0, op1, 0) == 0;
1086
1087     default:
1088       abort ();
1089     }
1090 }
1091
1092 /* Return floor log2(R).  */
1093
1094 int
1095 real_exponent (const REAL_VALUE_TYPE *r)
1096 {
1097   switch (r->class)
1098     {
1099     case rvc_zero:
1100       return 0;
1101     case rvc_inf:
1102     case rvc_nan:
1103       return (unsigned int)-1 >> 1;
1104     case rvc_normal:
1105       return r->exp;
1106     default:
1107       abort ();
1108     }
1109 }
1110
1111 /* R = OP0 * 2**EXP.  */
1112
1113 void
1114 real_ldexp (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *op0, int exp)
1115 {
1116   *r = *op0;
1117   switch (r->class)
1118     {
1119     case rvc_zero:
1120     case rvc_inf:
1121     case rvc_nan:
1122       break;
1123
1124     case rvc_normal:
1125       exp += op0->exp;
1126       if (exp > MAX_EXP)
1127         get_inf (r, r->sign);
1128       else if (exp < -MAX_EXP)
1129         get_zero (r, r->sign);
1130       else
1131         r->exp = exp;
1132       break;
1133
1134     default:
1135       abort ();
1136     }
1137 }
1138
1139 /* Determine whether a floating-point value X is infinite.  */
1140
1141 bool
1142 real_isinf (const REAL_VALUE_TYPE *r)
1143 {
1144   return (r->class == rvc_inf);
1145 }
1146
1147 /* Determine whether a floating-point value X is a NaN.  */
1148
1149 bool
1150 real_isnan (const REAL_VALUE_TYPE *r)
1151 {
1152   return (r->class == rvc_nan);
1153 }
1154
1155 /* Determine whether a floating-point value X is negative.  */
1156
1157 bool
1158 real_isneg (const REAL_VALUE_TYPE *r)
1159 {
1160   return r->sign;
1161 }
1162
1163 /* Determine whether a floating-point value X is minus zero.  */
1164
1165 bool
1166 real_isnegzero (const REAL_VALUE_TYPE *r)
1167 {
1168   return r->sign && r->class == rvc_zero;
1169 }
1170
1171 /* Compare two floating-point objects for bitwise identity.  */
1172
1173 bool
1174 real_identical (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
1175 {
1176   int i;
1177
1178   if (a->class != b->class)
1179     return false;
1180   if (a->sign != b->sign)
1181     return false;
1182
1183   switch (a->class)
1184     {
1185     case rvc_zero:
1186     case rvc_inf:
1187       return true;
1188
1189     case rvc_normal:
1190       if (a->exp != b->exp)
1191         return false;
1192       break;
1193
1194     case rvc_nan:
1195       if (a->signalling != b->signalling)
1196         return false;
1197       /* The significand is ignored for canonical NaNs.  */
1198       if (a->canonical || b->canonical)
1199         return a->canonical == b->canonical;
1200       break;
1201
1202     default:
1203       abort ();
1204     }
1205
1206   for (i = 0; i < SIGSZ; ++i)
1207     if (a->sig[i] != b->sig[i])
1208       return false;
1209
1210   return true;
1211 }
1212
1213 /* Try to change R into its exact multiplicative inverse in machine
1214    mode MODE.  Return true if successful.  */
1215
1216 bool
1217 exact_real_inverse (enum machine_mode mode, REAL_VALUE_TYPE *r)
1218 {
1219   const REAL_VALUE_TYPE *one = real_digit (1);
1220   REAL_VALUE_TYPE u;
1221   int i;
1222
1223   if (r->class != rvc_normal)
1224     return false;
1225
1226   /* Check for a power of two: all significand bits zero except the MSB.  */
1227   for (i = 0; i < SIGSZ-1; ++i)
1228     if (r->sig[i] != 0)
1229       return false;
1230   if (r->sig[SIGSZ-1] != SIG_MSB)
1231     return false;
1232
1233   /* Find the inverse and truncate to the required mode.  */
1234   do_divide (&u, one, r);
1235   real_convert (&u, mode, &u);
1236
1237   /* The rounding may have overflowed.  */
1238   if (u.class != rvc_normal)
1239     return false;
1240   for (i = 0; i < SIGSZ-1; ++i)
1241     if (u.sig[i] != 0)
1242       return false;
1243   if (u.sig[SIGSZ-1] != SIG_MSB)
1244     return false;
1245
1246   *r = u;
1247   return true;
1248 }
1249 \f
1250 /* Render R as an integer.  */
1251
1252 HOST_WIDE_INT
1253 real_to_integer (const REAL_VALUE_TYPE *r)
1254 {
1255   unsigned HOST_WIDE_INT i;
1256
1257   switch (r->class)
1258     {
1259     case rvc_zero:
1260     underflow:
1261       return 0;
1262
1263     case rvc_inf:
1264     case rvc_nan:
1265     overflow:
1266       i = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1267       if (!r->sign)
1268         i--;
1269       return i;
1270
1271     case rvc_normal:
1272       if (r->exp <= 0)
1273         goto underflow;
1274       /* Only force overflow for unsigned overflow.  Signed overflow is
1275          undefined, so it doesn't matter what we return, and some callers
1276          expect to be able to use this routine for both signed and
1277          unsigned conversions.  */
1278       if (r->exp > HOST_BITS_PER_WIDE_INT)
1279         goto overflow;
1280
1281       if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1282         i = r->sig[SIGSZ-1];
1283       else if (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG)
1284         {
1285           i = r->sig[SIGSZ-1];
1286           i = i << (HOST_BITS_PER_LONG - 1) << 1;
1287           i |= r->sig[SIGSZ-2];
1288         }
1289       else
1290         abort ();
1291
1292       i >>= HOST_BITS_PER_WIDE_INT - r->exp;
1293
1294       if (r->sign)
1295         i = -i;
1296       return i;
1297
1298     default:
1299       abort ();
1300     }
1301 }
1302
1303 /* Likewise, but to an integer pair, HI+LOW.  */
1304
1305 void
1306 real_to_integer2 (HOST_WIDE_INT *plow, HOST_WIDE_INT *phigh,
1307                   const REAL_VALUE_TYPE *r)
1308 {
1309   REAL_VALUE_TYPE t;
1310   HOST_WIDE_INT low, high;
1311   int exp;
1312
1313   switch (r->class)
1314     {
1315     case rvc_zero:
1316     underflow:
1317       low = high = 0;
1318       break;
1319
1320     case rvc_inf:
1321     case rvc_nan:
1322     overflow:
1323       high = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1324       if (r->sign)
1325         low = 0;
1326       else
1327         {
1328           high--;
1329           low = -1;
1330         }
1331       break;
1332
1333     case rvc_normal:
1334       exp = r->exp;
1335       if (exp <= 0)
1336         goto underflow;
1337       /* Only force overflow for unsigned overflow.  Signed overflow is
1338          undefined, so it doesn't matter what we return, and some callers
1339          expect to be able to use this routine for both signed and
1340          unsigned conversions.  */
1341       if (exp > 2*HOST_BITS_PER_WIDE_INT)
1342         goto overflow;
1343
1344       rshift_significand (&t, r, 2*HOST_BITS_PER_WIDE_INT - exp);
1345       if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1346         {
1347           high = t.sig[SIGSZ-1];
1348           low = t.sig[SIGSZ-2];
1349         }
1350       else if (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG)
1351         {
1352           high = t.sig[SIGSZ-1];
1353           high = high << (HOST_BITS_PER_LONG - 1) << 1;
1354           high |= t.sig[SIGSZ-2];
1355
1356           low = t.sig[SIGSZ-3];
1357           low = low << (HOST_BITS_PER_LONG - 1) << 1;
1358           low |= t.sig[SIGSZ-4];
1359         }
1360       else
1361         abort ();
1362
1363       if (r->sign)
1364         {
1365           if (low == 0)
1366             high = -high;
1367           else
1368             low = -low, high = ~high;
1369         }
1370       break;
1371
1372     default:
1373       abort ();
1374     }
1375
1376   *plow = low;
1377   *phigh = high;
1378 }
1379
1380 /* A subroutine of real_to_decimal.  Compute the quotient and remainder
1381    of NUM / DEN.  Return the quotient and place the remainder in NUM.
1382    It is expected that NUM / DEN are close enough that the quotient is
1383    small.  */
1384
1385 static unsigned long
1386 rtd_divmod (REAL_VALUE_TYPE *num, REAL_VALUE_TYPE *den)
1387 {
1388   unsigned long q, msb;
1389   int expn = num->exp, expd = den->exp;
1390
1391   if (expn < expd)
1392     return 0;
1393
1394   q = msb = 0;
1395   goto start;
1396   do
1397     {
1398       msb = num->sig[SIGSZ-1] & SIG_MSB;
1399       q <<= 1;
1400       lshift_significand_1 (num, num);
1401     start:
1402       if (msb || cmp_significands (num, den) >= 0)
1403         {
1404           sub_significands (num, num, den, 0);
1405           q |= 1;
1406         }
1407     }
1408   while (--expn >= expd);
1409
1410   num->exp = expd;
1411   normalize (num);
1412
1413   return q;
1414 }
1415
1416 /* Render R as a decimal floating point constant.  Emit DIGITS significant
1417    digits in the result, bounded by BUF_SIZE.  If DIGITS is 0, choose the
1418    maximum for the representation.  If CROP_TRAILING_ZEROS, strip trailing
1419    zeros.  */
1420
1421 #define M_LOG10_2       0.30102999566398119521
1422
1423 void
1424 real_to_decimal (char *str, const REAL_VALUE_TYPE *r_orig, size_t buf_size,
1425                  size_t digits, int crop_trailing_zeros)
1426 {
1427   const REAL_VALUE_TYPE *one, *ten;
1428   REAL_VALUE_TYPE r, pten, u, v;
1429   int dec_exp, cmp_one, digit;
1430   size_t max_digits;
1431   char *p, *first, *last;
1432   bool sign;
1433
1434   r = *r_orig;
1435   switch (r.class)
1436     {
1437     case rvc_zero:
1438       strcpy (str, (r.sign ? "-0.0" : "0.0"));
1439       return;
1440     case rvc_normal:
1441       break;
1442     case rvc_inf:
1443       strcpy (str, (r.sign ? "-Inf" : "+Inf"));
1444       return;
1445     case rvc_nan:
1446       /* ??? Print the significand as well, if not canonical?  */
1447       strcpy (str, (r.sign ? "-NaN" : "+NaN"));
1448       return;
1449     default:
1450       abort ();
1451     }
1452
1453   /* Bound the number of digits printed by the size of the representation.  */
1454   max_digits = SIGNIFICAND_BITS * M_LOG10_2;
1455   if (digits == 0 || digits > max_digits)
1456     digits = max_digits;
1457
1458   /* Estimate the decimal exponent, and compute the length of the string it
1459      will print as.  Be conservative and add one to account for possible
1460      overflow or rounding error.  */
1461   dec_exp = r.exp * M_LOG10_2;
1462   for (max_digits = 1; dec_exp ; max_digits++)
1463     dec_exp /= 10;
1464
1465   /* Bound the number of digits printed by the size of the output buffer.  */
1466   max_digits = buf_size - 1 - 1 - 2 - max_digits - 1;
1467   if (max_digits > buf_size)
1468     abort ();
1469   if (digits > max_digits)
1470     digits = max_digits;
1471
1472   one = real_digit (1);
1473   ten = ten_to_ptwo (0);
1474
1475   sign = r.sign;
1476   r.sign = 0;
1477
1478   dec_exp = 0;
1479   pten = *one;
1480
1481   cmp_one = do_compare (&r, one, 0);
1482   if (cmp_one > 0)
1483     {
1484       int m;
1485
1486       /* Number is greater than one.  Convert significand to an integer
1487          and strip trailing decimal zeros.  */
1488
1489       u = r;
1490       u.exp = SIGNIFICAND_BITS - 1;
1491
1492       /* Largest M, such that 10**2**M fits within SIGNIFICAND_BITS.  */
1493       m = floor_log2 (max_digits);
1494
1495       /* Iterate over the bits of the possible powers of 10 that might
1496          be present in U and eliminate them.  That is, if we find that
1497          10**2**M divides U evenly, keep the division and increase
1498          DEC_EXP by 2**M.  */
1499       do
1500         {
1501           REAL_VALUE_TYPE t;
1502
1503           do_divide (&t, &u, ten_to_ptwo (m));
1504           do_fix_trunc (&v, &t);
1505           if (cmp_significands (&v, &t) == 0)
1506             {
1507               u = t;
1508               dec_exp += 1 << m;
1509             }
1510         }
1511       while (--m >= 0);
1512
1513       /* Revert the scaling to integer that we performed earlier.  */
1514       u.exp += r.exp - (SIGNIFICAND_BITS - 1);
1515       r = u;
1516
1517       /* Find power of 10.  Do this by dividing out 10**2**M when
1518          this is larger than the current remainder.  Fill PTEN with
1519          the power of 10 that we compute.  */
1520       if (r.exp > 0)
1521         {
1522           m = floor_log2 ((int)(r.exp * M_LOG10_2)) + 1;
1523           do
1524             {
1525               const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1526               if (do_compare (&u, ptentwo, 0) >= 0)
1527                 {
1528                   do_divide (&u, &u, ptentwo);
1529                   do_multiply (&pten, &pten, ptentwo);
1530                   dec_exp += 1 << m;
1531                 }
1532             }
1533           while (--m >= 0);
1534         }
1535       else
1536         /* We managed to divide off enough tens in the above reduction
1537            loop that we've now got a negative exponent.  Fall into the
1538            less-than-one code to compute the proper value for PTEN.  */
1539         cmp_one = -1;
1540     }
1541   if (cmp_one < 0)
1542     {
1543       int m;
1544
1545       /* Number is less than one.  Pad significand with leading
1546          decimal zeros.  */
1547
1548       v = r;
1549       while (1)
1550         {
1551           /* Stop if we'd shift bits off the bottom.  */
1552           if (v.sig[0] & 7)
1553             break;
1554
1555           do_multiply (&u, &v, ten);
1556
1557           /* Stop if we're now >= 1.  */
1558           if (u.exp > 0)
1559             break;
1560
1561           v = u;
1562           dec_exp -= 1;
1563         }
1564       r = v;
1565
1566       /* Find power of 10.  Do this by multiplying in P=10**2**M when
1567          the current remainder is smaller than 1/P.  Fill PTEN with the
1568          power of 10 that we compute.  */
1569       m = floor_log2 ((int)(-r.exp * M_LOG10_2)) + 1;
1570       do
1571         {
1572           const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1573           const REAL_VALUE_TYPE *ptenmtwo = ten_to_mptwo (m);
1574
1575           if (do_compare (&v, ptenmtwo, 0) <= 0)
1576             {
1577               do_multiply (&v, &v, ptentwo);
1578               do_multiply (&pten, &pten, ptentwo);
1579               dec_exp -= 1 << m;
1580             }
1581         }
1582       while (--m >= 0);
1583
1584       /* Invert the positive power of 10 that we've collected so far.  */
1585       do_divide (&pten, one, &pten);
1586     }
1587
1588   p = str;
1589   if (sign)
1590     *p++ = '-';
1591   first = p++;
1592
1593   /* At this point, PTEN should contain the nearest power of 10 smaller
1594      than R, such that this division produces the first digit.
1595
1596      Using a divide-step primitive that returns the complete integral
1597      remainder avoids the rounding error that would be produced if
1598      we were to use do_divide here and then simply multiply by 10 for
1599      each subsequent digit.  */
1600
1601   digit = rtd_divmod (&r, &pten);
1602
1603   /* Be prepared for error in that division via underflow ...  */
1604   if (digit == 0 && cmp_significand_0 (&r))
1605     {
1606       /* Multiply by 10 and try again.  */
1607       do_multiply (&r, &r, ten);
1608       digit = rtd_divmod (&r, &pten);
1609       dec_exp -= 1;
1610       if (digit == 0)
1611         abort ();
1612     }
1613
1614   /* ... or overflow.  */
1615   if (digit == 10)
1616     {
1617       *p++ = '1';
1618       if (--digits > 0)
1619         *p++ = '0';
1620       dec_exp += 1;
1621     }
1622   else if (digit > 10)
1623     abort ();
1624   else
1625     *p++ = digit + '0';
1626
1627   /* Generate subsequent digits.  */
1628   while (--digits > 0)
1629     {
1630       do_multiply (&r, &r, ten);
1631       digit = rtd_divmod (&r, &pten);
1632       *p++ = digit + '0';
1633     }
1634   last = p;
1635
1636   /* Generate one more digit with which to do rounding.  */
1637   do_multiply (&r, &r, ten);
1638   digit = rtd_divmod (&r, &pten);
1639
1640   /* Round the result.  */
1641   if (digit == 5)
1642     {
1643       /* Round to nearest.  If R is nonzero there are additional
1644          nonzero digits to be extracted.  */
1645       if (cmp_significand_0 (&r))
1646         digit++;
1647       /* Round to even.  */
1648       else if ((p[-1] - '0') & 1)
1649         digit++;
1650     }
1651   if (digit > 5)
1652     {
1653       while (p > first)
1654         {
1655           digit = *--p;
1656           if (digit == '9')
1657             *p = '0';
1658           else
1659             {
1660               *p = digit + 1;
1661               break;
1662             }
1663         }
1664
1665       /* Carry out of the first digit.  This means we had all 9's and
1666          now have all 0's.  "Prepend" a 1 by overwriting the first 0.  */
1667       if (p == first)
1668         {
1669           first[1] = '1';
1670           dec_exp++;
1671         }
1672     }
1673
1674   /* Insert the decimal point.  */
1675   first[0] = first[1];
1676   first[1] = '.';
1677
1678   /* If requested, drop trailing zeros.  Never crop past "1.0".  */
1679   if (crop_trailing_zeros)
1680     while (last > first + 3 && last[-1] == '0')
1681       last--;
1682
1683   /* Append the exponent.  */
1684   sprintf (last, "e%+d", dec_exp);
1685 }
1686
1687 /* Render R as a hexadecimal floating point constant.  Emit DIGITS
1688    significant digits in the result, bounded by BUF_SIZE.  If DIGITS is 0,
1689    choose the maximum for the representation.  If CROP_TRAILING_ZEROS,
1690    strip trailing zeros.  */
1691
1692 void
1693 real_to_hexadecimal (char *str, const REAL_VALUE_TYPE *r, size_t buf_size,
1694                      size_t digits, int crop_trailing_zeros)
1695 {
1696   int i, j, exp = r->exp;
1697   char *p, *first;
1698   char exp_buf[16];
1699   size_t max_digits;
1700
1701   switch (r->class)
1702     {
1703     case rvc_zero:
1704       exp = 0;
1705       break;
1706     case rvc_normal:
1707       break;
1708     case rvc_inf:
1709       strcpy (str, (r->sign ? "-Inf" : "+Inf"));
1710       return;
1711     case rvc_nan:
1712       /* ??? Print the significand as well, if not canonical?  */
1713       strcpy (str, (r->sign ? "-NaN" : "+NaN"));
1714       return;
1715     default:
1716       abort ();
1717     }
1718
1719   if (digits == 0)
1720     digits = SIGNIFICAND_BITS / 4;
1721
1722   /* Bound the number of digits printed by the size of the output buffer.  */
1723
1724   sprintf (exp_buf, "p%+d", exp);
1725   max_digits = buf_size - strlen (exp_buf) - r->sign - 4 - 1;
1726   if (max_digits > buf_size)
1727     abort ();
1728   if (digits > max_digits)
1729     digits = max_digits;
1730
1731   p = str;
1732   if (r->sign)
1733     *p++ = '-';
1734   *p++ = '0';
1735   *p++ = 'x';
1736   *p++ = '0';
1737   *p++ = '.';
1738   first = p;
1739
1740   for (i = SIGSZ - 1; i >= 0; --i)
1741     for (j = HOST_BITS_PER_LONG - 4; j >= 0; j -= 4)
1742       {
1743         *p++ = "0123456789abcdef"[(r->sig[i] >> j) & 15];
1744         if (--digits == 0)
1745           goto out;
1746       }
1747
1748  out:
1749   if (crop_trailing_zeros)
1750     while (p > first + 1 && p[-1] == '0')
1751       p--;
1752
1753   sprintf (p, "p%+d", exp);
1754 }
1755
1756 /* Initialize R from a decimal or hexadecimal string.  The string is
1757    assumed to have been syntax checked already.  */
1758
1759 void
1760 real_from_string (REAL_VALUE_TYPE *r, const char *str)
1761 {
1762   int exp = 0;
1763   bool sign = false;
1764
1765   get_zero (r, 0);
1766
1767   if (*str == '-')
1768     {
1769       sign = true;
1770       str++;
1771     }
1772   else if (*str == '+')
1773     str++;
1774
1775   if (str[0] == '0' && (str[1] == 'x' || str[1] == 'X'))
1776     {
1777       /* Hexadecimal floating point.  */
1778       int pos = SIGNIFICAND_BITS - 4, d;
1779
1780       str += 2;
1781
1782       while (*str == '0')
1783         str++;
1784       while (1)
1785         {
1786           d = hex_value (*str);
1787           if (d == _hex_bad)
1788             break;
1789           if (pos >= 0)
1790             {
1791               r->sig[pos / HOST_BITS_PER_LONG]
1792                 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1793               pos -= 4;
1794             }
1795           exp += 4;
1796           str++;
1797         }
1798       if (*str == '.')
1799         {
1800           str++;
1801           if (pos == SIGNIFICAND_BITS - 4)
1802             {
1803               while (*str == '0')
1804                 str++, exp -= 4;
1805             }
1806           while (1)
1807             {
1808               d = hex_value (*str);
1809               if (d == _hex_bad)
1810                 break;
1811               if (pos >= 0)
1812                 {
1813                   r->sig[pos / HOST_BITS_PER_LONG]
1814                     |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1815                   pos -= 4;
1816                 }
1817               str++;
1818             }
1819         }
1820       if (*str == 'p' || *str == 'P')
1821         {
1822           bool exp_neg = false;
1823
1824           str++;
1825           if (*str == '-')
1826             {
1827               exp_neg = true;
1828               str++;
1829             }
1830           else if (*str == '+')
1831             str++;
1832
1833           d = 0;
1834           while (ISDIGIT (*str))
1835             {
1836               d *= 10;
1837               d += *str - '0';
1838               if (d > MAX_EXP)
1839                 {
1840                   /* Overflowed the exponent.  */
1841                   if (exp_neg)
1842                     goto underflow;
1843                   else
1844                     goto overflow;
1845                 }
1846               str++;
1847             }
1848           if (exp_neg)
1849             d = -d;
1850
1851           exp += d;
1852         }
1853
1854       r->class = rvc_normal;
1855       r->exp = exp;
1856
1857       normalize (r);
1858     }
1859   else
1860     {
1861       /* Decimal floating point.  */
1862       const REAL_VALUE_TYPE *ten = ten_to_ptwo (0);
1863       int d;
1864
1865       while (*str == '0')
1866         str++;
1867       while (ISDIGIT (*str))
1868         {
1869           d = *str++ - '0';
1870           do_multiply (r, r, ten);
1871           if (d)
1872             do_add (r, r, real_digit (d), 0);
1873         }
1874       if (*str == '.')
1875         {
1876           str++;
1877           if (r->class == rvc_zero)
1878             {
1879               while (*str == '0')
1880                 str++, exp--;
1881             }
1882           while (ISDIGIT (*str))
1883             {
1884               d = *str++ - '0';
1885               do_multiply (r, r, ten);
1886               if (d)
1887                 do_add (r, r, real_digit (d), 0);
1888               exp--;
1889             }
1890         }
1891
1892       if (*str == 'e' || *str == 'E')
1893         {
1894           bool exp_neg = false;
1895
1896           str++;
1897           if (*str == '-')
1898             {
1899               exp_neg = true;
1900               str++;
1901             }
1902           else if (*str == '+')
1903             str++;
1904
1905           d = 0;
1906           while (ISDIGIT (*str))
1907             {
1908               d *= 10;
1909               d += *str - '0';
1910               if (d > MAX_EXP)
1911                 {
1912                   /* Overflowed the exponent.  */
1913                   if (exp_neg)
1914                     goto underflow;
1915                   else
1916                     goto overflow;
1917                 }
1918               str++;
1919             }
1920           if (exp_neg)
1921             d = -d;
1922           exp += d;
1923         }
1924
1925       if (exp)
1926         times_pten (r, exp);
1927     }
1928
1929   r->sign = sign;
1930   return;
1931
1932  underflow:
1933   get_zero (r, sign);
1934   return;
1935
1936  overflow:
1937   get_inf (r, sign);
1938   return;
1939 }
1940
1941 /* Legacy.  Similar, but return the result directly.  */
1942
1943 REAL_VALUE_TYPE
1944 real_from_string2 (const char *s, enum machine_mode mode)
1945 {
1946   REAL_VALUE_TYPE r;
1947
1948   real_from_string (&r, s);
1949   if (mode != VOIDmode)
1950     real_convert (&r, mode, &r);
1951
1952   return r;
1953 }
1954
1955 /* Initialize R from the integer pair HIGH+LOW.  */
1956
1957 void
1958 real_from_integer (REAL_VALUE_TYPE *r, enum machine_mode mode,
1959                    unsigned HOST_WIDE_INT low, HOST_WIDE_INT high,
1960                    int unsigned_p)
1961 {
1962   if (low == 0 && high == 0)
1963     get_zero (r, 0);
1964   else
1965     {
1966       memset (r, 0, sizeof (*r));
1967       r->class = rvc_normal;
1968       r->sign = high < 0 && !unsigned_p;
1969       r->exp = 2 * HOST_BITS_PER_WIDE_INT;
1970
1971       if (r->sign)
1972         {
1973           high = ~high;
1974           if (low == 0)
1975             high += 1;
1976           else
1977             low = -low;
1978         }
1979
1980       if (HOST_BITS_PER_LONG == HOST_BITS_PER_WIDE_INT)
1981         {
1982           r->sig[SIGSZ-1] = high;
1983           r->sig[SIGSZ-2] = low;
1984         }
1985       else if (HOST_BITS_PER_LONG*2 == HOST_BITS_PER_WIDE_INT)
1986         {
1987           r->sig[SIGSZ-1] = high >> (HOST_BITS_PER_LONG - 1) >> 1;
1988           r->sig[SIGSZ-2] = high;
1989           r->sig[SIGSZ-3] = low >> (HOST_BITS_PER_LONG - 1) >> 1;
1990           r->sig[SIGSZ-4] = low;
1991         }
1992       else
1993         abort ();
1994
1995       normalize (r);
1996     }
1997
1998   if (mode != VOIDmode)
1999     real_convert (r, mode, r);
2000 }
2001
2002 /* Returns 10**2**N.  */
2003
2004 static const REAL_VALUE_TYPE *
2005 ten_to_ptwo (int n)
2006 {
2007   static REAL_VALUE_TYPE tens[EXP_BITS];
2008
2009   if (n < 0 || n >= EXP_BITS)
2010     abort ();
2011
2012   if (tens[n].class == rvc_zero)
2013     {
2014       if (n < (HOST_BITS_PER_WIDE_INT == 64 ? 5 : 4))
2015         {
2016           HOST_WIDE_INT t = 10;
2017           int i;
2018
2019           for (i = 0; i < n; ++i)
2020             t *= t;
2021
2022           real_from_integer (&tens[n], VOIDmode, t, 0, 1);
2023         }
2024       else
2025         {
2026           const REAL_VALUE_TYPE *t = ten_to_ptwo (n - 1);
2027           do_multiply (&tens[n], t, t);
2028         }
2029     }
2030
2031   return &tens[n];
2032 }
2033
2034 /* Returns 10**(-2**N).  */
2035
2036 static const REAL_VALUE_TYPE *
2037 ten_to_mptwo (int n)
2038 {
2039   static REAL_VALUE_TYPE tens[EXP_BITS];
2040
2041   if (n < 0 || n >= EXP_BITS)
2042     abort ();
2043
2044   if (tens[n].class == rvc_zero)
2045     do_divide (&tens[n], real_digit (1), ten_to_ptwo (n));
2046
2047   return &tens[n];
2048 }
2049
2050 /* Returns N.  */
2051
2052 static const REAL_VALUE_TYPE *
2053 real_digit (int n)
2054 {
2055   static REAL_VALUE_TYPE num[10];
2056
2057   if (n < 0 || n > 9)
2058     abort ();
2059
2060   if (n > 0 && num[n].class == rvc_zero)
2061     real_from_integer (&num[n], VOIDmode, n, 0, 1);
2062
2063   return &num[n];
2064 }
2065
2066 /* Multiply R by 10**EXP.  */
2067
2068 static void
2069 times_pten (REAL_VALUE_TYPE *r, int exp)
2070 {
2071   REAL_VALUE_TYPE pten, *rr;
2072   bool negative = (exp < 0);
2073   int i;
2074
2075   if (negative)
2076     {
2077       exp = -exp;
2078       pten = *real_digit (1);
2079       rr = &pten;
2080     }
2081   else
2082     rr = r;
2083
2084   for (i = 0; exp > 0; ++i, exp >>= 1)
2085     if (exp & 1)
2086       do_multiply (rr, rr, ten_to_ptwo (i));
2087
2088   if (negative)
2089     do_divide (r, r, &pten);
2090 }
2091
2092 /* Fills R with +Inf.  */
2093
2094 void
2095 real_inf (REAL_VALUE_TYPE *r)
2096 {
2097   get_inf (r, 0);
2098 }
2099
2100 /* Fills R with a NaN whose significand is described by STR.  If QUIET,
2101    we force a QNaN, else we force an SNaN.  The string, if not empty,
2102    is parsed as a number and placed in the significand.  Return true
2103    if the string was successfully parsed.  */
2104
2105 bool
2106 real_nan (REAL_VALUE_TYPE *r, const char *str, int quiet,
2107           enum machine_mode mode)
2108 {
2109   const struct real_format *fmt;
2110
2111   fmt = REAL_MODE_FORMAT (mode);
2112   if (fmt == NULL)
2113     abort ();
2114
2115   if (*str == 0)
2116     {
2117       if (quiet)
2118         get_canonical_qnan (r, 0);
2119       else
2120         get_canonical_snan (r, 0);
2121     }
2122   else
2123     {
2124       int base = 10, d;
2125       bool neg = false;
2126
2127       memset (r, 0, sizeof (*r));
2128       r->class = rvc_nan;
2129
2130       /* Parse akin to strtol into the significand of R.  */
2131
2132       while (ISSPACE (*str))
2133         str++;
2134       if (*str == '-')
2135         str++, neg = true;
2136       else if (*str == '+')
2137         str++;
2138       if (*str == '0')
2139         {
2140           if (*++str == 'x')
2141             str++, base = 16;
2142           else
2143             base = 8;
2144         }
2145
2146       while ((d = hex_value (*str)) < base)
2147         {
2148           REAL_VALUE_TYPE u;
2149
2150           switch (base)
2151             {
2152             case 8:
2153               lshift_significand (r, r, 3);
2154               break;
2155             case 16:
2156               lshift_significand (r, r, 4);
2157               break;
2158             case 10:
2159               lshift_significand_1 (&u, r);
2160               lshift_significand (r, r, 3);
2161               add_significands (r, r, &u);
2162               break;
2163             default:
2164               abort ();
2165             }
2166
2167           get_zero (&u, 0);
2168           u.sig[0] = d;
2169           add_significands (r, r, &u);
2170
2171           str++;
2172         }
2173
2174       /* Must have consumed the entire string for success.  */
2175       if (*str != 0)
2176         return false;
2177
2178       /* Shift the significand into place such that the bits
2179          are in the most significant bits for the format.  */
2180       lshift_significand (r, r, SIGNIFICAND_BITS - fmt->pnan);
2181
2182       /* Our MSB is always unset for NaNs.  */
2183       r->sig[SIGSZ-1] &= ~SIG_MSB;
2184
2185       /* Force quiet or signalling NaN.  */
2186       r->signalling = !quiet;
2187     }
2188
2189   return true;
2190 }
2191
2192 /* Fills R with the largest finite value representable in mode MODE.
2193    If SIGN is nonzero, R is set to the most negative finite value.  */
2194
2195 void
2196 real_maxval (REAL_VALUE_TYPE *r, int sign, enum machine_mode mode)
2197 {
2198   const struct real_format *fmt;
2199   int np2;
2200
2201   fmt = REAL_MODE_FORMAT (mode);
2202   if (fmt == NULL)
2203     abort ();
2204
2205   r->class = rvc_normal;
2206   r->sign = sign;
2207   r->signalling = 0;
2208   r->canonical = 0;
2209   r->exp = fmt->emax * fmt->log2_b;
2210
2211   np2 = SIGNIFICAND_BITS - fmt->p * fmt->log2_b;
2212   memset (r->sig, -1, SIGSZ * sizeof (unsigned long));
2213   clear_significand_below (r, np2);
2214 }
2215
2216 /* Fills R with 2**N.  */
2217
2218 void
2219 real_2expN (REAL_VALUE_TYPE *r, int n)
2220 {
2221   memset (r, 0, sizeof (*r));
2222
2223   n++;
2224   if (n > MAX_EXP)
2225     r->class = rvc_inf;
2226   else if (n < -MAX_EXP)
2227     ;
2228   else
2229     {
2230       r->class = rvc_normal;
2231       r->exp = n;
2232       r->sig[SIGSZ-1] = SIG_MSB;
2233     }
2234 }
2235
2236 \f
2237 static void
2238 round_for_format (const struct real_format *fmt, REAL_VALUE_TYPE *r)
2239 {
2240   int p2, np2, i, w;
2241   unsigned long sticky;
2242   bool guard, lsb;
2243   int emin2m1, emax2;
2244
2245   p2 = fmt->p * fmt->log2_b;
2246   emin2m1 = (fmt->emin - 1) * fmt->log2_b;
2247   emax2 = fmt->emax * fmt->log2_b;
2248
2249   np2 = SIGNIFICAND_BITS - p2;
2250   switch (r->class)
2251     {
2252     underflow:
2253       get_zero (r, r->sign);
2254     case rvc_zero:
2255       if (!fmt->has_signed_zero)
2256         r->sign = 0;
2257       return;
2258
2259     overflow:
2260       get_inf (r, r->sign);
2261     case rvc_inf:
2262       return;
2263
2264     case rvc_nan:
2265       clear_significand_below (r, np2);
2266       return;
2267
2268     case rvc_normal:
2269       break;
2270
2271     default:
2272       abort ();
2273     }
2274
2275   /* If we're not base2, normalize the exponent to a multiple of
2276      the true base.  */
2277   if (fmt->log2_b != 1)
2278     {
2279       int shift = r->exp & (fmt->log2_b - 1);
2280       if (shift)
2281         {
2282           shift = fmt->log2_b - shift;
2283           r->sig[0] |= sticky_rshift_significand (r, r, shift);
2284           r->exp += shift;
2285         }
2286     }
2287
2288   /* Check the range of the exponent.  If we're out of range,
2289      either underflow or overflow.  */
2290   if (r->exp > emax2)
2291     goto overflow;
2292   else if (r->exp <= emin2m1)
2293     {
2294       int diff;
2295
2296       if (!fmt->has_denorm)
2297         {
2298           /* Don't underflow completely until we've had a chance to round.  */
2299           if (r->exp < emin2m1)
2300             goto underflow;
2301         }
2302       else
2303         {
2304           diff = emin2m1 - r->exp + 1;
2305           if (diff > p2)
2306             goto underflow;
2307
2308           /* De-normalize the significand.  */
2309           r->sig[0] |= sticky_rshift_significand (r, r, diff);
2310           r->exp += diff;
2311         }
2312     }
2313
2314   /* There are P2 true significand bits, followed by one guard bit,
2315      followed by one sticky bit, followed by stuff.  Fold nonzero
2316      stuff into the sticky bit.  */
2317
2318   sticky = 0;
2319   for (i = 0, w = (np2 - 1) / HOST_BITS_PER_LONG; i < w; ++i)
2320     sticky |= r->sig[i];
2321   sticky |=
2322     r->sig[w] & (((unsigned long)1 << ((np2 - 1) % HOST_BITS_PER_LONG)) - 1);
2323
2324   guard = test_significand_bit (r, np2 - 1);
2325   lsb = test_significand_bit (r, np2);
2326
2327   /* Round to even.  */
2328   if (guard && (sticky || lsb))
2329     {
2330       REAL_VALUE_TYPE u;
2331       get_zero (&u, 0);
2332       set_significand_bit (&u, np2);
2333
2334       if (add_significands (r, r, &u))
2335         {
2336           /* Overflow.  Means the significand had been all ones, and
2337              is now all zeros.  Need to increase the exponent, and
2338              possibly re-normalize it.  */
2339           if (++r->exp > emax2)
2340             goto overflow;
2341           r->sig[SIGSZ-1] = SIG_MSB;
2342
2343           if (fmt->log2_b != 1)
2344             {
2345               int shift = r->exp & (fmt->log2_b - 1);
2346               if (shift)
2347                 {
2348                   shift = fmt->log2_b - shift;
2349                   rshift_significand (r, r, shift);
2350                   r->exp += shift;
2351                   if (r->exp > emax2)
2352                     goto overflow;
2353                 }
2354             }
2355         }
2356     }
2357
2358   /* Catch underflow that we deferred until after rounding.  */
2359   if (r->exp <= emin2m1)
2360     goto underflow;
2361
2362   /* Clear out trailing garbage.  */
2363   clear_significand_below (r, np2);
2364 }
2365
2366 /* Extend or truncate to a new mode.  */
2367
2368 void
2369 real_convert (REAL_VALUE_TYPE *r, enum machine_mode mode,
2370               const REAL_VALUE_TYPE *a)
2371 {
2372   const struct real_format *fmt;
2373
2374   fmt = REAL_MODE_FORMAT (mode);
2375   if (fmt == NULL)
2376     abort ();
2377
2378   *r = *a;
2379   round_for_format (fmt, r);
2380
2381   /* round_for_format de-normalizes denormals.  Undo just that part.  */
2382   if (r->class == rvc_normal)
2383     normalize (r);
2384 }
2385
2386 /* Legacy.  Likewise, except return the struct directly.  */
2387
2388 REAL_VALUE_TYPE
2389 real_value_truncate (enum machine_mode mode, REAL_VALUE_TYPE a)
2390 {
2391   REAL_VALUE_TYPE r;
2392   real_convert (&r, mode, &a);
2393   return r;
2394 }
2395
2396 /* Return true if truncating to MODE is exact.  */
2397
2398 bool
2399 exact_real_truncate (enum machine_mode mode, const REAL_VALUE_TYPE *a)
2400 {
2401   REAL_VALUE_TYPE t;
2402   real_convert (&t, mode, a);
2403   return real_identical (&t, a);
2404 }
2405
2406 /* Write R to the given target format.  Place the words of the result
2407    in target word order in BUF.  There are always 32 bits in each
2408    long, no matter the size of the host long.
2409
2410    Legacy: return word 0 for implementing REAL_VALUE_TO_TARGET_SINGLE.  */
2411
2412 long
2413 real_to_target_fmt (long *buf, const REAL_VALUE_TYPE *r_orig,
2414                     const struct real_format *fmt)
2415 {
2416   REAL_VALUE_TYPE r;
2417   long buf1;
2418
2419   r = *r_orig;
2420   round_for_format (fmt, &r);
2421
2422   if (!buf)
2423     buf = &buf1;
2424   (*fmt->encode) (fmt, buf, &r);
2425
2426   return *buf;
2427 }
2428
2429 /* Similar, but look up the format from MODE.  */
2430
2431 long
2432 real_to_target (long *buf, const REAL_VALUE_TYPE *r, enum machine_mode mode)
2433 {
2434   const struct real_format *fmt;
2435
2436   fmt = REAL_MODE_FORMAT (mode);
2437   if (fmt == NULL)
2438     abort ();
2439
2440   return real_to_target_fmt (buf, r, fmt);
2441 }
2442
2443 /* Read R from the given target format.  Read the words of the result
2444    in target word order in BUF.  There are always 32 bits in each
2445    long, no matter the size of the host long.  */
2446
2447 void
2448 real_from_target_fmt (REAL_VALUE_TYPE *r, const long *buf,
2449                       const struct real_format *fmt)
2450 {
2451   (*fmt->decode) (fmt, r, buf);
2452 }
2453
2454 /* Similar, but look up the format from MODE.  */
2455
2456 void
2457 real_from_target (REAL_VALUE_TYPE *r, const long *buf, enum machine_mode mode)
2458 {
2459   const struct real_format *fmt;
2460
2461   fmt = REAL_MODE_FORMAT (mode);
2462   if (fmt == NULL)
2463     abort ();
2464
2465   (*fmt->decode) (fmt, r, buf);
2466 }
2467
2468 /* Return the number of bits in the significand for MODE.  */
2469 /* ??? Legacy.  Should get access to real_format directly.  */
2470
2471 int
2472 significand_size (enum machine_mode mode)
2473 {
2474   const struct real_format *fmt;
2475
2476   fmt = REAL_MODE_FORMAT (mode);
2477   if (fmt == NULL)
2478     return 0;
2479
2480   return fmt->p * fmt->log2_b;
2481 }
2482
2483 /* Return a hash value for the given real value.  */
2484 /* ??? The "unsigned int" return value is intended to be hashval_t,
2485    but I didn't want to pull hashtab.h into real.h.  */
2486
2487 unsigned int
2488 real_hash (const REAL_VALUE_TYPE *r)
2489 {
2490   unsigned int h;
2491   size_t i;
2492
2493   h = r->class | (r->sign << 2);
2494   switch (r->class)
2495     {
2496     case rvc_zero:
2497     case rvc_inf:
2498       return h;
2499
2500     case rvc_normal:
2501       h |= r->exp << 3;
2502       break;
2503
2504     case rvc_nan:
2505       if (r->signalling)
2506         h ^= (unsigned int)-1;
2507       if (r->canonical)
2508         return h;
2509       break;
2510
2511     default:
2512       abort ();
2513     }
2514
2515   if (sizeof(unsigned long) > sizeof(unsigned int))
2516     for (i = 0; i < SIGSZ; ++i)
2517       {
2518         unsigned long s = r->sig[i];
2519         h ^= s ^ (s >> (HOST_BITS_PER_LONG / 2));
2520       }
2521   else
2522     for (i = 0; i < SIGSZ; ++i)
2523       h ^= r->sig[i];
2524
2525   return h;
2526 }
2527 \f
2528 /* IEEE single-precision format.  */
2529
2530 static void encode_ieee_single (const struct real_format *fmt,
2531                                 long *, const REAL_VALUE_TYPE *);
2532 static void decode_ieee_single (const struct real_format *,
2533                                 REAL_VALUE_TYPE *, const long *);
2534
2535 static void
2536 encode_ieee_single (const struct real_format *fmt, long *buf,
2537                     const REAL_VALUE_TYPE *r)
2538 {
2539   unsigned long image, sig, exp;
2540   unsigned long sign = r->sign;
2541   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2542
2543   image = sign << 31;
2544   sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
2545
2546   switch (r->class)
2547     {
2548     case rvc_zero:
2549       break;
2550
2551     case rvc_inf:
2552       if (fmt->has_inf)
2553         image |= 255 << 23;
2554       else
2555         image |= 0x7fffffff;
2556       break;
2557
2558     case rvc_nan:
2559       if (fmt->has_nans)
2560         {
2561           if (r->canonical)
2562             sig = 0;
2563           if (r->signalling == fmt->qnan_msb_set)
2564             sig &= ~(1 << 22);
2565           else
2566             sig |= 1 << 22;
2567           /* We overload qnan_msb_set here: it's only clear for
2568              mips_ieee_single, which wants all mantissa bits but the
2569              quiet/signalling one set in canonical NaNs (at least
2570              Quiet ones).  */
2571           if (r->canonical && !fmt->qnan_msb_set)
2572             sig |= (1 << 22) - 1;
2573           else if (sig == 0)
2574             sig = 1 << 21;
2575
2576           image |= 255 << 23;
2577           image |= sig;
2578         }
2579       else
2580         image |= 0x7fffffff;
2581       break;
2582
2583     case rvc_normal:
2584       /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2585          whereas the intermediate representation is 0.F x 2**exp.
2586          Which means we're off by one.  */
2587       if (denormal)
2588         exp = 0;
2589       else
2590       exp = r->exp + 127 - 1;
2591       image |= exp << 23;
2592       image |= sig;
2593       break;
2594
2595     default:
2596       abort ();
2597     }
2598
2599   buf[0] = image;
2600 }
2601
2602 static void
2603 decode_ieee_single (const struct real_format *fmt, REAL_VALUE_TYPE *r,
2604                     const long *buf)
2605 {
2606   unsigned long image = buf[0] & 0xffffffff;
2607   bool sign = (image >> 31) & 1;
2608   int exp = (image >> 23) & 0xff;
2609
2610   memset (r, 0, sizeof (*r));
2611   image <<= HOST_BITS_PER_LONG - 24;
2612   image &= ~SIG_MSB;
2613
2614   if (exp == 0)
2615     {
2616       if (image && fmt->has_denorm)
2617         {
2618           r->class = rvc_normal;
2619           r->sign = sign;
2620           r->exp = -126;
2621           r->sig[SIGSZ-1] = image << 1;
2622           normalize (r);
2623         }
2624       else if (fmt->has_signed_zero)
2625         r->sign = sign;
2626     }
2627   else if (exp == 255 && (fmt->has_nans || fmt->has_inf))
2628     {
2629       if (image)
2630         {
2631           r->class = rvc_nan;
2632           r->sign = sign;
2633           r->signalling = (((image >> (HOST_BITS_PER_LONG - 2)) & 1)
2634                            ^ fmt->qnan_msb_set);
2635           r->sig[SIGSZ-1] = image;
2636         }
2637       else
2638         {
2639           r->class = rvc_inf;
2640           r->sign = sign;
2641         }
2642     }
2643   else
2644     {
2645       r->class = rvc_normal;
2646       r->sign = sign;
2647       r->exp = exp - 127 + 1;
2648       r->sig[SIGSZ-1] = image | SIG_MSB;
2649     }
2650 }
2651
2652 const struct real_format ieee_single_format =
2653   {
2654     encode_ieee_single,
2655     decode_ieee_single,
2656     2,
2657     1,
2658     24,
2659     24,
2660     -125,
2661     128,
2662     31,
2663     true,
2664     true,
2665     true,
2666     true,
2667     true
2668   };
2669
2670 const struct real_format mips_single_format =
2671   {
2672     encode_ieee_single,
2673     decode_ieee_single,
2674     2,
2675     1,
2676     24,
2677     24,
2678     -125,
2679     128,
2680     31,
2681     true,
2682     true,
2683     true,
2684     true,
2685     false
2686   };
2687
2688 \f
2689 /* IEEE double-precision format.  */
2690
2691 static void encode_ieee_double (const struct real_format *fmt,
2692                                 long *, const REAL_VALUE_TYPE *);
2693 static void decode_ieee_double (const struct real_format *,
2694                                 REAL_VALUE_TYPE *, const long *);
2695
2696 static void
2697 encode_ieee_double (const struct real_format *fmt, long *buf,
2698                     const REAL_VALUE_TYPE *r)
2699 {
2700   unsigned long image_lo, image_hi, sig_lo, sig_hi, exp;
2701   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2702
2703   image_hi = r->sign << 31;
2704   image_lo = 0;
2705
2706   if (HOST_BITS_PER_LONG == 64)
2707     {
2708       sig_hi = r->sig[SIGSZ-1];
2709       sig_lo = (sig_hi >> (64 - 53)) & 0xffffffff;
2710       sig_hi = (sig_hi >> (64 - 53 + 1) >> 31) & 0xfffff;
2711     }
2712   else
2713     {
2714       sig_hi = r->sig[SIGSZ-1];
2715       sig_lo = r->sig[SIGSZ-2];
2716       sig_lo = (sig_hi << 21) | (sig_lo >> 11);
2717       sig_hi = (sig_hi >> 11) & 0xfffff;
2718     }
2719
2720   switch (r->class)
2721     {
2722     case rvc_zero:
2723       break;
2724
2725     case rvc_inf:
2726       if (fmt->has_inf)
2727         image_hi |= 2047 << 20;
2728       else
2729         {
2730           image_hi |= 0x7fffffff;
2731           image_lo = 0xffffffff;
2732         }
2733       break;
2734
2735     case rvc_nan:
2736       if (fmt->has_nans)
2737         {
2738           if (r->canonical)
2739             sig_hi = sig_lo = 0;
2740           if (r->signalling == fmt->qnan_msb_set)
2741             sig_hi &= ~(1 << 19);
2742           else
2743             sig_hi |= 1 << 19;
2744           /* We overload qnan_msb_set here: it's only clear for
2745              mips_ieee_single, which wants all mantissa bits but the
2746              quiet/signalling one set in canonical NaNs (at least
2747              Quiet ones).  */
2748           if (r->canonical && !fmt->qnan_msb_set)
2749             {
2750               sig_hi |= (1 << 19) - 1;
2751               sig_lo = 0xffffffff;
2752             }
2753           else if (sig_hi == 0 && sig_lo == 0)
2754             sig_hi = 1 << 18;
2755
2756           image_hi |= 2047 << 20;
2757           image_hi |= sig_hi;
2758           image_lo = sig_lo;
2759         }
2760       else
2761         {
2762           image_hi |= 0x7fffffff;
2763           image_lo = 0xffffffff;
2764         }
2765       break;
2766
2767     case rvc_normal:
2768       /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2769          whereas the intermediate representation is 0.F x 2**exp.
2770          Which means we're off by one.  */
2771       if (denormal)
2772         exp = 0;
2773       else
2774         exp = r->exp + 1023 - 1;
2775       image_hi |= exp << 20;
2776       image_hi |= sig_hi;
2777       image_lo = sig_lo;
2778       break;
2779
2780     default:
2781       abort ();
2782     }
2783
2784   if (FLOAT_WORDS_BIG_ENDIAN)
2785     buf[0] = image_hi, buf[1] = image_lo;
2786   else
2787     buf[0] = image_lo, buf[1] = image_hi;
2788 }
2789
2790 static void
2791 decode_ieee_double (const struct real_format *fmt, REAL_VALUE_TYPE *r,
2792                     const long *buf)
2793 {
2794   unsigned long image_hi, image_lo;
2795   bool sign;
2796   int exp;
2797
2798   if (FLOAT_WORDS_BIG_ENDIAN)
2799     image_hi = buf[0], image_lo = buf[1];
2800   else
2801     image_lo = buf[0], image_hi = buf[1];
2802   image_lo &= 0xffffffff;
2803   image_hi &= 0xffffffff;
2804
2805   sign = (image_hi >> 31) & 1;
2806   exp = (image_hi >> 20) & 0x7ff;
2807
2808   memset (r, 0, sizeof (*r));
2809
2810   image_hi <<= 32 - 21;
2811   image_hi |= image_lo >> 21;
2812   image_hi &= 0x7fffffff;
2813   image_lo <<= 32 - 21;
2814
2815   if (exp == 0)
2816     {
2817       if ((image_hi || image_lo) && fmt->has_denorm)
2818         {
2819           r->class = rvc_normal;
2820           r->sign = sign;
2821           r->exp = -1022;
2822           if (HOST_BITS_PER_LONG == 32)
2823             {
2824               image_hi = (image_hi << 1) | (image_lo >> 31);
2825               image_lo <<= 1;
2826               r->sig[SIGSZ-1] = image_hi;
2827               r->sig[SIGSZ-2] = image_lo;
2828             }
2829           else
2830             {
2831               image_hi = (image_hi << 31 << 2) | (image_lo << 1);
2832               r->sig[SIGSZ-1] = image_hi;
2833             }
2834           normalize (r);
2835         }
2836       else if (fmt->has_signed_zero)
2837         r->sign = sign;
2838     }
2839   else if (exp == 2047 && (fmt->has_nans || fmt->has_inf))
2840     {
2841       if (image_hi || image_lo)
2842         {
2843           r->class = rvc_nan;
2844           r->sign = sign;
2845           r->signalling = ((image_hi >> 30) & 1) ^ fmt->qnan_msb_set;
2846           if (HOST_BITS_PER_LONG == 32)
2847             {
2848               r->sig[SIGSZ-1] = image_hi;
2849               r->sig[SIGSZ-2] = image_lo;
2850             }
2851           else
2852             r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo;
2853         }
2854       else
2855         {
2856           r->class = rvc_inf;
2857           r->sign = sign;
2858         }
2859     }
2860   else
2861     {
2862       r->class = rvc_normal;
2863       r->sign = sign;
2864       r->exp = exp - 1023 + 1;
2865       if (HOST_BITS_PER_LONG == 32)
2866         {
2867           r->sig[SIGSZ-1] = image_hi | SIG_MSB;
2868           r->sig[SIGSZ-2] = image_lo;
2869         }
2870       else
2871         r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo | SIG_MSB;
2872     }
2873 }
2874
2875 const struct real_format ieee_double_format =
2876   {
2877     encode_ieee_double,
2878     decode_ieee_double,
2879     2,
2880     1,
2881     53,
2882     53,
2883     -1021,
2884     1024,
2885     63,
2886     true,
2887     true,
2888     true,
2889     true,
2890     true
2891   };
2892
2893 const struct real_format mips_double_format =
2894   {
2895     encode_ieee_double,
2896     decode_ieee_double,
2897     2,
2898     1,
2899     53,
2900     53,
2901     -1021,
2902     1024,
2903     63,
2904     true,
2905     true,
2906     true,
2907     true,
2908     false
2909   };
2910
2911 \f
2912 /* IEEE extended real format.  This comes in three flavors: Intel's as
2913    a 12 byte image, Intel's as a 16 byte image, and Motorola's.  Intel
2914    12- and 16-byte images may be big- or little endian; Motorola's is
2915    always big endian.  */
2916
2917 /* Helper subroutine which converts from the internal format to the
2918    12-byte little-endian Intel format.  Functions below adjust this
2919    for the other possible formats.  */
2920 static void
2921 encode_ieee_extended (const struct real_format *fmt, long *buf,
2922                       const REAL_VALUE_TYPE *r)
2923 {
2924   unsigned long image_hi, sig_hi, sig_lo;
2925   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2926
2927   image_hi = r->sign << 15;
2928   sig_hi = sig_lo = 0;
2929
2930   switch (r->class)
2931     {
2932     case rvc_zero:
2933       break;
2934
2935     case rvc_inf:
2936       if (fmt->has_inf)
2937         {
2938           image_hi |= 32767;
2939
2940           /* Intel requires the explicit integer bit to be set, otherwise
2941              it considers the value a "pseudo-infinity".  Motorola docs
2942              say it doesn't care.  */
2943           sig_hi = 0x80000000;
2944         }
2945       else
2946         {
2947           image_hi |= 32767;
2948           sig_lo = sig_hi = 0xffffffff;
2949         }
2950       break;
2951
2952     case rvc_nan:
2953       if (fmt->has_nans)
2954         {
2955           image_hi |= 32767;
2956           if (HOST_BITS_PER_LONG == 32)
2957             {
2958               sig_hi = r->sig[SIGSZ-1];
2959               sig_lo = r->sig[SIGSZ-2];
2960             }
2961           else
2962             {
2963               sig_lo = r->sig[SIGSZ-1];
2964               sig_hi = sig_lo >> 31 >> 1;
2965               sig_lo &= 0xffffffff;
2966             }
2967           if (r->signalling == fmt->qnan_msb_set)
2968             sig_hi &= ~(1 << 30);
2969           else
2970             sig_hi |= 1 << 30;
2971           if ((sig_hi & 0x7fffffff) == 0 && sig_lo == 0)
2972             sig_hi = 1 << 29;
2973
2974           /* Intel requires the explicit integer bit to be set, otherwise
2975              it considers the value a "pseudo-nan".  Motorola docs say it
2976              doesn't care.  */
2977           sig_hi |= 0x80000000;
2978         }
2979       else
2980         {
2981           image_hi |= 32767;
2982           sig_lo = sig_hi = 0xffffffff;
2983         }
2984       break;
2985
2986     case rvc_normal:
2987       {
2988         int exp = r->exp;
2989
2990         /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2991            whereas the intermediate representation is 0.F x 2**exp.
2992            Which means we're off by one.
2993
2994            Except for Motorola, which consider exp=0 and explicit
2995            integer bit set to continue to be normalized.  In theory
2996            this discrepancy has been taken care of by the difference
2997            in fmt->emin in round_for_format.  */
2998
2999         if (denormal)
3000           exp = 0;
3001         else
3002           {
3003             exp += 16383 - 1;
3004             if (exp < 0)
3005               abort ();
3006           }
3007         image_hi |= exp;
3008
3009         if (HOST_BITS_PER_LONG == 32)
3010           {
3011             sig_hi = r->sig[SIGSZ-1];
3012             sig_lo = r->sig[SIGSZ-2];
3013           }
3014         else
3015           {
3016             sig_lo = r->sig[SIGSZ-1];
3017             sig_hi = sig_lo >> 31 >> 1;
3018             sig_lo &= 0xffffffff;
3019           }
3020       }
3021       break;
3022
3023     default:
3024       abort ();
3025     }
3026
3027   buf[0] = sig_lo, buf[1] = sig_hi, buf[2] = image_hi;
3028 }
3029
3030 /* Convert from the internal format to the 12-byte Motorola format
3031    for an IEEE extended real.  */
3032 static void
3033 encode_ieee_extended_motorola (const struct real_format *fmt, long *buf,
3034                                const REAL_VALUE_TYPE *r)
3035 {
3036   long intermed[3];
3037   encode_ieee_extended (fmt, intermed, r);
3038
3039   /* Motorola chips are assumed always to be big-endian.  Also, the
3040      padding in a Motorola extended real goes between the exponent and
3041      the mantissa.  At this point the mantissa is entirely within
3042      elements 0 and 1 of intermed, and the exponent entirely within
3043      element 2, so all we have to do is swap the order around, and
3044      shift element 2 left 16 bits.  */
3045   buf[0] = intermed[2] << 16;
3046   buf[1] = intermed[1];
3047   buf[2] = intermed[0];
3048 }
3049
3050 /* Convert from the internal format to the 12-byte Intel format for
3051    an IEEE extended real.  */
3052 static void
3053 encode_ieee_extended_intel_96 (const struct real_format *fmt, long *buf,
3054                                const REAL_VALUE_TYPE *r)
3055 {
3056   if (FLOAT_WORDS_BIG_ENDIAN)
3057     {
3058       /* All the padding in an Intel-format extended real goes at the high
3059          end, which in this case is after the mantissa, not the exponent.
3060          Therefore we must shift everything down 16 bits.  */
3061       long intermed[3];
3062       encode_ieee_extended (fmt, intermed, r);
3063       buf[0] = ((intermed[2] << 16) | ((unsigned long)(intermed[1] & 0xFFFF0000) >> 16));
3064       buf[1] = ((intermed[1] << 16) | ((unsigned long)(intermed[0] & 0xFFFF0000) >> 16));
3065       buf[2] =  (intermed[0] << 16);
3066     }
3067   else
3068     /* encode_ieee_extended produces what we want directly.  */
3069     encode_ieee_extended (fmt, buf, r);
3070 }
3071
3072 /* Convert from the internal format to the 16-byte Intel format for
3073    an IEEE extended real.  */
3074 static void
3075 encode_ieee_extended_intel_128 (const struct real_format *fmt, long *buf,
3076                                 const REAL_VALUE_TYPE *r)
3077 {
3078   /* All the padding in an Intel-format extended real goes at the high end.  */
3079   encode_ieee_extended_intel_96 (fmt, buf, r);
3080   buf[3] = 0;
3081 }
3082
3083 /* As above, we have a helper function which converts from 12-byte
3084    little-endian Intel format to internal format.  Functions below
3085    adjust for the other possible formats.  */
3086 static void
3087 decode_ieee_extended (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3088                       const long *buf)
3089 {
3090   unsigned long image_hi, sig_hi, sig_lo;
3091   bool sign;
3092   int exp;
3093
3094   sig_lo = buf[0], sig_hi = buf[1], image_hi = buf[2];
3095   sig_lo &= 0xffffffff;
3096   sig_hi &= 0xffffffff;
3097   image_hi &= 0xffffffff;
3098
3099   sign = (image_hi >> 15) & 1;
3100   exp = image_hi & 0x7fff;
3101
3102   memset (r, 0, sizeof (*r));
3103
3104   if (exp == 0)
3105     {
3106       if ((sig_hi || sig_lo) && fmt->has_denorm)
3107         {
3108           r->class = rvc_normal;
3109           r->sign = sign;
3110
3111           /* When the IEEE format contains a hidden bit, we know that
3112              it's zero at this point, and so shift up the significand
3113              and decrease the exponent to match.  In this case, Motorola
3114              defines the explicit integer bit to be valid, so we don't
3115              know whether the msb is set or not.  */
3116           r->exp = fmt->emin;
3117           if (HOST_BITS_PER_LONG == 32)
3118             {
3119               r->sig[SIGSZ-1] = sig_hi;
3120               r->sig[SIGSZ-2] = sig_lo;
3121             }
3122           else
3123             r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3124
3125           normalize (r);
3126         }
3127       else if (fmt->has_signed_zero)
3128         r->sign = sign;
3129     }
3130   else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3131     {
3132       /* See above re "pseudo-infinities" and "pseudo-nans".
3133          Short summary is that the MSB will likely always be
3134          set, and that we don't care about it.  */
3135       sig_hi &= 0x7fffffff;
3136
3137       if (sig_hi || sig_lo)
3138         {
3139           r->class = rvc_nan;
3140           r->sign = sign;
3141           r->signalling = ((sig_hi >> 30) & 1) ^ fmt->qnan_msb_set;
3142           if (HOST_BITS_PER_LONG == 32)
3143             {
3144               r->sig[SIGSZ-1] = sig_hi;
3145               r->sig[SIGSZ-2] = sig_lo;
3146             }
3147           else
3148             r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3149         }
3150       else
3151         {
3152           r->class = rvc_inf;
3153           r->sign = sign;
3154         }
3155     }
3156   else
3157     {
3158       r->class = rvc_normal;
3159       r->sign = sign;
3160       r->exp = exp - 16383 + 1;
3161       if (HOST_BITS_PER_LONG == 32)
3162         {
3163           r->sig[SIGSZ-1] = sig_hi;
3164           r->sig[SIGSZ-2] = sig_lo;
3165         }
3166       else
3167         r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3168     }
3169 }
3170
3171 /* Convert from the internal format to the 12-byte Motorola format
3172    for an IEEE extended real.  */
3173 static void
3174 decode_ieee_extended_motorola (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3175                                const long *buf)
3176 {
3177   long intermed[3];
3178
3179   /* Motorola chips are assumed always to be big-endian.  Also, the
3180      padding in a Motorola extended real goes between the exponent and
3181      the mantissa; remove it.  */
3182   intermed[0] = buf[2];
3183   intermed[1] = buf[1];
3184   intermed[2] = (unsigned long)buf[0] >> 16;
3185
3186   decode_ieee_extended (fmt, r, intermed);
3187 }
3188
3189 /* Convert from the internal format to the 12-byte Intel format for
3190    an IEEE extended real.  */
3191 static void
3192 decode_ieee_extended_intel_96 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3193                                const long *buf)
3194 {
3195   if (FLOAT_WORDS_BIG_ENDIAN)
3196     {
3197       /* All the padding in an Intel-format extended real goes at the high
3198          end, which in this case is after the mantissa, not the exponent.
3199          Therefore we must shift everything up 16 bits.  */
3200       long intermed[3];
3201
3202       intermed[0] = (((unsigned long)buf[2] >> 16) | (buf[1] << 16));
3203       intermed[1] = (((unsigned long)buf[1] >> 16) | (buf[0] << 16));
3204       intermed[2] =  ((unsigned long)buf[0] >> 16);
3205
3206       decode_ieee_extended (fmt, r, intermed);
3207     }
3208   else
3209     /* decode_ieee_extended produces what we want directly.  */
3210     decode_ieee_extended (fmt, r, buf);
3211 }
3212
3213 /* Convert from the internal format to the 16-byte Intel format for
3214    an IEEE extended real.  */
3215 static void
3216 decode_ieee_extended_intel_128 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3217                                 const long *buf)
3218 {
3219   /* All the padding in an Intel-format extended real goes at the high end.  */
3220   decode_ieee_extended_intel_96 (fmt, r, buf);
3221 }
3222
3223 const struct real_format ieee_extended_motorola_format =
3224   {
3225     encode_ieee_extended_motorola,
3226     decode_ieee_extended_motorola,
3227     2,
3228     1,
3229     64,
3230     64,
3231     -16382,
3232     16384,
3233     95,
3234     true,
3235     true,
3236     true,
3237     true,
3238     true
3239   };
3240
3241 const struct real_format ieee_extended_intel_96_format =
3242   {
3243     encode_ieee_extended_intel_96,
3244     decode_ieee_extended_intel_96,
3245     2,
3246     1,
3247     64,
3248     64,
3249     -16381,
3250     16384,
3251     79,
3252     true,
3253     true,
3254     true,
3255     true,
3256     true
3257   };
3258
3259 const struct real_format ieee_extended_intel_128_format =
3260   {
3261     encode_ieee_extended_intel_128,
3262     decode_ieee_extended_intel_128,
3263     2,
3264     1,
3265     64,
3266     64,
3267     -16381,
3268     16384,
3269     79,
3270     true,
3271     true,
3272     true,
3273     true,
3274     true
3275   };
3276
3277 /* The following caters to i386 systems that set the rounding precision
3278    to 53 bits instead of 64, e.g. FreeBSD.  */
3279 const struct real_format ieee_extended_intel_96_round_53_format =
3280   {
3281     encode_ieee_extended_intel_96,
3282     decode_ieee_extended_intel_96,
3283     2,
3284     1,
3285     53,
3286     53,
3287     -16381,
3288     16384,
3289     79,
3290     true,
3291     true,
3292     true,
3293     true,
3294     true
3295   };
3296 \f
3297 /* IBM 128-bit extended precision format: a pair of IEEE double precision
3298    numbers whose sum is equal to the extended precision value.  The number
3299    with greater magnitude is first.  This format has the same magnitude
3300    range as an IEEE double precision value, but effectively 106 bits of
3301    significand precision.  Infinity and NaN are represented by their IEEE
3302    double precision value stored in the first number, the second number is
3303    ignored.  Zeroes, Infinities, and NaNs are set in both doubles
3304    due to precedent.  */
3305
3306 static void encode_ibm_extended (const struct real_format *fmt,
3307                                  long *, const REAL_VALUE_TYPE *);
3308 static void decode_ibm_extended (const struct real_format *,
3309                                  REAL_VALUE_TYPE *, const long *);
3310
3311 static void
3312 encode_ibm_extended (const struct real_format *fmt, long *buf,
3313                      const REAL_VALUE_TYPE *r)
3314 {
3315   REAL_VALUE_TYPE u, normr, v;
3316   const struct real_format *base_fmt;
3317
3318   base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3319
3320   /* Renormlize R before doing any arithmetic on it.  */
3321   normr = *r;
3322   if (normr.class == rvc_normal)
3323     normalize (&normr);
3324
3325   /* u = IEEE double precision portion of significand.  */
3326   u = normr;
3327   round_for_format (base_fmt, &u);
3328   encode_ieee_double (base_fmt, &buf[0], &u);
3329
3330   if (u.class == rvc_normal)
3331     {
3332       do_add (&v, &normr, &u, 1);
3333       /* Call round_for_format since we might need to denormalize.  */
3334       round_for_format (base_fmt, &v);
3335       encode_ieee_double (base_fmt, &buf[2], &v);
3336     }
3337   else
3338     {
3339       /* Inf, NaN, 0 are all representable as doubles, so the
3340          least-significant part can be 0.0.  */
3341       buf[2] = 0;
3342       buf[3] = 0;
3343     }
3344 }
3345
3346 static void
3347 decode_ibm_extended (const struct real_format *fmt ATTRIBUTE_UNUSED, REAL_VALUE_TYPE *r,
3348                      const long *buf)
3349 {
3350   REAL_VALUE_TYPE u, v;
3351   const struct real_format *base_fmt;
3352
3353   base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3354   decode_ieee_double (base_fmt, &u, &buf[0]);
3355
3356   if (u.class != rvc_zero && u.class != rvc_inf && u.class != rvc_nan)
3357     {
3358       decode_ieee_double (base_fmt, &v, &buf[2]);
3359       do_add (r, &u, &v, 0);
3360     }
3361   else
3362     *r = u;
3363 }
3364
3365 const struct real_format ibm_extended_format =
3366   {
3367     encode_ibm_extended,
3368     decode_ibm_extended,
3369     2,
3370     1,
3371     53 + 53,
3372     53,
3373     -1021 + 53,
3374     1024,
3375     -1,
3376     true,
3377     true,
3378     true,
3379     true,
3380     true
3381   };
3382
3383 const struct real_format mips_extended_format =
3384   {
3385     encode_ibm_extended,
3386     decode_ibm_extended,
3387     2,
3388     1,
3389     53 + 53,
3390     53,
3391     -1021 + 53,
3392     1024,
3393     -1,
3394     true,
3395     true,
3396     true,
3397     true,
3398     false
3399   };
3400
3401 \f
3402 /* IEEE quad precision format.  */
3403
3404 static void encode_ieee_quad (const struct real_format *fmt,
3405                               long *, const REAL_VALUE_TYPE *);
3406 static void decode_ieee_quad (const struct real_format *,
3407                               REAL_VALUE_TYPE *, const long *);
3408
3409 static void
3410 encode_ieee_quad (const struct real_format *fmt, long *buf,
3411                   const REAL_VALUE_TYPE *r)
3412 {
3413   unsigned long image3, image2, image1, image0, exp;
3414   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3415   REAL_VALUE_TYPE u;
3416
3417   image3 = r->sign << 31;
3418   image2 = 0;
3419   image1 = 0;
3420   image0 = 0;
3421
3422   rshift_significand (&u, r, SIGNIFICAND_BITS - 113);
3423
3424   switch (r->class)
3425     {
3426     case rvc_zero:
3427       break;
3428
3429     case rvc_inf:
3430       if (fmt->has_inf)
3431         image3 |= 32767 << 16;
3432       else
3433         {
3434           image3 |= 0x7fffffff;
3435           image2 = 0xffffffff;
3436           image1 = 0xffffffff;
3437           image0 = 0xffffffff;
3438         }
3439       break;
3440
3441     case rvc_nan:
3442       if (fmt->has_nans)
3443         {
3444           image3 |= 32767 << 16;
3445
3446           if (r->canonical)
3447             {
3448               /* Don't use bits from the significand.  The
3449                  initialization above is right.  */
3450             }
3451           else if (HOST_BITS_PER_LONG == 32)
3452             {
3453               image0 = u.sig[0];
3454               image1 = u.sig[1];
3455               image2 = u.sig[2];
3456               image3 |= u.sig[3] & 0xffff;
3457             }
3458           else
3459             {
3460               image0 = u.sig[0];
3461               image1 = image0 >> 31 >> 1;
3462               image2 = u.sig[1];
3463               image3 |= (image2 >> 31 >> 1) & 0xffff;
3464               image0 &= 0xffffffff;
3465               image2 &= 0xffffffff;
3466             }
3467           if (r->signalling == fmt->qnan_msb_set)
3468             image3 &= ~0x8000;
3469           else
3470             image3 |= 0x8000;
3471           /* We overload qnan_msb_set here: it's only clear for
3472              mips_ieee_single, which wants all mantissa bits but the
3473              quiet/signalling one set in canonical NaNs (at least
3474              Quiet ones).  */
3475           if (r->canonical && !fmt->qnan_msb_set)
3476             {
3477               image3 |= 0x7fff;
3478               image2 = image1 = image0 = 0xffffffff;
3479             }
3480           else if (((image3 & 0xffff) | image2 | image1 | image0) == 0)
3481             image3 |= 0x4000;
3482         }
3483       else
3484         {
3485           image3 |= 0x7fffffff;
3486           image2 = 0xffffffff;
3487           image1 = 0xffffffff;
3488           image0 = 0xffffffff;
3489         }
3490       break;
3491
3492     case rvc_normal:
3493       /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3494          whereas the intermediate representation is 0.F x 2**exp.
3495          Which means we're off by one.  */
3496       if (denormal)
3497         exp = 0;
3498       else
3499         exp = r->exp + 16383 - 1;
3500       image3 |= exp << 16;
3501
3502       if (HOST_BITS_PER_LONG == 32)
3503         {
3504           image0 = u.sig[0];
3505           image1 = u.sig[1];
3506           image2 = u.sig[2];
3507           image3 |= u.sig[3] & 0xffff;
3508         }
3509       else
3510         {
3511           image0 = u.sig[0];
3512           image1 = image0 >> 31 >> 1;
3513           image2 = u.sig[1];
3514           image3 |= (image2 >> 31 >> 1) & 0xffff;
3515           image0 &= 0xffffffff;
3516           image2 &= 0xffffffff;
3517         }
3518       break;
3519
3520     default:
3521       abort ();
3522     }
3523
3524   if (FLOAT_WORDS_BIG_ENDIAN)
3525     {
3526       buf[0] = image3;
3527       buf[1] = image2;
3528       buf[2] = image1;
3529       buf[3] = image0;
3530     }
3531   else
3532     {
3533       buf[0] = image0;
3534       buf[1] = image1;
3535       buf[2] = image2;
3536       buf[3] = image3;
3537     }
3538 }
3539
3540 static void
3541 decode_ieee_quad (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3542                   const long *buf)
3543 {
3544   unsigned long image3, image2, image1, image0;
3545   bool sign;
3546   int exp;
3547
3548   if (FLOAT_WORDS_BIG_ENDIAN)
3549     {
3550       image3 = buf[0];
3551       image2 = buf[1];
3552       image1 = buf[2];
3553       image0 = buf[3];
3554     }
3555   else
3556     {
3557       image0 = buf[0];
3558       image1 = buf[1];
3559       image2 = buf[2];
3560       image3 = buf[3];
3561     }
3562   image0 &= 0xffffffff;
3563   image1 &= 0xffffffff;
3564   image2 &= 0xffffffff;
3565
3566   sign = (image3 >> 31) & 1;
3567   exp = (image3 >> 16) & 0x7fff;
3568   image3 &= 0xffff;
3569
3570   memset (r, 0, sizeof (*r));
3571
3572   if (exp == 0)
3573     {
3574       if ((image3 | image2 | image1 | image0) && fmt->has_denorm)
3575         {
3576           r->class = rvc_normal;
3577           r->sign = sign;
3578
3579           r->exp = -16382 + (SIGNIFICAND_BITS - 112);
3580           if (HOST_BITS_PER_LONG == 32)
3581             {
3582               r->sig[0] = image0;
3583               r->sig[1] = image1;
3584               r->sig[2] = image2;
3585               r->sig[3] = image3;
3586             }
3587           else
3588             {
3589               r->sig[0] = (image1 << 31 << 1) | image0;
3590               r->sig[1] = (image3 << 31 << 1) | image2;
3591             }
3592
3593           normalize (r);
3594         }
3595       else if (fmt->has_signed_zero)
3596         r->sign = sign;
3597     }
3598   else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3599     {
3600       if (image3 | image2 | image1 | image0)
3601         {
3602           r->class = rvc_nan;
3603           r->sign = sign;
3604           r->signalling = ((image3 >> 15) & 1) ^ fmt->qnan_msb_set;
3605
3606           if (HOST_BITS_PER_LONG == 32)
3607             {
3608               r->sig[0] = image0;
3609               r->sig[1] = image1;
3610               r->sig[2] = image2;
3611               r->sig[3] = image3;
3612             }
3613           else
3614             {
3615               r->sig[0] = (image1 << 31 << 1) | image0;
3616               r->sig[1] = (image3 << 31 << 1) | image2;
3617             }
3618           lshift_significand (r, r, SIGNIFICAND_BITS - 113);
3619         }
3620       else
3621         {
3622           r->class = rvc_inf;
3623           r->sign = sign;
3624         }
3625     }
3626   else
3627     {
3628       r->class = rvc_normal;
3629       r->sign = sign;
3630       r->exp = exp - 16383 + 1;
3631
3632       if (HOST_BITS_PER_LONG == 32)
3633         {
3634           r->sig[0] = image0;
3635           r->sig[1] = image1;
3636           r->sig[2] = image2;
3637           r->sig[3] = image3;
3638         }
3639       else
3640         {
3641           r->sig[0] = (image1 << 31 << 1) | image0;
3642           r->sig[1] = (image3 << 31 << 1) | image2;
3643         }
3644       lshift_significand (r, r, SIGNIFICAND_BITS - 113);
3645       r->sig[SIGSZ-1] |= SIG_MSB;
3646     }
3647 }
3648
3649 const struct real_format ieee_quad_format =
3650   {
3651     encode_ieee_quad,
3652     decode_ieee_quad,
3653     2,
3654     1,
3655     113,
3656     113,
3657     -16381,
3658     16384,
3659     127,
3660     true,
3661     true,
3662     true,
3663     true,
3664     true
3665   };
3666
3667 const struct real_format mips_quad_format =
3668   {
3669     encode_ieee_quad,
3670     decode_ieee_quad,
3671     2,
3672     1,
3673     113,
3674     113,
3675     -16381,
3676     16384,
3677     127,
3678     true,
3679     true,
3680     true,
3681     true,
3682     false
3683   };
3684 \f
3685 /* Descriptions of VAX floating point formats can be found beginning at
3686
3687    http://h71000.www7.hp.com/doc/73FINAL/4515/4515pro_013.html#f_floating_point_format
3688
3689    The thing to remember is that they're almost IEEE, except for word
3690    order, exponent bias, and the lack of infinities, nans, and denormals.
3691
3692    We don't implement the H_floating format here, simply because neither
3693    the VAX or Alpha ports use it.  */
3694
3695 static void encode_vax_f (const struct real_format *fmt,
3696                           long *, const REAL_VALUE_TYPE *);
3697 static void decode_vax_f (const struct real_format *,
3698                           REAL_VALUE_TYPE *, const long *);
3699 static void encode_vax_d (const struct real_format *fmt,
3700                           long *, const REAL_VALUE_TYPE *);
3701 static void decode_vax_d (const struct real_format *,
3702                           REAL_VALUE_TYPE *, const long *);
3703 static void encode_vax_g (const struct real_format *fmt,
3704                           long *, const REAL_VALUE_TYPE *);
3705 static void decode_vax_g (const struct real_format *,
3706                           REAL_VALUE_TYPE *, const long *);
3707
3708 static void
3709 encode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3710               const REAL_VALUE_TYPE *r)
3711 {
3712   unsigned long sign, exp, sig, image;
3713
3714   sign = r->sign << 15;
3715
3716   switch (r->class)
3717     {
3718     case rvc_zero:
3719       image = 0;
3720       break;
3721
3722     case rvc_inf:
3723     case rvc_nan:
3724       image = 0xffff7fff | sign;
3725       break;
3726
3727     case rvc_normal:
3728       sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
3729       exp = r->exp + 128;
3730
3731       image = (sig << 16) & 0xffff0000;
3732       image |= sign;
3733       image |= exp << 7;
3734       image |= sig >> 16;
3735       break;
3736
3737     default:
3738       abort ();
3739     }
3740
3741   buf[0] = image;
3742 }
3743
3744 static void
3745 decode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED,
3746               REAL_VALUE_TYPE *r, const long *buf)
3747 {
3748   unsigned long image = buf[0] & 0xffffffff;
3749   int exp = (image >> 7) & 0xff;
3750
3751   memset (r, 0, sizeof (*r));
3752
3753   if (exp != 0)
3754     {
3755       r->class = rvc_normal;
3756       r->sign = (image >> 15) & 1;
3757       r->exp = exp - 128;
3758
3759       image = ((image & 0x7f) << 16) | ((image >> 16) & 0xffff);
3760       r->sig[SIGSZ-1] = (image << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
3761     }
3762 }
3763
3764 static void
3765 encode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3766               const REAL_VALUE_TYPE *r)
3767 {
3768   unsigned long image0, image1, sign = r->sign << 15;
3769
3770   switch (r->class)
3771     {
3772     case rvc_zero:
3773       image0 = image1 = 0;
3774       break;
3775
3776     case rvc_inf:
3777     case rvc_nan:
3778       image0 = 0xffff7fff | sign;
3779       image1 = 0xffffffff;
3780       break;
3781
3782     case rvc_normal:
3783       /* Extract the significand into straight hi:lo.  */
3784       if (HOST_BITS_PER_LONG == 64)
3785         {
3786           image0 = r->sig[SIGSZ-1];
3787           image1 = (image0 >> (64 - 56)) & 0xffffffff;
3788           image0 = (image0 >> (64 - 56 + 1) >> 31) & 0x7fffff;
3789         }
3790       else
3791         {
3792           image0 = r->sig[SIGSZ-1];
3793           image1 = r->sig[SIGSZ-2];
3794           image1 = (image0 << 24) | (image1 >> 8);
3795           image0 = (image0 >> 8) & 0xffffff;
3796         }
3797
3798       /* Rearrange the half-words of the significand to match the
3799          external format.  */
3800       image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff007f;
3801       image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
3802
3803       /* Add the sign and exponent.  */
3804       image0 |= sign;
3805       image0 |= (r->exp + 128) << 7;
3806       break;
3807
3808     default:
3809       abort ();
3810     }
3811
3812   if (FLOAT_WORDS_BIG_ENDIAN)
3813     buf[0] = image1, buf[1] = image0;
3814   else
3815     buf[0] = image0, buf[1] = image1;
3816 }
3817
3818 static void
3819 decode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED,
3820               REAL_VALUE_TYPE *r, const long *buf)
3821 {
3822   unsigned long image0, image1;
3823   int exp;
3824
3825   if (FLOAT_WORDS_BIG_ENDIAN)
3826     image1 = buf[0], image0 = buf[1];
3827   else
3828     image0 = buf[0], image1 = buf[1];
3829   image0 &= 0xffffffff;
3830   image1 &= 0xffffffff;
3831
3832   exp = (image0 >> 7) & 0xff;
3833
3834   memset (r, 0, sizeof (*r));
3835
3836   if (exp != 0)
3837     {
3838       r->class = rvc_normal;
3839       r->sign = (image0 >> 15) & 1;
3840       r->exp = exp - 128;
3841
3842       /* Rearrange the half-words of the external format into
3843          proper ascending order.  */
3844       image0 = ((image0 & 0x7f) << 16) | ((image0 >> 16) & 0xffff);
3845       image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
3846
3847       if (HOST_BITS_PER_LONG == 64)
3848         {
3849           image0 = (image0 << 31 << 1) | image1;
3850           image0 <<= 64 - 56;
3851           image0 |= SIG_MSB;
3852           r->sig[SIGSZ-1] = image0;
3853         }
3854       else
3855         {
3856           r->sig[SIGSZ-1] = image0;
3857           r->sig[SIGSZ-2] = image1;
3858           lshift_significand (r, r, 2*HOST_BITS_PER_LONG - 56);
3859           r->sig[SIGSZ-1] |= SIG_MSB;
3860         }
3861     }
3862 }
3863
3864 static void
3865 encode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
3866               const REAL_VALUE_TYPE *r)
3867 {
3868   unsigned long image0, image1, sign = r->sign << 15;
3869
3870   switch (r->class)
3871     {
3872     case rvc_zero:
3873       image0 = image1 = 0;
3874       break;
3875
3876     case rvc_inf:
3877     case rvc_nan:
3878       image0 = 0xffff7fff | sign;
3879       image1 = 0xffffffff;
3880       break;
3881
3882     case rvc_normal:
3883       /* Extract the significand into straight hi:lo.  */
3884       if (HOST_BITS_PER_LONG == 64)
3885         {
3886           image0 = r->sig[SIGSZ-1];
3887           image1 = (image0 >> (64 - 53)) & 0xffffffff;
3888           image0 = (image0 >> (64 - 53 + 1) >> 31) & 0xfffff;
3889         }
3890       else
3891         {
3892           image0 = r->sig[SIGSZ-1];
3893           image1 = r->sig[SIGSZ-2];
3894           image1 = (image0 << 21) | (image1 >> 11);
3895           image0 = (image0 >> 11) & 0xfffff;
3896         }
3897
3898       /* Rearrange the half-words of the significand to match the
3899          external format.  */
3900       image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff000f;
3901       image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
3902
3903       /* Add the sign and exponent.  */
3904       image0 |= sign;
3905       image0 |= (r->exp + 1024) << 4;
3906       break;
3907
3908     default:
3909       abort ();
3910     }
3911
3912   if (FLOAT_WORDS_BIG_ENDIAN)
3913     buf[0] = image1, buf[1] = image0;
3914   else
3915     buf[0] = image0, buf[1] = image1;
3916 }
3917
3918 static void
3919 decode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED,
3920               REAL_VALUE_TYPE *r, const long *buf)
3921 {
3922   unsigned long image0, image1;
3923   int exp;
3924
3925   if (FLOAT_WORDS_BIG_ENDIAN)
3926     image1 = buf[0], image0 = buf[1];
3927   else
3928     image0 = buf[0], image1 = buf[1];
3929   image0 &= 0xffffffff;
3930   image1 &= 0xffffffff;
3931
3932   exp = (image0 >> 4) & 0x7ff;
3933
3934   memset (r, 0, sizeof (*r));
3935
3936   if (exp != 0)
3937     {
3938       r->class = rvc_normal;
3939       r->sign = (image0 >> 15) & 1;
3940       r->exp = exp - 1024;
3941
3942       /* Rearrange the half-words of the external format into
3943          proper ascending order.  */
3944       image0 = ((image0 & 0xf) << 16) | ((image0 >> 16) & 0xffff);
3945       image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
3946
3947       if (HOST_BITS_PER_LONG == 64)
3948         {
3949           image0 = (image0 << 31 << 1) | image1;
3950           image0 <<= 64 - 53;
3951           image0 |= SIG_MSB;
3952           r->sig[SIGSZ-1] = image0;
3953         }
3954       else
3955         {
3956           r->sig[SIGSZ-1] = image0;
3957           r->sig[SIGSZ-2] = image1;
3958           lshift_significand (r, r, 64 - 53);
3959           r->sig[SIGSZ-1] |= SIG_MSB;
3960         }
3961     }
3962 }
3963
3964 const struct real_format vax_f_format =
3965   {
3966     encode_vax_f,
3967     decode_vax_f,
3968     2,
3969     1,
3970     24,
3971     24,
3972     -127,
3973     127,
3974     15,
3975     false,
3976     false,
3977     false,
3978     false,
3979     false
3980   };
3981
3982 const struct real_format vax_d_format =
3983   {
3984     encode_vax_d,
3985     decode_vax_d,
3986     2,
3987     1,
3988     56,
3989     56,
3990     -127,
3991     127,
3992     15,
3993     false,
3994     false,
3995     false,
3996     false,
3997     false
3998   };
3999
4000 const struct real_format vax_g_format =
4001   {
4002     encode_vax_g,
4003     decode_vax_g,
4004     2,
4005     1,
4006     53,
4007     53,
4008     -1023,
4009     1023,
4010     15,
4011     false,
4012     false,
4013     false,
4014     false,
4015     false
4016   };
4017 \f
4018 /* A good reference for these can be found in chapter 9 of
4019    "ESA/390 Principles of Operation", IBM document number SA22-7201-01.
4020    An on-line version can be found here:
4021
4022    http://publibz.boulder.ibm.com/cgi-bin/bookmgr_OS390/BOOKS/DZ9AR001/9.1?DT=19930923083613
4023 */
4024
4025 static void encode_i370_single (const struct real_format *fmt,
4026                                 long *, const REAL_VALUE_TYPE *);
4027 static void decode_i370_single (const struct real_format *,
4028                                 REAL_VALUE_TYPE *, const long *);
4029 static void encode_i370_double (const struct real_format *fmt,
4030                                 long *, const REAL_VALUE_TYPE *);
4031 static void decode_i370_double (const struct real_format *,
4032                                 REAL_VALUE_TYPE *, const long *);
4033
4034 static void
4035 encode_i370_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4036                     long *buf, const REAL_VALUE_TYPE *r)
4037 {
4038   unsigned long sign, exp, sig, image;
4039
4040   sign = r->sign << 31;
4041
4042   switch (r->class)
4043     {
4044     case rvc_zero:
4045       image = 0;
4046       break;
4047
4048     case rvc_inf:
4049     case rvc_nan:
4050       image = 0x7fffffff | sign;
4051       break;
4052
4053     case rvc_normal:
4054       sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0xffffff;
4055       exp = ((r->exp / 4) + 64) << 24;
4056       image = sign | exp | sig;
4057       break;
4058
4059     default:
4060       abort ();
4061     }
4062
4063   buf[0] = image;
4064 }
4065
4066 static void
4067 decode_i370_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4068                     REAL_VALUE_TYPE *r, const long *buf)
4069 {
4070   unsigned long sign, sig, image = buf[0];
4071   int exp;
4072
4073   sign = (image >> 31) & 1;
4074   exp = (image >> 24) & 0x7f;
4075   sig = image & 0xffffff;
4076
4077   memset (r, 0, sizeof (*r));
4078
4079   if (exp || sig)
4080     {
4081       r->class = rvc_normal;
4082       r->sign = sign;
4083       r->exp = (exp - 64) * 4;
4084       r->sig[SIGSZ-1] = sig << (HOST_BITS_PER_LONG - 24);
4085       normalize (r);
4086     }
4087 }
4088
4089 static void
4090 encode_i370_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4091                     long *buf, const REAL_VALUE_TYPE *r)
4092 {
4093   unsigned long sign, exp, image_hi, image_lo;
4094
4095   sign = r->sign << 31;
4096
4097   switch (r->class)
4098     {
4099     case rvc_zero:
4100       image_hi = image_lo = 0;
4101       break;
4102
4103     case rvc_inf:
4104     case rvc_nan:
4105       image_hi = 0x7fffffff | sign;
4106       image_lo = 0xffffffff;
4107       break;
4108
4109     case rvc_normal:
4110       if (HOST_BITS_PER_LONG == 64)
4111         {
4112           image_hi = r->sig[SIGSZ-1];
4113           image_lo = (image_hi >> (64 - 56)) & 0xffffffff;
4114           image_hi = (image_hi >> (64 - 56 + 1) >> 31) & 0xffffff;
4115         }
4116       else
4117         {
4118           image_hi = r->sig[SIGSZ-1];
4119           image_lo = r->sig[SIGSZ-2];
4120           image_lo = (image_lo >> 8) | (image_hi << 24);
4121           image_hi >>= 8;
4122         }
4123
4124       exp = ((r->exp / 4) + 64) << 24;
4125       image_hi |= sign | exp;
4126       break;
4127
4128     default:
4129       abort ();
4130     }
4131
4132   if (FLOAT_WORDS_BIG_ENDIAN)
4133     buf[0] = image_hi, buf[1] = image_lo;
4134   else
4135     buf[0] = image_lo, buf[1] = image_hi;
4136 }
4137
4138 static void
4139 decode_i370_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4140                     REAL_VALUE_TYPE *r, const long *buf)
4141 {
4142   unsigned long sign, image_hi, image_lo;
4143   int exp;
4144
4145   if (FLOAT_WORDS_BIG_ENDIAN)
4146     image_hi = buf[0], image_lo = buf[1];
4147   else
4148     image_lo = buf[0], image_hi = buf[1];
4149
4150   sign = (image_hi >> 31) & 1;
4151   exp = (image_hi >> 24) & 0x7f;
4152   image_hi &= 0xffffff;
4153   image_lo &= 0xffffffff;
4154
4155   memset (r, 0, sizeof (*r));
4156
4157   if (exp || image_hi || image_lo)
4158     {
4159       r->class = rvc_normal;
4160       r->sign = sign;
4161       r->exp = (exp - 64) * 4 + (SIGNIFICAND_BITS - 56);
4162
4163       if (HOST_BITS_PER_LONG == 32)
4164         {
4165           r->sig[0] = image_lo;
4166           r->sig[1] = image_hi;
4167         }
4168       else
4169         r->sig[0] = image_lo | (image_hi << 31 << 1);
4170
4171       normalize (r);
4172     }
4173 }
4174
4175 const struct real_format i370_single_format =
4176   {
4177     encode_i370_single,
4178     decode_i370_single,
4179     16,
4180     4,
4181     6,
4182     6,
4183     -64,
4184     63,
4185     31,
4186     false,
4187     false,
4188     false, /* ??? The encoding does allow for "unnormals".  */
4189     false, /* ??? The encoding does allow for "unnormals".  */
4190     false
4191   };
4192
4193 const struct real_format i370_double_format =
4194   {
4195     encode_i370_double,
4196     decode_i370_double,
4197     16,
4198     4,
4199     14,
4200     14,
4201     -64,
4202     63,
4203     63,
4204     false,
4205     false,
4206     false, /* ??? The encoding does allow for "unnormals".  */
4207     false, /* ??? The encoding does allow for "unnormals".  */
4208     false
4209   };
4210 \f
4211 /* The "twos-complement" c4x format is officially defined as
4212
4213         x = s(~s).f * 2**e
4214
4215    This is rather misleading.  One must remember that F is signed.
4216    A better description would be
4217
4218         x = -1**s * ((s + 1 + .f) * 2**e
4219
4220    So if we have a (4 bit) fraction of .1000 with a sign bit of 1,
4221    that's -1 * (1+1+(-.5)) == -1.5.  I think.
4222
4223    The constructions here are taken from Tables 5-1 and 5-2 of the
4224    TMS320C4x User's Guide wherein step-by-step instructions for
4225    conversion from IEEE are presented.  That's close enough to our
4226    internal representation so as to make things easy.
4227
4228    See http://www-s.ti.com/sc/psheets/spru063c/spru063c.pdf  */
4229
4230 static void encode_c4x_single (const struct real_format *fmt,
4231                                long *, const REAL_VALUE_TYPE *);
4232 static void decode_c4x_single (const struct real_format *,
4233                                REAL_VALUE_TYPE *, const long *);
4234 static void encode_c4x_extended (const struct real_format *fmt,
4235                                  long *, const REAL_VALUE_TYPE *);
4236 static void decode_c4x_extended (const struct real_format *,
4237                                  REAL_VALUE_TYPE *, const long *);
4238
4239 static void
4240 encode_c4x_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4241                    long *buf, const REAL_VALUE_TYPE *r)
4242 {
4243   unsigned long image, exp, sig;
4244
4245   switch (r->class)
4246     {
4247     case rvc_zero:
4248       exp = -128;
4249       sig = 0;
4250       break;
4251
4252     case rvc_inf:
4253     case rvc_nan:
4254       exp = 127;
4255       sig = 0x800000 - r->sign;
4256       break;
4257
4258     case rvc_normal:
4259       exp = r->exp - 1;
4260       sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
4261       if (r->sign)
4262         {
4263           if (sig)
4264             sig = -sig;
4265           else
4266             exp--;
4267           sig |= 0x800000;
4268         }
4269       break;
4270
4271     default:
4272       abort ();
4273     }
4274
4275   image = ((exp & 0xff) << 24) | (sig & 0xffffff);
4276   buf[0] = image;
4277 }
4278
4279 static void
4280 decode_c4x_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4281                    REAL_VALUE_TYPE *r, const long *buf)
4282 {
4283   unsigned long image = buf[0];
4284   unsigned long sig;
4285   int exp, sf;
4286
4287   exp = (((image >> 24) & 0xff) ^ 0x80) - 0x80;
4288   sf = ((image & 0xffffff) ^ 0x800000) - 0x800000;
4289
4290   memset (r, 0, sizeof (*r));
4291
4292   if (exp != -128)
4293     {
4294       r->class = rvc_normal;
4295
4296       sig = sf & 0x7fffff;
4297       if (sf < 0)
4298         {
4299           r->sign = 1;
4300           if (sig)
4301             sig = -sig;
4302           else
4303             exp++;
4304         }
4305       sig = (sig << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
4306
4307       r->exp = exp + 1;
4308       r->sig[SIGSZ-1] = sig;
4309     }
4310 }
4311
4312 static void
4313 encode_c4x_extended (const struct real_format *fmt ATTRIBUTE_UNUSED,
4314                      long *buf, const REAL_VALUE_TYPE *r)
4315 {
4316   unsigned long exp, sig;
4317
4318   switch (r->class)
4319     {
4320     case rvc_zero:
4321       exp = -128;
4322       sig = 0;
4323       break;
4324
4325     case rvc_inf:
4326     case rvc_nan:
4327       exp = 127;
4328       sig = 0x80000000 - r->sign;
4329       break;
4330
4331     case rvc_normal:
4332       exp = r->exp - 1;
4333
4334       sig = r->sig[SIGSZ-1];
4335       if (HOST_BITS_PER_LONG == 64)
4336         sig = sig >> 1 >> 31;
4337       sig &= 0x7fffffff;
4338
4339       if (r->sign)
4340         {
4341           if (sig)
4342             sig = -sig;
4343           else
4344             exp--;
4345           sig |= 0x80000000;
4346         }
4347       break;
4348
4349     default:
4350       abort ();
4351     }
4352
4353   exp = (exp & 0xff) << 24;
4354   sig &= 0xffffffff;
4355
4356   if (FLOAT_WORDS_BIG_ENDIAN)
4357     buf[0] = exp, buf[1] = sig;
4358   else
4359     buf[0] = sig, buf[0] = exp;
4360 }
4361
4362 static void
4363 decode_c4x_extended (const struct real_format *fmt ATTRIBUTE_UNUSED,
4364                      REAL_VALUE_TYPE *r, const long *buf)
4365 {
4366   unsigned long sig;
4367   int exp, sf;
4368
4369   if (FLOAT_WORDS_BIG_ENDIAN)
4370     exp = buf[0], sf = buf[1];
4371   else
4372     sf = buf[0], exp = buf[1];
4373
4374   exp = (((exp >> 24) & 0xff) & 0x80) - 0x80;
4375   sf = ((sf & 0xffffffff) ^ 0x80000000) - 0x80000000;
4376
4377   memset (r, 0, sizeof (*r));
4378
4379   if (exp != -128)
4380     {
4381       r->class = rvc_normal;
4382
4383       sig = sf & 0x7fffffff;
4384       if (sf < 0)
4385         {
4386           r->sign = 1;
4387           if (sig)
4388             sig = -sig;
4389           else
4390             exp++;
4391         }
4392       if (HOST_BITS_PER_LONG == 64)
4393         sig = sig << 1 << 31;
4394       sig |= SIG_MSB;
4395
4396       r->exp = exp + 1;
4397       r->sig[SIGSZ-1] = sig;
4398     }
4399 }
4400
4401 const struct real_format c4x_single_format =
4402   {
4403     encode_c4x_single,
4404     decode_c4x_single,
4405     2,
4406     1,
4407     24,
4408     24,
4409     -126,
4410     128,
4411     -1,
4412     false,
4413     false,
4414     false,
4415     false,
4416     false
4417   };
4418
4419 const struct real_format c4x_extended_format =
4420   {
4421     encode_c4x_extended,
4422     decode_c4x_extended,
4423     2,
4424     1,
4425     32,
4426     32,
4427     -126,
4428     128,
4429     -1,
4430     false,
4431     false,
4432     false,
4433     false,
4434     false
4435   };
4436
4437 \f
4438 /* A synthetic "format" for internal arithmetic.  It's the size of the
4439    internal significand minus the two bits needed for proper rounding.
4440    The encode and decode routines exist only to satisfy our paranoia
4441    harness.  */
4442
4443 static void encode_internal (const struct real_format *fmt,
4444                              long *, const REAL_VALUE_TYPE *);
4445 static void decode_internal (const struct real_format *,
4446                              REAL_VALUE_TYPE *, const long *);
4447
4448 static void
4449 encode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4450                  const REAL_VALUE_TYPE *r)
4451 {
4452   memcpy (buf, r, sizeof (*r));
4453 }
4454
4455 static void
4456 decode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED,
4457                  REAL_VALUE_TYPE *r, const long *buf)
4458 {
4459   memcpy (r, buf, sizeof (*r));
4460 }
4461
4462 const struct real_format real_internal_format =
4463   {
4464     encode_internal,
4465     decode_internal,
4466     2,
4467     1,
4468     SIGNIFICAND_BITS - 2,
4469     SIGNIFICAND_BITS - 2,
4470     -MAX_EXP,
4471     MAX_EXP,
4472     -1,
4473     true,
4474     true,
4475     false,
4476     true,
4477     true
4478   };
4479 \f
4480 /* Calculate the square root of X in mode MODE, and store the result
4481    in R.  Return TRUE if the operation does not raise an exception.
4482    For details see "High Precision Division and Square Root",
4483    Alan H. Karp and Peter Markstein, HP Lab Report 93-93-42, June
4484    1993.  http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf.  */
4485
4486 bool
4487 real_sqrt (REAL_VALUE_TYPE *r, enum machine_mode mode,
4488            const REAL_VALUE_TYPE *x)
4489 {
4490   static REAL_VALUE_TYPE halfthree;
4491   static bool init = false;
4492   REAL_VALUE_TYPE h, t, i;
4493   int iter, exp;
4494
4495   /* sqrt(-0.0) is -0.0.  */
4496   if (real_isnegzero (x))
4497     {
4498       *r = *x;
4499       return false;
4500     }
4501
4502   /* Negative arguments return NaN.  */
4503   if (real_isneg (x))
4504     {
4505       get_canonical_qnan (r, 0);
4506       return false;
4507     }
4508
4509   /* Infinity and NaN return themselves.  */
4510   if (real_isinf (x) || real_isnan (x))
4511     {
4512       *r = *x;
4513       return false;
4514     }
4515
4516   if (!init)
4517     {
4518       do_add (&halfthree, &dconst1, &dconsthalf, 0);
4519       init = true;
4520     }
4521
4522   /* Initial guess for reciprocal sqrt, i.  */
4523   exp = real_exponent (x);
4524   real_ldexp (&i, &dconst1, -exp/2);
4525
4526   /* Newton's iteration for reciprocal sqrt, i.  */
4527   for (iter = 0; iter < 16; iter++)
4528     {
4529       /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x).  */
4530       do_multiply (&t, x, &i);
4531       do_multiply (&h, &t, &i);
4532       do_multiply (&t, &h, &dconsthalf);
4533       do_add (&h, &halfthree, &t, 1);
4534       do_multiply (&t, &i, &h);
4535
4536       /* Check for early convergence.  */
4537       if (iter >= 6 && real_identical (&i, &t))
4538         break;
4539
4540       /* ??? Unroll loop to avoid copying.  */
4541       i = t;
4542     }
4543
4544   /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)).  */
4545   do_multiply (&t, x, &i);
4546   do_multiply (&h, &t, &i);
4547   do_add (&i, &dconst1, &h, 1);
4548   do_multiply (&h, &t, &i);
4549   do_multiply (&i, &dconsthalf, &h);
4550   do_add (&h, &t, &i, 0);
4551
4552   /* ??? We need a Tuckerman test to get the last bit.  */
4553
4554   real_convert (r, mode, &h);
4555   return true;
4556 }
4557
4558 /* Calculate X raised to the integer exponent N in mode MODE and store
4559    the result in R.  Return true if the result may be inexact due to
4560    loss of precision.  The algorithm is t