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