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