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