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 || category == fcNaN)
 251     copySignificand(rhs);
 252 }
 253
 254 void
 255 APFloat::copySignificand(const APFloat &rhs)
 256 {
 257   assert(category == fcNormal || category == fcNaN);
 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::bitwiseIsEqual(const APFloat &rhs) const {
 280   if (this == &rhs)
 281     return true;
 282   if (semantics != rhs.semantics ||
 283       category != rhs.category ||
 284       sign != rhs.sign)
 285     return false;
 286   if (category==fcZero || category==fcInfinity)
 287     return true;
 288   else if (category==fcNormal && exponent!=rhs.exponent)
 289     return false;
 290   else {
 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 || category == fcNaN);
 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   /* NaNs 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(fcNaN, fcZero):
 854   case convolve(fcNaN, fcNormal):
 855   case convolve(fcNaN, fcInfinity):
 856   case convolve(fcNaN, fcNaN):
 857   case convolve(fcNormal, fcZero):
 858   case convolve(fcInfinity, fcNormal):
 859   case convolve(fcInfinity, fcZero):
 860     return opOK;
 861
 862   case convolve(fcZero, fcNaN):
 863   case convolve(fcNormal, fcNaN):
 864   case convolve(fcInfinity, fcNaN):
 865     category = fcNaN;
 866     copySignificand(rhs);
 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 = fcNaN;
 889       // Arbitrary but deterministic value for significand
 890       APInt::tcSet(significandParts(), ~0U, partCount());
 891       return opInvalidOp;
 892     }
 893
 894     return opOK;
 895
 896   case convolve(fcNormal, fcNormal):
 897     return opDivByZero;
 898   }
 899 }
 900
 901 /* Add or subtract two normal numbers.  */
 902 lostFraction
 903 APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
 904 {
 905   integerPart carry;
 906   lostFraction lost_fraction;
 907   int bits;
 908
 909   /* Determine if the operation on the absolute values is effectively
 910      an addition or subtraction.  */
 911   subtract ^= (sign ^ rhs.sign);
 912
 913   /* Are we bigger exponent-wise than the RHS?  */
 914   bits = exponent - rhs.exponent;
 915
 916   /* Subtraction is more subtle than one might naively expect.  */
 917   if(subtract) {
 918     APFloat temp_rhs(rhs);
 919     bool reverse;
 920
 921     if (bits == 0) {
 922       reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
 923       lost_fraction = lfExactlyZero;
 924     } else if (bits > 0) {
 925       lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
 926       shiftSignificandLeft(1);
 927       reverse = false;
 928     } else {
 929       lost_fraction = shiftSignificandRight(-bits - 1);
 930       temp_rhs.shiftSignificandLeft(1);
 931       reverse = true;
 932     }
 933
 934     if (reverse) {
 935       carry = temp_rhs.subtractSignificand
 936         (*this, lost_fraction != lfExactlyZero);
 937       copySignificand(temp_rhs);
 938       sign = !sign;
 939     } else {
 940       carry = subtractSignificand
 941         (temp_rhs, lost_fraction != lfExactlyZero);
 942     }
 943
 944     /* Invert the lost fraction - it was on the RHS and
 945        subtracted.  */
 946     if(lost_fraction == lfLessThanHalf)
 947       lost_fraction = lfMoreThanHalf;
 948     else if(lost_fraction == lfMoreThanHalf)
 949       lost_fraction = lfLessThanHalf;
 950
 951     /* The code above is intended to ensure that no borrow is
 952        necessary.  */
 953     assert(!carry);
 954   } else {
 955     if(bits > 0) {
 956       APFloat temp_rhs(rhs);
 957
 958       lost_fraction = temp_rhs.shiftSignificandRight(bits);
 959       carry = addSignificand(temp_rhs);
 960     } else {
 961       lost_fraction = shiftSignificandRight(-bits);
 962       carry = addSignificand(rhs);
 963     }
 964
 965     /* We have a guard bit; generating a carry cannot happen.  */
 966     assert(!carry);
 967   }
 968
 969   return lost_fraction;
 970 }
 971
 972 APFloat::opStatus
 973 APFloat::multiplySpecials(const APFloat &rhs)
 974 {
 975   switch(convolve(category, rhs.category)) {
 976   default:
 977     assert(0);
 978
 979   case convolve(fcNaN, fcZero):
 980   case convolve(fcNaN, fcNormal):
 981   case convolve(fcNaN, fcInfinity):
 982   case convolve(fcNaN, fcNaN):
 983     return opOK;
 984
 985   case convolve(fcZero, fcNaN):
 986   case convolve(fcNormal, fcNaN):
 987   case convolve(fcInfinity, fcNaN):
 988     category = fcNaN;
 989     copySignificand(rhs);
 990     return opOK;
 991
 992   case convolve(fcNormal, fcInfinity):
 993   case convolve(fcInfinity, fcNormal):
 994   case convolve(fcInfinity, fcInfinity):
 995     category = fcInfinity;
 996     return opOK;
 997
 998   case convolve(fcZero, fcNormal):
 999   case convolve(fcNormal, fcZero):
1000   case convolve(fcZero, fcZero):
1001     category = fcZero;
1002     return opOK;
1003
1004   case convolve(fcZero, fcInfinity):
1005   case convolve(fcInfinity, fcZero):
1006     category = fcNaN;
1007     // Arbitrary but deterministic value for significand
1008     APInt::tcSet(significandParts(), ~0U, partCount());
1009     return opInvalidOp;
1010
1011   case convolve(fcNormal, fcNormal):
1012     return opOK;
1013   }
1014 }
1015
1016 APFloat::opStatus
1017 APFloat::divideSpecials(const APFloat &rhs)
1018 {
1019   switch(convolve(category, rhs.category)) {
1020   default:
1021     assert(0);
1022
1023   case convolve(fcNaN, fcZero):
1024   case convolve(fcNaN, fcNormal):
1025   case convolve(fcNaN, fcInfinity):
1026   case convolve(fcNaN, fcNaN):
1027   case convolve(fcInfinity, fcZero):
1028   case convolve(fcInfinity, fcNormal):
1029   case convolve(fcZero, fcInfinity):
1030   case convolve(fcZero, fcNormal):
1031     return opOK;
1032
1033   case convolve(fcZero, fcNaN):
1034   case convolve(fcNormal, fcNaN):
1035   case convolve(fcInfinity, fcNaN):
1036     category = fcNaN;
1037     copySignificand(rhs);
1038     return opOK;
1039
1040   case convolve(fcNormal, fcInfinity):
1041     category = fcZero;
1042     return opOK;
1043
1044   case convolve(fcNormal, fcZero):
1045     category = fcInfinity;
1046     return opDivByZero;
1047
1048   case convolve(fcInfinity, fcInfinity):
1049   case convolve(fcZero, fcZero):
1050     category = fcNaN;
1051     // Arbitrary but deterministic value for significand
1052     APInt::tcSet(significandParts(), ~0U, partCount());
1053     return opInvalidOp;
1054
1055   case convolve(fcNormal, fcNormal):
1056     return opOK;
1057   }
1058 }
1059
1060 /* Change sign.  */
1061 void
1062 APFloat::changeSign()
1063 {
1064   /* Look mummy, this one's easy.  */
1065   sign = !sign;
1066 }
1067
1068 void
1069 APFloat::clearSign()
1070 {
1071   /* So is this one. */
1072   sign = 0;
1073 }
1074
1075 void
1076 APFloat::copySign(const APFloat &rhs)
1077 {
1078   /* And this one. */
1079   sign = rhs.sign;
1080 }
1081
1082 /* Normalized addition or subtraction.  */
1083 APFloat::opStatus
1084 APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1085                        bool subtract)
1086 {
1087   opStatus fs;
1088
1089   fs = addOrSubtractSpecials(rhs, subtract);
1090
1091   /* This return code means it was not a simple case.  */
1092   if(fs == opDivByZero) {
1093     lostFraction lost_fraction;
1094
1095     lost_fraction = addOrSubtractSignificand(rhs, subtract);
1096     fs = normalize(rounding_mode, lost_fraction);
1097
1098     /* Can only be zero if we lost no fraction.  */
1099     assert(category != fcZero || lost_fraction == lfExactlyZero);
1100   }
1101
1102   /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1103      positive zero unless rounding to minus infinity, except that
1104      adding two like-signed zeroes gives that zero.  */
1105   if(category == fcZero) {
1106     if(rhs.category != fcZero || (sign == rhs.sign) == subtract)
1107       sign = (rounding_mode == rmTowardNegative);
1108   }
1109
1110   return fs;
1111 }
1112
1113 /* Normalized addition.  */
1114 APFloat::opStatus
1115 APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1116 {
1117   return addOrSubtract(rhs, rounding_mode, false);
1118 }
1119
1120 /* Normalized subtraction.  */
1121 APFloat::opStatus
1122 APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1123 {
1124   return addOrSubtract(rhs, rounding_mode, true);
1125 }
1126
1127 /* Normalized multiply.  */
1128 APFloat::opStatus
1129 APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1130 {
1131   opStatus fs;
1132
1133   sign ^= rhs.sign;
1134   fs = multiplySpecials(rhs);
1135
1136   if(category == fcNormal) {
1137     lostFraction lost_fraction = multiplySignificand(rhs, 0);
1138     fs = normalize(rounding_mode, lost_fraction);
1139     if(lost_fraction != lfExactlyZero)
1140       fs = (opStatus) (fs | opInexact);
1141   }
1142
1143   return fs;
1144 }
1145
1146 /* Normalized divide.  */
1147 APFloat::opStatus
1148 APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1149 {
1150   opStatus fs;
1151
1152   sign ^= rhs.sign;
1153   fs = divideSpecials(rhs);
1154
1155   if(category == fcNormal) {
1156     lostFraction lost_fraction = divideSignificand(rhs);
1157     fs = normalize(rounding_mode, lost_fraction);
1158     if(lost_fraction != lfExactlyZero)
1159       fs = (opStatus) (fs | opInexact);
1160   }
1161
1162   return fs;
1163 }
1164
1165 /* Normalized remainder. */
1166 APFloat::opStatus
1167 APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1168 {
1169   opStatus fs;
1170   APFloat V = *this;
1171   unsigned int origSign = sign;
1172   fs = V.divide(rhs, rmNearestTiesToEven);
1173   if (fs == opDivByZero)
1174     return fs;
1175
1176   int parts = partCount();
1177   integerPart *x = new integerPart[parts];
1178   fs = V.convertToInteger(x, parts * integerPartWidth, true,
1179                           rmNearestTiesToEven);
1180   if (fs==opInvalidOp)
1181     return fs;
1182
1183   fs = V.convertFromInteger(x, parts * integerPartWidth, true,
1184                             rmNearestTiesToEven);
1185   assert(fs==opOK);   // should always work
1186
1187   fs = V.multiply(rhs, rounding_mode);
1188   assert(fs==opOK || fs==opInexact);   // should not overflow or underflow
1189
1190   fs = subtract(V, rounding_mode);
1191   assert(fs==opOK || fs==opInexact);   // likewise
1192
1193   if (isZero())
1194     sign = origSign;    // IEEE754 requires this
1195   delete[] x;
1196   return fs;
1197 }
1198
1199 /* Normalized fused-multiply-add.  */
1200 APFloat::opStatus
1201 APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1202                           const APFloat &addend,
1203                           roundingMode rounding_mode)
1204 {
1205   opStatus fs;
1206
1207   /* Post-multiplication sign, before addition.  */
1208   sign ^= multiplicand.sign;
1209
1210   /* If and only if all arguments are normal do we need to do an
1211      extended-precision calculation.  */
1212   if(category == fcNormal
1213      && multiplicand.category == fcNormal
1214      && addend.category == fcNormal) {
1215     lostFraction lost_fraction;
1216
1217     lost_fraction = multiplySignificand(multiplicand, &addend);
1218     fs = normalize(rounding_mode, lost_fraction);
1219     if(lost_fraction != lfExactlyZero)
1220       fs = (opStatus) (fs | opInexact);
1221
1222     /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1223        positive zero unless rounding to minus infinity, except that
1224        adding two like-signed zeroes gives that zero.  */
1225     if(category == fcZero && sign != addend.sign)
1226       sign = (rounding_mode == rmTowardNegative);
1227   } else {
1228     fs = multiplySpecials(multiplicand);
1229
1230     /* FS can only be opOK or opInvalidOp.  There is no more work
1231        to do in the latter case.  The IEEE-754R standard says it is
1232        implementation-defined in this case whether, if ADDEND is a
1233        quiet NaN, we raise invalid op; this implementation does so.
1234
1235        If we need to do the addition we can do so with normal
1236        precision.  */
1237     if(fs == opOK)
1238       fs = addOrSubtract(addend, rounding_mode, false);
1239   }
1240
1241   return fs;
1242 }
1243
1244 /* Comparison requires normalized numbers.  */
1245 APFloat::cmpResult
1246 APFloat::compare(const APFloat &rhs) const
1247 {
1248   cmpResult result;
1249
1250   assert(semantics == rhs.semantics);
1251
1252   switch(convolve(category, rhs.category)) {
1253   default:
1254     assert(0);
1255
1256   case convolve(fcNaN, fcZero):
1257   case convolve(fcNaN, fcNormal):
1258   case convolve(fcNaN, fcInfinity):
1259   case convolve(fcNaN, fcNaN):
1260   case convolve(fcZero, fcNaN):
1261   case convolve(fcNormal, fcNaN):
1262   case convolve(fcInfinity, fcNaN):
1263     return cmpUnordered;
1264
1265   case convolve(fcInfinity, fcNormal):
1266   case convolve(fcInfinity, fcZero):
1267   case convolve(fcNormal, fcZero):
1268     if(sign)
1269       return cmpLessThan;
1270     else
1271       return cmpGreaterThan;
1272
1273   case convolve(fcNormal, fcInfinity):
1274   case convolve(fcZero, fcInfinity):
1275   case convolve(fcZero, fcNormal):
1276     if(rhs.sign)
1277       return cmpGreaterThan;
1278     else
1279       return cmpLessThan;
1280
1281   case convolve(fcInfinity, fcInfinity):
1282     if(sign == rhs.sign)
1283       return cmpEqual;
1284     else if(sign)
1285       return cmpLessThan;
1286     else
1287       return cmpGreaterThan;
1288
1289   case convolve(fcZero, fcZero):
1290     return cmpEqual;
1291
1292   case convolve(fcNormal, fcNormal):
1293     break;
1294   }
1295
1296   /* Two normal numbers.  Do they have the same sign?  */
1297   if(sign != rhs.sign) {
1298     if(sign)
1299       result = cmpLessThan;
1300     else
1301       result = cmpGreaterThan;
1302   } else {
1303     /* Compare absolute values; invert result if negative.  */
1304     result = compareAbsoluteValue(rhs);
1305
1306     if(sign) {
1307       if(result == cmpLessThan)
1308         result = cmpGreaterThan;
1309       else if(result == cmpGreaterThan)
1310         result = cmpLessThan;
1311     }
1312   }
1313
1314   return result;
1315 }
1316
1317 APFloat::opStatus
1318 APFloat::convert(const fltSemantics &toSemantics,
1319                  roundingMode rounding_mode)
1320 {
1321   lostFraction lostFraction;
1322   unsigned int newPartCount, oldPartCount;
1323   opStatus fs;
1324
1325   lostFraction = lfExactlyZero;
1326   newPartCount = partCountForBits(toSemantics.precision + 1);
1327   oldPartCount = partCount();
1328
1329   /* Handle storage complications.  If our new form is wider,
1330      re-allocate our bit pattern into wider storage.  If it is
1331      narrower, we ignore the excess parts, but if narrowing to a
1332      single part we need to free the old storage.
1333      Be careful not to reference significandParts for zeroes
1334      and infinities, since it aborts.  */
1335   if (newPartCount > oldPartCount) {
1336     integerPart *newParts;
1337     newParts = new integerPart[newPartCount];
1338     APInt::tcSet(newParts, 0, newPartCount);
1339     if (category==fcNormal || category==fcNaN)
1340       APInt::tcAssign(newParts, significandParts(), oldPartCount);
1341     freeSignificand();
1342     significand.parts = newParts;
1343   } else if (newPartCount < oldPartCount) {
1344     /* Capture any lost fraction through truncation of parts so we get
1345        correct rounding whilst normalizing.  */
1346     if (category==fcNormal)
1347       lostFraction = lostFractionThroughTruncation
1348         (significandParts(), oldPartCount, toSemantics.precision);
1349     if (newPartCount == 1) {
1350         integerPart newPart = 0;
1351         if (category==fcNormal || category==fcNaN)
1352           newPart = significandParts()[0];
1353         freeSignificand();
1354         significand.part = newPart;
1355     }
1356   }
1357
1358   if(category == fcNormal) {
1359     /* Re-interpret our bit-pattern.  */
1360     exponent += toSemantics.precision - semantics->precision;
1361     semantics = &toSemantics;
1362     fs = normalize(rounding_mode, lostFraction);
1363   } else if (category == fcNaN) {
1364     int shift = toSemantics.precision - semantics->precision;
1365     // No normalization here, just truncate
1366     if (shift>0)
1367       APInt::tcShiftLeft(significandParts(), newPartCount, shift);
1368     else if (shift < 0)
1369       APInt::tcShiftRight(significandParts(), newPartCount, -shift);
1370     // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
1371     // does not give you back the same bits.  This is dubious, and we
1372     // don't currently do it.  You're really supposed to get
1373     // an invalid operation signal at runtime, but nobody does that.
1374     semantics = &toSemantics;
1375     fs = opOK;
1376   } else {
1377     semantics = &toSemantics;
1378     fs = opOK;
1379   }
1380
1381   return fs;
1382 }
1383
1384 /* Convert a floating point number to an integer according to the
1385    rounding mode.  If the rounded integer value is out of range this
1386    returns an invalid operation exception.  If the rounded value is in
1387    range but the floating point number is not the exact integer, the C
1388    standard doesn't require an inexact exception to be raised.  IEEE
1389    854 does require it so we do that.
1390
1391    Note that for conversions to integer type the C standard requires
1392    round-to-zero to always be used.  */
1393 APFloat::opStatus
1394 APFloat::convertToInteger(integerPart *parts, unsigned int width,
1395                           bool isSigned,
1396                           roundingMode rounding_mode) const
1397 {
1398   lostFraction lost_fraction;
1399   unsigned int msb, partsCount;
1400   int bits;
1401
1402   partsCount = partCountForBits(width);
1403
1404   /* Handle the three special cases first.  We produce
1405      a deterministic result even for the Invalid cases. */
1406   if (category == fcNaN) {
1407     // Neither sign nor isSigned affects this.
1408     APInt::tcSet(parts, 0, partsCount);
1409     return opInvalidOp;
1410   }
1411   if (category == fcInfinity) {
1412     if (!sign && isSigned)
1413       APInt::tcSetLeastSignificantBits(parts, partsCount, width-1);
1414     else if (!sign && !isSigned)
1415       APInt::tcSetLeastSignificantBits(parts, partsCount, width);
1416     else if (sign && isSigned) {
1417       APInt::tcSetLeastSignificantBits(parts, partsCount, 1);
1418       APInt::tcShiftLeft(parts, partsCount, width-1);
1419     } else // sign && !isSigned
1420       APInt::tcSet(parts, 0, partsCount);
1421     return opInvalidOp;
1422   }
1423   if (category == fcZero) {
1424     APInt::tcSet(parts, 0, partsCount);
1425     return opOK;
1426   }
1427
1428   /* Shift the bit pattern so the fraction is lost.  */
1429   APFloat tmp(*this);
1430
1431   bits = (int) semantics->precision - 1 - exponent;
1432
1433   if(bits > 0) {
1434     lost_fraction = tmp.shiftSignificandRight(bits);
1435   } else {
1436     if (-bits >= semantics->precision) {
1437       // Unrepresentably large.
1438       if (!sign && isSigned)
1439         APInt::tcSetLeastSignificantBits(parts, partsCount, width-1);
1440       else if (!sign && !isSigned)
1441         APInt::tcSetLeastSignificantBits(parts, partsCount, width);
1442       else if (sign && isSigned) {
1443         APInt::tcSetLeastSignificantBits(parts, partsCount, 1);
1444         APInt::tcShiftLeft(parts, partsCount, width-1);
1445       } else // sign && !isSigned
1446         APInt::tcSet(parts, 0, partsCount);
1447       return (opStatus)(opOverflow | opInexact);
1448     }
1449     tmp.shiftSignificandLeft(-bits);
1450     lost_fraction = lfExactlyZero;
1451   }
1452
1453   if(lost_fraction != lfExactlyZero
1454      && tmp.roundAwayFromZero(rounding_mode, lost_fraction))
1455     tmp.incrementSignificand();
1456
1457   msb = tmp.significandMSB();
1458
1459   /* Negative numbers cannot be represented as unsigned.  */
1460   if(!isSigned && tmp.sign && msb != -1U)
1461     return opInvalidOp;
1462
1463   /* It takes exponent + 1 bits to represent the truncated floating
1464      point number without its sign.  We lose a bit for the sign, but
1465      the maximally negative integer is a special case.  */
1466   if(msb + 1 > width)                /* !! Not same as msb >= width !! */
1467     return opInvalidOp;
1468
1469   if(isSigned && msb + 1 == width
1470      && (!tmp.sign || tmp.significandLSB() != msb))
1471     return opInvalidOp;
1472
1473   APInt::tcAssign(parts, tmp.significandParts(), partsCount);
1474
1475   if(tmp.sign)
1476     APInt::tcNegate(parts, partsCount);
1477
1478   if(lost_fraction == lfExactlyZero)
1479     return opOK;
1480   else
1481     return opInexact;
1482 }
1483
1484 APFloat::opStatus
1485 APFloat::convertFromUnsignedInteger(integerPart *parts,
1486                                     unsigned int partCount,
1487                                     roundingMode rounding_mode)
1488 {
1489   unsigned int msb, precision;
1490   lostFraction lost_fraction;
1491
1492   msb = APInt::tcMSB(parts, partCount) + 1;
1493   precision = semantics->precision;
1494
1495   category = fcNormal;
1496   exponent = precision - 1;
1497
1498   if(msb > precision) {
1499     exponent += (msb - precision);
1500     lost_fraction = shiftRight(parts, partCount, msb - precision);
1501     msb = precision;
1502   } else
1503     lost_fraction = lfExactlyZero;
1504
1505   /* Copy the bit image.  */
1506   zeroSignificand();
1507   APInt::tcAssign(significandParts(), parts, partCountForBits(msb));
1508
1509   return normalize(rounding_mode, lost_fraction);
1510 }
1511
1512 APFloat::opStatus
1513 APFloat::convertFromInteger(const integerPart *parts, unsigned int width,
1514                             bool isSigned, roundingMode rounding_mode)
1515 {
1516   unsigned int partCount = partCountForBits(width);
1517   opStatus status;
1518   APInt api = APInt(width, partCount, parts);
1519   integerPart *copy = new integerPart[partCount];
1520
1521   sign = false;
1522   if(isSigned && APInt::tcExtractBit(parts, width - 1)) {
1523     sign = true;
1524     api = -api;
1525   }
1526
1527   APInt::tcAssign(copy, api.getRawData(), partCount);
1528   status = convertFromUnsignedInteger(copy, partCount, rounding_mode);
1529   return status;
1530 }
1531
1532 APFloat::opStatus
1533 APFloat::convertFromHexadecimalString(const char *p,
1534                                       roundingMode rounding_mode)
1535 {
1536   lostFraction lost_fraction;
1537   integerPart *significand;
1538   unsigned int bitPos, partsCount;
1539   const char *dot, *firstSignificantDigit;
1540
1541   zeroSignificand();
1542   exponent = 0;
1543   category = fcNormal;
1544
1545   significand = significandParts();
1546   partsCount = partCount();
1547   bitPos = partsCount * integerPartWidth;
1548
1549   /* Skip leading zeroes and any(hexa)decimal point.  */
1550   p = skipLeadingZeroesAndAnyDot(p, &dot);
1551   firstSignificantDigit = p;
1552
1553   for(;;) {
1554     integerPart hex_value;
1555
1556     if(*p == '.') {
1557       assert(dot == 0);
1558       dot = p++;
1559     }
1560
1561     hex_value = hexDigitValue(*p);
1562     if(hex_value == -1U) {
1563       lost_fraction = lfExactlyZero;
1564       break;
1565     }
1566
1567     p++;
1568
1569     /* Store the number whilst 4-bit nibbles remain.  */
1570     if(bitPos) {
1571       bitPos -= 4;
1572       hex_value <<= bitPos % integerPartWidth;
1573       significand[bitPos / integerPartWidth] |= hex_value;
1574     } else {
1575       lost_fraction = trailingHexadecimalFraction(p, hex_value);
1576       while(hexDigitValue(*p) != -1U)
1577         p++;
1578       break;
1579     }
1580   }
1581
1582   /* Hex floats require an exponent but not a hexadecimal point.  */
1583   assert(*p == 'p' || *p == 'P');
1584
1585   /* Ignore the exponent if we are zero.  */
1586   if(p != firstSignificantDigit) {
1587     int expAdjustment;
1588
1589     /* Implicit hexadecimal point?  */
1590     if(!dot)
1591       dot = p;
1592
1593     /* Calculate the exponent adjustment implicit in the number of
1594        significant digits.  */
1595     expAdjustment = dot - firstSignificantDigit;
1596     if(expAdjustment < 0)
1597       expAdjustment++;
1598     expAdjustment = expAdjustment * 4 - 1;
1599
1600     /* Adjust for writing the significand starting at the most
1601        significant nibble.  */
1602     expAdjustment += semantics->precision;
1603     expAdjustment -= partsCount * integerPartWidth;
1604
1605     /* Adjust for the given exponent.  */
1606     exponent = totalExponent(p, expAdjustment);
1607   }
1608
1609   return normalize(rounding_mode, lost_fraction);
1610 }
1611
1612 APFloat::opStatus
1613 APFloat::convertFromString(const char *p, roundingMode rounding_mode)
1614 {
1615   /* Handle a leading minus sign.  */
1616   if(*p == '-')
1617     sign = 1, p++;
1618   else
1619     sign = 0;
1620
1621   if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
1622     return convertFromHexadecimalString(p + 2, rounding_mode);
1623
1624   assert(0 && "Decimal to binary conversions not yet implemented");
1625   abort();
1626 }
1627
1628 // For good performance it is desirable for different APFloats
1629 // to produce different integers.
1630 uint32_t
1631 APFloat::getHashValue() const
1632 {
1633   if (category==fcZero) return sign<<8 | semantics->precision ;
1634   else if (category==fcInfinity) return sign<<9 | semantics->precision;
1635   else if (category==fcNaN) return 1<<10 | semantics->precision;
1636   else {
1637     uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
1638     const integerPart* p = significandParts();
1639     for (int i=partCount(); i>0; i--, p++)
1640       hash ^= ((uint32_t)*p) ^ (*p)>>32;
1641     return hash;
1642   }
1643 }
1644
1645 // Conversion from APFloat to/from host float/double.  It may eventually be
1646 // possible to eliminate these and have everybody deal with APFloats, but that
1647 // will take a while.  This approach will not easily extend to long double.
1648 // Current implementation requires integerPartWidth==64, which is correct at
1649 // the moment but could be made more general.
1650
1651 // Denormals have exponent minExponent in APFloat, but minExponent-1 in
1652 // the actual IEEE respresentations.  We compensate for that here.
1653
1654 APInt
1655 APFloat::convertF80LongDoubleAPFloatToAPInt() const
1656 {
1657   assert(semantics == (const llvm::fltSemantics* const)&x87DoubleExtended);
1658   assert (partCount()==2);
1659
1660   uint64_t myexponent, mysignificand;
1661
1662   if (category==fcNormal) {
1663     myexponent = exponent+16383; //bias
1664     mysignificand = significandParts()[0];
1665     if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
1666       myexponent = 0;   // denormal
1667   } else if (category==fcZero) {
1668     myexponent = 0;
1669     mysignificand = 0;
1670   } else if (category==fcInfinity) {
1671     myexponent = 0x7fff;
1672     mysignificand = 0x8000000000000000ULL;
1673   } else if (category==fcNaN) {
1674     myexponent = 0x7fff;
1675     mysignificand = significandParts()[0];
1676   } else
1677     assert(0);
1678
1679   uint64_t words[2];
1680   words[0] =  (((uint64_t)sign & 1) << 63) |
1681               ((myexponent & 0x7fff) <<  48) |
1682               ((mysignificand >>16) & 0xffffffffffffLL);
1683   words[1] = mysignificand & 0xffff;
1684   APInt api(80, 2, words);
1685   return api;
1686 }
1687
1688 APInt
1689 APFloat::convertDoubleAPFloatToAPInt() const
1690 {
1691   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
1692   assert (partCount()==1);
1693
1694   uint64_t myexponent, mysignificand;
1695
1696   if (category==fcNormal) {
1697     myexponent = exponent+1023; //bias
1698     mysignificand = *significandParts();
1699     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
1700       myexponent = 0;   // denormal
1701   } else if (category==fcZero) {
1702     myexponent = 0;
1703     mysignificand = 0;
1704   } else if (category==fcInfinity) {
1705     myexponent = 0x7ff;
1706     mysignificand = 0;
1707   } else if (category==fcNaN) {
1708     myexponent = 0x7ff;
1709     mysignificand = *significandParts();
1710   } else
1711     assert(0);
1712
1713   APInt api(64, (((((uint64_t)sign & 1) << 63) |
1714                  ((myexponent & 0x7ff) <<  52) |
1715                  (mysignificand & 0xfffffffffffffLL))));
1716   return api;
1717 }
1718
1719 APInt
1720 APFloat::convertFloatAPFloatToAPInt() const
1721 {
1722   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
1723   assert (partCount()==1);
1724
1725   uint32_t myexponent, mysignificand;
1726
1727   if (category==fcNormal) {
1728     myexponent = exponent+127; //bias
1729     mysignificand = *significandParts();
1730     if (myexponent == 1 && !(mysignificand & 0x400000))
1731       myexponent = 0;   // denormal
1732   } else if (category==fcZero) {
1733     myexponent = 0;
1734     mysignificand = 0;
1735   } else if (category==fcInfinity) {
1736     myexponent = 0xff;
1737     mysignificand = 0;
1738   } else if (category==fcNaN) {
1739     myexponent = 0xff;
1740     mysignificand = *significandParts();
1741   } else
1742     assert(0);
1743
1744   APInt api(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
1745                  (mysignificand & 0x7fffff)));
1746   return api;
1747 }
1748
1749 APInt
1750 APFloat::convertToAPInt() const
1751 {
1752   if (semantics == (const llvm::fltSemantics* const)&IEEEsingle)
1753     return convertFloatAPFloatToAPInt();
1754   else if (semantics == (const llvm::fltSemantics* const)&IEEEdouble)
1755     return convertDoubleAPFloatToAPInt();
1756   else if (semantics == (const llvm::fltSemantics* const)&x87DoubleExtended)
1757     return convertF80LongDoubleAPFloatToAPInt();
1758
1759   assert(0);
1760   abort();
1761 }
1762
1763 float
1764 APFloat::convertToFloat() const
1765 {
1766   assert(semantics == (const llvm::fltSemantics* const)&IEEEsingle);
1767   APInt api = convertToAPInt();
1768   return api.bitsToFloat();
1769 }
1770
1771 double
1772 APFloat::convertToDouble() const
1773 {
1774   assert(semantics == (const llvm::fltSemantics* const)&IEEEdouble);
1775   APInt api = convertToAPInt();
1776   return api.bitsToDouble();
1777 }
1778
1779 /// Integer bit is explicit in this format.  Current Intel book does not
1780 /// define meaning of:
1781 ///  exponent = all 1's, integer bit not set.
1782 ///  exponent = 0, integer bit set. (formerly "psuedodenormals")
1783 ///  exponent!=0 nor all 1's, integer bit not set. (formerly "unnormals")
1784 void
1785 APFloat::initFromF80LongDoubleAPInt(const APInt &api)
1786 {
1787   assert(api.getBitWidth()==80);
1788   uint64_t i1 = api.getRawData()[0];
1789   uint64_t i2 = api.getRawData()[1];
1790   uint64_t myexponent = (i1 >> 48) & 0x7fff;
1791   uint64_t mysignificand = ((i1 << 16) &  0xffffffffffff0000ULL) |
1792                           (i2 & 0xffff);
1793
1794   initialize(&APFloat::x87DoubleExtended);
1795   assert(partCount()==2);
1796
1797   sign = i1>>63;
1798   if (myexponent==0 && mysignificand==0) {
1799     // exponent, significand meaningless
1800     category = fcZero;
1801   } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
1802     // exponent, significand meaningless
1803     category = fcInfinity;
1804   } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
1805     // exponent meaningless
1806     category = fcNaN;
1807     significandParts()[0] = mysignificand;
1808     significandParts()[1] = 0;
1809   } else {
1810     category = fcNormal;
1811     exponent = myexponent - 16383;
1812     significandParts()[0] = mysignificand;
1813     significandParts()[1] = 0;
1814     if (myexponent==0)          // denormal
1815       exponent = -16382;
1816   }
1817 }
1818
1819 void
1820 APFloat::initFromDoubleAPInt(const APInt &api)
1821 {
1822   assert(api.getBitWidth()==64);
1823   uint64_t i = *api.getRawData();
1824   uint64_t myexponent = (i >> 52) & 0x7ff;
1825   uint64_t mysignificand = i & 0xfffffffffffffLL;
1826
1827   initialize(&APFloat::IEEEdouble);
1828   assert(partCount()==1);
1829
1830   sign = i>>63;
1831   if (myexponent==0 && mysignificand==0) {
1832     // exponent, significand meaningless
1833     category = fcZero;
1834   } else if (myexponent==0x7ff && mysignificand==0) {
1835     // exponent, significand meaningless
1836     category = fcInfinity;
1837   } else if (myexponent==0x7ff && mysignificand!=0) {
1838     // exponent meaningless
1839     category = fcNaN;
1840     *significandParts() = mysignificand;
1841   } else {
1842     category = fcNormal;
1843     exponent = myexponent - 1023;
1844     *significandParts() = mysignificand;
1845     if (myexponent==0)          // denormal
1846       exponent = -1022;
1847     else
1848       *significandParts() |= 0x10000000000000LL;  // integer bit
1849   }
1850 }
1851
1852 void
1853 APFloat::initFromFloatAPInt(const APInt & api)
1854 {
1855   assert(api.getBitWidth()==32);
1856   uint32_t i = (uint32_t)*api.getRawData();
1857   uint32_t myexponent = (i >> 23) & 0xff;
1858   uint32_t mysignificand = i & 0x7fffff;
1859
1860   initialize(&APFloat::IEEEsingle);
1861   assert(partCount()==1);
1862
1863   sign = i >> 31;
1864   if (myexponent==0 && mysignificand==0) {
1865     // exponent, significand meaningless
1866     category = fcZero;
1867   } else if (myexponent==0xff && mysignificand==0) {
1868     // exponent, significand meaningless
1869     category = fcInfinity;
1870   } else if (myexponent==0xff && mysignificand!=0) {
1871     // sign, exponent, significand meaningless
1872     category = fcNaN;
1873     *significandParts() = mysignificand;
1874   } else {
1875     category = fcNormal;
1876     exponent = myexponent - 127;  //bias
1877     *significandParts() = mysignificand;
1878     if (myexponent==0)    // denormal
1879       exponent = -126;
1880     else
1881       *significandParts() |= 0x800000; // integer bit
1882   }
1883 }
1884
1885 /// Treat api as containing the bits of a floating point number.  Currently
1886 /// we infer the floating point type from the size of the APInt.  FIXME: This
1887 /// breaks when we get to PPC128 and IEEE128 (but both cannot exist in the
1888 /// same compile...)
1889 void
1890 APFloat::initFromAPInt(const APInt& api)
1891 {
1892   if (api.getBitWidth() == 32)
1893     return initFromFloatAPInt(api);
1894   else if (api.getBitWidth()==64)
1895     return initFromDoubleAPInt(api);
1896   else if (api.getBitWidth()==80)
1897     return initFromF80LongDoubleAPInt(api);
1898   else
1899     assert(0);
1900 }
1901
1902 APFloat::APFloat(const APInt& api)
1903 {
1904   initFromAPInt(api);
1905 }
1906
1907 APFloat::APFloat(float f)
1908 {
1909   APInt api = APInt(32, 0);
1910   initFromAPInt(api.floatToBits(f));
1911 }
1912
1913 APFloat::APFloat(double d)
1914 {
1915   APInt api = APInt(64, 0);
1916   initFromAPInt(api.doubleToBits(d));
1917 }