remove unused files
[platform/upstream/gcc48.git] / gcc / real.c
1 /* real.c - software floating point emulation.
2    Copyright (C) 1993-2013 Free Software Foundation, Inc.
3    Contributed by Stephen L. Moshier (moshier@world.std.com).
4    Re-written by Richard Henderson <rth@redhat.com>
5
6    This file is part of GCC.
7
8    GCC is free software; you can redistribute it and/or modify it under
9    the terms of the GNU General Public License as published by the Free
10    Software Foundation; either version 3, or (at your option) any later
11    version.
12
13    GCC is distributed in the hope that it will be useful, but WITHOUT ANY
14    WARRANTY; without even the implied warranty of MERCHANTABILITY or
15    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
16    for more details.
17
18    You should have received a copy of the GNU General Public License
19    along with GCC; see the file COPYING3.  If not see
20    <http://www.gnu.org/licenses/>.  */
21
22 #include "config.h"
23 #include "system.h"
24 #include "coretypes.h"
25 #include "tm.h"
26 #include "tree.h"
27 #include "diagnostic-core.h"
28 #include "real.h"
29 #include "realmpfr.h"
30 #include "tm_p.h"
31 #include "dfp.h"
32
33 /* The floating point model used internally is not exactly IEEE 754
34    compliant, and close to the description in the ISO C99 standard,
35    section 5.2.4.2.2 Characteristics of floating types.
36
37    Specifically
38
39         x = s * b^e * \sum_{k=1}^p f_k * b^{-k}
40
41         where
42                 s = sign (+- 1)
43                 b = base or radix, here always 2
44                 e = exponent
45                 p = precision (the number of base-b digits in the significand)
46                 f_k = the digits of the significand.
47
48    We differ from typical IEEE 754 encodings in that the entire
49    significand is fractional.  Normalized significands are in the
50    range [0.5, 1.0).
51
52    A requirement of the model is that P be larger than the largest
53    supported target floating-point type by at least 2 bits.  This gives
54    us proper rounding when we truncate to the target type.  In addition,
55    E must be large enough to hold the smallest supported denormal number
56    in a normalized form.
57
58    Both of these requirements are easily satisfied.  The largest target
59    significand is 113 bits; we store at least 160.  The smallest
60    denormal number fits in 17 exponent bits; we store 26.
61
62    Note that the decimal string conversion routines are sensitive to
63    rounding errors.  Since the raw arithmetic routines do not themselves
64    have guard digits or rounding, the computation of 10**exp can
65    accumulate more than a few digits of error.  The previous incarnation
66    of real.c successfully used a 144-bit fraction; given the current
67    layout of REAL_VALUE_TYPE we're forced to expand to at least 160 bits.  */
68
69
70 /* Used to classify two numbers simultaneously.  */
71 #define CLASS2(A, B)  ((A) << 2 | (B))
72
73 #if HOST_BITS_PER_LONG != 64 && HOST_BITS_PER_LONG != 32
74  #error "Some constant folding done by hand to avoid shift count warnings"
75 #endif
76
77 static void get_zero (REAL_VALUE_TYPE *, int);
78 static void get_canonical_qnan (REAL_VALUE_TYPE *, int);
79 static void get_canonical_snan (REAL_VALUE_TYPE *, int);
80 static void get_inf (REAL_VALUE_TYPE *, int);
81 static bool sticky_rshift_significand (REAL_VALUE_TYPE *,
82                                        const REAL_VALUE_TYPE *, unsigned int);
83 static void rshift_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
84                                 unsigned int);
85 static void lshift_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
86                                 unsigned int);
87 static void lshift_significand_1 (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
88 static bool add_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *,
89                               const REAL_VALUE_TYPE *);
90 static bool sub_significands (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
91                               const REAL_VALUE_TYPE *, int);
92 static void neg_significand (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
93 static int cmp_significands (const REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
94 static int cmp_significand_0 (const REAL_VALUE_TYPE *);
95 static void set_significand_bit (REAL_VALUE_TYPE *, unsigned int);
96 static void clear_significand_bit (REAL_VALUE_TYPE *, unsigned int);
97 static bool test_significand_bit (REAL_VALUE_TYPE *, unsigned int);
98 static void clear_significand_below (REAL_VALUE_TYPE *, unsigned int);
99 static bool div_significands (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
100                               const REAL_VALUE_TYPE *);
101 static void normalize (REAL_VALUE_TYPE *);
102
103 static bool do_add (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
104                     const REAL_VALUE_TYPE *, int);
105 static bool do_multiply (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
106                          const REAL_VALUE_TYPE *);
107 static bool do_divide (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *,
108                        const REAL_VALUE_TYPE *);
109 static int do_compare (const REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *, int);
110 static void do_fix_trunc (REAL_VALUE_TYPE *, const REAL_VALUE_TYPE *);
111
112 static unsigned long rtd_divmod (REAL_VALUE_TYPE *, REAL_VALUE_TYPE *);
113 static void decimal_from_integer (REAL_VALUE_TYPE *);
114 static void decimal_integer_string (char *, const REAL_VALUE_TYPE *,
115                                     size_t);
116
117 static const REAL_VALUE_TYPE * ten_to_ptwo (int);
118 static const REAL_VALUE_TYPE * ten_to_mptwo (int);
119 static const REAL_VALUE_TYPE * real_digit (int);
120 static void times_pten (REAL_VALUE_TYPE *, int);
121
122 static void round_for_format (const struct real_format *, REAL_VALUE_TYPE *);
123 \f
124 /* Initialize R with a positive zero.  */
125
126 static inline void
127 get_zero (REAL_VALUE_TYPE *r, int sign)
128 {
129   memset (r, 0, sizeof (*r));
130   r->sign = sign;
131 }
132
133 /* Initialize R with the canonical quiet NaN.  */
134
135 static inline void
136 get_canonical_qnan (REAL_VALUE_TYPE *r, int sign)
137 {
138   memset (r, 0, sizeof (*r));
139   r->cl = rvc_nan;
140   r->sign = sign;
141   r->canonical = 1;
142 }
143
144 static inline void
145 get_canonical_snan (REAL_VALUE_TYPE *r, int sign)
146 {
147   memset (r, 0, sizeof (*r));
148   r->cl = rvc_nan;
149   r->sign = sign;
150   r->signalling = 1;
151   r->canonical = 1;
152 }
153
154 static inline void
155 get_inf (REAL_VALUE_TYPE *r, int sign)
156 {
157   memset (r, 0, sizeof (*r));
158   r->cl = rvc_inf;
159   r->sign = sign;
160 }
161
162 \f
163 /* Right-shift the significand of A by N bits; put the result in the
164    significand of R.  If any one bits are shifted out, return true.  */
165
166 static bool
167 sticky_rshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
168                            unsigned int n)
169 {
170   unsigned long sticky = 0;
171   unsigned int i, ofs = 0;
172
173   if (n >= HOST_BITS_PER_LONG)
174     {
175       for (i = 0, ofs = n / HOST_BITS_PER_LONG; i < ofs; ++i)
176         sticky |= a->sig[i];
177       n &= HOST_BITS_PER_LONG - 1;
178     }
179
180   if (n != 0)
181     {
182       sticky |= a->sig[ofs] & (((unsigned long)1 << n) - 1);
183       for (i = 0; i < SIGSZ; ++i)
184         {
185           r->sig[i]
186             = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
187                | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
188                   << (HOST_BITS_PER_LONG - n)));
189         }
190     }
191   else
192     {
193       for (i = 0; ofs + i < SIGSZ; ++i)
194         r->sig[i] = a->sig[ofs + i];
195       for (; i < SIGSZ; ++i)
196         r->sig[i] = 0;
197     }
198
199   return sticky != 0;
200 }
201
202 /* Right-shift the significand of A by N bits; put the result in the
203    significand of R.  */
204
205 static void
206 rshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
207                     unsigned int n)
208 {
209   unsigned int i, ofs = n / HOST_BITS_PER_LONG;
210
211   n &= HOST_BITS_PER_LONG - 1;
212   if (n != 0)
213     {
214       for (i = 0; i < SIGSZ; ++i)
215         {
216           r->sig[i]
217             = (((ofs + i >= SIGSZ ? 0 : a->sig[ofs + i]) >> n)
218                | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[ofs + i + 1])
219                   << (HOST_BITS_PER_LONG - n)));
220         }
221     }
222   else
223     {
224       for (i = 0; ofs + i < SIGSZ; ++i)
225         r->sig[i] = a->sig[ofs + i];
226       for (; i < SIGSZ; ++i)
227         r->sig[i] = 0;
228     }
229 }
230
231 /* Left-shift the significand of A by N bits; put the result in the
232    significand of R.  */
233
234 static void
235 lshift_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
236                     unsigned int n)
237 {
238   unsigned int i, ofs = n / HOST_BITS_PER_LONG;
239
240   n &= HOST_BITS_PER_LONG - 1;
241   if (n == 0)
242     {
243       for (i = 0; ofs + i < SIGSZ; ++i)
244         r->sig[SIGSZ-1-i] = a->sig[SIGSZ-1-i-ofs];
245       for (; i < SIGSZ; ++i)
246         r->sig[SIGSZ-1-i] = 0;
247     }
248   else
249     for (i = 0; i < SIGSZ; ++i)
250       {
251         r->sig[SIGSZ-1-i]
252           = (((ofs + i >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs]) << n)
253              | ((ofs + i + 1 >= SIGSZ ? 0 : a->sig[SIGSZ-1-i-ofs-1])
254                 >> (HOST_BITS_PER_LONG - n)));
255       }
256 }
257
258 /* Likewise, but N is specialized to 1.  */
259
260 static inline void
261 lshift_significand_1 (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
262 {
263   unsigned int i;
264
265   for (i = SIGSZ - 1; i > 0; --i)
266     r->sig[i] = (a->sig[i] << 1) | (a->sig[i-1] >> (HOST_BITS_PER_LONG - 1));
267   r->sig[0] = a->sig[0] << 1;
268 }
269
270 /* Add the significands of A and B, placing the result in R.  Return
271    true if there was carry out of the most significant word.  */
272
273 static inline bool
274 add_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
275                   const REAL_VALUE_TYPE *b)
276 {
277   bool carry = false;
278   int i;
279
280   for (i = 0; i < SIGSZ; ++i)
281     {
282       unsigned long ai = a->sig[i];
283       unsigned long ri = ai + b->sig[i];
284
285       if (carry)
286         {
287           carry = ri < ai;
288           carry |= ++ri == 0;
289         }
290       else
291         carry = ri < ai;
292
293       r->sig[i] = ri;
294     }
295
296   return carry;
297 }
298
299 /* Subtract the significands of A and B, placing the result in R.  CARRY is
300    true if there's a borrow incoming to the least significant word.
301    Return true if there was borrow out of the most significant word.  */
302
303 static inline bool
304 sub_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
305                   const REAL_VALUE_TYPE *b, int carry)
306 {
307   int i;
308
309   for (i = 0; i < SIGSZ; ++i)
310     {
311       unsigned long ai = a->sig[i];
312       unsigned long ri = ai - b->sig[i];
313
314       if (carry)
315         {
316           carry = ri > ai;
317           carry |= ~--ri == 0;
318         }
319       else
320         carry = ri > ai;
321
322       r->sig[i] = ri;
323     }
324
325   return carry;
326 }
327
328 /* Negate the significand A, placing the result in R.  */
329
330 static inline void
331 neg_significand (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
332 {
333   bool carry = true;
334   int i;
335
336   for (i = 0; i < SIGSZ; ++i)
337     {
338       unsigned long ri, ai = a->sig[i];
339
340       if (carry)
341         {
342           if (ai)
343             {
344               ri = -ai;
345               carry = false;
346             }
347           else
348             ri = ai;
349         }
350       else
351         ri = ~ai;
352
353       r->sig[i] = ri;
354     }
355 }
356
357 /* Compare significands.  Return tri-state vs zero.  */
358
359 static inline int
360 cmp_significands (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
361 {
362   int i;
363
364   for (i = SIGSZ - 1; i >= 0; --i)
365     {
366       unsigned long ai = a->sig[i];
367       unsigned long bi = b->sig[i];
368
369       if (ai > bi)
370         return 1;
371       if (ai < bi)
372         return -1;
373     }
374
375   return 0;
376 }
377
378 /* Return true if A is nonzero.  */
379
380 static inline int
381 cmp_significand_0 (const REAL_VALUE_TYPE *a)
382 {
383   int i;
384
385   for (i = SIGSZ - 1; i >= 0; --i)
386     if (a->sig[i])
387       return 1;
388
389   return 0;
390 }
391
392 /* Set bit N of the significand of R.  */
393
394 static inline void
395 set_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
396 {
397   r->sig[n / HOST_BITS_PER_LONG]
398     |= (unsigned long)1 << (n % HOST_BITS_PER_LONG);
399 }
400
401 /* Clear bit N of the significand of R.  */
402
403 static inline void
404 clear_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
405 {
406   r->sig[n / HOST_BITS_PER_LONG]
407     &= ~((unsigned long)1 << (n % HOST_BITS_PER_LONG));
408 }
409
410 /* Test bit N of the significand of R.  */
411
412 static inline bool
413 test_significand_bit (REAL_VALUE_TYPE *r, unsigned int n)
414 {
415   /* ??? Compiler bug here if we return this expression directly.
416      The conversion to bool strips the "&1" and we wind up testing
417      e.g. 2 != 0 -> true.  Seen in gcc version 3.2 20020520.  */
418   int t = (r->sig[n / HOST_BITS_PER_LONG] >> (n % HOST_BITS_PER_LONG)) & 1;
419   return t;
420 }
421
422 /* Clear bits 0..N-1 of the significand of R.  */
423
424 static void
425 clear_significand_below (REAL_VALUE_TYPE *r, unsigned int n)
426 {
427   int i, w = n / HOST_BITS_PER_LONG;
428
429   for (i = 0; i < w; ++i)
430     r->sig[i] = 0;
431
432   r->sig[w] &= ~(((unsigned long)1 << (n % HOST_BITS_PER_LONG)) - 1);
433 }
434
435 /* Divide the significands of A and B, placing the result in R.  Return
436    true if the division was inexact.  */
437
438 static inline bool
439 div_significands (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
440                   const REAL_VALUE_TYPE *b)
441 {
442   REAL_VALUE_TYPE u;
443   int i, bit = SIGNIFICAND_BITS - 1;
444   unsigned long msb, inexact;
445
446   u = *a;
447   memset (r->sig, 0, sizeof (r->sig));
448
449   msb = 0;
450   goto start;
451   do
452     {
453       msb = u.sig[SIGSZ-1] & SIG_MSB;
454       lshift_significand_1 (&u, &u);
455     start:
456       if (msb || cmp_significands (&u, b) >= 0)
457         {
458           sub_significands (&u, &u, b, 0);
459           set_significand_bit (r, bit);
460         }
461     }
462   while (--bit >= 0);
463
464   for (i = 0, inexact = 0; i < SIGSZ; i++)
465     inexact |= u.sig[i];
466
467   return inexact != 0;
468 }
469
470 /* Adjust the exponent and significand of R such that the most
471    significant bit is set.  We underflow to zero and overflow to
472    infinity here, without denormals.  (The intermediate representation
473    exponent is large enough to handle target denormals normalized.)  */
474
475 static void
476 normalize (REAL_VALUE_TYPE *r)
477 {
478   int shift = 0, exp;
479   int i, j;
480
481   if (r->decimal)
482     return;
483
484   /* Find the first word that is nonzero.  */
485   for (i = SIGSZ - 1; i >= 0; i--)
486     if (r->sig[i] == 0)
487       shift += HOST_BITS_PER_LONG;
488     else
489       break;
490
491   /* Zero significand flushes to zero.  */
492   if (i < 0)
493     {
494       r->cl = rvc_zero;
495       SET_REAL_EXP (r, 0);
496       return;
497     }
498
499   /* Find the first bit that is nonzero.  */
500   for (j = 0; ; j++)
501     if (r->sig[i] & ((unsigned long)1 << (HOST_BITS_PER_LONG - 1 - j)))
502       break;
503   shift += j;
504
505   if (shift > 0)
506     {
507       exp = REAL_EXP (r) - shift;
508       if (exp > MAX_EXP)
509         get_inf (r, r->sign);
510       else if (exp < -MAX_EXP)
511         get_zero (r, r->sign);
512       else
513         {
514           SET_REAL_EXP (r, exp);
515           lshift_significand (r, r, shift);
516         }
517     }
518 }
519 \f
520 /* Calculate R = A + (SUBTRACT_P ? -B : B).  Return true if the
521    result may be inexact due to a loss of precision.  */
522
523 static bool
524 do_add (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
525         const REAL_VALUE_TYPE *b, int subtract_p)
526 {
527   int dexp, sign, exp;
528   REAL_VALUE_TYPE t;
529   bool inexact = false;
530
531   /* Determine if we need to add or subtract.  */
532   sign = a->sign;
533   subtract_p = (sign ^ b->sign) ^ subtract_p;
534
535   switch (CLASS2 (a->cl, b->cl))
536     {
537     case CLASS2 (rvc_zero, rvc_zero):
538       /* -0 + -0 = -0, -0 - +0 = -0; all other cases yield +0.  */
539       get_zero (r, sign & !subtract_p);
540       return false;
541
542     case CLASS2 (rvc_zero, rvc_normal):
543     case CLASS2 (rvc_zero, rvc_inf):
544     case CLASS2 (rvc_zero, rvc_nan):
545       /* 0 + ANY = ANY.  */
546     case CLASS2 (rvc_normal, rvc_nan):
547     case CLASS2 (rvc_inf, rvc_nan):
548     case CLASS2 (rvc_nan, rvc_nan):
549       /* ANY + NaN = NaN.  */
550     case CLASS2 (rvc_normal, rvc_inf):
551       /* R + Inf = Inf.  */
552       *r = *b;
553       r->sign = sign ^ subtract_p;
554       return false;
555
556     case CLASS2 (rvc_normal, rvc_zero):
557     case CLASS2 (rvc_inf, rvc_zero):
558     case CLASS2 (rvc_nan, rvc_zero):
559       /* ANY + 0 = ANY.  */
560     case CLASS2 (rvc_nan, rvc_normal):
561     case CLASS2 (rvc_nan, rvc_inf):
562       /* NaN + ANY = NaN.  */
563     case CLASS2 (rvc_inf, rvc_normal):
564       /* Inf + R = Inf.  */
565       *r = *a;
566       return false;
567
568     case CLASS2 (rvc_inf, rvc_inf):
569       if (subtract_p)
570         /* Inf - Inf = NaN.  */
571         get_canonical_qnan (r, 0);
572       else
573         /* Inf + Inf = Inf.  */
574         *r = *a;
575       return false;
576
577     case CLASS2 (rvc_normal, rvc_normal):
578       break;
579
580     default:
581       gcc_unreachable ();
582     }
583
584   /* Swap the arguments such that A has the larger exponent.  */
585   dexp = REAL_EXP (a) - REAL_EXP (b);
586   if (dexp < 0)
587     {
588       const REAL_VALUE_TYPE *t;
589       t = a, a = b, b = t;
590       dexp = -dexp;
591       sign ^= subtract_p;
592     }
593   exp = REAL_EXP (a);
594
595   /* If the exponents are not identical, we need to shift the
596      significand of B down.  */
597   if (dexp > 0)
598     {
599       /* If the exponents are too far apart, the significands
600          do not overlap, which makes the subtraction a noop.  */
601       if (dexp >= SIGNIFICAND_BITS)
602         {
603           *r = *a;
604           r->sign = sign;
605           return true;
606         }
607
608       inexact |= sticky_rshift_significand (&t, b, dexp);
609       b = &t;
610     }
611
612   if (subtract_p)
613     {
614       if (sub_significands (r, a, b, inexact))
615         {
616           /* We got a borrow out of the subtraction.  That means that
617              A and B had the same exponent, and B had the larger
618              significand.  We need to swap the sign and negate the
619              significand.  */
620           sign ^= 1;
621           neg_significand (r, r);
622         }
623     }
624   else
625     {
626       if (add_significands (r, a, b))
627         {
628           /* We got carry out of the addition.  This means we need to
629              shift the significand back down one bit and increase the
630              exponent.  */
631           inexact |= sticky_rshift_significand (r, r, 1);
632           r->sig[SIGSZ-1] |= SIG_MSB;
633           if (++exp > MAX_EXP)
634             {
635               get_inf (r, sign);
636               return true;
637             }
638         }
639     }
640
641   r->cl = rvc_normal;
642   r->sign = sign;
643   SET_REAL_EXP (r, exp);
644   /* Zero out the remaining fields.  */
645   r->signalling = 0;
646   r->canonical = 0;
647   r->decimal = 0;
648
649   /* Re-normalize the result.  */
650   normalize (r);
651
652   /* Special case: if the subtraction results in zero, the result
653      is positive.  */
654   if (r->cl == rvc_zero)
655     r->sign = 0;
656   else
657     r->sig[0] |= inexact;
658
659   return inexact;
660 }
661
662 /* Calculate R = A * B.  Return true if the result may be inexact.  */
663
664 static bool
665 do_multiply (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
666              const REAL_VALUE_TYPE *b)
667 {
668   REAL_VALUE_TYPE u, t, *rr;
669   unsigned int i, j, k;
670   int sign = a->sign ^ b->sign;
671   bool inexact = false;
672
673   switch (CLASS2 (a->cl, b->cl))
674     {
675     case CLASS2 (rvc_zero, rvc_zero):
676     case CLASS2 (rvc_zero, rvc_normal):
677     case CLASS2 (rvc_normal, rvc_zero):
678       /* +-0 * ANY = 0 with appropriate sign.  */
679       get_zero (r, sign);
680       return false;
681
682     case CLASS2 (rvc_zero, rvc_nan):
683     case CLASS2 (rvc_normal, rvc_nan):
684     case CLASS2 (rvc_inf, rvc_nan):
685     case CLASS2 (rvc_nan, rvc_nan):
686       /* ANY * NaN = NaN.  */
687       *r = *b;
688       r->sign = sign;
689       return false;
690
691     case CLASS2 (rvc_nan, rvc_zero):
692     case CLASS2 (rvc_nan, rvc_normal):
693     case CLASS2 (rvc_nan, rvc_inf):
694       /* NaN * ANY = NaN.  */
695       *r = *a;
696       r->sign = sign;
697       return false;
698
699     case CLASS2 (rvc_zero, rvc_inf):
700     case CLASS2 (rvc_inf, rvc_zero):
701       /* 0 * Inf = NaN */
702       get_canonical_qnan (r, sign);
703       return false;
704
705     case CLASS2 (rvc_inf, rvc_inf):
706     case CLASS2 (rvc_normal, rvc_inf):
707     case CLASS2 (rvc_inf, rvc_normal):
708       /* Inf * Inf = Inf, R * Inf = Inf */
709       get_inf (r, sign);
710       return false;
711
712     case CLASS2 (rvc_normal, rvc_normal):
713       break;
714
715     default:
716       gcc_unreachable ();
717     }
718
719   if (r == a || r == b)
720     rr = &t;
721   else
722     rr = r;
723   get_zero (rr, 0);
724
725   /* Collect all the partial products.  Since we don't have sure access
726      to a widening multiply, we split each long into two half-words.
727
728      Consider the long-hand form of a four half-word multiplication:
729
730                  A  B  C  D
731               *  E  F  G  H
732              --------------
733                 DE DF DG DH
734              CE CF CG CH
735           BE BF BG BH
736        AE AF AG AH
737
738      We construct partial products of the widened half-word products
739      that are known to not overlap, e.g. DF+DH.  Each such partial
740      product is given its proper exponent, which allows us to sum them
741      and obtain the finished product.  */
742
743   for (i = 0; i < SIGSZ * 2; ++i)
744     {
745       unsigned long ai = a->sig[i / 2];
746       if (i & 1)
747         ai >>= HOST_BITS_PER_LONG / 2;
748       else
749         ai &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
750
751       if (ai == 0)
752         continue;
753
754       for (j = 0; j < 2; ++j)
755         {
756           int exp = (REAL_EXP (a) - (2*SIGSZ-1-i)*(HOST_BITS_PER_LONG/2)
757                      + (REAL_EXP (b) - (1-j)*(HOST_BITS_PER_LONG/2)));
758
759           if (exp > MAX_EXP)
760             {
761               get_inf (r, sign);
762               return true;
763             }
764           if (exp < -MAX_EXP)
765             {
766               /* Would underflow to zero, which we shouldn't bother adding.  */
767               inexact = true;
768               continue;
769             }
770
771           memset (&u, 0, sizeof (u));
772           u.cl = rvc_normal;
773           SET_REAL_EXP (&u, exp);
774
775           for (k = j; k < SIGSZ * 2; k += 2)
776             {
777               unsigned long bi = b->sig[k / 2];
778               if (k & 1)
779                 bi >>= HOST_BITS_PER_LONG / 2;
780               else
781                 bi &= ((unsigned long)1 << (HOST_BITS_PER_LONG / 2)) - 1;
782
783               u.sig[k / 2] = ai * bi;
784             }
785
786           normalize (&u);
787           inexact |= do_add (rr, rr, &u, 0);
788         }
789     }
790
791   rr->sign = sign;
792   if (rr != r)
793     *r = t;
794
795   return inexact;
796 }
797
798 /* Calculate R = A / B.  Return true if the result may be inexact.  */
799
800 static bool
801 do_divide (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a,
802            const REAL_VALUE_TYPE *b)
803 {
804   int exp, sign = a->sign ^ b->sign;
805   REAL_VALUE_TYPE t, *rr;
806   bool inexact;
807
808   switch (CLASS2 (a->cl, b->cl))
809     {
810     case CLASS2 (rvc_zero, rvc_zero):
811       /* 0 / 0 = NaN.  */
812     case CLASS2 (rvc_inf, rvc_inf):
813       /* Inf / Inf = NaN.  */
814       get_canonical_qnan (r, sign);
815       return false;
816
817     case CLASS2 (rvc_zero, rvc_normal):
818     case CLASS2 (rvc_zero, rvc_inf):
819       /* 0 / ANY = 0.  */
820     case CLASS2 (rvc_normal, rvc_inf):
821       /* R / Inf = 0.  */
822       get_zero (r, sign);
823       return false;
824
825     case CLASS2 (rvc_normal, rvc_zero):
826       /* R / 0 = Inf.  */
827     case CLASS2 (rvc_inf, rvc_zero):
828       /* Inf / 0 = Inf.  */
829       get_inf (r, sign);
830       return false;
831
832     case CLASS2 (rvc_zero, rvc_nan):
833     case CLASS2 (rvc_normal, rvc_nan):
834     case CLASS2 (rvc_inf, rvc_nan):
835     case CLASS2 (rvc_nan, rvc_nan):
836       /* ANY / NaN = NaN.  */
837       *r = *b;
838       r->sign = sign;
839       return false;
840
841     case CLASS2 (rvc_nan, rvc_zero):
842     case CLASS2 (rvc_nan, rvc_normal):
843     case CLASS2 (rvc_nan, rvc_inf):
844       /* NaN / ANY = NaN.  */
845       *r = *a;
846       r->sign = sign;
847       return false;
848
849     case CLASS2 (rvc_inf, rvc_normal):
850       /* Inf / R = Inf.  */
851       get_inf (r, sign);
852       return false;
853
854     case CLASS2 (rvc_normal, rvc_normal):
855       break;
856
857     default:
858       gcc_unreachable ();
859     }
860
861   if (r == a || r == b)
862     rr = &t;
863   else
864     rr = r;
865
866   /* Make sure all fields in the result are initialized.  */
867   get_zero (rr, 0);
868   rr->cl = rvc_normal;
869   rr->sign = sign;
870
871   exp = REAL_EXP (a) - REAL_EXP (b) + 1;
872   if (exp > MAX_EXP)
873     {
874       get_inf (r, sign);
875       return true;
876     }
877   if (exp < -MAX_EXP)
878     {
879       get_zero (r, sign);
880       return true;
881     }
882   SET_REAL_EXP (rr, exp);
883
884   inexact = div_significands (rr, a, b);
885
886   /* Re-normalize the result.  */
887   normalize (rr);
888   rr->sig[0] |= inexact;
889
890   if (rr != r)
891     *r = t;
892
893   return inexact;
894 }
895
896 /* Return a tri-state comparison of A vs B.  Return NAN_RESULT if
897    one of the two operands is a NaN.  */
898
899 static int
900 do_compare (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b,
901             int nan_result)
902 {
903   int ret;
904
905   switch (CLASS2 (a->cl, b->cl))
906     {
907     case CLASS2 (rvc_zero, rvc_zero):
908       /* Sign of zero doesn't matter for compares.  */
909       return 0;
910
911     case CLASS2 (rvc_normal, rvc_zero):
912       /* Decimal float zero is special and uses rvc_normal, not rvc_zero.  */
913       if (a->decimal)
914         return decimal_do_compare (a, b, nan_result);
915       /* Fall through.  */
916     case CLASS2 (rvc_inf, rvc_zero):
917     case CLASS2 (rvc_inf, rvc_normal):
918       return (a->sign ? -1 : 1);
919
920     case CLASS2 (rvc_inf, rvc_inf):
921       return -a->sign - -b->sign;
922
923     case CLASS2 (rvc_zero, rvc_normal):
924       /* Decimal float zero is special and uses rvc_normal, not rvc_zero.  */
925       if (b->decimal)
926         return decimal_do_compare (a, b, nan_result);
927       /* Fall through.  */
928     case CLASS2 (rvc_zero, rvc_inf):
929     case CLASS2 (rvc_normal, rvc_inf):
930       return (b->sign ? 1 : -1);
931
932     case CLASS2 (rvc_zero, rvc_nan):
933     case CLASS2 (rvc_normal, rvc_nan):
934     case CLASS2 (rvc_inf, rvc_nan):
935     case CLASS2 (rvc_nan, rvc_nan):
936     case CLASS2 (rvc_nan, rvc_zero):
937     case CLASS2 (rvc_nan, rvc_normal):
938     case CLASS2 (rvc_nan, rvc_inf):
939       return nan_result;
940
941     case CLASS2 (rvc_normal, rvc_normal):
942       break;
943
944     default:
945       gcc_unreachable ();
946     }
947
948   if (a->sign != b->sign)
949     return -a->sign - -b->sign;
950
951   if (a->decimal || b->decimal)
952     return decimal_do_compare (a, b, nan_result);
953
954   if (REAL_EXP (a) > REAL_EXP (b))
955     ret = 1;
956   else if (REAL_EXP (a) < REAL_EXP (b))
957     ret = -1;
958   else
959     ret = cmp_significands (a, b);
960
961   return (a->sign ? -ret : ret);
962 }
963
964 /* Return A truncated to an integral value toward zero.  */
965
966 static void
967 do_fix_trunc (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *a)
968 {
969   *r = *a;
970
971   switch (r->cl)
972     {
973     case rvc_zero:
974     case rvc_inf:
975     case rvc_nan:
976       break;
977
978     case rvc_normal:
979       if (r->decimal)
980         {
981           decimal_do_fix_trunc (r, a);
982           return;
983         }
984       if (REAL_EXP (r) <= 0)
985         get_zero (r, r->sign);
986       else if (REAL_EXP (r) < SIGNIFICAND_BITS)
987         clear_significand_below (r, SIGNIFICAND_BITS - REAL_EXP (r));
988       break;
989
990     default:
991       gcc_unreachable ();
992     }
993 }
994
995 /* Perform the binary or unary operation described by CODE.
996    For a unary operation, leave OP1 NULL.  This function returns
997    true if the result may be inexact due to loss of precision.  */
998
999 bool
1000 real_arithmetic (REAL_VALUE_TYPE *r, int icode, const REAL_VALUE_TYPE *op0,
1001                  const REAL_VALUE_TYPE *op1)
1002 {
1003   enum tree_code code = (enum tree_code) icode;
1004
1005   if (op0->decimal || (op1 && op1->decimal))
1006     return decimal_real_arithmetic (r, code, op0, op1);
1007
1008   switch (code)
1009     {
1010     case PLUS_EXPR:
1011       /* Clear any padding areas in *r if it isn't equal to one of the
1012          operands so that we can later do bitwise comparisons later on.  */
1013       if (r != op0 && r != op1)
1014         memset (r, '\0', sizeof (*r));
1015       return do_add (r, op0, op1, 0);
1016
1017     case MINUS_EXPR:
1018       if (r != op0 && r != op1)
1019         memset (r, '\0', sizeof (*r));
1020       return do_add (r, op0, op1, 1);
1021
1022     case MULT_EXPR:
1023       if (r != op0 && r != op1)
1024         memset (r, '\0', sizeof (*r));
1025       return do_multiply (r, op0, op1);
1026
1027     case RDIV_EXPR:
1028       if (r != op0 && r != op1)
1029         memset (r, '\0', sizeof (*r));
1030       return do_divide (r, op0, op1);
1031
1032     case MIN_EXPR:
1033       if (op1->cl == rvc_nan)
1034         *r = *op1;
1035       else if (do_compare (op0, op1, -1) < 0)
1036         *r = *op0;
1037       else
1038         *r = *op1;
1039       break;
1040
1041     case MAX_EXPR:
1042       if (op1->cl == rvc_nan)
1043         *r = *op1;
1044       else if (do_compare (op0, op1, 1) < 0)
1045         *r = *op1;
1046       else
1047         *r = *op0;
1048       break;
1049
1050     case NEGATE_EXPR:
1051       *r = *op0;
1052       r->sign ^= 1;
1053       break;
1054
1055     case ABS_EXPR:
1056       *r = *op0;
1057       r->sign = 0;
1058       break;
1059
1060     case FIX_TRUNC_EXPR:
1061       do_fix_trunc (r, op0);
1062       break;
1063
1064     default:
1065       gcc_unreachable ();
1066     }
1067   return false;
1068 }
1069
1070 REAL_VALUE_TYPE
1071 real_value_negate (const REAL_VALUE_TYPE *op0)
1072 {
1073   REAL_VALUE_TYPE r;
1074   real_arithmetic (&r, NEGATE_EXPR, op0, NULL);
1075   return r;
1076 }
1077
1078 REAL_VALUE_TYPE
1079 real_value_abs (const REAL_VALUE_TYPE *op0)
1080 {
1081   REAL_VALUE_TYPE r;
1082   real_arithmetic (&r, ABS_EXPR, op0, NULL);
1083   return r;
1084 }
1085
1086 bool
1087 real_compare (int icode, const REAL_VALUE_TYPE *op0,
1088               const REAL_VALUE_TYPE *op1)
1089 {
1090   enum tree_code code = (enum tree_code) icode;
1091
1092   switch (code)
1093     {
1094     case LT_EXPR:
1095       return do_compare (op0, op1, 1) < 0;
1096     case LE_EXPR:
1097       return do_compare (op0, op1, 1) <= 0;
1098     case GT_EXPR:
1099       return do_compare (op0, op1, -1) > 0;
1100     case GE_EXPR:
1101       return do_compare (op0, op1, -1) >= 0;
1102     case EQ_EXPR:
1103       return do_compare (op0, op1, -1) == 0;
1104     case NE_EXPR:
1105       return do_compare (op0, op1, -1) != 0;
1106     case UNORDERED_EXPR:
1107       return op0->cl == rvc_nan || op1->cl == rvc_nan;
1108     case ORDERED_EXPR:
1109       return op0->cl != rvc_nan && op1->cl != rvc_nan;
1110     case UNLT_EXPR:
1111       return do_compare (op0, op1, -1) < 0;
1112     case UNLE_EXPR:
1113       return do_compare (op0, op1, -1) <= 0;
1114     case UNGT_EXPR:
1115       return do_compare (op0, op1, 1) > 0;
1116     case UNGE_EXPR:
1117       return do_compare (op0, op1, 1) >= 0;
1118     case UNEQ_EXPR:
1119       return do_compare (op0, op1, 0) == 0;
1120     case LTGT_EXPR:
1121       return do_compare (op0, op1, 0) != 0;
1122
1123     default:
1124       gcc_unreachable ();
1125     }
1126 }
1127
1128 /* Return floor log2(R).  */
1129
1130 int
1131 real_exponent (const REAL_VALUE_TYPE *r)
1132 {
1133   switch (r->cl)
1134     {
1135     case rvc_zero:
1136       return 0;
1137     case rvc_inf:
1138     case rvc_nan:
1139       return (unsigned int)-1 >> 1;
1140     case rvc_normal:
1141       return REAL_EXP (r);
1142     default:
1143       gcc_unreachable ();
1144     }
1145 }
1146
1147 /* R = OP0 * 2**EXP.  */
1148
1149 void
1150 real_ldexp (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *op0, int exp)
1151 {
1152   *r = *op0;
1153   switch (r->cl)
1154     {
1155     case rvc_zero:
1156     case rvc_inf:
1157     case rvc_nan:
1158       break;
1159
1160     case rvc_normal:
1161       exp += REAL_EXP (op0);
1162       if (exp > MAX_EXP)
1163         get_inf (r, r->sign);
1164       else if (exp < -MAX_EXP)
1165         get_zero (r, r->sign);
1166       else
1167         SET_REAL_EXP (r, exp);
1168       break;
1169
1170     default:
1171       gcc_unreachable ();
1172     }
1173 }
1174
1175 /* Determine whether a floating-point value X is infinite.  */
1176
1177 bool
1178 real_isinf (const REAL_VALUE_TYPE *r)
1179 {
1180   return (r->cl == rvc_inf);
1181 }
1182
1183 /* Determine whether a floating-point value X is a NaN.  */
1184
1185 bool
1186 real_isnan (const REAL_VALUE_TYPE *r)
1187 {
1188   return (r->cl == rvc_nan);
1189 }
1190
1191 /* Determine whether a floating-point value X is finite.  */
1192
1193 bool
1194 real_isfinite (const REAL_VALUE_TYPE *r)
1195 {
1196   return (r->cl != rvc_nan) && (r->cl != rvc_inf);
1197 }
1198
1199 /* Determine whether a floating-point value X is negative.  */
1200
1201 bool
1202 real_isneg (const REAL_VALUE_TYPE *r)
1203 {
1204   return r->sign;
1205 }
1206
1207 /* Determine whether a floating-point value X is minus zero.  */
1208
1209 bool
1210 real_isnegzero (const REAL_VALUE_TYPE *r)
1211 {
1212   return r->sign && r->cl == rvc_zero;
1213 }
1214
1215 /* Compare two floating-point objects for bitwise identity.  */
1216
1217 bool
1218 real_identical (const REAL_VALUE_TYPE *a, const REAL_VALUE_TYPE *b)
1219 {
1220   int i;
1221
1222   if (a->cl != b->cl)
1223     return false;
1224   if (a->sign != b->sign)
1225     return false;
1226
1227   switch (a->cl)
1228     {
1229     case rvc_zero:
1230     case rvc_inf:
1231       return true;
1232
1233     case rvc_normal:
1234       if (a->decimal != b->decimal)
1235         return false;
1236       if (REAL_EXP (a) != REAL_EXP (b))
1237         return false;
1238       break;
1239
1240     case rvc_nan:
1241       if (a->signalling != b->signalling)
1242         return false;
1243       /* The significand is ignored for canonical NaNs.  */
1244       if (a->canonical || b->canonical)
1245         return a->canonical == b->canonical;
1246       break;
1247
1248     default:
1249       gcc_unreachable ();
1250     }
1251
1252   for (i = 0; i < SIGSZ; ++i)
1253     if (a->sig[i] != b->sig[i])
1254       return false;
1255
1256   return true;
1257 }
1258
1259 /* Try to change R into its exact multiplicative inverse in machine
1260    mode MODE.  Return true if successful.  */
1261
1262 bool
1263 exact_real_inverse (enum machine_mode mode, REAL_VALUE_TYPE *r)
1264 {
1265   const REAL_VALUE_TYPE *one = real_digit (1);
1266   REAL_VALUE_TYPE u;
1267   int i;
1268
1269   if (r->cl != rvc_normal)
1270     return false;
1271
1272   /* Check for a power of two: all significand bits zero except the MSB.  */
1273   for (i = 0; i < SIGSZ-1; ++i)
1274     if (r->sig[i] != 0)
1275       return false;
1276   if (r->sig[SIGSZ-1] != SIG_MSB)
1277     return false;
1278
1279   /* Find the inverse and truncate to the required mode.  */
1280   do_divide (&u, one, r);
1281   real_convert (&u, mode, &u);
1282
1283   /* The rounding may have overflowed.  */
1284   if (u.cl != rvc_normal)
1285     return false;
1286   for (i = 0; i < SIGSZ-1; ++i)
1287     if (u.sig[i] != 0)
1288       return false;
1289   if (u.sig[SIGSZ-1] != SIG_MSB)
1290     return false;
1291
1292   *r = u;
1293   return true;
1294 }
1295
1296 /* Return true if arithmetic on values in IMODE that were promoted
1297    from values in TMODE is equivalent to direct arithmetic on values
1298    in TMODE.  */
1299
1300 bool
1301 real_can_shorten_arithmetic (enum machine_mode imode, enum machine_mode tmode)
1302 {
1303   const struct real_format *tfmt, *ifmt;
1304   tfmt = REAL_MODE_FORMAT (tmode);
1305   ifmt = REAL_MODE_FORMAT (imode);
1306   /* These conditions are conservative rather than trying to catch the
1307      exact boundary conditions; the main case to allow is IEEE float
1308      and double.  */
1309   return (ifmt->b == tfmt->b
1310           && ifmt->p > 2 * tfmt->p
1311           && ifmt->emin < 2 * tfmt->emin - tfmt->p - 2
1312           && ifmt->emin < tfmt->emin - tfmt->emax - tfmt->p - 2
1313           && ifmt->emax > 2 * tfmt->emax + 2
1314           && ifmt->emax > tfmt->emax - tfmt->emin + tfmt->p + 2
1315           && ifmt->round_towards_zero == tfmt->round_towards_zero
1316           && (ifmt->has_sign_dependent_rounding
1317               == tfmt->has_sign_dependent_rounding)
1318           && ifmt->has_nans >= tfmt->has_nans
1319           && ifmt->has_inf >= tfmt->has_inf
1320           && ifmt->has_signed_zero >= tfmt->has_signed_zero
1321           && !MODE_COMPOSITE_P (tmode)
1322           && !MODE_COMPOSITE_P (imode));
1323 }
1324 \f
1325 /* Render R as an integer.  */
1326
1327 HOST_WIDE_INT
1328 real_to_integer (const REAL_VALUE_TYPE *r)
1329 {
1330   unsigned HOST_WIDE_INT i;
1331
1332   switch (r->cl)
1333     {
1334     case rvc_zero:
1335     underflow:
1336       return 0;
1337
1338     case rvc_inf:
1339     case rvc_nan:
1340     overflow:
1341       i = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1342       if (!r->sign)
1343         i--;
1344       return i;
1345
1346     case rvc_normal:
1347       if (r->decimal)
1348         return decimal_real_to_integer (r);
1349
1350       if (REAL_EXP (r) <= 0)
1351         goto underflow;
1352       /* Only force overflow for unsigned overflow.  Signed overflow is
1353          undefined, so it doesn't matter what we return, and some callers
1354          expect to be able to use this routine for both signed and
1355          unsigned conversions.  */
1356       if (REAL_EXP (r) > HOST_BITS_PER_WIDE_INT)
1357         goto overflow;
1358
1359       if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1360         i = r->sig[SIGSZ-1];
1361       else
1362         {
1363           gcc_assert (HOST_BITS_PER_WIDE_INT == 2 * HOST_BITS_PER_LONG);
1364           i = r->sig[SIGSZ-1];
1365           i = i << (HOST_BITS_PER_LONG - 1) << 1;
1366           i |= r->sig[SIGSZ-2];
1367         }
1368
1369       i >>= HOST_BITS_PER_WIDE_INT - REAL_EXP (r);
1370
1371       if (r->sign)
1372         i = -i;
1373       return i;
1374
1375     default:
1376       gcc_unreachable ();
1377     }
1378 }
1379
1380 /* Likewise, but to an integer pair, HI+LOW.  */
1381
1382 void
1383 real_to_integer2 (HOST_WIDE_INT *plow, HOST_WIDE_INT *phigh,
1384                   const REAL_VALUE_TYPE *r)
1385 {
1386   REAL_VALUE_TYPE t;
1387   HOST_WIDE_INT low, high;
1388   int exp;
1389
1390   switch (r->cl)
1391     {
1392     case rvc_zero:
1393     underflow:
1394       low = high = 0;
1395       break;
1396
1397     case rvc_inf:
1398     case rvc_nan:
1399     overflow:
1400       high = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
1401       if (r->sign)
1402         low = 0;
1403       else
1404         {
1405           high--;
1406           low = -1;
1407         }
1408       break;
1409
1410     case rvc_normal:
1411       if (r->decimal)
1412         {
1413           decimal_real_to_integer2 (plow, phigh, r);
1414           return;
1415         }
1416
1417       exp = REAL_EXP (r);
1418       if (exp <= 0)
1419         goto underflow;
1420       /* Only force overflow for unsigned overflow.  Signed overflow is
1421          undefined, so it doesn't matter what we return, and some callers
1422          expect to be able to use this routine for both signed and
1423          unsigned conversions.  */
1424       if (exp > HOST_BITS_PER_DOUBLE_INT)
1425         goto overflow;
1426
1427       rshift_significand (&t, r, HOST_BITS_PER_DOUBLE_INT - exp);
1428       if (HOST_BITS_PER_WIDE_INT == HOST_BITS_PER_LONG)
1429         {
1430           high = t.sig[SIGSZ-1];
1431           low = t.sig[SIGSZ-2];
1432         }
1433       else
1434         {
1435           gcc_assert (HOST_BITS_PER_WIDE_INT == 2*HOST_BITS_PER_LONG);
1436           high = t.sig[SIGSZ-1];
1437           high = high << (HOST_BITS_PER_LONG - 1) << 1;
1438           high |= t.sig[SIGSZ-2];
1439
1440           low = t.sig[SIGSZ-3];
1441           low = low << (HOST_BITS_PER_LONG - 1) << 1;
1442           low |= t.sig[SIGSZ-4];
1443         }
1444
1445       if (r->sign)
1446         {
1447           if (low == 0)
1448             high = -high;
1449           else
1450             low = -low, high = ~high;
1451         }
1452       break;
1453
1454     default:
1455       gcc_unreachable ();
1456     }
1457
1458   *plow = low;
1459   *phigh = high;
1460 }
1461
1462 /* A subroutine of real_to_decimal.  Compute the quotient and remainder
1463    of NUM / DEN.  Return the quotient and place the remainder in NUM.
1464    It is expected that NUM / DEN are close enough that the quotient is
1465    small.  */
1466
1467 static unsigned long
1468 rtd_divmod (REAL_VALUE_TYPE *num, REAL_VALUE_TYPE *den)
1469 {
1470   unsigned long q, msb;
1471   int expn = REAL_EXP (num), expd = REAL_EXP (den);
1472
1473   if (expn < expd)
1474     return 0;
1475
1476   q = msb = 0;
1477   goto start;
1478   do
1479     {
1480       msb = num->sig[SIGSZ-1] & SIG_MSB;
1481       q <<= 1;
1482       lshift_significand_1 (num, num);
1483     start:
1484       if (msb || cmp_significands (num, den) >= 0)
1485         {
1486           sub_significands (num, num, den, 0);
1487           q |= 1;
1488         }
1489     }
1490   while (--expn >= expd);
1491
1492   SET_REAL_EXP (num, expd);
1493   normalize (num);
1494
1495   return q;
1496 }
1497
1498 /* Render R as a decimal floating point constant.  Emit DIGITS significant
1499    digits in the result, bounded by BUF_SIZE.  If DIGITS is 0, choose the
1500    maximum for the representation.  If CROP_TRAILING_ZEROS, strip trailing
1501    zeros.  If MODE is VOIDmode, round to nearest value.  Otherwise, round
1502    to a string that, when parsed back in mode MODE, yields the same value.  */
1503
1504 #define M_LOG10_2       0.30102999566398119521
1505
1506 void
1507 real_to_decimal_for_mode (char *str, const REAL_VALUE_TYPE *r_orig,
1508                           size_t buf_size, size_t digits,
1509                           int crop_trailing_zeros, enum machine_mode mode)
1510 {
1511   const struct real_format *fmt = NULL;
1512   const REAL_VALUE_TYPE *one, *ten;
1513   REAL_VALUE_TYPE r, pten, u, v;
1514   int dec_exp, cmp_one, digit;
1515   size_t max_digits;
1516   char *p, *first, *last;
1517   bool sign;
1518   bool round_up;
1519
1520   if (mode != VOIDmode)
1521    {
1522      fmt = REAL_MODE_FORMAT (mode);
1523      gcc_assert (fmt);
1524    }
1525
1526   r = *r_orig;
1527   switch (r.cl)
1528     {
1529     case rvc_zero:
1530       strcpy (str, (r.sign ? "-0.0" : "0.0"));
1531       return;
1532     case rvc_normal:
1533       break;
1534     case rvc_inf:
1535       strcpy (str, (r.sign ? "-Inf" : "+Inf"));
1536       return;
1537     case rvc_nan:
1538       /* ??? Print the significand as well, if not canonical?  */
1539       sprintf (str, "%c%cNaN", (r_orig->sign ? '-' : '+'),
1540                (r_orig->signalling ? 'S' : 'Q'));
1541       return;
1542     default:
1543       gcc_unreachable ();
1544     }
1545
1546   if (r.decimal)
1547     {
1548       decimal_real_to_decimal (str, &r, buf_size, digits, crop_trailing_zeros);
1549       return;
1550     }
1551
1552   /* Bound the number of digits printed by the size of the representation.  */
1553   max_digits = SIGNIFICAND_BITS * M_LOG10_2;
1554   if (digits == 0 || digits > max_digits)
1555     digits = max_digits;
1556
1557   /* Estimate the decimal exponent, and compute the length of the string it
1558      will print as.  Be conservative and add one to account for possible
1559      overflow or rounding error.  */
1560   dec_exp = REAL_EXP (&r) * M_LOG10_2;
1561   for (max_digits = 1; dec_exp ; max_digits++)
1562     dec_exp /= 10;
1563
1564   /* Bound the number of digits printed by the size of the output buffer.  */
1565   max_digits = buf_size - 1 - 1 - 2 - max_digits - 1;
1566   gcc_assert (max_digits <= buf_size);
1567   if (digits > max_digits)
1568     digits = max_digits;
1569
1570   one = real_digit (1);
1571   ten = ten_to_ptwo (0);
1572
1573   sign = r.sign;
1574   r.sign = 0;
1575
1576   dec_exp = 0;
1577   pten = *one;
1578
1579   cmp_one = do_compare (&r, one, 0);
1580   if (cmp_one > 0)
1581     {
1582       int m;
1583
1584       /* Number is greater than one.  Convert significand to an integer
1585          and strip trailing decimal zeros.  */
1586
1587       u = r;
1588       SET_REAL_EXP (&u, SIGNIFICAND_BITS - 1);
1589
1590       /* Largest M, such that 10**2**M fits within SIGNIFICAND_BITS.  */
1591       m = floor_log2 (max_digits);
1592
1593       /* Iterate over the bits of the possible powers of 10 that might
1594          be present in U and eliminate them.  That is, if we find that
1595          10**2**M divides U evenly, keep the division and increase
1596          DEC_EXP by 2**M.  */
1597       do
1598         {
1599           REAL_VALUE_TYPE t;
1600
1601           do_divide (&t, &u, ten_to_ptwo (m));
1602           do_fix_trunc (&v, &t);
1603           if (cmp_significands (&v, &t) == 0)
1604             {
1605               u = t;
1606               dec_exp += 1 << m;
1607             }
1608         }
1609       while (--m >= 0);
1610
1611       /* Revert the scaling to integer that we performed earlier.  */
1612       SET_REAL_EXP (&u, REAL_EXP (&u) + REAL_EXP (&r)
1613                     - (SIGNIFICAND_BITS - 1));
1614       r = u;
1615
1616       /* Find power of 10.  Do this by dividing out 10**2**M when
1617          this is larger than the current remainder.  Fill PTEN with
1618          the power of 10 that we compute.  */
1619       if (REAL_EXP (&r) > 0)
1620         {
1621           m = floor_log2 ((int)(REAL_EXP (&r) * M_LOG10_2)) + 1;
1622           do
1623             {
1624               const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1625               if (do_compare (&u, ptentwo, 0) >= 0)
1626                 {
1627                   do_divide (&u, &u, ptentwo);
1628                   do_multiply (&pten, &pten, ptentwo);
1629                   dec_exp += 1 << m;
1630                 }
1631             }
1632           while (--m >= 0);
1633         }
1634       else
1635         /* We managed to divide off enough tens in the above reduction
1636            loop that we've now got a negative exponent.  Fall into the
1637            less-than-one code to compute the proper value for PTEN.  */
1638         cmp_one = -1;
1639     }
1640   if (cmp_one < 0)
1641     {
1642       int m;
1643
1644       /* Number is less than one.  Pad significand with leading
1645          decimal zeros.  */
1646
1647       v = r;
1648       while (1)
1649         {
1650           /* Stop if we'd shift bits off the bottom.  */
1651           if (v.sig[0] & 7)
1652             break;
1653
1654           do_multiply (&u, &v, ten);
1655
1656           /* Stop if we're now >= 1.  */
1657           if (REAL_EXP (&u) > 0)
1658             break;
1659
1660           v = u;
1661           dec_exp -= 1;
1662         }
1663       r = v;
1664
1665       /* Find power of 10.  Do this by multiplying in P=10**2**M when
1666          the current remainder is smaller than 1/P.  Fill PTEN with the
1667          power of 10 that we compute.  */
1668       m = floor_log2 ((int)(-REAL_EXP (&r) * M_LOG10_2)) + 1;
1669       do
1670         {
1671           const REAL_VALUE_TYPE *ptentwo = ten_to_ptwo (m);
1672           const REAL_VALUE_TYPE *ptenmtwo = ten_to_mptwo (m);
1673
1674           if (do_compare (&v, ptenmtwo, 0) <= 0)
1675             {
1676               do_multiply (&v, &v, ptentwo);
1677               do_multiply (&pten, &pten, ptentwo);
1678               dec_exp -= 1 << m;
1679             }
1680         }
1681       while (--m >= 0);
1682
1683       /* Invert the positive power of 10 that we've collected so far.  */
1684       do_divide (&pten, one, &pten);
1685     }
1686
1687   p = str;
1688   if (sign)
1689     *p++ = '-';
1690   first = p++;
1691
1692   /* At this point, PTEN should contain the nearest power of 10 smaller
1693      than R, such that this division produces the first digit.
1694
1695      Using a divide-step primitive that returns the complete integral
1696      remainder avoids the rounding error that would be produced if
1697      we were to use do_divide here and then simply multiply by 10 for
1698      each subsequent digit.  */
1699
1700   digit = rtd_divmod (&r, &pten);
1701
1702   /* Be prepared for error in that division via underflow ...  */
1703   if (digit == 0 && cmp_significand_0 (&r))
1704     {
1705       /* Multiply by 10 and try again.  */
1706       do_multiply (&r, &r, ten);
1707       digit = rtd_divmod (&r, &pten);
1708       dec_exp -= 1;
1709       gcc_assert (digit != 0);
1710     }
1711
1712   /* ... or overflow.  */
1713   if (digit == 10)
1714     {
1715       *p++ = '1';
1716       if (--digits > 0)
1717         *p++ = '0';
1718       dec_exp += 1;
1719     }
1720   else
1721     {
1722       gcc_assert (digit <= 10);
1723       *p++ = digit + '0';
1724     }
1725
1726   /* Generate subsequent digits.  */
1727   while (--digits > 0)
1728     {
1729       do_multiply (&r, &r, ten);
1730       digit = rtd_divmod (&r, &pten);
1731       *p++ = digit + '0';
1732     }
1733   last = p;
1734
1735   /* Generate one more digit with which to do rounding.  */
1736   do_multiply (&r, &r, ten);
1737   digit = rtd_divmod (&r, &pten);
1738
1739   /* Round the result.  */
1740   if (fmt && fmt->round_towards_zero)
1741     {
1742       /* If the format uses round towards zero when parsing the string
1743          back in, we need to always round away from zero here.  */
1744       if (cmp_significand_0 (&r))
1745         digit++;
1746       round_up = digit > 0;
1747     }
1748   else
1749     {
1750       if (digit == 5)
1751         {
1752           /* Round to nearest.  If R is nonzero there are additional
1753              nonzero digits to be extracted.  */
1754           if (cmp_significand_0 (&r))
1755             digit++;
1756           /* Round to even.  */
1757           else if ((p[-1] - '0') & 1)
1758             digit++;
1759         }
1760
1761       round_up = digit > 5;
1762     }
1763
1764   if (round_up)
1765     {
1766       while (p > first)
1767         {
1768           digit = *--p;
1769           if (digit == '9')
1770             *p = '0';
1771           else
1772             {
1773               *p = digit + 1;
1774               break;
1775             }
1776         }
1777
1778       /* Carry out of the first digit.  This means we had all 9's and
1779          now have all 0's.  "Prepend" a 1 by overwriting the first 0.  */
1780       if (p == first)
1781         {
1782           first[1] = '1';
1783           dec_exp++;
1784         }
1785     }
1786
1787   /* Insert the decimal point.  */
1788   first[0] = first[1];
1789   first[1] = '.';
1790
1791   /* If requested, drop trailing zeros.  Never crop past "1.0".  */
1792   if (crop_trailing_zeros)
1793     while (last > first + 3 && last[-1] == '0')
1794       last--;
1795
1796   /* Append the exponent.  */
1797   sprintf (last, "e%+d", dec_exp);
1798
1799 #ifdef ENABLE_CHECKING
1800   /* Verify that we can read the original value back in.  */
1801   if (mode != VOIDmode)
1802     {
1803       real_from_string (&r, str);
1804       real_convert (&r, mode, &r);
1805       gcc_assert (real_identical (&r, r_orig));
1806     }
1807 #endif
1808 }
1809
1810 /* Likewise, except always uses round-to-nearest.  */
1811
1812 void
1813 real_to_decimal (char *str, const REAL_VALUE_TYPE *r_orig, size_t buf_size,
1814                  size_t digits, int crop_trailing_zeros)
1815 {
1816   real_to_decimal_for_mode (str, r_orig, buf_size,
1817                             digits, crop_trailing_zeros, VOIDmode);
1818 }
1819
1820 /* Render R as a hexadecimal floating point constant.  Emit DIGITS
1821    significant digits in the result, bounded by BUF_SIZE.  If DIGITS is 0,
1822    choose the maximum for the representation.  If CROP_TRAILING_ZEROS,
1823    strip trailing zeros.  */
1824
1825 void
1826 real_to_hexadecimal (char *str, const REAL_VALUE_TYPE *r, size_t buf_size,
1827                      size_t digits, int crop_trailing_zeros)
1828 {
1829   int i, j, exp = REAL_EXP (r);
1830   char *p, *first;
1831   char exp_buf[16];
1832   size_t max_digits;
1833
1834   switch (r->cl)
1835     {
1836     case rvc_zero:
1837       exp = 0;
1838       break;
1839     case rvc_normal:
1840       break;
1841     case rvc_inf:
1842       strcpy (str, (r->sign ? "-Inf" : "+Inf"));
1843       return;
1844     case rvc_nan:
1845       /* ??? Print the significand as well, if not canonical?  */
1846       sprintf (str, "%c%cNaN", (r->sign ? '-' : '+'),
1847                (r->signalling ? 'S' : 'Q'));
1848       return;
1849     default:
1850       gcc_unreachable ();
1851     }
1852
1853   if (r->decimal)
1854     {
1855       /* Hexadecimal format for decimal floats is not interesting. */
1856       strcpy (str, "N/A");
1857       return;
1858     }
1859
1860   if (digits == 0)
1861     digits = SIGNIFICAND_BITS / 4;
1862
1863   /* Bound the number of digits printed by the size of the output buffer.  */
1864
1865   sprintf (exp_buf, "p%+d", exp);
1866   max_digits = buf_size - strlen (exp_buf) - r->sign - 4 - 1;
1867   gcc_assert (max_digits <= buf_size);
1868   if (digits > max_digits)
1869     digits = max_digits;
1870
1871   p = str;
1872   if (r->sign)
1873     *p++ = '-';
1874   *p++ = '0';
1875   *p++ = 'x';
1876   *p++ = '0';
1877   *p++ = '.';
1878   first = p;
1879
1880   for (i = SIGSZ - 1; i >= 0; --i)
1881     for (j = HOST_BITS_PER_LONG - 4; j >= 0; j -= 4)
1882       {
1883         *p++ = "0123456789abcdef"[(r->sig[i] >> j) & 15];
1884         if (--digits == 0)
1885           goto out;
1886       }
1887
1888  out:
1889   if (crop_trailing_zeros)
1890     while (p > first + 1 && p[-1] == '0')
1891       p--;
1892
1893   sprintf (p, "p%+d", exp);
1894 }
1895
1896 /* Initialize R from a decimal or hexadecimal string.  The string is
1897    assumed to have been syntax checked already.  Return -1 if the
1898    value underflows, +1 if overflows, and 0 otherwise. */
1899
1900 int
1901 real_from_string (REAL_VALUE_TYPE *r, const char *str)
1902 {
1903   int exp = 0;
1904   bool sign = false;
1905
1906   get_zero (r, 0);
1907
1908   if (*str == '-')
1909     {
1910       sign = true;
1911       str++;
1912     }
1913   else if (*str == '+')
1914     str++;
1915
1916   if (!strncmp (str, "QNaN", 4))
1917     {
1918       get_canonical_qnan (r, sign);
1919       return 0;
1920     }
1921   else if (!strncmp (str, "SNaN", 4))
1922     {
1923       get_canonical_snan (r, sign);
1924       return 0;
1925     }
1926   else if (!strncmp (str, "Inf", 3))
1927     {
1928       get_inf (r, sign);
1929       return 0;
1930     }
1931
1932   if (str[0] == '0' && (str[1] == 'x' || str[1] == 'X'))
1933     {
1934       /* Hexadecimal floating point.  */
1935       int pos = SIGNIFICAND_BITS - 4, d;
1936
1937       str += 2;
1938
1939       while (*str == '0')
1940         str++;
1941       while (1)
1942         {
1943           d = hex_value (*str);
1944           if (d == _hex_bad)
1945             break;
1946           if (pos >= 0)
1947             {
1948               r->sig[pos / HOST_BITS_PER_LONG]
1949                 |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1950               pos -= 4;
1951             }
1952           else if (d)
1953             /* Ensure correct rounding by setting last bit if there is
1954                a subsequent nonzero digit.  */
1955             r->sig[0] |= 1;
1956           exp += 4;
1957           str++;
1958         }
1959       if (*str == '.')
1960         {
1961           str++;
1962           if (pos == SIGNIFICAND_BITS - 4)
1963             {
1964               while (*str == '0')
1965                 str++, exp -= 4;
1966             }
1967           while (1)
1968             {
1969               d = hex_value (*str);
1970               if (d == _hex_bad)
1971                 break;
1972               if (pos >= 0)
1973                 {
1974                   r->sig[pos / HOST_BITS_PER_LONG]
1975                     |= (unsigned long) d << (pos % HOST_BITS_PER_LONG);
1976                   pos -= 4;
1977                 }
1978               else if (d)
1979                 /* Ensure correct rounding by setting last bit if there is
1980                    a subsequent nonzero digit.  */
1981                 r->sig[0] |= 1;
1982               str++;
1983             }
1984         }
1985
1986       /* If the mantissa is zero, ignore the exponent.  */
1987       if (!cmp_significand_0 (r))
1988         goto is_a_zero;
1989
1990       if (*str == 'p' || *str == 'P')
1991         {
1992           bool exp_neg = false;
1993
1994           str++;
1995           if (*str == '-')
1996             {
1997               exp_neg = true;
1998               str++;
1999             }
2000           else if (*str == '+')
2001             str++;
2002
2003           d = 0;
2004           while (ISDIGIT (*str))
2005             {
2006               d *= 10;
2007               d += *str - '0';
2008               if (d > MAX_EXP)
2009                 {
2010                   /* Overflowed the exponent.  */
2011                   if (exp_neg)
2012                     goto underflow;
2013                   else
2014                     goto overflow;
2015                 }
2016               str++;
2017             }
2018           if (exp_neg)
2019             d = -d;
2020
2021           exp += d;
2022         }
2023
2024       r->cl = rvc_normal;
2025       SET_REAL_EXP (r, exp);
2026
2027       normalize (r);
2028     }
2029   else
2030     {
2031       /* Decimal floating point.  */
2032       const REAL_VALUE_TYPE *ten = ten_to_ptwo (0);
2033       int d;
2034
2035       while (*str == '0')
2036         str++;
2037       while (ISDIGIT (*str))
2038         {
2039           d = *str++ - '0';
2040           do_multiply (r, r, ten);
2041           if (d)
2042             do_add (r, r, real_digit (d), 0);
2043         }
2044       if (*str == '.')
2045         {
2046           str++;
2047           if (r->cl == rvc_zero)
2048             {
2049               while (*str == '0')
2050                 str++, exp--;
2051             }
2052           while (ISDIGIT (*str))
2053             {
2054               d = *str++ - '0';
2055               do_multiply (r, r, ten);
2056               if (d)
2057                 do_add (r, r, real_digit (d), 0);
2058               exp--;
2059             }
2060         }
2061
2062       /* If the mantissa is zero, ignore the exponent.  */
2063       if (r->cl == rvc_zero)
2064         goto is_a_zero;
2065
2066       if (*str == 'e' || *str == 'E')
2067         {
2068           bool exp_neg = false;
2069
2070           str++;
2071           if (*str == '-')
2072             {
2073               exp_neg = true;
2074               str++;
2075             }
2076           else if (*str == '+')
2077             str++;
2078
2079           d = 0;
2080           while (ISDIGIT (*str))
2081             {
2082               d *= 10;
2083               d += *str - '0';
2084               if (d > MAX_EXP)
2085                 {
2086                   /* Overflowed the exponent.  */
2087                   if (exp_neg)
2088                     goto underflow;
2089                   else
2090                     goto overflow;
2091                 }
2092               str++;
2093             }
2094           if (exp_neg)
2095             d = -d;
2096           exp += d;
2097         }
2098
2099       if (exp)
2100         times_pten (r, exp);
2101     }
2102
2103   r->sign = sign;
2104   return 0;
2105
2106  is_a_zero:
2107   get_zero (r, sign);
2108   return 0;
2109
2110  underflow:
2111   get_zero (r, sign);
2112   return -1;
2113
2114  overflow:
2115   get_inf (r, sign);
2116   return 1;
2117 }
2118
2119 /* Legacy.  Similar, but return the result directly.  */
2120
2121 REAL_VALUE_TYPE
2122 real_from_string2 (const char *s, enum machine_mode mode)
2123 {
2124   REAL_VALUE_TYPE r;
2125
2126   real_from_string (&r, s);
2127   if (mode != VOIDmode)
2128     real_convert (&r, mode, &r);
2129
2130   return r;
2131 }
2132
2133 /* Initialize R from string S and desired MODE. */
2134
2135 void
2136 real_from_string3 (REAL_VALUE_TYPE *r, const char *s, enum machine_mode mode)
2137 {
2138   if (DECIMAL_FLOAT_MODE_P (mode))
2139     decimal_real_from_string (r, s);
2140   else
2141     real_from_string (r, s);
2142
2143   if (mode != VOIDmode)
2144     real_convert (r, mode, r);
2145 }
2146
2147 /* Initialize R from the integer pair HIGH+LOW.  */
2148
2149 void
2150 real_from_integer (REAL_VALUE_TYPE *r, enum machine_mode mode,
2151                    unsigned HOST_WIDE_INT low, HOST_WIDE_INT high,
2152                    int unsigned_p)
2153 {
2154   if (low == 0 && high == 0)
2155     get_zero (r, 0);
2156   else
2157     {
2158       memset (r, 0, sizeof (*r));
2159       r->cl = rvc_normal;
2160       r->sign = high < 0 && !unsigned_p;
2161       SET_REAL_EXP (r, HOST_BITS_PER_DOUBLE_INT);
2162
2163       if (r->sign)
2164         {
2165           high = ~high;
2166           if (low == 0)
2167             high += 1;
2168           else
2169             low = -low;
2170         }
2171
2172       if (HOST_BITS_PER_LONG == HOST_BITS_PER_WIDE_INT)
2173         {
2174           r->sig[SIGSZ-1] = high;
2175           r->sig[SIGSZ-2] = low;
2176         }
2177       else
2178         {
2179           gcc_assert (HOST_BITS_PER_LONG*2 == HOST_BITS_PER_WIDE_INT);
2180           r->sig[SIGSZ-1] = high >> (HOST_BITS_PER_LONG - 1) >> 1;
2181           r->sig[SIGSZ-2] = high;
2182           r->sig[SIGSZ-3] = low >> (HOST_BITS_PER_LONG - 1) >> 1;
2183           r->sig[SIGSZ-4] = low;
2184         }
2185
2186       normalize (r);
2187     }
2188
2189   if (DECIMAL_FLOAT_MODE_P (mode))
2190     decimal_from_integer (r);
2191   else if (mode != VOIDmode)
2192     real_convert (r, mode, r);
2193 }
2194
2195 /* Render R, an integral value, as a floating point constant with no
2196    specified exponent.  */
2197
2198 static void
2199 decimal_integer_string (char *str, const REAL_VALUE_TYPE *r_orig,
2200                         size_t buf_size)
2201 {
2202   int dec_exp, digit, digits;
2203   REAL_VALUE_TYPE r, pten;
2204   char *p;
2205   bool sign;
2206
2207   r = *r_orig;
2208
2209   if (r.cl == rvc_zero)
2210     {
2211       strcpy (str, "0.");
2212       return;
2213     }
2214
2215   sign = r.sign;
2216   r.sign = 0;
2217
2218   dec_exp = REAL_EXP (&r) * M_LOG10_2;
2219   digits = dec_exp + 1;
2220   gcc_assert ((digits + 2) < (int)buf_size);
2221
2222   pten = *real_digit (1);
2223   times_pten (&pten, dec_exp);
2224
2225   p = str;
2226   if (sign)
2227     *p++ = '-';
2228
2229   digit = rtd_divmod (&r, &pten);
2230   gcc_assert (digit >= 0 && digit <= 9);
2231   *p++ = digit + '0';
2232   while (--digits > 0)
2233     {
2234       times_pten (&r, 1);
2235       digit = rtd_divmod (&r, &pten);
2236       *p++ = digit + '0';
2237     }
2238   *p++ = '.';
2239   *p++ = '\0';
2240 }
2241
2242 /* Convert a real with an integral value to decimal float.  */
2243
2244 static void
2245 decimal_from_integer (REAL_VALUE_TYPE *r)
2246 {
2247   char str[256];
2248
2249   decimal_integer_string (str, r, sizeof (str) - 1);
2250   decimal_real_from_string (r, str);
2251 }
2252
2253 /* Returns 10**2**N.  */
2254
2255 static const REAL_VALUE_TYPE *
2256 ten_to_ptwo (int n)
2257 {
2258   static REAL_VALUE_TYPE tens[EXP_BITS];
2259
2260   gcc_assert (n >= 0);
2261   gcc_assert (n < EXP_BITS);
2262
2263   if (tens[n].cl == rvc_zero)
2264     {
2265       if (n < (HOST_BITS_PER_WIDE_INT == 64 ? 5 : 4))
2266         {
2267           HOST_WIDE_INT t = 10;
2268           int i;
2269
2270           for (i = 0; i < n; ++i)
2271             t *= t;
2272
2273           real_from_integer (&tens[n], VOIDmode, t, 0, 1);
2274         }
2275       else
2276         {
2277           const REAL_VALUE_TYPE *t = ten_to_ptwo (n - 1);
2278           do_multiply (&tens[n], t, t);
2279         }
2280     }
2281
2282   return &tens[n];
2283 }
2284
2285 /* Returns 10**(-2**N).  */
2286
2287 static const REAL_VALUE_TYPE *
2288 ten_to_mptwo (int n)
2289 {
2290   static REAL_VALUE_TYPE tens[EXP_BITS];
2291
2292   gcc_assert (n >= 0);
2293   gcc_assert (n < EXP_BITS);
2294
2295   if (tens[n].cl == rvc_zero)
2296     do_divide (&tens[n], real_digit (1), ten_to_ptwo (n));
2297
2298   return &tens[n];
2299 }
2300
2301 /* Returns N.  */
2302
2303 static const REAL_VALUE_TYPE *
2304 real_digit (int n)
2305 {
2306   static REAL_VALUE_TYPE num[10];
2307
2308   gcc_assert (n >= 0);
2309   gcc_assert (n <= 9);
2310
2311   if (n > 0 && num[n].cl == rvc_zero)
2312     real_from_integer (&num[n], VOIDmode, n, 0, 1);
2313
2314   return &num[n];
2315 }
2316
2317 /* Multiply R by 10**EXP.  */
2318
2319 static void
2320 times_pten (REAL_VALUE_TYPE *r, int exp)
2321 {
2322   REAL_VALUE_TYPE pten, *rr;
2323   bool negative = (exp < 0);
2324   int i;
2325
2326   if (negative)
2327     {
2328       exp = -exp;
2329       pten = *real_digit (1);
2330       rr = &pten;
2331     }
2332   else
2333     rr = r;
2334
2335   for (i = 0; exp > 0; ++i, exp >>= 1)
2336     if (exp & 1)
2337       do_multiply (rr, rr, ten_to_ptwo (i));
2338
2339   if (negative)
2340     do_divide (r, r, &pten);
2341 }
2342
2343 /* Returns the special REAL_VALUE_TYPE corresponding to 'e'.  */
2344
2345 const REAL_VALUE_TYPE *
2346 dconst_e_ptr (void)
2347 {
2348   static REAL_VALUE_TYPE value;
2349
2350   /* Initialize mathematical constants for constant folding builtins.
2351      These constants need to be given to at least 160 bits precision.  */
2352   if (value.cl == rvc_zero)
2353     {
2354       mpfr_t m;
2355       mpfr_init2 (m, SIGNIFICAND_BITS);
2356       mpfr_set_ui (m, 1, GMP_RNDN);
2357       mpfr_exp (m, m, GMP_RNDN);
2358       real_from_mpfr (&value, m, NULL_TREE, GMP_RNDN);
2359       mpfr_clear (m);
2360
2361     }
2362   return &value;
2363 }
2364
2365 /* Returns the special REAL_VALUE_TYPE corresponding to 1/3.  */
2366
2367 const REAL_VALUE_TYPE *
2368 dconst_third_ptr (void)
2369 {
2370   static REAL_VALUE_TYPE value;
2371
2372   /* Initialize mathematical constants for constant folding builtins.
2373      These constants need to be given to at least 160 bits precision.  */
2374   if (value.cl == rvc_zero)
2375     {
2376       real_arithmetic (&value, RDIV_EXPR, &dconst1, real_digit (3));
2377     }
2378   return &value;
2379 }
2380
2381 /* Returns the special REAL_VALUE_TYPE corresponding to sqrt(2).  */
2382
2383 const REAL_VALUE_TYPE *
2384 dconst_sqrt2_ptr (void)
2385 {
2386   static REAL_VALUE_TYPE value;
2387
2388   /* Initialize mathematical constants for constant folding builtins.
2389      These constants need to be given to at least 160 bits precision.  */
2390   if (value.cl == rvc_zero)
2391     {
2392       mpfr_t m;
2393       mpfr_init2 (m, SIGNIFICAND_BITS);
2394       mpfr_sqrt_ui (m, 2, GMP_RNDN);
2395       real_from_mpfr (&value, m, NULL_TREE, GMP_RNDN);
2396       mpfr_clear (m);
2397     }
2398   return &value;
2399 }
2400
2401 /* Fills R with +Inf.  */
2402
2403 void
2404 real_inf (REAL_VALUE_TYPE *r)
2405 {
2406   get_inf (r, 0);
2407 }
2408
2409 /* Fills R with a NaN whose significand is described by STR.  If QUIET,
2410    we force a QNaN, else we force an SNaN.  The string, if not empty,
2411    is parsed as a number and placed in the significand.  Return true
2412    if the string was successfully parsed.  */
2413
2414 bool
2415 real_nan (REAL_VALUE_TYPE *r, const char *str, int quiet,
2416           enum machine_mode mode)
2417 {
2418   const struct real_format *fmt;
2419
2420   fmt = REAL_MODE_FORMAT (mode);
2421   gcc_assert (fmt);
2422
2423   if (*str == 0)
2424     {
2425       if (quiet)
2426         get_canonical_qnan (r, 0);
2427       else
2428         get_canonical_snan (r, 0);
2429     }
2430   else
2431     {
2432       int base = 10, d;
2433
2434       memset (r, 0, sizeof (*r));
2435       r->cl = rvc_nan;
2436
2437       /* Parse akin to strtol into the significand of R.  */
2438
2439       while (ISSPACE (*str))
2440         str++;
2441       if (*str == '-')
2442         str++;
2443       else if (*str == '+')
2444         str++;
2445       if (*str == '0')
2446         {
2447           str++;
2448           if (*str == 'x' || *str == 'X')
2449             {
2450               base = 16;
2451               str++;
2452             }
2453           else
2454             base = 8;
2455         }
2456
2457       while ((d = hex_value (*str)) < base)
2458         {
2459           REAL_VALUE_TYPE u;
2460
2461           switch (base)
2462             {
2463             case 8:
2464               lshift_significand (r, r, 3);
2465               break;
2466             case 16:
2467               lshift_significand (r, r, 4);
2468               break;
2469             case 10:
2470               lshift_significand_1 (&u, r);
2471               lshift_significand (r, r, 3);
2472               add_significands (r, r, &u);
2473               break;
2474             default:
2475               gcc_unreachable ();
2476             }
2477
2478           get_zero (&u, 0);
2479           u.sig[0] = d;
2480           add_significands (r, r, &u);
2481
2482           str++;
2483         }
2484
2485       /* Must have consumed the entire string for success.  */
2486       if (*str != 0)
2487         return false;
2488
2489       /* Shift the significand into place such that the bits
2490          are in the most significant bits for the format.  */
2491       lshift_significand (r, r, SIGNIFICAND_BITS - fmt->pnan);
2492
2493       /* Our MSB is always unset for NaNs.  */
2494       r->sig[SIGSZ-1] &= ~SIG_MSB;
2495
2496       /* Force quiet or signalling NaN.  */
2497       r->signalling = !quiet;
2498     }
2499
2500   return true;
2501 }
2502
2503 /* Fills R with the largest finite value representable in mode MODE.
2504    If SIGN is nonzero, R is set to the most negative finite value.  */
2505
2506 void
2507 real_maxval (REAL_VALUE_TYPE *r, int sign, enum machine_mode mode)
2508 {
2509   const struct real_format *fmt;
2510   int np2;
2511
2512   fmt = REAL_MODE_FORMAT (mode);
2513   gcc_assert (fmt);
2514   memset (r, 0, sizeof (*r));
2515
2516   if (fmt->b == 10)
2517     decimal_real_maxval (r, sign, mode);
2518   else
2519     {
2520       r->cl = rvc_normal;
2521       r->sign = sign;
2522       SET_REAL_EXP (r, fmt->emax);
2523
2524       np2 = SIGNIFICAND_BITS - fmt->p;
2525       memset (r->sig, -1, SIGSZ * sizeof (unsigned long));
2526       clear_significand_below (r, np2);
2527
2528       if (fmt->pnan < fmt->p)
2529         /* This is an IBM extended double format made up of two IEEE
2530            doubles.  The value of the long double is the sum of the
2531            values of the two parts.  The most significant part is
2532            required to be the value of the long double rounded to the
2533            nearest double.  Rounding means we need a slightly smaller
2534            value for LDBL_MAX.  */
2535         clear_significand_bit (r, SIGNIFICAND_BITS - fmt->pnan - 1);
2536     }
2537 }
2538
2539 /* Fills R with 2**N.  */
2540
2541 void
2542 real_2expN (REAL_VALUE_TYPE *r, int n, enum machine_mode fmode)
2543 {
2544   memset (r, 0, sizeof (*r));
2545
2546   n++;
2547   if (n > MAX_EXP)
2548     r->cl = rvc_inf;
2549   else if (n < -MAX_EXP)
2550     ;
2551   else
2552     {
2553       r->cl = rvc_normal;
2554       SET_REAL_EXP (r, n);
2555       r->sig[SIGSZ-1] = SIG_MSB;
2556     }
2557   if (DECIMAL_FLOAT_MODE_P (fmode))
2558     decimal_real_convert (r, fmode, r);
2559 }
2560
2561 \f
2562 static void
2563 round_for_format (const struct real_format *fmt, REAL_VALUE_TYPE *r)
2564 {
2565   int p2, np2, i, w;
2566   int emin2m1, emax2;
2567   bool round_up = false;
2568
2569   if (r->decimal)
2570     {
2571       if (fmt->b == 10)
2572         {
2573           decimal_round_for_format (fmt, r);
2574           return;
2575         }
2576       /* FIXME. We can come here via fp_easy_constant
2577          (e.g. -O0 on '_Decimal32 x = 1.0 + 2.0dd'), but have not
2578          investigated whether this convert needs to be here, or
2579          something else is missing. */
2580       decimal_real_convert (r, DFmode, r);
2581     }
2582
2583   p2 = fmt->p;
2584   emin2m1 = fmt->emin - 1;
2585   emax2 = fmt->emax;
2586
2587   np2 = SIGNIFICAND_BITS - p2;
2588   switch (r->cl)
2589     {
2590     underflow:
2591       get_zero (r, r->sign);
2592     case rvc_zero:
2593       if (!fmt->has_signed_zero)
2594         r->sign = 0;
2595       return;
2596
2597     overflow:
2598       get_inf (r, r->sign);
2599     case rvc_inf:
2600       return;
2601
2602     case rvc_nan:
2603       clear_significand_below (r, np2);
2604       return;
2605
2606     case rvc_normal:
2607       break;
2608
2609     default:
2610       gcc_unreachable ();
2611     }
2612
2613   /* Check the range of the exponent.  If we're out of range,
2614      either underflow or overflow.  */
2615   if (REAL_EXP (r) > emax2)
2616     goto overflow;
2617   else if (REAL_EXP (r) <= emin2m1)
2618     {
2619       int diff;
2620
2621       if (!fmt->has_denorm)
2622         {
2623           /* Don't underflow completely until we've had a chance to round.  */
2624           if (REAL_EXP (r) < emin2m1)
2625             goto underflow;
2626         }
2627       else
2628         {
2629           diff = emin2m1 - REAL_EXP (r) + 1;
2630           if (diff > p2)
2631             goto underflow;
2632
2633           /* De-normalize the significand.  */
2634           r->sig[0] |= sticky_rshift_significand (r, r, diff);
2635           SET_REAL_EXP (r, REAL_EXP (r) + diff);
2636         }
2637     }
2638
2639   if (!fmt->round_towards_zero)
2640     {
2641       /* There are P2 true significand bits, followed by one guard bit,
2642          followed by one sticky bit, followed by stuff.  Fold nonzero
2643          stuff into the sticky bit.  */
2644       unsigned long sticky;
2645       bool guard, lsb;
2646
2647       sticky = 0;
2648       for (i = 0, w = (np2 - 1) / HOST_BITS_PER_LONG; i < w; ++i)
2649         sticky |= r->sig[i];
2650       sticky |= r->sig[w]
2651                 & (((unsigned long)1 << ((np2 - 1) % HOST_BITS_PER_LONG)) - 1);
2652
2653       guard = test_significand_bit (r, np2 - 1);
2654       lsb = test_significand_bit (r, np2);
2655
2656       /* Round to even.  */
2657       round_up = guard && (sticky || lsb);
2658     }
2659
2660   if (round_up)
2661     {
2662       REAL_VALUE_TYPE u;
2663       get_zero (&u, 0);
2664       set_significand_bit (&u, np2);
2665
2666       if (add_significands (r, r, &u))
2667         {
2668           /* Overflow.  Means the significand had been all ones, and
2669              is now all zeros.  Need to increase the exponent, and
2670              possibly re-normalize it.  */
2671           SET_REAL_EXP (r, REAL_EXP (r) + 1);
2672           if (REAL_EXP (r) > emax2)
2673             goto overflow;
2674           r->sig[SIGSZ-1] = SIG_MSB;
2675         }
2676     }
2677
2678   /* Catch underflow that we deferred until after rounding.  */
2679   if (REAL_EXP (r) <= emin2m1)
2680     goto underflow;
2681
2682   /* Clear out trailing garbage.  */
2683   clear_significand_below (r, np2);
2684 }
2685
2686 /* Extend or truncate to a new mode.  */
2687
2688 void
2689 real_convert (REAL_VALUE_TYPE *r, enum machine_mode mode,
2690               const REAL_VALUE_TYPE *a)
2691 {
2692   const struct real_format *fmt;
2693
2694   fmt = REAL_MODE_FORMAT (mode);
2695   gcc_assert (fmt);
2696
2697   *r = *a;
2698
2699   if (a->decimal || fmt->b == 10)
2700     decimal_real_convert (r, mode, a);
2701
2702   round_for_format (fmt, r);
2703
2704   /* round_for_format de-normalizes denormals.  Undo just that part.  */
2705   if (r->cl == rvc_normal)
2706     normalize (r);
2707 }
2708
2709 /* Legacy.  Likewise, except return the struct directly.  */
2710
2711 REAL_VALUE_TYPE
2712 real_value_truncate (enum machine_mode mode, REAL_VALUE_TYPE a)
2713 {
2714   REAL_VALUE_TYPE r;
2715   real_convert (&r, mode, &a);
2716   return r;
2717 }
2718
2719 /* Return true if truncating to MODE is exact.  */
2720
2721 bool
2722 exact_real_truncate (enum machine_mode mode, const REAL_VALUE_TYPE *a)
2723 {
2724   const struct real_format *fmt;
2725   REAL_VALUE_TYPE t;
2726   int emin2m1;
2727
2728   fmt = REAL_MODE_FORMAT (mode);
2729   gcc_assert (fmt);
2730
2731   /* Don't allow conversion to denormals.  */
2732   emin2m1 = fmt->emin - 1;
2733   if (REAL_EXP (a) <= emin2m1)
2734     return false;
2735
2736   /* After conversion to the new mode, the value must be identical.  */
2737   real_convert (&t, mode, a);
2738   return real_identical (&t, a);
2739 }
2740
2741 /* Write R to the given target format.  Place the words of the result
2742    in target word order in BUF.  There are always 32 bits in each
2743    long, no matter the size of the host long.
2744
2745    Legacy: return word 0 for implementing REAL_VALUE_TO_TARGET_SINGLE.  */
2746
2747 long
2748 real_to_target_fmt (long *buf, const REAL_VALUE_TYPE *r_orig,
2749                     const struct real_format *fmt)
2750 {
2751   REAL_VALUE_TYPE r;
2752   long buf1;
2753
2754   r = *r_orig;
2755   round_for_format (fmt, &r);
2756
2757   if (!buf)
2758     buf = &buf1;
2759   (*fmt->encode) (fmt, buf, &r);
2760
2761   return *buf;
2762 }
2763
2764 /* Similar, but look up the format from MODE.  */
2765
2766 long
2767 real_to_target (long *buf, const REAL_VALUE_TYPE *r, enum machine_mode mode)
2768 {
2769   const struct real_format *fmt;
2770
2771   fmt = REAL_MODE_FORMAT (mode);
2772   gcc_assert (fmt);
2773
2774   return real_to_target_fmt (buf, r, fmt);
2775 }
2776
2777 /* Read R from the given target format.  Read the words of the result
2778    in target word order in BUF.  There are always 32 bits in each
2779    long, no matter the size of the host long.  */
2780
2781 void
2782 real_from_target_fmt (REAL_VALUE_TYPE *r, const long *buf,
2783                       const struct real_format *fmt)
2784 {
2785   (*fmt->decode) (fmt, r, buf);
2786 }
2787
2788 /* Similar, but look up the format from MODE.  */
2789
2790 void
2791 real_from_target (REAL_VALUE_TYPE *r, const long *buf, enum machine_mode mode)
2792 {
2793   const struct real_format *fmt;
2794
2795   fmt = REAL_MODE_FORMAT (mode);
2796   gcc_assert (fmt);
2797
2798   (*fmt->decode) (fmt, r, buf);
2799 }
2800
2801 /* Return the number of bits of the largest binary value that the
2802    significand of MODE will hold.  */
2803 /* ??? Legacy.  Should get access to real_format directly.  */
2804
2805 int
2806 significand_size (enum machine_mode mode)
2807 {
2808   const struct real_format *fmt;
2809
2810   fmt = REAL_MODE_FORMAT (mode);
2811   if (fmt == NULL)
2812     return 0;
2813
2814   if (fmt->b == 10)
2815     {
2816       /* Return the size in bits of the largest binary value that can be
2817          held by the decimal coefficient for this mode.  This is one more
2818          than the number of bits required to hold the largest coefficient
2819          of this mode.  */
2820       double log2_10 = 3.3219281;
2821       return fmt->p * log2_10;
2822     }
2823   return fmt->p;
2824 }
2825
2826 /* Return a hash value for the given real value.  */
2827 /* ??? The "unsigned int" return value is intended to be hashval_t,
2828    but I didn't want to pull hashtab.h into real.h.  */
2829
2830 unsigned int
2831 real_hash (const REAL_VALUE_TYPE *r)
2832 {
2833   unsigned int h;
2834   size_t i;
2835
2836   h = r->cl | (r->sign << 2);
2837   switch (r->cl)
2838     {
2839     case rvc_zero:
2840     case rvc_inf:
2841       return h;
2842
2843     case rvc_normal:
2844       h |= REAL_EXP (r) << 3;
2845       break;
2846
2847     case rvc_nan:
2848       if (r->signalling)
2849         h ^= (unsigned int)-1;
2850       if (r->canonical)
2851         return h;
2852       break;
2853
2854     default:
2855       gcc_unreachable ();
2856     }
2857
2858   if (sizeof(unsigned long) > sizeof(unsigned int))
2859     for (i = 0; i < SIGSZ; ++i)
2860       {
2861         unsigned long s = r->sig[i];
2862         h ^= s ^ (s >> (HOST_BITS_PER_LONG / 2));
2863       }
2864   else
2865     for (i = 0; i < SIGSZ; ++i)
2866       h ^= r->sig[i];
2867
2868   return h;
2869 }
2870 \f
2871 /* IEEE single-precision format.  */
2872
2873 static void encode_ieee_single (const struct real_format *fmt,
2874                                 long *, const REAL_VALUE_TYPE *);
2875 static void decode_ieee_single (const struct real_format *,
2876                                 REAL_VALUE_TYPE *, const long *);
2877
2878 static void
2879 encode_ieee_single (const struct real_format *fmt, long *buf,
2880                     const REAL_VALUE_TYPE *r)
2881 {
2882   unsigned long image, sig, exp;
2883   unsigned long sign = r->sign;
2884   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
2885
2886   image = sign << 31;
2887   sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
2888
2889   switch (r->cl)
2890     {
2891     case rvc_zero:
2892       break;
2893
2894     case rvc_inf:
2895       if (fmt->has_inf)
2896         image |= 255 << 23;
2897       else
2898         image |= 0x7fffffff;
2899       break;
2900
2901     case rvc_nan:
2902       if (fmt->has_nans)
2903         {
2904           if (r->canonical)
2905             sig = (fmt->canonical_nan_lsbs_set ? (1 << 22) - 1 : 0);
2906           if (r->signalling == fmt->qnan_msb_set)
2907             sig &= ~(1 << 22);
2908           else
2909             sig |= 1 << 22;
2910           if (sig == 0)
2911             sig = 1 << 21;
2912
2913           image |= 255 << 23;
2914           image |= sig;
2915         }
2916       else
2917         image |= 0x7fffffff;
2918       break;
2919
2920     case rvc_normal:
2921       /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
2922          whereas the intermediate representation is 0.F x 2**exp.
2923          Which means we're off by one.  */
2924       if (denormal)
2925         exp = 0;
2926       else
2927       exp = REAL_EXP (r) + 127 - 1;
2928       image |= exp << 23;
2929       image |= sig;
2930       break;
2931
2932     default:
2933       gcc_unreachable ();
2934     }
2935
2936   buf[0] = image;
2937 }
2938
2939 static void
2940 decode_ieee_single (const struct real_format *fmt, REAL_VALUE_TYPE *r,
2941                     const long *buf)
2942 {
2943   unsigned long image = buf[0] & 0xffffffff;
2944   bool sign = (image >> 31) & 1;
2945   int exp = (image >> 23) & 0xff;
2946
2947   memset (r, 0, sizeof (*r));
2948   image <<= HOST_BITS_PER_LONG - 24;
2949   image &= ~SIG_MSB;
2950
2951   if (exp == 0)
2952     {
2953       if (image && fmt->has_denorm)
2954         {
2955           r->cl = rvc_normal;
2956           r->sign = sign;
2957           SET_REAL_EXP (r, -126);
2958           r->sig[SIGSZ-1] = image << 1;
2959           normalize (r);
2960         }
2961       else if (fmt->has_signed_zero)
2962         r->sign = sign;
2963     }
2964   else if (exp == 255 && (fmt->has_nans || fmt->has_inf))
2965     {
2966       if (image)
2967         {
2968           r->cl = rvc_nan;
2969           r->sign = sign;
2970           r->signalling = (((image >> (HOST_BITS_PER_LONG - 2)) & 1)
2971                            ^ fmt->qnan_msb_set);
2972           r->sig[SIGSZ-1] = image;
2973         }
2974       else
2975         {
2976           r->cl = rvc_inf;
2977           r->sign = sign;
2978         }
2979     }
2980   else
2981     {
2982       r->cl = rvc_normal;
2983       r->sign = sign;
2984       SET_REAL_EXP (r, exp - 127 + 1);
2985       r->sig[SIGSZ-1] = image | SIG_MSB;
2986     }
2987 }
2988
2989 const struct real_format ieee_single_format =
2990   {
2991     encode_ieee_single,
2992     decode_ieee_single,
2993     2,
2994     24,
2995     24,
2996     -125,
2997     128,
2998     31,
2999     31,
3000     false,
3001     true,
3002     true,
3003     true,
3004     true,
3005     true,
3006     true,
3007     false
3008   };
3009
3010 const struct real_format mips_single_format =
3011   {
3012     encode_ieee_single,
3013     decode_ieee_single,
3014     2,
3015     24,
3016     24,
3017     -125,
3018     128,
3019     31,
3020     31,
3021     false,
3022     true,
3023     true,
3024     true,
3025     true,
3026     true,
3027     false,
3028     true
3029   };
3030
3031 const struct real_format motorola_single_format =
3032   {
3033     encode_ieee_single,
3034     decode_ieee_single,
3035     2,
3036     24,
3037     24,
3038     -125,
3039     128,
3040     31,
3041     31,
3042     false,
3043     true,
3044     true,
3045     true,
3046     true,
3047     true,
3048     true,
3049     true
3050   };
3051
3052 /*  SPU Single Precision (Extended-Range Mode) format is the same as IEEE
3053     single precision with the following differences:
3054       - Infinities are not supported.  Instead MAX_FLOAT or MIN_FLOAT
3055         are generated.
3056       - NaNs are not supported.
3057       - The range of non-zero numbers in binary is
3058         (001)[1.]000...000 to (255)[1.]111...111.
3059       - Denormals can be represented, but are treated as +0.0 when
3060         used as an operand and are never generated as a result.
3061       - -0.0 can be represented, but a zero result is always +0.0.
3062       - the only supported rounding mode is trunction (towards zero).  */
3063 const struct real_format spu_single_format =
3064   {
3065     encode_ieee_single,
3066     decode_ieee_single,
3067     2,
3068     24,
3069     24,
3070     -125,
3071     129,
3072     31,
3073     31,
3074     true,
3075     false,
3076     false,
3077     false,
3078     true,
3079     true,
3080     false,
3081     false
3082   };
3083 \f
3084 /* IEEE double-precision format.  */
3085
3086 static void encode_ieee_double (const struct real_format *fmt,
3087                                 long *, const REAL_VALUE_TYPE *);
3088 static void decode_ieee_double (const struct real_format *,
3089                                 REAL_VALUE_TYPE *, const long *);
3090
3091 static void
3092 encode_ieee_double (const struct real_format *fmt, long *buf,
3093                     const REAL_VALUE_TYPE *r)
3094 {
3095   unsigned long image_lo, image_hi, sig_lo, sig_hi, exp;
3096   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3097
3098   image_hi = r->sign << 31;
3099   image_lo = 0;
3100
3101   if (HOST_BITS_PER_LONG == 64)
3102     {
3103       sig_hi = r->sig[SIGSZ-1];
3104       sig_lo = (sig_hi >> (64 - 53)) & 0xffffffff;
3105       sig_hi = (sig_hi >> (64 - 53 + 1) >> 31) & 0xfffff;
3106     }
3107   else
3108     {
3109       sig_hi = r->sig[SIGSZ-1];
3110       sig_lo = r->sig[SIGSZ-2];
3111       sig_lo = (sig_hi << 21) | (sig_lo >> 11);
3112       sig_hi = (sig_hi >> 11) & 0xfffff;
3113     }
3114
3115   switch (r->cl)
3116     {
3117     case rvc_zero:
3118       break;
3119
3120     case rvc_inf:
3121       if (fmt->has_inf)
3122         image_hi |= 2047 << 20;
3123       else
3124         {
3125           image_hi |= 0x7fffffff;
3126           image_lo = 0xffffffff;
3127         }
3128       break;
3129
3130     case rvc_nan:
3131       if (fmt->has_nans)
3132         {
3133           if (r->canonical)
3134             {
3135               if (fmt->canonical_nan_lsbs_set)
3136                 {
3137                   sig_hi = (1 << 19) - 1;
3138                   sig_lo = 0xffffffff;
3139                 }
3140               else
3141                 {
3142                   sig_hi = 0;
3143                   sig_lo = 0;
3144                 }
3145             }
3146           if (r->signalling == fmt->qnan_msb_set)
3147             sig_hi &= ~(1 << 19);
3148           else
3149             sig_hi |= 1 << 19;
3150           if (sig_hi == 0 && sig_lo == 0)
3151             sig_hi = 1 << 18;
3152
3153           image_hi |= 2047 << 20;
3154           image_hi |= sig_hi;
3155           image_lo = sig_lo;
3156         }
3157       else
3158         {
3159           image_hi |= 0x7fffffff;
3160           image_lo = 0xffffffff;
3161         }
3162       break;
3163
3164     case rvc_normal:
3165       /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3166          whereas the intermediate representation is 0.F x 2**exp.
3167          Which means we're off by one.  */
3168       if (denormal)
3169         exp = 0;
3170       else
3171         exp = REAL_EXP (r) + 1023 - 1;
3172       image_hi |= exp << 20;
3173       image_hi |= sig_hi;
3174       image_lo = sig_lo;
3175       break;
3176
3177     default:
3178       gcc_unreachable ();
3179     }
3180
3181   if (FLOAT_WORDS_BIG_ENDIAN)
3182     buf[0] = image_hi, buf[1] = image_lo;
3183   else
3184     buf[0] = image_lo, buf[1] = image_hi;
3185 }
3186
3187 static void
3188 decode_ieee_double (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3189                     const long *buf)
3190 {
3191   unsigned long image_hi, image_lo;
3192   bool sign;
3193   int exp;
3194
3195   if (FLOAT_WORDS_BIG_ENDIAN)
3196     image_hi = buf[0], image_lo = buf[1];
3197   else
3198     image_lo = buf[0], image_hi = buf[1];
3199   image_lo &= 0xffffffff;
3200   image_hi &= 0xffffffff;
3201
3202   sign = (image_hi >> 31) & 1;
3203   exp = (image_hi >> 20) & 0x7ff;
3204
3205   memset (r, 0, sizeof (*r));
3206
3207   image_hi <<= 32 - 21;
3208   image_hi |= image_lo >> 21;
3209   image_hi &= 0x7fffffff;
3210   image_lo <<= 32 - 21;
3211
3212   if (exp == 0)
3213     {
3214       if ((image_hi || image_lo) && fmt->has_denorm)
3215         {
3216           r->cl = rvc_normal;
3217           r->sign = sign;
3218           SET_REAL_EXP (r, -1022);
3219           if (HOST_BITS_PER_LONG == 32)
3220             {
3221               image_hi = (image_hi << 1) | (image_lo >> 31);
3222               image_lo <<= 1;
3223               r->sig[SIGSZ-1] = image_hi;
3224               r->sig[SIGSZ-2] = image_lo;
3225             }
3226           else
3227             {
3228               image_hi = (image_hi << 31 << 2) | (image_lo << 1);
3229               r->sig[SIGSZ-1] = image_hi;
3230             }
3231           normalize (r);
3232         }
3233       else if (fmt->has_signed_zero)
3234         r->sign = sign;
3235     }
3236   else if (exp == 2047 && (fmt->has_nans || fmt->has_inf))
3237     {
3238       if (image_hi || image_lo)
3239         {
3240           r->cl = rvc_nan;
3241           r->sign = sign;
3242           r->signalling = ((image_hi >> 30) & 1) ^ fmt->qnan_msb_set;
3243           if (HOST_BITS_PER_LONG == 32)
3244             {
3245               r->sig[SIGSZ-1] = image_hi;
3246               r->sig[SIGSZ-2] = image_lo;
3247             }
3248           else
3249             r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo;
3250         }
3251       else
3252         {
3253           r->cl = rvc_inf;
3254           r->sign = sign;
3255         }
3256     }
3257   else
3258     {
3259       r->cl = rvc_normal;
3260       r->sign = sign;
3261       SET_REAL_EXP (r, exp - 1023 + 1);
3262       if (HOST_BITS_PER_LONG == 32)
3263         {
3264           r->sig[SIGSZ-1] = image_hi | SIG_MSB;
3265           r->sig[SIGSZ-2] = image_lo;
3266         }
3267       else
3268         r->sig[SIGSZ-1] = (image_hi << 31 << 1) | image_lo | SIG_MSB;
3269     }
3270 }
3271
3272 const struct real_format ieee_double_format =
3273   {
3274     encode_ieee_double,
3275     decode_ieee_double,
3276     2,
3277     53,
3278     53,
3279     -1021,
3280     1024,
3281     63,
3282     63,
3283     false,
3284     true,
3285     true,
3286     true,
3287     true,
3288     true,
3289     true,
3290     false
3291   };
3292
3293 const struct real_format mips_double_format =
3294   {
3295     encode_ieee_double,
3296     decode_ieee_double,
3297     2,
3298     53,
3299     53,
3300     -1021,
3301     1024,
3302     63,
3303     63,
3304     false,
3305     true,
3306     true,
3307     true,
3308     true,
3309     true,
3310     false,
3311     true
3312   };
3313
3314 const struct real_format motorola_double_format =
3315   {
3316     encode_ieee_double,
3317     decode_ieee_double,
3318     2,
3319     53,
3320     53,
3321     -1021,
3322     1024,
3323     63,
3324     63,
3325     false,
3326     true,
3327     true,
3328     true,
3329     true,
3330     true,
3331     true,
3332     true
3333   };
3334 \f
3335 /* IEEE extended real format.  This comes in three flavors: Intel's as
3336    a 12 byte image, Intel's as a 16 byte image, and Motorola's.  Intel
3337    12- and 16-byte images may be big- or little endian; Motorola's is
3338    always big endian.  */
3339
3340 /* Helper subroutine which converts from the internal format to the
3341    12-byte little-endian Intel format.  Functions below adjust this
3342    for the other possible formats.  */
3343 static void
3344 encode_ieee_extended (const struct real_format *fmt, long *buf,
3345                       const REAL_VALUE_TYPE *r)
3346 {
3347   unsigned long image_hi, sig_hi, sig_lo;
3348   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3349
3350   image_hi = r->sign << 15;
3351   sig_hi = sig_lo = 0;
3352
3353   switch (r->cl)
3354     {
3355     case rvc_zero:
3356       break;
3357
3358     case rvc_inf:
3359       if (fmt->has_inf)
3360         {
3361           image_hi |= 32767;
3362
3363           /* Intel requires the explicit integer bit to be set, otherwise
3364              it considers the value a "pseudo-infinity".  Motorola docs
3365              say it doesn't care.  */
3366           sig_hi = 0x80000000;
3367         }
3368       else
3369         {
3370           image_hi |= 32767;
3371           sig_lo = sig_hi = 0xffffffff;
3372         }
3373       break;
3374
3375     case rvc_nan:
3376       if (fmt->has_nans)
3377         {
3378           image_hi |= 32767;
3379           if (r->canonical)
3380             {
3381               if (fmt->canonical_nan_lsbs_set)
3382                 {
3383                   sig_hi = (1 << 30) - 1;
3384                   sig_lo = 0xffffffff;
3385                 }
3386             }
3387           else if (HOST_BITS_PER_LONG == 32)
3388             {
3389               sig_hi = r->sig[SIGSZ-1];
3390               sig_lo = r->sig[SIGSZ-2];
3391             }
3392           else
3393             {
3394               sig_lo = r->sig[SIGSZ-1];
3395               sig_hi = sig_lo >> 31 >> 1;
3396               sig_lo &= 0xffffffff;
3397             }
3398           if (r->signalling == fmt->qnan_msb_set)
3399             sig_hi &= ~(1 << 30);
3400           else
3401             sig_hi |= 1 << 30;
3402           if ((sig_hi & 0x7fffffff) == 0 && sig_lo == 0)
3403             sig_hi = 1 << 29;
3404
3405           /* Intel requires the explicit integer bit to be set, otherwise
3406              it considers the value a "pseudo-nan".  Motorola docs say it
3407              doesn't care.  */
3408           sig_hi |= 0x80000000;
3409         }
3410       else
3411         {
3412           image_hi |= 32767;
3413           sig_lo = sig_hi = 0xffffffff;
3414         }
3415       break;
3416
3417     case rvc_normal:
3418       {
3419         int exp = REAL_EXP (r);
3420
3421         /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3422            whereas the intermediate representation is 0.F x 2**exp.
3423            Which means we're off by one.
3424
3425            Except for Motorola, which consider exp=0 and explicit
3426            integer bit set to continue to be normalized.  In theory
3427            this discrepancy has been taken care of by the difference
3428            in fmt->emin in round_for_format.  */
3429
3430         if (denormal)
3431           exp = 0;
3432         else
3433           {
3434             exp += 16383 - 1;
3435             gcc_assert (exp >= 0);
3436           }
3437         image_hi |= exp;
3438
3439         if (HOST_BITS_PER_LONG == 32)
3440           {
3441             sig_hi = r->sig[SIGSZ-1];
3442             sig_lo = r->sig[SIGSZ-2];
3443           }
3444         else
3445           {
3446             sig_lo = r->sig[SIGSZ-1];
3447             sig_hi = sig_lo >> 31 >> 1;
3448             sig_lo &= 0xffffffff;
3449           }
3450       }
3451       break;
3452
3453     default:
3454       gcc_unreachable ();
3455     }
3456
3457   buf[0] = sig_lo, buf[1] = sig_hi, buf[2] = image_hi;
3458 }
3459
3460 /* Convert from the internal format to the 12-byte Motorola format
3461    for an IEEE extended real.  */
3462 static void
3463 encode_ieee_extended_motorola (const struct real_format *fmt, long *buf,
3464                                const REAL_VALUE_TYPE *r)
3465 {
3466   long intermed[3];
3467   encode_ieee_extended (fmt, intermed, r);
3468
3469   /* Motorola chips are assumed always to be big-endian.  Also, the
3470      padding in a Motorola extended real goes between the exponent and
3471      the mantissa.  At this point the mantissa is entirely within
3472      elements 0 and 1 of intermed, and the exponent entirely within
3473      element 2, so all we have to do is swap the order around, and
3474      shift element 2 left 16 bits.  */
3475   buf[0] = intermed[2] << 16;
3476   buf[1] = intermed[1];
3477   buf[2] = intermed[0];
3478 }
3479
3480 /* Convert from the internal format to the 12-byte Intel format for
3481    an IEEE extended real.  */
3482 static void
3483 encode_ieee_extended_intel_96 (const struct real_format *fmt, long *buf,
3484                                const REAL_VALUE_TYPE *r)
3485 {
3486   if (FLOAT_WORDS_BIG_ENDIAN)
3487     {
3488       /* All the padding in an Intel-format extended real goes at the high
3489          end, which in this case is after the mantissa, not the exponent.
3490          Therefore we must shift everything down 16 bits.  */
3491       long intermed[3];
3492       encode_ieee_extended (fmt, intermed, r);
3493       buf[0] = ((intermed[2] << 16) | ((unsigned long)(intermed[1] & 0xFFFF0000) >> 16));
3494       buf[1] = ((intermed[1] << 16) | ((unsigned long)(intermed[0] & 0xFFFF0000) >> 16));
3495       buf[2] =  (intermed[0] << 16);
3496     }
3497   else
3498     /* encode_ieee_extended produces what we want directly.  */
3499     encode_ieee_extended (fmt, buf, r);
3500 }
3501
3502 /* Convert from the internal format to the 16-byte Intel format for
3503    an IEEE extended real.  */
3504 static void
3505 encode_ieee_extended_intel_128 (const struct real_format *fmt, long *buf,
3506                                 const REAL_VALUE_TYPE *r)
3507 {
3508   /* All the padding in an Intel-format extended real goes at the high end.  */
3509   encode_ieee_extended_intel_96 (fmt, buf, r);
3510   buf[3] = 0;
3511 }
3512
3513 /* As above, we have a helper function which converts from 12-byte
3514    little-endian Intel format to internal format.  Functions below
3515    adjust for the other possible formats.  */
3516 static void
3517 decode_ieee_extended (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3518                       const long *buf)
3519 {
3520   unsigned long image_hi, sig_hi, sig_lo;
3521   bool sign;
3522   int exp;
3523
3524   sig_lo = buf[0], sig_hi = buf[1], image_hi = buf[2];
3525   sig_lo &= 0xffffffff;
3526   sig_hi &= 0xffffffff;
3527   image_hi &= 0xffffffff;
3528
3529   sign = (image_hi >> 15) & 1;
3530   exp = image_hi & 0x7fff;
3531
3532   memset (r, 0, sizeof (*r));
3533
3534   if (exp == 0)
3535     {
3536       if ((sig_hi || sig_lo) && fmt->has_denorm)
3537         {
3538           r->cl = rvc_normal;
3539           r->sign = sign;
3540
3541           /* When the IEEE format contains a hidden bit, we know that
3542              it's zero at this point, and so shift up the significand
3543              and decrease the exponent to match.  In this case, Motorola
3544              defines the explicit integer bit to be valid, so we don't
3545              know whether the msb is set or not.  */
3546           SET_REAL_EXP (r, fmt->emin);
3547           if (HOST_BITS_PER_LONG == 32)
3548             {
3549               r->sig[SIGSZ-1] = sig_hi;
3550               r->sig[SIGSZ-2] = sig_lo;
3551             }
3552           else
3553             r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3554
3555           normalize (r);
3556         }
3557       else if (fmt->has_signed_zero)
3558         r->sign = sign;
3559     }
3560   else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
3561     {
3562       /* See above re "pseudo-infinities" and "pseudo-nans".
3563          Short summary is that the MSB will likely always be
3564          set, and that we don't care about it.  */
3565       sig_hi &= 0x7fffffff;
3566
3567       if (sig_hi || sig_lo)
3568         {
3569           r->cl = rvc_nan;
3570           r->sign = sign;
3571           r->signalling = ((sig_hi >> 30) & 1) ^ fmt->qnan_msb_set;
3572           if (HOST_BITS_PER_LONG == 32)
3573             {
3574               r->sig[SIGSZ-1] = sig_hi;
3575               r->sig[SIGSZ-2] = sig_lo;
3576             }
3577           else
3578             r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3579         }
3580       else
3581         {
3582           r->cl = rvc_inf;
3583           r->sign = sign;
3584         }
3585     }
3586   else
3587     {
3588       r->cl = rvc_normal;
3589       r->sign = sign;
3590       SET_REAL_EXP (r, exp - 16383 + 1);
3591       if (HOST_BITS_PER_LONG == 32)
3592         {
3593           r->sig[SIGSZ-1] = sig_hi;
3594           r->sig[SIGSZ-2] = sig_lo;
3595         }
3596       else
3597         r->sig[SIGSZ-1] = (sig_hi << 31 << 1) | sig_lo;
3598     }
3599 }
3600
3601 /* Convert from the internal format to the 12-byte Motorola format
3602    for an IEEE extended real.  */
3603 static void
3604 decode_ieee_extended_motorola (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3605                                const long *buf)
3606 {
3607   long intermed[3];
3608
3609   /* Motorola chips are assumed always to be big-endian.  Also, the
3610      padding in a Motorola extended real goes between the exponent and
3611      the mantissa; remove it.  */
3612   intermed[0] = buf[2];
3613   intermed[1] = buf[1];
3614   intermed[2] = (unsigned long)buf[0] >> 16;
3615
3616   decode_ieee_extended (fmt, r, intermed);
3617 }
3618
3619 /* Convert from the internal format to the 12-byte Intel format for
3620    an IEEE extended real.  */
3621 static void
3622 decode_ieee_extended_intel_96 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3623                                const long *buf)
3624 {
3625   if (FLOAT_WORDS_BIG_ENDIAN)
3626     {
3627       /* All the padding in an Intel-format extended real goes at the high
3628          end, which in this case is after the mantissa, not the exponent.
3629          Therefore we must shift everything up 16 bits.  */
3630       long intermed[3];
3631
3632       intermed[0] = (((unsigned long)buf[2] >> 16) | (buf[1] << 16));
3633       intermed[1] = (((unsigned long)buf[1] >> 16) | (buf[0] << 16));
3634       intermed[2] =  ((unsigned long)buf[0] >> 16);
3635
3636       decode_ieee_extended (fmt, r, intermed);
3637     }
3638   else
3639     /* decode_ieee_extended produces what we want directly.  */
3640     decode_ieee_extended (fmt, r, buf);
3641 }
3642
3643 /* Convert from the internal format to the 16-byte Intel format for
3644    an IEEE extended real.  */
3645 static void
3646 decode_ieee_extended_intel_128 (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3647                                 const long *buf)
3648 {
3649   /* All the padding in an Intel-format extended real goes at the high end.  */
3650   decode_ieee_extended_intel_96 (fmt, r, buf);
3651 }
3652
3653 const struct real_format ieee_extended_motorola_format =
3654   {
3655     encode_ieee_extended_motorola,
3656     decode_ieee_extended_motorola,
3657     2,
3658     64,
3659     64,
3660     -16382,
3661     16384,
3662     95,
3663     95,
3664     false,
3665     true,
3666     true,
3667     true,
3668     true,
3669     true,
3670     true,
3671     true
3672   };
3673
3674 const struct real_format ieee_extended_intel_96_format =
3675   {
3676     encode_ieee_extended_intel_96,
3677     decode_ieee_extended_intel_96,
3678     2,
3679     64,
3680     64,
3681     -16381,
3682     16384,
3683     79,
3684     79,
3685     false,
3686     true,
3687     true,
3688     true,
3689     true,
3690     true,
3691     true,
3692     false
3693   };
3694
3695 const struct real_format ieee_extended_intel_128_format =
3696   {
3697     encode_ieee_extended_intel_128,
3698     decode_ieee_extended_intel_128,
3699     2,
3700     64,
3701     64,
3702     -16381,
3703     16384,
3704     79,
3705     79,
3706     false,
3707     true,
3708     true,
3709     true,
3710     true,
3711     true,
3712     true,
3713     false
3714   };
3715
3716 /* The following caters to i386 systems that set the rounding precision
3717    to 53 bits instead of 64, e.g. FreeBSD.  */
3718 const struct real_format ieee_extended_intel_96_round_53_format =
3719   {
3720     encode_ieee_extended_intel_96,
3721     decode_ieee_extended_intel_96,
3722     2,
3723     53,
3724     53,
3725     -16381,
3726     16384,
3727     79,
3728     79,
3729     false,
3730     true,
3731     true,
3732     true,
3733     true,
3734     true,
3735     true,
3736     false
3737   };
3738 \f
3739 /* IBM 128-bit extended precision format: a pair of IEEE double precision
3740    numbers whose sum is equal to the extended precision value.  The number
3741    with greater magnitude is first.  This format has the same magnitude
3742    range as an IEEE double precision value, but effectively 106 bits of
3743    significand precision.  Infinity and NaN are represented by their IEEE
3744    double precision value stored in the first number, the second number is
3745    +0.0 or -0.0 for Infinity and don't-care for NaN.  */
3746
3747 static void encode_ibm_extended (const struct real_format *fmt,
3748                                  long *, const REAL_VALUE_TYPE *);
3749 static void decode_ibm_extended (const struct real_format *,
3750                                  REAL_VALUE_TYPE *, const long *);
3751
3752 static void
3753 encode_ibm_extended (const struct real_format *fmt, long *buf,
3754                      const REAL_VALUE_TYPE *r)
3755 {
3756   REAL_VALUE_TYPE u, normr, v;
3757   const struct real_format *base_fmt;
3758
3759   base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3760
3761   /* Renormalize R before doing any arithmetic on it.  */
3762   normr = *r;
3763   if (normr.cl == rvc_normal)
3764     normalize (&normr);
3765
3766   /* u = IEEE double precision portion of significand.  */
3767   u = normr;
3768   round_for_format (base_fmt, &u);
3769   encode_ieee_double (base_fmt, &buf[0], &u);
3770
3771   if (u.cl == rvc_normal)
3772     {
3773       do_add (&v, &normr, &u, 1);
3774       /* Call round_for_format since we might need to denormalize.  */
3775       round_for_format (base_fmt, &v);
3776       encode_ieee_double (base_fmt, &buf[2], &v);
3777     }
3778   else
3779     {
3780       /* Inf, NaN, 0 are all representable as doubles, so the
3781          least-significant part can be 0.0.  */
3782       buf[2] = 0;
3783       buf[3] = 0;
3784     }
3785 }
3786
3787 static void
3788 decode_ibm_extended (const struct real_format *fmt ATTRIBUTE_UNUSED, REAL_VALUE_TYPE *r,
3789                      const long *buf)
3790 {
3791   REAL_VALUE_TYPE u, v;
3792   const struct real_format *base_fmt;
3793
3794   base_fmt = fmt->qnan_msb_set ? &ieee_double_format : &mips_double_format;
3795   decode_ieee_double (base_fmt, &u, &buf[0]);
3796
3797   if (u.cl != rvc_zero && u.cl != rvc_inf && u.cl != rvc_nan)
3798     {
3799       decode_ieee_double (base_fmt, &v, &buf[2]);
3800       do_add (r, &u, &v, 0);
3801     }
3802   else
3803     *r = u;
3804 }
3805
3806 const struct real_format ibm_extended_format =
3807   {
3808     encode_ibm_extended,
3809     decode_ibm_extended,
3810     2,
3811     53 + 53,
3812     53,
3813     -1021 + 53,
3814     1024,
3815     127,
3816     -1,
3817     false,
3818     true,
3819     true,
3820     true,
3821     true,
3822     true,
3823     true,
3824     false
3825   };
3826
3827 const struct real_format mips_extended_format =
3828   {
3829     encode_ibm_extended,
3830     decode_ibm_extended,
3831     2,
3832     53 + 53,
3833     53,
3834     -1021 + 53,
3835     1024,
3836     127,
3837     -1,
3838     false,
3839     true,
3840     true,
3841     true,
3842     true,
3843     true,
3844     false,
3845     true
3846   };
3847
3848 \f
3849 /* IEEE quad precision format.  */
3850
3851 static void encode_ieee_quad (const struct real_format *fmt,
3852                               long *, const REAL_VALUE_TYPE *);
3853 static void decode_ieee_quad (const struct real_format *,
3854                               REAL_VALUE_TYPE *, const long *);
3855
3856 static void
3857 encode_ieee_quad (const struct real_format *fmt, long *buf,
3858                   const REAL_VALUE_TYPE *r)
3859 {
3860   unsigned long image3, image2, image1, image0, exp;
3861   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
3862   REAL_VALUE_TYPE u;
3863
3864   image3 = r->sign << 31;
3865   image2 = 0;
3866   image1 = 0;
3867   image0 = 0;
3868
3869   rshift_significand (&u, r, SIGNIFICAND_BITS - 113);
3870
3871   switch (r->cl)
3872     {
3873     case rvc_zero:
3874       break;
3875
3876     case rvc_inf:
3877       if (fmt->has_inf)
3878         image3 |= 32767 << 16;
3879       else
3880         {
3881           image3 |= 0x7fffffff;
3882           image2 = 0xffffffff;
3883           image1 = 0xffffffff;
3884           image0 = 0xffffffff;
3885         }
3886       break;
3887
3888     case rvc_nan:
3889       if (fmt->has_nans)
3890         {
3891           image3 |= 32767 << 16;
3892
3893           if (r->canonical)
3894             {
3895               if (fmt->canonical_nan_lsbs_set)
3896                 {
3897                   image3 |= 0x7fff;
3898                   image2 = image1 = image0 = 0xffffffff;
3899                 }
3900             }
3901           else if (HOST_BITS_PER_LONG == 32)
3902             {
3903               image0 = u.sig[0];
3904               image1 = u.sig[1];
3905               image2 = u.sig[2];
3906               image3 |= u.sig[3] & 0xffff;
3907             }
3908           else
3909             {
3910               image0 = u.sig[0];
3911               image1 = image0 >> 31 >> 1;
3912               image2 = u.sig[1];
3913               image3 |= (image2 >> 31 >> 1) & 0xffff;
3914               image0 &= 0xffffffff;
3915               image2 &= 0xffffffff;
3916             }
3917           if (r->signalling == fmt->qnan_msb_set)
3918             image3 &= ~0x8000;
3919           else
3920             image3 |= 0x8000;
3921           if (((image3 & 0xffff) | image2 | image1 | image0) == 0)
3922             image3 |= 0x4000;
3923         }
3924       else
3925         {
3926           image3 |= 0x7fffffff;
3927           image2 = 0xffffffff;
3928           image1 = 0xffffffff;
3929           image0 = 0xffffffff;
3930         }
3931       break;
3932
3933     case rvc_normal:
3934       /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
3935          whereas the intermediate representation is 0.F x 2**exp.
3936          Which means we're off by one.  */
3937       if (denormal)
3938         exp = 0;
3939       else
3940         exp = REAL_EXP (r) + 16383 - 1;
3941       image3 |= exp << 16;
3942
3943       if (HOST_BITS_PER_LONG == 32)
3944         {
3945           image0 = u.sig[0];
3946           image1 = u.sig[1];
3947           image2 = u.sig[2];
3948           image3 |= u.sig[3] & 0xffff;
3949         }
3950       else
3951         {
3952           image0 = u.sig[0];
3953           image1 = image0 >> 31 >> 1;
3954           image2 = u.sig[1];
3955           image3 |= (image2 >> 31 >> 1) & 0xffff;
3956           image0 &= 0xffffffff;
3957           image2 &= 0xffffffff;
3958         }
3959       break;
3960
3961     default:
3962       gcc_unreachable ();
3963     }
3964
3965   if (FLOAT_WORDS_BIG_ENDIAN)
3966     {
3967       buf[0] = image3;
3968       buf[1] = image2;
3969       buf[2] = image1;
3970       buf[3] = image0;
3971     }
3972   else
3973     {
3974       buf[0] = image0;
3975       buf[1] = image1;
3976       buf[2] = image2;
3977       buf[3] = image3;
3978     }
3979 }
3980
3981 static void
3982 decode_ieee_quad (const struct real_format *fmt, REAL_VALUE_TYPE *r,
3983                   const long *buf)
3984 {
3985   unsigned long image3, image2, image1, image0;
3986   bool sign;
3987   int exp;
3988
3989   if (FLOAT_WORDS_BIG_ENDIAN)
3990     {
3991       image3 = buf[0];
3992       image2 = buf[1];
3993       image1 = buf[2];
3994       image0 = buf[3];
3995     }
3996   else
3997     {
3998       image0 = buf[0];
3999       image1 = buf[1];
4000       image2 = buf[2];
4001       image3 = buf[3];
4002     }
4003   image0 &= 0xffffffff;
4004   image1 &= 0xffffffff;
4005   image2 &= 0xffffffff;
4006
4007   sign = (image3 >> 31) & 1;
4008   exp = (image3 >> 16) & 0x7fff;
4009   image3 &= 0xffff;
4010
4011   memset (r, 0, sizeof (*r));
4012
4013   if (exp == 0)
4014     {
4015       if ((image3 | image2 | image1 | image0) && fmt->has_denorm)
4016         {
4017           r->cl = rvc_normal;
4018           r->sign = sign;
4019
4020           SET_REAL_EXP (r, -16382 + (SIGNIFICAND_BITS - 112));
4021           if (HOST_BITS_PER_LONG == 32)
4022             {
4023               r->sig[0] = image0;
4024               r->sig[1] = image1;
4025               r->sig[2] = image2;
4026               r->sig[3] = image3;
4027             }
4028           else
4029             {
4030               r->sig[0] = (image1 << 31 << 1) | image0;
4031               r->sig[1] = (image3 << 31 << 1) | image2;
4032             }
4033
4034           normalize (r);
4035         }
4036       else if (fmt->has_signed_zero)
4037         r->sign = sign;
4038     }
4039   else if (exp == 32767 && (fmt->has_nans || fmt->has_inf))
4040     {
4041       if (image3 | image2 | image1 | image0)
4042         {
4043           r->cl = rvc_nan;
4044           r->sign = sign;
4045           r->signalling = ((image3 >> 15) & 1) ^ fmt->qnan_msb_set;
4046
4047           if (HOST_BITS_PER_LONG == 32)
4048             {
4049               r->sig[0] = image0;
4050               r->sig[1] = image1;
4051               r->sig[2] = image2;
4052               r->sig[3] = image3;
4053             }
4054           else
4055             {
4056               r->sig[0] = (image1 << 31 << 1) | image0;
4057               r->sig[1] = (image3 << 31 << 1) | image2;
4058             }
4059           lshift_significand (r, r, SIGNIFICAND_BITS - 113);
4060         }
4061       else
4062         {
4063           r->cl = rvc_inf;
4064           r->sign = sign;
4065         }
4066     }
4067   else
4068     {
4069       r->cl = rvc_normal;
4070       r->sign = sign;
4071       SET_REAL_EXP (r, exp - 16383 + 1);
4072
4073       if (HOST_BITS_PER_LONG == 32)
4074         {
4075           r->sig[0] = image0;
4076           r->sig[1] = image1;
4077           r->sig[2] = image2;
4078           r->sig[3] = image3;
4079         }
4080       else
4081         {
4082           r->sig[0] = (image1 << 31 << 1) | image0;
4083           r->sig[1] = (image3 << 31 << 1) | image2;
4084         }
4085       lshift_significand (r, r, SIGNIFICAND_BITS - 113);
4086       r->sig[SIGSZ-1] |= SIG_MSB;
4087     }
4088 }
4089
4090 const struct real_format ieee_quad_format =
4091   {
4092     encode_ieee_quad,
4093     decode_ieee_quad,
4094     2,
4095     113,
4096     113,
4097     -16381,
4098     16384,
4099     127,
4100     127,
4101     false,
4102     true,
4103     true,
4104     true,
4105     true,
4106     true,
4107     true,
4108     false
4109   };
4110
4111 const struct real_format mips_quad_format =
4112   {
4113     encode_ieee_quad,
4114     decode_ieee_quad,
4115     2,
4116     113,
4117     113,
4118     -16381,
4119     16384,
4120     127,
4121     127,
4122     false,
4123     true,
4124     true,
4125     true,
4126     true,
4127     true,
4128     false,
4129     true
4130   };
4131 \f
4132 /* Descriptions of VAX floating point formats can be found beginning at
4133
4134    http://h71000.www7.hp.com/doc/73FINAL/4515/4515pro_013.html#f_floating_point_format
4135
4136    The thing to remember is that they're almost IEEE, except for word
4137    order, exponent bias, and the lack of infinities, nans, and denormals.
4138
4139    We don't implement the H_floating format here, simply because neither
4140    the VAX or Alpha ports use it.  */
4141
4142 static void encode_vax_f (const struct real_format *fmt,
4143                           long *, const REAL_VALUE_TYPE *);
4144 static void decode_vax_f (const struct real_format *,
4145                           REAL_VALUE_TYPE *, const long *);
4146 static void encode_vax_d (const struct real_format *fmt,
4147                           long *, const REAL_VALUE_TYPE *);
4148 static void decode_vax_d (const struct real_format *,
4149                           REAL_VALUE_TYPE *, const long *);
4150 static void encode_vax_g (const struct real_format *fmt,
4151                           long *, const REAL_VALUE_TYPE *);
4152 static void decode_vax_g (const struct real_format *,
4153                           REAL_VALUE_TYPE *, const long *);
4154
4155 static void
4156 encode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4157               const REAL_VALUE_TYPE *r)
4158 {
4159   unsigned long sign, exp, sig, image;
4160
4161   sign = r->sign << 15;
4162
4163   switch (r->cl)
4164     {
4165     case rvc_zero:
4166       image = 0;
4167       break;
4168
4169     case rvc_inf:
4170     case rvc_nan:
4171       image = 0xffff7fff | sign;
4172       break;
4173
4174     case rvc_normal:
4175       sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 24)) & 0x7fffff;
4176       exp = REAL_EXP (r) + 128;
4177
4178       image = (sig << 16) & 0xffff0000;
4179       image |= sign;
4180       image |= exp << 7;
4181       image |= sig >> 16;
4182       break;
4183
4184     default:
4185       gcc_unreachable ();
4186     }
4187
4188   buf[0] = image;
4189 }
4190
4191 static void
4192 decode_vax_f (const struct real_format *fmt ATTRIBUTE_UNUSED,
4193               REAL_VALUE_TYPE *r, const long *buf)
4194 {
4195   unsigned long image = buf[0] & 0xffffffff;
4196   int exp = (image >> 7) & 0xff;
4197
4198   memset (r, 0, sizeof (*r));
4199
4200   if (exp != 0)
4201     {
4202       r->cl = rvc_normal;
4203       r->sign = (image >> 15) & 1;
4204       SET_REAL_EXP (r, exp - 128);
4205
4206       image = ((image & 0x7f) << 16) | ((image >> 16) & 0xffff);
4207       r->sig[SIGSZ-1] = (image << (HOST_BITS_PER_LONG - 24)) | SIG_MSB;
4208     }
4209 }
4210
4211 static void
4212 encode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4213               const REAL_VALUE_TYPE *r)
4214 {
4215   unsigned long image0, image1, sign = r->sign << 15;
4216
4217   switch (r->cl)
4218     {
4219     case rvc_zero:
4220       image0 = image1 = 0;
4221       break;
4222
4223     case rvc_inf:
4224     case rvc_nan:
4225       image0 = 0xffff7fff | sign;
4226       image1 = 0xffffffff;
4227       break;
4228
4229     case rvc_normal:
4230       /* Extract the significand into straight hi:lo.  */
4231       if (HOST_BITS_PER_LONG == 64)
4232         {
4233           image0 = r->sig[SIGSZ-1];
4234           image1 = (image0 >> (64 - 56)) & 0xffffffff;
4235           image0 = (image0 >> (64 - 56 + 1) >> 31) & 0x7fffff;
4236         }
4237       else
4238         {
4239           image0 = r->sig[SIGSZ-1];
4240           image1 = r->sig[SIGSZ-2];
4241           image1 = (image0 << 24) | (image1 >> 8);
4242           image0 = (image0 >> 8) & 0xffffff;
4243         }
4244
4245       /* Rearrange the half-words of the significand to match the
4246          external format.  */
4247       image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff007f;
4248       image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
4249
4250       /* Add the sign and exponent.  */
4251       image0 |= sign;
4252       image0 |= (REAL_EXP (r) + 128) << 7;
4253       break;
4254
4255     default:
4256       gcc_unreachable ();
4257     }
4258
4259   if (FLOAT_WORDS_BIG_ENDIAN)
4260     buf[0] = image1, buf[1] = image0;
4261   else
4262     buf[0] = image0, buf[1] = image1;
4263 }
4264
4265 static void
4266 decode_vax_d (const struct real_format *fmt ATTRIBUTE_UNUSED,
4267               REAL_VALUE_TYPE *r, const long *buf)
4268 {
4269   unsigned long image0, image1;
4270   int exp;
4271
4272   if (FLOAT_WORDS_BIG_ENDIAN)
4273     image1 = buf[0], image0 = buf[1];
4274   else
4275     image0 = buf[0], image1 = buf[1];
4276   image0 &= 0xffffffff;
4277   image1 &= 0xffffffff;
4278
4279   exp = (image0 >> 7) & 0xff;
4280
4281   memset (r, 0, sizeof (*r));
4282
4283   if (exp != 0)
4284     {
4285       r->cl = rvc_normal;
4286       r->sign = (image0 >> 15) & 1;
4287       SET_REAL_EXP (r, exp - 128);
4288
4289       /* Rearrange the half-words of the external format into
4290          proper ascending order.  */
4291       image0 = ((image0 & 0x7f) << 16) | ((image0 >> 16) & 0xffff);
4292       image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
4293
4294       if (HOST_BITS_PER_LONG == 64)
4295         {
4296           image0 = (image0 << 31 << 1) | image1;
4297           image0 <<= 64 - 56;
4298           image0 |= SIG_MSB;
4299           r->sig[SIGSZ-1] = image0;
4300         }
4301       else
4302         {
4303           r->sig[SIGSZ-1] = image0;
4304           r->sig[SIGSZ-2] = image1;
4305           lshift_significand (r, r, 2*HOST_BITS_PER_LONG - 56);
4306           r->sig[SIGSZ-1] |= SIG_MSB;
4307         }
4308     }
4309 }
4310
4311 static void
4312 encode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4313               const REAL_VALUE_TYPE *r)
4314 {
4315   unsigned long image0, image1, sign = r->sign << 15;
4316
4317   switch (r->cl)
4318     {
4319     case rvc_zero:
4320       image0 = image1 = 0;
4321       break;
4322
4323     case rvc_inf:
4324     case rvc_nan:
4325       image0 = 0xffff7fff | sign;
4326       image1 = 0xffffffff;
4327       break;
4328
4329     case rvc_normal:
4330       /* Extract the significand into straight hi:lo.  */
4331       if (HOST_BITS_PER_LONG == 64)
4332         {
4333           image0 = r->sig[SIGSZ-1];
4334           image1 = (image0 >> (64 - 53)) & 0xffffffff;
4335           image0 = (image0 >> (64 - 53 + 1) >> 31) & 0xfffff;
4336         }
4337       else
4338         {
4339           image0 = r->sig[SIGSZ-1];
4340           image1 = r->sig[SIGSZ-2];
4341           image1 = (image0 << 21) | (image1 >> 11);
4342           image0 = (image0 >> 11) & 0xfffff;
4343         }
4344
4345       /* Rearrange the half-words of the significand to match the
4346          external format.  */
4347       image0 = ((image0 << 16) | (image0 >> 16)) & 0xffff000f;
4348       image1 = ((image1 << 16) | (image1 >> 16)) & 0xffffffff;
4349
4350       /* Add the sign and exponent.  */
4351       image0 |= sign;
4352       image0 |= (REAL_EXP (r) + 1024) << 4;
4353       break;
4354
4355     default:
4356       gcc_unreachable ();
4357     }
4358
4359   if (FLOAT_WORDS_BIG_ENDIAN)
4360     buf[0] = image1, buf[1] = image0;
4361   else
4362     buf[0] = image0, buf[1] = image1;
4363 }
4364
4365 static void
4366 decode_vax_g (const struct real_format *fmt ATTRIBUTE_UNUSED,
4367               REAL_VALUE_TYPE *r, const long *buf)
4368 {
4369   unsigned long image0, image1;
4370   int exp;
4371
4372   if (FLOAT_WORDS_BIG_ENDIAN)
4373     image1 = buf[0], image0 = buf[1];
4374   else
4375     image0 = buf[0], image1 = buf[1];
4376   image0 &= 0xffffffff;
4377   image1 &= 0xffffffff;
4378
4379   exp = (image0 >> 4) & 0x7ff;
4380
4381   memset (r, 0, sizeof (*r));
4382
4383   if (exp != 0)
4384     {
4385       r->cl = rvc_normal;
4386       r->sign = (image0 >> 15) & 1;
4387       SET_REAL_EXP (r, exp - 1024);
4388
4389       /* Rearrange the half-words of the external format into
4390          proper ascending order.  */
4391       image0 = ((image0 & 0xf) << 16) | ((image0 >> 16) & 0xffff);
4392       image1 = ((image1 & 0xffff) << 16) | ((image1 >> 16) & 0xffff);
4393
4394       if (HOST_BITS_PER_LONG == 64)
4395         {
4396           image0 = (image0 << 31 << 1) | image1;
4397           image0 <<= 64 - 53;
4398           image0 |= SIG_MSB;
4399           r->sig[SIGSZ-1] = image0;
4400         }
4401       else
4402         {
4403           r->sig[SIGSZ-1] = image0;
4404           r->sig[SIGSZ-2] = image1;
4405           lshift_significand (r, r, 64 - 53);
4406           r->sig[SIGSZ-1] |= SIG_MSB;
4407         }
4408     }
4409 }
4410
4411 const struct real_format vax_f_format =
4412   {
4413     encode_vax_f,
4414     decode_vax_f,
4415     2,
4416     24,
4417     24,
4418     -127,
4419     127,
4420     15,
4421     15,
4422     false,
4423     false,
4424     false,
4425     false,
4426     false,
4427     false,
4428     false,
4429     false
4430   };
4431
4432 const struct real_format vax_d_format =
4433   {
4434     encode_vax_d,
4435     decode_vax_d,
4436     2,
4437     56,
4438     56,
4439     -127,
4440     127,
4441     15,
4442     15,
4443     false,
4444     false,
4445     false,
4446     false,
4447     false,
4448     false,
4449     false,
4450     false
4451   };
4452
4453 const struct real_format vax_g_format =
4454   {
4455     encode_vax_g,
4456     decode_vax_g,
4457     2,
4458     53,
4459     53,
4460     -1023,
4461     1023,
4462     15,
4463     15,
4464     false,
4465     false,
4466     false,
4467     false,
4468     false,
4469     false,
4470     false,
4471     false
4472   };
4473 \f
4474 /* Encode real R into a single precision DFP value in BUF.  */
4475 static void
4476 encode_decimal_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4477                        long *buf ATTRIBUTE_UNUSED,
4478                        const REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED)
4479 {
4480   encode_decimal32 (fmt, buf, r);
4481 }
4482
4483 /* Decode a single precision DFP value in BUF into a real R.  */
4484 static void
4485 decode_decimal_single (const struct real_format *fmt ATTRIBUTE_UNUSED,
4486                        REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED,
4487                        const long *buf ATTRIBUTE_UNUSED)
4488 {
4489   decode_decimal32 (fmt, r, buf);
4490 }
4491
4492 /* Encode real R into a double precision DFP value in BUF.  */
4493 static void
4494 encode_decimal_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4495                        long *buf ATTRIBUTE_UNUSED,
4496                        const REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED)
4497 {
4498   encode_decimal64 (fmt, buf, r);
4499 }
4500
4501 /* Decode a double precision DFP value in BUF into a real R.  */
4502 static void
4503 decode_decimal_double (const struct real_format *fmt ATTRIBUTE_UNUSED,
4504                        REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED,
4505                        const long *buf ATTRIBUTE_UNUSED)
4506 {
4507   decode_decimal64 (fmt, r, buf);
4508 }
4509
4510 /* Encode real R into a quad precision DFP value in BUF.  */
4511 static void
4512 encode_decimal_quad (const struct real_format *fmt ATTRIBUTE_UNUSED,
4513                      long *buf ATTRIBUTE_UNUSED,
4514                      const REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED)
4515 {
4516   encode_decimal128 (fmt, buf, r);
4517 }
4518
4519 /* Decode a quad precision DFP value in BUF into a real R.  */
4520 static void
4521 decode_decimal_quad (const struct real_format *fmt ATTRIBUTE_UNUSED,
4522                      REAL_VALUE_TYPE *r ATTRIBUTE_UNUSED,
4523                      const long *buf ATTRIBUTE_UNUSED)
4524 {
4525   decode_decimal128 (fmt, r, buf);
4526 }
4527
4528 /* Single precision decimal floating point (IEEE 754). */
4529 const struct real_format decimal_single_format =
4530   {
4531     encode_decimal_single,
4532     decode_decimal_single,
4533     10,
4534     7,
4535     7,
4536     -94,
4537     97,
4538     31,
4539     31,
4540     false,
4541     true,
4542     true,
4543     true,
4544     true,
4545     true,
4546     true,
4547     false
4548   };
4549
4550 /* Double precision decimal floating point (IEEE 754). */
4551 const struct real_format decimal_double_format =
4552   {
4553     encode_decimal_double,
4554     decode_decimal_double,
4555     10,
4556     16,
4557     16,
4558     -382,
4559     385,
4560     63,
4561     63,
4562     false,
4563     true,
4564     true,
4565     true,
4566     true,
4567     true,
4568     true,
4569     false
4570   };
4571
4572 /* Quad precision decimal floating point (IEEE 754). */
4573 const struct real_format decimal_quad_format =
4574   {
4575     encode_decimal_quad,
4576     decode_decimal_quad,
4577     10,
4578     34,
4579     34,
4580     -6142,
4581     6145,
4582     127,
4583     127,
4584     false,
4585     true,
4586     true,
4587     true,
4588     true,
4589     true,
4590     true,
4591     false
4592   };
4593 \f
4594 /* Encode half-precision floats.  This routine is used both for the IEEE
4595    ARM alternative encodings.  */
4596 static void
4597 encode_ieee_half (const struct real_format *fmt, long *buf,
4598                   const REAL_VALUE_TYPE *r)
4599 {
4600   unsigned long image, sig, exp;
4601   unsigned long sign = r->sign;
4602   bool denormal = (r->sig[SIGSZ-1] & SIG_MSB) == 0;
4603
4604   image = sign << 15;
4605   sig = (r->sig[SIGSZ-1] >> (HOST_BITS_PER_LONG - 11)) & 0x3ff;
4606
4607   switch (r->cl)
4608     {
4609     case rvc_zero:
4610       break;
4611
4612     case rvc_inf:
4613       if (fmt->has_inf)
4614         image |= 31 << 10;
4615       else
4616         image |= 0x7fff;
4617       break;
4618
4619     case rvc_nan:
4620       if (fmt->has_nans)
4621         {
4622           if (r->canonical)
4623             sig = (fmt->canonical_nan_lsbs_set ? (1 << 9) - 1 : 0);
4624           if (r->signalling == fmt->qnan_msb_set)
4625             sig &= ~(1 << 9);
4626           else
4627             sig |= 1 << 9;
4628           if (sig == 0)
4629             sig = 1 << 8;
4630
4631           image |= 31 << 10;
4632           image |= sig;
4633         }
4634       else
4635         image |= 0x3ff;
4636       break;
4637
4638     case rvc_normal:
4639       /* Recall that IEEE numbers are interpreted as 1.F x 2**exp,
4640          whereas the intermediate representation is 0.F x 2**exp.
4641          Which means we're off by one.  */
4642       if (denormal)
4643         exp = 0;
4644       else
4645         exp = REAL_EXP (r) + 15 - 1;
4646       image |= exp << 10;
4647       image |= sig;
4648       break;
4649
4650     default:
4651       gcc_unreachable ();
4652     }
4653
4654   buf[0] = image;
4655 }
4656
4657 /* Decode half-precision floats.  This routine is used both for the IEEE
4658    ARM alternative encodings.  */
4659 static void
4660 decode_ieee_half (const struct real_format *fmt, REAL_VALUE_TYPE *r,
4661                   const long *buf)
4662 {
4663   unsigned long image = buf[0] & 0xffff;
4664   bool sign = (image >> 15) & 1;
4665   int exp = (image >> 10) & 0x1f;
4666
4667   memset (r, 0, sizeof (*r));
4668   image <<= HOST_BITS_PER_LONG - 11;
4669   image &= ~SIG_MSB;
4670
4671   if (exp == 0)
4672     {
4673       if (image && fmt->has_denorm)
4674         {
4675           r->cl = rvc_normal;
4676           r->sign = sign;
4677           SET_REAL_EXP (r, -14);
4678           r->sig[SIGSZ-1] = image << 1;
4679           normalize (r);
4680         }
4681       else if (fmt->has_signed_zero)
4682         r->sign = sign;
4683     }
4684   else if (exp == 31 && (fmt->has_nans || fmt->has_inf))
4685     {
4686       if (image)
4687         {
4688           r->cl = rvc_nan;
4689           r->sign = sign;
4690           r->signalling = (((image >> (HOST_BITS_PER_LONG - 2)) & 1)
4691                            ^ fmt->qnan_msb_set);
4692           r->sig[SIGSZ-1] = image;
4693         }
4694       else
4695         {
4696           r->cl = rvc_inf;
4697           r->sign = sign;
4698         }
4699     }
4700   else
4701     {
4702       r->cl = rvc_normal;
4703       r->sign = sign;
4704       SET_REAL_EXP (r, exp - 15 + 1);
4705       r->sig[SIGSZ-1] = image | SIG_MSB;
4706     }
4707 }
4708
4709 /* Half-precision format, as specified in IEEE 754R.  */
4710 const struct real_format ieee_half_format =
4711   {
4712     encode_ieee_half,
4713     decode_ieee_half,
4714     2,
4715     11,
4716     11,
4717     -13,
4718     16,
4719     15,
4720     15,
4721     false,
4722     true,
4723     true,
4724     true,
4725     true,
4726     true,
4727     true,
4728     false
4729   };
4730
4731 /* ARM's alternative half-precision format, similar to IEEE but with
4732    no reserved exponent value for NaNs and infinities; rather, it just
4733    extends the range of exponents by one.  */
4734 const struct real_format arm_half_format =
4735   {
4736     encode_ieee_half,
4737     decode_ieee_half,
4738     2,
4739     11,
4740     11,
4741     -13,
4742     17,
4743     15,
4744     15,
4745     false,
4746     true,
4747     false,
4748     false,
4749     true,
4750     true,
4751     false,
4752     false
4753   };
4754 \f
4755 /* A synthetic "format" for internal arithmetic.  It's the size of the
4756    internal significand minus the two bits needed for proper rounding.
4757    The encode and decode routines exist only to satisfy our paranoia
4758    harness.  */
4759
4760 static void encode_internal (const struct real_format *fmt,
4761                              long *, const REAL_VALUE_TYPE *);
4762 static void decode_internal (const struct real_format *,
4763                              REAL_VALUE_TYPE *, const long *);
4764
4765 static void
4766 encode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED, long *buf,
4767                  const REAL_VALUE_TYPE *r)
4768 {
4769   memcpy (buf, r, sizeof (*r));
4770 }
4771
4772 static void
4773 decode_internal (const struct real_format *fmt ATTRIBUTE_UNUSED,
4774                  REAL_VALUE_TYPE *r, const long *buf)
4775 {
4776   memcpy (r, buf, sizeof (*r));
4777 }
4778
4779 const struct real_format real_internal_format =
4780   {
4781     encode_internal,
4782     decode_internal,
4783     2,
4784     SIGNIFICAND_BITS - 2,
4785     SIGNIFICAND_BITS - 2,
4786     -MAX_EXP,
4787     MAX_EXP,
4788     -1,
4789     -1,
4790     false,
4791     false,
4792     true,
4793     true,
4794     false,
4795     true,
4796     true,
4797     false
4798   };
4799 \f
4800 /* Calculate the square root of X in mode MODE, and store the result
4801    in R.  Return TRUE if the operation does not raise an exception.
4802    For details see "High Precision Division and Square Root",
4803    Alan H. Karp and Peter Markstein, HP Lab Report 93-93-42, June
4804    1993.  http://www.hpl.hp.com/techreports/93/HPL-93-42.pdf.  */
4805
4806 bool
4807 real_sqrt (REAL_VALUE_TYPE *r, enum machine_mode mode,
4808            const REAL_VALUE_TYPE *x)
4809 {
4810   static REAL_VALUE_TYPE halfthree;
4811   static bool init = false;
4812   REAL_VALUE_TYPE h, t, i;
4813   int iter, exp;
4814
4815   /* sqrt(-0.0) is -0.0.  */
4816   if (real_isnegzero (x))
4817     {
4818       *r = *x;
4819       return false;
4820     }
4821
4822   /* Negative arguments return NaN.  */
4823   if (real_isneg (x))
4824     {
4825       get_canonical_qnan (r, 0);
4826       return false;
4827     }
4828
4829   /* Infinity and NaN return themselves.  */
4830   if (!real_isfinite (x))
4831     {
4832       *r = *x;
4833       return false;
4834     }
4835
4836   if (!init)
4837     {
4838       do_add (&halfthree, &dconst1, &dconsthalf, 0);
4839       init = true;
4840     }
4841
4842   /* Initial guess for reciprocal sqrt, i.  */
4843   exp = real_exponent (x);
4844   real_ldexp (&i, &dconst1, -exp/2);
4845
4846   /* Newton's iteration for reciprocal sqrt, i.  */
4847   for (iter = 0; iter < 16; iter++)
4848     {
4849       /* i(n+1) = i(n) * (1.5 - 0.5*i(n)*i(n)*x).  */
4850       do_multiply (&t, x, &i);
4851       do_multiply (&h, &t, &i);
4852       do_multiply (&t, &h, &dconsthalf);
4853       do_add (&h, &halfthree, &t, 1);
4854       do_multiply (&t, &i, &h);
4855
4856       /* Check for early convergence.  */
4857       if (iter >= 6 && real_identical (&i, &t))
4858         break;
4859
4860       /* ??? Unroll loop to avoid copying.  */
4861       i = t;
4862     }
4863
4864   /* Final iteration: r = i*x + 0.5*i*x*(1.0 - i*(i*x)).  */
4865   do_multiply (&t, x, &i);
4866   do_multiply (&h, &t, &i);
4867   do_add (&i, &dconst1, &h, 1);
4868   do_multiply (&h, &t, &i);
4869   do_multiply (&i, &dconsthalf, &h);
4870   do_add (&h, &t, &i, 0);
4871
4872   /* ??? We need a Tuckerman test to get the last bit.  */
4873
4874   real_convert (r, mode, &h);
4875   return true;
4876 }
4877
4878 /* Calculate X raised to the integer exponent N in mode MODE and store
4879    the result in R.  Return true if the result may be inexact due to
4880    loss of precision.  The algorithm is the classic "left-to-right binary
4881    method" described in section 4.6.3 of Donald Knuth's "Seminumerical
4882    Algorithms", "The Art of Computer Programming", Volume 2.  */
4883
4884 bool
4885 real_powi (REAL_VALUE_TYPE *r, enum machine_mode mode,
4886            const REAL_VALUE_TYPE *x, HOST_WIDE_INT n)
4887 {
4888   unsigned HOST_WIDE_INT bit;
4889   REAL_VALUE_TYPE t;
4890   bool inexact = false;
4891   bool init = false;
4892   bool neg;
4893   int i;
4894
4895   if (n == 0)
4896     {
4897       *r = dconst1;
4898       return false;
4899     }
4900   else if (n < 0)
4901     {
4902       /* Don't worry about overflow, from now on n is unsigned.  */
4903       neg = true;
4904       n = -n;
4905     }
4906   else
4907     neg = false;
4908
4909   t = *x;
4910   bit = (unsigned HOST_WIDE_INT) 1 << (HOST_BITS_PER_WIDE_INT - 1);
4911   for (i = 0; i < HOST_BITS_PER_WIDE_INT; i++)
4912     {
4913       if (init)
4914         {
4915           inexact |= do_multiply (&t, &t, &t);
4916           if (n & bit)
4917             inexact |= do_multiply (&t, &t, x);
4918         }
4919       else if (n & bit)
4920         init = true;
4921       bit >>= 1;
4922     }
4923
4924   if (neg)
4925     inexact |= do_divide (&t, &dconst1, &t);
4926
4927   real_convert (r, mode, &t);
4928   return inexact;
4929 }
4930
4931 /* Round X to the nearest integer not larger in absolute value, i.e.
4932    towards zero, placing the result in R in mode MODE.  */
4933
4934 void
4935 real_trunc (REAL_VALUE_TYPE *r, enum machine_mode mode,
4936             const REAL_VALUE_TYPE *x)
4937 {
4938   do_fix_trunc (r, x);
4939   if (mode != VOIDmode)
4940     real_convert (r, mode, r);
4941 }
4942
4943 /* Round X to the largest integer not greater in value, i.e. round
4944    down, placing the result in R in mode MODE.  */
4945
4946 void
4947 real_floor (REAL_VALUE_TYPE *r, enum machine_mode mode,
4948             const REAL_VALUE_TYPE *x)
4949 {
4950   REAL_VALUE_TYPE t;
4951
4952   do_fix_trunc (&t, x);
4953   if (! real_identical (&t, x) && x->sign)
4954     do_add (&t, &t, &dconstm1, 0);
4955   if (mode != VOIDmode)
4956     real_convert (r, mode, &t);
4957   else
4958     *r = t;
4959 }
4960
4961 /* Round X to the smallest integer not less then argument, i.e. round
4962    up, placing the result in R in mode MODE.  */
4963
4964 void
4965 real_ceil (REAL_VALUE_TYPE *r, enum machine_mode mode,
4966            const REAL_VALUE_TYPE *x)
4967 {
4968   REAL_VALUE_TYPE t;
4969
4970   do_fix_trunc (&t, x);
4971   if (! real_identical (&t, x) && ! x->sign)
4972     do_add (&t, &t, &dconst1, 0);
4973   if (mode != VOIDmode)
4974     real_convert (r, mode, &t);
4975   else
4976     *r = t;
4977 }
4978
4979 /* Round X to the nearest integer, but round halfway cases away from
4980    zero.  */
4981
4982 void
4983 real_round (REAL_VALUE_TYPE *r, enum machine_mode mode,
4984             const REAL_VALUE_TYPE *x)
4985 {
4986   do_add (r, x, &dconsthalf, x->sign);
4987   do_fix_trunc (r, r);
4988   if (mode != VOIDmode)
4989     real_convert (r, mode, r);
4990 }
4991
4992 /* Set the sign of R to the sign of X.  */
4993
4994 void
4995 real_copysign (REAL_VALUE_TYPE *r, const REAL_VALUE_TYPE *x)
4996 {
4997   r->sign = x->sign;
4998 }
4999
5000 /* Check whether the real constant value given is an integer.  */
5001
5002 bool
5003 real_isinteger (const REAL_VALUE_TYPE *c, enum machine_mode mode)
5004 {
5005   REAL_VALUE_TYPE cint;
5006
5007   real_trunc (&cint, mode, c);
5008   return real_identical (c, &cint);
5009 }
5010
5011 /* Write into BUF the maximum representable finite floating-point
5012    number, (1 - b**-p) * b**emax for a given FP format FMT as a hex
5013    float string.  LEN is the size of BUF, and the buffer must be large
5014    enough to contain the resulting string.  */
5015
5016 void
5017 get_max_float (const struct real_format *fmt, char *buf, size_t len)
5018 {
5019   int i, n;
5020   char *p;
5021
5022   strcpy (buf, "0x0.");
5023   n = fmt->p;
5024   for (i = 0, p = buf + 4; i + 3 < n; i += 4)
5025     *p++ = 'f';
5026   if (i < n)
5027     *p++ = "08ce"[n - i];
5028   sprintf (p, "p%d", fmt->emax);
5029   if (fmt->pnan < fmt->p)
5030     {
5031       /* This is an IBM extended double format made up of two IEEE
5032          doubles.  The value of the long double is the sum of the
5033          values of the two parts.  The most significant part is
5034          required to be the value of the long double rounded to the
5035          nearest double.  Rounding means we need a slightly smaller
5036          value for LDBL_MAX.  */
5037       buf[4 + fmt->pnan / 4] = "7bde"[fmt->pnan % 4];
5038     }
5039
5040   gcc_assert (strlen (buf) < len);
5041 }