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