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