gcc/real.c

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