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