gcc/real.c

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