APFloat.cpp revision 193323
1175061Sobrien//===-- APFloat.cpp - Implement APFloat class -----------------------------===// 21590Srgrimes// 31590Srgrimes// The LLVM Compiler Infrastructure 41590Srgrimes// 51590Srgrimes// This file is distributed under the University of Illinois Open Source 61590Srgrimes// License. See LICENSE.TXT for details. 71590Srgrimes// 81590Srgrimes//===----------------------------------------------------------------------===// 91590Srgrimes// 101590Srgrimes// This file implements a class to represent arbitrary precision floating 111590Srgrimes// point values and provide a variety of arithmetic operations on them. 121590Srgrimes// 131590Srgrimes//===----------------------------------------------------------------------===// 141590Srgrimes 151590Srgrimes#include "llvm/ADT/APFloat.h" 161590Srgrimes#include "llvm/ADT/FoldingSet.h" 171590Srgrimes#include "llvm/Support/MathExtras.h" 181590Srgrimes#include <cstring> 191590Srgrimes 201590Srgrimesusing namespace llvm; 211590Srgrimes 221590Srgrimes#define convolve(lhs, rhs) ((lhs) * 4 + (rhs)) 231590Srgrimes 241590Srgrimes/* Assumed in hexadecimal significand parsing, and conversion to 251590Srgrimes hexadecimal strings. */ 261590Srgrimes#define COMPILE_TIME_ASSERT(cond) extern int CTAssert[(cond) ? 1 : -1] 271590SrgrimesCOMPILE_TIME_ASSERT(integerPartWidth % 4 == 0); 281590Srgrimes 291590Srgrimesnamespace llvm { 30132671Scharnier 311590Srgrimes /* Represents floating point arithmetic semantics. */ 3213430Speter struct fltSemantics { 33132671Scharnier /* The largest E such that 2^E is representable; this matches the 345811Swollman definition of IEEE 754. */ 351590Srgrimes exponent_t maxExponent; 36132671Scharnier 37132671Scharnier /* The smallest E such that 2^E is a normalized number; this 38132671Scharnier matches the definition of IEEE 754. */ 391590Srgrimes exponent_t minExponent; 401590Srgrimes 411590Srgrimes /* Number of bits in the significand. This includes the integer 42171465Sjhb bit. */ 4320287Swollman unsigned int precision; 441590Srgrimes 4593957Sru /* True if arithmetic is supported. */ 461590Srgrimes unsigned int arithmeticOK; 4721263Swollman }; 481590Srgrimes 491590Srgrimes const fltSemantics APFloat::IEEEsingle = { 127, -126, 24, true }; 50120716Ssam const fltSemantics APFloat::IEEEdouble = { 1023, -1022, 53, true }; 511590Srgrimes const fltSemantics APFloat::IEEEquad = { 16383, -16382, 113, true }; 5216178Sjulian const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64, true }; 531590Srgrimes const fltSemantics APFloat::Bogus = { 0, 0, 0, true }; 5411819Sjulian 5516178Sjulian // The PowerPC format consists of two doubles. It does not map cleanly 5652419Sjulian // onto the usual format above. For now only storage of constants of 5711819Sjulian // this type is supported, no arithmetic. 581590Srgrimes const fltSemantics APFloat::PPCDoubleDouble = { 1023, -1022, 106, false }; 591590Srgrimes 6054263Sshin /* A tight upper bound on number of parts required to hold the value 6174262Sbrian pow(5, power) is 621590Srgrimes 63160787Syar power * 815 / (351 * integerPartWidth) + 1 641590Srgrimes 651590Srgrimes However, whilst the result may require only this many parts, 661590Srgrimes because we are multiplying two values to get it, the 67160130Soleg multiplication may require an extra part with the excess part 68200462Sdelphij being zero (consider the trivial case of 1 * 1, tcFullMultiply 695811Swollman requires two parts to hold the single-part result). So we add an 701590Srgrimes extra one to guarantee enough space whilst multiplying. */ 711590Srgrimes const unsigned int maxExponent = 16383; 72175061Sobrien const unsigned int maxPrecision = 113; 731590Srgrimes const unsigned int maxPowerOfFiveExponent = maxExponent + maxPrecision - 1; 741590Srgrimes const unsigned int maxPowerOfFiveParts = 2 + ((maxPowerOfFiveExponent * 815) 751590Srgrimes / (351 * integerPartWidth)); 761590Srgrimes} 771590Srgrimes 785103Swollman/* A bunch of private, handy routines. */ 791590Srgrimes 801590Srgrimesstatic inline unsigned int 811590SrgrimespartCountForBits(unsigned int bits) 821590Srgrimes{ 831590Srgrimes return ((bits) + integerPartWidth - 1) / integerPartWidth; 841590Srgrimes} 851590Srgrimes 861590Srgrimes/* Returns 0U-9U. Return values >= 10U are not digits. */ 871590Srgrimesstatic inline unsigned int 881590SrgrimesdecDigitValue(unsigned int c) 891590Srgrimes{ 901590Srgrimes return c - '0'; 911590Srgrimes} 925103Swollman 935103Swollmanstatic unsigned int 9416392SpeterhexDigitValue(unsigned int c) 9516392Speter{ 96186119Sqingli unsigned int r; 97186119Sqingli 98186119Sqingli r = c - '0'; 99186119Sqingli if(r <= 9) 100186119Sqingli return r; 101186119Sqingli 102186119Sqingli r = c - 'A'; 103186119Sqingli if(r <= 5) 104186119Sqingli return r + 10; 105102975Sdwmalone 1061590Srgrimes r = c - 'a'; 1071590Srgrimes if(r <= 5) 10835383Sphk return r + 10; 10938428Sjb 1101590Srgrimes return -1U; 1111590Srgrimes} 11235383Sphk 1131590Srgrimesstatic inline void 11435383SphkassertArithmeticOK(const llvm::fltSemantics &semantics) { 11535383Sphk assert(semantics.arithmeticOK 1161590Srgrimes && "Compile-time arithmetic does not support these semantics"); 1171590Srgrimes} 1181590Srgrimes 1191590Srgrimes/* Return the value of a decimal exponent of the form 120193232Sbz [+-]ddddddd. 1211590Srgrimes 1221590Srgrimes If the exponent overflows, returns a large exponent with the 1231590Srgrimes appropriate sign. */ 124160130Solegstatic int 125160130SolegreadExponent(const char *p) 126175061Sobrien{ 127175061Sobrien bool isNegative; 128175061Sobrien unsigned int absExponent; 129175061Sobrien const unsigned int overlargeExponent = 24000; /* FIXME. */ 130175061Sobrien 131175061Sobrien isNegative = (*p == '-'); 132175061Sobrien if (*p == '-' || *p == '+') 133175061Sobrien p++; 134175061Sobrien 135175061Sobrien absExponent = decDigitValue(*p++); 136182602Sobrien assert (absExponent < 10U); 137175061Sobrien 13897878Skbyanc for (;;) { 139175061Sobrien unsigned int value; 140176194Smarius 1411590Srgrimes value = decDigitValue(*p); 1421590Srgrimes if (value >= 10U) 1431590Srgrimes break; 1441590Srgrimes 1451590Srgrimes p++; 146253275Shrs value += absExponent * 10; 1471590Srgrimes if (absExponent >= overlargeExponent) { 148193232Sbz absExponent = overlargeExponent; 149182602Sobrien break; 150253275Shrs } 1511590Srgrimes absExponent = value; 152178912Sdelphij } 153253275Shrs 154253275Shrs if (isNegative) 155178887Sjulian return -(int) absExponent; 156178912Sdelphij else 157178887Sjulian return (int) absExponent; 158253275Shrs} 159253275Shrs 160193232Sbz/* This is ugly and needs cleaning up, but I don't immediately see 161193232Sbz how whilst remaining safe. */ 162178887Sjulianstatic int 163178887SjuliantotalExponent(const char *p, int exponentAdjustment) 164160130Soleg{ 165160130Soleg int unsignedExponent; 166160130Soleg bool negative, overflow; 167160130Soleg int exponent; 168160130Soleg 169160130Soleg /* Move past the exponent letter and sign to the digits. */ 170160130Soleg p++; 171160130Soleg negative = *p == '-'; 172253275Shrs if(*p == '-' || *p == '+') 173253275Shrs p++; 174253275Shrs 175253275Shrs unsignedExponent = 0; 1761590Srgrimes overflow = false; 1771590Srgrimes for(;;) { 1781590Srgrimes unsigned int value; 1791590Srgrimes 1801590Srgrimes value = decDigitValue(*p); 1811590Srgrimes if(value >= 10U) 1821590Srgrimes break; 1831590Srgrimes 1841590Srgrimes p++; 185193232Sbz unsignedExponent = unsignedExponent * 10 + value; 186193232Sbz if(unsignedExponent > 65535) 187176289Sjhb overflow = true; 188231852Sbz } 189178887Sjulian 190231852Sbz if(exponentAdjustment > 65535 || exponentAdjustment < -65536) 191231852Sbz overflow = true; 192231852Sbz 193231852Sbz if(!overflow) { 194231852Sbz exponent = unsignedExponent; 195231852Sbz if(negative) 196231852Sbz exponent = -exponent; 197178887Sjulian exponent += exponentAdjustment; 198231852Sbz if(exponent > 65535 || exponent < -65536) 199193232Sbz overflow = true; 200193232Sbz } 201231852Sbz 202193232Sbz if(overflow) 203193232Sbz exponent = negative ? -65536: 65535; 2041590Srgrimes 205193232Sbz return exponent; 206193232Sbz} 207193232Sbz 208176289Sjhbstatic const char * 209176289SjhbskipLeadingZeroesAndAnyDot(const char *p, const char **dot) 210231852Sbz{ 2111590Srgrimes *dot = 0; 2121590Srgrimes while(*p == '0') 2131590Srgrimes p++; 2141590Srgrimes 215231852Sbz if(*p == '.') { 216231852Sbz *dot = p++; 217231852Sbz while(*p == '0') 2181590Srgrimes p++; 219231852Sbz } 2201590Srgrimes 2211590Srgrimes return p; 2221590Srgrimes} 2231590Srgrimes 2241590Srgrimes/* Given a normal decimal floating point number of the form 2251590Srgrimes 2261590Srgrimes dddd.dddd[eE][+-]ddd 2271590Srgrimes 2281590Srgrimes where the decimal point and exponent are optional, fill out the 2291590Srgrimes structure D. Exponent is appropriate if the significand is 230102975Sdwmalone treated as an integer, and normalizedExponent if the significand 2311590Srgrimes is taken to have the decimal point after a single leading 232102975Sdwmalone non-zero digit. 2331590Srgrimes 234102975Sdwmalone If the value is zero, V->firstSigDigit points to a non-digit, and 2351590Srgrimes the return exponent is zero. 2361590Srgrimes*/ 2371590Srgrimesstruct decimalInfo { 23854263Sshin const char *firstSigDigit; 23954263Sshin const char *lastSigDigit; 24054263Sshin int exponent; 24154263Sshin int normalizedExponent; 24254263Sshin}; 24311819Sjulian 24411819Sjulianstatic void 24511819SjulianinterpretDecimal(const char *p, decimalInfo *D) 2461590Srgrimes{ 2471590Srgrimes const char *dot; 2481590Srgrimes 24916178Sjulian p = skipLeadingZeroesAndAnyDot (p, &dot); 25017254Sjulian 25116178Sjulian D->firstSigDigit = p; 2521590Srgrimes D->exponent = 0; 2531590Srgrimes D->normalizedExponent = 0; 2541590Srgrimes 25552419Sjulian for (;;) { 25652419Sjulian if (*p == '.') { 25752419Sjulian assert(dot == 0); 2581590Srgrimes dot = p++; 2591590Srgrimes } 2601590Srgrimes if (decDigitValue(*p) >= 10U) 2611590Srgrimes break; 2621590Srgrimes p++; 2631590Srgrimes } 2641590Srgrimes 265102975Sdwmalone /* If number is all zerooes accept any exponent. */ 2661590Srgrimes if (p != D->firstSigDigit) { 2671590Srgrimes if (*p == 'e' || *p == 'E') 2681590Srgrimes D->exponent = readExponent(p + 1); 26962584Sitojun 27097878Skbyanc /* Implied decimal point? */ 27197878Skbyanc if (!dot) 272123030Sbms dot = p; 27362584Sitojun 27497878Skbyanc /* Drop insignificant trailing zeroes. */ 27597878Skbyanc do 27697878Skbyanc do 27797878Skbyanc p--; 278123030Sbms while (*p == '0'); 27954263Sshin while (*p == '.'); 2801590Srgrimes 28197878Skbyanc /* Adjust the exponents for any decimal point. */ 28297878Skbyanc D->exponent += static_cast<exponent_t>((dot - p) - (dot > p)); 28397878Skbyanc D->normalizedExponent = (D->exponent + 28497878Skbyanc static_cast<exponent_t>((p - D->firstSigDigit) 28597878Skbyanc - (dot > D->firstSigDigit && dot < p))); 28697878Skbyanc } 28797878Skbyanc 28897878Skbyanc D->lastSigDigit = p; 28997878Skbyanc} 29097878Skbyanc 291175207Ssam/* Return the trailing fraction of a hexadecimal number. 29297878Skbyanc DIGITVALUE is the first hex digit of the fraction, P points to 29397878Skbyanc the next digit. */ 29497878Skbyancstatic lostFraction 29597878SkbyanctrailingHexadecimalFraction(const char *p, unsigned int digitValue) 29697878Skbyanc{ 29797878Skbyanc unsigned int hexDigit; 29897878Skbyanc 29997878Skbyanc /* If the first trailing digit isn't 0 or 8 we can work out the 30097878Skbyanc fraction immediately. */ 30197878Skbyanc if(digitValue > 8) 30297878Skbyanc return lfMoreThanHalf; 30397878Skbyanc else if(digitValue < 8 && digitValue > 0) 30497878Skbyanc return lfLessThanHalf; 30597878Skbyanc 30697878Skbyanc /* Otherwise we need to find the first non-zero digit. */ 30797878Skbyanc while(*p == '0') 30897878Skbyanc p++; 30997878Skbyanc 310176289Sjhb hexDigit = hexDigitValue(*p); 311176289Sjhb 312176289Sjhb /* If we ran off the end it is exactly zero or one-half, otherwise 313176289Sjhb a little more. */ 31497878Skbyanc if(hexDigit == -1U) 31597878Skbyanc return digitValue == 0 ? lfExactlyZero: lfExactlyHalf; 316176289Sjhb else 317176289Sjhb return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf; 31897878Skbyanc} 31997878Skbyanc 32097878Skbyanc/* Return the fraction lost were a bignum truncated losing the least 32197878Skbyanc significant BITS bits. */ 32297878Skbyancstatic lostFraction 32397878SkbyanclostFractionThroughTruncation(const integerPart *parts, 32497878Skbyanc unsigned int partCount, 32597878Skbyanc unsigned int bits) 32697878Skbyanc{ 32797878Skbyanc unsigned int lsb; 32897878Skbyanc 32997878Skbyanc lsb = APInt::tcLSB(parts, partCount); 33097878Skbyanc 33197878Skbyanc /* Note this is guaranteed true if bits == 0, or LSB == -1U. */ 33297878Skbyanc if(bits <= lsb) 33397878Skbyanc return lfExactlyZero; 33497878Skbyanc if(bits == lsb + 1) 33597878Skbyanc return lfExactlyHalf; 33697878Skbyanc if(bits <= partCount * integerPartWidth 33797878Skbyanc && APInt::tcExtractBit(parts, bits - 1)) 33897878Skbyanc return lfMoreThanHalf; 33997878Skbyanc 34097878Skbyanc return lfLessThanHalf; 34197878Skbyanc} 34297878Skbyanc 34397878Skbyanc/* Shift DST right BITS bits noting lost fraction. */ 34497878Skbyancstatic lostFraction 34597878SkbyancshiftRight(integerPart *dst, unsigned int parts, unsigned int bits) 34697878Skbyanc{ 34797878Skbyanc lostFraction lost_fraction; 34897878Skbyanc 34997878Skbyanc lost_fraction = lostFractionThroughTruncation(dst, parts, bits); 35097878Skbyanc 35197878Skbyanc APInt::tcShiftRight(dst, parts, bits); 35297878Skbyanc 35397878Skbyanc return lost_fraction; 35497878Skbyanc} 35597878Skbyanc 35697878Skbyanc/* Combine the effect of two lost fractions. */ 35797878Skbyancstatic lostFraction 358186119SqinglicombineLostFractions(lostFraction moreSignificant, 35997878Skbyanc lostFraction lessSignificant) 36097878Skbyanc{ 36197878Skbyanc if(lessSignificant != lfExactlyZero) { 36297878Skbyanc if(moreSignificant == lfExactlyZero) 36397878Skbyanc moreSignificant = lfLessThanHalf; 36497878Skbyanc else if(moreSignificant == lfExactlyHalf) 36597878Skbyanc moreSignificant = lfMoreThanHalf; 36697878Skbyanc } 36797878Skbyanc 36897878Skbyanc return moreSignificant; 36997878Skbyanc} 370176289Sjhb 371176289Sjhb/* The error from the true value, in half-ulps, on multiplying two 372176289Sjhb floating point numbers, which differ from the value they 373176289Sjhb approximate by at most HUE1 and HUE2 half-ulps, is strictly less 37497878Skbyanc than the returned value. 37597878Skbyanc 37697878Skbyanc See "How to Read Floating Point Numbers Accurately" by William D 37797878Skbyanc Clinger. */ 37897878Skbyancstatic unsigned int 37997878SkbyancHUerrBound(bool inexactMultiply, unsigned int HUerr1, unsigned int HUerr2) 38097878Skbyanc{ 381160130Soleg assert(HUerr1 < 2 || HUerr2 < 2 || (HUerr1 + HUerr2 < 8)); 38297878Skbyanc 38397878Skbyanc if (HUerr1 + HUerr2 == 0) 38497878Skbyanc return inexactMultiply * 2; /* <= inexactMultiply half-ulps. */ 38597878Skbyanc else 38697878Skbyanc return inexactMultiply + 2 * (HUerr1 + HUerr2); 38797878Skbyanc} 38897878Skbyanc 38997878Skbyanc/* The number of ulps from the boundary (zero, or half if ISNEAREST) 39097878Skbyanc when the least significant BITS are truncated. BITS cannot be 3911590Srgrimes zero. */ 3921590Srgrimesstatic integerPart 3931590SrgrimesulpsFromBoundary(const integerPart *parts, unsigned int bits, bool isNearest) 3941590Srgrimes{ 395102975Sdwmalone unsigned int count, partBits; 3961590Srgrimes integerPart part, boundary; 39754946Sshin 3981590Srgrimes assert (bits != 0); 3991590Srgrimes 400102975Sdwmalone bits--; 40197878Skbyanc count = bits / integerPartWidth; 40297878Skbyanc partBits = bits % integerPartWidth + 1; 40397878Skbyanc 40497878Skbyanc part = parts[count] & (~(integerPart) 0 >> (integerPartWidth - partBits)); 40597878Skbyanc 40697878Skbyanc if (isNearest) 40797878Skbyanc boundary = (integerPart) 1 << (partBits - 1); 40897878Skbyanc else 40997878Skbyanc boundary = 0; 41097878Skbyanc 41197878Skbyanc if (count == 0) { 41297878Skbyanc if (part - boundary <= boundary - part) 41397878Skbyanc return part - boundary; 41497878Skbyanc else 41597878Skbyanc return boundary - part; 41697878Skbyanc } 41797878Skbyanc 41897878Skbyanc if (part == boundary) { 41997878Skbyanc while (--count) 42097878Skbyanc if (parts[count]) 42197878Skbyanc return ~(integerPart) 0; /* A lot. */ 42297878Skbyanc 42397878Skbyanc return parts[0]; 42497878Skbyanc } else if (part == boundary - 1) { 42597878Skbyanc while (--count) 42697878Skbyanc if (~parts[count]) 42797878Skbyanc return ~(integerPart) 0; /* A lot. */ 42897878Skbyanc 4291590Srgrimes return -parts[0]; 4301590Srgrimes } 4311590Srgrimes 43278314Sassar return ~(integerPart) 0; /* A lot. */ 4331590Srgrimes} 4341590Srgrimes 435176289Sjhb/* Place pow(5, power) in DST, and return the number of parts used. 436176289Sjhb DST must be at least one part larger than size of the answer. */ 4371590Srgrimesstatic unsigned int 4381590SrgrimespowerOf5(integerPart *dst, unsigned int power) 4391590Srgrimes{ 4401590Srgrimes static const integerPart firstEightPowers[] = { 1, 5, 25, 125, 625, 3125, 4411590Srgrimes 15625, 78125 }; 4421590Srgrimes integerPart pow5s[maxPowerOfFiveParts * 2 + 5]; 44378314Sassar pow5s[0] = 78125 * 5; 4441590Srgrimes 4451590Srgrimes unsigned int partsCount[16] = { 1 }; 4461590Srgrimes integerPart scratch[maxPowerOfFiveParts], *p1, *p2, *pow5; 447176289Sjhb unsigned int result; 448176289Sjhb assert(power <= maxExponent); 449176289Sjhb 450176289Sjhb p1 = dst; 45159540Smarkm p2 = scratch; 4521590Srgrimes 45338428Sjb *p1 = firstEightPowers[power & 7]; 4541590Srgrimes power >>= 3; 4551590Srgrimes 4561590Srgrimes result = 1; 4571590Srgrimes pow5 = pow5s; 4581590Srgrimes 459176289Sjhb for (unsigned int n = 0; power; power >>= 1, n++) { 460176289Sjhb unsigned int pc; 461176289Sjhb 462176289Sjhb pc = partsCount[n]; 463176289Sjhb 4641590Srgrimes /* Calculate pow(5,pow(2,n+3)) if we haven't yet. */ 4651590Srgrimes if (pc == 0) { 46613430Speter pc = partsCount[n - 1]; 4671590Srgrimes APInt::tcFullMultiply(pow5, pow5 - pc, pow5 - pc, pc, pc); 4681590Srgrimes pc *= 2; 46916080Salex if (pow5[pc - 1] == 0) 4701590Srgrimes pc--; 4711590Srgrimes partsCount[n] = pc; 4721590Srgrimes } 47338428Sjb 4741590Srgrimes if (power & 1) { 4751590Srgrimes integerPart *tmp; 47659540Smarkm 47759540Smarkm APInt::tcFullMultiply(p2, p1, pow5, result, pc); 4781590Srgrimes result += pc; 4791590Srgrimes if (p2[result - 1] == 0) 4801590Srgrimes result--; 4811590Srgrimes 4821590Srgrimes /* Now result is in p1 with partsCount parts and p2 is scratch 4831590Srgrimes space. */ 4841590Srgrimes tmp = p1, p1 = p2, p2 = tmp; 48578314Sassar } 4861590Srgrimes 4871590Srgrimes pow5 += pc; 4881590Srgrimes } 48959540Smarkm 4901590Srgrimes if (p1 != dst) 4911590Srgrimes APInt::tcAssign(dst, p1, result); 4921590Srgrimes 49313430Speter return result; 4941590Srgrimes} 4951590Srgrimes 4961590Srgrimes/* Zero at the end to avoid modular arithmetic when adding one; used 49759540Smarkm when rounding up during hexadecimal output. */ 49859540Smarkmstatic const char hexDigitsLower[] = "0123456789abcdef0"; 4991590Srgrimesstatic const char hexDigitsUpper[] = "0123456789ABCDEF0"; 5001590Srgrimesstatic const char infinityL[] = "infinity"; 501176289Sjhbstatic const char infinityU[] = "INFINITY"; 502176289Sjhbstatic const char NaNL[] = "nan"; 5031590Srgrimesstatic const char NaNU[] = "NAN"; 50438428Sjb 50559540Smarkm/* Write out an integerPart in hexadecimal, starting with the most 50613430Speter significant nibble. Write out exactly COUNT hexdigits, return 50713430Speter COUNT. */ 50813430Speterstatic unsigned int 509176289SjhbpartAsHex (char *dst, integerPart part, unsigned int count, 510176289Sjhb const char *hexDigitChars) 51113430Speter{ 512176289Sjhb unsigned int result = count; 513176289Sjhb 51413430Speter assert (count != 0 && count <= integerPartWidth / 4); 51513430Speter 51613430Speter part >>= (integerPartWidth - 4 * count); 5171590Srgrimes while (count--) { 51816080Salex dst[count] = hexDigitChars[part & 0xf]; 5191590Srgrimes part >>= 4; 5201590Srgrimes } 5211590Srgrimes 5221590Srgrimes return result; 5231590Srgrimes} 5241590Srgrimes 52578314Sassar/* Write out an unsigned decimal integer. */ 5261590Srgrimesstatic char * 5271590SrgrimeswriteUnsignedDecimal (char *dst, unsigned int n) 5281590Srgrimes{ 5291590Srgrimes char buff[40], *p; 530102975Sdwmalone 5311590Srgrimes p = buff; 53235308Sphk do 53335308Sphk *p++ = '0' + n % 10; 53435308Sphk while (n /= 10); 53535308Sphk 53635308Sphk do 53735308Sphk *dst++ = *--p; 53835308Sphk while (p != buff); 5395811Swollman 5405811Swollman return dst; 5415811Swollman} 5425811Swollman 543132671Scharnier/* Write out a signed decimal integer. */ 5445811Swollmanstatic char * 54535308SphkwriteSignedDecimal (char *dst, int value) 5465811Swollman{ 5475811Swollman if (value < 0) { 5481590Srgrimes *dst++ = '-'; 5491590Srgrimes dst = writeUnsignedDecimal(dst, -(unsigned) value); 5501590Srgrimes } else 5511590Srgrimes dst = writeUnsignedDecimal(dst, value); 5521590Srgrimes 5531590Srgrimes return dst; 5541590Srgrimes} 5551590Srgrimes 55678314Sassar/* Constructors. */ 5571590Srgrimesvoid 558102975SdwmaloneAPFloat::initialize(const fltSemantics *ourSemantics) 5591590Srgrimes{ 5601590Srgrimes unsigned int count; 5611590Srgrimes 5621590Srgrimes semantics = ourSemantics; 563102975Sdwmalone count = partCount(); 5641590Srgrimes if(count > 1) 5651590Srgrimes significand.parts = new integerPart[count]; 5661590Srgrimes} 5671590Srgrimes 5681590Srgrimesvoid 5691590SrgrimesAPFloat::freeSignificand() 5701590Srgrimes{ 5711590Srgrimes if(partCount() > 1) 5721590Srgrimes delete [] significand.parts; 5731590Srgrimes} 574102975Sdwmalone 5751590Srgrimesvoid 5761590SrgrimesAPFloat::assign(const APFloat &rhs) 5771590Srgrimes{ 578102975Sdwmalone assert(semantics == rhs.semantics); 579102975Sdwmalone 580102975Sdwmalone sign = rhs.sign; 581102975Sdwmalone category = rhs.category; 5821590Srgrimes exponent = rhs.exponent; 5831590Srgrimes sign2 = rhs.sign2; 58413430Speter exponent2 = rhs.exponent2; 5851590Srgrimes if(category == fcNormal || category == fcNaN) 58613430Speter copySignificand(rhs); 587128186Sluigi} 58813430Speter 5891590Srgrimesvoid 5901590SrgrimesAPFloat::copySignificand(const APFloat &rhs) 5911590Srgrimes{ 5921590Srgrimes assert(category == fcNormal || category == fcNaN); 5931590Srgrimes assert(rhs.partCount() >= partCount()); 5941590Srgrimes 59578314Sassar APInt::tcAssign(significandParts(), rhs.significandParts(), 5961590Srgrimes partCount()); 59797878Skbyanc} 5981590Srgrimes 59997878Skbyanc/* Make this number a NaN, with an arbitrary but deterministic value 60097878Skbyanc for the significand. If double or longer, this is a signalling NaN, 60197878Skbyanc which may not be ideal. If float, this is QNaN(0). */ 60297878Skbyancvoid 60397878SkbyancAPFloat::makeNaN(unsigned type) 60497878Skbyanc{ 60597878Skbyanc category = fcNaN; 60697878Skbyanc // FIXME: Add double and long double support for QNaN(0). 60797878Skbyanc if (semantics->precision == 24 && semantics->maxExponent == 127) { 60897878Skbyanc type |= 0x7fc00000U; 60997878Skbyanc type &= ~0x80000000U; 61097878Skbyanc } else 61197878Skbyanc type = ~0U; 61297878Skbyanc APInt::tcSet(significandParts(), type, partCount()); 61397878Skbyanc} 61497878Skbyanc 615102975SdwmaloneAPFloat & 61697878SkbyancAPFloat::operator=(const APFloat &rhs) 617176289Sjhb{ 618176289Sjhb if(this != &rhs) { 619176289Sjhb if(semantics != rhs.semantics) { 6201590Srgrimes freeSignificand(); 6211590Srgrimes initialize(rhs.semantics); 6221590Srgrimes } 623132671Scharnier assign(rhs); 6241590Srgrimes } 625132671Scharnier 62655575Srgrimes return *this; 62755575Srgrimes} 62855575Srgrimes 62955575Srgrimesbool 63013430SpeterAPFloat::bitwiseIsEqual(const APFloat &rhs) const { 631132671Scharnier if (this == &rhs) 63213430Speter return true; 633132671Scharnier if (semantics != rhs.semantics || 63413430Speter category != rhs.category || 63513430Speter sign != rhs.sign) 63613430Speter return false; 637132671Scharnier if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble && 6381590Srgrimes sign2 != rhs.sign2) 6391590Srgrimes return false; 6401590Srgrimes if (category==fcZero || category==fcInfinity) 64154263Sshin return true; 64254263Sshin else if (category==fcNormal && exponent!=rhs.exponent) 64354263Sshin return false; 64454263Sshin else if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble && 64554263Sshin exponent2!=rhs.exponent2) 646243187Shrs return false; 647243187Shrs else { 648243187Shrs int i= partCount(); 649243187Shrs const integerPart* p=significandParts(); 650243187Shrs const integerPart* q=rhs.significandParts(); 651243187Shrs for (; i>0; i--, p++, q++) { 652217642Sume if (*p != *q) 65354263Sshin return false; 65454263Sshin } 65554263Sshin return true; 65654263Sshin } 65754263Sshin} 65854263Sshin 65954263SshinAPFloat::APFloat(const fltSemantics &ourSemantics, integerPart value) 66054263Sshin{ 66154263Sshin assertArithmeticOK(ourSemantics); 66254263Sshin initialize(&ourSemantics); 66354263Sshin sign = 0; 66454263Sshin zeroSignificand(); 66554263Sshin exponent = ourSemantics.precision - 1; 66611819Sjulian significandParts()[0] = value; 66711819Sjulian normalize(rmNearestTiesToEven, lfExactlyZero); 66812632Sjulian} 66912632Sjulian 67011819SjulianAPFloat::APFloat(const fltSemantics &ourSemantics, 67111819Sjulian fltCategory ourCategory, bool negative, unsigned type) 67211819Sjulian{ 67311819Sjulian assertArithmeticOK(ourSemantics); 67411819Sjulian initialize(&ourSemantics); 67516178Sjulian category = ourCategory; 67616178Sjulian sign = negative; 67717265Sjulian if (category == fcNormal) 67818066Sjulian category = fcZero; 67917254Sjulian else if (ourCategory == fcNaN) 68017254Sjulian makeNaN(type); 68116178Sjulian} 68216178Sjulian 68352419SjulianAPFloat::APFloat(const fltSemantics &ourSemantics, const char *text) 68452419Sjulian{ 685164686Syar assertArithmeticOK(ourSemantics); 686164686Syar initialize(&ourSemantics); 687164686Syar convertFromString(text, rmNearestTiesToEven); 68852419Sjulian} 68952419Sjulian 6901590SrgrimesAPFloat::APFloat(const APFloat &rhs) 6911590Srgrimes{ 6921590Srgrimes initialize(rhs.semantics); 693102975Sdwmalone assign(rhs); 6941590Srgrimes} 6951590Srgrimes 696102975SdwmaloneAPFloat::~APFloat() 6971590Srgrimes{ 698102975Sdwmalone freeSignificand(); 699102975Sdwmalone} 70035308Sphk 7011590Srgrimes// Profile - This method 'profiles' an APFloat for use with FoldingSet. 70235308Sphkvoid APFloat::Profile(FoldingSetNodeID& ID) const { 703132803Sglebius ID.Add(bitcastToAPInt()); 704175217Sthompsa} 70593957Sru 70693957Sruunsigned int 707100161SkbyancAPFloat::partCount() const 70893957Sru{ 70935308Sphk return partCountForBits(semantics->precision + 1); 71093957Sru} 71135308Sphk 71235308Sphkunsigned int 71335308SphkAPFloat::semanticsPrecision(const fltSemantics &semantics) 7141590Srgrimes{ 7151590Srgrimes return semantics.precision; 7161590Srgrimes} 7171590Srgrimes 7181590Srgrimesconst integerPart * 7191590SrgrimesAPFloat::significandParts() const 720102975Sdwmalone{ 721102975Sdwmalone return const_cast<APFloat *>(this)->significandParts(); 7221590Srgrimes} 723102975Sdwmalone 7241590SrgrimesintegerPart * 725102975SdwmaloneAPFloat::significandParts() 726102975Sdwmalone{ 727102975Sdwmalone assert(category == fcNormal || category == fcNaN); 728102975Sdwmalone 7291590Srgrimes if(partCount() > 1) 730102975Sdwmalone return significand.parts; 7311590Srgrimes else 7321590Srgrimes return &significand.part; 7331590Srgrimes} 7341590Srgrimes 73597878Skbyancvoid 73697878SkbyancAPFloat::zeroSignificand() 7371590Srgrimes{ 7381590Srgrimes category = fcNormal; 7391590Srgrimes APInt::tcSet(significandParts(), 0, partCount()); 740102975Sdwmalone} 7411590Srgrimes 74297878Skbyanc/* Increment an fcNormal floating point number's significand. */ 74397878Skbyancvoid 7441590SrgrimesAPFloat::incrementSignificand() 74597878Skbyanc{ 74697878Skbyanc integerPart carry; 74797878Skbyanc 74897878Skbyanc carry = APInt::tcIncrement(significandParts(), partCount()); 74997878Skbyanc 75097878Skbyanc /* Our callers should never cause us to overflow. */ 75197878Skbyanc assert(carry == 0); 7521590Srgrimes} 7531590Srgrimes 7541590Srgrimes/* Add the significand of the RHS. Returns the carry flag. */ 7551590SrgrimesintegerPart 75697878SkbyancAPFloat::addSignificand(const APFloat &rhs) 7571590Srgrimes{ 7581590Srgrimes integerPart *parts; 7591590Srgrimes 76078314Sassar parts = significandParts(); 7611590Srgrimes 7621590Srgrimes assert(semantics == rhs.semantics); 76397878Skbyanc assert(exponent == rhs.exponent); 76497878Skbyanc 76513430Speter return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount()); 76635383Sphk} 7671590Srgrimes 76837452Sbde/* Subtract the significand of the RHS with a borrow flag. Returns 76935383Sphk the borrow flag. */ 77037452SbdeintegerPart 77137452SbdeAPFloat::subtractSignificand(const APFloat &rhs, integerPart borrow) 77235383Sphk{ 77337452Sbde integerPart *parts; 77497878Skbyanc 77597878Skbyanc parts = significandParts(); 77697878Skbyanc 77797878Skbyanc assert(semantics == rhs.semantics); 77883200Sru assert(exponent == rhs.exponent); 779186119Sqingli 78097878Skbyanc return APInt::tcSubtract(parts, rhs.significandParts(), borrow, 78183200Sru partCount()); 78278064Sume} 78397878Skbyanc 78478064Sume/* Multiply the significand of the RHS. If ADDEND is non-NULL, add it 78597878Skbyanc on to the full-precision result of the multiplication. Returns the 78678064Sume lost fraction. */ 78778064SumelostFraction 7881590SrgrimesAPFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend) 7891590Srgrimes{ 790176289Sjhb unsigned int omsb; // One, not zero, based MSB. 791176289Sjhb unsigned int partsCount, newPartsCount, precision; 792176289Sjhb integerPart *lhsSignificand; 793176289Sjhb integerPart scratch[4]; 794176289Sjhb integerPart *fullSignificand; 7951590Srgrimes lostFraction lost_fraction; 7961590Srgrimes bool ignored; 79797878Skbyanc 79835308Sphk assert(semantics == rhs.semantics); 7997642Sjkh 8005811Swollman precision = semantics->precision; 80135308Sphk newPartsCount = partCountForBits(precision * 2); 802160130Soleg 80397878Skbyanc if(newPartsCount > 4) 80478957Sru fullSignificand = new integerPart[newPartsCount]; 80578957Sru else 80635308Sphk fullSignificand = scratch; 8071590Srgrimes 8081590Srgrimes lhsSignificand = significandParts(); 8091590Srgrimes partsCount = partCount(); 8101590Srgrimes 8111590Srgrimes APInt::tcFullMultiply(fullSignificand, lhsSignificand, 812176099Smarius rhs.significandParts(), partsCount, partsCount); 8131590Srgrimes 814102975Sdwmalone lost_fraction = lfExactlyZero; 81574262Sbrian omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1; 8161590Srgrimes exponent += rhs.exponent; 8171590Srgrimes 8181590Srgrimes if(addend) { 81978238Sassar Significand savedSignificand = significand; 820176099Smarius const fltSemantics *savedSemantics = semantics; 8211590Srgrimes fltSemantics extendedSemantics; 8221590Srgrimes opStatus status; 82374262Sbrian unsigned int extendedPrecision; 8241590Srgrimes 8251590Srgrimes /* Normalize our MSB. */ 82636788Simp extendedPrecision = precision + precision - 1; 827183988Sdelphij if(omsb != extendedPrecision) 82836788Simp { 829175061Sobrien APInt::tcShiftLeft(fullSignificand, newPartsCount, 8301590Srgrimes extendedPrecision - omsb); 831176099Smarius exponent -= extendedPrecision - omsb; 8321590Srgrimes } 8331590Srgrimes 8341590Srgrimes /* Create new semantics. */ 8351590Srgrimes extendedSemantics = *semantics; 8361590Srgrimes extendedSemantics.precision = extendedPrecision; 837166711Sbms 838166711Sbms if(newPartsCount == 1) 839166711Sbms significand.part = fullSignificand[0]; 840166711Sbms else 841166711Sbms significand.parts = fullSignificand; 84213430Speter semantics = &extendedSemantics; 84313430Speter 844176194Smarius APFloat extendedAddend(*addend); 84513430Speter status = extendedAddend.convert(extendedSemantics, rmTowardZero, &ignored); 846102975Sdwmalone assert(status == opOK); 84713430Speter lost_fraction = addOrSubtractSignificand(extendedAddend, false); 848166711Sbms 84913430Speter /* Restore our state. */ 85013430Speter if(newPartsCount == 1) 85113430Speter fullSignificand[0] = significand.part; 85213430Speter significand = savedSignificand; 85313430Speter semantics = savedSemantics; 85413430Speter 855102975Sdwmalone omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1; 85613430Speter } 85713430Speter 85813430Speter exponent -= (precision - 1); 85913430Speter 86013430Speter if(omsb > precision) { 86113430Speter unsigned int bits, significantParts; 86213430Speter lostFraction lf; 86313430Speter 86413430Speter bits = omsb - precision; 86513430Speter significantParts = partCountForBits(omsb); 86613430Speter lf = shiftRight(fullSignificand, significantParts, bits); 86713430Speter lost_fraction = combineLostFractions(lf, lost_fraction); 86813430Speter exponent += bits; 86913430Speter } 87013430Speter 8711590Srgrimes APInt::tcAssign(lhsSignificand, fullSignificand, partsCount); 8721590Srgrimes 8731590Srgrimes if(newPartsCount > 4) 8741590Srgrimes delete [] fullSignificand; 8751590Srgrimes 876176099Smarius return lost_fraction; 8771590Srgrimes} 8781590Srgrimes 87974262Sbrian/* Multiply the significands of LHS and RHS to DST. */ 8801590SrgrimeslostFraction 881176194SmariusAPFloat::divideSignificand(const APFloat &rhs) 8821590Srgrimes{ 8831590Srgrimes unsigned int bit, i, partsCount; 88478238Sassar const integerPart *rhsSignificand; 88584803Sru integerPart *lhsSignificand, *dividend, *divisor; 88684803Sru integerPart scratch[4]; 8871590Srgrimes lostFraction lost_fraction; 88874262Sbrian 88913433Speter assert(semantics == rhs.semantics); 8901590Srgrimes 89184803Sru lhsSignificand = significandParts(); 892183988Sdelphij rhsSignificand = rhs.significandParts(); 89384803Sru partsCount = partCount(); 894176099Smarius 89577911Sru if(partsCount > 2) 89684803Sru dividend = new integerPart[partsCount * 2]; 8971590Srgrimes else 8981590Srgrimes dividend = scratch; 8991590Srgrimes 900166711Sbms divisor = dividend + partsCount; 901166711Sbms 90254263Sshin /* Copy the dividend and divisor as they will be modified in-place. */ 903217642Sume for(i = 0; i < partsCount; i++) { 904217642Sume dividend[i] = lhsSignificand[i]; 905217642Sume divisor[i] = rhsSignificand[i]; 906217642Sume lhsSignificand[i] = 0; 907217642Sume } 908217642Sume 909217642Sume exponent -= rhs.exponent; 910217642Sume 911217642Sume unsigned int precision = semantics->precision; 912217642Sume 913217642Sume /* Normalize the divisor. */ 914217642Sume bit = precision - APInt::tcMSB(divisor, partsCount) - 1; 915217642Sume if(bit) { 916217642Sume exponent += bit; 917217642Sume APInt::tcShiftLeft(divisor, partsCount, bit); 918217642Sume } 919217642Sume 920217642Sume /* Normalize the dividend. */ 921217642Sume bit = precision - APInt::tcMSB(dividend, partsCount) - 1; 922102975Sdwmalone if(bit) { 92378314Sassar exponent -= bit; 92454263Sshin APInt::tcShiftLeft(dividend, partsCount, bit); 92574262Sbrian } 92654263Sshin 92754263Sshin /* Ensure the dividend >= divisor initially for the loop below. 928146187Sume Incidentally, this means that the division loop below is 92954263Sshin guaranteed to set the integer bit to one. */ 93054263Sshin if(APInt::tcCompare(dividend, divisor, partsCount) < 0) { 93154263Sshin exponent--; 93254263Sshin APInt::tcShiftLeft(dividend, partsCount, 1); 93362584Sitojun assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0); 93462584Sitojun } 93562584Sitojun 93654263Sshin /* Long division. */ 93754263Sshin for(bit = precision; bit; bit -= 1) { 93854263Sshin if(APInt::tcCompare(dividend, divisor, partsCount) >= 0) { 93954263Sshin APInt::tcSubtract(dividend, divisor, 0, partsCount); 94054263Sshin APInt::tcSetBit(lhsSignificand, bit - 1); 94154263Sshin } 94254263Sshin 94354263Sshin APInt::tcShiftLeft(dividend, partsCount, 1); 94454263Sshin } 94554263Sshin 94654263Sshin /* Figure out the lost fraction. */ 94754263Sshin int cmp = APInt::tcCompare(dividend, divisor, partsCount); 94854263Sshin 94954263Sshin if(cmp > 0) 95054263Sshin lost_fraction = lfMoreThanHalf; 95154263Sshin else if(cmp == 0) 95254263Sshin lost_fraction = lfExactlyHalf; 95354263Sshin else if(APInt::tcIsZero(dividend, partsCount)) 95454263Sshin lost_fraction = lfExactlyZero; 95554263Sshin else 95654263Sshin lost_fraction = lfLessThanHalf; 95754263Sshin 95854263Sshin if(partsCount > 2) 95954263Sshin delete [] dividend; 96054263Sshin 96154263Sshin return lost_fraction; 96254263Sshin} 96354263Sshin 96454263Sshinunsigned int 96554263SshinAPFloat::significandMSB() const 96654263Sshin{ 96754263Sshin return APInt::tcMSB(significandParts(), partCount()); 96854263Sshin} 96954263Sshin 97054263Sshinunsigned int 97154263SshinAPFloat::significandLSB() const 97254263Sshin{ 97378238Sassar return APInt::tcLSB(significandParts(), partCount()); 97454263Sshin} 97554263Sshin 97654263Sshin/* Note that a zero result is NOT normalized to fcZero. */ 97754263SshinlostFraction 97878238SassarAPFloat::shiftSignificandRight(unsigned int bits) 97954263Sshin{ 98054263Sshin /* Our exponent should not overflow. */ 98154263Sshin assert((exponent_t) (exponent + bits) >= exponent); 98254263Sshin 98354263Sshin exponent += bits; 98454263Sshin 98578314Sassar return shiftRight(significandParts(), partCount(), bits); 98654263Sshin} 98774262Sbrian 988146187Sume/* Shift the significand left BITS bits, subtract BITS from its exponent. */ 98962584Sitojunvoid 990160789SyarAPFloat::shiftSignificandLeft(unsigned int bits) 99154263Sshin{ 992160789Syar assert(bits < semantics->precision); 993160789Syar 99462584Sitojun if(bits) { 99562584Sitojun unsigned int partsCount = partCount(); 99662584Sitojun 99778238Sassar APInt::tcShiftLeft(significandParts(), partsCount, bits); 99854263Sshin exponent -= bits; 99954263Sshin 100062584Sitojun assert(!APInt::tcIsZero(significandParts(), partsCount)); 100162584Sitojun } 100254263Sshin} 100354263Sshin 100454263SshinAPFloat::cmpResult 100554263SshinAPFloat::compareAbsoluteValue(const APFloat &rhs) const 100654263Sshin{ 10071590Srgrimes int compare; 10081590Srgrimes 10091590Srgrimes assert(semantics == rhs.semantics); 10101590Srgrimes assert(category == fcNormal); 101178958Sru assert(rhs.category == fcNormal); 10121590Srgrimes 10131590Srgrimes compare = exponent - rhs.exponent; 101478958Sru 10151590Srgrimes /* If exponents are equal, do an unsigned bignum comparison of the 101678958Sru significands. */ 10171590Srgrimes if(compare == 0) 10181590Srgrimes compare = APInt::tcCompare(significandParts(), rhs.significandParts(), 10191590Srgrimes partCount()); 102078958Sru 102178958Sru if(compare > 0) 102278958Sru return cmpGreaterThan; 102378958Sru else if(compare < 0) 102478958Sru return cmpLessThan; 102578958Sru else 10261590Srgrimes return cmpEqual; 102778659Sru} 102878659Sru 102978659Sru/* Handle overflow. Sign is preserved. We either become infinity or 103078659Sru the largest finite number. */ 1031198118SrwatsonAPFloat::opStatus 1032198118SrwatsonAPFloat::handleOverflow(roundingMode rounding_mode) 1033198118Srwatson{ 1034198118Srwatson /* Infinity? */ 1035198118Srwatson if(rounding_mode == rmNearestTiesToEven 103678659Sru || rounding_mode == rmNearestTiesToAway 103778958Sru || (rounding_mode == rmTowardPositive && !sign) 103878958Sru || (rounding_mode == rmTowardNegative && sign)) 103978958Sru { 104078958Sru category = fcInfinity; 10411590Srgrimes return (opStatus) (opOverflow | opInexact); 104211819Sjulian } 104311819Sjulian 104478314Sassar /* Otherwise we become the largest finite number. */ 104511819Sjulian category = fcNormal; 104611819Sjulian exponent = semantics->maxExponent; 104735308Sphk APInt::tcSetLeastSignificantBits(significandParts(), partCount(), 1048102975Sdwmalone semantics->precision); 1049102975Sdwmalone 1050102975Sdwmalone return opInexact; 105111819Sjulian} 105235308Sphk 105311819Sjulian/* Returns TRUE if, when truncating the current number, with BIT the 105411819Sjulian new LSB, with the given lost fraction and rounding mode, the result 105511819Sjulian would need to be rounded away from zero (i.e., by increasing the 105611819Sjulian signficand). This routine must work for fcZero of both signs, and 105711819Sjulian fcNormal numbers. */ 105811819Sjulianbool 105911819SjulianAPFloat::roundAwayFromZero(roundingMode rounding_mode, 106035308Sphk lostFraction lost_fraction, 106135308Sphk unsigned int bit) const 106235308Sphk{ 106335308Sphk /* NaNs and infinities should not have lost fractions. */ 106411819Sjulian assert(category == fcNormal || category == fcZero); 106511819Sjulian 106611819Sjulian /* Current callers never pass this so we don't handle it. */ 106711819Sjulian assert(lost_fraction != lfExactlyZero); 106811819Sjulian 106911819Sjulian switch (rounding_mode) { 107025654Sjhay default: 107111819Sjulian assert(0); 107211819Sjulian 107311819Sjulian case rmNearestTiesToAway: 107411819Sjulian return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf; 107511819Sjulian 107611819Sjulian case rmNearestTiesToEven: 107711819Sjulian if(lost_fraction == lfMoreThanHalf) 107811819Sjulian return true; 107911819Sjulian 108011819Sjulian /* Our zeroes don't have a significand to test. */ 108111819Sjulian if(lost_fraction == lfExactlyHalf && category != fcZero) 108211819Sjulian return APInt::tcExtractBit(significandParts(), bit); 108325654Sjhay 108411819Sjulian return false; 108511819Sjulian 108611819Sjulian case rmTowardZero: 108711819Sjulian return false; 108811819Sjulian 108911819Sjulian case rmTowardPositive: 109011819Sjulian return sign == false; 109111819Sjulian 109211819Sjulian case rmTowardNegative: 109311819Sjulian return sign == true; 109411819Sjulian } 109511819Sjulian} 109611819Sjulian 109735308SphkAPFloat::opStatus 109835308SphkAPFloat::normalize(roundingMode rounding_mode, 1099175061Sobrien lostFraction lost_fraction) 110036788Simp{ 110136788Simp unsigned int omsb; /* One, not zero, based MSB. */ 1102175061Sobrien int exponentChange; 110336788Simp 110436788Simp if(category != fcNormal) 110511819Sjulian return opOK; 110611819Sjulian 110711819Sjulian /* Before rounding normalize the exponent of fcNormal numbers. */ 110836788Simp omsb = significandMSB() + 1; 110911819Sjulian 111011819Sjulian if(omsb) { 111111819Sjulian /* OMSB is numbered from 1. We want to place it in the integer 111211819Sjulian bit numbered PRECISON if possible, with a compensating change in 111378314Sassar the exponent. */ 111411819Sjulian exponentChange = omsb - semantics->precision; 1115102975Sdwmalone 111611819Sjulian /* If the resulting exponent is too high, overflow according to 111735308Sphk the rounding mode. */ 111811819Sjulian if(exponent + exponentChange > semantics->maxExponent) 111911819Sjulian return handleOverflow(rounding_mode); 112011819Sjulian 112111819Sjulian /* Subnormal numbers have exponent minExponent, and their MSB 112211819Sjulian is forced based on that. */ 112311819Sjulian if(exponent + exponentChange < semantics->minExponent) 112411819Sjulian exponentChange = semantics->minExponent - exponent; 112511819Sjulian 112611819Sjulian /* Shifting left is easy as we don't lose precision. */ 112711819Sjulian if(exponentChange < 0) { 112811819Sjulian assert(lost_fraction == lfExactlyZero); 112911819Sjulian 11301590Srgrimes shiftSignificandLeft(-exponentChange); 113178314Sassar 11321590Srgrimes return opOK; 1133102975Sdwmalone } 11341590Srgrimes 113535308Sphk if(exponentChange > 0) { 113635308Sphk lostFraction lf; 113735308Sphk 113835308Sphk /* Shift right and capture any new lost fraction. */ 113935308Sphk lf = shiftSignificandRight(exponentChange); 114035308Sphk 114135308Sphk lost_fraction = combineLostFractions(lf, lost_fraction); 114235308Sphk 114335308Sphk /* Keep OMSB up-to-date. */ 114435308Sphk if(omsb > (unsigned) exponentChange) 114535308Sphk omsb -= exponentChange; 114635308Sphk else 11471590Srgrimes omsb = 0; 1148 } 1149 } 1150 1151 /* Now round the number according to rounding_mode given the lost 1152 fraction. */ 1153 1154 /* As specified in IEEE 754, since we do not trap we do not report 1155 underflow for exact results. */ 1156 if(lost_fraction == lfExactlyZero) { 1157 /* Canonicalize zeroes. */ 1158 if(omsb == 0) 1159 category = fcZero; 1160 1161 return opOK; 1162 } 1163 1164 /* Increment the significand if we're rounding away from zero. */ 1165 if(roundAwayFromZero(rounding_mode, lost_fraction, 0)) { 1166 if(omsb == 0) 1167 exponent = semantics->minExponent; 1168 1169 incrementSignificand(); 1170 omsb = significandMSB() + 1; 1171 1172 /* Did the significand increment overflow? */ 1173 if(omsb == (unsigned) semantics->precision + 1) { 1174 /* Renormalize by incrementing the exponent and shifting our 1175 significand right one. However if we already have the 1176 maximum exponent we overflow to infinity. */ 1177 if(exponent == semantics->maxExponent) { 1178 category = fcInfinity; 1179 1180 return (opStatus) (opOverflow | opInexact); 1181 } 1182 1183 shiftSignificandRight(1); 1184 1185 return opInexact; 1186 } 1187 } 1188 1189 /* The normal case - we were and are not denormal, and any 1190 significand increment above didn't overflow. */ 1191 if(omsb == semantics->precision) 1192 return opInexact; 1193 1194 /* We have a non-zero denormal. */ 1195 assert(omsb < semantics->precision); 1196 1197 /* Canonicalize zeroes. */ 1198 if(omsb == 0) 1199 category = fcZero; 1200 1201 /* The fcZero case is a denormal that underflowed to zero. */ 1202 return (opStatus) (opUnderflow | opInexact); 1203} 1204 1205APFloat::opStatus 1206APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract) 1207{ 1208 switch (convolve(category, rhs.category)) { 1209 default: 1210 assert(0); 1211 1212 case convolve(fcNaN, fcZero): 1213 case convolve(fcNaN, fcNormal): 1214 case convolve(fcNaN, fcInfinity): 1215 case convolve(fcNaN, fcNaN): 1216 case convolve(fcNormal, fcZero): 1217 case convolve(fcInfinity, fcNormal): 1218 case convolve(fcInfinity, fcZero): 1219 return opOK; 1220 1221 case convolve(fcZero, fcNaN): 1222 case convolve(fcNormal, fcNaN): 1223 case convolve(fcInfinity, fcNaN): 1224 category = fcNaN; 1225 copySignificand(rhs); 1226 return opOK; 1227 1228 case convolve(fcNormal, fcInfinity): 1229 case convolve(fcZero, fcInfinity): 1230 category = fcInfinity; 1231 sign = rhs.sign ^ subtract; 1232 return opOK; 1233 1234 case convolve(fcZero, fcNormal): 1235 assign(rhs); 1236 sign = rhs.sign ^ subtract; 1237 return opOK; 1238 1239 case convolve(fcZero, fcZero): 1240 /* Sign depends on rounding mode; handled by caller. */ 1241 return opOK; 1242 1243 case convolve(fcInfinity, fcInfinity): 1244 /* Differently signed infinities can only be validly 1245 subtracted. */ 1246 if(((sign ^ rhs.sign)!=0) != subtract) { 1247 makeNaN(); 1248 return opInvalidOp; 1249 } 1250 1251 return opOK; 1252 1253 case convolve(fcNormal, fcNormal): 1254 return opDivByZero; 1255 } 1256} 1257 1258/* Add or subtract two normal numbers. */ 1259lostFraction 1260APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract) 1261{ 1262 integerPart carry; 1263 lostFraction lost_fraction; 1264 int bits; 1265 1266 /* Determine if the operation on the absolute values is effectively 1267 an addition or subtraction. */ 1268 subtract ^= (sign ^ rhs.sign) ? true : false; 1269 1270 /* Are we bigger exponent-wise than the RHS? */ 1271 bits = exponent - rhs.exponent; 1272 1273 /* Subtraction is more subtle than one might naively expect. */ 1274 if(subtract) { 1275 APFloat temp_rhs(rhs); 1276 bool reverse; 1277 1278 if (bits == 0) { 1279 reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan; 1280 lost_fraction = lfExactlyZero; 1281 } else if (bits > 0) { 1282 lost_fraction = temp_rhs.shiftSignificandRight(bits - 1); 1283 shiftSignificandLeft(1); 1284 reverse = false; 1285 } else { 1286 lost_fraction = shiftSignificandRight(-bits - 1); 1287 temp_rhs.shiftSignificandLeft(1); 1288 reverse = true; 1289 } 1290 1291 if (reverse) { 1292 carry = temp_rhs.subtractSignificand 1293 (*this, lost_fraction != lfExactlyZero); 1294 copySignificand(temp_rhs); 1295 sign = !sign; 1296 } else { 1297 carry = subtractSignificand 1298 (temp_rhs, lost_fraction != lfExactlyZero); 1299 } 1300 1301 /* Invert the lost fraction - it was on the RHS and 1302 subtracted. */ 1303 if(lost_fraction == lfLessThanHalf) 1304 lost_fraction = lfMoreThanHalf; 1305 else if(lost_fraction == lfMoreThanHalf) 1306 lost_fraction = lfLessThanHalf; 1307 1308 /* The code above is intended to ensure that no borrow is 1309 necessary. */ 1310 assert(!carry); 1311 } else { 1312 if(bits > 0) { 1313 APFloat temp_rhs(rhs); 1314 1315 lost_fraction = temp_rhs.shiftSignificandRight(bits); 1316 carry = addSignificand(temp_rhs); 1317 } else { 1318 lost_fraction = shiftSignificandRight(-bits); 1319 carry = addSignificand(rhs); 1320 } 1321 1322 /* We have a guard bit; generating a carry cannot happen. */ 1323 assert(!carry); 1324 } 1325 1326 return lost_fraction; 1327} 1328 1329APFloat::opStatus 1330APFloat::multiplySpecials(const APFloat &rhs) 1331{ 1332 switch (convolve(category, rhs.category)) { 1333 default: 1334 assert(0); 1335 1336 case convolve(fcNaN, fcZero): 1337 case convolve(fcNaN, fcNormal): 1338 case convolve(fcNaN, fcInfinity): 1339 case convolve(fcNaN, fcNaN): 1340 return opOK; 1341 1342 case convolve(fcZero, fcNaN): 1343 case convolve(fcNormal, fcNaN): 1344 case convolve(fcInfinity, fcNaN): 1345 category = fcNaN; 1346 copySignificand(rhs); 1347 return opOK; 1348 1349 case convolve(fcNormal, fcInfinity): 1350 case convolve(fcInfinity, fcNormal): 1351 case convolve(fcInfinity, fcInfinity): 1352 category = fcInfinity; 1353 return opOK; 1354 1355 case convolve(fcZero, fcNormal): 1356 case convolve(fcNormal, fcZero): 1357 case convolve(fcZero, fcZero): 1358 category = fcZero; 1359 return opOK; 1360 1361 case convolve(fcZero, fcInfinity): 1362 case convolve(fcInfinity, fcZero): 1363 makeNaN(); 1364 return opInvalidOp; 1365 1366 case convolve(fcNormal, fcNormal): 1367 return opOK; 1368 } 1369} 1370 1371APFloat::opStatus 1372APFloat::divideSpecials(const APFloat &rhs) 1373{ 1374 switch (convolve(category, rhs.category)) { 1375 default: 1376 assert(0); 1377 1378 case convolve(fcNaN, fcZero): 1379 case convolve(fcNaN, fcNormal): 1380 case convolve(fcNaN, fcInfinity): 1381 case convolve(fcNaN, fcNaN): 1382 case convolve(fcInfinity, fcZero): 1383 case convolve(fcInfinity, fcNormal): 1384 case convolve(fcZero, fcInfinity): 1385 case convolve(fcZero, fcNormal): 1386 return opOK; 1387 1388 case convolve(fcZero, fcNaN): 1389 case convolve(fcNormal, fcNaN): 1390 case convolve(fcInfinity, fcNaN): 1391 category = fcNaN; 1392 copySignificand(rhs); 1393 return opOK; 1394 1395 case convolve(fcNormal, fcInfinity): 1396 category = fcZero; 1397 return opOK; 1398 1399 case convolve(fcNormal, fcZero): 1400 category = fcInfinity; 1401 return opDivByZero; 1402 1403 case convolve(fcInfinity, fcInfinity): 1404 case convolve(fcZero, fcZero): 1405 makeNaN(); 1406 return opInvalidOp; 1407 1408 case convolve(fcNormal, fcNormal): 1409 return opOK; 1410 } 1411} 1412 1413APFloat::opStatus 1414APFloat::modSpecials(const APFloat &rhs) 1415{ 1416 switch (convolve(category, rhs.category)) { 1417 default: 1418 assert(0); 1419 1420 case convolve(fcNaN, fcZero): 1421 case convolve(fcNaN, fcNormal): 1422 case convolve(fcNaN, fcInfinity): 1423 case convolve(fcNaN, fcNaN): 1424 case convolve(fcZero, fcInfinity): 1425 case convolve(fcZero, fcNormal): 1426 case convolve(fcNormal, fcInfinity): 1427 return opOK; 1428 1429 case convolve(fcZero, fcNaN): 1430 case convolve(fcNormal, fcNaN): 1431 case convolve(fcInfinity, fcNaN): 1432 category = fcNaN; 1433 copySignificand(rhs); 1434 return opOK; 1435 1436 case convolve(fcNormal, fcZero): 1437 case convolve(fcInfinity, fcZero): 1438 case convolve(fcInfinity, fcNormal): 1439 case convolve(fcInfinity, fcInfinity): 1440 case convolve(fcZero, fcZero): 1441 makeNaN(); 1442 return opInvalidOp; 1443 1444 case convolve(fcNormal, fcNormal): 1445 return opOK; 1446 } 1447} 1448 1449/* Change sign. */ 1450void 1451APFloat::changeSign() 1452{ 1453 /* Look mummy, this one's easy. */ 1454 sign = !sign; 1455} 1456 1457void 1458APFloat::clearSign() 1459{ 1460 /* So is this one. */ 1461 sign = 0; 1462} 1463 1464void 1465APFloat::copySign(const APFloat &rhs) 1466{ 1467 /* And this one. */ 1468 sign = rhs.sign; 1469} 1470 1471/* Normalized addition or subtraction. */ 1472APFloat::opStatus 1473APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode, 1474 bool subtract) 1475{ 1476 opStatus fs; 1477 1478 assertArithmeticOK(*semantics); 1479 1480 fs = addOrSubtractSpecials(rhs, subtract); 1481 1482 /* This return code means it was not a simple case. */ 1483 if(fs == opDivByZero) { 1484 lostFraction lost_fraction; 1485 1486 lost_fraction = addOrSubtractSignificand(rhs, subtract); 1487 fs = normalize(rounding_mode, lost_fraction); 1488 1489 /* Can only be zero if we lost no fraction. */ 1490 assert(category != fcZero || lost_fraction == lfExactlyZero); 1491 } 1492 1493 /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a 1494 positive zero unless rounding to minus infinity, except that 1495 adding two like-signed zeroes gives that zero. */ 1496 if(category == fcZero) { 1497 if(rhs.category != fcZero || (sign == rhs.sign) == subtract) 1498 sign = (rounding_mode == rmTowardNegative); 1499 } 1500 1501 return fs; 1502} 1503 1504/* Normalized addition. */ 1505APFloat::opStatus 1506APFloat::add(const APFloat &rhs, roundingMode rounding_mode) 1507{ 1508 return addOrSubtract(rhs, rounding_mode, false); 1509} 1510 1511/* Normalized subtraction. */ 1512APFloat::opStatus 1513APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode) 1514{ 1515 return addOrSubtract(rhs, rounding_mode, true); 1516} 1517 1518/* Normalized multiply. */ 1519APFloat::opStatus 1520APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode) 1521{ 1522 opStatus fs; 1523 1524 assertArithmeticOK(*semantics); 1525 sign ^= rhs.sign; 1526 fs = multiplySpecials(rhs); 1527 1528 if(category == fcNormal) { 1529 lostFraction lost_fraction = multiplySignificand(rhs, 0); 1530 fs = normalize(rounding_mode, lost_fraction); 1531 if(lost_fraction != lfExactlyZero) 1532 fs = (opStatus) (fs | opInexact); 1533 } 1534 1535 return fs; 1536} 1537 1538/* Normalized divide. */ 1539APFloat::opStatus 1540APFloat::divide(const APFloat &rhs, roundingMode rounding_mode) 1541{ 1542 opStatus fs; 1543 1544 assertArithmeticOK(*semantics); 1545 sign ^= rhs.sign; 1546 fs = divideSpecials(rhs); 1547 1548 if(category == fcNormal) { 1549 lostFraction lost_fraction = divideSignificand(rhs); 1550 fs = normalize(rounding_mode, lost_fraction); 1551 if(lost_fraction != lfExactlyZero) 1552 fs = (opStatus) (fs | opInexact); 1553 } 1554 1555 return fs; 1556} 1557 1558/* Normalized remainder. This is not currently correct in all cases. */ 1559APFloat::opStatus 1560APFloat::remainder(const APFloat &rhs) 1561{ 1562 opStatus fs; 1563 APFloat V = *this; 1564 unsigned int origSign = sign; 1565 1566 assertArithmeticOK(*semantics); 1567 fs = V.divide(rhs, rmNearestTiesToEven); 1568 if (fs == opDivByZero) 1569 return fs; 1570 1571 int parts = partCount(); 1572 integerPart *x = new integerPart[parts]; 1573 bool ignored; 1574 fs = V.convertToInteger(x, parts * integerPartWidth, true, 1575 rmNearestTiesToEven, &ignored); 1576 if (fs==opInvalidOp) 1577 return fs; 1578 1579 fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true, 1580 rmNearestTiesToEven); 1581 assert(fs==opOK); // should always work 1582 1583 fs = V.multiply(rhs, rmNearestTiesToEven); 1584 assert(fs==opOK || fs==opInexact); // should not overflow or underflow 1585 1586 fs = subtract(V, rmNearestTiesToEven); 1587 assert(fs==opOK || fs==opInexact); // likewise 1588 1589 if (isZero()) 1590 sign = origSign; // IEEE754 requires this 1591 delete[] x; 1592 return fs; 1593} 1594 1595/* Normalized llvm frem (C fmod). 1596 This is not currently correct in all cases. */ 1597APFloat::opStatus 1598APFloat::mod(const APFloat &rhs, roundingMode rounding_mode) 1599{ 1600 opStatus fs; 1601 assertArithmeticOK(*semantics); 1602 fs = modSpecials(rhs); 1603 1604 if (category == fcNormal && rhs.category == fcNormal) { 1605 APFloat V = *this; 1606 unsigned int origSign = sign; 1607 1608 fs = V.divide(rhs, rmNearestTiesToEven); 1609 if (fs == opDivByZero) 1610 return fs; 1611 1612 int parts = partCount(); 1613 integerPart *x = new integerPart[parts]; 1614 bool ignored; 1615 fs = V.convertToInteger(x, parts * integerPartWidth, true, 1616 rmTowardZero, &ignored); 1617 if (fs==opInvalidOp) 1618 return fs; 1619 1620 fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true, 1621 rmNearestTiesToEven); 1622 assert(fs==opOK); // should always work 1623 1624 fs = V.multiply(rhs, rounding_mode); 1625 assert(fs==opOK || fs==opInexact); // should not overflow or underflow 1626 1627 fs = subtract(V, rounding_mode); 1628 assert(fs==opOK || fs==opInexact); // likewise 1629 1630 if (isZero()) 1631 sign = origSign; // IEEE754 requires this 1632 delete[] x; 1633 } 1634 return fs; 1635} 1636 1637/* Normalized fused-multiply-add. */ 1638APFloat::opStatus 1639APFloat::fusedMultiplyAdd(const APFloat &multiplicand, 1640 const APFloat &addend, 1641 roundingMode rounding_mode) 1642{ 1643 opStatus fs; 1644 1645 assertArithmeticOK(*semantics); 1646 1647 /* Post-multiplication sign, before addition. */ 1648 sign ^= multiplicand.sign; 1649 1650 /* If and only if all arguments are normal do we need to do an 1651 extended-precision calculation. */ 1652 if(category == fcNormal 1653 && multiplicand.category == fcNormal 1654 && addend.category == fcNormal) { 1655 lostFraction lost_fraction; 1656 1657 lost_fraction = multiplySignificand(multiplicand, &addend); 1658 fs = normalize(rounding_mode, lost_fraction); 1659 if(lost_fraction != lfExactlyZero) 1660 fs = (opStatus) (fs | opInexact); 1661 1662 /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a 1663 positive zero unless rounding to minus infinity, except that 1664 adding two like-signed zeroes gives that zero. */ 1665 if(category == fcZero && sign != addend.sign) 1666 sign = (rounding_mode == rmTowardNegative); 1667 } else { 1668 fs = multiplySpecials(multiplicand); 1669 1670 /* FS can only be opOK or opInvalidOp. There is no more work 1671 to do in the latter case. The IEEE-754R standard says it is 1672 implementation-defined in this case whether, if ADDEND is a 1673 quiet NaN, we raise invalid op; this implementation does so. 1674 1675 If we need to do the addition we can do so with normal 1676 precision. */ 1677 if(fs == opOK) 1678 fs = addOrSubtract(addend, rounding_mode, false); 1679 } 1680 1681 return fs; 1682} 1683 1684/* Comparison requires normalized numbers. */ 1685APFloat::cmpResult 1686APFloat::compare(const APFloat &rhs) const 1687{ 1688 cmpResult result; 1689 1690 assertArithmeticOK(*semantics); 1691 assert(semantics == rhs.semantics); 1692 1693 switch (convolve(category, rhs.category)) { 1694 default: 1695 assert(0); 1696 1697 case convolve(fcNaN, fcZero): 1698 case convolve(fcNaN, fcNormal): 1699 case convolve(fcNaN, fcInfinity): 1700 case convolve(fcNaN, fcNaN): 1701 case convolve(fcZero, fcNaN): 1702 case convolve(fcNormal, fcNaN): 1703 case convolve(fcInfinity, fcNaN): 1704 return cmpUnordered; 1705 1706 case convolve(fcInfinity, fcNormal): 1707 case convolve(fcInfinity, fcZero): 1708 case convolve(fcNormal, fcZero): 1709 if(sign) 1710 return cmpLessThan; 1711 else 1712 return cmpGreaterThan; 1713 1714 case convolve(fcNormal, fcInfinity): 1715 case convolve(fcZero, fcInfinity): 1716 case convolve(fcZero, fcNormal): 1717 if(rhs.sign) 1718 return cmpGreaterThan; 1719 else 1720 return cmpLessThan; 1721 1722 case convolve(fcInfinity, fcInfinity): 1723 if(sign == rhs.sign) 1724 return cmpEqual; 1725 else if(sign) 1726 return cmpLessThan; 1727 else 1728 return cmpGreaterThan; 1729 1730 case convolve(fcZero, fcZero): 1731 return cmpEqual; 1732 1733 case convolve(fcNormal, fcNormal): 1734 break; 1735 } 1736 1737 /* Two normal numbers. Do they have the same sign? */ 1738 if(sign != rhs.sign) { 1739 if(sign) 1740 result = cmpLessThan; 1741 else 1742 result = cmpGreaterThan; 1743 } else { 1744 /* Compare absolute values; invert result if negative. */ 1745 result = compareAbsoluteValue(rhs); 1746 1747 if(sign) { 1748 if(result == cmpLessThan) 1749 result = cmpGreaterThan; 1750 else if(result == cmpGreaterThan) 1751 result = cmpLessThan; 1752 } 1753 } 1754 1755 return result; 1756} 1757 1758/// APFloat::convert - convert a value of one floating point type to another. 1759/// The return value corresponds to the IEEE754 exceptions. *losesInfo 1760/// records whether the transformation lost information, i.e. whether 1761/// converting the result back to the original type will produce the 1762/// original value (this is almost the same as return value==fsOK, but there 1763/// are edge cases where this is not so). 1764 1765APFloat::opStatus 1766APFloat::convert(const fltSemantics &toSemantics, 1767 roundingMode rounding_mode, bool *losesInfo) 1768{ 1769 lostFraction lostFraction; 1770 unsigned int newPartCount, oldPartCount; 1771 opStatus fs; 1772 1773 assertArithmeticOK(*semantics); 1774 assertArithmeticOK(toSemantics); 1775 lostFraction = lfExactlyZero; 1776 newPartCount = partCountForBits(toSemantics.precision + 1); 1777 oldPartCount = partCount(); 1778 1779 /* Handle storage complications. If our new form is wider, 1780 re-allocate our bit pattern into wider storage. If it is 1781 narrower, we ignore the excess parts, but if narrowing to a 1782 single part we need to free the old storage. 1783 Be careful not to reference significandParts for zeroes 1784 and infinities, since it aborts. */ 1785 if (newPartCount > oldPartCount) { 1786 integerPart *newParts; 1787 newParts = new integerPart[newPartCount]; 1788 APInt::tcSet(newParts, 0, newPartCount); 1789 if (category==fcNormal || category==fcNaN) 1790 APInt::tcAssign(newParts, significandParts(), oldPartCount); 1791 freeSignificand(); 1792 significand.parts = newParts; 1793 } else if (newPartCount < oldPartCount) { 1794 /* Capture any lost fraction through truncation of parts so we get 1795 correct rounding whilst normalizing. */ 1796 if (category==fcNormal) 1797 lostFraction = lostFractionThroughTruncation 1798 (significandParts(), oldPartCount, toSemantics.precision); 1799 if (newPartCount == 1) { 1800 integerPart newPart = 0; 1801 if (category==fcNormal || category==fcNaN) 1802 newPart = significandParts()[0]; 1803 freeSignificand(); 1804 significand.part = newPart; 1805 } 1806 } 1807 1808 if(category == fcNormal) { 1809 /* Re-interpret our bit-pattern. */ 1810 exponent += toSemantics.precision - semantics->precision; 1811 semantics = &toSemantics; 1812 fs = normalize(rounding_mode, lostFraction); 1813 *losesInfo = (fs != opOK); 1814 } else if (category == fcNaN) { 1815 int shift = toSemantics.precision - semantics->precision; 1816 // Do this now so significandParts gets the right answer 1817 const fltSemantics *oldSemantics = semantics; 1818 semantics = &toSemantics; 1819 *losesInfo = false; 1820 // No normalization here, just truncate 1821 if (shift>0) 1822 APInt::tcShiftLeft(significandParts(), newPartCount, shift); 1823 else if (shift < 0) { 1824 unsigned ushift = -shift; 1825 // Figure out if we are losing information. This happens 1826 // if are shifting out something other than 0s, or if the x87 long 1827 // double input did not have its integer bit set (pseudo-NaN), or if the 1828 // x87 long double input did not have its QNan bit set (because the x87 1829 // hardware sets this bit when converting a lower-precision NaN to 1830 // x87 long double). 1831 if (APInt::tcLSB(significandParts(), newPartCount) < ushift) 1832 *losesInfo = true; 1833 if (oldSemantics == &APFloat::x87DoubleExtended && 1834 (!(*significandParts() & 0x8000000000000000ULL) || 1835 !(*significandParts() & 0x4000000000000000ULL))) 1836 *losesInfo = true; 1837 APInt::tcShiftRight(significandParts(), newPartCount, ushift); 1838 } 1839 // gcc forces the Quiet bit on, which means (float)(double)(float_sNan) 1840 // does not give you back the same bits. This is dubious, and we 1841 // don't currently do it. You're really supposed to get 1842 // an invalid operation signal at runtime, but nobody does that. 1843 fs = opOK; 1844 } else { 1845 semantics = &toSemantics; 1846 fs = opOK; 1847 *losesInfo = false; 1848 } 1849 1850 return fs; 1851} 1852 1853/* Convert a floating point number to an integer according to the 1854 rounding mode. If the rounded integer value is out of range this 1855 returns an invalid operation exception and the contents of the 1856 destination parts are unspecified. If the rounded value is in 1857 range but the floating point number is not the exact integer, the C 1858 standard doesn't require an inexact exception to be raised. IEEE 1859 854 does require it so we do that. 1860 1861 Note that for conversions to integer type the C standard requires 1862 round-to-zero to always be used. */ 1863APFloat::opStatus 1864APFloat::convertToSignExtendedInteger(integerPart *parts, unsigned int width, 1865 bool isSigned, 1866 roundingMode rounding_mode, 1867 bool *isExact) const 1868{ 1869 lostFraction lost_fraction; 1870 const integerPart *src; 1871 unsigned int dstPartsCount, truncatedBits; 1872 1873 assertArithmeticOK(*semantics); 1874 1875 *isExact = false; 1876 1877 /* Handle the three special cases first. */ 1878 if(category == fcInfinity || category == fcNaN) 1879 return opInvalidOp; 1880 1881 dstPartsCount = partCountForBits(width); 1882 1883 if(category == fcZero) { 1884 APInt::tcSet(parts, 0, dstPartsCount); 1885 // Negative zero can't be represented as an int. 1886 *isExact = !sign; 1887 return opOK; 1888 } 1889 1890 src = significandParts(); 1891 1892 /* Step 1: place our absolute value, with any fraction truncated, in 1893 the destination. */ 1894 if (exponent < 0) { 1895 /* Our absolute value is less than one; truncate everything. */ 1896 APInt::tcSet(parts, 0, dstPartsCount); 1897 /* For exponent -1 the integer bit represents .5, look at that. 1898 For smaller exponents leftmost truncated bit is 0. */ 1899 truncatedBits = semantics->precision -1U - exponent; 1900 } else { 1901 /* We want the most significant (exponent + 1) bits; the rest are 1902 truncated. */ 1903 unsigned int bits = exponent + 1U; 1904 1905 /* Hopelessly large in magnitude? */ 1906 if (bits > width) 1907 return opInvalidOp; 1908 1909 if (bits < semantics->precision) { 1910 /* We truncate (semantics->precision - bits) bits. */ 1911 truncatedBits = semantics->precision - bits; 1912 APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits); 1913 } else { 1914 /* We want at least as many bits as are available. */ 1915 APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0); 1916 APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision); 1917 truncatedBits = 0; 1918 } 1919 } 1920 1921 /* Step 2: work out any lost fraction, and increment the absolute 1922 value if we would round away from zero. */ 1923 if (truncatedBits) { 1924 lost_fraction = lostFractionThroughTruncation(src, partCount(), 1925 truncatedBits); 1926 if (lost_fraction != lfExactlyZero 1927 && roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) { 1928 if (APInt::tcIncrement(parts, dstPartsCount)) 1929 return opInvalidOp; /* Overflow. */ 1930 } 1931 } else { 1932 lost_fraction = lfExactlyZero; 1933 } 1934 1935 /* Step 3: check if we fit in the destination. */ 1936 unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1; 1937 1938 if (sign) { 1939 if (!isSigned) { 1940 /* Negative numbers cannot be represented as unsigned. */ 1941 if (omsb != 0) 1942 return opInvalidOp; 1943 } else { 1944 /* It takes omsb bits to represent the unsigned integer value. 1945 We lose a bit for the sign, but care is needed as the 1946 maximally negative integer is a special case. */ 1947 if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb) 1948 return opInvalidOp; 1949 1950 /* This case can happen because of rounding. */ 1951 if (omsb > width) 1952 return opInvalidOp; 1953 } 1954 1955 APInt::tcNegate (parts, dstPartsCount); 1956 } else { 1957 if (omsb >= width + !isSigned) 1958 return opInvalidOp; 1959 } 1960 1961 if (lost_fraction == lfExactlyZero) { 1962 *isExact = true; 1963 return opOK; 1964 } else 1965 return opInexact; 1966} 1967 1968/* Same as convertToSignExtendedInteger, except we provide 1969 deterministic values in case of an invalid operation exception, 1970 namely zero for NaNs and the minimal or maximal value respectively 1971 for underflow or overflow. 1972 The *isExact output tells whether the result is exact, in the sense 1973 that converting it back to the original floating point type produces 1974 the original value. This is almost equivalent to result==opOK, 1975 except for negative zeroes. 1976*/ 1977APFloat::opStatus 1978APFloat::convertToInteger(integerPart *parts, unsigned int width, 1979 bool isSigned, 1980 roundingMode rounding_mode, bool *isExact) const 1981{ 1982 opStatus fs; 1983 1984 fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode, 1985 isExact); 1986 1987 if (fs == opInvalidOp) { 1988 unsigned int bits, dstPartsCount; 1989 1990 dstPartsCount = partCountForBits(width); 1991 1992 if (category == fcNaN) 1993 bits = 0; 1994 else if (sign) 1995 bits = isSigned; 1996 else 1997 bits = width - isSigned; 1998 1999 APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits); 2000 if (sign && isSigned) 2001 APInt::tcShiftLeft(parts, dstPartsCount, width - 1); 2002 } 2003 2004 return fs; 2005} 2006 2007/* Convert an unsigned integer SRC to a floating point number, 2008 rounding according to ROUNDING_MODE. The sign of the floating 2009 point number is not modified. */ 2010APFloat::opStatus 2011APFloat::convertFromUnsignedParts(const integerPart *src, 2012 unsigned int srcCount, 2013 roundingMode rounding_mode) 2014{ 2015 unsigned int omsb, precision, dstCount; 2016 integerPart *dst; 2017 lostFraction lost_fraction; 2018 2019 assertArithmeticOK(*semantics); 2020 category = fcNormal; 2021 omsb = APInt::tcMSB(src, srcCount) + 1; 2022 dst = significandParts(); 2023 dstCount = partCount(); 2024 precision = semantics->precision; 2025 2026 /* We want the most significant PRECISON bits of SRC. There may not 2027 be that many; extract what we can. */ 2028 if (precision <= omsb) { 2029 exponent = omsb - 1; 2030 lost_fraction = lostFractionThroughTruncation(src, srcCount, 2031 omsb - precision); 2032 APInt::tcExtract(dst, dstCount, src, precision, omsb - precision); 2033 } else { 2034 exponent = precision - 1; 2035 lost_fraction = lfExactlyZero; 2036 APInt::tcExtract(dst, dstCount, src, omsb, 0); 2037 } 2038 2039 return normalize(rounding_mode, lost_fraction); 2040} 2041 2042APFloat::opStatus 2043APFloat::convertFromAPInt(const APInt &Val, 2044 bool isSigned, 2045 roundingMode rounding_mode) 2046{ 2047 unsigned int partCount = Val.getNumWords(); 2048 APInt api = Val; 2049 2050 sign = false; 2051 if (isSigned && api.isNegative()) { 2052 sign = true; 2053 api = -api; 2054 } 2055 2056 return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode); 2057} 2058 2059/* Convert a two's complement integer SRC to a floating point number, 2060 rounding according to ROUNDING_MODE. ISSIGNED is true if the 2061 integer is signed, in which case it must be sign-extended. */ 2062APFloat::opStatus 2063APFloat::convertFromSignExtendedInteger(const integerPart *src, 2064 unsigned int srcCount, 2065 bool isSigned, 2066 roundingMode rounding_mode) 2067{ 2068 opStatus status; 2069 2070 assertArithmeticOK(*semantics); 2071 if (isSigned 2072 && APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) { 2073 integerPart *copy; 2074 2075 /* If we're signed and negative negate a copy. */ 2076 sign = true; 2077 copy = new integerPart[srcCount]; 2078 APInt::tcAssign(copy, src, srcCount); 2079 APInt::tcNegate(copy, srcCount); 2080 status = convertFromUnsignedParts(copy, srcCount, rounding_mode); 2081 delete [] copy; 2082 } else { 2083 sign = false; 2084 status = convertFromUnsignedParts(src, srcCount, rounding_mode); 2085 } 2086 2087 return status; 2088} 2089 2090/* FIXME: should this just take a const APInt reference? */ 2091APFloat::opStatus 2092APFloat::convertFromZeroExtendedInteger(const integerPart *parts, 2093 unsigned int width, bool isSigned, 2094 roundingMode rounding_mode) 2095{ 2096 unsigned int partCount = partCountForBits(width); 2097 APInt api = APInt(width, partCount, parts); 2098 2099 sign = false; 2100 if(isSigned && APInt::tcExtractBit(parts, width - 1)) { 2101 sign = true; 2102 api = -api; 2103 } 2104 2105 return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode); 2106} 2107 2108APFloat::opStatus 2109APFloat::convertFromHexadecimalString(const char *p, 2110 roundingMode rounding_mode) 2111{ 2112 lostFraction lost_fraction; 2113 integerPart *significand; 2114 unsigned int bitPos, partsCount; 2115 const char *dot, *firstSignificantDigit; 2116 2117 zeroSignificand(); 2118 exponent = 0; 2119 category = fcNormal; 2120 2121 significand = significandParts(); 2122 partsCount = partCount(); 2123 bitPos = partsCount * integerPartWidth; 2124 2125 /* Skip leading zeroes and any (hexa)decimal point. */ 2126 p = skipLeadingZeroesAndAnyDot(p, &dot); 2127 firstSignificantDigit = p; 2128 2129 for(;;) { 2130 integerPart hex_value; 2131 2132 if(*p == '.') { 2133 assert(dot == 0); 2134 dot = p++; 2135 } 2136 2137 hex_value = hexDigitValue(*p); 2138 if(hex_value == -1U) { 2139 lost_fraction = lfExactlyZero; 2140 break; 2141 } 2142 2143 p++; 2144 2145 /* Store the number whilst 4-bit nibbles remain. */ 2146 if(bitPos) { 2147 bitPos -= 4; 2148 hex_value <<= bitPos % integerPartWidth; 2149 significand[bitPos / integerPartWidth] |= hex_value; 2150 } else { 2151 lost_fraction = trailingHexadecimalFraction(p, hex_value); 2152 while(hexDigitValue(*p) != -1U) 2153 p++; 2154 break; 2155 } 2156 } 2157 2158 /* Hex floats require an exponent but not a hexadecimal point. */ 2159 assert(*p == 'p' || *p == 'P'); 2160 2161 /* Ignore the exponent if we are zero. */ 2162 if(p != firstSignificantDigit) { 2163 int expAdjustment; 2164 2165 /* Implicit hexadecimal point? */ 2166 if(!dot) 2167 dot = p; 2168 2169 /* Calculate the exponent adjustment implicit in the number of 2170 significant digits. */ 2171 expAdjustment = static_cast<int>(dot - firstSignificantDigit); 2172 if(expAdjustment < 0) 2173 expAdjustment++; 2174 expAdjustment = expAdjustment * 4 - 1; 2175 2176 /* Adjust for writing the significand starting at the most 2177 significant nibble. */ 2178 expAdjustment += semantics->precision; 2179 expAdjustment -= partsCount * integerPartWidth; 2180 2181 /* Adjust for the given exponent. */ 2182 exponent = totalExponent(p, expAdjustment); 2183 } 2184 2185 return normalize(rounding_mode, lost_fraction); 2186} 2187 2188APFloat::opStatus 2189APFloat::roundSignificandWithExponent(const integerPart *decSigParts, 2190 unsigned sigPartCount, int exp, 2191 roundingMode rounding_mode) 2192{ 2193 unsigned int parts, pow5PartCount; 2194 fltSemantics calcSemantics = { 32767, -32767, 0, true }; 2195 integerPart pow5Parts[maxPowerOfFiveParts]; 2196 bool isNearest; 2197 2198 isNearest = (rounding_mode == rmNearestTiesToEven 2199 || rounding_mode == rmNearestTiesToAway); 2200 2201 parts = partCountForBits(semantics->precision + 11); 2202 2203 /* Calculate pow(5, abs(exp)). */ 2204 pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp); 2205 2206 for (;; parts *= 2) { 2207 opStatus sigStatus, powStatus; 2208 unsigned int excessPrecision, truncatedBits; 2209 2210 calcSemantics.precision = parts * integerPartWidth - 1; 2211 excessPrecision = calcSemantics.precision - semantics->precision; 2212 truncatedBits = excessPrecision; 2213 2214 APFloat decSig(calcSemantics, fcZero, sign); 2215 APFloat pow5(calcSemantics, fcZero, false); 2216 2217 sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount, 2218 rmNearestTiesToEven); 2219 powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount, 2220 rmNearestTiesToEven); 2221 /* Add exp, as 10^n = 5^n * 2^n. */ 2222 decSig.exponent += exp; 2223 2224 lostFraction calcLostFraction; 2225 integerPart HUerr, HUdistance; 2226 unsigned int powHUerr; 2227 2228 if (exp >= 0) { 2229 /* multiplySignificand leaves the precision-th bit set to 1. */ 2230 calcLostFraction = decSig.multiplySignificand(pow5, NULL); 2231 powHUerr = powStatus != opOK; 2232 } else { 2233 calcLostFraction = decSig.divideSignificand(pow5); 2234 /* Denormal numbers have less precision. */ 2235 if (decSig.exponent < semantics->minExponent) { 2236 excessPrecision += (semantics->minExponent - decSig.exponent); 2237 truncatedBits = excessPrecision; 2238 if (excessPrecision > calcSemantics.precision) 2239 excessPrecision = calcSemantics.precision; 2240 } 2241 /* Extra half-ulp lost in reciprocal of exponent. */ 2242 powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2; 2243 } 2244 2245 /* Both multiplySignificand and divideSignificand return the 2246 result with the integer bit set. */ 2247 assert (APInt::tcExtractBit 2248 (decSig.significandParts(), calcSemantics.precision - 1) == 1); 2249 2250 HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK, 2251 powHUerr); 2252 HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(), 2253 excessPrecision, isNearest); 2254 2255 /* Are we guaranteed to round correctly if we truncate? */ 2256 if (HUdistance >= HUerr) { 2257 APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(), 2258 calcSemantics.precision - excessPrecision, 2259 excessPrecision); 2260 /* Take the exponent of decSig. If we tcExtract-ed less bits 2261 above we must adjust our exponent to compensate for the 2262 implicit right shift. */ 2263 exponent = (decSig.exponent + semantics->precision 2264 - (calcSemantics.precision - excessPrecision)); 2265 calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(), 2266 decSig.partCount(), 2267 truncatedBits); 2268 return normalize(rounding_mode, calcLostFraction); 2269 } 2270 } 2271} 2272 2273APFloat::opStatus 2274APFloat::convertFromDecimalString(const char *p, roundingMode rounding_mode) 2275{ 2276 decimalInfo D; 2277 opStatus fs; 2278 2279 /* Scan the text. */ 2280 interpretDecimal(p, &D); 2281 2282 /* Handle the quick cases. First the case of no significant digits, 2283 i.e. zero, and then exponents that are obviously too large or too 2284 small. Writing L for log 10 / log 2, a number d.ddddd*10^exp 2285 definitely overflows if 2286 2287 (exp - 1) * L >= maxExponent 2288 2289 and definitely underflows to zero where 2290 2291 (exp + 1) * L <= minExponent - precision 2292 2293 With integer arithmetic the tightest bounds for L are 2294 2295 93/28 < L < 196/59 [ numerator <= 256 ] 2296 42039/12655 < L < 28738/8651 [ numerator <= 65536 ] 2297 */ 2298 2299 if (decDigitValue(*D.firstSigDigit) >= 10U) { 2300 category = fcZero; 2301 fs = opOK; 2302 } else if ((D.normalizedExponent + 1) * 28738 2303 <= 8651 * (semantics->minExponent - (int) semantics->precision)) { 2304 /* Underflow to zero and round. */ 2305 zeroSignificand(); 2306 fs = normalize(rounding_mode, lfLessThanHalf); 2307 } else if ((D.normalizedExponent - 1) * 42039 2308 >= 12655 * semantics->maxExponent) { 2309 /* Overflow and round. */ 2310 fs = handleOverflow(rounding_mode); 2311 } else { 2312 integerPart *decSignificand; 2313 unsigned int partCount; 2314 2315 /* A tight upper bound on number of bits required to hold an 2316 N-digit decimal integer is N * 196 / 59. Allocate enough space 2317 to hold the full significand, and an extra part required by 2318 tcMultiplyPart. */ 2319 partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1; 2320 partCount = partCountForBits(1 + 196 * partCount / 59); 2321 decSignificand = new integerPart[partCount + 1]; 2322 partCount = 0; 2323 2324 /* Convert to binary efficiently - we do almost all multiplication 2325 in an integerPart. When this would overflow do we do a single 2326 bignum multiplication, and then revert again to multiplication 2327 in an integerPart. */ 2328 do { 2329 integerPart decValue, val, multiplier; 2330 2331 val = 0; 2332 multiplier = 1; 2333 2334 do { 2335 if (*p == '.') 2336 p++; 2337 2338 decValue = decDigitValue(*p++); 2339 multiplier *= 10; 2340 val = val * 10 + decValue; 2341 /* The maximum number that can be multiplied by ten with any 2342 digit added without overflowing an integerPart. */ 2343 } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10); 2344 2345 /* Multiply out the current part. */ 2346 APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val, 2347 partCount, partCount + 1, false); 2348 2349 /* If we used another part (likely but not guaranteed), increase 2350 the count. */ 2351 if (decSignificand[partCount]) 2352 partCount++; 2353 } while (p <= D.lastSigDigit); 2354 2355 category = fcNormal; 2356 fs = roundSignificandWithExponent(decSignificand, partCount, 2357 D.exponent, rounding_mode); 2358 2359 delete [] decSignificand; 2360 } 2361 2362 return fs; 2363} 2364 2365APFloat::opStatus 2366APFloat::convertFromString(const char *p, roundingMode rounding_mode) 2367{ 2368 assertArithmeticOK(*semantics); 2369 2370 /* Handle a leading minus sign. */ 2371 if(*p == '-') 2372 sign = 1, p++; 2373 else 2374 sign = 0; 2375 2376 if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X')) 2377 return convertFromHexadecimalString(p + 2, rounding_mode); 2378 2379 return convertFromDecimalString(p, rounding_mode); 2380} 2381 2382/* Write out a hexadecimal representation of the floating point value 2383 to DST, which must be of sufficient size, in the C99 form 2384 [-]0xh.hhhhp[+-]d. Return the number of characters written, 2385 excluding the terminating NUL. 2386 2387 If UPPERCASE, the output is in upper case, otherwise in lower case. 2388 2389 HEXDIGITS digits appear altogether, rounding the value if 2390 necessary. If HEXDIGITS is 0, the minimal precision to display the 2391 number precisely is used instead. If nothing would appear after 2392 the decimal point it is suppressed. 2393 2394 The decimal exponent is always printed and has at least one digit. 2395 Zero values display an exponent of zero. Infinities and NaNs 2396 appear as "infinity" or "nan" respectively. 2397 2398 The above rules are as specified by C99. There is ambiguity about 2399 what the leading hexadecimal digit should be. This implementation 2400 uses whatever is necessary so that the exponent is displayed as 2401 stored. This implies the exponent will fall within the IEEE format 2402 range, and the leading hexadecimal digit will be 0 (for denormals), 2403 1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with 2404 any other digits zero). 2405*/ 2406unsigned int 2407APFloat::convertToHexString(char *dst, unsigned int hexDigits, 2408 bool upperCase, roundingMode rounding_mode) const 2409{ 2410 char *p; 2411 2412 assertArithmeticOK(*semantics); 2413 2414 p = dst; 2415 if (sign) 2416 *dst++ = '-'; 2417 2418 switch (category) { 2419 case fcInfinity: 2420 memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1); 2421 dst += sizeof infinityL - 1; 2422 break; 2423 2424 case fcNaN: 2425 memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1); 2426 dst += sizeof NaNU - 1; 2427 break; 2428 2429 case fcZero: 2430 *dst++ = '0'; 2431 *dst++ = upperCase ? 'X': 'x'; 2432 *dst++ = '0'; 2433 if (hexDigits > 1) { 2434 *dst++ = '.'; 2435 memset (dst, '0', hexDigits - 1); 2436 dst += hexDigits - 1; 2437 } 2438 *dst++ = upperCase ? 'P': 'p'; 2439 *dst++ = '0'; 2440 break; 2441 2442 case fcNormal: 2443 dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode); 2444 break; 2445 } 2446 2447 *dst = 0; 2448 2449 return static_cast<unsigned int>(dst - p); 2450} 2451 2452/* Does the hard work of outputting the correctly rounded hexadecimal 2453 form of a normal floating point number with the specified number of 2454 hexadecimal digits. If HEXDIGITS is zero the minimum number of 2455 digits necessary to print the value precisely is output. */ 2456char * 2457APFloat::convertNormalToHexString(char *dst, unsigned int hexDigits, 2458 bool upperCase, 2459 roundingMode rounding_mode) const 2460{ 2461 unsigned int count, valueBits, shift, partsCount, outputDigits; 2462 const char *hexDigitChars; 2463 const integerPart *significand; 2464 char *p; 2465 bool roundUp; 2466 2467 *dst++ = '0'; 2468 *dst++ = upperCase ? 'X': 'x'; 2469 2470 roundUp = false; 2471 hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower; 2472 2473 significand = significandParts(); 2474 partsCount = partCount(); 2475 2476 /* +3 because the first digit only uses the single integer bit, so 2477 we have 3 virtual zero most-significant-bits. */ 2478 valueBits = semantics->precision + 3; 2479 shift = integerPartWidth - valueBits % integerPartWidth; 2480 2481 /* The natural number of digits required ignoring trailing 2482 insignificant zeroes. */ 2483 outputDigits = (valueBits - significandLSB () + 3) / 4; 2484 2485 /* hexDigits of zero means use the required number for the 2486 precision. Otherwise, see if we are truncating. If we are, 2487 find out if we need to round away from zero. */ 2488 if (hexDigits) { 2489 if (hexDigits < outputDigits) { 2490 /* We are dropping non-zero bits, so need to check how to round. 2491 "bits" is the number of dropped bits. */ 2492 unsigned int bits; 2493 lostFraction fraction; 2494 2495 bits = valueBits - hexDigits * 4; 2496 fraction = lostFractionThroughTruncation (significand, partsCount, bits); 2497 roundUp = roundAwayFromZero(rounding_mode, fraction, bits); 2498 } 2499 outputDigits = hexDigits; 2500 } 2501 2502 /* Write the digits consecutively, and start writing in the location 2503 of the hexadecimal point. We move the most significant digit 2504 left and add the hexadecimal point later. */ 2505 p = ++dst; 2506 2507 count = (valueBits + integerPartWidth - 1) / integerPartWidth; 2508 2509 while (outputDigits && count) { 2510 integerPart part; 2511 2512 /* Put the most significant integerPartWidth bits in "part". */ 2513 if (--count == partsCount) 2514 part = 0; /* An imaginary higher zero part. */ 2515 else 2516 part = significand[count] << shift; 2517 2518 if (count && shift) 2519 part |= significand[count - 1] >> (integerPartWidth - shift); 2520 2521 /* Convert as much of "part" to hexdigits as we can. */ 2522 unsigned int curDigits = integerPartWidth / 4; 2523 2524 if (curDigits > outputDigits) 2525 curDigits = outputDigits; 2526 dst += partAsHex (dst, part, curDigits, hexDigitChars); 2527 outputDigits -= curDigits; 2528 } 2529 2530 if (roundUp) { 2531 char *q = dst; 2532 2533 /* Note that hexDigitChars has a trailing '0'. */ 2534 do { 2535 q--; 2536 *q = hexDigitChars[hexDigitValue (*q) + 1]; 2537 } while (*q == '0'); 2538 assert (q >= p); 2539 } else { 2540 /* Add trailing zeroes. */ 2541 memset (dst, '0', outputDigits); 2542 dst += outputDigits; 2543 } 2544 2545 /* Move the most significant digit to before the point, and if there 2546 is something after the decimal point add it. This must come 2547 after rounding above. */ 2548 p[-1] = p[0]; 2549 if (dst -1 == p) 2550 dst--; 2551 else 2552 p[0] = '.'; 2553 2554 /* Finally output the exponent. */ 2555 *dst++ = upperCase ? 'P': 'p'; 2556 2557 return writeSignedDecimal (dst, exponent); 2558} 2559 2560// For good performance it is desirable for different APFloats 2561// to produce different integers. 2562uint32_t 2563APFloat::getHashValue() const 2564{ 2565 if (category==fcZero) return sign<<8 | semantics->precision ; 2566 else if (category==fcInfinity) return sign<<9 | semantics->precision; 2567 else if (category==fcNaN) return 1<<10 | semantics->precision; 2568 else { 2569 uint32_t hash = sign<<11 | semantics->precision | exponent<<12; 2570 const integerPart* p = significandParts(); 2571 for (int i=partCount(); i>0; i--, p++) 2572 hash ^= ((uint32_t)*p) ^ (uint32_t)((*p)>>32); 2573 return hash; 2574 } 2575} 2576 2577// Conversion from APFloat to/from host float/double. It may eventually be 2578// possible to eliminate these and have everybody deal with APFloats, but that 2579// will take a while. This approach will not easily extend to long double. 2580// Current implementation requires integerPartWidth==64, which is correct at 2581// the moment but could be made more general. 2582 2583// Denormals have exponent minExponent in APFloat, but minExponent-1 in 2584// the actual IEEE respresentations. We compensate for that here. 2585 2586APInt 2587APFloat::convertF80LongDoubleAPFloatToAPInt() const 2588{ 2589 assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended); 2590 assert (partCount()==2); 2591 2592 uint64_t myexponent, mysignificand; 2593 2594 if (category==fcNormal) { 2595 myexponent = exponent+16383; //bias 2596 mysignificand = significandParts()[0]; 2597 if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL)) 2598 myexponent = 0; // denormal 2599 } else if (category==fcZero) { 2600 myexponent = 0; 2601 mysignificand = 0; 2602 } else if (category==fcInfinity) { 2603 myexponent = 0x7fff; 2604 mysignificand = 0x8000000000000000ULL; 2605 } else { 2606 assert(category == fcNaN && "Unknown category"); 2607 myexponent = 0x7fff; 2608 mysignificand = significandParts()[0]; 2609 } 2610 2611 uint64_t words[2]; 2612 words[0] = mysignificand; 2613 words[1] = ((uint64_t)(sign & 1) << 15) | 2614 (myexponent & 0x7fffLL); 2615 return APInt(80, 2, words); 2616} 2617 2618APInt 2619APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const 2620{ 2621 assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble); 2622 assert (partCount()==2); 2623 2624 uint64_t myexponent, mysignificand, myexponent2, mysignificand2; 2625 2626 if (category==fcNormal) { 2627 myexponent = exponent + 1023; //bias 2628 myexponent2 = exponent2 + 1023; 2629 mysignificand = significandParts()[0]; 2630 mysignificand2 = significandParts()[1]; 2631 if (myexponent==1 && !(mysignificand & 0x10000000000000LL)) 2632 myexponent = 0; // denormal 2633 if (myexponent2==1 && !(mysignificand2 & 0x10000000000000LL)) 2634 myexponent2 = 0; // denormal 2635 } else if (category==fcZero) { 2636 myexponent = 0; 2637 mysignificand = 0; 2638 myexponent2 = 0; 2639 mysignificand2 = 0; 2640 } else if (category==fcInfinity) { 2641 myexponent = 0x7ff; 2642 myexponent2 = 0; 2643 mysignificand = 0; 2644 mysignificand2 = 0; 2645 } else { 2646 assert(category == fcNaN && "Unknown category"); 2647 myexponent = 0x7ff; 2648 mysignificand = significandParts()[0]; 2649 myexponent2 = exponent2; 2650 mysignificand2 = significandParts()[1]; 2651 } 2652 2653 uint64_t words[2]; 2654 words[0] = ((uint64_t)(sign & 1) << 63) | 2655 ((myexponent & 0x7ff) << 52) | 2656 (mysignificand & 0xfffffffffffffLL); 2657 words[1] = ((uint64_t)(sign2 & 1) << 63) | 2658 ((myexponent2 & 0x7ff) << 52) | 2659 (mysignificand2 & 0xfffffffffffffLL); 2660 return APInt(128, 2, words); 2661} 2662 2663APInt 2664APFloat::convertDoubleAPFloatToAPInt() const 2665{ 2666 assert(semantics == (const llvm::fltSemantics*)&IEEEdouble); 2667 assert (partCount()==1); 2668 2669 uint64_t myexponent, mysignificand; 2670 2671 if (category==fcNormal) { 2672 myexponent = exponent+1023; //bias 2673 mysignificand = *significandParts(); 2674 if (myexponent==1 && !(mysignificand & 0x10000000000000LL)) 2675 myexponent = 0; // denormal 2676 } else if (category==fcZero) { 2677 myexponent = 0; 2678 mysignificand = 0; 2679 } else if (category==fcInfinity) { 2680 myexponent = 0x7ff; 2681 mysignificand = 0; 2682 } else { 2683 assert(category == fcNaN && "Unknown category!"); 2684 myexponent = 0x7ff; 2685 mysignificand = *significandParts(); 2686 } 2687 2688 return APInt(64, ((((uint64_t)(sign & 1) << 63) | 2689 ((myexponent & 0x7ff) << 52) | 2690 (mysignificand & 0xfffffffffffffLL)))); 2691} 2692 2693APInt 2694APFloat::convertFloatAPFloatToAPInt() const 2695{ 2696 assert(semantics == (const llvm::fltSemantics*)&IEEEsingle); 2697 assert (partCount()==1); 2698 2699 uint32_t myexponent, mysignificand; 2700 2701 if (category==fcNormal) { 2702 myexponent = exponent+127; //bias 2703 mysignificand = (uint32_t)*significandParts(); 2704 if (myexponent == 1 && !(mysignificand & 0x800000)) 2705 myexponent = 0; // denormal 2706 } else if (category==fcZero) { 2707 myexponent = 0; 2708 mysignificand = 0; 2709 } else if (category==fcInfinity) { 2710 myexponent = 0xff; 2711 mysignificand = 0; 2712 } else { 2713 assert(category == fcNaN && "Unknown category!"); 2714 myexponent = 0xff; 2715 mysignificand = (uint32_t)*significandParts(); 2716 } 2717 2718 return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) | 2719 (mysignificand & 0x7fffff))); 2720} 2721 2722// This function creates an APInt that is just a bit map of the floating 2723// point constant as it would appear in memory. It is not a conversion, 2724// and treating the result as a normal integer is unlikely to be useful. 2725 2726APInt 2727APFloat::bitcastToAPInt() const 2728{ 2729 if (semantics == (const llvm::fltSemantics*)&IEEEsingle) 2730 return convertFloatAPFloatToAPInt(); 2731 2732 if (semantics == (const llvm::fltSemantics*)&IEEEdouble) 2733 return convertDoubleAPFloatToAPInt(); 2734 2735 if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble) 2736 return convertPPCDoubleDoubleAPFloatToAPInt(); 2737 2738 assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended && 2739 "unknown format!"); 2740 return convertF80LongDoubleAPFloatToAPInt(); 2741} 2742 2743float 2744APFloat::convertToFloat() const 2745{ 2746 assert(semantics == (const llvm::fltSemantics*)&IEEEsingle); 2747 APInt api = bitcastToAPInt(); 2748 return api.bitsToFloat(); 2749} 2750 2751double 2752APFloat::convertToDouble() const 2753{ 2754 assert(semantics == (const llvm::fltSemantics*)&IEEEdouble); 2755 APInt api = bitcastToAPInt(); 2756 return api.bitsToDouble(); 2757} 2758 2759/// Integer bit is explicit in this format. Intel hardware (387 and later) 2760/// does not support these bit patterns: 2761/// exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity") 2762/// exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN") 2763/// exponent = 0, integer bit 1 ("pseudodenormal") 2764/// exponent!=0 nor all 1's, integer bit 0 ("unnormal") 2765/// At the moment, the first two are treated as NaNs, the second two as Normal. 2766void 2767APFloat::initFromF80LongDoubleAPInt(const APInt &api) 2768{ 2769 assert(api.getBitWidth()==80); 2770 uint64_t i1 = api.getRawData()[0]; 2771 uint64_t i2 = api.getRawData()[1]; 2772 uint64_t myexponent = (i2 & 0x7fff); 2773 uint64_t mysignificand = i1; 2774 2775 initialize(&APFloat::x87DoubleExtended); 2776 assert(partCount()==2); 2777 2778 sign = static_cast<unsigned int>(i2>>15); 2779 if (myexponent==0 && mysignificand==0) { 2780 // exponent, significand meaningless 2781 category = fcZero; 2782 } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) { 2783 // exponent, significand meaningless 2784 category = fcInfinity; 2785 } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) { 2786 // exponent meaningless 2787 category = fcNaN; 2788 significandParts()[0] = mysignificand; 2789 significandParts()[1] = 0; 2790 } else { 2791 category = fcNormal; 2792 exponent = myexponent - 16383; 2793 significandParts()[0] = mysignificand; 2794 significandParts()[1] = 0; 2795 if (myexponent==0) // denormal 2796 exponent = -16382; 2797 } 2798} 2799 2800void 2801APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api) 2802{ 2803 assert(api.getBitWidth()==128); 2804 uint64_t i1 = api.getRawData()[0]; 2805 uint64_t i2 = api.getRawData()[1]; 2806 uint64_t myexponent = (i1 >> 52) & 0x7ff; 2807 uint64_t mysignificand = i1 & 0xfffffffffffffLL; 2808 uint64_t myexponent2 = (i2 >> 52) & 0x7ff; 2809 uint64_t mysignificand2 = i2 & 0xfffffffffffffLL; 2810 2811 initialize(&APFloat::PPCDoubleDouble); 2812 assert(partCount()==2); 2813 2814 sign = static_cast<unsigned int>(i1>>63); 2815 sign2 = static_cast<unsigned int>(i2>>63); 2816 if (myexponent==0 && mysignificand==0) { 2817 // exponent, significand meaningless 2818 // exponent2 and significand2 are required to be 0; we don't check 2819 category = fcZero; 2820 } else if (myexponent==0x7ff && mysignificand==0) { 2821 // exponent, significand meaningless 2822 // exponent2 and significand2 are required to be 0; we don't check 2823 category = fcInfinity; 2824 } else if (myexponent==0x7ff && mysignificand!=0) { 2825 // exponent meaningless. So is the whole second word, but keep it 2826 // for determinism. 2827 category = fcNaN; 2828 exponent2 = myexponent2; 2829 significandParts()[0] = mysignificand; 2830 significandParts()[1] = mysignificand2; 2831 } else { 2832 category = fcNormal; 2833 // Note there is no category2; the second word is treated as if it is 2834 // fcNormal, although it might be something else considered by itself. 2835 exponent = myexponent - 1023; 2836 exponent2 = myexponent2 - 1023; 2837 significandParts()[0] = mysignificand; 2838 significandParts()[1] = mysignificand2; 2839 if (myexponent==0) // denormal 2840 exponent = -1022; 2841 else 2842 significandParts()[0] |= 0x10000000000000LL; // integer bit 2843 if (myexponent2==0) 2844 exponent2 = -1022; 2845 else 2846 significandParts()[1] |= 0x10000000000000LL; // integer bit 2847 } 2848} 2849 2850void 2851APFloat::initFromDoubleAPInt(const APInt &api) 2852{ 2853 assert(api.getBitWidth()==64); 2854 uint64_t i = *api.getRawData(); 2855 uint64_t myexponent = (i >> 52) & 0x7ff; 2856 uint64_t mysignificand = i & 0xfffffffffffffLL; 2857 2858 initialize(&APFloat::IEEEdouble); 2859 assert(partCount()==1); 2860 2861 sign = static_cast<unsigned int>(i>>63); 2862 if (myexponent==0 && mysignificand==0) { 2863 // exponent, significand meaningless 2864 category = fcZero; 2865 } else if (myexponent==0x7ff && mysignificand==0) { 2866 // exponent, significand meaningless 2867 category = fcInfinity; 2868 } else if (myexponent==0x7ff && mysignificand!=0) { 2869 // exponent meaningless 2870 category = fcNaN; 2871 *significandParts() = mysignificand; 2872 } else { 2873 category = fcNormal; 2874 exponent = myexponent - 1023; 2875 *significandParts() = mysignificand; 2876 if (myexponent==0) // denormal 2877 exponent = -1022; 2878 else 2879 *significandParts() |= 0x10000000000000LL; // integer bit 2880 } 2881} 2882 2883void 2884APFloat::initFromFloatAPInt(const APInt & api) 2885{ 2886 assert(api.getBitWidth()==32); 2887 uint32_t i = (uint32_t)*api.getRawData(); 2888 uint32_t myexponent = (i >> 23) & 0xff; 2889 uint32_t mysignificand = i & 0x7fffff; 2890 2891 initialize(&APFloat::IEEEsingle); 2892 assert(partCount()==1); 2893 2894 sign = i >> 31; 2895 if (myexponent==0 && mysignificand==0) { 2896 // exponent, significand meaningless 2897 category = fcZero; 2898 } else if (myexponent==0xff && mysignificand==0) { 2899 // exponent, significand meaningless 2900 category = fcInfinity; 2901 } else if (myexponent==0xff && mysignificand!=0) { 2902 // sign, exponent, significand meaningless 2903 category = fcNaN; 2904 *significandParts() = mysignificand; 2905 } else { 2906 category = fcNormal; 2907 exponent = myexponent - 127; //bias 2908 *significandParts() = mysignificand; 2909 if (myexponent==0) // denormal 2910 exponent = -126; 2911 else 2912 *significandParts() |= 0x800000; // integer bit 2913 } 2914} 2915 2916/// Treat api as containing the bits of a floating point number. Currently 2917/// we infer the floating point type from the size of the APInt. The 2918/// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful 2919/// when the size is anything else). 2920void 2921APFloat::initFromAPInt(const APInt& api, bool isIEEE) 2922{ 2923 if (api.getBitWidth() == 32) 2924 return initFromFloatAPInt(api); 2925 else if (api.getBitWidth()==64) 2926 return initFromDoubleAPInt(api); 2927 else if (api.getBitWidth()==80) 2928 return initFromF80LongDoubleAPInt(api); 2929 else if (api.getBitWidth()==128 && !isIEEE) 2930 return initFromPPCDoubleDoubleAPInt(api); 2931 else 2932 assert(0); 2933} 2934 2935APFloat::APFloat(const APInt& api, bool isIEEE) 2936{ 2937 initFromAPInt(api, isIEEE); 2938} 2939 2940APFloat::APFloat(float f) 2941{ 2942 APInt api = APInt(32, 0); 2943 initFromAPInt(api.floatToBits(f)); 2944} 2945 2946APFloat::APFloat(double d) 2947{ 2948 APInt api = APInt(64, 0); 2949 initFromAPInt(api.doubleToBits(d)); 2950} 2951