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 correct in all cases.  */
1518 APFloat::opStatus
1519 APFloat::remainder(const APFloat &rhs)
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                           rmNearestTiesToEven, &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, rmNearestTiesToEven);
1543   assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
1544
1545   fs = subtract(V, rmNearestTiesToEven);
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 llvm frem (C fmod).
1555    This is not currently correct in all cases.  */
1556 APFloat::opStatus
1557 APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1558 {
1559   opStatus fs;
1560   APFloat V = *this;
1561   unsigned int origSign = sign;
1562
1563   assertArithmeticOK(*semantics);
1564   fs = V.divide(rhs, rmNearestTiesToEven);
1565   if (fs == opDivByZero)
1566     return fs;
1567
1568   int parts = partCount();
1569   integerPart *x = new integerPart[parts];
1570   bool ignored;
1571   fs = V.convertToInteger(x, parts * integerPartWidth, true,
1572                           rmTowardZero, &ignored);
1573   if (fs==opInvalidOp)
1574     return fs;
1575
1576   fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1577                                         rmNearestTiesToEven);
1578   assert(fs==opOK);   // should always work
1579
1580   fs = V.multiply(rhs, rounding_mode);
1581   assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
1582
1583   fs = subtract(V, rounding_mode);
1584   assert(fs==opOK || fs==opInexact);   // likewise
1585
1586   if (isZero())
1587     sign = origSign;    // IEEE754 requires this
1588   delete[] x;
1589   return fs;
1590 }
1591
1592 /* Normalized fused-multiply-add.  */
1593 APFloat::opStatus
1594 APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1595                           const APFloat &addend,
1596                           roundingMode rounding_mode)
1597 {
1598   opStatus fs;
1599
1600   assertArithmeticOK(*semantics);
1601
1602   /* Post-multiplication sign, before addition.  */
1603   sign ^= multiplicand.sign;
1604
1605   /* If and only if all arguments are normal do we need to do an
1606      extended-precision calculation.  */
1607   if(category == fcNormal
1608      && multiplicand.category == fcNormal
1609      && addend.category == fcNormal) {
1610     lostFraction lost_fraction;
1611
1612     lost_fraction = multiplySignificand(multiplicand, &addend);
1613     fs = normalize(rounding_mode, lost_fraction);
1614     if(lost_fraction != lfExactlyZero)
1615       fs = (opStatus) (fs | opInexact);
1616
1617     /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1618        positive zero unless rounding to minus infinity, except that
1619        adding two like-signed zeroes gives that zero.  */
1620     if(category == fcZero && sign != addend.sign)
1621       sign = (rounding_mode == rmTowardNegative);
1622   } else {
1623     fs = multiplySpecials(multiplicand);
1624
1625     /* FS can only be opOK or opInvalidOp.  There is no more work
1626        to do in the latter case.  The IEEE-754R standard says it is
1627        implementation-defined in this case whether, if ADDEND is a
1628        quiet NaN, we raise invalid op; this implementation does so.
1629
1630        If we need to do the addition we can do so with normal
1631        precision.  */
1632     if(fs == opOK)
1633       fs = addOrSubtract(addend, rounding_mode, false);
1634   }
1635
1636   return fs;
1637 }
1638
1639 /* Comparison requires normalized numbers.  */
1640 APFloat::cmpResult
1641 APFloat::compare(const APFloat &rhs) const
1642 {
1643   cmpResult result;
1644
1645   assertArithmeticOK(*semantics);
1646   assert(semantics == rhs.semantics);
1647
1648   switch(convolve(category, rhs.category)) {
1649   default:
1650     assert(0);
1651
1652   case convolve(fcNaN, fcZero):
1653   case convolve(fcNaN, fcNormal):
1654   case convolve(fcNaN, fcInfinity):
1655   case convolve(fcNaN, fcNaN):
1656   case convolve(fcZero, fcNaN):
1657   case convolve(fcNormal, fcNaN):
1658   case convolve(fcInfinity, fcNaN):
1659     return cmpUnordered;
1660
1661   case convolve(fcInfinity, fcNormal):
1662   case convolve(fcInfinity, fcZero):
1663   case convolve(fcNormal, fcZero):
1664     if(sign)
1665       return cmpLessThan;
1666     else
1667       return cmpGreaterThan;
1668
1669   case convolve(fcNormal, fcInfinity):
1670   case convolve(fcZero, fcInfinity):
1671   case convolve(fcZero, fcNormal):
1672     if(rhs.sign)
1673       return cmpGreaterThan;
1674     else
1675       return cmpLessThan;
1676
1677   case convolve(fcInfinity, fcInfinity):
1678     if(sign == rhs.sign)
1679       return cmpEqual;
1680     else if(sign)
1681       return cmpLessThan;
1682     else
1683       return cmpGreaterThan;
1684
1685   case convolve(fcZero, fcZero):
1686     return cmpEqual;
1687
1688   case convolve(fcNormal, fcNormal):
1689     break;
1690   }
1691
1692   /* Two normal numbers.  Do they have the same sign?  */
1693   if(sign != rhs.sign) {
1694     if(sign)
1695       result = cmpLessThan;
1696     else
1697       result = cmpGreaterThan;
1698   } else {
1699     /* Compare absolute values; invert result if negative.  */
1700     result = compareAbsoluteValue(rhs);
1701
1702     if(sign) {
1703       if(result == cmpLessThan)
1704         result = cmpGreaterThan;
1705       else if(result == cmpGreaterThan)
1706         result = cmpLessThan;
1707     }
1708   }
1709
1710   return result;
1711 }
1712
1713 /// APFloat::convert - convert a value of one floating point type to another.
1714 /// The return value corresponds to the IEEE754 exceptions.  *losesInfo
1715 /// records whether the transformation lost information, i.e. whether
1716 /// converting the result back to the original type will produce the
1717 /// original value (this is almost the same as return value==fsOK, but there
1718 /// are edge cases where this is not so).
1719
1720 APFloat::opStatus
1721 APFloat::convert(const fltSemantics &toSemantics,
1722                  roundingMode rounding_mode, bool *losesInfo)
1723 {
1724   lostFraction lostFraction;
1725   unsigned int newPartCount, oldPartCount;
1726   opStatus fs;
1727
1728   assertArithmeticOK(*semantics);
1729   assertArithmeticOK(toSemantics);
1730   lostFraction = lfExactlyZero;
1731   newPartCount = partCountForBits(toSemantics.precision + 1);
1732   oldPartCount = partCount();
1733
1734   /* Handle storage complications.  If our new form is wider,
1735      re-allocate our bit pattern into wider storage.  If it is
1736      narrower, we ignore the excess parts, but if narrowing to a
1737      single part we need to free the old storage.
1738      Be careful not to reference significandParts for zeroes
1739      and infinities, since it aborts.  */
1740   if (newPartCount > oldPartCount) {
1741     integerPart *newParts;
1742     newParts = new integerPart[newPartCount];
1743     APInt::tcSet(newParts, 0, newPartCount);
1744     if (category==fcNormal || category==fcNaN)
1745       APInt::tcAssign(newParts, significandParts(), oldPartCount);
1746     freeSignificand();
1747     significand.parts = newParts;
1748   } else if (newPartCount < oldPartCount) {
1749     /* Capture any lost fraction through truncation of parts so we get
1750        correct rounding whilst normalizing.  */
1751     if (category==fcNormal)
1752       lostFraction = lostFractionThroughTruncation
1753         (significandParts(), oldPartCount, toSemantics.precision);
1754     if (newPartCount == 1) {
1755         integerPart newPart = 0;
1756         if (category==fcNormal || category==fcNaN)
1757           newPart = significandParts()[0];
1758         freeSignificand();
1759         significand.part = newPart;
1760     }
1761   }
1762
1763   if(category == fcNormal) {
1764     /* Re-interpret our bit-pattern.  */
1765     exponent += toSemantics.precision - semantics->precision;
1766     semantics = &toSemantics;
1767     fs = normalize(rounding_mode, lostFraction);
1768     *losesInfo = (fs != opOK);
1769   } else if (category == fcNaN) {
1770     int shift = toSemantics.precision - semantics->precision;
1771     // Do this now so significandParts gets the right answer
1772     const fltSemantics *oldSemantics = semantics;
1773     semantics = &toSemantics;
1774     *losesInfo = false;
1775     // No normalization here, just truncate
1776     if (shift>0)
1777       APInt::tcShiftLeft(significandParts(), newPartCount, shift);
1778     else if (shift < 0) {
1779       unsigned ushift = -shift;
1780       // Figure out if we are losing information.  This happens
1781       // if are shifting out something other than 0s, or if the x87 long
1782       // double input did not have its integer bit set (pseudo-NaN), or if the
1783       // x87 long double input did not have its QNan bit set (because the x87
1784       // hardware sets this bit when converting a lower-precision NaN to
1785       // x87 long double).
1786       if (APInt::tcLSB(significandParts(), newPartCount) < ushift)
1787         *losesInfo = true;
1788       if (oldSemantics == &APFloat::x87DoubleExtended &&
1789           (!(*significandParts() & 0x8000000000000000ULL) ||
1790            !(*significandParts() & 0x4000000000000000ULL)))
1791         *losesInfo = true;
1792       APInt::tcShiftRight(significandParts(), newPartCount, ushift);
1793     }
1794     // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
1795     // does not give you back the same bits.  This is dubious, and we
1796     // don't currently do it.  You're really supposed to get
1797     // an invalid operation signal at runtime, but nobody does that.
1798     fs = opOK;
1799   } else {
1800     semantics = &toSemantics;
1801     fs = opOK;
1802     *losesInfo = false;
1803   }
1804
1805   return fs;
1806 }
1807
1808 /* Convert a floating point number to an integer according to the
1809    rounding mode.  If the rounded integer value is out of range this
1810    returns an invalid operation exception and the contents of the
1811    destination parts are unspecified.  If the rounded value is in
1812    range but the floating point number is not the exact integer, the C
1813    standard doesn't require an inexact exception to be raised.  IEEE
1814    854 does require it so we do that.
1815
1816    Note that for conversions to integer type the C standard requires
1817    round-to-zero to always be used.  */
1818 APFloat::opStatus
1819 APFloat::convertToSignExtendedInteger(integerPart *parts, unsigned int width,
1820                                       bool isSigned,
1821                                       roundingMode rounding_mode,
1822                                       bool *isExact) const
1823 {
1824   lostFraction lost_fraction;
1825   const integerPart *src;
1826   unsigned int dstPartsCount, truncatedBits;
1827
1828   assertArithmeticOK(*semantics);
1829
1830   *isExact = false;
1831
1832   /* Handle the three special cases first.  */
1833   if(category == fcInfinity || category == fcNaN)
1834     return opInvalidOp;
1835
1836   dstPartsCount = partCountForBits(width);
1837
1838   if(category == fcZero) {
1839     APInt::tcSet(parts, 0, dstPartsCount);
1840     // Negative zero can't be represented as an int.
1841     *isExact = !sign;
1842     return opOK;
1843   }
1844
1845   src = significandParts();
1846
1847   /* Step 1: place our absolute value, with any fraction truncated, in
1848      the destination.  */
1849   if (exponent < 0) {
1850     /* Our absolute value is less than one; truncate everything.  */
1851     APInt::tcSet(parts, 0, dstPartsCount);
1852     /* For exponent -1 the integer bit represents .5, look at that.
1853        For smaller exponents leftmost truncated bit is 0. */
1854     truncatedBits = semantics->precision -1U - exponent;
1855   } else {
1856     /* We want the most significant (exponent + 1) bits; the rest are
1857        truncated.  */
1858     unsigned int bits = exponent + 1U;
1859
1860     /* Hopelessly large in magnitude?  */
1861     if (bits > width)
1862       return opInvalidOp;
1863
1864     if (bits < semantics->precision) {
1865       /* We truncate (semantics->precision - bits) bits.  */
1866       truncatedBits = semantics->precision - bits;
1867       APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits);
1868     } else {
1869       /* We want at least as many bits as are available.  */
1870       APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0);
1871       APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision);
1872       truncatedBits = 0;
1873     }
1874   }
1875
1876   /* Step 2: work out any lost fraction, and increment the absolute
1877      value if we would round away from zero.  */
1878   if (truncatedBits) {
1879     lost_fraction = lostFractionThroughTruncation(src, partCount(),
1880                                                   truncatedBits);
1881     if (lost_fraction != lfExactlyZero
1882         && roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) {
1883       if (APInt::tcIncrement(parts, dstPartsCount))
1884         return opInvalidOp;     /* Overflow.  */
1885     }
1886   } else {
1887     lost_fraction = lfExactlyZero;
1888   }
1889
1890   /* Step 3: check if we fit in the destination.  */
1891   unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1;
1892
1893   if (sign) {
1894     if (!isSigned) {
1895       /* Negative numbers cannot be represented as unsigned.  */
1896       if (omsb != 0)
1897         return opInvalidOp;
1898     } else {
1899       /* It takes omsb bits to represent the unsigned integer value.
1900          We lose a bit for the sign, but care is needed as the
1901          maximally negative integer is a special case.  */
1902       if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb)
1903         return opInvalidOp;
1904
1905       /* This case can happen because of rounding.  */
1906       if (omsb > width)
1907         return opInvalidOp;
1908     }
1909
1910     APInt::tcNegate (parts, dstPartsCount);
1911   } else {
1912     if (omsb >= width + !isSigned)
1913       return opInvalidOp;
1914   }
1915
1916   if (lost_fraction == lfExactlyZero) {
1917     *isExact = true;
1918     return opOK;
1919   } else
1920     return opInexact;
1921 }
1922
1923 /* Same as convertToSignExtendedInteger, except we provide
1924    deterministic values in case of an invalid operation exception,
1925    namely zero for NaNs and the minimal or maximal value respectively
1926    for underflow or overflow.
1927    The *isExact output tells whether the result is exact, in the sense
1928    that converting it back to the original floating point type produces
1929    the original value.  This is almost equivalent to result==opOK,
1930    except for negative zeroes.
1931 */
1932 APFloat::opStatus
1933 APFloat::convertToInteger(integerPart *parts, unsigned int width,
1934                           bool isSigned,
1935                           roundingMode rounding_mode, bool *isExact) const
1936 {
1937   opStatus fs;
1938
1939   fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode,
1940                                     isExact);
1941
1942   if (fs == opInvalidOp) {
1943     unsigned int bits, dstPartsCount;
1944
1945     dstPartsCount = partCountForBits(width);
1946
1947     if (category == fcNaN)
1948       bits = 0;
1949     else if (sign)
1950       bits = isSigned;
1951     else
1952       bits = width - isSigned;
1953
1954     APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits);
1955     if (sign && isSigned)
1956       APInt::tcShiftLeft(parts, dstPartsCount, width - 1);
1957   }
1958
1959   return fs;
1960 }
1961
1962 /* Convert an unsigned integer SRC to a floating point number,
1963    rounding according to ROUNDING_MODE.  The sign of the floating
1964    point number is not modified.  */
1965 APFloat::opStatus
1966 APFloat::convertFromUnsignedParts(const integerPart *src,
1967                                   unsigned int srcCount,
1968                                   roundingMode rounding_mode)
1969 {
1970   unsigned int omsb, precision, dstCount;
1971   integerPart *dst;
1972   lostFraction lost_fraction;
1973
1974   assertArithmeticOK(*semantics);
1975   category = fcNormal;
1976   omsb = APInt::tcMSB(src, srcCount) + 1;
1977   dst = significandParts();
1978   dstCount = partCount();
1979   precision = semantics->precision;
1980
1981   /* We want the most significant PRECISON bits of SRC.  There may not
1982      be that many; extract what we can.  */
1983   if (precision <= omsb) {
1984     exponent = omsb - 1;
1985     lost_fraction = lostFractionThroughTruncation(src, srcCount,
1986                                                   omsb - precision);
1987     APInt::tcExtract(dst, dstCount, src, precision, omsb - precision);
1988   } else {
1989     exponent = precision - 1;
1990     lost_fraction = lfExactlyZero;
1991     APInt::tcExtract(dst, dstCount, src, omsb, 0);
1992   }
1993
1994   return normalize(rounding_mode, lost_fraction);
1995 }
1996
1997 APFloat::opStatus
1998 APFloat::convertFromAPInt(const APInt &Val,
1999                           bool isSigned,
2000                           roundingMode rounding_mode)
2001 {
2002   unsigned int partCount = Val.getNumWords();
2003   APInt api = Val;
2004
2005   sign = false;
2006   if (isSigned && api.isNegative()) {
2007     sign = true;
2008     api = -api;
2009   }
2010
2011   return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2012 }
2013
2014 /* Convert a two's complement integer SRC to a floating point number,
2015    rounding according to ROUNDING_MODE.  ISSIGNED is true if the
2016    integer is signed, in which case it must be sign-extended.  */
2017 APFloat::opStatus
2018 APFloat::convertFromSignExtendedInteger(const integerPart *src,
2019                                         unsigned int srcCount,
2020                                         bool isSigned,
2021                                         roundingMode rounding_mode)
2022 {
2023   opStatus status;
2024
2025   assertArithmeticOK(*semantics);
2026   if (isSigned
2027       && APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) {
2028     integerPart *copy;
2029
2030     /* If we're signed and negative negate a copy.  */
2031     sign = true;
2032     copy = new integerPart[srcCount];
2033     APInt::tcAssign(copy, src, srcCount);
2034     APInt::tcNegate(copy, srcCount);
2035     status = convertFromUnsignedParts(copy, srcCount, rounding_mode);
2036     delete [] copy;
2037   } else {
2038     sign = false;
2039     status = convertFromUnsignedParts(src, srcCount, rounding_mode);
2040   }
2041
2042   return status;
2043 }
2044
2045 /* FIXME: should this just take a const APInt reference?  */
2046 APFloat::opStatus
2047 APFloat::convertFromZeroExtendedInteger(const integerPart *parts,
2048                                         unsigned int width, bool isSigned,
2049                                         roundingMode rounding_mode)
2050 {
2051   unsigned int partCount = partCountForBits(width);
2052   APInt api = APInt(width, partCount, parts);
2053
2054   sign = false;
2055   if(isSigned && APInt::tcExtractBit(parts, width - 1)) {
2056     sign = true;
2057     api = -api;
2058   }
2059
2060   return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2061 }
2062
2063 APFloat::opStatus
2064 APFloat::convertFromHexadecimalString(const char *p,
2065                                       roundingMode rounding_mode)
2066 {
2067   lostFraction lost_fraction;
2068   integerPart *significand;
2069   unsigned int bitPos, partsCount;
2070   const char *dot, *firstSignificantDigit;
2071
2072   zeroSignificand();
2073   exponent = 0;
2074   category = fcNormal;
2075
2076   significand = significandParts();
2077   partsCount = partCount();
2078   bitPos = partsCount * integerPartWidth;
2079
2080   /* Skip leading zeroes and any (hexa)decimal point.  */
2081   p = skipLeadingZeroesAndAnyDot(p, &dot);
2082   firstSignificantDigit = p;
2083
2084   for(;;) {
2085     integerPart hex_value;
2086
2087     if(*p == '.') {
2088       assert(dot == 0);
2089       dot = p++;
2090     }
2091
2092     hex_value = hexDigitValue(*p);
2093     if(hex_value == -1U) {
2094       lost_fraction = lfExactlyZero;
2095       break;
2096     }
2097
2098     p++;
2099
2100     /* Store the number whilst 4-bit nibbles remain.  */
2101     if(bitPos) {
2102       bitPos -= 4;
2103       hex_value <<= bitPos % integerPartWidth;
2104       significand[bitPos / integerPartWidth] |= hex_value;
2105     } else {
2106       lost_fraction = trailingHexadecimalFraction(p, hex_value);
2107       while(hexDigitValue(*p) != -1U)
2108         p++;
2109       break;
2110     }
2111   }
2112
2113   /* Hex floats require an exponent but not a hexadecimal point.  */
2114   assert(*p == 'p' || *p == 'P');
2115
2116   /* Ignore the exponent if we are zero.  */
2117   if(p != firstSignificantDigit) {
2118     int expAdjustment;
2119
2120     /* Implicit hexadecimal point?  */
2121     if(!dot)
2122       dot = p;
2123
2124     /* Calculate the exponent adjustment implicit in the number of
2125        significant digits.  */
2126     expAdjustment = static_cast<int>(dot - firstSignificantDigit);
2127     if(expAdjustment < 0)
2128       expAdjustment++;
2129     expAdjustment = expAdjustment * 4 - 1;
2130
2131     /* Adjust for writing the significand starting at the most
2132        significant nibble.  */
2133     expAdjustment += semantics->precision;
2134     expAdjustment -= partsCount * integerPartWidth;
2135
2136     /* Adjust for the given exponent.  */
2137     exponent = totalExponent(p, expAdjustment);
2138   }
2139
2140   return normalize(rounding_mode, lost_fraction);
2141 }
2142
2143 APFloat::opStatus
2144 APFloat::roundSignificandWithExponent(const integerPart *decSigParts,
2145                                       unsigned sigPartCount, int exp,
2146                                       roundingMode rounding_mode)
2147 {
2148   unsigned int parts, pow5PartCount;
2149   fltSemantics calcSemantics = { 32767, -32767, 0, true };
2150   integerPart pow5Parts[maxPowerOfFiveParts];
2151   bool isNearest;
2152
2153   isNearest = (rounding_mode == rmNearestTiesToEven
2154                || rounding_mode == rmNearestTiesToAway);
2155
2156   parts = partCountForBits(semantics->precision + 11);
2157
2158   /* Calculate pow(5, abs(exp)).  */
2159   pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp);
2160
2161   for (;; parts *= 2) {
2162     opStatus sigStatus, powStatus;
2163     unsigned int excessPrecision, truncatedBits;
2164
2165     calcSemantics.precision = parts * integerPartWidth - 1;
2166     excessPrecision = calcSemantics.precision - semantics->precision;
2167     truncatedBits = excessPrecision;
2168
2169     APFloat decSig(calcSemantics, fcZero, sign);
2170     APFloat pow5(calcSemantics, fcZero, false);
2171
2172     sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount,
2173                                                 rmNearestTiesToEven);
2174     powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount,
2175                                               rmNearestTiesToEven);
2176     /* Add exp, as 10^n = 5^n * 2^n.  */
2177     decSig.exponent += exp;
2178
2179     lostFraction calcLostFraction;
2180     integerPart HUerr, HUdistance;
2181     unsigned int powHUerr;
2182
2183     if (exp >= 0) {
2184       /* multiplySignificand leaves the precision-th bit set to 1.  */
2185       calcLostFraction = decSig.multiplySignificand(pow5, NULL);
2186       powHUerr = powStatus != opOK;
2187     } else {
2188       calcLostFraction = decSig.divideSignificand(pow5);
2189       /* Denormal numbers have less precision.  */
2190       if (decSig.exponent < semantics->minExponent) {
2191         excessPrecision += (semantics->minExponent - decSig.exponent);
2192         truncatedBits = excessPrecision;
2193         if (excessPrecision > calcSemantics.precision)
2194           excessPrecision = calcSemantics.precision;
2195       }
2196       /* Extra half-ulp lost in reciprocal of exponent.  */
2197       powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2;
2198     }
2199
2200     /* Both multiplySignificand and divideSignificand return the
2201        result with the integer bit set.  */
2202     assert (APInt::tcExtractBit
2203             (decSig.significandParts(), calcSemantics.precision - 1) == 1);
2204
2205     HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK,
2206                        powHUerr);
2207     HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(),
2208                                       excessPrecision, isNearest);
2209
2210     /* Are we guaranteed to round correctly if we truncate?  */
2211     if (HUdistance >= HUerr) {
2212       APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(),
2213                        calcSemantics.precision - excessPrecision,
2214                        excessPrecision);
2215       /* Take the exponent of decSig.  If we tcExtract-ed less bits
2216          above we must adjust our exponent to compensate for the
2217          implicit right shift.  */
2218       exponent = (decSig.exponent + semantics->precision
2219                   - (calcSemantics.precision - excessPrecision));
2220       calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(),
2221                                                        decSig.partCount(),
2222                                                        truncatedBits);
2223       return normalize(rounding_mode, calcLostFraction);
2224     }
2225   }
2226 }
2227
2228 APFloat::opStatus
2229 APFloat::convertFromDecimalString(const char *p, roundingMode rounding_mode)
2230 {
2231   decimalInfo D;
2232   opStatus fs;
2233
2234   /* Scan the text.  */
2235   interpretDecimal(p, &D);
2236
2237   /* Handle the quick cases.  First the case of no significant digits,
2238      i.e. zero, and then exponents that are obviously too large or too
2239      small.  Writing L for log 10 / log 2, a number d.ddddd*10^exp
2240      definitely overflows if
2241
2242            (exp - 1) * L >= maxExponent
2243
2244      and definitely underflows to zero where
2245
2246            (exp + 1) * L <= minExponent - precision
2247
2248      With integer arithmetic the tightest bounds for L are
2249
2250            93/28 < L < 196/59            [ numerator <= 256 ]
2251            42039/12655 < L < 28738/8651  [ numerator <= 65536 ]
2252   */
2253
2254   if (decDigitValue(*D.firstSigDigit) >= 10U) {
2255     category = fcZero;
2256     fs = opOK;
2257   } else if ((D.normalizedExponent + 1) * 28738
2258              <= 8651 * (semantics->minExponent - (int) semantics->precision)) {
2259     /* Underflow to zero and round.  */
2260     zeroSignificand();
2261     fs = normalize(rounding_mode, lfLessThanHalf);
2262   } else if ((D.normalizedExponent - 1) * 42039
2263              >= 12655 * semantics->maxExponent) {
2264     /* Overflow and round.  */
2265     fs = handleOverflow(rounding_mode);
2266   } else {
2267     integerPart *decSignificand;
2268     unsigned int partCount;
2269
2270     /* A tight upper bound on number of bits required to hold an
2271        N-digit decimal integer is N * 196 / 59.  Allocate enough space
2272        to hold the full significand, and an extra part required by
2273        tcMultiplyPart.  */
2274     partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1;
2275     partCount = partCountForBits(1 + 196 * partCount / 59);
2276     decSignificand = new integerPart[partCount + 1];
2277     partCount = 0;
2278
2279     /* Convert to binary efficiently - we do almost all multiplication
2280        in an integerPart.  When this would overflow do we do a single
2281        bignum multiplication, and then revert again to multiplication
2282        in an integerPart.  */
2283     do {
2284       integerPart decValue, val, multiplier;
2285
2286       val = 0;
2287       multiplier = 1;
2288
2289       do {
2290         if (*p == '.')
2291           p++;
2292
2293         decValue = decDigitValue(*p++);
2294         multiplier *= 10;
2295         val = val * 10 + decValue;
2296         /* The maximum number that can be multiplied by ten with any
2297            digit added without overflowing an integerPart.  */
2298       } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10);
2299
2300       /* Multiply out the current part.  */
2301       APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val,
2302                             partCount, partCount + 1, false);
2303
2304       /* If we used another part (likely but not guaranteed), increase
2305          the count.  */
2306       if (decSignificand[partCount])
2307         partCount++;
2308     } while (p <= D.lastSigDigit);
2309
2310     category = fcNormal;
2311     fs = roundSignificandWithExponent(decSignificand, partCount,
2312                                       D.exponent, rounding_mode);
2313
2314     delete [] decSignificand;
2315   }
2316
2317   return fs;
2318 }
2319
2320 APFloat::opStatus
2321 APFloat::convertFromString(const char *p, roundingMode rounding_mode)
2322 {
2323   assertArithmeticOK(*semantics);
2324
2325   /* Handle a leading minus sign.  */
2326   if(*p == '-')
2327     sign = 1, p++;
2328   else
2329     sign = 0;
2330
2331   if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
2332     return convertFromHexadecimalString(p + 2, rounding_mode);
2333
2334   return convertFromDecimalString(p, rounding_mode);
2335 }
2336
2337 /* Write out a hexadecimal representation of the floating point value
2338    to DST, which must be of sufficient size, in the C99 form
2339    [-]0xh.hhhhp[+-]d.  Return the number of characters written,
2340    excluding the terminating NUL.
2341
2342    If UPPERCASE, the output is in upper case, otherwise in lower case.
2343
2344    HEXDIGITS digits appear altogether, rounding the value if
2345    necessary.  If HEXDIGITS is 0, the minimal precision to display the
2346    number precisely is used instead.  If nothing would appear after
2347    the decimal point it is suppressed.
2348
2349    The decimal exponent is always printed and has at least one digit.
2350    Zero values display an exponent of zero.  Infinities and NaNs
2351    appear as "infinity" or "nan" respectively.
2352
2353    The above rules are as specified by C99.  There is ambiguity about
2354    what the leading hexadecimal digit should be.  This implementation
2355    uses whatever is necessary so that the exponent is displayed as
2356    stored.  This implies the exponent will fall within the IEEE format
2357    range, and the leading hexadecimal digit will be 0 (for denormals),
2358    1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with
2359    any other digits zero).
2360 */
2361 unsigned int
2362 APFloat::convertToHexString(char *dst, unsigned int hexDigits,
2363                             bool upperCase, roundingMode rounding_mode) const
2364 {
2365   char *p;
2366
2367   assertArithmeticOK(*semantics);
2368
2369   p = dst;
2370   if (sign)
2371     *dst++ = '-';
2372
2373   switch (category) {
2374   case fcInfinity:
2375     memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1);
2376     dst += sizeof infinityL - 1;
2377     break;
2378
2379   case fcNaN:
2380     memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1);
2381     dst += sizeof NaNU - 1;
2382     break;
2383
2384   case fcZero:
2385     *dst++ = '0';
2386     *dst++ = upperCase ? 'X': 'x';
2387     *dst++ = '0';
2388     if (hexDigits > 1) {
2389       *dst++ = '.';
2390       memset (dst, '0', hexDigits - 1);
2391       dst += hexDigits - 1;
2392     }
2393     *dst++ = upperCase ? 'P': 'p';
2394     *dst++ = '0';
2395     break;
2396
2397   case fcNormal:
2398     dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode);
2399     break;
2400   }
2401
2402   *dst = 0;
2403
2404   return static_cast<unsigned int>(dst - p);
2405 }
2406
2407 /* Does the hard work of outputting the correctly rounded hexadecimal
2408    form of a normal floating point number with the specified number of
2409    hexadecimal digits.  If HEXDIGITS is zero the minimum number of
2410    digits necessary to print the value precisely is output.  */
2411 char *
2412 APFloat::convertNormalToHexString(char *dst, unsigned int hexDigits,
2413                                   bool upperCase,
2414                                   roundingMode rounding_mode) const
2415 {
2416   unsigned int count, valueBits, shift, partsCount, outputDigits;
2417   const char *hexDigitChars;
2418   const integerPart *significand;
2419   char *p;
2420   bool roundUp;
2421
2422   *dst++ = '0';
2423   *dst++ = upperCase ? 'X': 'x';
2424
2425   roundUp = false;
2426   hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower;
2427
2428   significand = significandParts();
2429   partsCount = partCount();
2430
2431   /* +3 because the first digit only uses the single integer bit, so
2432      we have 3 virtual zero most-significant-bits.  */
2433   valueBits = semantics->precision + 3;
2434   shift = integerPartWidth - valueBits % integerPartWidth;
2435
2436   /* The natural number of digits required ignoring trailing
2437      insignificant zeroes.  */
2438   outputDigits = (valueBits - significandLSB () + 3) / 4;
2439
2440   /* hexDigits of zero means use the required number for the
2441      precision.  Otherwise, see if we are truncating.  If we are,
2442      find out if we need to round away from zero.  */
2443   if (hexDigits) {
2444     if (hexDigits < outputDigits) {
2445       /* We are dropping non-zero bits, so need to check how to round.
2446          "bits" is the number of dropped bits.  */
2447       unsigned int bits;
2448       lostFraction fraction;
2449
2450       bits = valueBits - hexDigits * 4;
2451       fraction = lostFractionThroughTruncation (significand, partsCount, bits);
2452       roundUp = roundAwayFromZero(rounding_mode, fraction, bits);
2453     }
2454     outputDigits = hexDigits;
2455   }
2456
2457   /* Write the digits consecutively, and start writing in the location
2458      of the hexadecimal point.  We move the most significant digit
2459      left and add the hexadecimal point later.  */
2460   p = ++dst;
2461
2462   count = (valueBits + integerPartWidth - 1) / integerPartWidth;
2463
2464   while (outputDigits && count) {
2465     integerPart part;
2466
2467     /* Put the most significant integerPartWidth bits in "part".  */
2468     if (--count == partsCount)
2469       part = 0;  /* An imaginary higher zero part.  */
2470     else
2471       part = significand[count] << shift;
2472
2473     if (count && shift)
2474       part |= significand[count - 1] >> (integerPartWidth - shift);
2475
2476     /* Convert as much of "part" to hexdigits as we can.  */
2477     unsigned int curDigits = integerPartWidth / 4;
2478
2479     if (curDigits > outputDigits)
2480       curDigits = outputDigits;
2481     dst += partAsHex (dst, part, curDigits, hexDigitChars);
2482     outputDigits -= curDigits;
2483   }
2484
2485   if (roundUp) {
2486     char *q = dst;
2487
2488     /* Note that hexDigitChars has a trailing '0'.  */
2489     do {
2490       q--;
2491       *q = hexDigitChars[hexDigitValue (*q) + 1];
2492     } while (*q == '0');
2493     assert (q >= p);
2494   } else {
2495     /* Add trailing zeroes.  */
2496     memset (dst, '0', outputDigits);
2497     dst += outputDigits;
2498   }
2499
2500   /* Move the most significant digit to before the point, and if there
2501      is something after the decimal point add it.  This must come
2502      after rounding above.  */
2503   p[-1] = p[0];
2504   if (dst -1 == p)
2505     dst--;
2506   else
2507     p[0] = '.';
2508
2509   /* Finally output the exponent.  */
2510   *dst++ = upperCase ? 'P': 'p';
2511
2512   return writeSignedDecimal (dst, exponent);
2513 }
2514
2515 // For good performance it is desirable for different APFloats
2516 // to produce different integers.
2517 uint32_t
2518 APFloat::getHashValue() const
2519 {
2520   if (category==fcZero) return sign<<8 | semantics->precision ;
2521   else if (category==fcInfinity) return sign<<9 | semantics->precision;
2522   else if (category==fcNaN) return 1<<10 | semantics->precision;
2523   else {
2524     uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
2525     const integerPart* p = significandParts();
2526     for (int i=partCount(); i>0; i--, p++)
2527       hash ^= ((uint32_t)*p) ^ (uint32_t)((*p)>>32);
2528     return hash;
2529   }
2530 }
2531
2532 // Conversion from APFloat to/from host float/double.  It may eventually be
2533 // possible to eliminate these and have everybody deal with APFloats, but that
2534 // will take a while.  This approach will not easily extend to long double.
2535 // Current implementation requires integerPartWidth==64, which is correct at
2536 // the moment but could be made more general.
2537
2538 // Denormals have exponent minExponent in APFloat, but minExponent-1 in
2539 // the actual IEEE respresentations.  We compensate for that here.
2540
2541 APInt
2542 APFloat::convertF80LongDoubleAPFloatToAPInt() const
2543 {
2544   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended);
2545   assert (partCount()==2);
2546
2547   uint64_t myexponent, mysignificand;
2548
2549   if (category==fcNormal) {
2550     myexponent = exponent+16383; //bias
2551     mysignificand = significandParts()[0];
2552     if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
2553       myexponent = 0;   // denormal
2554   } else if (category==fcZero) {
2555     myexponent = 0;
2556     mysignificand = 0;
2557   } else if (category==fcInfinity) {
2558     myexponent = 0x7fff;
2559     mysignificand = 0x8000000000000000ULL;
2560   } else {
2561     assert(category == fcNaN && "Unknown category");
2562     myexponent = 0x7fff;
2563     mysignificand = significandParts()[0];
2564   }
2565
2566   uint64_t words[2];
2567   words[0] =  ((uint64_t)(sign & 1) << 63) |
2568               ((myexponent & 0x7fffLL) <<  48) |
2569               ((mysignificand >>16) & 0xffffffffffffLL);
2570   words[1] = mysignificand & 0xffff;
2571   return APInt(80, 2, words);
2572 }
2573
2574 APInt
2575 APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const
2576 {
2577   assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble);
2578   assert (partCount()==2);
2579
2580   uint64_t myexponent, mysignificand, myexponent2, mysignificand2;
2581
2582   if (category==fcNormal) {
2583     myexponent = exponent + 1023; //bias
2584     myexponent2 = exponent2 + 1023;
2585     mysignificand = significandParts()[0];
2586     mysignificand2 = significandParts()[1];
2587     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2588       myexponent = 0;   // denormal
2589     if (myexponent2==1 && !(mysignificand2 & 0x10000000000000LL))
2590       myexponent2 = 0;   // denormal
2591   } else if (category==fcZero) {
2592     myexponent = 0;
2593     mysignificand = 0;
2594     myexponent2 = 0;
2595     mysignificand2 = 0;
2596   } else if (category==fcInfinity) {
2597     myexponent = 0x7ff;
2598     myexponent2 = 0;
2599     mysignificand = 0;
2600     mysignificand2 = 0;
2601   } else {
2602     assert(category == fcNaN && "Unknown category");
2603     myexponent = 0x7ff;
2604     mysignificand = significandParts()[0];
2605     myexponent2 = exponent2;
2606     mysignificand2 = significandParts()[1];
2607   }
2608
2609   uint64_t words[2];
2610   words[0] =  ((uint64_t)(sign & 1) << 63) |
2611               ((myexponent & 0x7ff) <<  52) |
2612               (mysignificand & 0xfffffffffffffLL);
2613   words[1] =  ((uint64_t)(sign2 & 1) << 63) |
2614               ((myexponent2 & 0x7ff) <<  52) |
2615               (mysignificand2 & 0xfffffffffffffLL);
2616   return APInt(128, 2, words);
2617 }
2618
2619 APInt
2620 APFloat::convertDoubleAPFloatToAPInt() const
2621 {
2622   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2623   assert (partCount()==1);
2624
2625   uint64_t myexponent, mysignificand;
2626
2627   if (category==fcNormal) {
2628     myexponent = exponent+1023; //bias
2629     mysignificand = *significandParts();
2630     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2631       myexponent = 0;   // denormal
2632   } else if (category==fcZero) {
2633     myexponent = 0;
2634     mysignificand = 0;
2635   } else if (category==fcInfinity) {
2636     myexponent = 0x7ff;
2637     mysignificand = 0;
2638   } else {
2639     assert(category == fcNaN && "Unknown category!");
2640     myexponent = 0x7ff;
2641     mysignificand = *significandParts();
2642   }
2643
2644   return APInt(64, ((((uint64_t)(sign & 1) << 63) |
2645                      ((myexponent & 0x7ff) <<  52) |
2646                      (mysignificand & 0xfffffffffffffLL))));
2647 }
2648
2649 APInt
2650 APFloat::convertFloatAPFloatToAPInt() const
2651 {
2652   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2653   assert (partCount()==1);
2654
2655   uint32_t myexponent, mysignificand;
2656
2657   if (category==fcNormal) {
2658     myexponent = exponent+127; //bias
2659     mysignificand = (uint32_t)*significandParts();
2660     if (myexponent == 1 && !(mysignificand & 0x800000))
2661       myexponent = 0;   // denormal
2662   } else if (category==fcZero) {
2663     myexponent = 0;
2664     mysignificand = 0;
2665   } else if (category==fcInfinity) {
2666     myexponent = 0xff;
2667     mysignificand = 0;
2668   } else {
2669     assert(category == fcNaN && "Unknown category!");
2670     myexponent = 0xff;
2671     mysignificand = (uint32_t)*significandParts();
2672   }
2673
2674   return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
2675                     (mysignificand & 0x7fffff)));
2676 }
2677
2678 // This function creates an APInt that is just a bit map of the floating
2679 // point constant as it would appear in memory.  It is not a conversion,
2680 // and treating the result as a normal integer is unlikely to be useful.
2681
2682 APInt
2683 APFloat::bitcastToAPInt() const
2684 {
2685   if (semantics == (const llvm::fltSemantics*)&IEEEsingle)
2686     return convertFloatAPFloatToAPInt();
2687
2688   if (semantics == (const llvm::fltSemantics*)&IEEEdouble)
2689     return convertDoubleAPFloatToAPInt();
2690
2691   if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble)
2692     return convertPPCDoubleDoubleAPFloatToAPInt();
2693
2694   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended &&
2695          "unknown format!");
2696   return convertF80LongDoubleAPFloatToAPInt();
2697 }
2698
2699 float
2700 APFloat::convertToFloat() const
2701 {
2702   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2703   APInt api = bitcastToAPInt();
2704   return api.bitsToFloat();
2705 }
2706
2707 double
2708 APFloat::convertToDouble() const
2709 {
2710   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2711   APInt api = bitcastToAPInt();
2712   return api.bitsToDouble();
2713 }
2714
2715 /// Integer bit is explicit in this format.  Intel hardware (387 and later)
2716 /// does not support these bit patterns:
2717 ///  exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
2718 ///  exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
2719 ///  exponent = 0, integer bit 1 ("pseudodenormal")
2720 ///  exponent!=0 nor all 1's, integer bit 0 ("unnormal")
2721 /// At the moment, the first two are treated as NaNs, the second two as Normal.
2722 void
2723 APFloat::initFromF80LongDoubleAPInt(const APInt &api)
2724 {
2725   assert(api.getBitWidth()==80);
2726   uint64_t i1 = api.getRawData()[0];
2727   uint64_t i2 = api.getRawData()[1];
2728   uint64_t myexponent = (i1 >> 48) & 0x7fff;
2729   uint64_t mysignificand = ((i1 << 16) &  0xffffffffffff0000ULL) |
2730                           (i2 & 0xffff);
2731
2732   initialize(&APFloat::x87DoubleExtended);
2733   assert(partCount()==2);
2734
2735   sign = static_cast<unsigned int>(i1>>63);
2736   if (myexponent==0 && mysignificand==0) {
2737     // exponent, significand meaningless
2738     category = fcZero;
2739   } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
2740     // exponent, significand meaningless
2741     category = fcInfinity;
2742   } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
2743     // exponent meaningless
2744     category = fcNaN;
2745     significandParts()[0] = mysignificand;
2746     significandParts()[1] = 0;
2747   } else {
2748     category = fcNormal;
2749     exponent = myexponent - 16383;
2750     significandParts()[0] = mysignificand;
2751     significandParts()[1] = 0;
2752     if (myexponent==0)          // denormal
2753       exponent = -16382;
2754   }
2755 }
2756
2757 void
2758 APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api)
2759 {
2760   assert(api.getBitWidth()==128);
2761   uint64_t i1 = api.getRawData()[0];
2762   uint64_t i2 = api.getRawData()[1];
2763   uint64_t myexponent = (i1 >> 52) & 0x7ff;
2764   uint64_t mysignificand = i1 & 0xfffffffffffffLL;
2765   uint64_t myexponent2 = (i2 >> 52) & 0x7ff;
2766   uint64_t mysignificand2 = i2 & 0xfffffffffffffLL;
2767
2768   initialize(&APFloat::PPCDoubleDouble);
2769   assert(partCount()==2);
2770
2771   sign = static_cast<unsigned int>(i1>>63);
2772   sign2 = static_cast<unsigned int>(i2>>63);
2773   if (myexponent==0 && mysignificand==0) {
2774     // exponent, significand meaningless
2775     // exponent2 and significand2 are required to be 0; we don't check
2776     category = fcZero;
2777   } else if (myexponent==0x7ff && mysignificand==0) {
2778     // exponent, significand meaningless
2779     // exponent2 and significand2 are required to be 0; we don't check
2780     category = fcInfinity;
2781   } else if (myexponent==0x7ff && mysignificand!=0) {
2782     // exponent meaningless.  So is the whole second word, but keep it
2783     // for determinism.
2784     category = fcNaN;
2785     exponent2 = myexponent2;
2786     significandParts()[0] = mysignificand;
2787     significandParts()[1] = mysignificand2;
2788   } else {
2789     category = fcNormal;
2790     // Note there is no category2; the second word is treated as if it is
2791     // fcNormal, although it might be something else considered by itself.
2792     exponent = myexponent - 1023;
2793     exponent2 = myexponent2 - 1023;
2794     significandParts()[0] = mysignificand;
2795     significandParts()[1] = mysignificand2;
2796     if (myexponent==0)          // denormal
2797       exponent = -1022;
2798     else
2799       significandParts()[0] |= 0x10000000000000LL;  // integer bit
2800     if (myexponent2==0)
2801       exponent2 = -1022;
2802     else
2803       significandParts()[1] |= 0x10000000000000LL;  // integer bit
2804   }
2805 }
2806
2807 void
2808 APFloat::initFromDoubleAPInt(const APInt &api)
2809 {
2810   assert(api.getBitWidth()==64);
2811   uint64_t i = *api.getRawData();
2812   uint64_t myexponent = (i >> 52) & 0x7ff;
2813   uint64_t mysignificand = i & 0xfffffffffffffLL;
2814
2815   initialize(&APFloat::IEEEdouble);
2816   assert(partCount()==1);
2817
2818   sign = static_cast<unsigned int>(i>>63);
2819   if (myexponent==0 && mysignificand==0) {
2820     // exponent, significand meaningless
2821     category = fcZero;
2822   } else if (myexponent==0x7ff && mysignificand==0) {
2823     // exponent, significand meaningless
2824     category = fcInfinity;
2825   } else if (myexponent==0x7ff && mysignificand!=0) {
2826     // exponent meaningless
2827     category = fcNaN;
2828     *significandParts() = mysignificand;
2829   } else {
2830     category = fcNormal;
2831     exponent = myexponent - 1023;
2832     *significandParts() = mysignificand;
2833     if (myexponent==0)          // denormal
2834       exponent = -1022;
2835     else
2836       *significandParts() |= 0x10000000000000LL;  // integer bit
2837   }
2838 }
2839
2840 void
2841 APFloat::initFromFloatAPInt(const APInt & api)
2842 {
2843   assert(api.getBitWidth()==32);
2844   uint32_t i = (uint32_t)*api.getRawData();
2845   uint32_t myexponent = (i >> 23) & 0xff;
2846   uint32_t mysignificand = i & 0x7fffff;
2847
2848   initialize(&APFloat::IEEEsingle);
2849   assert(partCount()==1);
2850
2851   sign = i >> 31;
2852   if (myexponent==0 && mysignificand==0) {
2853     // exponent, significand meaningless
2854     category = fcZero;
2855   } else if (myexponent==0xff && mysignificand==0) {
2856     // exponent, significand meaningless
2857     category = fcInfinity;
2858   } else if (myexponent==0xff && mysignificand!=0) {
2859     // sign, exponent, significand meaningless
2860     category = fcNaN;
2861     *significandParts() = mysignificand;
2862   } else {
2863     category = fcNormal;
2864     exponent = myexponent - 127;  //bias
2865     *significandParts() = mysignificand;
2866     if (myexponent==0)    // denormal
2867       exponent = -126;
2868     else
2869       *significandParts() |= 0x800000; // integer bit
2870   }
2871 }
2872
2873 /// Treat api as containing the bits of a floating point number.  Currently
2874 /// we infer the floating point type from the size of the APInt.  The
2875 /// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful
2876 /// when the size is anything else).
2877 void
2878 APFloat::initFromAPInt(const APInt& api, bool isIEEE)
2879 {
2880   if (api.getBitWidth() == 32)
2881     return initFromFloatAPInt(api);
2882   else if (api.getBitWidth()==64)
2883     return initFromDoubleAPInt(api);
2884   else if (api.getBitWidth()==80)
2885     return initFromF80LongDoubleAPInt(api);
2886   else if (api.getBitWidth()==128 && !isIEEE)
2887     return initFromPPCDoubleDoubleAPInt(api);
2888   else
2889     assert(0);
2890 }
2891
2892 APFloat::APFloat(const APInt& api, bool isIEEE)
2893 {
2894   initFromAPInt(api, isIEEE);
2895 }
2896
2897 APFloat::APFloat(float f)
2898 {
2899   APInt api = APInt(32, 0);
2900   initFromAPInt(api.floatToBits(f));
2901 }
2902
2903 APFloat::APFloat(double d)
2904 {
2905   APInt api = APInt(64, 0);
2906   initFromAPInt(api.doubleToBits(d));
2907 }