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