gcc/real.c

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