lib/Support/APFloat.cpp

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