gcc/real.c

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