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