lib/Support/APFloat.cpp

   1 //===-- APFloat.cpp - Implement APFloat class -----------------------------===//
   2 //
   3 //                     The LLVM Compiler Infrastructure
   4 //
   5 // This file is distributed under the University of Illinois Open Source
   6 // License. See LICENSE.TXT for details.
   7 //
   8 //===----------------------------------------------------------------------===//
   9 //
  10 // This file implements a class to represent arbitrary precision floating
  11 // point values and provide a variety of arithmetic operations on them.
  12 //
  13 //===----------------------------------------------------------------------===//
  14
  15 #include "llvm/ADT/APFloat.h"
  16 #include "llvm/ADT/FoldingSet.h"
  17 #include "llvm/Support/MathExtras.h"
  18 #include <cstring>
  19
  20 using namespace llvm;
  21
  22 #define convolve(lhs, rhs) ((lhs) * 4 + (rhs))
  23
  24 /* Assumed in hexadecimal significand parsing, and conversion to
  25    hexadecimal strings.  */
  26 #define COMPILE_TIME_ASSERT(cond) extern int CTAssert[(cond) ? 1 : -1]
  27 COMPILE_TIME_ASSERT(integerPartWidth % 4 == 0);
  28
  29 namespace llvm {
  30
  31   /* Represents floating point arithmetic semantics.  */
  32   struct fltSemantics {
  33     /* The largest E such that 2^E is representable; this matches the
  34        definition of IEEE 754.  */
  35     exponent_t maxExponent;
  36
  37     /* The smallest E such that 2^E is a normalized number; this
  38        matches the definition of IEEE 754.  */
  39     exponent_t minExponent;
  40
  41     /* Number of bits in the significand.  This includes the integer
  42        bit.  */
  43     unsigned int precision;
  44
  45     /* True if arithmetic is supported.  */
  46     unsigned int arithmeticOK;
  47   };
  48
  49   const fltSemantics APFloat::IEEEsingle = { 127, -126, 24, true };
  50   const fltSemantics APFloat::IEEEdouble = { 1023, -1022, 53, true };
  51   const fltSemantics APFloat::IEEEquad = { 16383, -16382, 113, true };
  52   const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64, true };
  53   const fltSemantics APFloat::Bogus = { 0, 0, 0, true };
  54
  55   // The PowerPC format consists of two doubles.  It does not map cleanly
  56   // onto the usual format above.  For now only storage of constants of
  57   // this type is supported, no arithmetic.
  58   const fltSemantics APFloat::PPCDoubleDouble = { 1023, -1022, 106, false };
  59
  60   /* A tight upper bound on number of parts required to hold the value
  61      pow(5, power) is
  62
  63        power * 815 / (351 * integerPartWidth) + 1
  64
  65      However, whilst the result may require only this many parts,
  66      because we are multiplying two values to get it, the
  67      multiplication may require an extra part with the excess part
  68      being zero (consider the trivial case of 1 * 1, tcFullMultiply
  69      requires two parts to hold the single-part result).  So we add an
  70      extra one to guarantee enough space whilst multiplying.  */
  71   const unsigned int maxExponent = 16383;
  72   const unsigned int maxPrecision = 113;
  73   const unsigned int maxPowerOfFiveExponent = maxExponent + maxPrecision - 1;
  74   const unsigned int maxPowerOfFiveParts = 2 + ((maxPowerOfFiveExponent * 815)
  75                                                 / (351 * integerPartWidth));
  76 }
  77
  78 /* Put a bunch of private, handy routines in an anonymous namespace.  */
  79 namespace {
  80
  81   static inline unsigned int
  82   partCountForBits(unsigned int bits)
  83   {
  84     return ((bits) + integerPartWidth - 1) / integerPartWidth;
  85   }
  86
  87   /* Returns 0U-9U.  Return values >= 10U are not digits.  */
  88   static inline unsigned int
  89   decDigitValue(unsigned int c)
  90   {
  91     return c - '0';
  92   }
  93
  94   static unsigned int
  95   hexDigitValue(unsigned int c)
  96   {
  97     unsigned int r;
  98
  99     r = c - '0';
 100     if(r <= 9)
 101       return r;
 102
 103     r = c - 'A';
 104     if(r <= 5)
 105       return r + 10;
 106
 107     r = c - 'a';
 108     if(r <= 5)
 109       return r + 10;
 110
 111     return -1U;
 112   }
 113
 114   static inline void
 115   assertArithmeticOK(const llvm::fltSemantics &semantics) {
 116     assert(semantics.arithmeticOK
 117            && "Compile-time arithmetic does not support these semantics");
 118   }
 119
 120   /* Return the value of a decimal exponent of the form
 121      [+-]ddddddd.
 122
 123      If the exponent overflows, returns a large exponent with the
 124      appropriate sign.  */
 125   static int
 126   readExponent(const char *p)
 127   {
 128     bool isNegative;
 129     unsigned int absExponent;
 130     const unsigned int overlargeExponent = 24000;  /* FIXME.  */
 131
 132     isNegative = (*p == '-');
 133     if (*p == '-' || *p == '+')
 134       p++;
 135
 136     absExponent = decDigitValue(*p++);
 137     assert (absExponent < 10U);
 138
 139     for (;;) {
 140       unsigned int value;
 141
 142       value = decDigitValue(*p);
 143       if (value >= 10U)
 144         break;
 145
 146       p++;
 147       value += absExponent * 10;
 148       if (absExponent >= overlargeExponent) {
 149         absExponent = overlargeExponent;
 150         break;
 151       }
 152       absExponent = value;
 153     }
 154
 155     if (isNegative)
 156       return -(int) absExponent;
 157     else
 158       return (int) absExponent;
 159   }
 160
 161   /* This is ugly and needs cleaning up, but I don't immediately see
 162      how whilst remaining safe.  */
 163   static int
 164   totalExponent(const char *p, int exponentAdjustment)
 165   {
 166     int unsignedExponent;
 167     bool negative, overflow;
 168     int exponent;
 169
 170     /* Move past the exponent letter and sign to the digits.  */
 171     p++;
 172     negative = *p == '-';
 173     if(*p == '-' || *p == '+')
 174       p++;
 175
 176     unsignedExponent = 0;
 177     overflow = false;
 178     for(;;) {
 179       unsigned int value;
 180
 181       value = decDigitValue(*p);
 182       if(value >= 10U)
 183         break;
 184
 185       p++;
 186       unsignedExponent = unsignedExponent * 10 + value;
 187       if(unsignedExponent > 65535)
 188         overflow = true;
 189     }
 190
 191     if(exponentAdjustment > 65535 || exponentAdjustment < -65536)
 192       overflow = true;
 193
 194     if(!overflow) {
 195       exponent = unsignedExponent;
 196       if(negative)
 197         exponent = -exponent;
 198       exponent += exponentAdjustment;
 199       if(exponent > 65535 || exponent < -65536)
 200         overflow = true;
 201     }
 202
 203     if(overflow)
 204       exponent = negative ? -65536: 65535;
 205
 206     return exponent;
 207   }
 208
 209   static const char *
 210   skipLeadingZeroesAndAnyDot(const char *p, const char **dot)
 211   {
 212     *dot = 0;
 213     while(*p == '0')
 214       p++;
 215
 216     if(*p == '.') {
 217       *dot = p++;
 218       while(*p == '0')
 219         p++;
 220     }
 221
 222     return p;
 223   }
 224
 225   /* Given a normal decimal floating point number of the form
 226
 227        dddd.dddd[eE][+-]ddd
 228
 229      where the decimal point and exponent are optional, fill out the
 230      structure D.  Exponent is appropriate if the significand is
 231      treated as an integer, and normalizedExponent if the significand
 232      is taken to have the decimal point after a single leading
 233      non-zero digit.
 234
 235      If the value is zero, V->firstSigDigit points to a non-digit, and
 236      the return exponent is zero.
 237   */
 238   struct decimalInfo {
 239     const char *firstSigDigit;
 240     const char *lastSigDigit;
 241     int exponent;
 242     int normalizedExponent;
 243   };
 244
 245   static void
 246   interpretDecimal(const char *p, decimalInfo *D)
 247   {
 248     const char *dot;
 249
 250     p = skipLeadingZeroesAndAnyDot (p, &dot);
 251
 252     D->firstSigDigit = p;
 253     D->exponent = 0;
 254     D->normalizedExponent = 0;
 255
 256     for (;;) {
 257       if (*p == '.') {
 258         assert(dot == 0);
 259         dot = p++;
 260       }
 261       if (decDigitValue(*p) >= 10U)
 262         break;
 263       p++;
 264     }
 265
 266     /* If number is all zerooes accept any exponent.  */
 267     if (p != D->firstSigDigit) {
 268       if (*p == 'e' || *p == 'E')
 269         D->exponent = readExponent(p + 1);
 270
 271       /* Implied decimal point?  */
 272       if (!dot)
 273         dot = p;
 274
 275       /* Drop insignificant trailing zeroes.  */
 276       do
 277         do
 278           p--;
 279         while (*p == '0');
 280       while (*p == '.');
 281
 282       /* Adjust the exponents for any decimal point.  */
 283       D->exponent += static_cast<exponent_t>((dot - p) - (dot > p));
 284       D->normalizedExponent = (D->exponent +
 285                 static_cast<exponent_t>((p - D->firstSigDigit)
 286                                         - (dot > D->firstSigDigit && dot < p)));
 287     }
 288
 289     D->lastSigDigit = p;
 290   }
 291
 292   /* Return the trailing fraction of a hexadecimal number.
 293      DIGITVALUE is the first hex digit of the fraction, P points to
 294      the next digit.  */
 295   static lostFraction
 296   trailingHexadecimalFraction(const char *p, unsigned int digitValue)
 297   {
 298     unsigned int hexDigit;
 299
 300     /* If the first trailing digit isn't 0 or 8 we can work out the
 301        fraction immediately.  */
 302     if(digitValue > 8)
 303       return lfMoreThanHalf;
 304     else if(digitValue < 8 && digitValue > 0)
 305       return lfLessThanHalf;
 306
 307     /* Otherwise we need to find the first non-zero digit.  */
 308     while(*p == '0')
 309       p++;
 310
 311     hexDigit = hexDigitValue(*p);
 312
 313     /* If we ran off the end it is exactly zero or one-half, otherwise
 314        a little more.  */
 315     if(hexDigit == -1U)
 316       return digitValue == 0 ? lfExactlyZero: lfExactlyHalf;
 317     else
 318       return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf;
 319   }
 320
 321   /* Return the fraction lost were a bignum truncated losing the least
 322      significant BITS bits.  */
 323   static lostFraction
 324   lostFractionThroughTruncation(const integerPart *parts,
 325                                 unsigned int partCount,
 326                                 unsigned int bits)
 327   {
 328     unsigned int lsb;
 329
 330     lsb = APInt::tcLSB(parts, partCount);
 331
 332     /* Note this is guaranteed true if bits == 0, or LSB == -1U.  */
 333     if(bits <= lsb)
 334       return lfExactlyZero;
 335     if(bits == lsb + 1)
 336       return lfExactlyHalf;
 337     if(bits <= partCount * integerPartWidth
 338        && APInt::tcExtractBit(parts, bits - 1))
 339       return lfMoreThanHalf;
 340
 341     return lfLessThanHalf;
 342   }
 343
 344   /* Shift DST right BITS bits noting lost fraction.  */
 345   static lostFraction
 346   shiftRight(integerPart *dst, unsigned int parts, unsigned int bits)
 347   {
 348     lostFraction lost_fraction;
 349
 350     lost_fraction = lostFractionThroughTruncation(dst, parts, bits);
 351
 352     APInt::tcShiftRight(dst, parts, bits);
 353
 354     return lost_fraction;
 355   }
 356
 357   /* Combine the effect of two lost fractions.  */
 358   static lostFraction
 359   combineLostFractions(lostFraction moreSignificant,
 360                        lostFraction lessSignificant)
 361   {
 362     if(lessSignificant != lfExactlyZero) {
 363       if(moreSignificant == lfExactlyZero)
 364         moreSignificant = lfLessThanHalf;
 365       else if(moreSignificant == lfExactlyHalf)
 366         moreSignificant = lfMoreThanHalf;
 367     }
 368
 369     return moreSignificant;
 370   }
 371
 372   /* The error from the true value, in half-ulps, on multiplying two
 373      floating point numbers, which differ from the value they
 374      approximate by at most HUE1 and HUE2 half-ulps, is strictly less
 375      than the returned value.
 376
 377      See "How to Read Floating Point Numbers Accurately" by William D
 378      Clinger.  */
 379   static unsigned int
 380   HUerrBound(bool inexactMultiply, unsigned int HUerr1, unsigned int HUerr2)
 381   {
 382     assert(HUerr1 < 2 || HUerr2 < 2 || (HUerr1 + HUerr2 < 8));
 383
 384     if (HUerr1 + HUerr2 == 0)
 385       return inexactMultiply * 2;  /* <= inexactMultiply half-ulps.  */
 386     else
 387       return inexactMultiply + 2 * (HUerr1 + HUerr2);
 388   }
 389
 390   /* The number of ulps from the boundary (zero, or half if ISNEAREST)
 391      when the least significant BITS are truncated.  BITS cannot be
 392      zero.  */
 393   static integerPart
 394   ulpsFromBoundary(const integerPart *parts, unsigned int bits, bool isNearest)
 395   {
 396     unsigned int count, partBits;
 397     integerPart part, boundary;
 398
 399     assert (bits != 0);
 400
 401     bits--;
 402     count = bits / integerPartWidth;
 403     partBits = bits % integerPartWidth + 1;
 404
 405     part = parts[count] & (~(integerPart) 0 >> (integerPartWidth - partBits));
 406
 407     if (isNearest)
 408       boundary = (integerPart) 1 << (partBits - 1);
 409     else
 410       boundary = 0;
 411
 412     if (count == 0) {
 413       if (part - boundary <= boundary - part)
 414         return part - boundary;
 415       else
 416         return boundary - part;
 417     }
 418
 419     if (part == boundary) {
 420       while (--count)
 421         if (parts[count])
 422           return ~(integerPart) 0; /* A lot.  */
 423
 424       return parts[0];
 425     } else if (part == boundary - 1) {
 426       while (--count)
 427         if (~parts[count])
 428           return ~(integerPart) 0; /* A lot.  */
 429
 430       return -parts[0];
 431     }
 432
 433     return ~(integerPart) 0; /* A lot.  */
 434   }
 435
 436   /* Place pow(5, power) in DST, and return the number of parts used.
 437      DST must be at least one part larger than size of the answer.  */
 438   static unsigned int
 439   powerOf5(integerPart *dst, unsigned int power)
 440   {
 441     static const integerPart firstEightPowers[] = { 1, 5, 25, 125, 625, 3125,
 442                                                     15625, 78125 };
 443     static integerPart pow5s[maxPowerOfFiveParts * 2 + 5] = { 78125 * 5 };
 444     static unsigned int partsCount[16] = { 1 };
 445
 446     integerPart scratch[maxPowerOfFiveParts], *p1, *p2, *pow5;
 447     unsigned int result;
 448
 449     assert(power <= maxExponent);
 450
 451     p1 = dst;
 452     p2 = scratch;
 453
 454     *p1 = firstEightPowers[power & 7];
 455     power >>= 3;
 456
 457     result = 1;
 458     pow5 = pow5s;
 459
 460     for (unsigned int n = 0; power; power >>= 1, n++) {
 461       unsigned int pc;
 462
 463       pc = partsCount[n];
 464
 465       /* Calculate pow(5,pow(2,n+3)) if we haven't yet.  */
 466       if (pc == 0) {
 467         pc = partsCount[n - 1];
 468         APInt::tcFullMultiply(pow5, pow5 - pc, pow5 - pc, pc, pc);
 469         pc *= 2;
 470         if (pow5[pc - 1] == 0)
 471           pc--;
 472         partsCount[n] = pc;
 473       }
 474
 475       if (power & 1) {
 476         integerPart *tmp;
 477
 478         APInt::tcFullMultiply(p2, p1, pow5, result, pc);
 479         result += pc;
 480         if (p2[result - 1] == 0)
 481           result--;
 482
 483         /* Now result is in p1 with partsCount parts and p2 is scratch
 484            space.  */
 485         tmp = p1, p1 = p2, p2 = tmp;
 486       }
 487
 488       pow5 += pc;
 489     }
 490
 491     if (p1 != dst)
 492       APInt::tcAssign(dst, p1, result);
 493
 494     return result;
 495   }
 496
 497   /* Zero at the end to avoid modular arithmetic when adding one; used
 498      when rounding up during hexadecimal output.  */
 499   static const char hexDigitsLower[] = "0123456789abcdef0";
 500   static const char hexDigitsUpper[] = "0123456789ABCDEF0";
 501   static const char infinityL[] = "infinity";
 502   static const char infinityU[] = "INFINITY";
 503   static const char NaNL[] = "nan";
 504   static const char NaNU[] = "NAN";
 505
 506   /* Write out an integerPart in hexadecimal, starting with the most
 507      significant nibble.  Write out exactly COUNT hexdigits, return
 508      COUNT.  */
 509   static unsigned int
 510   partAsHex (char *dst, integerPart part, unsigned int count,
 511              const char *hexDigitChars)
 512   {
 513     unsigned int result = count;
 514
 515     assert (count != 0 && count <= integerPartWidth / 4);
 516
 517     part >>= (integerPartWidth - 4 * count);
 518     while (count--) {
 519       dst[count] = hexDigitChars[part & 0xf];
 520       part >>= 4;
 521     }
 522
 523     return result;
 524   }
 525
 526   /* Write out an unsigned decimal integer.  */
 527   static char *
 528   writeUnsignedDecimal (char *dst, unsigned int n)
 529   {
 530     char buff[40], *p;
 531
 532     p = buff;
 533     do
 534       *p++ = '0' + n % 10;
 535     while (n /= 10);
 536
 537     do
 538       *dst++ = *--p;
 539     while (p != buff);
 540
 541     return dst;
 542   }
 543
 544   /* Write out a signed decimal integer.  */
 545   static char *
 546   writeSignedDecimal (char *dst, int value)
 547   {
 548     if (value < 0) {
 549       *dst++ = '-';
 550       dst = writeUnsignedDecimal(dst, -(unsigned) value);
 551     } else
 552       dst = writeUnsignedDecimal(dst, value);
 553
 554     return dst;
 555   }
 556 }
 557
 558 /* Constructors.  */
 559 void
 560 APFloat::initialize(const fltSemantics *ourSemantics)
 561 {
 562   unsigned int count;
 563
 564   semantics = ourSemantics;
 565   count = partCount();
 566   if(count > 1)
 567     significand.parts = new integerPart[count];
 568 }
 569
 570 void
 571 APFloat::freeSignificand()
 572 {
 573   if(partCount() > 1)
 574     delete [] significand.parts;
 575 }
 576
 577 void
 578 APFloat::assign(const APFloat &rhs)
 579 {
 580   assert(semantics == rhs.semantics);
 581
 582   sign = rhs.sign;
 583   category = rhs.category;
 584   exponent = rhs.exponent;
 585   sign2 = rhs.sign2;
 586   exponent2 = rhs.exponent2;
 587   if(category == fcNormal || category == fcNaN)
 588     copySignificand(rhs);
 589 }
 590
 591 void
 592 APFloat::copySignificand(const APFloat &rhs)
 593 {
 594   assert(category == fcNormal || category == fcNaN);
 595   assert(rhs.partCount() >= partCount());
 596
 597   APInt::tcAssign(significandParts(), rhs.significandParts(),
 598                   partCount());
 599 }
 600
 601 /* Make this number a NaN, with an arbitrary but deterministic value
 602    for the significand.  */
 603 void
 604 APFloat::makeNaN(void)
 605 {
 606   category = fcNaN;
 607   APInt::tcSet(significandParts(), ~0U, partCount());
 608 }
 609
 610 APFloat &
 611 APFloat::operator=(const APFloat &rhs)
 612 {
 613   if(this != &rhs) {
 614     if(semantics != rhs.semantics) {
 615       freeSignificand();
 616       initialize(rhs.semantics);
 617     }
 618     assign(rhs);
 619   }
 620
 621   return *this;
 622 }
 623
 624 bool
 625 APFloat::bitwiseIsEqual(const APFloat &rhs) const {
 626   if (this == &rhs)
 627     return true;
 628   if (semantics != rhs.semantics ||
 629       category != rhs.category ||
 630       sign != rhs.sign)
 631     return false;
 632   if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble &&
 633       sign2 != rhs.sign2)
 634     return false;
 635   if (category==fcZero || category==fcInfinity)
 636     return true;
 637   else if (category==fcNormal && exponent!=rhs.exponent)
 638     return false;
 639   else if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble &&
 640            exponent2!=rhs.exponent2)
 641     return false;
 642   else {
 643     int i= partCount();
 644     const integerPart* p=significandParts();
 645     const integerPart* q=rhs.significandParts();
 646     for (; i>0; i--, p++, q++) {
 647       if (*p != *q)
 648         return false;
 649     }
 650     return true;
 651   }
 652 }
 653
 654 APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value)
 655 {
 656   assertArithmeticOK(ourSemantics);
 657   initialize(&ourSemantics);
 658   sign = 0;
 659   zeroSignificand();
 660   exponent = ourSemantics.precision - 1;
 661   significandParts()[0] = value;
 662   normalize(rmNearestTiesToEven, lfExactlyZero);
 663 }
 664
 665 APFloat::APFloat(const fltSemantics &ourSemantics,
 666                  fltCategory ourCategory, bool negative)
 667 {
 668   assertArithmeticOK(ourSemantics);
 669   initialize(&ourSemantics);
 670   category = ourCategory;
 671   sign = negative;
 672   if(category == fcNormal)
 673     category = fcZero;
 674   else if (ourCategory == fcNaN)
 675     makeNaN();
 676 }
 677
 678 APFloat::APFloat(const fltSemantics &ourSemantics, const char *text)
 679 {
 680   assertArithmeticOK(ourSemantics);
 681   initialize(&ourSemantics);
 682   convertFromString(text, rmNearestTiesToEven);
 683 }
 684
 685 APFloat::APFloat(const APFloat &rhs)
 686 {
 687   initialize(rhs.semantics);
 688   assign(rhs);
 689 }
 690
 691 APFloat::~APFloat()
 692 {
 693   freeSignificand();
 694 }
 695
 696 // Profile - This method 'profiles' an APFloat for use with FoldingSet.
 697 void APFloat::Profile(FoldingSetNodeID& ID) const {
 698   ID.Add(bitcastToAPInt());
 699 }
 700
 701 unsigned int
 702 APFloat::partCount() const
 703 {
 704   return partCountForBits(semantics->precision + 1);
 705 }
 706
 707 unsigned int
 708 APFloat::semanticsPrecision(const fltSemantics &semantics)
 709 {
 710   return semantics.precision;
 711 }
 712
 713 const integerPart *
 714 APFloat::significandParts() const
 715 {
 716   return const_cast<APFloat *>(this)->significandParts();
 717 }
 718
 719 integerPart *
 720 APFloat::significandParts()
 721 {
 722   assert(category == fcNormal || category == fcNaN);
 723
 724   if(partCount() > 1)
 725     return significand.parts;
 726   else
 727     return &significand.part;
 728 }
 729
 730 void
 731 APFloat::zeroSignificand()
 732 {
 733   category = fcNormal;
 734   APInt::tcSet(significandParts(), 0, partCount());
 735 }
 736
 737 /* Increment an fcNormal floating point number's significand.  */
 738 void
 739 APFloat::incrementSignificand()
 740 {
 741   integerPart carry;
 742
 743   carry = APInt::tcIncrement(significandParts(), partCount());
 744
 745   /* Our callers should never cause us to overflow.  */
 746   assert(carry == 0);
 747 }
 748
 749 /* Add the significand of the RHS.  Returns the carry flag.  */
 750 integerPart
 751 APFloat::addSignificand(const APFloat &rhs)
 752 {
 753   integerPart *parts;
 754
 755   parts = significandParts();
 756
 757   assert(semantics == rhs.semantics);
 758   assert(exponent == rhs.exponent);
 759
 760   return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
 761 }
 762
 763 /* Subtract the significand of the RHS with a borrow flag.  Returns
 764    the borrow flag.  */
 765 integerPart
 766 APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
 767 {
 768   integerPart *parts;
 769
 770   parts = significandParts();
 771
 772   assert(semantics == rhs.semantics);
 773   assert(exponent == rhs.exponent);
 774
 775   return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
 776                            partCount());
 777 }
 778
 779 /* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
 780    on to the full-precision result of the multiplication.  Returns the
 781    lost fraction.  */
 782 lostFraction
 783 APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
 784 {
 785   unsigned int omsb;        // One, not zero, based MSB.
 786   unsigned int partsCount, newPartsCount, precision;
 787   integerPart *lhsSignificand;
 788   integerPart scratch[4];
 789   integerPart *fullSignificand;
 790   lostFraction lost_fraction;
 791   bool ignored;
 792
 793   assert(semantics == rhs.semantics);
 794
 795   precision = semantics->precision;
 796   newPartsCount = partCountForBits(precision * 2);
 797
 798   if(newPartsCount > 4)
 799     fullSignificand = new integerPart[newPartsCount];
 800   else
 801     fullSignificand = scratch;
 802
 803   lhsSignificand = significandParts();
 804   partsCount = partCount();
 805
 806   APInt::tcFullMultiply(fullSignificand, lhsSignificand,
 807                         rhs.significandParts(), partsCount, partsCount);
 808
 809   lost_fraction = lfExactlyZero;
 810   omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
 811   exponent += rhs.exponent;
 812
 813   if(addend) {
 814     Significand savedSignificand = significand;
 815     const fltSemantics *savedSemantics = semantics;
 816     fltSemantics extendedSemantics;
 817     opStatus status;
 818     unsigned int extendedPrecision;
 819
 820     /* Normalize our MSB.  */
 821     extendedPrecision = precision + precision - 1;
 822     if(omsb != extendedPrecision)
 823       {
 824         APInt::tcShiftLeft(fullSignificand, newPartsCount,
 825                            extendedPrecision - omsb);
 826         exponent -= extendedPrecision - omsb;
 827       }
 828
 829     /* Create new semantics.  */
 830     extendedSemantics = *semantics;
 831     extendedSemantics.precision = extendedPrecision;
 832
 833     if(newPartsCount == 1)
 834       significand.part = fullSignificand[0];
 835     else
 836       significand.parts = fullSignificand;
 837     semantics = &extendedSemantics;
 838
 839     APFloat extendedAddend(*addend);
 840     status = extendedAddend.convert(extendedSemantics, rmTowardZero, &ignored);
 841     assert(status == opOK);
 842     lost_fraction = addOrSubtractSignificand(extendedAddend, false);
 843
 844     /* Restore our state.  */
 845     if(newPartsCount == 1)
 846       fullSignificand[0] = significand.part;
 847     significand = savedSignificand;
 848     semantics = savedSemantics;
 849
 850     omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
 851   }
 852
 853   exponent -= (precision - 1);
 854
 855   if(omsb > precision) {
 856     unsigned int bits, significantParts;
 857     lostFraction lf;
 858
 859     bits = omsb - precision;
 860     significantParts = partCountForBits(omsb);
 861     lf = shiftRight(fullSignificand, significantParts, bits);
 862     lost_fraction = combineLostFractions(lf, lost_fraction);
 863     exponent += bits;
 864   }
 865
 866   APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
 867
 868   if(newPartsCount > 4)
 869     delete [] fullSignificand;
 870
 871   return lost_fraction;
 872 }
 873
 874 /* Multiply the significands of LHS and RHS to DST.  */
 875 lostFraction
 876 APFloat::divideSignificand(const APFloat &rhs)
 877 {
 878   unsigned int bit, i, partsCount;
 879   const integerPart *rhsSignificand;
 880   integerPart *lhsSignificand, *dividend, *divisor;
 881   integerPart scratch[4];
 882   lostFraction lost_fraction;
 883
 884   assert(semantics == rhs.semantics);
 885
 886   lhsSignificand = significandParts();
 887   rhsSignificand = rhs.significandParts();
 888   partsCount = partCount();
 889
 890   if(partsCount > 2)
 891     dividend = new integerPart[partsCount * 2];
 892   else
 893     dividend = scratch;
 894
 895   divisor = dividend + partsCount;
 896
 897   /* Copy the dividend and divisor as they will be modified in-place.  */
 898   for(i = 0; i < partsCount; i++) {
 899     dividend[i] = lhsSignificand[i];
 900     divisor[i] = rhsSignificand[i];
 901     lhsSignificand[i] = 0;
 902   }
 903
 904   exponent -= rhs.exponent;
 905
 906   unsigned int precision = semantics->precision;
 907
 908   /* Normalize the divisor.  */
 909   bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
 910   if(bit) {
 911     exponent += bit;
 912     APInt::tcShiftLeft(divisor, partsCount, bit);
 913   }
 914
 915   /* Normalize the dividend.  */
 916   bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
 917   if(bit) {
 918     exponent -= bit;
 919     APInt::tcShiftLeft(dividend, partsCount, bit);
 920   }
 921
 922   /* Ensure the dividend >= divisor initially for the loop below.
 923      Incidentally, this means that the division loop below is
 924      guaranteed to set the integer bit to one.  */
 925   if(APInt::tcCompare(dividend, divisor, partsCount) < 0) {
 926     exponent--;
 927     APInt::tcShiftLeft(dividend, partsCount, 1);
 928     assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
 929   }
 930
 931   /* Long division.  */
 932   for(bit = precision; bit; bit -= 1) {
 933     if(APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
 934       APInt::tcSubtract(dividend, divisor, 0, partsCount);
 935       APInt::tcSetBit(lhsSignificand, bit - 1);
 936     }
 937
 938     APInt::tcShiftLeft(dividend, partsCount, 1);
 939   }
 940
 941   /* Figure out the lost fraction.  */
 942   int cmp = APInt::tcCompare(dividend, divisor, partsCount);
 943
 944   if(cmp > 0)
 945     lost_fraction = lfMoreThanHalf;
 946   else if(cmp == 0)
 947     lost_fraction = lfExactlyHalf;
 948   else if(APInt::tcIsZero(dividend, partsCount))
 949     lost_fraction = lfExactlyZero;
 950   else
 951     lost_fraction = lfLessThanHalf;
 952
 953   if(partsCount > 2)
 954     delete [] dividend;
 955
 956   return lost_fraction;
 957 }
 958
 959 unsigned int
 960 APFloat::significandMSB() const
 961 {
 962   return APInt::tcMSB(significandParts(), partCount());
 963 }
 964
 965 unsigned int
 966 APFloat::significandLSB() const
 967 {
 968   return APInt::tcLSB(significandParts(), partCount());
 969 }
 970
 971 /* Note that a zero result is NOT normalized to fcZero.  */
 972 lostFraction
 973 APFloat::shiftSignificandRight(unsigned int bits)
 974 {
 975   /* Our exponent should not overflow.  */
 976   assert((exponent_t) (exponent + bits) >= exponent);
 977
 978   exponent += bits;
 979
 980   return shiftRight(significandParts(), partCount(), bits);
 981 }
 982
 983 /* Shift the significand left BITS bits, subtract BITS from its exponent.  */
 984 void
 985 APFloat::shiftSignificandLeft(unsigned int bits)
 986 {
 987   assert(bits < semantics->precision);
 988
 989   if(bits) {
 990     unsigned int partsCount = partCount();
 991
 992     APInt::tcShiftLeft(significandParts(), partsCount, bits);
 993     exponent -= bits;
 994
 995     assert(!APInt::tcIsZero(significandParts(), partsCount));
 996   }
 997 }
 998
 999 APFloat::cmpResult
1000 APFloat::compareAbsoluteValue(const APFloat &rhs) const
1001 {
1002   int compare;
1003
1004   assert(semantics == rhs.semantics);
1005   assert(category == fcNormal);
1006   assert(rhs.category == fcNormal);
1007
1008   compare = exponent - rhs.exponent;
1009
1010   /* If exponents are equal, do an unsigned bignum comparison of the
1011      significands.  */
1012   if(compare == 0)
1013     compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
1014                                partCount());
1015
1016   if(compare > 0)
1017     return cmpGreaterThan;
1018   else if(compare < 0)
1019     return cmpLessThan;
1020   else
1021     return cmpEqual;
1022 }
1023
1024 /* Handle overflow.  Sign is preserved.  We either become infinity or
1025    the largest finite number.  */
1026 APFloat::opStatus
1027 APFloat::handleOverflow(roundingMode rounding_mode)
1028 {
1029   /* Infinity?  */
1030   if(rounding_mode == rmNearestTiesToEven
1031      || rounding_mode == rmNearestTiesToAway
1032      || (rounding_mode == rmTowardPositive && !sign)
1033      || (rounding_mode == rmTowardNegative && sign))
1034     {
1035       category = fcInfinity;
1036       return (opStatus) (opOverflow | opInexact);
1037     }
1038
1039   /* Otherwise we become the largest finite number.  */
1040   category = fcNormal;
1041   exponent = semantics->maxExponent;
1042   APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
1043                                    semantics->precision);
1044
1045   return opInexact;
1046 }
1047
1048 /* Returns TRUE if, when truncating the current number, with BIT the
1049    new LSB, with the given lost fraction and rounding mode, the result
1050    would need to be rounded away from zero (i.e., by increasing the
1051    signficand).  This routine must work for fcZero of both signs, and
1052    fcNormal numbers.  */
1053 bool
1054 APFloat::roundAwayFromZero(roundingMode rounding_mode,
1055                            lostFraction lost_fraction,
1056                            unsigned int bit) const
1057 {
1058   /* NaNs and infinities should not have lost fractions.  */
1059   assert(category == fcNormal || category == fcZero);
1060
1061   /* Current callers never pass this so we don't handle it.  */
1062   assert(lost_fraction != lfExactlyZero);
1063
1064   switch(rounding_mode) {
1065   default:
1066     assert(0);
1067
1068   case rmNearestTiesToAway:
1069     return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
1070
1071   case rmNearestTiesToEven:
1072     if(lost_fraction == lfMoreThanHalf)
1073       return true;
1074
1075     /* Our zeroes don't have a significand to test.  */
1076     if(lost_fraction == lfExactlyHalf && category != fcZero)
1077       return APInt::tcExtractBit(significandParts(), bit);
1078
1079     return false;
1080
1081   case rmTowardZero:
1082     return false;
1083
1084   case rmTowardPositive:
1085     return sign == false;
1086
1087   case rmTowardNegative:
1088     return sign == true;
1089   }
1090 }
1091
1092 APFloat::opStatus
1093 APFloat::normalize(roundingMode rounding_mode,
1094                    lostFraction lost_fraction)
1095 {
1096   unsigned int omsb;                /* One, not zero, based MSB.  */
1097   int exponentChange;
1098
1099   if(category != fcNormal)
1100     return opOK;
1101
1102   /* Before rounding normalize the exponent of fcNormal numbers.  */
1103   omsb = significandMSB() + 1;
1104
1105   if(omsb) {
1106     /* OMSB is numbered from 1.  We want to place it in the integer
1107        bit numbered PRECISON if possible, with a compensating change in
1108        the exponent.  */
1109     exponentChange = omsb - semantics->precision;
1110
1111     /* If the resulting exponent is too high, overflow according to
1112        the rounding mode.  */
1113     if(exponent + exponentChange > semantics->maxExponent)
1114       return handleOverflow(rounding_mode);
1115
1116     /* Subnormal numbers have exponent minExponent, and their MSB
1117        is forced based on that.  */
1118     if(exponent + exponentChange < semantics->minExponent)
1119       exponentChange = semantics->minExponent - exponent;
1120
1121     /* Shifting left is easy as we don't lose precision.  */
1122     if(exponentChange < 0) {
1123       assert(lost_fraction == lfExactlyZero);
1124
1125       shiftSignificandLeft(-exponentChange);
1126
1127       return opOK;
1128     }
1129
1130     if(exponentChange > 0) {
1131       lostFraction lf;
1132
1133       /* Shift right and capture any new lost fraction.  */
1134       lf = shiftSignificandRight(exponentChange);
1135
1136       lost_fraction = combineLostFractions(lf, lost_fraction);
1137
1138       /* Keep OMSB up-to-date.  */
1139       if(omsb > (unsigned) exponentChange)
1140         omsb -= exponentChange;
1141       else
1142         omsb = 0;
1143     }
1144   }
1145
1146   /* Now round the number according to rounding_mode given the lost
1147      fraction.  */
1148
1149   /* As specified in IEEE 754, since we do not trap we do not report
1150      underflow for exact results.  */
1151   if(lost_fraction == lfExactlyZero) {
1152     /* Canonicalize zeroes.  */
1153     if(omsb == 0)
1154       category = fcZero;
1155
1156     return opOK;
1157   }
1158
1159   /* Increment the significand if we're rounding away from zero.  */
1160   if(roundAwayFromZero(rounding_mode, lost_fraction, 0)) {
1161     if(omsb == 0)
1162       exponent = semantics->minExponent;
1163
1164     incrementSignificand();
1165     omsb = significandMSB() + 1;
1166
1167     /* Did the significand increment overflow?  */
1168     if(omsb == (unsigned) semantics->precision + 1) {
1169       /* Renormalize by incrementing the exponent and shifting our
1170          significand right one.  However if we already have the
1171          maximum exponent we overflow to infinity.  */
1172       if(exponent == semantics->maxExponent) {
1173         category = fcInfinity;
1174
1175         return (opStatus) (opOverflow | opInexact);
1176       }
1177
1178       shiftSignificandRight(1);
1179
1180       return opInexact;
1181     }
1182   }
1183
1184   /* The normal case - we were and are not denormal, and any
1185      significand increment above didn't overflow.  */
1186   if(omsb == semantics->precision)
1187     return opInexact;
1188
1189   /* We have a non-zero denormal.  */
1190   assert(omsb < semantics->precision);
1191
1192   /* Canonicalize zeroes.  */
1193   if(omsb == 0)
1194     category = fcZero;
1195
1196   /* The fcZero case is a denormal that underflowed to zero.  */
1197   return (opStatus) (opUnderflow | opInexact);
1198 }
1199
1200 APFloat::opStatus
1201 APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
1202 {
1203   switch(convolve(category, rhs.category)) {
1204   default:
1205     assert(0);
1206
1207   case convolve(fcNaN, fcZero):
1208   case convolve(fcNaN, fcNormal):
1209   case convolve(fcNaN, fcInfinity):
1210   case convolve(fcNaN, fcNaN):
1211   case convolve(fcNormal, fcZero):
1212   case convolve(fcInfinity, fcNormal):
1213   case convolve(fcInfinity, fcZero):
1214     return opOK;
1215
1216   case convolve(fcZero, fcNaN):
1217   case convolve(fcNormal, fcNaN):
1218   case convolve(fcInfinity, fcNaN):
1219     category = fcNaN;
1220     copySignificand(rhs);
1221     return opOK;
1222
1223   case convolve(fcNormal, fcInfinity):
1224   case convolve(fcZero, fcInfinity):
1225     category = fcInfinity;
1226     sign = rhs.sign ^ subtract;
1227     return opOK;
1228
1229   case convolve(fcZero, fcNormal):
1230     assign(rhs);
1231     sign = rhs.sign ^ subtract;
1232     return opOK;
1233
1234   case convolve(fcZero, fcZero):
1235     /* Sign depends on rounding mode; handled by caller.  */
1236     return opOK;
1237
1238   case convolve(fcInfinity, fcInfinity):
1239     /* Differently signed infinities can only be validly
1240        subtracted.  */
1241     if((sign ^ rhs.sign) != subtract) {
1242       makeNaN();
1243       return opInvalidOp;
1244     }
1245
1246     return opOK;
1247
1248   case convolve(fcNormal, fcNormal):
1249     return opDivByZero;
1250   }
1251 }
1252
1253 /* Add or subtract two normal numbers.  */
1254 lostFraction
1255 APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
1256 {
1257   integerPart carry;
1258   lostFraction lost_fraction;
1259   int bits;
1260
1261   /* Determine if the operation on the absolute values is effectively
1262      an addition or subtraction.  */
1263   subtract ^= (sign ^ rhs.sign) ? true : false;
1264
1265   /* Are we bigger exponent-wise than the RHS?  */
1266   bits = exponent - rhs.exponent;
1267
1268   /* Subtraction is more subtle than one might naively expect.  */
1269   if(subtract) {
1270     APFloat temp_rhs(rhs);
1271     bool reverse;
1272
1273     if (bits == 0) {
1274       reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
1275       lost_fraction = lfExactlyZero;
1276     } else if (bits > 0) {
1277       lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
1278       shiftSignificandLeft(1);
1279       reverse = false;
1280     } else {
1281       lost_fraction = shiftSignificandRight(-bits - 1);
1282       temp_rhs.shiftSignificandLeft(1);
1283       reverse = true;
1284     }
1285
1286     if (reverse) {
1287       carry = temp_rhs.subtractSignificand
1288         (*this, lost_fraction != lfExactlyZero);
1289       copySignificand(temp_rhs);
1290       sign = !sign;
1291     } else {
1292       carry = subtractSignificand
1293         (temp_rhs, lost_fraction != lfExactlyZero);
1294     }
1295
1296     /* Invert the lost fraction - it was on the RHS and
1297        subtracted.  */
1298     if(lost_fraction == lfLessThanHalf)
1299       lost_fraction = lfMoreThanHalf;
1300     else if(lost_fraction == lfMoreThanHalf)
1301       lost_fraction = lfLessThanHalf;
1302
1303     /* The code above is intended to ensure that no borrow is
1304        necessary.  */
1305     assert(!carry);
1306   } else {
1307     if(bits > 0) {
1308       APFloat temp_rhs(rhs);
1309
1310       lost_fraction = temp_rhs.shiftSignificandRight(bits);
1311       carry = addSignificand(temp_rhs);
1312     } else {
1313       lost_fraction = shiftSignificandRight(-bits);
1314       carry = addSignificand(rhs);
1315     }
1316
1317     /* We have a guard bit; generating a carry cannot happen.  */
1318     assert(!carry);
1319   }
1320
1321   return lost_fraction;
1322 }
1323
1324 APFloat::opStatus
1325 APFloat::multiplySpecials(const APFloat &rhs)
1326 {
1327   switch(convolve(category, rhs.category)) {
1328   default:
1329     assert(0);
1330
1331   case convolve(fcNaN, fcZero):
1332   case convolve(fcNaN, fcNormal):
1333   case convolve(fcNaN, fcInfinity):
1334   case convolve(fcNaN, fcNaN):
1335     return opOK;
1336
1337   case convolve(fcZero, fcNaN):
1338   case convolve(fcNormal, fcNaN):
1339   case convolve(fcInfinity, fcNaN):
1340     category = fcNaN;
1341     copySignificand(rhs);
1342     return opOK;
1343
1344   case convolve(fcNormal, fcInfinity):
1345   case convolve(fcInfinity, fcNormal):
1346   case convolve(fcInfinity, fcInfinity):
1347     category = fcInfinity;
1348     return opOK;
1349
1350   case convolve(fcZero, fcNormal):
1351   case convolve(fcNormal, fcZero):
1352   case convolve(fcZero, fcZero):
1353     category = fcZero;
1354     return opOK;
1355
1356   case convolve(fcZero, fcInfinity):
1357   case convolve(fcInfinity, fcZero):
1358     makeNaN();
1359     return opInvalidOp;
1360
1361   case convolve(fcNormal, fcNormal):
1362     return opOK;
1363   }
1364 }
1365
1366 APFloat::opStatus
1367 APFloat::divideSpecials(const APFloat &rhs)
1368 {
1369   switch(convolve(category, rhs.category)) {
1370   default:
1371     assert(0);
1372
1373   case convolve(fcNaN, fcZero):
1374   case convolve(fcNaN, fcNormal):
1375   case convolve(fcNaN, fcInfinity):
1376   case convolve(fcNaN, fcNaN):
1377   case convolve(fcInfinity, fcZero):
1378   case convolve(fcInfinity, fcNormal):
1379   case convolve(fcZero, fcInfinity):
1380   case convolve(fcZero, fcNormal):
1381     return opOK;
1382
1383   case convolve(fcZero, fcNaN):
1384   case convolve(fcNormal, fcNaN):
1385   case convolve(fcInfinity, fcNaN):
1386     category = fcNaN;
1387     copySignificand(rhs);
1388     return opOK;
1389
1390   case convolve(fcNormal, fcInfinity):
1391     category = fcZero;
1392     return opOK;
1393
1394   case convolve(fcNormal, fcZero):
1395     category = fcInfinity;
1396     return opDivByZero;
1397
1398   case convolve(fcInfinity, fcInfinity):
1399   case convolve(fcZero, fcZero):
1400     makeNaN();
1401     return opInvalidOp;
1402
1403   case convolve(fcNormal, fcNormal):
1404     return opOK;
1405   }
1406 }
1407
1408 /* Change sign.  */
1409 void
1410 APFloat::changeSign()
1411 {
1412   /* Look mummy, this one's easy.  */
1413   sign = !sign;
1414 }
1415
1416 void
1417 APFloat::clearSign()
1418 {
1419   /* So is this one. */
1420   sign = 0;
1421 }
1422
1423 void
1424 APFloat::copySign(const APFloat &rhs)
1425 {
1426   /* And this one. */
1427   sign = rhs.sign;
1428 }
1429
1430 /* Normalized addition or subtraction.  */
1431 APFloat::opStatus
1432 APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1433                        bool subtract)
1434 {
1435   opStatus fs;
1436
1437   assertArithmeticOK(*semantics);
1438
1439   fs = addOrSubtractSpecials(rhs, subtract);
1440
1441   /* This return code means it was not a simple case.  */
1442   if(fs == opDivByZero) {
1443     lostFraction lost_fraction;
1444
1445     lost_fraction = addOrSubtractSignificand(rhs, subtract);
1446     fs = normalize(rounding_mode, lost_fraction);
1447
1448     /* Can only be zero if we lost no fraction.  */
1449     assert(category != fcZero || lost_fraction == lfExactlyZero);
1450   }
1451
1452   /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1453      positive zero unless rounding to minus infinity, except that
1454      adding two like-signed zeroes gives that zero.  */
1455   if(category == fcZero) {
1456     if(rhs.category != fcZero || (sign == rhs.sign) == subtract)
1457       sign = (rounding_mode == rmTowardNegative);
1458   }
1459
1460   return fs;
1461 }
1462
1463 /* Normalized addition.  */
1464 APFloat::opStatus
1465 APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1466 {
1467   return addOrSubtract(rhs, rounding_mode, false);
1468 }
1469
1470 /* Normalized subtraction.  */
1471 APFloat::opStatus
1472 APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1473 {
1474   return addOrSubtract(rhs, rounding_mode, true);
1475 }
1476
1477 /* Normalized multiply.  */
1478 APFloat::opStatus
1479 APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1480 {
1481   opStatus fs;
1482
1483   assertArithmeticOK(*semantics);
1484   sign ^= rhs.sign;
1485   fs = multiplySpecials(rhs);
1486
1487   if(category == fcNormal) {
1488     lostFraction lost_fraction = multiplySignificand(rhs, 0);
1489     fs = normalize(rounding_mode, lost_fraction);
1490     if(lost_fraction != lfExactlyZero)
1491       fs = (opStatus) (fs | opInexact);
1492   }
1493
1494   return fs;
1495 }
1496
1497 /* Normalized divide.  */
1498 APFloat::opStatus
1499 APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1500 {
1501   opStatus fs;
1502
1503   assertArithmeticOK(*semantics);
1504   sign ^= rhs.sign;
1505   fs = divideSpecials(rhs);
1506
1507   if(category == fcNormal) {
1508     lostFraction lost_fraction = divideSignificand(rhs);
1509     fs = normalize(rounding_mode, lost_fraction);
1510     if(lost_fraction != lfExactlyZero)
1511       fs = (opStatus) (fs | opInexact);
1512   }
1513
1514   return fs;
1515 }
1516
1517 /* Normalized remainder.  This is not currently doing TRT.  */
1518 APFloat::opStatus
1519 APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1520 {
1521   opStatus fs;
1522   APFloat V = *this;
1523   unsigned int origSign = sign;
1524
1525   assertArithmeticOK(*semantics);
1526   fs = V.divide(rhs, rmNearestTiesToEven);
1527   if (fs == opDivByZero)
1528     return fs;
1529
1530   int parts = partCount();
1531   integerPart *x = new integerPart[parts];
1532   bool ignored;
1533   fs = V.convertToInteger(x, parts * integerPartWidth, true,
1534                           rmTowardZero, &ignored);
1535   if (fs==opInvalidOp)
1536     return fs;
1537
1538   fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1539                                         rmNearestTiesToEven);
1540   assert(fs==opOK);   // should always work
1541
1542   fs = V.multiply(rhs, rounding_mode);
1543   assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
1544
1545   fs = subtract(V, rounding_mode);
1546   assert(fs==opOK || fs==opInexact);   // likewise
1547
1548   if (isZero())
1549     sign = origSign;    // IEEE754 requires this
1550   delete[] x;
1551   return fs;
1552 }
1553
1554 /* Normalized fused-multiply-add.  */
1555 APFloat::opStatus
1556 APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1557                           const APFloat &addend,
1558                           roundingMode rounding_mode)
1559 {
1560   opStatus fs;
1561
1562   assertArithmeticOK(*semantics);
1563
1564   /* Post-multiplication sign, before addition.  */
1565   sign ^= multiplicand.sign;
1566
1567   /* If and only if all arguments are normal do we need to do an
1568      extended-precision calculation.  */
1569   if(category == fcNormal
1570      && multiplicand.category == fcNormal
1571      && addend.category == fcNormal) {
1572     lostFraction lost_fraction;
1573
1574     lost_fraction = multiplySignificand(multiplicand, &addend);
1575     fs = normalize(rounding_mode, lost_fraction);
1576     if(lost_fraction != lfExactlyZero)
1577       fs = (opStatus) (fs | opInexact);
1578
1579     /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1580        positive zero unless rounding to minus infinity, except that
1581        adding two like-signed zeroes gives that zero.  */
1582     if(category == fcZero && sign != addend.sign)
1583       sign = (rounding_mode == rmTowardNegative);
1584   } else {
1585     fs = multiplySpecials(multiplicand);
1586
1587     /* FS can only be opOK or opInvalidOp.  There is no more work
1588        to do in the latter case.  The IEEE-754R standard says it is
1589        implementation-defined in this case whether, if ADDEND is a
1590        quiet NaN, we raise invalid op; this implementation does so.
1591
1592        If we need to do the addition we can do so with normal
1593        precision.  */
1594     if(fs == opOK)
1595       fs = addOrSubtract(addend, rounding_mode, false);
1596   }
1597
1598   return fs;
1599 }
1600
1601 /* Comparison requires normalized numbers.  */
1602 APFloat::cmpResult
1603 APFloat::compare(const APFloat &rhs) const
1604 {
1605   cmpResult result;
1606
1607   assertArithmeticOK(*semantics);
1608   assert(semantics == rhs.semantics);
1609
1610   switch(convolve(category, rhs.category)) {
1611   default:
1612     assert(0);
1613
1614   case convolve(fcNaN, fcZero):
1615   case convolve(fcNaN, fcNormal):
1616   case convolve(fcNaN, fcInfinity):
1617   case convolve(fcNaN, fcNaN):
1618   case convolve(fcZero, fcNaN):
1619   case convolve(fcNormal, fcNaN):
1620   case convolve(fcInfinity, fcNaN):
1621     return cmpUnordered;
1622
1623   case convolve(fcInfinity, fcNormal):
1624   case convolve(fcInfinity, fcZero):
1625   case convolve(fcNormal, fcZero):
1626     if(sign)
1627       return cmpLessThan;
1628     else
1629       return cmpGreaterThan;
1630
1631   case convolve(fcNormal, fcInfinity):
1632   case convolve(fcZero, fcInfinity):
1633   case convolve(fcZero, fcNormal):
1634     if(rhs.sign)
1635       return cmpGreaterThan;
1636     else
1637       return cmpLessThan;
1638
1639   case convolve(fcInfinity, fcInfinity):
1640     if(sign == rhs.sign)
1641       return cmpEqual;
1642     else if(sign)
1643       return cmpLessThan;
1644     else
1645       return cmpGreaterThan;
1646
1647   case convolve(fcZero, fcZero):
1648     return cmpEqual;
1649
1650   case convolve(fcNormal, fcNormal):
1651     break;
1652   }
1653
1654   /* Two normal numbers.  Do they have the same sign?  */
1655   if(sign != rhs.sign) {
1656     if(sign)
1657       result = cmpLessThan;
1658     else
1659       result = cmpGreaterThan;
1660   } else {
1661     /* Compare absolute values; invert result if negative.  */
1662     result = compareAbsoluteValue(rhs);
1663
1664     if(sign) {
1665       if(result == cmpLessThan)
1666         result = cmpGreaterThan;
1667       else if(result == cmpGreaterThan)
1668         result = cmpLessThan;
1669     }
1670   }
1671
1672   return result;
1673 }
1674
1675 /// APFloat::convert - convert a value of one floating point type to another.
1676 /// The return value corresponds to the IEEE754 exceptions.  *losesInfo
1677 /// records whether the transformation lost information, i.e. whether
1678 /// converting the result back to the original type will produce the
1679 /// original value (this is almost the same as return value==fsOK, but there
1680 /// are edge cases where this is not so).
1681
1682 APFloat::opStatus
1683 APFloat::convert(const fltSemantics &toSemantics,
1684                  roundingMode rounding_mode, bool *losesInfo)
1685 {
1686   lostFraction lostFraction;
1687   unsigned int newPartCount, oldPartCount;
1688   opStatus fs;
1689
1690   assertArithmeticOK(*semantics);
1691   assertArithmeticOK(toSemantics);
1692   lostFraction = lfExactlyZero;
1693   newPartCount = partCountForBits(toSemantics.precision + 1);
1694   oldPartCount = partCount();
1695
1696   /* Handle storage complications.  If our new form is wider,
1697      re-allocate our bit pattern into wider storage.  If it is
1698      narrower, we ignore the excess parts, but if narrowing to a
1699      single part we need to free the old storage.
1700      Be careful not to reference significandParts for zeroes
1701      and infinities, since it aborts.  */
1702   if (newPartCount > oldPartCount) {
1703     integerPart *newParts;
1704     newParts = new integerPart[newPartCount];
1705     APInt::tcSet(newParts, 0, newPartCount);
1706     if (category==fcNormal || category==fcNaN)
1707       APInt::tcAssign(newParts, significandParts(), oldPartCount);
1708     freeSignificand();
1709     significand.parts = newParts;
1710   } else if (newPartCount < oldPartCount) {
1711     /* Capture any lost fraction through truncation of parts so we get
1712        correct rounding whilst normalizing.  */
1713     if (category==fcNormal)
1714       lostFraction = lostFractionThroughTruncation
1715         (significandParts(), oldPartCount, toSemantics.precision);
1716     if (newPartCount == 1) {
1717         integerPart newPart = 0;
1718         if (category==fcNormal || category==fcNaN)
1719           newPart = significandParts()[0];
1720         freeSignificand();
1721         significand.part = newPart;
1722     }
1723   }
1724
1725   if(category == fcNormal) {
1726     /* Re-interpret our bit-pattern.  */
1727     exponent += toSemantics.precision - semantics->precision;
1728     semantics = &toSemantics;
1729     fs = normalize(rounding_mode, lostFraction);
1730     *losesInfo = (fs != opOK);
1731   } else if (category == fcNaN) {
1732     int shift = toSemantics.precision - semantics->precision;
1733     // Do this now so significandParts gets the right answer
1734     const fltSemantics *oldSemantics = semantics;
1735     semantics = &toSemantics;
1736     *losesInfo = false;
1737     // No normalization here, just truncate
1738     if (shift>0)
1739       APInt::tcShiftLeft(significandParts(), newPartCount, shift);
1740     else if (shift < 0) {
1741       unsigned ushift = -shift;
1742       // Figure out if we are losing information.  This happens
1743       // if are shifting out something other than 0s, or if the x87 long
1744       // double input did not have its integer bit set (pseudo-NaN), or if the
1745       // x87 long double input did not have its QNan bit set (because the x87
1746       // hardware sets this bit when converting a lower-precision NaN to
1747       // x87 long double).
1748       if (APInt::tcLSB(significandParts(), newPartCount) < ushift)
1749         *losesInfo = true;
1750       if (oldSemantics == &APFloat::x87DoubleExtended &&
1751           (!(*significandParts() & 0x8000000000000000ULL) ||
1752            !(*significandParts() & 0x4000000000000000ULL)))
1753         *losesInfo = true;
1754       APInt::tcShiftRight(significandParts(), newPartCount, ushift);
1755     }
1756     // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
1757     // does not give you back the same bits.  This is dubious, and we
1758     // don't currently do it.  You're really supposed to get
1759     // an invalid operation signal at runtime, but nobody does that.
1760     fs = opOK;
1761   } else {
1762     semantics = &toSemantics;
1763     fs = opOK;
1764     *losesInfo = false;
1765   }
1766
1767   return fs;
1768 }
1769
1770 /* Convert a floating point number to an integer according to the
1771    rounding mode.  If the rounded integer value is out of range this
1772    returns an invalid operation exception and the contents of the
1773    destination parts are unspecified.  If the rounded value is in
1774    range but the floating point number is not the exact integer, the C
1775    standard doesn't require an inexact exception to be raised.  IEEE
1776    854 does require it so we do that.
1777
1778    Note that for conversions to integer type the C standard requires
1779    round-to-zero to always be used.  */
1780 APFloat::opStatus
1781 APFloat::convertToSignExtendedInteger(integerPart *parts, unsigned int width,
1782                                       bool isSigned,
1783                                       roundingMode rounding_mode,
1784                                       bool *isExact) const
1785 {
1786   lostFraction lost_fraction;
1787   const integerPart *src;
1788   unsigned int dstPartsCount, truncatedBits;
1789
1790   assertArithmeticOK(*semantics);
1791
1792   *isExact = false;
1793
1794   /* Handle the three special cases first.  */
1795   if(category == fcInfinity || category == fcNaN)
1796     return opInvalidOp;
1797
1798   dstPartsCount = partCountForBits(width);
1799
1800   if(category == fcZero) {
1801     APInt::tcSet(parts, 0, dstPartsCount);
1802     // Negative zero can't be represented as an int.
1803     *isExact = !sign;
1804     return opOK;
1805   }
1806
1807   src = significandParts();
1808
1809   /* Step 1: place our absolute value, with any fraction truncated, in
1810      the destination.  */
1811   if (exponent < 0) {
1812     /* Our absolute value is less than one; truncate everything.  */
1813     APInt::tcSet(parts, 0, dstPartsCount);
1814     /* For exponent -1 the integer bit represents .5, look at that.
1815        For smaller exponents leftmost truncated bit is 0. */
1816     truncatedBits = semantics->precision -1U - exponent;
1817   } else {
1818     /* We want the most significant (exponent + 1) bits; the rest are
1819        truncated.  */
1820     unsigned int bits = exponent + 1U;
1821
1822     /* Hopelessly large in magnitude?  */
1823     if (bits > width)
1824       return opInvalidOp;
1825
1826     if (bits < semantics->precision) {
1827       /* We truncate (semantics->precision - bits) bits.  */
1828       truncatedBits = semantics->precision - bits;
1829       APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits);
1830     } else {
1831       /* We want at least as many bits as are available.  */
1832       APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0);
1833       APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision);
1834       truncatedBits = 0;
1835     }
1836   }
1837
1838   /* Step 2: work out any lost fraction, and increment the absolute
1839      value if we would round away from zero.  */
1840   if (truncatedBits) {
1841     lost_fraction = lostFractionThroughTruncation(src, partCount(),
1842                                                   truncatedBits);
1843     if (lost_fraction != lfExactlyZero
1844         && roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) {
1845       if (APInt::tcIncrement(parts, dstPartsCount))
1846         return opInvalidOp;     /* Overflow.  */
1847     }
1848   } else {
1849     lost_fraction = lfExactlyZero;
1850   }
1851
1852   /* Step 3: check if we fit in the destination.  */
1853   unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1;
1854
1855   if (sign) {
1856     if (!isSigned) {
1857       /* Negative numbers cannot be represented as unsigned.  */
1858       if (omsb != 0)
1859         return opInvalidOp;
1860     } else {
1861       /* It takes omsb bits to represent the unsigned integer value.
1862          We lose a bit for the sign, but care is needed as the
1863          maximally negative integer is a special case.  */
1864       if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb)
1865         return opInvalidOp;
1866
1867       /* This case can happen because of rounding.  */
1868       if (omsb > width)
1869         return opInvalidOp;
1870     }
1871
1872     APInt::tcNegate (parts, dstPartsCount);
1873   } else {
1874     if (omsb >= width + !isSigned)
1875       return opInvalidOp;
1876   }
1877
1878   if (lost_fraction == lfExactlyZero) {
1879     *isExact = true;
1880     return opOK;
1881   } else
1882     return opInexact;
1883 }
1884
1885 /* Same as convertToSignExtendedInteger, except we provide
1886    deterministic values in case of an invalid operation exception,
1887    namely zero for NaNs and the minimal or maximal value respectively
1888    for underflow or overflow.
1889    The *isExact output tells whether the result is exact, in the sense
1890    that converting it back to the original floating point type produces
1891    the original value.  This is almost equivalent to result==opOK,
1892    except for negative zeroes.
1893 */
1894 APFloat::opStatus
1895 APFloat::convertToInteger(integerPart *parts, unsigned int width,
1896                           bool isSigned,
1897                           roundingMode rounding_mode, bool *isExact) const
1898 {
1899   opStatus fs;
1900
1901   fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode,
1902                                     isExact);
1903
1904   if (fs == opInvalidOp) {
1905     unsigned int bits, dstPartsCount;
1906
1907     dstPartsCount = partCountForBits(width);
1908
1909     if (category == fcNaN)
1910       bits = 0;
1911     else if (sign)
1912       bits = isSigned;
1913     else
1914       bits = width - isSigned;
1915
1916     APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits);
1917     if (sign && isSigned)
1918       APInt::tcShiftLeft(parts, dstPartsCount, width - 1);
1919   }
1920
1921   return fs;
1922 }
1923
1924 /* Convert an unsigned integer SRC to a floating point number,
1925    rounding according to ROUNDING_MODE.  The sign of the floating
1926    point number is not modified.  */
1927 APFloat::opStatus
1928 APFloat::convertFromUnsignedParts(const integerPart *src,
1929                                   unsigned int srcCount,
1930                                   roundingMode rounding_mode)
1931 {
1932   unsigned int omsb, precision, dstCount;
1933   integerPart *dst;
1934   lostFraction lost_fraction;
1935
1936   assertArithmeticOK(*semantics);
1937   category = fcNormal;
1938   omsb = APInt::tcMSB(src, srcCount) + 1;
1939   dst = significandParts();
1940   dstCount = partCount();
1941   precision = semantics->precision;
1942
1943   /* We want the most significant PRECISON bits of SRC.  There may not
1944      be that many; extract what we can.  */
1945   if (precision <= omsb) {
1946     exponent = omsb - 1;
1947     lost_fraction = lostFractionThroughTruncation(src, srcCount,
1948                                                   omsb - precision);
1949     APInt::tcExtract(dst, dstCount, src, precision, omsb - precision);
1950   } else {
1951     exponent = precision - 1;
1952     lost_fraction = lfExactlyZero;
1953     APInt::tcExtract(dst, dstCount, src, omsb, 0);
1954   }
1955
1956   return normalize(rounding_mode, lost_fraction);
1957 }
1958
1959 APFloat::opStatus
1960 APFloat::convertFromAPInt(const APInt &Val,
1961                           bool isSigned,
1962                           roundingMode rounding_mode)
1963 {
1964   unsigned int partCount = Val.getNumWords();
1965   APInt api = Val;
1966
1967   sign = false;
1968   if (isSigned && api.isNegative()) {
1969     sign = true;
1970     api = -api;
1971   }
1972
1973   return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
1974 }
1975
1976 /* Convert a two's complement integer SRC to a floating point number,
1977    rounding according to ROUNDING_MODE.  ISSIGNED is true if the
1978    integer is signed, in which case it must be sign-extended.  */
1979 APFloat::opStatus
1980 APFloat::convertFromSignExtendedInteger(const integerPart *src,
1981                                         unsigned int srcCount,
1982                                         bool isSigned,
1983                                         roundingMode rounding_mode)
1984 {
1985   opStatus status;
1986
1987   assertArithmeticOK(*semantics);
1988   if (isSigned
1989       && APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) {
1990     integerPart *copy;
1991
1992     /* If we're signed and negative negate a copy.  */
1993     sign = true;
1994     copy = new integerPart[srcCount];
1995     APInt::tcAssign(copy, src, srcCount);
1996     APInt::tcNegate(copy, srcCount);
1997     status = convertFromUnsignedParts(copy, srcCount, rounding_mode);
1998     delete [] copy;
1999   } else {
2000     sign = false;
2001     status = convertFromUnsignedParts(src, srcCount, rounding_mode);
2002   }
2003
2004   return status;
2005 }
2006
2007 /* FIXME: should this just take a const APInt reference?  */
2008 APFloat::opStatus
2009 APFloat::convertFromZeroExtendedInteger(const integerPart *parts,
2010                                         unsigned int width, bool isSigned,
2011                                         roundingMode rounding_mode)
2012 {
2013   unsigned int partCount = partCountForBits(width);
2014   APInt api = APInt(width, partCount, parts);
2015
2016   sign = false;
2017   if(isSigned && APInt::tcExtractBit(parts, width - 1)) {
2018     sign = true;
2019     api = -api;
2020   }
2021
2022   return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2023 }
2024
2025 APFloat::opStatus
2026 APFloat::convertFromHexadecimalString(const char *p,
2027                                       roundingMode rounding_mode)
2028 {
2029   lostFraction lost_fraction;
2030   integerPart *significand;
2031   unsigned int bitPos, partsCount;
2032   const char *dot, *firstSignificantDigit;
2033
2034   zeroSignificand();
2035   exponent = 0;
2036   category = fcNormal;
2037
2038   significand = significandParts();
2039   partsCount = partCount();
2040   bitPos = partsCount * integerPartWidth;
2041
2042   /* Skip leading zeroes and any (hexa)decimal point.  */
2043   p = skipLeadingZeroesAndAnyDot(p, &dot);
2044   firstSignificantDigit = p;
2045
2046   for(;;) {
2047     integerPart hex_value;
2048
2049     if(*p == '.') {
2050       assert(dot == 0);
2051       dot = p++;
2052     }
2053
2054     hex_value = hexDigitValue(*p);
2055     if(hex_value == -1U) {
2056       lost_fraction = lfExactlyZero;
2057       break;
2058     }
2059
2060     p++;
2061
2062     /* Store the number whilst 4-bit nibbles remain.  */
2063     if(bitPos) {
2064       bitPos -= 4;
2065       hex_value <<= bitPos % integerPartWidth;
2066       significand[bitPos / integerPartWidth] |= hex_value;
2067     } else {
2068       lost_fraction = trailingHexadecimalFraction(p, hex_value);
2069       while(hexDigitValue(*p) != -1U)
2070         p++;
2071       break;
2072     }
2073   }
2074
2075   /* Hex floats require an exponent but not a hexadecimal point.  */
2076   assert(*p == 'p' || *p == 'P');
2077
2078   /* Ignore the exponent if we are zero.  */
2079   if(p != firstSignificantDigit) {
2080     int expAdjustment;
2081
2082     /* Implicit hexadecimal point?  */
2083     if(!dot)
2084       dot = p;
2085
2086     /* Calculate the exponent adjustment implicit in the number of
2087        significant digits.  */
2088     expAdjustment = static_cast<int>(dot - firstSignificantDigit);
2089     if(expAdjustment < 0)
2090       expAdjustment++;
2091     expAdjustment = expAdjustment * 4 - 1;
2092
2093     /* Adjust for writing the significand starting at the most
2094        significant nibble.  */
2095     expAdjustment += semantics->precision;
2096     expAdjustment -= partsCount * integerPartWidth;
2097
2098     /* Adjust for the given exponent.  */
2099     exponent = totalExponent(p, expAdjustment);
2100   }
2101
2102   return normalize(rounding_mode, lost_fraction);
2103 }
2104
2105 APFloat::opStatus
2106 APFloat::roundSignificandWithExponent(const integerPart *decSigParts,
2107                                       unsigned sigPartCount, int exp,
2108                                       roundingMode rounding_mode)
2109 {
2110   unsigned int parts, pow5PartCount;
2111   fltSemantics calcSemantics = { 32767, -32767, 0, true };
2112   integerPart pow5Parts[maxPowerOfFiveParts];
2113   bool isNearest;
2114
2115   isNearest = (rounding_mode == rmNearestTiesToEven
2116                || rounding_mode == rmNearestTiesToAway);
2117
2118   parts = partCountForBits(semantics->precision + 11);
2119
2120   /* Calculate pow(5, abs(exp)).  */
2121   pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp);
2122
2123   for (;; parts *= 2) {
2124     opStatus sigStatus, powStatus;
2125     unsigned int excessPrecision, truncatedBits;
2126
2127     calcSemantics.precision = parts * integerPartWidth - 1;
2128     excessPrecision = calcSemantics.precision - semantics->precision;
2129     truncatedBits = excessPrecision;
2130
2131     APFloat decSig(calcSemantics, fcZero, sign);
2132     APFloat pow5(calcSemantics, fcZero, false);
2133
2134     sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount,
2135                                                 rmNearestTiesToEven);
2136     powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount,
2137                                               rmNearestTiesToEven);
2138     /* Add exp, as 10^n = 5^n * 2^n.  */
2139     decSig.exponent += exp;
2140
2141     lostFraction calcLostFraction;
2142     integerPart HUerr, HUdistance;
2143     unsigned int powHUerr;
2144
2145     if (exp >= 0) {
2146       /* multiplySignificand leaves the precision-th bit set to 1.  */
2147       calcLostFraction = decSig.multiplySignificand(pow5, NULL);
2148       powHUerr = powStatus != opOK;
2149     } else {
2150       calcLostFraction = decSig.divideSignificand(pow5);
2151       /* Denormal numbers have less precision.  */
2152       if (decSig.exponent < semantics->minExponent) {
2153         excessPrecision += (semantics->minExponent - decSig.exponent);
2154         truncatedBits = excessPrecision;
2155         if (excessPrecision > calcSemantics.precision)
2156           excessPrecision = calcSemantics.precision;
2157       }
2158       /* Extra half-ulp lost in reciprocal of exponent.  */
2159       powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2;
2160     }
2161
2162     /* Both multiplySignificand and divideSignificand return the
2163        result with the integer bit set.  */
2164     assert (APInt::tcExtractBit
2165             (decSig.significandParts(), calcSemantics.precision - 1) == 1);
2166
2167     HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK,
2168                        powHUerr);
2169     HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(),
2170                                       excessPrecision, isNearest);
2171
2172     /* Are we guaranteed to round correctly if we truncate?  */
2173     if (HUdistance >= HUerr) {
2174       APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(),
2175                        calcSemantics.precision - excessPrecision,
2176                        excessPrecision);
2177       /* Take the exponent of decSig.  If we tcExtract-ed less bits
2178          above we must adjust our exponent to compensate for the
2179          implicit right shift.  */
2180       exponent = (decSig.exponent + semantics->precision
2181                   - (calcSemantics.precision - excessPrecision));
2182       calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(),
2183                                                        decSig.partCount(),
2184                                                        truncatedBits);
2185       return normalize(rounding_mode, calcLostFraction);
2186     }
2187   }
2188 }
2189
2190 APFloat::opStatus
2191 APFloat::convertFromDecimalString(const char *p, roundingMode rounding_mode)
2192 {
2193   decimalInfo D;
2194   opStatus fs;
2195
2196   /* Scan the text.  */
2197   interpretDecimal(p, &D);
2198
2199   /* Handle the quick cases.  First the case of no significant digits,
2200      i.e. zero, and then exponents that are obviously too large or too
2201      small.  Writing L for log 10 / log 2, a number d.ddddd*10^exp
2202      definitely overflows if
2203
2204            (exp - 1) * L >= maxExponent
2205
2206      and definitely underflows to zero where
2207
2208            (exp + 1) * L <= minExponent - precision
2209
2210      With integer arithmetic the tightest bounds for L are
2211
2212            93/28 < L < 196/59            [ numerator <= 256 ]
2213            42039/12655 < L < 28738/8651  [ numerator <= 65536 ]
2214   */
2215
2216   if (decDigitValue(*D.firstSigDigit) >= 10U) {
2217     category = fcZero;
2218     fs = opOK;
2219   } else if ((D.normalizedExponent + 1) * 28738
2220              <= 8651 * (semantics->minExponent - (int) semantics->precision)) {
2221     /* Underflow to zero and round.  */
2222     zeroSignificand();
2223     fs = normalize(rounding_mode, lfLessThanHalf);
2224   } else if ((D.normalizedExponent - 1) * 42039
2225              >= 12655 * semantics->maxExponent) {
2226     /* Overflow and round.  */
2227     fs = handleOverflow(rounding_mode);
2228   } else {
2229     integerPart *decSignificand;
2230     unsigned int partCount;
2231
2232     /* A tight upper bound on number of bits required to hold an
2233        N-digit decimal integer is N * 196 / 59.  Allocate enough space
2234        to hold the full significand, and an extra part required by
2235        tcMultiplyPart.  */
2236     partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1;
2237     partCount = partCountForBits(1 + 196 * partCount / 59);
2238     decSignificand = new integerPart[partCount + 1];
2239     partCount = 0;
2240
2241     /* Convert to binary efficiently - we do almost all multiplication
2242        in an integerPart.  When this would overflow do we do a single
2243        bignum multiplication, and then revert again to multiplication
2244        in an integerPart.  */
2245     do {
2246       integerPart decValue, val, multiplier;
2247
2248       val = 0;
2249       multiplier = 1;
2250
2251       do {
2252         if (*p == '.')
2253           p++;
2254
2255         decValue = decDigitValue(*p++);
2256         multiplier *= 10;
2257         val = val * 10 + decValue;
2258         /* The maximum number that can be multiplied by ten with any
2259            digit added without overflowing an integerPart.  */
2260       } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10);
2261
2262       /* Multiply out the current part.  */
2263       APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val,
2264                             partCount, partCount + 1, false);
2265
2266       /* If we used another part (likely but not guaranteed), increase
2267          the count.  */
2268       if (decSignificand[partCount])
2269         partCount++;
2270     } while (p <= D.lastSigDigit);
2271
2272     category = fcNormal;
2273     fs = roundSignificandWithExponent(decSignificand, partCount,
2274                                       D.exponent, rounding_mode);
2275
2276     delete [] decSignificand;
2277   }
2278
2279   return fs;
2280 }
2281
2282 APFloat::opStatus
2283 APFloat::convertFromString(const char *p, roundingMode rounding_mode)
2284 {
2285   assertArithmeticOK(*semantics);
2286
2287   /* Handle a leading minus sign.  */
2288   if(*p == '-')
2289     sign = 1, p++;
2290   else
2291     sign = 0;
2292
2293   if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
2294     return convertFromHexadecimalString(p + 2, rounding_mode);
2295
2296   return convertFromDecimalString(p, rounding_mode);
2297 }
2298
2299 /* Write out a hexadecimal representation of the floating point value
2300    to DST, which must be of sufficient size, in the C99 form
2301    [-]0xh.hhhhp[+-]d.  Return the number of characters written,
2302    excluding the terminating NUL.
2303
2304    If UPPERCASE, the output is in upper case, otherwise in lower case.
2305
2306    HEXDIGITS digits appear altogether, rounding the value if
2307    necessary.  If HEXDIGITS is 0, the minimal precision to display the
2308    number precisely is used instead.  If nothing would appear after
2309    the decimal point it is suppressed.
2310
2311    The decimal exponent is always printed and has at least one digit.
2312    Zero values display an exponent of zero.  Infinities and NaNs
2313    appear as "infinity" or "nan" respectively.
2314
2315    The above rules are as specified by C99.  There is ambiguity about
2316    what the leading hexadecimal digit should be.  This implementation
2317    uses whatever is necessary so that the exponent is displayed as
2318    stored.  This implies the exponent will fall within the IEEE format
2319    range, and the leading hexadecimal digit will be 0 (for denormals),
2320    1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with
2321    any other digits zero).
2322 */
2323 unsigned int
2324 APFloat::convertToHexString(char *dst, unsigned int hexDigits,
2325                             bool upperCase, roundingMode rounding_mode) const
2326 {
2327   char *p;
2328
2329   assertArithmeticOK(*semantics);
2330
2331   p = dst;
2332   if (sign)
2333     *dst++ = '-';
2334
2335   switch (category) {
2336   case fcInfinity:
2337     memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1);
2338     dst += sizeof infinityL - 1;
2339     break;
2340
2341   case fcNaN:
2342     memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1);
2343     dst += sizeof NaNU - 1;
2344     break;
2345
2346   case fcZero:
2347     *dst++ = '0';
2348     *dst++ = upperCase ? 'X': 'x';
2349     *dst++ = '0';
2350     if (hexDigits > 1) {
2351       *dst++ = '.';
2352       memset (dst, '0', hexDigits - 1);
2353       dst += hexDigits - 1;
2354     }
2355     *dst++ = upperCase ? 'P': 'p';
2356     *dst++ = '0';
2357     break;
2358
2359   case fcNormal:
2360     dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode);
2361     break;
2362   }
2363
2364   *dst = 0;
2365
2366   return static_cast<unsigned int>(dst - p);
2367 }
2368
2369 /* Does the hard work of outputting the correctly rounded hexadecimal
2370    form of a normal floating point number with the specified number of
2371    hexadecimal digits.  If HEXDIGITS is zero the minimum number of
2372    digits necessary to print the value precisely is output.  */
2373 char *
2374 APFloat::convertNormalToHexString(char *dst, unsigned int hexDigits,
2375                                   bool upperCase,
2376                                   roundingMode rounding_mode) const
2377 {
2378   unsigned int count, valueBits, shift, partsCount, outputDigits;
2379   const char *hexDigitChars;
2380   const integerPart *significand;
2381   char *p;
2382   bool roundUp;
2383
2384   *dst++ = '0';
2385   *dst++ = upperCase ? 'X': 'x';
2386
2387   roundUp = false;
2388   hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower;
2389
2390   significand = significandParts();
2391   partsCount = partCount();
2392
2393   /* +3 because the first digit only uses the single integer bit, so
2394      we have 3 virtual zero most-significant-bits.  */
2395   valueBits = semantics->precision + 3;
2396   shift = integerPartWidth - valueBits % integerPartWidth;
2397
2398   /* The natural number of digits required ignoring trailing
2399      insignificant zeroes.  */
2400   outputDigits = (valueBits - significandLSB () + 3) / 4;
2401
2402   /* hexDigits of zero means use the required number for the
2403      precision.  Otherwise, see if we are truncating.  If we are,
2404      find out if we need to round away from zero.  */
2405   if (hexDigits) {
2406     if (hexDigits < outputDigits) {
2407       /* We are dropping non-zero bits, so need to check how to round.
2408          "bits" is the number of dropped bits.  */
2409       unsigned int bits;
2410       lostFraction fraction;
2411
2412       bits = valueBits - hexDigits * 4;
2413       fraction = lostFractionThroughTruncation (significand, partsCount, bits);
2414       roundUp = roundAwayFromZero(rounding_mode, fraction, bits);
2415     }
2416     outputDigits = hexDigits;
2417   }
2418
2419   /* Write the digits consecutively, and start writing in the location
2420      of the hexadecimal point.  We move the most significant digit
2421      left and add the hexadecimal point later.  */
2422   p = ++dst;
2423
2424   count = (valueBits + integerPartWidth - 1) / integerPartWidth;
2425
2426   while (outputDigits && count) {
2427     integerPart part;
2428
2429     /* Put the most significant integerPartWidth bits in "part".  */
2430     if (--count == partsCount)
2431       part = 0;  /* An imaginary higher zero part.  */
2432     else
2433       part = significand[count] << shift;
2434
2435     if (count && shift)
2436       part |= significand[count - 1] >> (integerPartWidth - shift);
2437
2438     /* Convert as much of "part" to hexdigits as we can.  */
2439     unsigned int curDigits = integerPartWidth / 4;
2440
2441     if (curDigits > outputDigits)
2442       curDigits = outputDigits;
2443     dst += partAsHex (dst, part, curDigits, hexDigitChars);
2444     outputDigits -= curDigits;
2445   }
2446
2447   if (roundUp) {
2448     char *q = dst;
2449
2450     /* Note that hexDigitChars has a trailing '0'.  */
2451     do {
2452       q--;
2453       *q = hexDigitChars[hexDigitValue (*q) + 1];
2454     } while (*q == '0');
2455     assert (q >= p);
2456   } else {
2457     /* Add trailing zeroes.  */
2458     memset (dst, '0', outputDigits);
2459     dst += outputDigits;
2460   }
2461
2462   /* Move the most significant digit to before the point, and if there
2463      is something after the decimal point add it.  This must come
2464      after rounding above.  */
2465   p[-1] = p[0];
2466   if (dst -1 == p)
2467     dst--;
2468   else
2469     p[0] = '.';
2470
2471   /* Finally output the exponent.  */
2472   *dst++ = upperCase ? 'P': 'p';
2473
2474   return writeSignedDecimal (dst, exponent);
2475 }
2476
2477 // For good performance it is desirable for different APFloats
2478 // to produce different integers.
2479 uint32_t
2480 APFloat::getHashValue() const
2481 {
2482   if (category==fcZero) return sign<<8 | semantics->precision ;
2483   else if (category==fcInfinity) return sign<<9 | semantics->precision;
2484   else if (category==fcNaN) return 1<<10 | semantics->precision;
2485   else {
2486     uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
2487     const integerPart* p = significandParts();
2488     for (int i=partCount(); i>0; i--, p++)
2489       hash ^= ((uint32_t)*p) ^ (uint32_t)((*p)>>32);
2490     return hash;
2491   }
2492 }
2493
2494 // Conversion from APFloat to/from host float/double.  It may eventually be
2495 // possible to eliminate these and have everybody deal with APFloats, but that
2496 // will take a while.  This approach will not easily extend to long double.
2497 // Current implementation requires integerPartWidth==64, which is correct at
2498 // the moment but could be made more general.
2499
2500 // Denormals have exponent minExponent in APFloat, but minExponent-1 in
2501 // the actual IEEE respresentations.  We compensate for that here.
2502
2503 APInt
2504 APFloat::convertF80LongDoubleAPFloatToAPInt() const
2505 {
2506   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended);
2507   assert (partCount()==2);
2508
2509   uint64_t myexponent, mysignificand;
2510
2511   if (category==fcNormal) {
2512     myexponent = exponent+16383; //bias
2513     mysignificand = significandParts()[0];
2514     if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
2515       myexponent = 0;   // denormal
2516   } else if (category==fcZero) {
2517     myexponent = 0;
2518     mysignificand = 0;
2519   } else if (category==fcInfinity) {
2520     myexponent = 0x7fff;
2521     mysignificand = 0x8000000000000000ULL;
2522   } else {
2523     assert(category == fcNaN && "Unknown category");
2524     myexponent = 0x7fff;
2525     mysignificand = significandParts()[0];
2526   }
2527
2528   uint64_t words[2];
2529   words[0] =  ((uint64_t)(sign & 1) << 63) |
2530               ((myexponent & 0x7fffLL) <<  48) |
2531               ((mysignificand >>16) & 0xffffffffffffLL);
2532   words[1] = mysignificand & 0xffff;
2533   return APInt(80, 2, words);
2534 }
2535
2536 APInt
2537 APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const
2538 {
2539   assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble);
2540   assert (partCount()==2);
2541
2542   uint64_t myexponent, mysignificand, myexponent2, mysignificand2;
2543
2544   if (category==fcNormal) {
2545     myexponent = exponent + 1023; //bias
2546     myexponent2 = exponent2 + 1023;
2547     mysignificand = significandParts()[0];
2548     mysignificand2 = significandParts()[1];
2549     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2550       myexponent = 0;   // denormal
2551     if (myexponent2==1 && !(mysignificand2 & 0x10000000000000LL))
2552       myexponent2 = 0;   // denormal
2553   } else if (category==fcZero) {
2554     myexponent = 0;
2555     mysignificand = 0;
2556     myexponent2 = 0;
2557     mysignificand2 = 0;
2558   } else if (category==fcInfinity) {
2559     myexponent = 0x7ff;
2560     myexponent2 = 0;
2561     mysignificand = 0;
2562     mysignificand2 = 0;
2563   } else {
2564     assert(category == fcNaN && "Unknown category");
2565     myexponent = 0x7ff;
2566     mysignificand = significandParts()[0];
2567     myexponent2 = exponent2;
2568     mysignificand2 = significandParts()[1];
2569   }
2570
2571   uint64_t words[2];
2572   words[0] =  ((uint64_t)(sign & 1) << 63) |
2573               ((myexponent & 0x7ff) <<  52) |
2574               (mysignificand & 0xfffffffffffffLL);
2575   words[1] =  ((uint64_t)(sign2 & 1) << 63) |
2576               ((myexponent2 & 0x7ff) <<  52) |
2577               (mysignificand2 & 0xfffffffffffffLL);
2578   return APInt(128, 2, words);
2579 }
2580
2581 APInt
2582 APFloat::convertDoubleAPFloatToAPInt() const
2583 {
2584   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2585   assert (partCount()==1);
2586
2587   uint64_t myexponent, mysignificand;
2588
2589   if (category==fcNormal) {
2590     myexponent = exponent+1023; //bias
2591     mysignificand = *significandParts();
2592     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2593       myexponent = 0;   // denormal
2594   } else if (category==fcZero) {
2595     myexponent = 0;
2596     mysignificand = 0;
2597   } else if (category==fcInfinity) {
2598     myexponent = 0x7ff;
2599     mysignificand = 0;
2600   } else {
2601     assert(category == fcNaN && "Unknown category!");
2602     myexponent = 0x7ff;
2603     mysignificand = *significandParts();
2604   }
2605
2606   return APInt(64, ((((uint64_t)(sign & 1) << 63) |
2607                      ((myexponent & 0x7ff) <<  52) |
2608                      (mysignificand & 0xfffffffffffffLL))));
2609 }
2610
2611 APInt
2612 APFloat::convertFloatAPFloatToAPInt() const
2613 {
2614   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2615   assert (partCount()==1);
2616
2617   uint32_t myexponent, mysignificand;
2618
2619   if (category==fcNormal) {
2620     myexponent = exponent+127; //bias
2621     mysignificand = (uint32_t)*significandParts();
2622     if (myexponent == 1 && !(mysignificand & 0x800000))
2623       myexponent = 0;   // denormal
2624   } else if (category==fcZero) {
2625     myexponent = 0;
2626     mysignificand = 0;
2627   } else if (category==fcInfinity) {
2628     myexponent = 0xff;
2629     mysignificand = 0;
2630   } else {
2631     assert(category == fcNaN && "Unknown category!");
2632     myexponent = 0xff;
2633     mysignificand = (uint32_t)*significandParts();
2634   }
2635
2636   return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
2637                     (mysignificand & 0x7fffff)));
2638 }
2639
2640 // This function creates an APInt that is just a bit map of the floating
2641 // point constant as it would appear in memory.  It is not a conversion,
2642 // and treating the result as a normal integer is unlikely to be useful.
2643
2644 APInt
2645 APFloat::bitcastToAPInt() const
2646 {
2647   if (semantics == (const llvm::fltSemantics*)&IEEEsingle)
2648     return convertFloatAPFloatToAPInt();
2649
2650   if (semantics == (const llvm::fltSemantics*)&IEEEdouble)
2651     return convertDoubleAPFloatToAPInt();
2652
2653   if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble)
2654     return convertPPCDoubleDoubleAPFloatToAPInt();
2655
2656   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended &&
2657          "unknown format!");
2658   return convertF80LongDoubleAPFloatToAPInt();
2659 }
2660
2661 float
2662 APFloat::convertToFloat() const
2663 {
2664   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2665   APInt api = bitcastToAPInt();
2666   return api.bitsToFloat();
2667 }
2668
2669 double
2670 APFloat::convertToDouble() const
2671 {
2672   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2673   APInt api = bitcastToAPInt();
2674   return api.bitsToDouble();
2675 }
2676
2677 /// Integer bit is explicit in this format.  Intel hardware (387 and later)
2678 /// does not support these bit patterns:
2679 ///  exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
2680 ///  exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
2681 ///  exponent = 0, integer bit 1 ("pseudodenormal")
2682 ///  exponent!=0 nor all 1's, integer bit 0 ("unnormal")
2683 /// At the moment, the first two are treated as NaNs, the second two as Normal.
2684 void
2685 APFloat::initFromF80LongDoubleAPInt(const APInt &api)
2686 {
2687   assert(api.getBitWidth()==80);
2688   uint64_t i1 = api.getRawData()[0];
2689   uint64_t i2 = api.getRawData()[1];
2690   uint64_t myexponent = (i1 >> 48) & 0x7fff;
2691   uint64_t mysignificand = ((i1 << 16) &  0xffffffffffff0000ULL) |
2692                           (i2 & 0xffff);
2693
2694   initialize(&APFloat::x87DoubleExtended);
2695   assert(partCount()==2);
2696
2697   sign = static_cast<unsigned int>(i1>>63);
2698   if (myexponent==0 && mysignificand==0) {
2699     // exponent, significand meaningless
2700     category = fcZero;
2701   } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
2702     // exponent, significand meaningless
2703     category = fcInfinity;
2704   } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
2705     // exponent meaningless
2706     category = fcNaN;
2707     significandParts()[0] = mysignificand;
2708     significandParts()[1] = 0;
2709   } else {
2710     category = fcNormal;
2711     exponent = myexponent - 16383;
2712     significandParts()[0] = mysignificand;
2713     significandParts()[1] = 0;
2714     if (myexponent==0)          // denormal
2715       exponent = -16382;
2716   }
2717 }
2718
2719 void
2720 APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api)
2721 {
2722   assert(api.getBitWidth()==128);
2723   uint64_t i1 = api.getRawData()[0];
2724   uint64_t i2 = api.getRawData()[1];
2725   uint64_t myexponent = (i1 >> 52) & 0x7ff;
2726   uint64_t mysignificand = i1 & 0xfffffffffffffLL;
2727   uint64_t myexponent2 = (i2 >> 52) & 0x7ff;
2728   uint64_t mysignificand2 = i2 & 0xfffffffffffffLL;
2729
2730   initialize(&APFloat::PPCDoubleDouble);
2731   assert(partCount()==2);
2732
2733   sign = static_cast<unsigned int>(i1>>63);
2734   sign2 = static_cast<unsigned int>(i2>>63);
2735   if (myexponent==0 && mysignificand==0) {
2736     // exponent, significand meaningless
2737     // exponent2 and significand2 are required to be 0; we don't check
2738     category = fcZero;
2739   } else if (myexponent==0x7ff && mysignificand==0) {
2740     // exponent, significand meaningless
2741     // exponent2 and significand2 are required to be 0; we don't check
2742     category = fcInfinity;
2743   } else if (myexponent==0x7ff && mysignificand!=0) {
2744     // exponent meaningless.  So is the whole second word, but keep it
2745     // for determinism.
2746     category = fcNaN;
2747     exponent2 = myexponent2;
2748     significandParts()[0] = mysignificand;
2749     significandParts()[1] = mysignificand2;
2750   } else {
2751     category = fcNormal;
2752     // Note there is no category2; the second word is treated as if it is
2753     // fcNormal, although it might be something else considered by itself.
2754     exponent = myexponent - 1023;
2755     exponent2 = myexponent2 - 1023;
2756     significandParts()[0] = mysignificand;
2757     significandParts()[1] = mysignificand2;
2758     if (myexponent==0)          // denormal
2759       exponent = -1022;
2760     else
2761       significandParts()[0] |= 0x10000000000000LL;  // integer bit
2762     if (myexponent2==0)
2763       exponent2 = -1022;
2764     else
2765       significandParts()[1] |= 0x10000000000000LL;  // integer bit
2766   }
2767 }
2768
2769 void
2770 APFloat::initFromDoubleAPInt(const APInt &api)
2771 {
2772   assert(api.getBitWidth()==64);
2773   uint64_t i = *api.getRawData();
2774   uint64_t myexponent = (i >> 52) & 0x7ff;
2775   uint64_t mysignificand = i & 0xfffffffffffffLL;
2776
2777   initialize(&APFloat::IEEEdouble);
2778   assert(partCount()==1);
2779
2780   sign = static_cast<unsigned int>(i>>63);
2781   if (myexponent==0 && mysignificand==0) {
2782     // exponent, significand meaningless
2783     category = fcZero;
2784   } else if (myexponent==0x7ff && mysignificand==0) {
2785     // exponent, significand meaningless
2786     category = fcInfinity;
2787   } else if (myexponent==0x7ff && mysignificand!=0) {
2788     // exponent meaningless
2789     category = fcNaN;
2790     *significandParts() = mysignificand;
2791   } else {
2792     category = fcNormal;
2793     exponent = myexponent - 1023;
2794     *significandParts() = mysignificand;
2795     if (myexponent==0)          // denormal
2796       exponent = -1022;
2797     else
2798       *significandParts() |= 0x10000000000000LL;  // integer bit
2799   }
2800 }
2801
2802 void
2803 APFloat::initFromFloatAPInt(const APInt & api)
2804 {
2805   assert(api.getBitWidth()==32);
2806   uint32_t i = (uint32_t)*api.getRawData();
2807   uint32_t myexponent = (i >> 23) & 0xff;
2808   uint32_t mysignificand = i & 0x7fffff;
2809
2810   initialize(&APFloat::IEEEsingle);
2811   assert(partCount()==1);
2812
2813   sign = i >> 31;
2814   if (myexponent==0 && mysignificand==0) {
2815     // exponent, significand meaningless
2816     category = fcZero;
2817   } else if (myexponent==0xff && mysignificand==0) {
2818     // exponent, significand meaningless
2819     category = fcInfinity;
2820   } else if (myexponent==0xff && mysignificand!=0) {
2821     // sign, exponent, significand meaningless
2822     category = fcNaN;
2823     *significandParts() = mysignificand;
2824   } else {
2825     category = fcNormal;
2826     exponent = myexponent - 127;  //bias
2827     *significandParts() = mysignificand;
2828     if (myexponent==0)    // denormal
2829       exponent = -126;
2830     else
2831       *significandParts() |= 0x800000; // integer bit
2832   }
2833 }
2834
2835 /// Treat api as containing the bits of a floating point number.  Currently
2836 /// we infer the floating point type from the size of the APInt.  The
2837 /// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful
2838 /// when the size is anything else).
2839 void
2840 APFloat::initFromAPInt(const APInt& api, bool isIEEE)
2841 {
2842   if (api.getBitWidth() == 32)
2843     return initFromFloatAPInt(api);
2844   else if (api.getBitWidth()==64)
2845     return initFromDoubleAPInt(api);
2846   else if (api.getBitWidth()==80)
2847     return initFromF80LongDoubleAPInt(api);
2848   else if (api.getBitWidth()==128 && !isIEEE)
2849     return initFromPPCDoubleDoubleAPInt(api);
2850   else
2851     assert(0);
2852 }
2853
2854 APFloat::APFloat(const APInt& api, bool isIEEE)
2855 {
2856   initFromAPInt(api, isIEEE);
2857 }
2858
2859 APFloat::APFloat(float f)
2860 {
2861   APInt api = APInt(32, 0);
2862   initFromAPInt(api.floatToBits(f));
2863 }
2864
2865 APFloat::APFloat(double d)
2866 {
2867   APInt api = APInt(64, 0);
2868   initFromAPInt(api.doubleToBits(d));
2869 }