OSDN Git Service

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