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