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(convertToAPInt());
 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
 792   assert(semantics == rhs.semantics);
 793
 794   precision = semantics->precision;
 795   newPartsCount = partCountForBits(precision * 2);
 796
 797   if(newPartsCount > 4)
 798     fullSignificand = new integerPart[newPartsCount];
 799   else
 800     fullSignificand = scratch;
 801
 802   lhsSignificand = significandParts();
 803   partsCount = partCount();
 804
 805   APInt::tcFullMultiply(fullSignificand, lhsSignificand,
 806                         rhs.significandParts(), partsCount, partsCount);
 807
 808   lost_fraction = lfExactlyZero;
 809   omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
 810   exponent += rhs.exponent;
 811
 812   if(addend) {
 813     Significand savedSignificand = significand;
 814     const fltSemantics *savedSemantics = semantics;
 815     fltSemantics extendedSemantics;
 816     opStatus status;
 817     unsigned int extendedPrecision;
 818
 819     /* Normalize our MSB.  */
 820     extendedPrecision = precision + precision - 1;
 821     if(omsb != extendedPrecision)
 822       {
 823         APInt::tcShiftLeft(fullSignificand, newPartsCount,
 824                            extendedPrecision - omsb);
 825         exponent -= extendedPrecision - omsb;
 826       }
 827
 828     /* Create new semantics.  */
 829     extendedSemantics = *semantics;
 830     extendedSemantics.precision = extendedPrecision;
 831
 832     if(newPartsCount == 1)
 833       significand.part = fullSignificand[0];
 834     else
 835       significand.parts = fullSignificand;
 836     semantics = &extendedSemantics;
 837
 838     APFloat extendedAddend(*addend);
 839     status = extendedAddend.convert(extendedSemantics, rmTowardZero);
 840     assert(status == opOK);
 841     lost_fraction = addOrSubtractSignificand(extendedAddend, false);
 842
 843     /* Restore our state.  */
 844     if(newPartsCount == 1)
 845       fullSignificand[0] = significand.part;
 846     significand = savedSignificand;
 847     semantics = savedSemantics;
 848
 849     omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
 850   }
 851
 852   exponent -= (precision - 1);
 853
 854   if(omsb > precision) {
 855     unsigned int bits, significantParts;
 856     lostFraction lf;
 857
 858     bits = omsb - precision;
 859     significantParts = partCountForBits(omsb);
 860     lf = shiftRight(fullSignificand, significantParts, bits);
 861     lost_fraction = combineLostFractions(lf, lost_fraction);
 862     exponent += bits;
 863   }
 864
 865   APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
 866
 867   if(newPartsCount > 4)
 868     delete [] fullSignificand;
 869
 870   return lost_fraction;
 871 }
 872
 873 /* Multiply the significands of LHS and RHS to DST.  */
 874 lostFraction
 875 APFloat::divideSignificand(const APFloat &rhs)
 876 {
 877   unsigned int bit, i, partsCount;
 878   const integerPart *rhsSignificand;
 879   integerPart *lhsSignificand, *dividend, *divisor;
 880   integerPart scratch[4];
 881   lostFraction lost_fraction;
 882
 883   assert(semantics == rhs.semantics);
 884
 885   lhsSignificand = significandParts();
 886   rhsSignificand = rhs.significandParts();
 887   partsCount = partCount();
 888
 889   if(partsCount > 2)
 890     dividend = new integerPart[partsCount * 2];
 891   else
 892     dividend = scratch;
 893
 894   divisor = dividend + partsCount;
 895
 896   /* Copy the dividend and divisor as they will be modified in-place.  */
 897   for(i = 0; i < partsCount; i++) {
 898     dividend[i] = lhsSignificand[i];
 899     divisor[i] = rhsSignificand[i];
 900     lhsSignificand[i] = 0;
 901   }
 902
 903   exponent -= rhs.exponent;
 904
 905   unsigned int precision = semantics->precision;
 906
 907   /* Normalize the divisor.  */
 908   bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
 909   if(bit) {
 910     exponent += bit;
 911     APInt::tcShiftLeft(divisor, partsCount, bit);
 912   }
 913
 914   /* Normalize the dividend.  */
 915   bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
 916   if(bit) {
 917     exponent -= bit;
 918     APInt::tcShiftLeft(dividend, partsCount, bit);
 919   }
 920
 921   /* Ensure the dividend >= divisor initially for the loop below.
 922      Incidentally, this means that the division loop below is
 923      guaranteed to set the integer bit to one.  */
 924   if(APInt::tcCompare(dividend, divisor, partsCount) < 0) {
 925     exponent--;
 926     APInt::tcShiftLeft(dividend, partsCount, 1);
 927     assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
 928   }
 929
 930   /* Long division.  */
 931   for(bit = precision; bit; bit -= 1) {
 932     if(APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
 933       APInt::tcSubtract(dividend, divisor, 0, partsCount);
 934       APInt::tcSetBit(lhsSignificand, bit - 1);
 935     }
 936
 937     APInt::tcShiftLeft(dividend, partsCount, 1);
 938   }
 939
 940   /* Figure out the lost fraction.  */
 941   int cmp = APInt::tcCompare(dividend, divisor, partsCount);
 942
 943   if(cmp > 0)
 944     lost_fraction = lfMoreThanHalf;
 945   else if(cmp == 0)
 946     lost_fraction = lfExactlyHalf;
 947   else if(APInt::tcIsZero(dividend, partsCount))
 948     lost_fraction = lfExactlyZero;
 949   else
 950     lost_fraction = lfLessThanHalf;
 951
 952   if(partsCount > 2)
 953     delete [] dividend;
 954
 955   return lost_fraction;
 956 }
 957
 958 unsigned int
 959 APFloat::significandMSB() const
 960 {
 961   return APInt::tcMSB(significandParts(), partCount());
 962 }
 963
 964 unsigned int
 965 APFloat::significandLSB() const
 966 {
 967   return APInt::tcLSB(significandParts(), partCount());
 968 }
 969
 970 /* Note that a zero result is NOT normalized to fcZero.  */
 971 lostFraction
 972 APFloat::shiftSignificandRight(unsigned int bits)
 973 {
 974   /* Our exponent should not overflow.  */
 975   assert((exponent_t) (exponent + bits) >= exponent);
 976
 977   exponent += bits;
 978
 979   return shiftRight(significandParts(), partCount(), bits);
 980 }
 981
 982 /* Shift the significand left BITS bits, subtract BITS from its exponent.  */
 983 void
 984 APFloat::shiftSignificandLeft(unsigned int bits)
 985 {
 986   assert(bits < semantics->precision);
 987
 988   if(bits) {
 989     unsigned int partsCount = partCount();
 990
 991     APInt::tcShiftLeft(significandParts(), partsCount, bits);
 992     exponent -= bits;
 993
 994     assert(!APInt::tcIsZero(significandParts(), partsCount));
 995   }
 996 }
 997
 998 APFloat::cmpResult
 999 APFloat::compareAbsoluteValue(const APFloat &rhs) const
1000 {
1001   int compare;
1002
1003   assert(semantics == rhs.semantics);
1004   assert(category == fcNormal);
1005   assert(rhs.category == fcNormal);
1006
1007   compare = exponent - rhs.exponent;
1008
1009   /* If exponents are equal, do an unsigned bignum comparison of the
1010      significands.  */
1011   if(compare == 0)
1012     compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
1013                                partCount());
1014
1015   if(compare > 0)
1016     return cmpGreaterThan;
1017   else if(compare < 0)
1018     return cmpLessThan;
1019   else
1020     return cmpEqual;
1021 }
1022
1023 /* Handle overflow.  Sign is preserved.  We either become infinity or
1024    the largest finite number.  */
1025 APFloat::opStatus
1026 APFloat::handleOverflow(roundingMode rounding_mode)
1027 {
1028   /* Infinity?  */
1029   if(rounding_mode == rmNearestTiesToEven
1030      || rounding_mode == rmNearestTiesToAway
1031      || (rounding_mode == rmTowardPositive && !sign)
1032      || (rounding_mode == rmTowardNegative && sign))
1033     {
1034       category = fcInfinity;
1035       return (opStatus) (opOverflow | opInexact);
1036     }
1037
1038   /* Otherwise we become the largest finite number.  */
1039   category = fcNormal;
1040   exponent = semantics->maxExponent;
1041   APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
1042                                    semantics->precision);
1043
1044   return opInexact;
1045 }
1046
1047 /* Returns TRUE if, when truncating the current number, with BIT the
1048    new LSB, with the given lost fraction and rounding mode, the result
1049    would need to be rounded away from zero (i.e., by increasing the
1050    signficand).  This routine must work for fcZero of both signs, and
1051    fcNormal numbers.  */
1052 bool
1053 APFloat::roundAwayFromZero(roundingMode rounding_mode,
1054                            lostFraction lost_fraction,
1055                            unsigned int bit) const
1056 {
1057   /* NaNs and infinities should not have lost fractions.  */
1058   assert(category == fcNormal || category == fcZero);
1059
1060   /* Current callers never pass this so we don't handle it.  */
1061   assert(lost_fraction != lfExactlyZero);
1062
1063   switch(rounding_mode) {
1064   default:
1065     assert(0);
1066
1067   case rmNearestTiesToAway:
1068     return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
1069
1070   case rmNearestTiesToEven:
1071     if(lost_fraction == lfMoreThanHalf)
1072       return true;
1073
1074     /* Our zeroes don't have a significand to test.  */
1075     if(lost_fraction == lfExactlyHalf && category != fcZero)
1076       return APInt::tcExtractBit(significandParts(), bit);
1077
1078     return false;
1079
1080   case rmTowardZero:
1081     return false;
1082
1083   case rmTowardPositive:
1084     return sign == false;
1085
1086   case rmTowardNegative:
1087     return sign == true;
1088   }
1089 }
1090
1091 APFloat::opStatus
1092 APFloat::normalize(roundingMode rounding_mode,
1093                    lostFraction lost_fraction)
1094 {
1095   unsigned int omsb;                /* One, not zero, based MSB.  */
1096   int exponentChange;
1097
1098   if(category != fcNormal)
1099     return opOK;
1100
1101   /* Before rounding normalize the exponent of fcNormal numbers.  */
1102   omsb = significandMSB() + 1;
1103
1104   if(omsb) {
1105     /* OMSB is numbered from 1.  We want to place it in the integer
1106        bit numbered PRECISON if possible, with a compensating change in
1107        the exponent.  */
1108     exponentChange = omsb - semantics->precision;
1109
1110     /* If the resulting exponent is too high, overflow according to
1111        the rounding mode.  */
1112     if(exponent + exponentChange > semantics->maxExponent)
1113       return handleOverflow(rounding_mode);
1114
1115     /* Subnormal numbers have exponent minExponent, and their MSB
1116        is forced based on that.  */
1117     if(exponent + exponentChange < semantics->minExponent)
1118       exponentChange = semantics->minExponent - exponent;
1119
1120     /* Shifting left is easy as we don't lose precision.  */
1121     if(exponentChange < 0) {
1122       assert(lost_fraction == lfExactlyZero);
1123
1124       shiftSignificandLeft(-exponentChange);
1125
1126       return opOK;
1127     }
1128
1129     if(exponentChange > 0) {
1130       lostFraction lf;
1131
1132       /* Shift right and capture any new lost fraction.  */
1133       lf = shiftSignificandRight(exponentChange);
1134
1135       lost_fraction = combineLostFractions(lf, lost_fraction);
1136
1137       /* Keep OMSB up-to-date.  */
1138       if(omsb > (unsigned) exponentChange)
1139         omsb -= exponentChange;
1140       else
1141         omsb = 0;
1142     }
1143   }
1144
1145   /* Now round the number according to rounding_mode given the lost
1146      fraction.  */
1147
1148   /* As specified in IEEE 754, since we do not trap we do not report
1149      underflow for exact results.  */
1150   if(lost_fraction == lfExactlyZero) {
1151     /* Canonicalize zeroes.  */
1152     if(omsb == 0)
1153       category = fcZero;
1154
1155     return opOK;
1156   }
1157
1158   /* Increment the significand if we're rounding away from zero.  */
1159   if(roundAwayFromZero(rounding_mode, lost_fraction, 0)) {
1160     if(omsb == 0)
1161       exponent = semantics->minExponent;
1162
1163     incrementSignificand();
1164     omsb = significandMSB() + 1;
1165
1166     /* Did the significand increment overflow?  */
1167     if(omsb == (unsigned) semantics->precision + 1) {
1168       /* Renormalize by incrementing the exponent and shifting our
1169          significand right one.  However if we already have the
1170          maximum exponent we overflow to infinity.  */
1171       if(exponent == semantics->maxExponent) {
1172         category = fcInfinity;
1173
1174         return (opStatus) (opOverflow | opInexact);
1175       }
1176
1177       shiftSignificandRight(1);
1178
1179       return opInexact;
1180     }
1181   }
1182
1183   /* The normal case - we were and are not denormal, and any
1184      significand increment above didn't overflow.  */
1185   if(omsb == semantics->precision)
1186     return opInexact;
1187
1188   /* We have a non-zero denormal.  */
1189   assert(omsb < semantics->precision);
1190
1191   /* Canonicalize zeroes.  */
1192   if(omsb == 0)
1193     category = fcZero;
1194
1195   /* The fcZero case is a denormal that underflowed to zero.  */
1196   return (opStatus) (opUnderflow | opInexact);
1197 }
1198
1199 APFloat::opStatus
1200 APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
1201 {
1202   switch(convolve(category, rhs.category)) {
1203   default:
1204     assert(0);
1205
1206   case convolve(fcNaN, fcZero):
1207   case convolve(fcNaN, fcNormal):
1208   case convolve(fcNaN, fcInfinity):
1209   case convolve(fcNaN, fcNaN):
1210   case convolve(fcNormal, fcZero):
1211   case convolve(fcInfinity, fcNormal):
1212   case convolve(fcInfinity, fcZero):
1213     return opOK;
1214
1215   case convolve(fcZero, fcNaN):
1216   case convolve(fcNormal, fcNaN):
1217   case convolve(fcInfinity, fcNaN):
1218     category = fcNaN;
1219     copySignificand(rhs);
1220     return opOK;
1221
1222   case convolve(fcNormal, fcInfinity):
1223   case convolve(fcZero, fcInfinity):
1224     category = fcInfinity;
1225     sign = rhs.sign ^ subtract;
1226     return opOK;
1227
1228   case convolve(fcZero, fcNormal):
1229     assign(rhs);
1230     sign = rhs.sign ^ subtract;
1231     return opOK;
1232
1233   case convolve(fcZero, fcZero):
1234     /* Sign depends on rounding mode; handled by caller.  */
1235     return opOK;
1236
1237   case convolve(fcInfinity, fcInfinity):
1238     /* Differently signed infinities can only be validly
1239        subtracted.  */
1240     if((sign ^ rhs.sign) != subtract) {
1241       makeNaN();
1242       return opInvalidOp;
1243     }
1244
1245     return opOK;
1246
1247   case convolve(fcNormal, fcNormal):
1248     return opDivByZero;
1249   }
1250 }
1251
1252 /* Add or subtract two normal numbers.  */
1253 lostFraction
1254 APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
1255 {
1256   integerPart carry;
1257   lostFraction lost_fraction;
1258   int bits;
1259
1260   /* Determine if the operation on the absolute values is effectively
1261      an addition or subtraction.  */
1262   subtract ^= (sign ^ rhs.sign) ? true : false;
1263
1264   /* Are we bigger exponent-wise than the RHS?  */
1265   bits = exponent - rhs.exponent;
1266
1267   /* Subtraction is more subtle than one might naively expect.  */
1268   if(subtract) {
1269     APFloat temp_rhs(rhs);
1270     bool reverse;
1271
1272     if (bits == 0) {
1273       reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
1274       lost_fraction = lfExactlyZero;
1275     } else if (bits > 0) {
1276       lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
1277       shiftSignificandLeft(1);
1278       reverse = false;
1279     } else {
1280       lost_fraction = shiftSignificandRight(-bits - 1);
1281       temp_rhs.shiftSignificandLeft(1);
1282       reverse = true;
1283     }
1284
1285     if (reverse) {
1286       carry = temp_rhs.subtractSignificand
1287         (*this, lost_fraction != lfExactlyZero);
1288       copySignificand(temp_rhs);
1289       sign = !sign;
1290     } else {
1291       carry = subtractSignificand
1292         (temp_rhs, lost_fraction != lfExactlyZero);
1293     }
1294
1295     /* Invert the lost fraction - it was on the RHS and
1296        subtracted.  */
1297     if(lost_fraction == lfLessThanHalf)
1298       lost_fraction = lfMoreThanHalf;
1299     else if(lost_fraction == lfMoreThanHalf)
1300       lost_fraction = lfLessThanHalf;
1301
1302     /* The code above is intended to ensure that no borrow is
1303        necessary.  */
1304     assert(!carry);
1305   } else {
1306     if(bits > 0) {
1307       APFloat temp_rhs(rhs);
1308
1309       lost_fraction = temp_rhs.shiftSignificandRight(bits);
1310       carry = addSignificand(temp_rhs);
1311     } else {
1312       lost_fraction = shiftSignificandRight(-bits);
1313       carry = addSignificand(rhs);
1314     }
1315
1316     /* We have a guard bit; generating a carry cannot happen.  */
1317     assert(!carry);
1318   }
1319
1320   return lost_fraction;
1321 }
1322
1323 APFloat::opStatus
1324 APFloat::multiplySpecials(const APFloat &rhs)
1325 {
1326   switch(convolve(category, rhs.category)) {
1327   default:
1328     assert(0);
1329
1330   case convolve(fcNaN, fcZero):
1331   case convolve(fcNaN, fcNormal):
1332   case convolve(fcNaN, fcInfinity):
1333   case convolve(fcNaN, fcNaN):
1334     return opOK;
1335
1336   case convolve(fcZero, fcNaN):
1337   case convolve(fcNormal, fcNaN):
1338   case convolve(fcInfinity, fcNaN):
1339     category = fcNaN;
1340     copySignificand(rhs);
1341     return opOK;
1342
1343   case convolve(fcNormal, fcInfinity):
1344   case convolve(fcInfinity, fcNormal):
1345   case convolve(fcInfinity, fcInfinity):
1346     category = fcInfinity;
1347     return opOK;
1348
1349   case convolve(fcZero, fcNormal):
1350   case convolve(fcNormal, fcZero):
1351   case convolve(fcZero, fcZero):
1352     category = fcZero;
1353     return opOK;
1354
1355   case convolve(fcZero, fcInfinity):
1356   case convolve(fcInfinity, fcZero):
1357     makeNaN();
1358     return opInvalidOp;
1359
1360   case convolve(fcNormal, fcNormal):
1361     return opOK;
1362   }
1363 }
1364
1365 APFloat::opStatus
1366 APFloat::divideSpecials(const APFloat &rhs)
1367 {
1368   switch(convolve(category, rhs.category)) {
1369   default:
1370     assert(0);
1371
1372   case convolve(fcNaN, fcZero):
1373   case convolve(fcNaN, fcNormal):
1374   case convolve(fcNaN, fcInfinity):
1375   case convolve(fcNaN, fcNaN):
1376   case convolve(fcInfinity, fcZero):
1377   case convolve(fcInfinity, fcNormal):
1378   case convolve(fcZero, fcInfinity):
1379   case convolve(fcZero, fcNormal):
1380     return opOK;
1381
1382   case convolve(fcZero, fcNaN):
1383   case convolve(fcNormal, fcNaN):
1384   case convolve(fcInfinity, fcNaN):
1385     category = fcNaN;
1386     copySignificand(rhs);
1387     return opOK;
1388
1389   case convolve(fcNormal, fcInfinity):
1390     category = fcZero;
1391     return opOK;
1392
1393   case convolve(fcNormal, fcZero):
1394     category = fcInfinity;
1395     return opDivByZero;
1396
1397   case convolve(fcInfinity, fcInfinity):
1398   case convolve(fcZero, fcZero):
1399     makeNaN();
1400     return opInvalidOp;
1401
1402   case convolve(fcNormal, fcNormal):
1403     return opOK;
1404   }
1405 }
1406
1407 /* Change sign.  */
1408 void
1409 APFloat::changeSign()
1410 {
1411   /* Look mummy, this one's easy.  */
1412   sign = !sign;
1413 }
1414
1415 void
1416 APFloat::clearSign()
1417 {
1418   /* So is this one. */
1419   sign = 0;
1420 }
1421
1422 void
1423 APFloat::copySign(const APFloat &rhs)
1424 {
1425   /* And this one. */
1426   sign = rhs.sign;
1427 }
1428
1429 /* Normalized addition or subtraction.  */
1430 APFloat::opStatus
1431 APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1432                        bool subtract)
1433 {
1434   opStatus fs;
1435
1436   assertArithmeticOK(*semantics);
1437
1438   fs = addOrSubtractSpecials(rhs, subtract);
1439
1440   /* This return code means it was not a simple case.  */
1441   if(fs == opDivByZero) {
1442     lostFraction lost_fraction;
1443
1444     lost_fraction = addOrSubtractSignificand(rhs, subtract);
1445     fs = normalize(rounding_mode, lost_fraction);
1446
1447     /* Can only be zero if we lost no fraction.  */
1448     assert(category != fcZero || lost_fraction == lfExactlyZero);
1449   }
1450
1451   /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1452      positive zero unless rounding to minus infinity, except that
1453      adding two like-signed zeroes gives that zero.  */
1454   if(category == fcZero) {
1455     if(rhs.category != fcZero || (sign == rhs.sign) == subtract)
1456       sign = (rounding_mode == rmTowardNegative);
1457   }
1458
1459   return fs;
1460 }
1461
1462 /* Normalized addition.  */
1463 APFloat::opStatus
1464 APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1465 {
1466   return addOrSubtract(rhs, rounding_mode, false);
1467 }
1468
1469 /* Normalized subtraction.  */
1470 APFloat::opStatus
1471 APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1472 {
1473   return addOrSubtract(rhs, rounding_mode, true);
1474 }
1475
1476 /* Normalized multiply.  */
1477 APFloat::opStatus
1478 APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1479 {
1480   opStatus fs;
1481
1482   assertArithmeticOK(*semantics);
1483   sign ^= rhs.sign;
1484   fs = multiplySpecials(rhs);
1485
1486   if(category == fcNormal) {
1487     lostFraction lost_fraction = multiplySignificand(rhs, 0);
1488     fs = normalize(rounding_mode, lost_fraction);
1489     if(lost_fraction != lfExactlyZero)
1490       fs = (opStatus) (fs | opInexact);
1491   }
1492
1493   return fs;
1494 }
1495
1496 /* Normalized divide.  */
1497 APFloat::opStatus
1498 APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1499 {
1500   opStatus fs;
1501
1502   assertArithmeticOK(*semantics);
1503   sign ^= rhs.sign;
1504   fs = divideSpecials(rhs);
1505
1506   if(category == fcNormal) {
1507     lostFraction lost_fraction = divideSignificand(rhs);
1508     fs = normalize(rounding_mode, lost_fraction);
1509     if(lost_fraction != lfExactlyZero)
1510       fs = (opStatus) (fs | opInexact);
1511   }
1512
1513   return fs;
1514 }
1515
1516 /* Normalized remainder.  This is not currently doing TRT.  */
1517 APFloat::opStatus
1518 APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1519 {
1520   opStatus fs;
1521   APFloat V = *this;
1522   unsigned int origSign = sign;
1523
1524   assertArithmeticOK(*semantics);
1525   fs = V.divide(rhs, rmNearestTiesToEven);
1526   if (fs == opDivByZero)
1527     return fs;
1528
1529   int parts = partCount();
1530   integerPart *x = new integerPart[parts];
1531   fs = V.convertToInteger(x, parts * integerPartWidth, true,
1532                           rmNearestTiesToEven);
1533   if (fs==opInvalidOp)
1534     return fs;
1535
1536   fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1537                                         rmNearestTiesToEven);
1538   assert(fs==opOK);   // should always work
1539
1540   fs = V.multiply(rhs, rounding_mode);
1541   assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
1542
1543   fs = subtract(V, rounding_mode);
1544   assert(fs==opOK || fs==opInexact);   // likewise
1545
1546   if (isZero())
1547     sign = origSign;    // IEEE754 requires this
1548   delete[] x;
1549   return fs;
1550 }
1551
1552 /* Normalized fused-multiply-add.  */
1553 APFloat::opStatus
1554 APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1555                           const APFloat &addend,
1556                           roundingMode rounding_mode)
1557 {
1558   opStatus fs;
1559
1560   assertArithmeticOK(*semantics);
1561
1562   /* Post-multiplication sign, before addition.  */
1563   sign ^= multiplicand.sign;
1564
1565   /* If and only if all arguments are normal do we need to do an
1566      extended-precision calculation.  */
1567   if(category == fcNormal
1568      && multiplicand.category == fcNormal
1569      && addend.category == fcNormal) {
1570     lostFraction lost_fraction;
1571
1572     lost_fraction = multiplySignificand(multiplicand, &addend);
1573     fs = normalize(rounding_mode, lost_fraction);
1574     if(lost_fraction != lfExactlyZero)
1575       fs = (opStatus) (fs | opInexact);
1576
1577     /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1578        positive zero unless rounding to minus infinity, except that
1579        adding two like-signed zeroes gives that zero.  */
1580     if(category == fcZero && sign != addend.sign)
1581       sign = (rounding_mode == rmTowardNegative);
1582   } else {
1583     fs = multiplySpecials(multiplicand);
1584
1585     /* FS can only be opOK or opInvalidOp.  There is no more work
1586        to do in the latter case.  The IEEE-754R standard says it is
1587        implementation-defined in this case whether, if ADDEND is a
1588        quiet NaN, we raise invalid op; this implementation does so.
1589
1590        If we need to do the addition we can do so with normal
1591        precision.  */
1592     if(fs == opOK)
1593       fs = addOrSubtract(addend, rounding_mode, false);
1594   }
1595
1596   return fs;
1597 }
1598
1599 /* Comparison requires normalized numbers.  */
1600 APFloat::cmpResult
1601 APFloat::compare(const APFloat &rhs) const
1602 {
1603   cmpResult result;
1604
1605   assertArithmeticOK(*semantics);
1606   assert(semantics == rhs.semantics);
1607
1608   switch(convolve(category, rhs.category)) {
1609   default:
1610     assert(0);
1611
1612   case convolve(fcNaN, fcZero):
1613   case convolve(fcNaN, fcNormal):
1614   case convolve(fcNaN, fcInfinity):
1615   case convolve(fcNaN, fcNaN):
1616   case convolve(fcZero, fcNaN):
1617   case convolve(fcNormal, fcNaN):
1618   case convolve(fcInfinity, fcNaN):
1619     return cmpUnordered;
1620
1621   case convolve(fcInfinity, fcNormal):
1622   case convolve(fcInfinity, fcZero):
1623   case convolve(fcNormal, fcZero):
1624     if(sign)
1625       return cmpLessThan;
1626     else
1627       return cmpGreaterThan;
1628
1629   case convolve(fcNormal, fcInfinity):
1630   case convolve(fcZero, fcInfinity):
1631   case convolve(fcZero, fcNormal):
1632     if(rhs.sign)
1633       return cmpGreaterThan;
1634     else
1635       return cmpLessThan;
1636
1637   case convolve(fcInfinity, fcInfinity):
1638     if(sign == rhs.sign)
1639       return cmpEqual;
1640     else if(sign)
1641       return cmpLessThan;
1642     else
1643       return cmpGreaterThan;
1644
1645   case convolve(fcZero, fcZero):
1646     return cmpEqual;
1647
1648   case convolve(fcNormal, fcNormal):
1649     break;
1650   }
1651
1652   /* Two normal numbers.  Do they have the same sign?  */
1653   if(sign != rhs.sign) {
1654     if(sign)
1655       result = cmpLessThan;
1656     else
1657       result = cmpGreaterThan;
1658   } else {
1659     /* Compare absolute values; invert result if negative.  */
1660     result = compareAbsoluteValue(rhs);
1661
1662     if(sign) {
1663       if(result == cmpLessThan)
1664         result = cmpGreaterThan;
1665       else if(result == cmpGreaterThan)
1666         result = cmpLessThan;
1667     }
1668   }
1669
1670   return result;
1671 }
1672
1673 APFloat::opStatus
1674 APFloat::convert(const fltSemantics &toSemantics,
1675                  roundingMode rounding_mode)
1676 {
1677   lostFraction lostFraction;
1678   unsigned int newPartCount, oldPartCount;
1679   opStatus fs;
1680
1681   assertArithmeticOK(*semantics);
1682   assertArithmeticOK(toSemantics);
1683   lostFraction = lfExactlyZero;
1684   newPartCount = partCountForBits(toSemantics.precision + 1);
1685   oldPartCount = partCount();
1686
1687   /* Handle storage complications.  If our new form is wider,
1688      re-allocate our bit pattern into wider storage.  If it is
1689      narrower, we ignore the excess parts, but if narrowing to a
1690      single part we need to free the old storage.
1691      Be careful not to reference significandParts for zeroes
1692      and infinities, since it aborts.  */
1693   if (newPartCount > oldPartCount) {
1694     integerPart *newParts;
1695     newParts = new integerPart[newPartCount];
1696     APInt::tcSet(newParts, 0, newPartCount);
1697     if (category==fcNormal || category==fcNaN)
1698       APInt::tcAssign(newParts, significandParts(), oldPartCount);
1699     freeSignificand();
1700     significand.parts = newParts;
1701   } else if (newPartCount < oldPartCount) {
1702     /* Capture any lost fraction through truncation of parts so we get
1703        correct rounding whilst normalizing.  */
1704     if (category==fcNormal)
1705       lostFraction = lostFractionThroughTruncation
1706         (significandParts(), oldPartCount, toSemantics.precision);
1707     if (newPartCount == 1) {
1708         integerPart newPart = 0;
1709         if (category==fcNormal || category==fcNaN)
1710           newPart = significandParts()[0];
1711         freeSignificand();
1712         significand.part = newPart;
1713     }
1714   }
1715
1716   if(category == fcNormal) {
1717     /* Re-interpret our bit-pattern.  */
1718     exponent += toSemantics.precision - semantics->precision;
1719     semantics = &toSemantics;
1720     fs = normalize(rounding_mode, lostFraction);
1721   } else if (category == fcNaN) {
1722     int shift = toSemantics.precision - semantics->precision;
1723     // Do this now so significandParts gets the right answer
1724     const fltSemantics *oldSemantics = semantics;
1725     semantics = &toSemantics;
1726     fs = opOK;
1727     // No normalization here, just truncate
1728     if (shift>0)
1729       APInt::tcShiftLeft(significandParts(), newPartCount, shift);
1730     else if (shift < 0) {
1731       unsigned ushift = -shift;
1732       // We mark this as Inexact if we are losing information.  This happens
1733       // if are shifting out something other than 0s, or if the x87 long
1734       // double input did not have its integer bit set (pseudo-NaN), or if the
1735       // x87 long double input did not have its QNan bit set (because the x87
1736       // hardware sets this bit when converting a lower-precision NaN to
1737       // x87 long double).
1738       if (APInt::tcLSB(significandParts(), newPartCount) < ushift)
1739         fs = opInexact;
1740       if (oldSemantics == &APFloat::x87DoubleExtended &&
1741           (!(*significandParts() & 0x8000000000000000ULL) ||
1742            !(*significandParts() & 0x4000000000000000ULL)))
1743         fs = opInexact;
1744       APInt::tcShiftRight(significandParts(), newPartCount, ushift);
1745     }
1746     // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
1747     // does not give you back the same bits.  This is dubious, and we
1748     // don't currently do it.  You're really supposed to get
1749     // an invalid operation signal at runtime, but nobody does that.
1750   } else {
1751     semantics = &toSemantics;
1752     fs = opOK;
1753   }
1754
1755   return fs;
1756 }
1757
1758 /* Convert a floating point number to an integer according to the
1759    rounding mode.  If the rounded integer value is out of range this
1760    returns an invalid operation exception and the contents of the
1761    destination parts are unspecified.  If the rounded value is in
1762    range but the floating point number is not the exact integer, the C
1763    standard doesn't require an inexact exception to be raised.  IEEE
1764    854 does require it so we do that.
1765
1766    Note that for conversions to integer type the C standard requires
1767    round-to-zero to always be used.  */
1768 APFloat::opStatus
1769 APFloat::convertToSignExtendedInteger(integerPart *parts, unsigned int width,
1770                                       bool isSigned,
1771                                       roundingMode rounding_mode) const
1772 {
1773   lostFraction lost_fraction;
1774   const integerPart *src;
1775   unsigned int dstPartsCount, truncatedBits;
1776
1777   assertArithmeticOK(*semantics);
1778
1779   /* Handle the three special cases first.  */
1780   if(category == fcInfinity || category == fcNaN)
1781     return opInvalidOp;
1782
1783   dstPartsCount = partCountForBits(width);
1784
1785   if(category == fcZero) {
1786     APInt::tcSet(parts, 0, dstPartsCount);
1787     return opOK;
1788   }
1789
1790   src = significandParts();
1791
1792   /* Step 1: place our absolute value, with any fraction truncated, in
1793      the destination.  */
1794   if (exponent < 0) {
1795     /* Our absolute value is less than one; truncate everything.  */
1796     APInt::tcSet(parts, 0, dstPartsCount);
1797     truncatedBits = semantics->precision;
1798   } else {
1799     /* We want the most significant (exponent + 1) bits; the rest are
1800        truncated.  */
1801     unsigned int bits = exponent + 1U;
1802
1803     /* Hopelessly large in magnitude?  */
1804     if (bits > width)
1805       return opInvalidOp;
1806
1807     if (bits < semantics->precision) {
1808       /* We truncate (semantics->precision - bits) bits.  */
1809       truncatedBits = semantics->precision - bits;
1810       APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits);
1811     } else {
1812       /* We want at least as many bits as are available.  */
1813       APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0);
1814       APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision);
1815       truncatedBits = 0;
1816     }
1817   }
1818
1819   /* Step 2: work out any lost fraction, and increment the absolute
1820      value if we would round away from zero.  */
1821   if (truncatedBits) {
1822     lost_fraction = lostFractionThroughTruncation(src, partCount(),
1823                                                   truncatedBits);
1824     if (lost_fraction != lfExactlyZero
1825         && roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) {
1826       if (APInt::tcIncrement(parts, dstPartsCount))
1827         return opInvalidOp;     /* Overflow.  */
1828     }
1829   } else {
1830     lost_fraction = lfExactlyZero;
1831   }
1832
1833   /* Step 3: check if we fit in the destination.  */
1834   unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1;
1835
1836   if (sign) {
1837     if (!isSigned) {
1838       /* Negative numbers cannot be represented as unsigned.  */
1839       if (omsb != 0)
1840         return opInvalidOp;
1841     } else {
1842       /* It takes omsb bits to represent the unsigned integer value.
1843          We lose a bit for the sign, but care is needed as the
1844          maximally negative integer is a special case.  */
1845       if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb)
1846         return opInvalidOp;
1847
1848       /* This case can happen because of rounding.  */
1849       if (omsb > width)
1850         return opInvalidOp;
1851     }
1852
1853     APInt::tcNegate (parts, dstPartsCount);
1854   } else {
1855     if (omsb >= width + !isSigned)
1856       return opInvalidOp;
1857   }
1858
1859   if (lost_fraction == lfExactlyZero)
1860     return opOK;
1861   else
1862     return opInexact;
1863 }
1864
1865 /* Same as convertToSignExtendedInteger, except we provide
1866    deterministic values in case of an invalid operation exception,
1867    namely zero for NaNs and the minimal or maximal value respectively
1868    for underflow or overflow.  */
1869 APFloat::opStatus
1870 APFloat::convertToInteger(integerPart *parts, unsigned int width,
1871                           bool isSigned,
1872                           roundingMode rounding_mode) const
1873 {
1874   opStatus fs;
1875
1876   fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode);
1877
1878   if (fs == opInvalidOp) {
1879     unsigned int bits, dstPartsCount;
1880
1881     dstPartsCount = partCountForBits(width);
1882
1883     if (category == fcNaN)
1884       bits = 0;
1885     else if (sign)
1886       bits = isSigned;
1887     else
1888       bits = width - isSigned;
1889
1890     APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits);
1891     if (sign && isSigned)
1892       APInt::tcShiftLeft(parts, dstPartsCount, width - 1);
1893   }
1894
1895   return fs;
1896 }
1897
1898 /* Convert an unsigned integer SRC to a floating point number,
1899    rounding according to ROUNDING_MODE.  The sign of the floating
1900    point number is not modified.  */
1901 APFloat::opStatus
1902 APFloat::convertFromUnsignedParts(const integerPart *src,
1903                                   unsigned int srcCount,
1904                                   roundingMode rounding_mode)
1905 {
1906   unsigned int omsb, precision, dstCount;
1907   integerPart *dst;
1908   lostFraction lost_fraction;
1909
1910   assertArithmeticOK(*semantics);
1911   category = fcNormal;
1912   omsb = APInt::tcMSB(src, srcCount) + 1;
1913   dst = significandParts();
1914   dstCount = partCount();
1915   precision = semantics->precision;
1916
1917   /* We want the most significant PRECISON bits of SRC.  There may not
1918      be that many; extract what we can.  */
1919   if (precision <= omsb) {
1920     exponent = omsb - 1;
1921     lost_fraction = lostFractionThroughTruncation(src, srcCount,
1922                                                   omsb - precision);
1923     APInt::tcExtract(dst, dstCount, src, precision, omsb - precision);
1924   } else {
1925     exponent = precision - 1;
1926     lost_fraction = lfExactlyZero;
1927     APInt::tcExtract(dst, dstCount, src, omsb, 0);
1928   }
1929
1930   return normalize(rounding_mode, lost_fraction);
1931 }
1932
1933 APFloat::opStatus
1934 APFloat::convertFromAPInt(const APInt &Val,
1935                           bool isSigned,
1936                           roundingMode rounding_mode)
1937 {
1938   unsigned int partCount = Val.getNumWords();
1939   APInt api = Val;
1940
1941   sign = false;
1942   if (isSigned && api.isNegative()) {
1943     sign = true;
1944     api = -api;
1945   }
1946
1947   return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
1948 }
1949
1950 /* Convert a two's complement integer SRC to a floating point number,
1951    rounding according to ROUNDING_MODE.  ISSIGNED is true if the
1952    integer is signed, in which case it must be sign-extended.  */
1953 APFloat::opStatus
1954 APFloat::convertFromSignExtendedInteger(const integerPart *src,
1955                                         unsigned int srcCount,
1956                                         bool isSigned,
1957                                         roundingMode rounding_mode)
1958 {
1959   opStatus status;
1960
1961   assertArithmeticOK(*semantics);
1962   if (isSigned
1963       && APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) {
1964     integerPart *copy;
1965
1966     /* If we're signed and negative negate a copy.  */
1967     sign = true;
1968     copy = new integerPart[srcCount];
1969     APInt::tcAssign(copy, src, srcCount);
1970     APInt::tcNegate(copy, srcCount);
1971     status = convertFromUnsignedParts(copy, srcCount, rounding_mode);
1972     delete [] copy;
1973   } else {
1974     sign = false;
1975     status = convertFromUnsignedParts(src, srcCount, rounding_mode);
1976   }
1977
1978   return status;
1979 }
1980
1981 /* FIXME: should this just take a const APInt reference?  */
1982 APFloat::opStatus
1983 APFloat::convertFromZeroExtendedInteger(const integerPart *parts,
1984                                         unsigned int width, bool isSigned,
1985                                         roundingMode rounding_mode)
1986 {
1987   unsigned int partCount = partCountForBits(width);
1988   APInt api = APInt(width, partCount, parts);
1989
1990   sign = false;
1991   if(isSigned && APInt::tcExtractBit(parts, width - 1)) {
1992     sign = true;
1993     api = -api;
1994   }
1995
1996   return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
1997 }
1998
1999 APFloat::opStatus
2000 APFloat::convertFromHexadecimalString(const char *p,
2001                                       roundingMode rounding_mode)
2002 {
2003   lostFraction lost_fraction;
2004   integerPart *significand;
2005   unsigned int bitPos, partsCount;
2006   const char *dot, *firstSignificantDigit;
2007
2008   zeroSignificand();
2009   exponent = 0;
2010   category = fcNormal;
2011
2012   significand = significandParts();
2013   partsCount = partCount();
2014   bitPos = partsCount * integerPartWidth;
2015
2016   /* Skip leading zeroes and any (hexa)decimal point.  */
2017   p = skipLeadingZeroesAndAnyDot(p, &dot);
2018   firstSignificantDigit = p;
2019
2020   for(;;) {
2021     integerPart hex_value;
2022
2023     if(*p == '.') {
2024       assert(dot == 0);
2025       dot = p++;
2026     }
2027
2028     hex_value = hexDigitValue(*p);
2029     if(hex_value == -1U) {
2030       lost_fraction = lfExactlyZero;
2031       break;
2032     }
2033
2034     p++;
2035
2036     /* Store the number whilst 4-bit nibbles remain.  */
2037     if(bitPos) {
2038       bitPos -= 4;
2039       hex_value <<= bitPos % integerPartWidth;
2040       significand[bitPos / integerPartWidth] |= hex_value;
2041     } else {
2042       lost_fraction = trailingHexadecimalFraction(p, hex_value);
2043       while(hexDigitValue(*p) != -1U)
2044         p++;
2045       break;
2046     }
2047   }
2048
2049   /* Hex floats require an exponent but not a hexadecimal point.  */
2050   assert(*p == 'p' || *p == 'P');
2051
2052   /* Ignore the exponent if we are zero.  */
2053   if(p != firstSignificantDigit) {
2054     int expAdjustment;
2055
2056     /* Implicit hexadecimal point?  */
2057     if(!dot)
2058       dot = p;
2059
2060     /* Calculate the exponent adjustment implicit in the number of
2061        significant digits.  */
2062     expAdjustment = static_cast<int>(dot - firstSignificantDigit);
2063     if(expAdjustment < 0)
2064       expAdjustment++;
2065     expAdjustment = expAdjustment * 4 - 1;
2066
2067     /* Adjust for writing the significand starting at the most
2068        significant nibble.  */
2069     expAdjustment += semantics->precision;
2070     expAdjustment -= partsCount * integerPartWidth;
2071
2072     /* Adjust for the given exponent.  */
2073     exponent = totalExponent(p, expAdjustment);
2074   }
2075
2076   return normalize(rounding_mode, lost_fraction);
2077 }
2078
2079 APFloat::opStatus
2080 APFloat::roundSignificandWithExponent(const integerPart *decSigParts,
2081                                       unsigned sigPartCount, int exp,
2082                                       roundingMode rounding_mode)
2083 {
2084   unsigned int parts, pow5PartCount;
2085   fltSemantics calcSemantics = { 32767, -32767, 0, true };
2086   integerPart pow5Parts[maxPowerOfFiveParts];
2087   bool isNearest;
2088
2089   isNearest = (rounding_mode == rmNearestTiesToEven
2090                || rounding_mode == rmNearestTiesToAway);
2091
2092   parts = partCountForBits(semantics->precision + 11);
2093
2094   /* Calculate pow(5, abs(exp)).  */
2095   pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp);
2096
2097   for (;; parts *= 2) {
2098     opStatus sigStatus, powStatus;
2099     unsigned int excessPrecision, truncatedBits;
2100
2101     calcSemantics.precision = parts * integerPartWidth - 1;
2102     excessPrecision = calcSemantics.precision - semantics->precision;
2103     truncatedBits = excessPrecision;
2104
2105     APFloat decSig(calcSemantics, fcZero, sign);
2106     APFloat pow5(calcSemantics, fcZero, false);
2107
2108     sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount,
2109                                                 rmNearestTiesToEven);
2110     powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount,
2111                                               rmNearestTiesToEven);
2112     /* Add exp, as 10^n = 5^n * 2^n.  */
2113     decSig.exponent += exp;
2114
2115     lostFraction calcLostFraction;
2116     integerPart HUerr, HUdistance;
2117     unsigned int powHUerr;
2118
2119     if (exp >= 0) {
2120       /* multiplySignificand leaves the precision-th bit set to 1.  */
2121       calcLostFraction = decSig.multiplySignificand(pow5, NULL);
2122       powHUerr = powStatus != opOK;
2123     } else {
2124       calcLostFraction = decSig.divideSignificand(pow5);
2125       /* Denormal numbers have less precision.  */
2126       if (decSig.exponent < semantics->minExponent) {
2127         excessPrecision += (semantics->minExponent - decSig.exponent);
2128         truncatedBits = excessPrecision;
2129         if (excessPrecision > calcSemantics.precision)
2130           excessPrecision = calcSemantics.precision;
2131       }
2132       /* Extra half-ulp lost in reciprocal of exponent.  */
2133       powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2;
2134     }
2135
2136     /* Both multiplySignificand and divideSignificand return the
2137        result with the integer bit set.  */
2138     assert (APInt::tcExtractBit
2139             (decSig.significandParts(), calcSemantics.precision - 1) == 1);
2140
2141     HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK,
2142                        powHUerr);
2143     HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(),
2144                                       excessPrecision, isNearest);
2145
2146     /* Are we guaranteed to round correctly if we truncate?  */
2147     if (HUdistance >= HUerr) {
2148       APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(),
2149                        calcSemantics.precision - excessPrecision,
2150                        excessPrecision);
2151       /* Take the exponent of decSig.  If we tcExtract-ed less bits
2152          above we must adjust our exponent to compensate for the
2153          implicit right shift.  */
2154       exponent = (decSig.exponent + semantics->precision
2155                   - (calcSemantics.precision - excessPrecision));
2156       calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(),
2157                                                        decSig.partCount(),
2158                                                        truncatedBits);
2159       return normalize(rounding_mode, calcLostFraction);
2160     }
2161   }
2162 }
2163
2164 APFloat::opStatus
2165 APFloat::convertFromDecimalString(const char *p, roundingMode rounding_mode)
2166 {
2167   decimalInfo D;
2168   opStatus fs;
2169
2170   /* Scan the text.  */
2171   interpretDecimal(p, &D);
2172
2173   /* Handle the quick cases.  First the case of no significant digits,
2174      i.e. zero, and then exponents that are obviously too large or too
2175      small.  Writing L for log 10 / log 2, a number d.ddddd*10^exp
2176      definitely overflows if
2177
2178            (exp - 1) * L >= maxExponent
2179
2180      and definitely underflows to zero where
2181
2182            (exp + 1) * L <= minExponent - precision
2183
2184      With integer arithmetic the tightest bounds for L are
2185
2186            93/28 < L < 196/59            [ numerator <= 256 ]
2187            42039/12655 < L < 28738/8651  [ numerator <= 65536 ]
2188   */
2189
2190   if (decDigitValue(*D.firstSigDigit) >= 10U) {
2191     category = fcZero;
2192     fs = opOK;
2193   } else if ((D.normalizedExponent + 1) * 28738
2194              <= 8651 * (semantics->minExponent - (int) semantics->precision)) {
2195     /* Underflow to zero and round.  */
2196     zeroSignificand();
2197     fs = normalize(rounding_mode, lfLessThanHalf);
2198   } else if ((D.normalizedExponent - 1) * 42039
2199              >= 12655 * semantics->maxExponent) {
2200     /* Overflow and round.  */
2201     fs = handleOverflow(rounding_mode);
2202   } else {
2203     integerPart *decSignificand;
2204     unsigned int partCount;
2205
2206     /* A tight upper bound on number of bits required to hold an
2207        N-digit decimal integer is N * 196 / 59.  Allocate enough space
2208        to hold the full significand, and an extra part required by
2209        tcMultiplyPart.  */
2210     partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1;
2211     partCount = partCountForBits(1 + 196 * partCount / 59);
2212     decSignificand = new integerPart[partCount + 1];
2213     partCount = 0;
2214
2215     /* Convert to binary efficiently - we do almost all multiplication
2216        in an integerPart.  When this would overflow do we do a single
2217        bignum multiplication, and then revert again to multiplication
2218        in an integerPart.  */
2219     do {
2220       integerPart decValue, val, multiplier;
2221
2222       val = 0;
2223       multiplier = 1;
2224
2225       do {
2226         if (*p == '.')
2227           p++;
2228
2229         decValue = decDigitValue(*p++);
2230         multiplier *= 10;
2231         val = val * 10 + decValue;
2232         /* The maximum number that can be multiplied by ten with any
2233            digit added without overflowing an integerPart.  */
2234       } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10);
2235
2236       /* Multiply out the current part.  */
2237       APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val,
2238                             partCount, partCount + 1, false);
2239
2240       /* If we used another part (likely but not guaranteed), increase
2241          the count.  */
2242       if (decSignificand[partCount])
2243         partCount++;
2244     } while (p <= D.lastSigDigit);
2245
2246     category = fcNormal;
2247     fs = roundSignificandWithExponent(decSignificand, partCount,
2248                                       D.exponent, rounding_mode);
2249
2250     delete [] decSignificand;
2251   }
2252
2253   return fs;
2254 }
2255
2256 APFloat::opStatus
2257 APFloat::convertFromString(const char *p, roundingMode rounding_mode)
2258 {
2259   assertArithmeticOK(*semantics);
2260
2261   /* Handle a leading minus sign.  */
2262   if(*p == '-')
2263     sign = 1, p++;
2264   else
2265     sign = 0;
2266
2267   if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
2268     return convertFromHexadecimalString(p + 2, rounding_mode);
2269   else
2270     return convertFromDecimalString(p, rounding_mode);
2271 }
2272
2273 /* Write out a hexadecimal representation of the floating point value
2274    to DST, which must be of sufficient size, in the C99 form
2275    [-]0xh.hhhhp[+-]d.  Return the number of characters written,
2276    excluding the terminating NUL.
2277
2278    If UPPERCASE, the output is in upper case, otherwise in lower case.
2279
2280    HEXDIGITS digits appear altogether, rounding the value if
2281    necessary.  If HEXDIGITS is 0, the minimal precision to display the
2282    number precisely is used instead.  If nothing would appear after
2283    the decimal point it is suppressed.
2284
2285    The decimal exponent is always printed and has at least one digit.
2286    Zero values display an exponent of zero.  Infinities and NaNs
2287    appear as "infinity" or "nan" respectively.
2288
2289    The above rules are as specified by C99.  There is ambiguity about
2290    what the leading hexadecimal digit should be.  This implementation
2291    uses whatever is necessary so that the exponent is displayed as
2292    stored.  This implies the exponent will fall within the IEEE format
2293    range, and the leading hexadecimal digit will be 0 (for denormals),
2294    1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with
2295    any other digits zero).
2296 */
2297 unsigned int
2298 APFloat::convertToHexString(char *dst, unsigned int hexDigits,
2299                             bool upperCase, roundingMode rounding_mode) const
2300 {
2301   char *p;
2302
2303   assertArithmeticOK(*semantics);
2304
2305   p = dst;
2306   if (sign)
2307     *dst++ = '-';
2308
2309   switch (category) {
2310   case fcInfinity:
2311     memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1);
2312     dst += sizeof infinityL - 1;
2313     break;
2314
2315   case fcNaN:
2316     memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1);
2317     dst += sizeof NaNU - 1;
2318     break;
2319
2320   case fcZero:
2321     *dst++ = '0';
2322     *dst++ = upperCase ? 'X': 'x';
2323     *dst++ = '0';
2324     if (hexDigits > 1) {
2325       *dst++ = '.';
2326       memset (dst, '0', hexDigits - 1);
2327       dst += hexDigits - 1;
2328     }
2329     *dst++ = upperCase ? 'P': 'p';
2330     *dst++ = '0';
2331     break;
2332
2333   case fcNormal:
2334     dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode);
2335     break;
2336   }
2337
2338   *dst = 0;
2339
2340   return static_cast<unsigned int>(dst - p);
2341 }
2342
2343 /* Does the hard work of outputting the correctly rounded hexadecimal
2344    form of a normal floating point number with the specified number of
2345    hexadecimal digits.  If HEXDIGITS is zero the minimum number of
2346    digits necessary to print the value precisely is output.  */
2347 char *
2348 APFloat::convertNormalToHexString(char *dst, unsigned int hexDigits,
2349                                   bool upperCase,
2350                                   roundingMode rounding_mode) const
2351 {
2352   unsigned int count, valueBits, shift, partsCount, outputDigits;
2353   const char *hexDigitChars;
2354   const integerPart *significand;
2355   char *p;
2356   bool roundUp;
2357
2358   *dst++ = '0';
2359   *dst++ = upperCase ? 'X': 'x';
2360
2361   roundUp = false;
2362   hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower;
2363
2364   significand = significandParts();
2365   partsCount = partCount();
2366
2367   /* +3 because the first digit only uses the single integer bit, so
2368      we have 3 virtual zero most-significant-bits.  */
2369   valueBits = semantics->precision + 3;
2370   shift = integerPartWidth - valueBits % integerPartWidth;
2371
2372   /* The natural number of digits required ignoring trailing
2373      insignificant zeroes.  */
2374   outputDigits = (valueBits - significandLSB () + 3) / 4;
2375
2376   /* hexDigits of zero means use the required number for the
2377      precision.  Otherwise, see if we are truncating.  If we are,
2378      find out if we need to round away from zero.  */
2379   if (hexDigits) {
2380     if (hexDigits < outputDigits) {
2381       /* We are dropping non-zero bits, so need to check how to round.
2382          "bits" is the number of dropped bits.  */
2383       unsigned int bits;
2384       lostFraction fraction;
2385
2386       bits = valueBits - hexDigits * 4;
2387       fraction = lostFractionThroughTruncation (significand, partsCount, bits);
2388       roundUp = roundAwayFromZero(rounding_mode, fraction, bits);
2389     }
2390     outputDigits = hexDigits;
2391   }
2392
2393   /* Write the digits consecutively, and start writing in the location
2394      of the hexadecimal point.  We move the most significant digit
2395      left and add the hexadecimal point later.  */
2396   p = ++dst;
2397
2398   count = (valueBits + integerPartWidth - 1) / integerPartWidth;
2399
2400   while (outputDigits && count) {
2401     integerPart part;
2402
2403     /* Put the most significant integerPartWidth bits in "part".  */
2404     if (--count == partsCount)
2405       part = 0;  /* An imaginary higher zero part.  */
2406     else
2407       part = significand[count] << shift;
2408
2409     if (count && shift)
2410       part |= significand[count - 1] >> (integerPartWidth - shift);
2411
2412     /* Convert as much of "part" to hexdigits as we can.  */
2413     unsigned int curDigits = integerPartWidth / 4;
2414
2415     if (curDigits > outputDigits)
2416       curDigits = outputDigits;
2417     dst += partAsHex (dst, part, curDigits, hexDigitChars);
2418     outputDigits -= curDigits;
2419   }
2420
2421   if (roundUp) {
2422     char *q = dst;
2423
2424     /* Note that hexDigitChars has a trailing '0'.  */
2425     do {
2426       q--;
2427       *q = hexDigitChars[hexDigitValue (*q) + 1];
2428     } while (*q == '0');
2429     assert (q >= p);
2430   } else {
2431     /* Add trailing zeroes.  */
2432     memset (dst, '0', outputDigits);
2433     dst += outputDigits;
2434   }
2435
2436   /* Move the most significant digit to before the point, and if there
2437      is something after the decimal point add it.  This must come
2438      after rounding above.  */
2439   p[-1] = p[0];
2440   if (dst -1 == p)
2441     dst--;
2442   else
2443     p[0] = '.';
2444
2445   /* Finally output the exponent.  */
2446   *dst++ = upperCase ? 'P': 'p';
2447
2448   return writeSignedDecimal (dst, exponent);
2449 }
2450
2451 // For good performance it is desirable for different APFloats
2452 // to produce different integers.
2453 uint32_t
2454 APFloat::getHashValue() const
2455 {
2456   if (category==fcZero) return sign<<8 | semantics->precision ;
2457   else if (category==fcInfinity) return sign<<9 | semantics->precision;
2458   else if (category==fcNaN) return 1<<10 | semantics->precision;
2459   else {
2460     uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
2461     const integerPart* p = significandParts();
2462     for (int i=partCount(); i>0; i--, p++)
2463       hash ^= ((uint32_t)*p) ^ (uint32_t)((*p)>>32);
2464     return hash;
2465   }
2466 }
2467
2468 // Conversion from APFloat to/from host float/double.  It may eventually be
2469 // possible to eliminate these and have everybody deal with APFloats, but that
2470 // will take a while.  This approach will not easily extend to long double.
2471 // Current implementation requires integerPartWidth==64, which is correct at
2472 // the moment but could be made more general.
2473
2474 // Denormals have exponent minExponent in APFloat, but minExponent-1 in
2475 // the actual IEEE respresentations.  We compensate for that here.
2476
2477 APInt
2478 APFloat::convertF80LongDoubleAPFloatToAPInt() const
2479 {
2480   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended);
2481   assert (partCount()==2);
2482
2483   uint64_t myexponent, mysignificand;
2484
2485   if (category==fcNormal) {
2486     myexponent = exponent+16383; //bias
2487     mysignificand = significandParts()[0];
2488     if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
2489       myexponent = 0;   // denormal
2490   } else if (category==fcZero) {
2491     myexponent = 0;
2492     mysignificand = 0;
2493   } else if (category==fcInfinity) {
2494     myexponent = 0x7fff;
2495     mysignificand = 0x8000000000000000ULL;
2496   } else {
2497     assert(category == fcNaN && "Unknown category");
2498     myexponent = 0x7fff;
2499     mysignificand = significandParts()[0];
2500   }
2501
2502   uint64_t words[2];
2503   words[0] =  ((uint64_t)(sign & 1) << 63) |
2504               ((myexponent & 0x7fffLL) <<  48) |
2505               ((mysignificand >>16) & 0xffffffffffffLL);
2506   words[1] = mysignificand & 0xffff;
2507   return APInt(80, 2, words);
2508 }
2509
2510 APInt
2511 APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const
2512 {
2513   assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble);
2514   assert (partCount()==2);
2515
2516   uint64_t myexponent, mysignificand, myexponent2, mysignificand2;
2517
2518   if (category==fcNormal) {
2519     myexponent = exponent + 1023; //bias
2520     myexponent2 = exponent2 + 1023;
2521     mysignificand = significandParts()[0];
2522     mysignificand2 = significandParts()[1];
2523     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2524       myexponent = 0;   // denormal
2525     if (myexponent2==1 && !(mysignificand2 & 0x10000000000000LL))
2526       myexponent2 = 0;   // denormal
2527   } else if (category==fcZero) {
2528     myexponent = 0;
2529     mysignificand = 0;
2530     myexponent2 = 0;
2531     mysignificand2 = 0;
2532   } else if (category==fcInfinity) {
2533     myexponent = 0x7ff;
2534     myexponent2 = 0;
2535     mysignificand = 0;
2536     mysignificand2 = 0;
2537   } else {
2538     assert(category == fcNaN && "Unknown category");
2539     myexponent = 0x7ff;
2540     mysignificand = significandParts()[0];
2541     myexponent2 = exponent2;
2542     mysignificand2 = significandParts()[1];
2543   }
2544
2545   uint64_t words[2];
2546   words[0] =  ((uint64_t)(sign & 1) << 63) |
2547               ((myexponent & 0x7ff) <<  52) |
2548               (mysignificand & 0xfffffffffffffLL);
2549   words[1] =  ((uint64_t)(sign2 & 1) << 63) |
2550               ((myexponent2 & 0x7ff) <<  52) |
2551               (mysignificand2 & 0xfffffffffffffLL);
2552   return APInt(128, 2, words);
2553 }
2554
2555 APInt
2556 APFloat::convertDoubleAPFloatToAPInt() const
2557 {
2558   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2559   assert (partCount()==1);
2560
2561   uint64_t myexponent, mysignificand;
2562
2563   if (category==fcNormal) {
2564     myexponent = exponent+1023; //bias
2565     mysignificand = *significandParts();
2566     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2567       myexponent = 0;   // denormal
2568   } else if (category==fcZero) {
2569     myexponent = 0;
2570     mysignificand = 0;
2571   } else if (category==fcInfinity) {
2572     myexponent = 0x7ff;
2573     mysignificand = 0;
2574   } else {
2575     assert(category == fcNaN && "Unknown category!");
2576     myexponent = 0x7ff;
2577     mysignificand = *significandParts();
2578   }
2579
2580   return APInt(64, ((((uint64_t)(sign & 1) << 63) |
2581                      ((myexponent & 0x7ff) <<  52) |
2582                      (mysignificand & 0xfffffffffffffLL))));
2583 }
2584
2585 APInt
2586 APFloat::convertFloatAPFloatToAPInt() const
2587 {
2588   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2589   assert (partCount()==1);
2590
2591   uint32_t myexponent, mysignificand;
2592
2593   if (category==fcNormal) {
2594     myexponent = exponent+127; //bias
2595     mysignificand = (uint32_t)*significandParts();
2596     if (myexponent == 1 && !(mysignificand & 0x800000))
2597       myexponent = 0;   // denormal
2598   } else if (category==fcZero) {
2599     myexponent = 0;
2600     mysignificand = 0;
2601   } else if (category==fcInfinity) {
2602     myexponent = 0xff;
2603     mysignificand = 0;
2604   } else {
2605     assert(category == fcNaN && "Unknown category!");
2606     myexponent = 0xff;
2607     mysignificand = (uint32_t)*significandParts();
2608   }
2609
2610   return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
2611                     (mysignificand & 0x7fffff)));
2612 }
2613
2614 // This function creates an APInt that is just a bit map of the floating
2615 // point constant as it would appear in memory.  It is not a conversion,
2616 // and treating the result as a normal integer is unlikely to be useful.
2617
2618 APInt
2619 APFloat::convertToAPInt() const
2620 {
2621   if (semantics == (const llvm::fltSemantics*)&IEEEsingle)
2622     return convertFloatAPFloatToAPInt();
2623
2624   if (semantics == (const llvm::fltSemantics*)&IEEEdouble)
2625     return convertDoubleAPFloatToAPInt();
2626
2627   if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble)
2628     return convertPPCDoubleDoubleAPFloatToAPInt();
2629
2630   assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended &&
2631          "unknown format!");
2632   return convertF80LongDoubleAPFloatToAPInt();
2633 }
2634
2635 float
2636 APFloat::convertToFloat() const
2637 {
2638   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2639   APInt api = convertToAPInt();
2640   return api.bitsToFloat();
2641 }
2642
2643 double
2644 APFloat::convertToDouble() const
2645 {
2646   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2647   APInt api = convertToAPInt();
2648   return api.bitsToDouble();
2649 }
2650
2651 /// Integer bit is explicit in this format.  Intel hardware (387 and later)
2652 /// does not support these bit patterns:
2653 ///  exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
2654 ///  exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
2655 ///  exponent = 0, integer bit 1 ("pseudodenormal")
2656 ///  exponent!=0 nor all 1's, integer bit 0 ("unnormal")
2657 /// At the moment, the first two are treated as NaNs, the second two as Normal.
2658 void
2659 APFloat::initFromF80LongDoubleAPInt(const APInt &api)
2660 {
2661   assert(api.getBitWidth()==80);
2662   uint64_t i1 = api.getRawData()[0];
2663   uint64_t i2 = api.getRawData()[1];
2664   uint64_t myexponent = (i1 >> 48) & 0x7fff;
2665   uint64_t mysignificand = ((i1 << 16) &  0xffffffffffff0000ULL) |
2666                           (i2 & 0xffff);
2667
2668   initialize(&APFloat::x87DoubleExtended);
2669   assert(partCount()==2);
2670
2671   sign = static_cast<unsigned int>(i1>>63);
2672   if (myexponent==0 && mysignificand==0) {
2673     // exponent, significand meaningless
2674     category = fcZero;
2675   } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
2676     // exponent, significand meaningless
2677     category = fcInfinity;
2678   } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
2679     // exponent meaningless
2680     category = fcNaN;
2681     significandParts()[0] = mysignificand;
2682     significandParts()[1] = 0;
2683   } else {
2684     category = fcNormal;
2685     exponent = myexponent - 16383;
2686     significandParts()[0] = mysignificand;
2687     significandParts()[1] = 0;
2688     if (myexponent==0)          // denormal
2689       exponent = -16382;
2690   }
2691 }
2692
2693 void
2694 APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api)
2695 {
2696   assert(api.getBitWidth()==128);
2697   uint64_t i1 = api.getRawData()[0];
2698   uint64_t i2 = api.getRawData()[1];
2699   uint64_t myexponent = (i1 >> 52) & 0x7ff;
2700   uint64_t mysignificand = i1 & 0xfffffffffffffLL;
2701   uint64_t myexponent2 = (i2 >> 52) & 0x7ff;
2702   uint64_t mysignificand2 = i2 & 0xfffffffffffffLL;
2703
2704   initialize(&APFloat::PPCDoubleDouble);
2705   assert(partCount()==2);
2706
2707   sign = static_cast<unsigned int>(i1>>63);
2708   sign2 = static_cast<unsigned int>(i2>>63);
2709   if (myexponent==0 && mysignificand==0) {
2710     // exponent, significand meaningless
2711     // exponent2 and significand2 are required to be 0; we don't check
2712     category = fcZero;
2713   } else if (myexponent==0x7ff && mysignificand==0) {
2714     // exponent, significand meaningless
2715     // exponent2 and significand2 are required to be 0; we don't check
2716     category = fcInfinity;
2717   } else if (myexponent==0x7ff && mysignificand!=0) {
2718     // exponent meaningless.  So is the whole second word, but keep it
2719     // for determinism.
2720     category = fcNaN;
2721     exponent2 = myexponent2;
2722     significandParts()[0] = mysignificand;
2723     significandParts()[1] = mysignificand2;
2724   } else {
2725     category = fcNormal;
2726     // Note there is no category2; the second word is treated as if it is
2727     // fcNormal, although it might be something else considered by itself.
2728     exponent = myexponent - 1023;
2729     exponent2 = myexponent2 - 1023;
2730     significandParts()[0] = mysignificand;
2731     significandParts()[1] = mysignificand2;
2732     if (myexponent==0)          // denormal
2733       exponent = -1022;
2734     else
2735       significandParts()[0] |= 0x10000000000000LL;  // integer bit
2736     if (myexponent2==0)
2737       exponent2 = -1022;
2738     else
2739       significandParts()[1] |= 0x10000000000000LL;  // integer bit
2740   }
2741 }
2742
2743 void
2744 APFloat::initFromDoubleAPInt(const APInt &api)
2745 {
2746   assert(api.getBitWidth()==64);
2747   uint64_t i = *api.getRawData();
2748   uint64_t myexponent = (i >> 52) & 0x7ff;
2749   uint64_t mysignificand = i & 0xfffffffffffffLL;
2750
2751   initialize(&APFloat::IEEEdouble);
2752   assert(partCount()==1);
2753
2754   sign = static_cast<unsigned int>(i>>63);
2755   if (myexponent==0 && mysignificand==0) {
2756     // exponent, significand meaningless
2757     category = fcZero;
2758   } else if (myexponent==0x7ff && mysignificand==0) {
2759     // exponent, significand meaningless
2760     category = fcInfinity;
2761   } else if (myexponent==0x7ff && mysignificand!=0) {
2762     // exponent meaningless
2763     category = fcNaN;
2764     *significandParts() = mysignificand;
2765   } else {
2766     category = fcNormal;
2767     exponent = myexponent - 1023;
2768     *significandParts() = mysignificand;
2769     if (myexponent==0)          // denormal
2770       exponent = -1022;
2771     else
2772       *significandParts() |= 0x10000000000000LL;  // integer bit
2773   }
2774 }
2775
2776 void
2777 APFloat::initFromFloatAPInt(const APInt & api)
2778 {
2779   assert(api.getBitWidth()==32);
2780   uint32_t i = (uint32_t)*api.getRawData();
2781   uint32_t myexponent = (i >> 23) & 0xff;
2782   uint32_t mysignificand = i & 0x7fffff;
2783
2784   initialize(&APFloat::IEEEsingle);
2785   assert(partCount()==1);
2786
2787   sign = i >> 31;
2788   if (myexponent==0 && mysignificand==0) {
2789     // exponent, significand meaningless
2790     category = fcZero;
2791   } else if (myexponent==0xff && mysignificand==0) {
2792     // exponent, significand meaningless
2793     category = fcInfinity;
2794   } else if (myexponent==0xff && mysignificand!=0) {
2795     // sign, exponent, significand meaningless
2796     category = fcNaN;
2797     *significandParts() = mysignificand;
2798   } else {
2799     category = fcNormal;
2800     exponent = myexponent - 127;  //bias
2801     *significandParts() = mysignificand;
2802     if (myexponent==0)    // denormal
2803       exponent = -126;
2804     else
2805       *significandParts() |= 0x800000; // integer bit
2806   }
2807 }
2808
2809 /// Treat api as containing the bits of a floating point number.  Currently
2810 /// we infer the floating point type from the size of the APInt.  The
2811 /// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful
2812 /// when the size is anything else).
2813 void
2814 APFloat::initFromAPInt(const APInt& api, bool isIEEE)
2815 {
2816   if (api.getBitWidth() == 32)
2817     return initFromFloatAPInt(api);
2818   else if (api.getBitWidth()==64)
2819     return initFromDoubleAPInt(api);
2820   else if (api.getBitWidth()==80)
2821     return initFromF80LongDoubleAPInt(api);
2822   else if (api.getBitWidth()==128 && !isIEEE)
2823     return initFromPPCDoubleDoubleAPInt(api);
2824   else
2825     assert(0);
2826 }
2827
2828 APFloat::APFloat(const APInt& api, bool isIEEE)
2829 {
2830   initFromAPInt(api, isIEEE);
2831 }
2832
2833 APFloat::APFloat(float f)
2834 {
2835   APInt api = APInt(32, 0);
2836   initFromAPInt(api.floatToBits(f));
2837 }
2838
2839 APFloat::APFloat(double d)
2840 {
2841   APInt api = APInt(64, 0);
2842   initFromAPInt(api.doubleToBits(d));
2843 }