gcc/real.c

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