OSDN Git Service

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