gcc/real.c

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