gcc/real.c

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