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