lib/Support/APFloat.cpp

   1 //===-- APFloat.cpp - Implement APFloat class -----------------------------===//
   2 //
   3 //                     The LLVM Compiler Infrastructure
   4 //
   5 // This file was developed by Neil Booth and is distributed under the
   6 // University of Illinois Open Source 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 <cassert>
  16 #include "llvm/ADT/APFloat.h"
  17 #include "llvm/Support/MathExtras.h"
  18
  19 using namespace llvm;
  20
  21 #define convolve(lhs, rhs) ((lhs) * 4 + (rhs))
  22
  23 /* Assumed in hexadecimal significand parsing.  */
  24 COMPILE_TIME_ASSERT(integerPartWidth % 4 == 0);
  25
  26 namespace llvm {
  27
  28   /* Represents floating point arithmetic semantics.  */
  29   struct fltSemantics {
  30     /* The largest E such that 2^E is representable; this matches the
  31        definition of IEEE 754.  */
  32     exponent_t maxExponent;
  33
  34     /* The smallest E such that 2^E is a normalized number; this
  35        matches the definition of IEEE 754.  */
  36     exponent_t minExponent;
  37
  38     /* Number of bits in the significand.  This includes the integer
  39        bit.  */
  40     unsigned char precision;
  41
  42     /* If the target format has an implicit integer bit.  */
  43     bool implicitIntegerBit;
  44   };
  45
  46   const fltSemantics APFloat::IEEEsingle = { 127, -126, 24, true };
  47   const fltSemantics APFloat::IEEEdouble = { 1023, -1022, 53, true };
  48   const fltSemantics APFloat::IEEEquad = { 16383, -16382, 113, true };
  49   const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64, false };
  50   const fltSemantics APFloat::Bogus = { 0, 0, 0, false };
  51 }
  52
  53 /* Put a bunch of private, handy routines in an anonymous namespace.  */
  54 namespace {
  55
  56   inline unsigned int
  57   partCountForBits(unsigned int bits)
  58   {
  59     return ((bits) + integerPartWidth - 1) / integerPartWidth;
  60   }
  61
  62   unsigned int
  63   digitValue(unsigned int c)
  64   {
  65     unsigned int r;
  66
  67     r = c - '0';
  68     if(r <= 9)
  69       return r;
  70
  71     return -1U;
  72   }
  73
  74   unsigned int
  75   hexDigitValue (unsigned int c)
  76   {
  77     unsigned int r;
  78
  79     r = c - '0';
  80     if(r <= 9)
  81       return r;
  82
  83     r = c - 'A';
  84     if(r <= 5)
  85       return r + 10;
  86
  87     r = c - 'a';
  88     if(r <= 5)
  89       return r + 10;
  90
  91     return -1U;
  92   }
  93
  94   /* This is ugly and needs cleaning up, but I don't immediately see
  95      how whilst remaining safe.  */
  96   static int
  97   totalExponent(const char *p, int exponentAdjustment)
  98   {
  99     integerPart unsignedExponent;
 100     bool negative, overflow;
 101     long exponent;
 102
 103     /* Move past the exponent letter and sign to the digits.  */
 104     p++;
 105     negative = *p == '-';
 106     if(*p == '-' || *p == '+')
 107       p++;
 108
 109     unsignedExponent = 0;
 110     overflow = false;
 111     for(;;) {
 112       unsigned int value;
 113
 114       value = digitValue(*p);
 115       if(value == -1U)
 116         break;
 117
 118       p++;
 119       unsignedExponent = unsignedExponent * 10 + value;
 120       if(unsignedExponent > 65535)
 121         overflow = true;
 122     }
 123
 124     if(exponentAdjustment > 65535 || exponentAdjustment < -65536)
 125       overflow = true;
 126
 127     if(!overflow) {
 128       exponent = unsignedExponent;
 129       if(negative)
 130         exponent = -exponent;
 131       exponent += exponentAdjustment;
 132       if(exponent > 65535 || exponent < -65536)
 133         overflow = true;
 134     }
 135
 136     if(overflow)
 137       exponent = negative ? -65536: 65535;
 138
 139     return exponent;
 140   }
 141
 142   const char *
 143   skipLeadingZeroesAndAnyDot(const char *p, const char **dot)
 144   {
 145     *dot = 0;
 146     while(*p == '0')
 147       p++;
 148
 149     if(*p == '.') {
 150       *dot = p++;
 151       while(*p == '0')
 152         p++;
 153     }
 154
 155     return p;
 156   }
 157
 158   /* Return the trailing fraction of a hexadecimal number.
 159      DIGITVALUE is the first hex digit of the fraction, P points to
 160      the next digit.  */
 161   lostFraction
 162   trailingHexadecimalFraction(const char *p, unsigned int digitValue)
 163   {
 164     unsigned int hexDigit;
 165
 166     /* If the first trailing digit isn't 0 or 8 we can work out the
 167        fraction immediately.  */
 168     if(digitValue > 8)
 169       return lfMoreThanHalf;
 170     else if(digitValue < 8 && digitValue > 0)
 171       return lfLessThanHalf;
 172
 173     /* Otherwise we need to find the first non-zero digit.  */
 174     while(*p == '0')
 175       p++;
 176
 177     hexDigit = hexDigitValue(*p);
 178
 179     /* If we ran off the end it is exactly zero or one-half, otherwise
 180        a little more.  */
 181     if(hexDigit == -1U)
 182       return digitValue == 0 ? lfExactlyZero: lfExactlyHalf;
 183     else
 184       return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf;
 185   }
 186
 187   /* Return the fraction lost were a bignum truncated.  */
 188   lostFraction
 189   lostFractionThroughTruncation(integerPart *parts,
 190                                 unsigned int partCount,
 191                                 unsigned int bits)
 192   {
 193     unsigned int lsb;
 194
 195     lsb = APInt::tcLSB(parts, partCount);
 196
 197     /* Note this is guaranteed true if bits == 0, or LSB == -1U.  */
 198     if(bits <= lsb)
 199       return lfExactlyZero;
 200     if(bits == lsb + 1)
 201       return lfExactlyHalf;
 202     if(bits <= partCount * integerPartWidth
 203        && APInt::tcExtractBit(parts, bits - 1))
 204       return lfMoreThanHalf;
 205
 206     return lfLessThanHalf;
 207   }
 208
 209   /* Shift DST right BITS bits noting lost fraction.  */
 210   lostFraction
 211   shiftRight(integerPart *dst, unsigned int parts, unsigned int bits)
 212   {
 213     lostFraction lost_fraction;
 214
 215     lost_fraction = lostFractionThroughTruncation(dst, parts, bits);
 216
 217     APInt::tcShiftRight(dst, parts, bits);
 218
 219     return lost_fraction;
 220   }
 221 }
 222
 223 /* Constructors.  */
 224 void
 225 APFloat::initialize(const fltSemantics *ourSemantics)
 226 {
 227   unsigned int count;
 228
 229   semantics = ourSemantics;
 230   count = partCount();
 231   if(count > 1)
 232     significand.parts = new integerPart[count];
 233 }
 234
 235 void
 236 APFloat::freeSignificand()
 237 {
 238   if(partCount() > 1)
 239     delete [] significand.parts;
 240 }
 241
 242 void
 243 APFloat::assign(const APFloat &rhs)
 244 {
 245   assert(semantics == rhs.semantics);
 246
 247   sign = rhs.sign;
 248   category = rhs.category;
 249   exponent = rhs.exponent;
 250   if(category == fcNormal)
 251     copySignificand(rhs);
 252 }
 253
 254 void
 255 APFloat::copySignificand(const APFloat &rhs)
 256 {
 257   assert(category == fcNormal);
 258   assert(rhs.partCount() >= partCount());
 259
 260   APInt::tcAssign(significandParts(), rhs.significandParts(),
 261                   partCount());
 262 }
 263
 264 APFloat &
 265 APFloat::operator=(const APFloat &rhs)
 266 {
 267   if(this != &rhs) {
 268     if(semantics != rhs.semantics) {
 269       freeSignificand();
 270       initialize(rhs.semantics);
 271     }
 272     assign(rhs);
 273   }
 274
 275   return *this;
 276 }
 277
 278 bool
 279 APFloat::operator==(const APFloat &rhs) const {
 280   if (this == &rhs)
 281     return true;
 282   if (semantics != rhs.semantics ||
 283       category != rhs.category)
 284     return false;
 285   if (category==fcQNaN)
 286     return true;
 287   else if (category==fcZero || category==fcInfinity)
 288     return sign==rhs.sign;
 289   else {
 290     if (sign!=rhs.sign || exponent!=rhs.exponent)
 291       return false;
 292     int i= partCount();
 293     const integerPart* p=significandParts();
 294     const integerPart* q=rhs.significandParts();
 295     for (; i>0; i--, p++, q++) {
 296       if (*p != *q)
 297         return false;
 298     }
 299     return true;
 300   }
 301 }
 302
 303 APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value)
 304 {
 305   initialize(&ourSemantics);
 306   sign = 0;
 307   zeroSignificand();
 308   exponent = ourSemantics.precision - 1;
 309   significandParts()[0] = value;
 310   normalize(rmNearestTiesToEven, lfExactlyZero);
 311 }
 312
 313 APFloat::APFloat(const fltSemantics &ourSemantics,
 314                  fltCategory ourCategory, bool negative)
 315 {
 316   initialize(&ourSemantics);
 317   category = ourCategory;
 318   sign = negative;
 319   if(category == fcNormal)
 320     category = fcZero;
 321 }
 322
 323 APFloat::APFloat(const fltSemantics &ourSemantics, const char *text)
 324 {
 325   initialize(&ourSemantics);
 326   convertFromString(text, rmNearestTiesToEven);
 327 }
 328
 329 APFloat::APFloat(const APFloat &rhs)
 330 {
 331   initialize(rhs.semantics);
 332   assign(rhs);
 333 }
 334
 335 APFloat::~APFloat()
 336 {
 337   freeSignificand();
 338 }
 339
 340 unsigned int
 341 APFloat::partCount() const
 342 {
 343   return partCountForBits(semantics->precision + 1);
 344 }
 345
 346 unsigned int
 347 APFloat::semanticsPrecision(const fltSemantics &semantics)
 348 {
 349   return semantics.precision;
 350 }
 351
 352 const integerPart *
 353 APFloat::significandParts() const
 354 {
 355   return const_cast<APFloat *>(this)->significandParts();
 356 }
 357
 358 integerPart *
 359 APFloat::significandParts()
 360 {
 361   assert(category == fcNormal);
 362
 363   if(partCount() > 1)
 364     return significand.parts;
 365   else
 366     return &significand.part;
 367 }
 368
 369 /* Combine the effect of two lost fractions.  */
 370 lostFraction
 371 APFloat::combineLostFractions(lostFraction moreSignificant,
 372                               lostFraction lessSignificant)
 373 {
 374   if(lessSignificant != lfExactlyZero) {
 375     if(moreSignificant == lfExactlyZero)
 376       moreSignificant = lfLessThanHalf;
 377     else if(moreSignificant == lfExactlyHalf)
 378       moreSignificant = lfMoreThanHalf;
 379   }
 380
 381   return moreSignificant;
 382 }
 383
 384 void
 385 APFloat::zeroSignificand()
 386 {
 387   category = fcNormal;
 388   APInt::tcSet(significandParts(), 0, partCount());
 389 }
 390
 391 /* Increment an fcNormal floating point number's significand.  */
 392 void
 393 APFloat::incrementSignificand()
 394 {
 395   integerPart carry;
 396
 397   carry = APInt::tcIncrement(significandParts(), partCount());
 398
 399   /* Our callers should never cause us to overflow.  */
 400   assert(carry == 0);
 401 }
 402
 403 /* Add the significand of the RHS.  Returns the carry flag.  */
 404 integerPart
 405 APFloat::addSignificand(const APFloat &rhs)
 406 {
 407   integerPart *parts;
 408
 409   parts = significandParts();
 410
 411   assert(semantics == rhs.semantics);
 412   assert(exponent == rhs.exponent);
 413
 414   return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
 415 }
 416
 417 /* Subtract the significand of the RHS with a borrow flag.  Returns
 418    the borrow flag.  */
 419 integerPart
 420 APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
 421 {
 422   integerPart *parts;
 423
 424   parts = significandParts();
 425
 426   assert(semantics == rhs.semantics);
 427   assert(exponent == rhs.exponent);
 428
 429   return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
 430                            partCount());
 431 }
 432
 433 /* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
 434    on to the full-precision result of the multiplication.  Returns the
 435    lost fraction.  */
 436 lostFraction
 437 APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
 438 {
 439   unsigned int omsb;    // One, not zero, based MSB.
 440   unsigned int partsCount, newPartsCount, precision;
 441   integerPart *lhsSignificand;
 442   integerPart scratch[4];
 443   integerPart *fullSignificand;
 444   lostFraction lost_fraction;
 445
 446   assert(semantics == rhs.semantics);
 447
 448   precision = semantics->precision;
 449   newPartsCount = partCountForBits(precision * 2);
 450
 451   if(newPartsCount > 4)
 452     fullSignificand = new integerPart[newPartsCount];
 453   else
 454     fullSignificand = scratch;
 455
 456   lhsSignificand = significandParts();
 457   partsCount = partCount();
 458
 459   APInt::tcFullMultiply(fullSignificand, lhsSignificand,
 460                         rhs.significandParts(), partsCount);
 461
 462   lost_fraction = lfExactlyZero;
 463   omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
 464   exponent += rhs.exponent;
 465
 466   if(addend) {
 467     Significand savedSignificand = significand;
 468     const fltSemantics *savedSemantics = semantics;
 469     fltSemantics extendedSemantics;
 470     opStatus status;
 471     unsigned int extendedPrecision;
 472
 473     /* Normalize our MSB.  */
 474     extendedPrecision = precision + precision - 1;
 475     if(omsb != extendedPrecision)
 476       {
 477         APInt::tcShiftLeft(fullSignificand, newPartsCount,
 478                            extendedPrecision - omsb);
 479         exponent -= extendedPrecision - omsb;
 480       }
 481
 482     /* Create new semantics.  */
 483     extendedSemantics = *semantics;
 484     extendedSemantics.precision = extendedPrecision;
 485
 486     if(newPartsCount == 1)
 487       significand.part = fullSignificand[0];
 488     else
 489       significand.parts = fullSignificand;
 490     semantics = &extendedSemantics;
 491
 492     APFloat extendedAddend(*addend);
 493     status = extendedAddend.convert(extendedSemantics, rmTowardZero);
 494     assert(status == opOK);
 495     lost_fraction = addOrSubtractSignificand(extendedAddend, false);
 496
 497     /* Restore our state.  */
 498     if(newPartsCount == 1)
 499       fullSignificand[0] = significand.part;
 500     significand = savedSignificand;
 501     semantics = savedSemantics;
 502
 503     omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
 504   }
 505
 506   exponent -= (precision - 1);
 507
 508   if(omsb > precision) {
 509     unsigned int bits, significantParts;
 510     lostFraction lf;
 511
 512     bits = omsb - precision;
 513     significantParts = partCountForBits(omsb);
 514     lf = shiftRight(fullSignificand, significantParts, bits);
 515     lost_fraction = combineLostFractions(lf, lost_fraction);
 516     exponent += bits;
 517   }
 518
 519   APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
 520
 521   if(newPartsCount > 4)
 522     delete [] fullSignificand;
 523
 524   return lost_fraction;
 525 }
 526
 527 /* Multiply the significands of LHS and RHS to DST.  */
 528 lostFraction
 529 APFloat::divideSignificand(const APFloat &rhs)
 530 {
 531   unsigned int bit, i, partsCount;
 532   const integerPart *rhsSignificand;
 533   integerPart *lhsSignificand, *dividend, *divisor;
 534   integerPart scratch[4];
 535   lostFraction lost_fraction;
 536
 537   assert(semantics == rhs.semantics);
 538
 539   lhsSignificand = significandParts();
 540   rhsSignificand = rhs.significandParts();
 541   partsCount = partCount();
 542
 543   if(partsCount > 2)
 544     dividend = new integerPart[partsCount * 2];
 545   else
 546     dividend = scratch;
 547
 548   divisor = dividend + partsCount;
 549
 550   /* Copy the dividend and divisor as they will be modified in-place.  */
 551   for(i = 0; i < partsCount; i++) {
 552     dividend[i] = lhsSignificand[i];
 553     divisor[i] = rhsSignificand[i];
 554     lhsSignificand[i] = 0;
 555   }
 556
 557   exponent -= rhs.exponent;
 558
 559   unsigned int precision = semantics->precision;
 560
 561   /* Normalize the divisor.  */
 562   bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
 563   if(bit) {
 564     exponent += bit;
 565     APInt::tcShiftLeft(divisor, partsCount, bit);
 566   }
 567
 568   /* Normalize the dividend.  */
 569   bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
 570   if(bit) {
 571     exponent -= bit;
 572     APInt::tcShiftLeft(dividend, partsCount, bit);
 573   }
 574
 575   if(APInt::tcCompare(dividend, divisor, partsCount) < 0) {
 576     exponent--;
 577     APInt::tcShiftLeft(dividend, partsCount, 1);
 578     assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
 579   }
 580
 581   /* Long division.  */
 582   for(bit = precision; bit; bit -= 1) {
 583     if(APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
 584       APInt::tcSubtract(dividend, divisor, 0, partsCount);
 585       APInt::tcSetBit(lhsSignificand, bit - 1);
 586     }
 587
 588     APInt::tcShiftLeft(dividend, partsCount, 1);
 589   }
 590
 591   /* Figure out the lost fraction.  */
 592   int cmp = APInt::tcCompare(dividend, divisor, partsCount);
 593
 594   if(cmp > 0)
 595     lost_fraction = lfMoreThanHalf;
 596   else if(cmp == 0)
 597     lost_fraction = lfExactlyHalf;
 598   else if(APInt::tcIsZero(dividend, partsCount))
 599     lost_fraction = lfExactlyZero;
 600   else
 601     lost_fraction = lfLessThanHalf;
 602
 603   if(partsCount > 2)
 604     delete [] dividend;
 605
 606   return lost_fraction;
 607 }
 608
 609 unsigned int
 610 APFloat::significandMSB() const
 611 {
 612   return APInt::tcMSB(significandParts(), partCount());
 613 }
 614
 615 unsigned int
 616 APFloat::significandLSB() const
 617 {
 618   return APInt::tcLSB(significandParts(), partCount());
 619 }
 620
 621 /* Note that a zero result is NOT normalized to fcZero.  */
 622 lostFraction
 623 APFloat::shiftSignificandRight(unsigned int bits)
 624 {
 625   /* Our exponent should not overflow.  */
 626   assert((exponent_t) (exponent + bits) >= exponent);
 627
 628   exponent += bits;
 629
 630   return shiftRight(significandParts(), partCount(), bits);
 631 }
 632
 633 /* Shift the significand left BITS bits, subtract BITS from its exponent.  */
 634 void
 635 APFloat::shiftSignificandLeft(unsigned int bits)
 636 {
 637   assert(bits < semantics->precision);
 638
 639   if(bits) {
 640     unsigned int partsCount = partCount();
 641
 642     APInt::tcShiftLeft(significandParts(), partsCount, bits);
 643     exponent -= bits;
 644
 645     assert(!APInt::tcIsZero(significandParts(), partsCount));
 646   }
 647 }
 648
 649 APFloat::cmpResult
 650 APFloat::compareAbsoluteValue(const APFloat &rhs) const
 651 {
 652   int compare;
 653
 654   assert(semantics == rhs.semantics);
 655   assert(category == fcNormal);
 656   assert(rhs.category == fcNormal);
 657
 658   compare = exponent - rhs.exponent;
 659
 660   /* If exponents are equal, do an unsigned bignum comparison of the
 661      significands.  */
 662   if(compare == 0)
 663     compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
 664                                partCount());
 665
 666   if(compare > 0)
 667     return cmpGreaterThan;
 668   else if(compare < 0)
 669     return cmpLessThan;
 670   else
 671     return cmpEqual;
 672 }
 673
 674 /* Handle overflow.  Sign is preserved.  We either become infinity or
 675    the largest finite number.  */
 676 APFloat::opStatus
 677 APFloat::handleOverflow(roundingMode rounding_mode)
 678 {
 679   /* Infinity?  */
 680   if(rounding_mode == rmNearestTiesToEven
 681      || rounding_mode == rmNearestTiesToAway
 682      || (rounding_mode == rmTowardPositive && !sign)
 683      || (rounding_mode == rmTowardNegative && sign))
 684     {
 685       category = fcInfinity;
 686       return (opStatus) (opOverflow | opInexact);
 687     }
 688
 689   /* Otherwise we become the largest finite number.  */
 690   category = fcNormal;
 691   exponent = semantics->maxExponent;
 692   APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
 693                                    semantics->precision);
 694
 695   return opInexact;
 696 }
 697
 698 /* This routine must work for fcZero of both signs, and fcNormal
 699    numbers.  */
 700 bool
 701 APFloat::roundAwayFromZero(roundingMode rounding_mode,
 702                            lostFraction lost_fraction)
 703 {
 704   /* QNaNs and infinities should not have lost fractions.  */
 705   assert(category == fcNormal || category == fcZero);
 706
 707   /* Our caller has already handled this case.  */
 708   assert(lost_fraction != lfExactlyZero);
 709
 710   switch(rounding_mode) {
 711   default:
 712     assert(0);
 713
 714   case rmNearestTiesToAway:
 715     return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
 716
 717   case rmNearestTiesToEven:
 718     if(lost_fraction == lfMoreThanHalf)
 719       return true;
 720
 721     /* Our zeroes don't have a significand to test.  */
 722     if(lost_fraction == lfExactlyHalf && category != fcZero)
 723       return significandParts()[0] & 1;
 724
 725     return false;
 726
 727   case rmTowardZero:
 728     return false;
 729
 730   case rmTowardPositive:
 731     return sign == false;
 732
 733   case rmTowardNegative:
 734     return sign == true;
 735   }
 736 }
 737
 738 APFloat::opStatus
 739 APFloat::normalize(roundingMode rounding_mode,
 740                    lostFraction lost_fraction)
 741 {
 742   unsigned int omsb;            /* One, not zero, based MSB.  */
 743   int exponentChange;
 744
 745   if(category != fcNormal)
 746     return opOK;
 747
 748   /* Before rounding normalize the exponent of fcNormal numbers.  */
 749   omsb = significandMSB() + 1;
 750
 751   if(omsb) {
 752     /* OMSB is numbered from 1.  We want to place it in the integer
 753        bit numbered PRECISON if possible, with a compensating change in
 754        the exponent.  */
 755     exponentChange = omsb - semantics->precision;
 756
 757     /* If the resulting exponent is too high, overflow according to
 758        the rounding mode.  */
 759     if(exponent + exponentChange > semantics->maxExponent)
 760       return handleOverflow(rounding_mode);
 761
 762     /* Subnormal numbers have exponent minExponent, and their MSB
 763        is forced based on that.  */
 764     if(exponent + exponentChange < semantics->minExponent)
 765       exponentChange = semantics->minExponent - exponent;
 766
 767     /* Shifting left is easy as we don't lose precision.  */
 768     if(exponentChange < 0) {
 769       assert(lost_fraction == lfExactlyZero);
 770
 771       shiftSignificandLeft(-exponentChange);
 772
 773       return opOK;
 774     }
 775
 776     if(exponentChange > 0) {
 777       lostFraction lf;
 778
 779       /* Shift right and capture any new lost fraction.  */
 780       lf = shiftSignificandRight(exponentChange);
 781
 782       lost_fraction = combineLostFractions(lf, lost_fraction);
 783
 784       /* Keep OMSB up-to-date.  */
 785       if(omsb > (unsigned) exponentChange)
 786         omsb -= (unsigned) exponentChange;
 787       else
 788         omsb = 0;
 789     }
 790   }
 791
 792   /* Now round the number according to rounding_mode given the lost
 793      fraction.  */
 794
 795   /* As specified in IEEE 754, since we do not trap we do not report
 796      underflow for exact results.  */
 797   if(lost_fraction == lfExactlyZero) {
 798     /* Canonicalize zeroes.  */
 799     if(omsb == 0)
 800       category = fcZero;
 801
 802     return opOK;
 803   }
 804
 805   /* Increment the significand if we're rounding away from zero.  */
 806   if(roundAwayFromZero(rounding_mode, lost_fraction)) {
 807     if(omsb == 0)
 808       exponent = semantics->minExponent;
 809
 810     incrementSignificand();
 811     omsb = significandMSB() + 1;
 812
 813     /* Did the significand increment overflow?  */
 814     if(omsb == (unsigned) semantics->precision + 1) {
 815       /* Renormalize by incrementing the exponent and shifting our
 816          significand right one.  However if we already have the
 817          maximum exponent we overflow to infinity.  */
 818       if(exponent == semantics->maxExponent) {
 819         category = fcInfinity;
 820
 821         return (opStatus) (opOverflow | opInexact);
 822       }
 823
 824       shiftSignificandRight(1);
 825
 826       return opInexact;
 827     }
 828   }
 829
 830   /* The normal case - we were and are not denormal, and any
 831      significand increment above didn't overflow.  */
 832   if(omsb == semantics->precision)
 833     return opInexact;
 834
 835   /* We have a non-zero denormal.  */
 836   assert(omsb < semantics->precision);
 837   assert(exponent == semantics->minExponent);
 838
 839   /* Canonicalize zeroes.  */
 840   if(omsb == 0)
 841     category = fcZero;
 842
 843   /* The fcZero case is a denormal that underflowed to zero.  */
 844   return (opStatus) (opUnderflow | opInexact);
 845 }
 846
 847 APFloat::opStatus
 848 APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
 849 {
 850   switch(convolve(category, rhs.category)) {
 851   default:
 852     assert(0);
 853
 854   case convolve(fcQNaN, fcZero):
 855   case convolve(fcQNaN, fcNormal):
 856   case convolve(fcQNaN, fcInfinity):
 857   case convolve(fcQNaN, fcQNaN):
 858   case convolve(fcNormal, fcZero):
 859   case convolve(fcInfinity, fcNormal):
 860   case convolve(fcInfinity, fcZero):
 861     return opOK;
 862
 863   case convolve(fcZero, fcQNaN):
 864   case convolve(fcNormal, fcQNaN):
 865   case convolve(fcInfinity, fcQNaN):
 866     category = fcQNaN;
 867     return opOK;
 868
 869   case convolve(fcNormal, fcInfinity):
 870   case convolve(fcZero, fcInfinity):
 871     category = fcInfinity;
 872     sign = rhs.sign ^ subtract;
 873     return opOK;
 874
 875   case convolve(fcZero, fcNormal):
 876     assign(rhs);
 877     sign = rhs.sign ^ subtract;
 878     return opOK;
 879
 880   case convolve(fcZero, fcZero):
 881     /* Sign depends on rounding mode; handled by caller.  */
 882     return opOK;
 883
 884   case convolve(fcInfinity, fcInfinity):
 885     /* Differently signed infinities can only be validly
 886        subtracted.  */
 887     if(sign ^ rhs.sign != subtract) {
 888       category = fcQNaN;
 889       return opInvalidOp;
 890     }
 891
 892     return opOK;
 893
 894   case convolve(fcNormal, fcNormal):
 895     return opDivByZero;
 896   }
 897 }
 898
 899 /* Add or subtract two normal numbers.  */
 900 lostFraction
 901 APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
 902 {
 903   integerPart carry;
 904   lostFraction lost_fraction;
 905   int bits;
 906
 907   /* Determine if the operation on the absolute values is effectively
 908      an addition or subtraction.  */
 909   subtract ^= (sign ^ rhs.sign);
 910
 911   /* Are we bigger exponent-wise than the RHS?  */
 912   bits = exponent - rhs.exponent;
 913
 914   /* Subtraction is more subtle than one might naively expect.  */
 915   if(subtract) {
 916     APFloat temp_rhs(rhs);
 917     bool reverse;
 918
 919     if (bits == 0) {
 920       reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
 921       lost_fraction = lfExactlyZero;
 922     } else if (bits > 0) {
 923       lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
 924       shiftSignificandLeft(1);
 925       reverse = false;
 926     } else {
 927       lost_fraction = shiftSignificandRight(-bits - 1);
 928       temp_rhs.shiftSignificandLeft(1);
 929       reverse = true;
 930     }
 931
 932     if (reverse) {
 933       carry = temp_rhs.subtractSignificand
 934         (*this, lost_fraction != lfExactlyZero);
 935       copySignificand(temp_rhs);
 936       sign = !sign;
 937     } else {
 938       carry = subtractSignificand
 939         (temp_rhs, lost_fraction != lfExactlyZero);
 940     }
 941
 942     /* Invert the lost fraction - it was on the RHS and
 943        subtracted.  */
 944     if(lost_fraction == lfLessThanHalf)
 945       lost_fraction = lfMoreThanHalf;
 946     else if(lost_fraction == lfMoreThanHalf)
 947       lost_fraction = lfLessThanHalf;
 948
 949     /* The code above is intended to ensure that no borrow is
 950        necessary.  */
 951     assert(!carry);
 952   } else {
 953     if(bits > 0) {
 954       APFloat temp_rhs(rhs);
 955
 956       lost_fraction = temp_rhs.shiftSignificandRight(bits);
 957       carry = addSignificand(temp_rhs);
 958     } else {
 959       lost_fraction = shiftSignificandRight(-bits);
 960       carry = addSignificand(rhs);
 961     }
 962
 963     /* We have a guard bit; generating a carry cannot happen.  */
 964     assert(!carry);
 965   }
 966
 967   return lost_fraction;
 968 }
 969
 970 APFloat::opStatus
 971 APFloat::multiplySpecials(const APFloat &rhs)
 972 {
 973   switch(convolve(category, rhs.category)) {
 974   default:
 975     assert(0);
 976
 977   case convolve(fcQNaN, fcZero):
 978   case convolve(fcQNaN, fcNormal):
 979   case convolve(fcQNaN, fcInfinity):
 980   case convolve(fcQNaN, fcQNaN):
 981   case convolve(fcZero, fcQNaN):
 982   case convolve(fcNormal, fcQNaN):
 983   case convolve(fcInfinity, fcQNaN):
 984     category = fcQNaN;
 985     return opOK;
 986
 987   case convolve(fcNormal, fcInfinity):
 988   case convolve(fcInfinity, fcNormal):
 989   case convolve(fcInfinity, fcInfinity):
 990     category = fcInfinity;
 991     return opOK;
 992
 993   case convolve(fcZero, fcNormal):
 994   case convolve(fcNormal, fcZero):
 995   case convolve(fcZero, fcZero):
 996     category = fcZero;
 997     return opOK;
 998
 999   case convolve(fcZero, fcInfinity):
1000   case convolve(fcInfinity, fcZero):
1001     category = fcQNaN;
1002     return opInvalidOp;
1003
1004   case convolve(fcNormal, fcNormal):
1005     return opOK;
1006   }
1007 }
1008
1009 APFloat::opStatus
1010 APFloat::divideSpecials(const APFloat &rhs)
1011 {
1012   switch(convolve(category, rhs.category)) {
1013   default:
1014     assert(0);
1015
1016   case convolve(fcQNaN, fcZero):
1017   case convolve(fcQNaN, fcNormal):
1018   case convolve(fcQNaN, fcInfinity):
1019   case convolve(fcQNaN, fcQNaN):
1020   case convolve(fcInfinity, fcZero):
1021   case convolve(fcInfinity, fcNormal):
1022   case convolve(fcZero, fcInfinity):
1023   case convolve(fcZero, fcNormal):
1024     return opOK;
1025
1026   case convolve(fcZero, fcQNaN):
1027   case convolve(fcNormal, fcQNaN):
1028   case convolve(fcInfinity, fcQNaN):
1029     category = fcQNaN;
1030     return opOK;
1031
1032   case convolve(fcNormal, fcInfinity):
1033     category = fcZero;
1034     return opOK;
1035
1036   case convolve(fcNormal, fcZero):
1037     category = fcInfinity;
1038     return opDivByZero;
1039
1040   case convolve(fcInfinity, fcInfinity):
1041   case convolve(fcZero, fcZero):
1042     category = fcQNaN;
1043     return opInvalidOp;
1044
1045   case convolve(fcNormal, fcNormal):
1046     return opOK;
1047   }
1048 }
1049
1050 /* Change sign.  */
1051 void
1052 APFloat::changeSign()
1053 {
1054   /* Look mummy, this one's easy.  */
1055   sign = !sign;
1056 }
1057
1058 /* Normalized addition or subtraction.  */
1059 APFloat::opStatus
1060 APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1061                        bool subtract)
1062 {
1063   opStatus fs;
1064
1065   fs = addOrSubtractSpecials(rhs, subtract);
1066
1067   /* This return code means it was not a simple case.  */
1068   if(fs == opDivByZero) {
1069     lostFraction lost_fraction;
1070
1071     lost_fraction = addOrSubtractSignificand(rhs, subtract);
1072     fs = normalize(rounding_mode, lost_fraction);
1073
1074     /* Can only be zero if we lost no fraction.  */
1075     assert(category != fcZero || lost_fraction == lfExactlyZero);
1076   }
1077
1078   /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1079      positive zero unless rounding to minus infinity, except that
1080      adding two like-signed zeroes gives that zero.  */
1081   if(category == fcZero) {
1082     if(rhs.category != fcZero || (sign == rhs.sign) == subtract)
1083       sign = (rounding_mode == rmTowardNegative);
1084   }
1085
1086   return fs;
1087 }
1088
1089 /* Normalized addition.  */
1090 APFloat::opStatus
1091 APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1092 {
1093   return addOrSubtract(rhs, rounding_mode, false);
1094 }
1095
1096 /* Normalized subtraction.  */
1097 APFloat::opStatus
1098 APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1099 {
1100   return addOrSubtract(rhs, rounding_mode, true);
1101 }
1102
1103 /* Normalized multiply.  */
1104 APFloat::opStatus
1105 APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1106 {
1107   opStatus fs;
1108
1109   sign ^= rhs.sign;
1110   fs = multiplySpecials(rhs);
1111
1112   if(category == fcNormal) {
1113     lostFraction lost_fraction = multiplySignificand(rhs, 0);
1114     fs = normalize(rounding_mode, lost_fraction);
1115     if(lost_fraction != lfExactlyZero)
1116       fs = (opStatus) (fs | opInexact);
1117   }
1118
1119   return fs;
1120 }
1121
1122 /* Normalized divide.  */
1123 APFloat::opStatus
1124 APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1125 {
1126   opStatus fs;
1127
1128   sign ^= rhs.sign;
1129   fs = divideSpecials(rhs);
1130
1131   if(category == fcNormal) {
1132     lostFraction lost_fraction = divideSignificand(rhs);
1133     fs = normalize(rounding_mode, lost_fraction);
1134     if(lost_fraction != lfExactlyZero)
1135       fs = (opStatus) (fs | opInexact);
1136   }
1137
1138   return fs;
1139 }
1140
1141 /* Normalized fused-multiply-add.  */
1142 APFloat::opStatus
1143 APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1144                           const APFloat &addend,
1145                           roundingMode rounding_mode)
1146 {
1147   opStatus fs;
1148
1149   /* Post-multiplication sign, before addition.  */
1150   sign ^= multiplicand.sign;
1151
1152   /* If and only if all arguments are normal do we need to do an
1153      extended-precision calculation.  */
1154   if(category == fcNormal
1155      && multiplicand.category == fcNormal
1156      && addend.category == fcNormal) {
1157     lostFraction lost_fraction;
1158
1159     lost_fraction = multiplySignificand(multiplicand, &addend);
1160     fs = normalize(rounding_mode, lost_fraction);
1161     if(lost_fraction != lfExactlyZero)
1162       fs = (opStatus) (fs | opInexact);
1163
1164     /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1165        positive zero unless rounding to minus infinity, except that
1166        adding two like-signed zeroes gives that zero.  */
1167     if(category == fcZero && sign != addend.sign)
1168       sign = (rounding_mode == rmTowardNegative);
1169   } else {
1170     fs = multiplySpecials(multiplicand);
1171
1172     /* FS can only be opOK or opInvalidOp.  There is no more work
1173        to do in the latter case.  The IEEE-754R standard says it is
1174        implementation-defined in this case whether, if ADDEND is a
1175        quiet QNaN, we raise invalid op; this implementation does so.
1176
1177        If we need to do the addition we can do so with normal
1178        precision.  */
1179     if(fs == opOK)
1180       fs = addOrSubtract(addend, rounding_mode, false);
1181   }
1182
1183   return fs;
1184 }
1185
1186 /* Comparison requires normalized numbers.  */
1187 APFloat::cmpResult
1188 APFloat::compare(const APFloat &rhs) const
1189 {
1190   cmpResult result;
1191
1192   assert(semantics == rhs.semantics);
1193
1194   switch(convolve(category, rhs.category)) {
1195   default:
1196     assert(0);
1197
1198   case convolve(fcQNaN, fcZero):
1199   case convolve(fcQNaN, fcNormal):
1200   case convolve(fcQNaN, fcInfinity):
1201   case convolve(fcQNaN, fcQNaN):
1202   case convolve(fcZero, fcQNaN):
1203   case convolve(fcNormal, fcQNaN):
1204   case convolve(fcInfinity, fcQNaN):
1205     return cmpUnordered;
1206
1207   case convolve(fcInfinity, fcNormal):
1208   case convolve(fcInfinity, fcZero):
1209   case convolve(fcNormal, fcZero):
1210     if(sign)
1211       return cmpLessThan;
1212     else
1213       return cmpGreaterThan;
1214
1215   case convolve(fcNormal, fcInfinity):
1216   case convolve(fcZero, fcInfinity):
1217   case convolve(fcZero, fcNormal):
1218     if(rhs.sign)
1219       return cmpGreaterThan;
1220     else
1221       return cmpLessThan;
1222
1223   case convolve(fcInfinity, fcInfinity):
1224     if(sign == rhs.sign)
1225       return cmpEqual;
1226     else if(sign)
1227       return cmpLessThan;
1228     else
1229       return cmpGreaterThan;
1230
1231   case convolve(fcZero, fcZero):
1232     return cmpEqual;
1233
1234   case convolve(fcNormal, fcNormal):
1235     break;
1236   }
1237
1238   /* Two normal numbers.  Do they have the same sign?  */
1239   if(sign != rhs.sign) {
1240     if(sign)
1241       result = cmpLessThan;
1242     else
1243       result = cmpGreaterThan;
1244   } else {
1245     /* Compare absolute values; invert result if negative.  */
1246     result = compareAbsoluteValue(rhs);
1247
1248     if(sign) {
1249       if(result == cmpLessThan)
1250         result = cmpGreaterThan;
1251       else if(result == cmpGreaterThan)
1252         result = cmpLessThan;
1253     }
1254   }
1255
1256   return result;
1257 }
1258
1259 APFloat::opStatus
1260 APFloat::convert(const fltSemantics &toSemantics,
1261                  roundingMode rounding_mode)
1262 {
1263   unsigned int newPartCount;
1264   opStatus fs;
1265
1266   newPartCount = partCountForBits(toSemantics.precision + 1);
1267
1268   /* If our new form is wider, re-allocate our bit pattern into wider
1269      storage.  */
1270   if(newPartCount > partCount()) {
1271     integerPart *newParts;
1272
1273     newParts = new integerPart[newPartCount];
1274     APInt::tcSet(newParts, 0, newPartCount);
1275     APInt::tcAssign(newParts, significandParts(), partCount());
1276     freeSignificand();
1277     significand.parts = newParts;
1278   }
1279
1280   if(category == fcNormal) {
1281     /* Re-interpret our bit-pattern.  */
1282     exponent += toSemantics.precision - semantics->precision;
1283     semantics = &toSemantics;
1284     fs = normalize(rounding_mode, lfExactlyZero);
1285   } else {
1286     semantics = &toSemantics;
1287     fs = opOK;
1288   }
1289
1290   return fs;
1291 }
1292
1293 /* Convert a floating point number to an integer according to the
1294    rounding mode.  If the rounded integer value is out of range this
1295    returns an invalid operation exception.  If the rounded value is in
1296    range but the floating point number is not the exact integer, the C
1297    standard doesn't require an inexact exception to be raised.  IEEE
1298    854 does require it so we do that.
1299
1300    Note that for conversions to integer type the C standard requires
1301    round-to-zero to always be used.  */
1302 APFloat::opStatus
1303 APFloat::convertToInteger(integerPart *parts, unsigned int width,
1304                           bool isSigned,
1305                           roundingMode rounding_mode) const
1306 {
1307   lostFraction lost_fraction;
1308   unsigned int msb, partsCount;
1309   int bits;
1310
1311   /* Handle the three special cases first.  */
1312   if(category == fcInfinity || category == fcQNaN)
1313     return opInvalidOp;
1314
1315   partsCount = partCountForBits(width);
1316
1317   if(category == fcZero) {
1318     APInt::tcSet(parts, 0, partsCount);
1319     return opOK;
1320   }
1321
1322   /* Shift the bit pattern so the fraction is lost.  */
1323   APFloat tmp(*this);
1324
1325   bits = (int) semantics->precision - 1 - exponent;
1326
1327   if(bits > 0) {
1328     lost_fraction = tmp.shiftSignificandRight(bits);
1329   } else {
1330     tmp.shiftSignificandLeft(-bits);
1331     lost_fraction = lfExactlyZero;
1332   }
1333
1334   if(lost_fraction != lfExactlyZero
1335      && tmp.roundAwayFromZero(rounding_mode, lost_fraction))
1336     tmp.incrementSignificand();
1337
1338   msb = tmp.significandMSB();
1339
1340   /* Negative numbers cannot be represented as unsigned.  */
1341   if(!isSigned && tmp.sign && msb != -1U)
1342     return opInvalidOp;
1343
1344   /* It takes exponent + 1 bits to represent the truncated floating
1345      point number without its sign.  We lose a bit for the sign, but
1346      the maximally negative integer is a special case.  */
1347   if(msb + 1 > width)           /* !! Not same as msb >= width !! */
1348     return opInvalidOp;
1349
1350   if(isSigned && msb + 1 == width
1351      && (!tmp.sign || tmp.significandLSB() != msb))
1352     return opInvalidOp;
1353
1354   APInt::tcAssign(parts, tmp.significandParts(), partsCount);
1355
1356   if(tmp.sign)
1357     APInt::tcNegate(parts, partsCount);
1358
1359   if(lost_fraction == lfExactlyZero)
1360     return opOK;
1361   else
1362     return opInexact;
1363 }
1364
1365 APFloat::opStatus
1366 APFloat::convertFromUnsignedInteger(integerPart *parts,
1367                                     unsigned int partCount,
1368                                     roundingMode rounding_mode)
1369 {
1370   unsigned int msb, precision;
1371   lostFraction lost_fraction;
1372
1373   msb = APInt::tcMSB(parts, partCount) + 1;
1374   precision = semantics->precision;
1375
1376   category = fcNormal;
1377   exponent = precision - 1;
1378
1379   if(msb > precision) {
1380     exponent += (msb - precision);
1381     lost_fraction = shiftRight(parts, partCount, msb - precision);
1382     msb = precision;
1383   } else
1384     lost_fraction = lfExactlyZero;
1385
1386   /* Copy the bit image.  */
1387   zeroSignificand();
1388   APInt::tcAssign(significandParts(), parts, partCountForBits(msb));
1389
1390   return normalize(rounding_mode, lost_fraction);
1391 }
1392
1393 APFloat::opStatus
1394 APFloat::convertFromInteger(const integerPart *parts,
1395                             unsigned int partCount, bool isSigned,
1396                             roundingMode rounding_mode)
1397 {
1398   unsigned int width;
1399   opStatus status;
1400   integerPart *copy;
1401
1402   copy = new integerPart[partCount];
1403   APInt::tcAssign(copy, parts, partCount);
1404
1405   width = partCount * integerPartWidth;
1406
1407   sign = false;
1408   if(isSigned && APInt::tcExtractBit(parts, width - 1)) {
1409     sign = true;
1410     APInt::tcNegate(copy, partCount);
1411   }
1412
1413   status = convertFromUnsignedInteger(copy, partCount, rounding_mode);
1414   delete [] copy;
1415
1416   return status;
1417 }
1418
1419 APFloat::opStatus
1420 APFloat::convertFromHexadecimalString(const char *p,
1421                                       roundingMode rounding_mode)
1422 {
1423   lostFraction lost_fraction;
1424   integerPart *significand;
1425   unsigned int bitPos, partsCount;
1426   const char *dot, *firstSignificantDigit;
1427
1428   zeroSignificand();
1429   exponent = 0;
1430   category = fcNormal;
1431
1432   significand = significandParts();
1433   partsCount = partCount();
1434   bitPos = partsCount * integerPartWidth;
1435
1436   /* Skip leading zeroes and any(hexa)decimal point.  */
1437   p = skipLeadingZeroesAndAnyDot(p, &dot);
1438   firstSignificantDigit = p;
1439
1440   for(;;) {
1441     integerPart hex_value;
1442
1443     if(*p == '.') {
1444       assert(dot == 0);
1445       dot = p++;
1446     }
1447
1448     hex_value = hexDigitValue(*p);
1449     if(hex_value == -1U) {
1450       lost_fraction = lfExactlyZero;
1451       break;
1452     }
1453
1454     p++;
1455
1456     /* Store the number whilst 4-bit nibbles remain.  */
1457     if(bitPos) {
1458       bitPos -= 4;
1459       hex_value <<= bitPos % integerPartWidth;
1460       significand[bitPos / integerPartWidth] |= hex_value;
1461     } else {
1462       lost_fraction = trailingHexadecimalFraction(p, hex_value);
1463       while(hexDigitValue(*p) != -1U)
1464         p++;
1465       break;
1466     }
1467   }
1468
1469   /* Hex floats require an exponent but not a hexadecimal point.  */
1470   assert(*p == 'p' || *p == 'P');
1471
1472   /* Ignore the exponent if we are zero.  */
1473   if(p != firstSignificantDigit) {
1474     int expAdjustment;
1475
1476     /* Implicit hexadecimal point?  */
1477     if(!dot)
1478       dot = p;
1479
1480     /* Calculate the exponent adjustment implicit in the number of
1481        significant digits.  */
1482     expAdjustment = dot - firstSignificantDigit;
1483     if(expAdjustment < 0)
1484       expAdjustment++;
1485     expAdjustment = expAdjustment * 4 - 1;
1486
1487     /* Adjust for writing the significand starting at the most
1488        significant nibble.  */
1489     expAdjustment += semantics->precision;
1490     expAdjustment -= partsCount * integerPartWidth;
1491
1492     /* Adjust for the given exponent.  */
1493     exponent = totalExponent(p, expAdjustment);
1494   }
1495
1496   return normalize(rounding_mode, lost_fraction);
1497 }
1498
1499 APFloat::opStatus
1500 APFloat::convertFromString(const char *p, roundingMode rounding_mode) {
1501   /* Handle a leading minus sign.  */
1502   if(*p == '-')
1503     sign = 1, p++;
1504   else
1505     sign = 0;
1506
1507   if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
1508     return convertFromHexadecimalString(p + 2, rounding_mode);
1509
1510   assert(0 && "Decimal to binary conversions not yet implemented");
1511   abort();
1512 }
1513
1514 // For good performance it is desirable for different APFloats
1515 // to produce different integers.
1516 uint32_t
1517 APFloat::getHashValue() const {
1518   if (category==fcZero) return sign<<8 | semantics->precision ;
1519   else if (category==fcInfinity) return sign<<9 | semantics->precision;
1520   else if (category==fcQNaN) return 1<<10 | semantics->precision;
1521   else {
1522     uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
1523     const integerPart* p = significandParts();
1524     for (int i=partCount(); i>0; i--, p++)
1525       hash ^= ((uint32_t)*p) ^ (*p)>>32;
1526     return hash;
1527   }
1528 }
1529
1530 // Conversion from APFloat to/from host float/double.  It may eventually be
1531 // possible to eliminate these and have everybody deal with APFloats, but that
1532 // will take a while.  This approach will not easily extend to long double.
1533 // Current implementation requires partCount()==1, which is correct at the
1534 // moment but could be made more general.
1535
1536 double
1537 APFloat::convertToDouble() const {
1538   assert(semantics == (const llvm::fltSemantics* const)&IEEEdouble);
1539   assert (partCount()==1);
1540
1541   uint64_t myexponent, mysign, mysignificand;
1542
1543   if (category==fcNormal) {
1544     mysign = sign;
1545     mysignificand = *significandParts();
1546     myexponent = exponent+1023; //bias
1547   } else if (category==fcZero) {
1548     mysign = sign;
1549     myexponent = 0;
1550     mysignificand = 0;
1551   } else if (category==fcInfinity) {
1552     mysign = sign;
1553     myexponent = 0x7ff;
1554     mysignificand = 0;
1555   } else if (category==fcQNaN) {
1556     mysign = 0;
1557     myexponent = 0x7ff;
1558     mysignificand = 0xfffffffffffffLL;
1559   } else
1560     assert(0);
1561
1562   return BitsToDouble(((mysign & 1) << 63) | ((myexponent & 0x7ff) <<  52) |
1563         (mysignificand & 0xfffffffffffffLL));
1564 }
1565
1566 float
1567 APFloat::convertToFloat() const {
1568   assert(semantics == (const llvm::fltSemantics* const)&IEEEsingle);
1569   assert (partCount()==1);
1570
1571   uint32_t mysign, myexponent, mysignificand;
1572
1573   if (category==fcNormal) {
1574     mysign = sign;
1575     myexponent = exponent+127; //bias
1576     mysignificand = *significandParts();
1577   } else if (category==fcZero) {
1578     mysign = sign;
1579     myexponent = 0;
1580     mysignificand = 0;
1581   } else if (category==fcInfinity) {
1582     mysign = sign;
1583     myexponent = 0xff;
1584     mysignificand = 0;
1585   } else if (category==fcQNaN) {
1586     mysign = sign;
1587     myexponent = 0x7ff;
1588     mysignificand = 0x7fffff;
1589   } else
1590     assert(0);
1591
1592   return BitsToFloat(((mysign&1) << 31) | ((myexponent&0xff) << 23) |
1593         (mysignificand & 0x7fffff));
1594 }
1595
1596 APFloat::APFloat(double d) {
1597   uint64_t i = DoubleToBits(d);
1598   uint64_t mysign = i >> 63;
1599   uint64_t myexponent = (i >> 52) & 0x7ff;
1600   uint64_t mysignificand = i & 0xfffffffffffffLL;
1601
1602   initialize(&APFloat::IEEEdouble);
1603   assert(partCount()==1);
1604
1605   if (myexponent==0 && mysignificand==0) {
1606     // exponent, significand meaningless
1607     category = fcZero;
1608     sign = mysign;
1609   } else if (myexponent==0x7ff && mysignificand==0) {
1610     // exponent, significand meaningless
1611     category = fcInfinity;
1612     sign = mysign;
1613   } else if (myexponent==0x7ff && (mysignificand & 0x8000000000000LL)) {
1614     // sign, exponent, significand meaningless
1615     category = fcQNaN;
1616   } else {
1617     sign = mysign;
1618     category = fcNormal;
1619     exponent = myexponent - 1023;
1620     *significandParts() = mysignificand | 0x100000000000000LL;
1621   }
1622 }
1623
1624 APFloat::APFloat(float f) {
1625   uint32_t i = FloatToBits(f);
1626   uint32_t mysign = i >> 31;
1627   uint32_t myexponent = (i >> 23) & 0xff;
1628   uint32_t mysignificand = i & 0x7fffff;
1629
1630   initialize(&APFloat::IEEEsingle);
1631   assert(partCount()==1);
1632
1633   if (myexponent==0 && mysignificand==0) {
1634     // exponent, significand meaningless
1635     category = fcZero;
1636     sign = mysign;
1637   } else if (myexponent==0xff && mysignificand==0) {
1638     // exponent, significand meaningless
1639     category = fcInfinity;
1640     sign = mysign;
1641   } else if (myexponent==0xff && (mysignificand & 0x400000)) {
1642     // sign, exponent, significand meaningless
1643     category = fcQNaN;
1644   } else {
1645     category = fcNormal;
1646     sign = mysign;
1647     exponent = myexponent - 127;  //bias
1648     *significandParts() = mysignificand | 0x800000; // integer bit
1649   }
1650 }