OSDN Git Service

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