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, true };
  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 +
 343                           semantics->implicitIntegerBit ? 1 : 0);
 344 }
 345
 346 unsigned int
 347 APFloat::semanticsPrecision(const fltSemantics &semantics)
 348 {
 349   return semantics.precision;
 350 }
 351
 352 const integerPart *
 353 APFloat::significandParts() const
 354 {
 355   return const_cast<APFloat *>(this)->significandParts();
 356 }
 357
 358 integerPart *
 359 APFloat::significandParts()
 360 {
 361   assert(category == fcNormal || category == fcNaN);
 362
 363   if(partCount() > 1)
 364     return significand.parts;
 365   else
 366     return &significand.part;
 367 }
 368
 369 /* Combine the effect of two lost fractions.  */
 370 lostFraction
 371 APFloat::combineLostFractions(lostFraction moreSignificant,
 372                               lostFraction lessSignificant)
 373 {
 374   if(lessSignificant != lfExactlyZero) {
 375     if(moreSignificant == lfExactlyZero)
 376       moreSignificant = lfLessThanHalf;
 377     else if(moreSignificant == lfExactlyHalf)
 378       moreSignificant = lfMoreThanHalf;
 379   }
 380
 381   return moreSignificant;
 382 }
 383
 384 void
 385 APFloat::zeroSignificand()
 386 {
 387   category = fcNormal;
 388   APInt::tcSet(significandParts(), 0, partCount());
 389 }
 390
 391 /* Increment an fcNormal floating point number's significand.  */
 392 void
 393 APFloat::incrementSignificand()
 394 {
 395   integerPart carry;
 396
 397   carry = APInt::tcIncrement(significandParts(), partCount());
 398
 399   /* Our callers should never cause us to overflow.  */
 400   assert(carry == 0);
 401 }
 402
 403 /* Add the significand of the RHS.  Returns the carry flag.  */
 404 integerPart
 405 APFloat::addSignificand(const APFloat &rhs)
 406 {
 407   integerPart *parts;
 408
 409   parts = significandParts();
 410
 411   assert(semantics == rhs.semantics);
 412   assert(exponent == rhs.exponent);
 413
 414   return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
 415 }
 416
 417 /* Subtract the significand of the RHS with a borrow flag.  Returns
 418    the borrow flag.  */
 419 integerPart
 420 APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
 421 {
 422   integerPart *parts;
 423
 424   parts = significandParts();
 425
 426   assert(semantics == rhs.semantics);
 427   assert(exponent == rhs.exponent);
 428
 429   return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
 430                            partCount());
 431 }
 432
 433 /* Multiply the significand of the RHS.  If ADDEND is non-NULL, add it
 434    on to the full-precision result of the multiplication.  Returns the
 435    lost fraction.  */
 436 lostFraction
 437 APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
 438 {
 439   unsigned int omsb;    // One, not zero, based MSB.
 440   unsigned int partsCount, newPartsCount, precision;
 441   integerPart *lhsSignificand;
 442   integerPart scratch[4];
 443   integerPart *fullSignificand;
 444   lostFraction lost_fraction;
 445
 446   assert(semantics == rhs.semantics);
 447
 448   precision = semantics->precision;
 449   newPartsCount = partCountForBits(precision * 2);
 450
 451   if(newPartsCount > 4)
 452     fullSignificand = new integerPart[newPartsCount];
 453   else
 454     fullSignificand = scratch;
 455
 456   lhsSignificand = significandParts();
 457   partsCount = partCount();
 458
 459   APInt::tcFullMultiply(fullSignificand, lhsSignificand,
 460                         rhs.significandParts(), partsCount);
 461
 462   lost_fraction = lfExactlyZero;
 463   omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
 464   exponent += rhs.exponent;
 465
 466   if(addend) {
 467     Significand savedSignificand = significand;
 468     const fltSemantics *savedSemantics = semantics;
 469     fltSemantics extendedSemantics;
 470     opStatus status;
 471     unsigned int extendedPrecision;
 472
 473     /* Normalize our MSB.  */
 474     extendedPrecision = precision + precision - 1;
 475     if(omsb != extendedPrecision)
 476       {
 477         APInt::tcShiftLeft(fullSignificand, newPartsCount,
 478                            extendedPrecision - omsb);
 479         exponent -= extendedPrecision - omsb;
 480       }
 481
 482     /* Create new semantics.  */
 483     extendedSemantics = *semantics;
 484     extendedSemantics.precision = extendedPrecision;
 485
 486     if(newPartsCount == 1)
 487       significand.part = fullSignificand[0];
 488     else
 489       significand.parts = fullSignificand;
 490     semantics = &extendedSemantics;
 491
 492     APFloat extendedAddend(*addend);
 493     status = extendedAddend.convert(extendedSemantics, rmTowardZero);
 494     assert(status == opOK);
 495     lost_fraction = addOrSubtractSignificand(extendedAddend, false);
 496
 497     /* Restore our state.  */
 498     if(newPartsCount == 1)
 499       fullSignificand[0] = significand.part;
 500     significand = savedSignificand;
 501     semantics = savedSemantics;
 502
 503     omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
 504   }
 505
 506   exponent -= (precision - 1);
 507
 508   if(omsb > precision) {
 509     unsigned int bits, significantParts;
 510     lostFraction lf;
 511
 512     bits = omsb - precision;
 513     significantParts = partCountForBits(omsb);
 514     lf = shiftRight(fullSignificand, significantParts, bits);
 515     lost_fraction = combineLostFractions(lf, lost_fraction);
 516     exponent += bits;
 517   }
 518
 519   APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
 520
 521   if(newPartsCount > 4)
 522     delete [] fullSignificand;
 523
 524   return lost_fraction;
 525 }
 526
 527 /* Multiply the significands of LHS and RHS to DST.  */
 528 lostFraction
 529 APFloat::divideSignificand(const APFloat &rhs)
 530 {
 531   unsigned int bit, i, partsCount;
 532   const integerPart *rhsSignificand;
 533   integerPart *lhsSignificand, *dividend, *divisor;
 534   integerPart scratch[4];
 535   lostFraction lost_fraction;
 536
 537   assert(semantics == rhs.semantics);
 538
 539   lhsSignificand = significandParts();
 540   rhsSignificand = rhs.significandParts();
 541   partsCount = partCount();
 542
 543   if(partsCount > 2)
 544     dividend = new integerPart[partsCount * 2];
 545   else
 546     dividend = scratch;
 547
 548   divisor = dividend + partsCount;
 549
 550   /* Copy the dividend and divisor as they will be modified in-place.  */
 551   for(i = 0; i < partsCount; i++) {
 552     dividend[i] = lhsSignificand[i];
 553     divisor[i] = rhsSignificand[i];
 554     lhsSignificand[i] = 0;
 555   }
 556
 557   exponent -= rhs.exponent;
 558
 559   unsigned int precision = semantics->precision;
 560
 561   /* Normalize the divisor.  */
 562   bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
 563   if(bit) {
 564     exponent += bit;
 565     APInt::tcShiftLeft(divisor, partsCount, bit);
 566   }
 567
 568   /* Normalize the dividend.  */
 569   bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
 570   if(bit) {
 571     exponent -= bit;
 572     APInt::tcShiftLeft(dividend, partsCount, bit);
 573   }
 574
 575   if(APInt::tcCompare(dividend, divisor, partsCount) < 0) {
 576     exponent--;
 577     APInt::tcShiftLeft(dividend, partsCount, 1);
 578     assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
 579   }
 580
 581   /* Long division.  */
 582   for(bit = precision; bit; bit -= 1) {
 583     if(APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
 584       APInt::tcSubtract(dividend, divisor, 0, partsCount);
 585       APInt::tcSetBit(lhsSignificand, bit - 1);
 586     }
 587
 588     APInt::tcShiftLeft(dividend, partsCount, 1);
 589   }
 590
 591   /* Figure out the lost fraction.  */
 592   int cmp = APInt::tcCompare(dividend, divisor, partsCount);
 593
 594   if(cmp > 0)
 595     lost_fraction = lfMoreThanHalf;
 596   else if(cmp == 0)
 597     lost_fraction = lfExactlyHalf;
 598   else if(APInt::tcIsZero(dividend, partsCount))
 599     lost_fraction = lfExactlyZero;
 600   else
 601     lost_fraction = lfLessThanHalf;
 602
 603   if(partsCount > 2)
 604     delete [] dividend;
 605
 606   return lost_fraction;
 607 }
 608
 609 unsigned int
 610 APFloat::significandMSB() const
 611 {
 612   return APInt::tcMSB(significandParts(), partCount());
 613 }
 614
 615 unsigned int
 616 APFloat::significandLSB() const
 617 {
 618   return APInt::tcLSB(significandParts(), partCount());
 619 }
 620
 621 /* Note that a zero result is NOT normalized to fcZero.  */
 622 lostFraction
 623 APFloat::shiftSignificandRight(unsigned int bits)
 624 {
 625   /* Our exponent should not overflow.  */
 626   assert((exponent_t) (exponent + bits) >= exponent);
 627
 628   exponent += bits;
 629
 630   return shiftRight(significandParts(), partCount(), bits);
 631 }
 632
 633 /* Shift the significand left BITS bits, subtract BITS from its exponent.  */
 634 void
 635 APFloat::shiftSignificandLeft(unsigned int bits)
 636 {
 637   assert(bits < semantics->precision);
 638
 639   if(bits) {
 640     unsigned int partsCount = partCount();
 641
 642     APInt::tcShiftLeft(significandParts(), partsCount, bits);
 643     exponent -= bits;
 644
 645     assert(!APInt::tcIsZero(significandParts(), partsCount));
 646   }
 647 }
 648
 649 APFloat::cmpResult
 650 APFloat::compareAbsoluteValue(const APFloat &rhs) const
 651 {
 652   int compare;
 653
 654   assert(semantics == rhs.semantics);
 655   assert(category == fcNormal);
 656   assert(rhs.category == fcNormal);
 657
 658   compare = exponent - rhs.exponent;
 659
 660   /* If exponents are equal, do an unsigned bignum comparison of the
 661      significands.  */
 662   if(compare == 0)
 663     compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
 664                                partCount());
 665
 666   if(compare > 0)
 667     return cmpGreaterThan;
 668   else if(compare < 0)
 669     return cmpLessThan;
 670   else
 671     return cmpEqual;
 672 }
 673
 674 /* Handle overflow.  Sign is preserved.  We either become infinity or
 675    the largest finite number.  */
 676 APFloat::opStatus
 677 APFloat::handleOverflow(roundingMode rounding_mode)
 678 {
 679   /* Infinity?  */
 680   if(rounding_mode == rmNearestTiesToEven
 681      || rounding_mode == rmNearestTiesToAway
 682      || (rounding_mode == rmTowardPositive && !sign)
 683      || (rounding_mode == rmTowardNegative && sign))
 684     {
 685       category = fcInfinity;
 686       return (opStatus) (opOverflow | opInexact);
 687     }
 688
 689   /* Otherwise we become the largest finite number.  */
 690   category = fcNormal;
 691   exponent = semantics->maxExponent;
 692   APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
 693                                    semantics->precision);
 694
 695   return opInexact;
 696 }
 697
 698 /* This routine must work for fcZero of both signs, and fcNormal
 699    numbers.  */
 700 bool
 701 APFloat::roundAwayFromZero(roundingMode rounding_mode,
 702                            lostFraction lost_fraction)
 703 {
 704   /* NaNs and infinities should not have lost fractions.  */
 705   assert(category == fcNormal || category == fcZero);
 706
 707   /* Our caller has already handled this case.  */
 708   assert(lost_fraction != lfExactlyZero);
 709
 710   switch(rounding_mode) {
 711   default:
 712     assert(0);
 713
 714   case rmNearestTiesToAway:
 715     return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
 716
 717   case rmNearestTiesToEven:
 718     if(lost_fraction == lfMoreThanHalf)
 719       return true;
 720
 721     /* Our zeroes don't have a significand to test.  */
 722     if(lost_fraction == lfExactlyHalf && category != fcZero)
 723       return significandParts()[0] & 1;
 724
 725     return false;
 726
 727   case rmTowardZero:
 728     return false;
 729
 730   case rmTowardPositive:
 731     return sign == false;
 732
 733   case rmTowardNegative:
 734     return sign == true;
 735   }
 736 }
 737
 738 APFloat::opStatus
 739 APFloat::normalize(roundingMode rounding_mode,
 740                    lostFraction lost_fraction)
 741 {
 742   unsigned int omsb;            /* One, not zero, based MSB.  */
 743   int exponentChange;
 744
 745   if(category != fcNormal)
 746     return opOK;
 747
 748   /* Before rounding normalize the exponent of fcNormal numbers.  */
 749   omsb = significandMSB() + 1;
 750
 751   if(omsb) {
 752     /* OMSB is numbered from 1.  We want to place it in the integer
 753        bit numbered PRECISON if possible, with a compensating change in
 754        the exponent.  */
 755     exponentChange = omsb - semantics->precision;
 756
 757     /* If the resulting exponent is too high, overflow according to
 758        the rounding mode.  */
 759     if(exponent + exponentChange > semantics->maxExponent)
 760       return handleOverflow(rounding_mode);
 761
 762     /* Subnormal numbers have exponent minExponent, and their MSB
 763        is forced based on that.  */
 764     if(exponent + exponentChange < semantics->minExponent)
 765       exponentChange = semantics->minExponent - exponent;
 766
 767     /* Shifting left is easy as we don't lose precision.  */
 768     if(exponentChange < 0) {
 769       assert(lost_fraction == lfExactlyZero);
 770
 771       shiftSignificandLeft(-exponentChange);
 772
 773       return opOK;
 774     }
 775
 776     if(exponentChange > 0) {
 777       lostFraction lf;
 778
 779       /* Shift right and capture any new lost fraction.  */
 780       lf = shiftSignificandRight(exponentChange);
 781
 782       lost_fraction = combineLostFractions(lf, lost_fraction);
 783
 784       /* Keep OMSB up-to-date.  */
 785       if(omsb > (unsigned) exponentChange)
 786         omsb -= (unsigned) exponentChange;
 787       else
 788         omsb = 0;
 789     }
 790   }
 791
 792   /* Now round the number according to rounding_mode given the lost
 793      fraction.  */
 794
 795   /* As specified in IEEE 754, since we do not trap we do not report
 796      underflow for exact results.  */
 797   if(lost_fraction == lfExactlyZero) {
 798     /* Canonicalize zeroes.  */
 799     if(omsb == 0)
 800       category = fcZero;
 801
 802     return opOK;
 803   }
 804
 805   /* Increment the significand if we're rounding away from zero.  */
 806   if(roundAwayFromZero(rounding_mode, lost_fraction)) {
 807     if(omsb == 0)
 808       exponent = semantics->minExponent;
 809
 810     incrementSignificand();
 811     omsb = significandMSB() + 1;
 812
 813     /* Did the significand increment overflow?  */
 814     if(omsb == (unsigned) semantics->precision + 1) {
 815       /* Renormalize by incrementing the exponent and shifting our
 816          significand right one.  However if we already have the
 817          maximum exponent we overflow to infinity.  */
 818       if(exponent == semantics->maxExponent) {
 819         category = fcInfinity;
 820
 821         return (opStatus) (opOverflow | opInexact);
 822       }
 823
 824       shiftSignificandRight(1);
 825
 826       return opInexact;
 827     }
 828   }
 829
 830   /* The normal case - we were and are not denormal, and any
 831      significand increment above didn't overflow.  */
 832   if(omsb == semantics->precision)
 833     return opInexact;
 834
 835   /* We have a non-zero denormal.  */
 836   assert(omsb < semantics->precision);
 837   assert(exponent == semantics->minExponent);
 838
 839   /* Canonicalize zeroes.  */
 840   if(omsb == 0)
 841     category = fcZero;
 842
 843   /* The fcZero case is a denormal that underflowed to zero.  */
 844   return (opStatus) (opUnderflow | opInexact);
 845 }
 846
 847 APFloat::opStatus
 848 APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
 849 {
 850   switch(convolve(category, rhs.category)) {
 851   default:
 852     assert(0);
 853
 854   case convolve(fcNaN, fcZero):
 855   case convolve(fcNaN, fcNormal):
 856   case convolve(fcNaN, fcInfinity):
 857   case convolve(fcNaN, fcNaN):
 858   case convolve(fcNormal, fcZero):
 859   case convolve(fcInfinity, fcNormal):
 860   case convolve(fcInfinity, fcZero):
 861     return opOK;
 862
 863   case convolve(fcZero, fcNaN):
 864   case convolve(fcNormal, fcNaN):
 865   case convolve(fcInfinity, fcNaN):
 866     category = fcNaN;
 867     copySignificand(rhs);
 868     return opOK;
 869
 870   case convolve(fcNormal, fcInfinity):
 871   case convolve(fcZero, fcInfinity):
 872     category = fcInfinity;
 873     sign = rhs.sign ^ subtract;
 874     return opOK;
 875
 876   case convolve(fcZero, fcNormal):
 877     assign(rhs);
 878     sign = rhs.sign ^ subtract;
 879     return opOK;
 880
 881   case convolve(fcZero, fcZero):
 882     /* Sign depends on rounding mode; handled by caller.  */
 883     return opOK;
 884
 885   case convolve(fcInfinity, fcInfinity):
 886     /* Differently signed infinities can only be validly
 887        subtracted.  */
 888     if(sign ^ rhs.sign != subtract) {
 889       category = fcNaN;
 890       // Arbitrary but deterministic value for significand
 891       APInt::tcSet(significandParts(), ~0U, partCount());
 892       return opInvalidOp;
 893     }
 894
 895     return opOK;
 896
 897   case convolve(fcNormal, fcNormal):
 898     return opDivByZero;
 899   }
 900 }
 901
 902 /* Add or subtract two normal numbers.  */
 903 lostFraction
 904 APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
 905 {
 906   integerPart carry;
 907   lostFraction lost_fraction;
 908   int bits;
 909
 910   /* Determine if the operation on the absolute values is effectively
 911      an addition or subtraction.  */
 912   subtract ^= (sign ^ rhs.sign);
 913
 914   /* Are we bigger exponent-wise than the RHS?  */
 915   bits = exponent - rhs.exponent;
 916
 917   /* Subtraction is more subtle than one might naively expect.  */
 918   if(subtract) {
 919     APFloat temp_rhs(rhs);
 920     bool reverse;
 921
 922     if (bits == 0) {
 923       reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
 924       lost_fraction = lfExactlyZero;
 925     } else if (bits > 0) {
 926       lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
 927       shiftSignificandLeft(1);
 928       reverse = false;
 929     } else {
 930       lost_fraction = shiftSignificandRight(-bits - 1);
 931       temp_rhs.shiftSignificandLeft(1);
 932       reverse = true;
 933     }
 934
 935     if (reverse) {
 936       carry = temp_rhs.subtractSignificand
 937         (*this, lost_fraction != lfExactlyZero);
 938       copySignificand(temp_rhs);
 939       sign = !sign;
 940     } else {
 941       carry = subtractSignificand
 942         (temp_rhs, lost_fraction != lfExactlyZero);
 943     }
 944
 945     /* Invert the lost fraction - it was on the RHS and
 946        subtracted.  */
 947     if(lost_fraction == lfLessThanHalf)
 948       lost_fraction = lfMoreThanHalf;
 949     else if(lost_fraction == lfMoreThanHalf)
 950       lost_fraction = lfLessThanHalf;
 951
 952     /* The code above is intended to ensure that no borrow is
 953        necessary.  */
 954     assert(!carry);
 955   } else {
 956     if(bits > 0) {
 957       APFloat temp_rhs(rhs);
 958
 959       lost_fraction = temp_rhs.shiftSignificandRight(bits);
 960       carry = addSignificand(temp_rhs);
 961     } else {
 962       lost_fraction = shiftSignificandRight(-bits);
 963       carry = addSignificand(rhs);
 964     }
 965
 966     /* We have a guard bit; generating a carry cannot happen.  */
 967     assert(!carry);
 968   }
 969
 970   return lost_fraction;
 971 }
 972
 973 APFloat::opStatus
 974 APFloat::multiplySpecials(const APFloat &rhs)
 975 {
 976   switch(convolve(category, rhs.category)) {
 977   default:
 978     assert(0);
 979
 980   case convolve(fcNaN, fcZero):
 981   case convolve(fcNaN, fcNormal):
 982   case convolve(fcNaN, fcInfinity):
 983   case convolve(fcNaN, fcNaN):
 984     return opOK;
 985
 986   case convolve(fcZero, fcNaN):
 987   case convolve(fcNormal, fcNaN):
 988   case convolve(fcInfinity, fcNaN):
 989     category = fcNaN;
 990     copySignificand(rhs);
 991     return opOK;
 992
 993   case convolve(fcNormal, fcInfinity):
 994   case convolve(fcInfinity, fcNormal):
 995   case convolve(fcInfinity, fcInfinity):
 996     category = fcInfinity;
 997     return opOK;
 998
 999   case convolve(fcZero, fcNormal):
1000   case convolve(fcNormal, fcZero):
1001   case convolve(fcZero, fcZero):
1002     category = fcZero;
1003     return opOK;
1004
1005   case convolve(fcZero, fcInfinity):
1006   case convolve(fcInfinity, fcZero):
1007     category = fcNaN;
1008     // Arbitrary but deterministic value for significand
1009     APInt::tcSet(significandParts(), ~0U, partCount());
1010     return opInvalidOp;
1011
1012   case convolve(fcNormal, fcNormal):
1013     return opOK;
1014   }
1015 }
1016
1017 APFloat::opStatus
1018 APFloat::divideSpecials(const APFloat &rhs)
1019 {
1020   switch(convolve(category, rhs.category)) {
1021   default:
1022     assert(0);
1023
1024   case convolve(fcNaN, fcZero):
1025   case convolve(fcNaN, fcNormal):
1026   case convolve(fcNaN, fcInfinity):
1027   case convolve(fcNaN, fcNaN):
1028   case convolve(fcInfinity, fcZero):
1029   case convolve(fcInfinity, fcNormal):
1030   case convolve(fcZero, fcInfinity):
1031   case convolve(fcZero, fcNormal):
1032     return opOK;
1033
1034   case convolve(fcZero, fcNaN):
1035   case convolve(fcNormal, fcNaN):
1036   case convolve(fcInfinity, fcNaN):
1037     category = fcNaN;
1038     copySignificand(rhs);
1039     return opOK;
1040
1041   case convolve(fcNormal, fcInfinity):
1042     category = fcZero;
1043     return opOK;
1044
1045   case convolve(fcNormal, fcZero):
1046     category = fcInfinity;
1047     return opDivByZero;
1048
1049   case convolve(fcInfinity, fcInfinity):
1050   case convolve(fcZero, fcZero):
1051     category = fcNaN;
1052     // Arbitrary but deterministic value for significand
1053     APInt::tcSet(significandParts(), ~0U, partCount());
1054     return opInvalidOp;
1055
1056   case convolve(fcNormal, fcNormal):
1057     return opOK;
1058   }
1059 }
1060
1061 /* Change sign.  */
1062 void
1063 APFloat::changeSign()
1064 {
1065   /* Look mummy, this one's easy.  */
1066   sign = !sign;
1067 }
1068
1069 void
1070 APFloat::clearSign()
1071 {
1072   /* So is this one. */
1073   sign = 0;
1074 }
1075
1076 void
1077 APFloat::copySign(const APFloat &rhs)
1078 {
1079   /* And this one. */
1080   sign = rhs.sign;
1081 }
1082
1083 /* Normalized addition or subtraction.  */
1084 APFloat::opStatus
1085 APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1086                        bool subtract)
1087 {
1088   opStatus fs;
1089
1090   fs = addOrSubtractSpecials(rhs, subtract);
1091
1092   /* This return code means it was not a simple case.  */
1093   if(fs == opDivByZero) {
1094     lostFraction lost_fraction;
1095
1096     lost_fraction = addOrSubtractSignificand(rhs, subtract);
1097     fs = normalize(rounding_mode, lost_fraction);
1098
1099     /* Can only be zero if we lost no fraction.  */
1100     assert(category != fcZero || lost_fraction == lfExactlyZero);
1101   }
1102
1103   /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1104      positive zero unless rounding to minus infinity, except that
1105      adding two like-signed zeroes gives that zero.  */
1106   if(category == fcZero) {
1107     if(rhs.category != fcZero || (sign == rhs.sign) == subtract)
1108       sign = (rounding_mode == rmTowardNegative);
1109   }
1110
1111   return fs;
1112 }
1113
1114 /* Normalized addition.  */
1115 APFloat::opStatus
1116 APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1117 {
1118   return addOrSubtract(rhs, rounding_mode, false);
1119 }
1120
1121 /* Normalized subtraction.  */
1122 APFloat::opStatus
1123 APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1124 {
1125   return addOrSubtract(rhs, rounding_mode, true);
1126 }
1127
1128 /* Normalized multiply.  */
1129 APFloat::opStatus
1130 APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1131 {
1132   opStatus fs;
1133
1134   sign ^= rhs.sign;
1135   fs = multiplySpecials(rhs);
1136
1137   if(category == fcNormal) {
1138     lostFraction lost_fraction = multiplySignificand(rhs, 0);
1139     fs = normalize(rounding_mode, lost_fraction);
1140     if(lost_fraction != lfExactlyZero)
1141       fs = (opStatus) (fs | opInexact);
1142   }
1143
1144   return fs;
1145 }
1146
1147 /* Normalized divide.  */
1148 APFloat::opStatus
1149 APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1150 {
1151   opStatus fs;
1152
1153   sign ^= rhs.sign;
1154   fs = divideSpecials(rhs);
1155
1156   if(category == fcNormal) {
1157     lostFraction lost_fraction = divideSignificand(rhs);
1158     fs = normalize(rounding_mode, lost_fraction);
1159     if(lost_fraction != lfExactlyZero)
1160       fs = (opStatus) (fs | opInexact);
1161   }
1162
1163   return fs;
1164 }
1165
1166 /* Normalized remainder. */
1167 APFloat::opStatus
1168 APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1169 {
1170   opStatus fs;
1171   APFloat V = *this;
1172   unsigned int origSign = sign;
1173   fs = V.divide(rhs, rmNearestTiesToEven);
1174   if (fs == opDivByZero)
1175     return fs;
1176
1177   int parts = partCount();
1178   integerPart *x = new integerPart[parts];
1179   fs = V.convertToInteger(x, parts * integerPartWidth, true,
1180                           rmNearestTiesToEven);
1181   if (fs==opInvalidOp)
1182     return fs;
1183
1184   fs = V.convertFromInteger(x, parts, true, 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   unsigned int newPartCount;
1322   opStatus fs;
1323
1324   newPartCount = partCountForBits(toSemantics.precision + 1);
1325
1326   /* If our new form is wider, re-allocate our bit pattern into wider
1327      storage.  */
1328   if(newPartCount > partCount()) {
1329     integerPart *newParts;
1330
1331     newParts = new integerPart[newPartCount];
1332     APInt::tcSet(newParts, 0, newPartCount);
1333     APInt::tcAssign(newParts, significandParts(), partCount());
1334     freeSignificand();
1335     significand.parts = newParts;
1336   }
1337
1338   if(category == fcNormal) {
1339     /* Re-interpret our bit-pattern.  */
1340     exponent += toSemantics.precision - semantics->precision;
1341     semantics = &toSemantics;
1342     fs = normalize(rounding_mode, lfExactlyZero);
1343   } else {
1344     semantics = &toSemantics;
1345     fs = opOK;
1346   }
1347
1348   return fs;
1349 }
1350
1351 /* Convert a floating point number to an integer according to the
1352    rounding mode.  If the rounded integer value is out of range this
1353    returns an invalid operation exception.  If the rounded value is in
1354    range but the floating point number is not the exact integer, the C
1355    standard doesn't require an inexact exception to be raised.  IEEE
1356    854 does require it so we do that.
1357
1358    Note that for conversions to integer type the C standard requires
1359    round-to-zero to always be used.  */
1360 APFloat::opStatus
1361 APFloat::convertToInteger(integerPart *parts, unsigned int width,
1362                           bool isSigned,
1363                           roundingMode rounding_mode) const
1364 {
1365   lostFraction lost_fraction;
1366   unsigned int msb, partsCount;
1367   int bits;
1368
1369   /* Handle the three special cases first.  */
1370   if(category == fcInfinity || category == fcNaN)
1371     return opInvalidOp;
1372
1373   partsCount = partCountForBits(width);
1374
1375   if(category == fcZero) {
1376     APInt::tcSet(parts, 0, partsCount);
1377     return opOK;
1378   }
1379
1380   /* Shift the bit pattern so the fraction is lost.  */
1381   APFloat tmp(*this);
1382
1383   bits = (int) semantics->precision - 1 - exponent;
1384
1385   if(bits > 0) {
1386     lost_fraction = tmp.shiftSignificandRight(bits);
1387   } else {
1388     tmp.shiftSignificandLeft(-bits);
1389     lost_fraction = lfExactlyZero;
1390   }
1391
1392   if(lost_fraction != lfExactlyZero
1393      && tmp.roundAwayFromZero(rounding_mode, lost_fraction))
1394     tmp.incrementSignificand();
1395
1396   msb = tmp.significandMSB();
1397
1398   /* Negative numbers cannot be represented as unsigned.  */
1399   if(!isSigned && tmp.sign && msb != -1U)
1400     return opInvalidOp;
1401
1402   /* It takes exponent + 1 bits to represent the truncated floating
1403      point number without its sign.  We lose a bit for the sign, but
1404      the maximally negative integer is a special case.  */
1405   if(msb + 1 > width)           /* !! Not same as msb >= width !! */
1406     return opInvalidOp;
1407
1408   if(isSigned && msb + 1 == width
1409      && (!tmp.sign || tmp.significandLSB() != msb))
1410     return opInvalidOp;
1411
1412   APInt::tcAssign(parts, tmp.significandParts(), partsCount);
1413
1414   if(tmp.sign)
1415     APInt::tcNegate(parts, partsCount);
1416
1417   if(lost_fraction == lfExactlyZero)
1418     return opOK;
1419   else
1420     return opInexact;
1421 }
1422
1423 APFloat::opStatus
1424 APFloat::convertFromUnsignedInteger(integerPart *parts,
1425                                     unsigned int partCount,
1426                                     roundingMode rounding_mode)
1427 {
1428   unsigned int msb, precision;
1429   lostFraction lost_fraction;
1430
1431   msb = APInt::tcMSB(parts, partCount) + 1;
1432   precision = semantics->precision;
1433
1434   category = fcNormal;
1435   exponent = precision - 1;
1436
1437   if(msb > precision) {
1438     exponent += (msb - precision);
1439     lost_fraction = shiftRight(parts, partCount, msb - precision);
1440     msb = precision;
1441   } else
1442     lost_fraction = lfExactlyZero;
1443
1444   /* Copy the bit image.  */
1445   zeroSignificand();
1446   APInt::tcAssign(significandParts(), parts, partCountForBits(msb));
1447
1448   return normalize(rounding_mode, lost_fraction);
1449 }
1450
1451 APFloat::opStatus
1452 APFloat::convertFromInteger(const integerPart *parts,
1453                             unsigned int partCount, bool isSigned,
1454                             roundingMode rounding_mode)
1455 {
1456   unsigned int width;
1457   opStatus status;
1458   integerPart *copy;
1459
1460   copy = new integerPart[partCount];
1461   APInt::tcAssign(copy, parts, partCount);
1462
1463   width = partCount * integerPartWidth;
1464
1465   sign = false;
1466   if(isSigned && APInt::tcExtractBit(parts, width - 1)) {
1467     sign = true;
1468     APInt::tcNegate(copy, partCount);
1469   }
1470
1471   status = convertFromUnsignedInteger(copy, partCount, rounding_mode);
1472   delete [] copy;
1473
1474   return status;
1475 }
1476
1477 APFloat::opStatus
1478 APFloat::convertFromHexadecimalString(const char *p,
1479                                       roundingMode rounding_mode)
1480 {
1481   lostFraction lost_fraction;
1482   integerPart *significand;
1483   unsigned int bitPos, partsCount;
1484   const char *dot, *firstSignificantDigit;
1485
1486   zeroSignificand();
1487   exponent = 0;
1488   category = fcNormal;
1489
1490   significand = significandParts();
1491   partsCount = partCount();
1492   bitPos = partsCount * integerPartWidth;
1493
1494   /* Skip leading zeroes and any(hexa)decimal point.  */
1495   p = skipLeadingZeroesAndAnyDot(p, &dot);
1496   firstSignificantDigit = p;
1497
1498   for(;;) {
1499     integerPart hex_value;
1500
1501     if(*p == '.') {
1502       assert(dot == 0);
1503       dot = p++;
1504     }
1505
1506     hex_value = hexDigitValue(*p);
1507     if(hex_value == -1U) {
1508       lost_fraction = lfExactlyZero;
1509       break;
1510     }
1511
1512     p++;
1513
1514     /* Store the number whilst 4-bit nibbles remain.  */
1515     if(bitPos) {
1516       bitPos -= 4;
1517       hex_value <<= bitPos % integerPartWidth;
1518       significand[bitPos / integerPartWidth] |= hex_value;
1519     } else {
1520       lost_fraction = trailingHexadecimalFraction(p, hex_value);
1521       while(hexDigitValue(*p) != -1U)
1522         p++;
1523       break;
1524     }
1525   }
1526
1527   /* Hex floats require an exponent but not a hexadecimal point.  */
1528   assert(*p == 'p' || *p == 'P');
1529
1530   /* Ignore the exponent if we are zero.  */
1531   if(p != firstSignificantDigit) {
1532     int expAdjustment;
1533
1534     /* Implicit hexadecimal point?  */
1535     if(!dot)
1536       dot = p;
1537
1538     /* Calculate the exponent adjustment implicit in the number of
1539        significant digits.  */
1540     expAdjustment = dot - firstSignificantDigit;
1541     if(expAdjustment < 0)
1542       expAdjustment++;
1543     expAdjustment = expAdjustment * 4 - 1;
1544
1545     /* Adjust for writing the significand starting at the most
1546        significant nibble.  */
1547     expAdjustment += semantics->precision;
1548     expAdjustment -= partsCount * integerPartWidth;
1549
1550     /* Adjust for the given exponent.  */
1551     exponent = totalExponent(p, expAdjustment);
1552   }
1553
1554   return normalize(rounding_mode, lost_fraction);
1555 }
1556
1557 APFloat::opStatus
1558 APFloat::convertFromString(const char *p, roundingMode rounding_mode) {
1559   /* Handle a leading minus sign.  */
1560   if(*p == '-')
1561     sign = 1, p++;
1562   else
1563     sign = 0;
1564
1565   if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
1566     return convertFromHexadecimalString(p + 2, rounding_mode);
1567
1568   assert(0 && "Decimal to binary conversions not yet implemented");
1569   abort();
1570 }
1571
1572 // For good performance it is desirable for different APFloats
1573 // to produce different integers.
1574 uint32_t
1575 APFloat::getHashValue() const {
1576   if (category==fcZero) return sign<<8 | semantics->precision ;
1577   else if (category==fcInfinity) return sign<<9 | semantics->precision;
1578   else if (category==fcNaN) return 1<<10 | semantics->precision;
1579   else {
1580     uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
1581     const integerPart* p = significandParts();
1582     for (int i=partCount(); i>0; i--, p++)
1583       hash ^= ((uint32_t)*p) ^ (*p)>>32;
1584     return hash;
1585   }
1586 }
1587
1588 // Conversion from APFloat to/from host float/double.  It may eventually be
1589 // possible to eliminate these and have everybody deal with APFloats, but that
1590 // will take a while.  This approach will not easily extend to long double.
1591 // Current implementation requires partCount()==1, which is correct at the
1592 // moment but could be made more general.
1593
1594 // Denormals have exponent minExponent in APFloat, but minExponent-1 in
1595 // the actual IEEE respresentation.  We compensate for that here.
1596
1597 APInt
1598 APFloat::convertF80LongDoubleAPFloatToAPInt() const {
1599   assert(semantics == (const llvm::fltSemantics* const)&x87DoubleExtended);
1600   assert (partCount()==1);
1601
1602   uint64_t myexponent, mysignificand;
1603
1604   if (category==fcNormal) {
1605     myexponent = exponent+16383; //bias
1606     mysignificand = *significandParts();
1607     if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
1608       myexponent = 0;   // denormal
1609   } else if (category==fcZero) {
1610     myexponent = 0;
1611     mysignificand = 0;
1612   } else if (category==fcInfinity) {
1613     myexponent = 0x7fff;
1614     mysignificand = 0x8000000000000000ULL;
1615   } else if (category==fcNaN) {
1616     myexponent = 0x7fff;
1617     mysignificand = *significandParts();
1618   } else
1619     assert(0);
1620
1621   uint64_t words[2];
1622   words[0] =  (((uint64_t)sign & 1) << 63) |
1623               ((myexponent & 0x7fff) <<  48) |
1624               ((mysignificand >>16) & 0xffffffffffffLL);
1625   words[1] = mysignificand & 0xffff;
1626   APInt api(80, 2, words);
1627   return api;
1628 }
1629
1630 APInt
1631 APFloat::convertDoubleAPFloatToAPInt() const {
1632   assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
1633   assert (partCount()==1);
1634
1635   uint64_t myexponent, mysignificand;
1636
1637   if (category==fcNormal) {
1638     myexponent = exponent+1023; //bias
1639     mysignificand = *significandParts();
1640     if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
1641       myexponent = 0;   // denormal
1642   } else if (category==fcZero) {
1643     myexponent = 0;
1644     mysignificand = 0;
1645   } else if (category==fcInfinity) {
1646     myexponent = 0x7ff;
1647     mysignificand = 0;
1648   } else if (category==fcNaN) {
1649     myexponent = 0x7ff;
1650     mysignificand = *significandParts();
1651   } else
1652     assert(0);
1653
1654   APInt api(64, (((((uint64_t)sign & 1) << 63) |
1655                  ((myexponent & 0x7ff) <<  52) |
1656                  (mysignificand & 0xfffffffffffffLL))));
1657   return api;
1658 }
1659
1660 APInt
1661 APFloat::convertFloatAPFloatToAPInt() const {
1662   assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
1663   assert (partCount()==1);
1664
1665   uint32_t myexponent, mysignificand;
1666
1667   if (category==fcNormal) {
1668     myexponent = exponent+127; //bias
1669     mysignificand = *significandParts();
1670     if (myexponent == 1 && !(mysignificand & 0x400000))
1671       myexponent = 0;   // denormal
1672   } else if (category==fcZero) {
1673     myexponent = 0;
1674     mysignificand = 0;
1675   } else if (category==fcInfinity) {
1676     myexponent = 0xff;
1677     mysignificand = 0;
1678   } else if (category==fcNaN) {
1679     myexponent = 0xff;
1680     mysignificand = *significandParts();
1681   } else
1682     assert(0);
1683
1684   APInt api(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
1685                  (mysignificand & 0x7fffff)));
1686   return api;
1687 }
1688
1689 APInt
1690 APFloat::convertToAPInt() const {
1691   if (semantics == (const llvm::fltSemantics* const)&IEEEsingle)
1692     return convertFloatAPFloatToAPInt();
1693   else if (semantics == (const llvm::fltSemantics* const)&IEEEdouble)
1694     return convertDoubleAPFloatToAPInt();
1695   else if (semantics == (const llvm::fltSemantics* const)&x87DoubleExtended)
1696     return convertF80LongDoubleAPFloatToAPInt();
1697   else
1698     assert(0);
1699 }
1700
1701 float
1702 APFloat::convertToFloat() const {
1703   assert(semantics == (const llvm::fltSemantics* const)&IEEEsingle);
1704   APInt api = convertToAPInt();
1705   return api.bitsToFloat();
1706 }
1707
1708 double
1709 APFloat::convertToDouble() const {
1710   assert(semantics == (const llvm::fltSemantics* const)&IEEEdouble);
1711   APInt api = convertToAPInt();
1712   return api.bitsToDouble();
1713 }
1714
1715 /// Integer bit is explicit in this format.  Current Intel book does not
1716 /// define meaning of:
1717 ///  exponent = all 1's, integer bit not set.
1718 ///  exponent = 0, integer bit set. (formerly "psuedodenormals")
1719 ///  exponent!=0 nor all 1's, integer bit not set. (formerly "unnormals")
1720 void
1721 APFloat::initFromF80LongDoubleAPInt(const APInt &api) {
1722   assert(api.getBitWidth()==80);
1723   uint64_t i1 = api.getRawData()[0];
1724   uint64_t i2 = api.getRawData()[1];
1725   uint64_t myexponent = (i1 >> 48) & 0x7fff;
1726   uint64_t mysignificand = ((i1 << 16) &  0xffffffffffff0000ULL) |
1727                           (i2 & 0xffff);
1728
1729   initialize(&APFloat::x87DoubleExtended);
1730   assert(partCount()==1);
1731
1732   sign = i1>>63;
1733   if (myexponent==0 && mysignificand==0) {
1734     // exponent, significand meaningless
1735     category = fcZero;
1736   } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
1737     // exponent, significand meaningless
1738     category = fcInfinity;
1739   } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
1740     // exponent meaningless
1741     category = fcNaN;
1742     *significandParts() = mysignificand;
1743   } else {
1744     category = fcNormal;
1745     exponent = myexponent - 16383;
1746     *significandParts() = mysignificand;
1747     if (myexponent==0)          // denormal
1748       exponent = -16382;
1749  }
1750 }
1751
1752 void
1753 APFloat::initFromDoubleAPInt(const APInt &api) {
1754   assert(api.getBitWidth()==64);
1755   uint64_t i = *api.getRawData();
1756   uint64_t myexponent = (i >> 52) & 0x7ff;
1757   uint64_t mysignificand = i & 0xfffffffffffffLL;
1758
1759   initialize(&APFloat::IEEEdouble);
1760   assert(partCount()==1);
1761
1762   sign = i>>63;
1763   if (myexponent==0 && mysignificand==0) {
1764     // exponent, significand meaningless
1765     category = fcZero;
1766   } else if (myexponent==0x7ff && mysignificand==0) {
1767     // exponent, significand meaningless
1768     category = fcInfinity;
1769   } else if (myexponent==0x7ff && mysignificand!=0) {
1770     // exponent meaningless
1771     category = fcNaN;
1772     *significandParts() = mysignificand;
1773   } else {
1774     category = fcNormal;
1775     exponent = myexponent - 1023;
1776     *significandParts() = mysignificand;
1777     if (myexponent==0)          // denormal
1778       exponent = -1022;
1779     else
1780       *significandParts() |= 0x10000000000000LL;  // integer bit
1781  }
1782 }
1783
1784 void
1785 APFloat::initFromFloatAPInt(const APInt & api) {
1786   assert(api.getBitWidth()==32);
1787   uint32_t i = (uint32_t)*api.getRawData();
1788   uint32_t myexponent = (i >> 23) & 0xff;
1789   uint32_t mysignificand = i & 0x7fffff;
1790
1791   initialize(&APFloat::IEEEsingle);
1792   assert(partCount()==1);
1793
1794   sign = i >> 31;
1795   if (myexponent==0 && mysignificand==0) {
1796     // exponent, significand meaningless
1797     category = fcZero;
1798   } else if (myexponent==0xff && mysignificand==0) {
1799     // exponent, significand meaningless
1800     category = fcInfinity;
1801   } else if (myexponent==0xff && (mysignificand & 0x400000)) {
1802     // sign, exponent, significand meaningless
1803     category = fcNaN;
1804     *significandParts() = mysignificand;
1805   } else {
1806     category = fcNormal;
1807     exponent = myexponent - 127;  //bias
1808     *significandParts() = mysignificand;
1809     if (myexponent==0)    // denormal
1810       exponent = -126;
1811     else
1812       *significandParts() |= 0x800000; // integer bit
1813   }
1814 }
1815
1816 /// Treat api as containing the bits of a floating point number.  Currently
1817 /// we infer the floating point type from the size of the APInt.  FIXME: This
1818 /// breaks when we get to PPC128 and IEEE128 (but both cannot exist in the
1819 /// same compile...)
1820 void
1821 APFloat::initFromAPInt(const APInt& api) {
1822   if (api.getBitWidth() == 32)
1823     return initFromFloatAPInt(api);
1824   else if (api.getBitWidth()==64)
1825     return initFromDoubleAPInt(api);
1826   else if (api.getBitWidth()==80)
1827     return initFromF80LongDoubleAPInt(api);
1828   else
1829     assert(0);
1830 }
1831
1832 APFloat::APFloat(const APInt& api) {
1833   initFromAPInt(api);
1834 }
1835
1836 APFloat::APFloat(float f) {
1837   APInt api = APInt(32, 0);
1838   initFromAPInt(api.floatToBits(f));
1839 }
1840
1841 APFloat::APFloat(double d) {
1842   APInt api = APInt(64, 0);
1843   initFromAPInt(api.doubleToBits(d));
1844 }
1845