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