X-Git-Url: http://demsky.eecs.uci.edu/git/?a=blobdiff_plain;f=lib%2FSupport%2FAPFloat.cpp;h=cc86e795e7927e6616b05b0441bd3a66f062a099;hb=e6be34a53ecbe8c2ff9f0793b13d847e94c0de91;hp=2037210e05ba91ee4e97d5189614ac2214437433;hpb=33d4c92e78a32a8e07ab3ebc1487b07304875ddd;p=oota-llvm.git diff --git a/lib/Support/APFloat.cpp b/lib/Support/APFloat.cpp index 2037210e05b..cc86e795e79 100644 --- a/lib/Support/APFloat.cpp +++ b/lib/Support/APFloat.cpp @@ -2,8 +2,8 @@ // // The LLVM Compiler Infrastructure // -// This file was developed by Neil Booth and is distributed under the -// University of Illinois Open Source License. See LICENSE.TXT for details. +// This file is distributed under the University of Illinois Open Source +// License. See LICENSE.TXT for details. // //===----------------------------------------------------------------------===// // @@ -12,9 +12,10 @@ // //===----------------------------------------------------------------------===// +#include "llvm/ADT/APFloat.h" +#include "llvm/ADT/FoldingSet.h" #include #include -#include "llvm/ADT/APFloat.h" #include "llvm/Support/MathExtras.h" using namespace llvm; @@ -39,17 +40,39 @@ namespace llvm { /* Number of bits in the significand. This includes the integer bit. */ - unsigned char precision; + unsigned int precision; - /* If the target format has an implicit integer bit. */ - bool implicitIntegerBit; + /* True if arithmetic is supported. */ + unsigned int arithmeticOK; }; const fltSemantics APFloat::IEEEsingle = { 127, -126, 24, true }; const fltSemantics APFloat::IEEEdouble = { 1023, -1022, 53, true }; const fltSemantics APFloat::IEEEquad = { 16383, -16382, 113, true }; - const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64, false }; - const fltSemantics APFloat::Bogus = { 0, 0, 0, false }; + const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64, true }; + const fltSemantics APFloat::Bogus = { 0, 0, 0, true }; + + // The PowerPC format consists of two doubles. It does not map cleanly + // onto the usual format above. For now only storage of constants of + // this type is supported, no arithmetic. + const fltSemantics APFloat::PPCDoubleDouble = { 1023, -1022, 106, false }; + + /* A tight upper bound on number of parts required to hold the value + pow(5, power) is + + power * 815 / (351 * integerPartWidth) + 1 + + However, whilst the result may require only this many parts, + because we are multiplying two values to get it, the + multiplication may require an extra part with the excess part + being zero (consider the trivial case of 1 * 1, tcFullMultiply + requires two parts to hold the single-part result). So we add an + extra one to guarantee enough space whilst multiplying. */ + const unsigned int maxExponent = 16383; + const unsigned int maxPrecision = 113; + const unsigned int maxPowerOfFiveExponent = maxExponent + maxPrecision - 1; + const unsigned int maxPowerOfFiveParts = 2 + ((maxPowerOfFiveExponent * 815) + / (351 * integerPartWidth)); } /* Put a bunch of private, handy routines in an anonymous namespace. */ @@ -61,20 +84,15 @@ namespace { return ((bits) + integerPartWidth - 1) / integerPartWidth; } - unsigned int - digitValue(unsigned int c) + /* Returns 0U-9U. Return values >= 10U are not digits. */ + inline unsigned int + decDigitValue(unsigned int c) { - unsigned int r; - - r = c - '0'; - if(r <= 9) - return r; - - return -1U; + return c - '0'; } unsigned int - hexDigitValue (unsigned int c) + hexDigitValue(unsigned int c) { unsigned int r; @@ -93,9 +111,56 @@ namespace { return -1U; } + inline void + assertArithmeticOK(const llvm::fltSemantics &semantics) { + assert(semantics.arithmeticOK + && "Compile-time arithmetic does not support these semantics"); + } + + /* Return the value of a decimal exponent of the form + [+-]ddddddd. + + If the exponent overflows, returns a large exponent with the + appropriate sign. */ + int + readExponent(const char *p) + { + bool isNegative; + unsigned int absExponent; + const unsigned int overlargeExponent = 24000; /* FIXME. */ + + isNegative = (*p == '-'); + if (*p == '-' || *p == '+') + p++; + + absExponent = decDigitValue(*p++); + assert (absExponent < 10U); + + for (;;) { + unsigned int value; + + value = decDigitValue(*p); + if (value >= 10U) + break; + + p++; + value += absExponent * 10; + if (absExponent >= overlargeExponent) { + absExponent = overlargeExponent; + break; + } + absExponent = value; + } + + if (isNegative) + return -(int) absExponent; + else + return (int) absExponent; + } + /* This is ugly and needs cleaning up, but I don't immediately see how whilst remaining safe. */ - static int + int totalExponent(const char *p, int exponentAdjustment) { integerPart unsignedExponent; @@ -113,8 +178,8 @@ namespace { for(;;) { unsigned int value; - value = digitValue(*p); - if(value == -1U) + value = decDigitValue(*p); + if(value >= 10U) break; p++; @@ -157,6 +222,72 @@ namespace { return p; } + /* Given a normal decimal floating point number of the form + + dddd.dddd[eE][+-]ddd + + where the decimal point and exponent are optional, fill out the + structure D. Exponent is appropriate if the significand is + treated as an integer, and normalizedExponent if the significand + is taken to have the decimal point after a single leading + non-zero digit. + + If the value is zero, V->firstSigDigit points to a non-digit, and + the return exponent is zero. + */ + struct decimalInfo { + const char *firstSigDigit; + const char *lastSigDigit; + int exponent; + int normalizedExponent; + }; + + void + interpretDecimal(const char *p, decimalInfo *D) + { + const char *dot; + + p = skipLeadingZeroesAndAnyDot (p, &dot); + + D->firstSigDigit = p; + D->exponent = 0; + D->normalizedExponent = 0; + + for (;;) { + if (*p == '.') { + assert(dot == 0); + dot = p++; + } + if (decDigitValue(*p) >= 10U) + break; + p++; + } + + /* If number is all zerooes accept any exponent. */ + if (p != D->firstSigDigit) { + if (*p == 'e' || *p == 'E') + D->exponent = readExponent(p + 1); + + /* Implied decimal point? */ + if (!dot) + dot = p; + + /* Drop insignificant trailing zeroes. */ + do + do + p--; + while (*p == '0'); + while (*p == '.'); + + /* Adjust the exponents for any decimal point. */ + D->exponent += (dot - p) - (dot > p); + D->normalizedExponent = (D->exponent + (p - D->firstSigDigit) + - (dot > D->firstSigDigit && dot < p)); + } + + D->lastSigDigit = p; + } + /* Return the trailing fraction of a hexadecimal number. DIGITVALUE is the first hex digit of the fraction, P points to the next digit. */ @@ -237,6 +368,131 @@ namespace { return moreSignificant; } + /* The error from the true value, in half-ulps, on multiplying two + floating point numbers, which differ from the value they + approximate by at most HUE1 and HUE2 half-ulps, is strictly less + than the returned value. + + See "How to Read Floating Point Numbers Accurately" by William D + Clinger. */ + unsigned int + HUerrBound(bool inexactMultiply, unsigned int HUerr1, unsigned int HUerr2) + { + assert(HUerr1 < 2 || HUerr2 < 2 || (HUerr1 + HUerr2 < 8)); + + if (HUerr1 + HUerr2 == 0) + return inexactMultiply * 2; /* <= inexactMultiply half-ulps. */ + else + return inexactMultiply + 2 * (HUerr1 + HUerr2); + } + + /* The number of ulps from the boundary (zero, or half if ISNEAREST) + when the least significant BITS are truncated. BITS cannot be + zero. */ + integerPart + ulpsFromBoundary(const integerPart *parts, unsigned int bits, bool isNearest) + { + unsigned int count, partBits; + integerPart part, boundary; + + assert (bits != 0); + + bits--; + count = bits / integerPartWidth; + partBits = bits % integerPartWidth + 1; + + part = parts[count] & (~(integerPart) 0 >> (integerPartWidth - partBits)); + + if (isNearest) + boundary = (integerPart) 1 << (partBits - 1); + else + boundary = 0; + + if (count == 0) { + if (part - boundary <= boundary - part) + return part - boundary; + else + return boundary - part; + } + + if (part == boundary) { + while (--count) + if (parts[count]) + return ~(integerPart) 0; /* A lot. */ + + return parts[0]; + } else if (part == boundary - 1) { + while (--count) + if (~parts[count]) + return ~(integerPart) 0; /* A lot. */ + + return -parts[0]; + } + + return ~(integerPart) 0; /* A lot. */ + } + + /* Place pow(5, power) in DST, and return the number of parts used. + DST must be at least one part larger than size of the answer. */ + unsigned int + powerOf5(integerPart *dst, unsigned int power) + { + static integerPart firstEightPowers[] = { 1, 5, 25, 125, 625, 3125, + 15625, 78125 }; + static integerPart pow5s[maxPowerOfFiveParts * 2 + 5] = { 78125 * 5 }; + static unsigned int partsCount[16] = { 1 }; + + integerPart scratch[maxPowerOfFiveParts], *p1, *p2, *pow5; + unsigned int result; + + assert(power <= maxExponent); + + p1 = dst; + p2 = scratch; + + *p1 = firstEightPowers[power & 7]; + power >>= 3; + + result = 1; + pow5 = pow5s; + + for (unsigned int n = 0; power; power >>= 1, n++) { + unsigned int pc; + + pc = partsCount[n]; + + /* Calculate pow(5,pow(2,n+3)) if we haven't yet. */ + if (pc == 0) { + pc = partsCount[n - 1]; + APInt::tcFullMultiply(pow5, pow5 - pc, pow5 - pc, pc, pc); + pc *= 2; + if (pow5[pc - 1] == 0) + pc--; + partsCount[n] = pc; + } + + if (power & 1) { + integerPart *tmp; + + APInt::tcFullMultiply(p2, p1, pow5, result, pc); + result += pc; + if (p2[result - 1] == 0) + result--; + + /* Now result is in p1 with partsCount parts and p2 is scratch + space. */ + tmp = p1, p1 = p2, p2 = tmp; + } + + pow5 += pc; + } + + if (p1 != dst) + APInt::tcAssign(dst, p1, result); + + return result; + } + /* Zero at the end to avoid modular arithmetic when adding one; used when rounding up during hexadecimal output. */ static const char hexDigitsLower[] = "0123456789abcdef0"; @@ -249,7 +505,7 @@ namespace { /* Write out an integerPart in hexadecimal, starting with the most significant nibble. Write out exactly COUNT hexdigits, return COUNT. */ - static unsigned int + unsigned int partAsHex (char *dst, integerPart part, unsigned int count, const char *hexDigitChars) { @@ -267,7 +523,7 @@ namespace { } /* Write out an unsigned decimal integer. */ - static char * + char * writeUnsignedDecimal (char *dst, unsigned int n) { char buff[40], *p; @@ -285,7 +541,7 @@ namespace { } /* Write out a signed decimal integer. */ - static char * + char * writeSignedDecimal (char *dst, int value) { if (value < 0) { @@ -325,6 +581,8 @@ APFloat::assign(const APFloat &rhs) sign = rhs.sign; category = rhs.category; exponent = rhs.exponent; + sign2 = rhs.sign2; + exponent2 = rhs.exponent2; if(category == fcNormal || category == fcNaN) copySignificand(rhs); } @@ -339,6 +597,15 @@ APFloat::copySignificand(const APFloat &rhs) partCount()); } +/* Make this number a NaN, with an arbitrary but deterministic value + for the significand. */ +void +APFloat::makeNaN(void) +{ + category = fcNaN; + APInt::tcSet(significandParts(), ~0U, partCount()); +} + APFloat & APFloat::operator=(const APFloat &rhs) { @@ -361,10 +628,16 @@ APFloat::bitwiseIsEqual(const APFloat &rhs) const { category != rhs.category || sign != rhs.sign) return false; + if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble && + sign2 != rhs.sign2) + return false; if (category==fcZero || category==fcInfinity) return true; else if (category==fcNormal && exponent!=rhs.exponent) return false; + else if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble && + exponent2!=rhs.exponent2) + return false; else { int i= partCount(); const integerPart* p=significandParts(); @@ -379,6 +652,7 @@ APFloat::bitwiseIsEqual(const APFloat &rhs) const { APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value) { + assertArithmeticOK(ourSemantics); initialize(&ourSemantics); sign = 0; zeroSignificand(); @@ -390,15 +664,19 @@ APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value) APFloat::APFloat(const fltSemantics &ourSemantics, fltCategory ourCategory, bool negative) { + assertArithmeticOK(ourSemantics); initialize(&ourSemantics); category = ourCategory; sign = negative; if(category == fcNormal) category = fcZero; + else if (ourCategory == fcNaN) + makeNaN(); } APFloat::APFloat(const fltSemantics &ourSemantics, const char *text) { + assertArithmeticOK(ourSemantics); initialize(&ourSemantics); convertFromString(text, rmNearestTiesToEven); } @@ -414,6 +692,11 @@ APFloat::~APFloat() freeSignificand(); } +// Profile - This method 'profiles' an APFloat for use with FoldingSet. +void APFloat::Profile(FoldingSetNodeID& ID) const { + ID.Add(convertToAPInt()); +} + unsigned int APFloat::partCount() const { @@ -634,6 +917,9 @@ APFloat::divideSignificand(const APFloat &rhs) APInt::tcShiftLeft(dividend, partsCount, bit); } + /* Ensure the dividend >= divisor initially for the loop below. + Incidentally, this means that the division loop below is + guaranteed to set the integer bit to one. */ if(APInt::tcCompare(dividend, divisor, partsCount) < 0) { exponent--; APInt::tcShiftLeft(dividend, partsCount, 1); @@ -849,7 +1135,7 @@ APFloat::normalize(roundingMode rounding_mode, /* Keep OMSB up-to-date. */ if(omsb > (unsigned) exponentChange) - omsb -= (unsigned) exponentChange; + omsb -= exponentChange; else omsb = 0; } @@ -900,7 +1186,6 @@ APFloat::normalize(roundingMode rounding_mode, /* We have a non-zero denormal. */ assert(omsb < semantics->precision); - assert(exponent == semantics->minExponent); /* Canonicalize zeroes. */ if(omsb == 0) @@ -951,10 +1236,8 @@ APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract) case convolve(fcInfinity, fcInfinity): /* Differently signed infinities can only be validly subtracted. */ - if(sign ^ rhs.sign != subtract) { - category = fcNaN; - // Arbitrary but deterministic value for significand - APInt::tcSet(significandParts(), ~0U, partCount()); + if((sign ^ rhs.sign) != subtract) { + makeNaN(); return opInvalidOp; } @@ -975,7 +1258,7 @@ APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract) /* Determine if the operation on the absolute values is effectively an addition or subtraction. */ - subtract ^= (sign ^ rhs.sign); + subtract ^= (sign ^ rhs.sign) ? true : false; /* Are we bigger exponent-wise than the RHS? */ bits = exponent - rhs.exponent; @@ -1070,9 +1353,7 @@ APFloat::multiplySpecials(const APFloat &rhs) case convolve(fcZero, fcInfinity): case convolve(fcInfinity, fcZero): - category = fcNaN; - // Arbitrary but deterministic value for significand - APInt::tcSet(significandParts(), ~0U, partCount()); + makeNaN(); return opInvalidOp; case convolve(fcNormal, fcNormal): @@ -1114,9 +1395,7 @@ APFloat::divideSpecials(const APFloat &rhs) case convolve(fcInfinity, fcInfinity): case convolve(fcZero, fcZero): - category = fcNaN; - // Arbitrary but deterministic value for significand - APInt::tcSet(significandParts(), ~0U, partCount()); + makeNaN(); return opInvalidOp; case convolve(fcNormal, fcNormal): @@ -1153,6 +1432,8 @@ APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode, { opStatus fs; + assertArithmeticOK(*semantics); + fs = addOrSubtractSpecials(rhs, subtract); /* This return code means it was not a simple case. */ @@ -1197,6 +1478,7 @@ APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode) { opStatus fs; + assertArithmeticOK(*semantics); sign ^= rhs.sign; fs = multiplySpecials(rhs); @@ -1216,6 +1498,7 @@ APFloat::divide(const APFloat &rhs, roundingMode rounding_mode) { opStatus fs; + assertArithmeticOK(*semantics); sign ^= rhs.sign; fs = divideSpecials(rhs); @@ -1236,6 +1519,8 @@ APFloat::mod(const APFloat &rhs, roundingMode rounding_mode) opStatus fs; APFloat V = *this; unsigned int origSign = sign; + + assertArithmeticOK(*semantics); fs = V.divide(rhs, rmNearestTiesToEven); if (fs == opDivByZero) return fs; @@ -1247,8 +1532,8 @@ APFloat::mod(const APFloat &rhs, roundingMode rounding_mode) if (fs==opInvalidOp) return fs; - fs = V.convertFromInteger(x, parts * integerPartWidth, true, - rmNearestTiesToEven); + fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true, + rmNearestTiesToEven); assert(fs==opOK); // should always work fs = V.multiply(rhs, rounding_mode); @@ -1271,6 +1556,8 @@ APFloat::fusedMultiplyAdd(const APFloat &multiplicand, { opStatus fs; + assertArithmeticOK(*semantics); + /* Post-multiplication sign, before addition. */ sign ^= multiplicand.sign; @@ -1314,6 +1601,7 @@ APFloat::compare(const APFloat &rhs) const { cmpResult result; + assertArithmeticOK(*semantics); assert(semantics == rhs.semantics); switch(convolve(category, rhs.category)) { @@ -1389,6 +1677,7 @@ APFloat::convert(const fltSemantics &toSemantics, unsigned int newPartCount, oldPartCount; opStatus fs; + assertArithmeticOK(*semantics); lostFraction = lfExactlyZero; newPartCount = partCountForBits(toSemantics.precision + 1); oldPartCount = partCount(); @@ -1429,6 +1718,8 @@ APFloat::convert(const fltSemantics &toSemantics, fs = normalize(rounding_mode, lostFraction); } else if (category == fcNaN) { int shift = toSemantics.precision - semantics->precision; + // Do this now so significandParts gets the right answer + semantics = &toSemantics; // No normalization here, just truncate if (shift>0) APInt::tcShiftLeft(significandParts(), newPartCount, shift); @@ -1438,7 +1729,6 @@ APFloat::convert(const fltSemantics &toSemantics, // does not give you back the same bits. This is dubious, and we // don't currently do it. You're really supposed to get // an invalid operation signal at runtime, but nobody does that. - semantics = &toSemantics; fs = opOK; } else { semantics = &toSemantics; @@ -1450,7 +1740,8 @@ APFloat::convert(const fltSemantics &toSemantics, /* Convert a floating point number to an integer according to the rounding mode. If the rounded integer value is out of range this - returns an invalid operation exception. If the rounded value is in + returns an invalid operation exception and the contents of the + destination parts are unspecified. If the rounded value is in range but the floating point number is not the exact integer, the C standard doesn't require an inexact exception to be raised. IEEE 854 does require it so we do that. @@ -1458,132 +1749,209 @@ APFloat::convert(const fltSemantics &toSemantics, Note that for conversions to integer type the C standard requires round-to-zero to always be used. */ APFloat::opStatus -APFloat::convertToInteger(integerPart *parts, unsigned int width, - bool isSigned, - roundingMode rounding_mode) const +APFloat::convertToSignExtendedInteger(integerPart *parts, unsigned int width, + bool isSigned, + roundingMode rounding_mode) const { lostFraction lost_fraction; - unsigned int msb, partsCount; - int bits; + const integerPart *src; + unsigned int dstPartsCount, truncatedBits; - partsCount = partCountForBits(width); + assertArithmeticOK(*semantics); - /* Handle the three special cases first. We produce - a deterministic result even for the Invalid cases. */ - if (category == fcNaN) { - // Neither sign nor isSigned affects this. - APInt::tcSet(parts, 0, partsCount); - return opInvalidOp; - } - if (category == fcInfinity) { - if (!sign && isSigned) - APInt::tcSetLeastSignificantBits(parts, partsCount, width-1); - else if (!sign && !isSigned) - APInt::tcSetLeastSignificantBits(parts, partsCount, width); - else if (sign && isSigned) { - APInt::tcSetLeastSignificantBits(parts, partsCount, 1); - APInt::tcShiftLeft(parts, partsCount, width-1); - } else // sign && !isSigned - APInt::tcSet(parts, 0, partsCount); + /* Handle the three special cases first. */ + if(category == fcInfinity || category == fcNaN) return opInvalidOp; - } - if (category == fcZero) { - APInt::tcSet(parts, 0, partsCount); + + dstPartsCount = partCountForBits(width); + + if(category == fcZero) { + APInt::tcSet(parts, 0, dstPartsCount); return opOK; } - /* Shift the bit pattern so the fraction is lost. */ - APFloat tmp(*this); + src = significandParts(); - bits = (int) semantics->precision - 1 - exponent; - - if(bits > 0) { - lost_fraction = tmp.shiftSignificandRight(bits); + /* Step 1: place our absolute value, with any fraction truncated, in + the destination. */ + if (exponent < 0) { + /* Our absolute value is less than one; truncate everything. */ + APInt::tcSet(parts, 0, dstPartsCount); + truncatedBits = semantics->precision; } else { - if (-bits >= semantics->precision) { - // Unrepresentably large. - if (!sign && isSigned) - APInt::tcSetLeastSignificantBits(parts, partsCount, width-1); - else if (!sign && !isSigned) - APInt::tcSetLeastSignificantBits(parts, partsCount, width); - else if (sign && isSigned) { - APInt::tcSetLeastSignificantBits(parts, partsCount, 1); - APInt::tcShiftLeft(parts, partsCount, width-1); - } else // sign && !isSigned - APInt::tcSet(parts, 0, partsCount); - return (opStatus)(opOverflow | opInexact); + /* We want the most significant (exponent + 1) bits; the rest are + truncated. */ + unsigned int bits = exponent + 1U; + + /* Hopelessly large in magnitude? */ + if (bits > width) + return opInvalidOp; + + if (bits < semantics->precision) { + /* We truncate (semantics->precision - bits) bits. */ + truncatedBits = semantics->precision - bits; + APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits); + } else { + /* We want at least as many bits as are available. */ + APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0); + APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision); + truncatedBits = 0; } - tmp.shiftSignificandLeft(-bits); + } + + /* Step 2: work out any lost fraction, and increment the absolute + value if we would round away from zero. */ + if (truncatedBits) { + lost_fraction = lostFractionThroughTruncation(src, partCount(), + truncatedBits); + if (lost_fraction != lfExactlyZero + && roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) { + if (APInt::tcIncrement(parts, dstPartsCount)) + return opInvalidOp; /* Overflow. */ + } + } else { lost_fraction = lfExactlyZero; } - if(lost_fraction != lfExactlyZero - && tmp.roundAwayFromZero(rounding_mode, lost_fraction, 0)) - tmp.incrementSignificand(); + /* Step 3: check if we fit in the destination. */ + unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1; - msb = tmp.significandMSB(); + if (sign) { + if (!isSigned) { + /* Negative numbers cannot be represented as unsigned. */ + if (omsb != 0) + return opInvalidOp; + } else { + /* It takes omsb bits to represent the unsigned integer value. + We lose a bit for the sign, but care is needed as the + maximally negative integer is a special case. */ + if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb) + return opInvalidOp; + + /* This case can happen because of rounding. */ + if (omsb > width) + return opInvalidOp; + } - /* Negative numbers cannot be represented as unsigned. */ - if(!isSigned && tmp.sign && msb != -1U) - return opInvalidOp; + APInt::tcNegate (parts, dstPartsCount); + } else { + if (omsb >= width + !isSigned) + return opInvalidOp; + } - /* It takes exponent + 1 bits to represent the truncated floating - point number without its sign. We lose a bit for the sign, but - the maximally negative integer is a special case. */ - if(msb + 1 > width) /* !! Not same as msb >= width !! */ - return opInvalidOp; + if (lost_fraction == lfExactlyZero) + return opOK; + else + return opInexact; +} - if(isSigned && msb + 1 == width - && (!tmp.sign || tmp.significandLSB() != msb)) - return opInvalidOp; +/* Same as convertToSignExtendedInteger, except we provide + deterministic values in case of an invalid operation exception, + namely zero for NaNs and the minimal or maximal value respectively + for underflow or overflow. */ +APFloat::opStatus +APFloat::convertToInteger(integerPart *parts, unsigned int width, + bool isSigned, + roundingMode rounding_mode) const +{ + opStatus fs; - APInt::tcAssign(parts, tmp.significandParts(), partsCount); + fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode); - if(tmp.sign) - APInt::tcNegate(parts, partsCount); + if (fs == opInvalidOp) { + unsigned int bits, dstPartsCount; - if(lost_fraction == lfExactlyZero) - return opOK; - else - return opInexact; + dstPartsCount = partCountForBits(width); + + if (category == fcNaN) + bits = 0; + else if (sign) + bits = isSigned; + else + bits = width - isSigned; + + APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits); + if (sign && isSigned) + APInt::tcShiftLeft(parts, dstPartsCount, width - 1); + } + + return fs; } +/* Convert an unsigned integer SRC to a floating point number, + rounding according to ROUNDING_MODE. The sign of the floating + point number is not modified. */ APFloat::opStatus -APFloat::convertFromUnsignedInteger(integerPart *parts, - unsigned int partCount, - roundingMode rounding_mode) +APFloat::convertFromUnsignedParts(const integerPart *src, + unsigned int srcCount, + roundingMode rounding_mode) { - unsigned int msb, precision; + unsigned int omsb, precision, dstCount; + integerPart *dst; lostFraction lost_fraction; - msb = APInt::tcMSB(parts, partCount) + 1; - precision = semantics->precision; - + assertArithmeticOK(*semantics); category = fcNormal; - exponent = precision - 1; + omsb = APInt::tcMSB(src, srcCount) + 1; + dst = significandParts(); + dstCount = partCount(); + precision = semantics->precision; - if(msb > precision) { - exponent += (msb - precision); - lost_fraction = shiftRight(parts, partCount, msb - precision); - msb = precision; - } else + /* We want the most significant PRECISON bits of SRC. There may not + be that many; extract what we can. */ + if (precision <= omsb) { + exponent = omsb - 1; + lost_fraction = lostFractionThroughTruncation(src, srcCount, + omsb - precision); + APInt::tcExtract(dst, dstCount, src, precision, omsb - precision); + } else { + exponent = precision - 1; lost_fraction = lfExactlyZero; - - /* Copy the bit image. */ - zeroSignificand(); - APInt::tcAssign(significandParts(), parts, partCountForBits(msb)); + APInt::tcExtract(dst, dstCount, src, omsb, 0); + } return normalize(rounding_mode, lost_fraction); } +/* Convert a two's complement integer SRC to a floating point number, + rounding according to ROUNDING_MODE. ISSIGNED is true if the + integer is signed, in which case it must be sign-extended. */ APFloat::opStatus -APFloat::convertFromInteger(const integerPart *parts, unsigned int width, - bool isSigned, roundingMode rounding_mode) +APFloat::convertFromSignExtendedInteger(const integerPart *src, + unsigned int srcCount, + bool isSigned, + roundingMode rounding_mode) { - unsigned int partCount = partCountForBits(width); opStatus status; + + assertArithmeticOK(*semantics); + if (isSigned + && APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) { + integerPart *copy; + + /* If we're signed and negative negate a copy. */ + sign = true; + copy = new integerPart[srcCount]; + APInt::tcAssign(copy, src, srcCount); + APInt::tcNegate(copy, srcCount); + status = convertFromUnsignedParts(copy, srcCount, rounding_mode); + delete [] copy; + } else { + sign = false; + status = convertFromUnsignedParts(src, srcCount, rounding_mode); + } + + return status; +} + +/* FIXME: should this just take a const APInt reference? */ +APFloat::opStatus +APFloat::convertFromZeroExtendedInteger(const integerPart *parts, + unsigned int width, bool isSigned, + roundingMode rounding_mode) +{ + unsigned int partCount = partCountForBits(width); APInt api = APInt(width, partCount, parts); - integerPart *copy = new integerPart[partCount]; sign = false; if(isSigned && APInt::tcExtractBit(parts, width - 1)) { @@ -1591,9 +1959,7 @@ APFloat::convertFromInteger(const integerPart *parts, unsigned int width, api = -api; } - APInt::tcAssign(copy, api.getRawData(), partCount); - status = convertFromUnsignedInteger(copy, partCount, rounding_mode); - return status; + return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode); } APFloat::opStatus @@ -1676,9 +2042,187 @@ APFloat::convertFromHexadecimalString(const char *p, return normalize(rounding_mode, lost_fraction); } +APFloat::opStatus +APFloat::roundSignificandWithExponent(const integerPart *decSigParts, + unsigned sigPartCount, int exp, + roundingMode rounding_mode) +{ + unsigned int parts, pow5PartCount; + fltSemantics calcSemantics = { 32767, -32767, 0, true }; + integerPart pow5Parts[maxPowerOfFiveParts]; + bool isNearest; + + isNearest = (rounding_mode == rmNearestTiesToEven + || rounding_mode == rmNearestTiesToAway); + + parts = partCountForBits(semantics->precision + 11); + + /* Calculate pow(5, abs(exp)). */ + pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp); + + for (;; parts *= 2) { + opStatus sigStatus, powStatus; + unsigned int excessPrecision, truncatedBits; + + calcSemantics.precision = parts * integerPartWidth - 1; + excessPrecision = calcSemantics.precision - semantics->precision; + truncatedBits = excessPrecision; + + APFloat decSig(calcSemantics, fcZero, sign); + APFloat pow5(calcSemantics, fcZero, false); + + sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount, + rmNearestTiesToEven); + powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount, + rmNearestTiesToEven); + /* Add exp, as 10^n = 5^n * 2^n. */ + decSig.exponent += exp; + + lostFraction calcLostFraction; + integerPart HUerr, HUdistance, powHUerr; + + if (exp >= 0) { + /* multiplySignificand leaves the precision-th bit set to 1. */ + calcLostFraction = decSig.multiplySignificand(pow5, NULL); + powHUerr = powStatus != opOK; + } else { + calcLostFraction = decSig.divideSignificand(pow5); + /* Denormal numbers have less precision. */ + if (decSig.exponent < semantics->minExponent) { + excessPrecision += (semantics->minExponent - decSig.exponent); + truncatedBits = excessPrecision; + if (excessPrecision > calcSemantics.precision) + excessPrecision = calcSemantics.precision; + } + /* Extra half-ulp lost in reciprocal of exponent. */ + powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0: 2; + } + + /* Both multiplySignificand and divideSignificand return the + result with the integer bit set. */ + assert (APInt::tcExtractBit + (decSig.significandParts(), calcSemantics.precision - 1) == 1); + + HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK, + powHUerr); + HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(), + excessPrecision, isNearest); + + /* Are we guaranteed to round correctly if we truncate? */ + if (HUdistance >= HUerr) { + APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(), + calcSemantics.precision - excessPrecision, + excessPrecision); + /* Take the exponent of decSig. If we tcExtract-ed less bits + above we must adjust our exponent to compensate for the + implicit right shift. */ + exponent = (decSig.exponent + semantics->precision + - (calcSemantics.precision - excessPrecision)); + calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(), + decSig.partCount(), + truncatedBits); + return normalize(rounding_mode, calcLostFraction); + } + } +} + +APFloat::opStatus +APFloat::convertFromDecimalString(const char *p, roundingMode rounding_mode) +{ + decimalInfo D; + opStatus fs; + + /* Scan the text. */ + interpretDecimal(p, &D); + + /* Handle the quick cases. First the case of no significant digits, + i.e. zero, and then exponents that are obviously too large or too + small. Writing L for log 10 / log 2, a number d.ddddd*10^exp + definitely overflows if + + (exp - 1) * L >= maxExponent + + and definitely underflows to zero where + + (exp + 1) * L <= minExponent - precision + + With integer arithmetic the tightest bounds for L are + + 93/28 < L < 196/59 [ numerator <= 256 ] + 42039/12655 < L < 28738/8651 [ numerator <= 65536 ] + */ + + if (decDigitValue(*D.firstSigDigit) >= 10U) { + category = fcZero; + fs = opOK; + } else if ((D.normalizedExponent + 1) * 28738 + <= 8651 * (semantics->minExponent - (int) semantics->precision)) { + /* Underflow to zero and round. */ + zeroSignificand(); + fs = normalize(rounding_mode, lfLessThanHalf); + } else if ((D.normalizedExponent - 1) * 42039 + >= 12655 * semantics->maxExponent) { + /* Overflow and round. */ + fs = handleOverflow(rounding_mode); + } else { + integerPart *decSignificand; + unsigned int partCount; + + /* A tight upper bound on number of bits required to hold an + N-digit decimal integer is N * 196 / 59. Allocate enough space + to hold the full significand, and an extra part required by + tcMultiplyPart. */ + partCount = (D.lastSigDigit - D.firstSigDigit) + 1; + partCount = partCountForBits(1 + 196 * partCount / 59); + decSignificand = new integerPart[partCount + 1]; + partCount = 0; + + /* Convert to binary efficiently - we do almost all multiplication + in an integerPart. When this would overflow do we do a single + bignum multiplication, and then revert again to multiplication + in an integerPart. */ + do { + integerPart decValue, val, multiplier; + + val = 0; + multiplier = 1; + + do { + if (*p == '.') + p++; + + decValue = decDigitValue(*p++); + multiplier *= 10; + val = val * 10 + decValue; + /* The maximum number that can be multiplied by ten with any + digit added without overflowing an integerPart. */ + } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10); + + /* Multiply out the current part. */ + APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val, + partCount, partCount + 1, false); + + /* If we used another part (likely but not guaranteed), increase + the count. */ + if (decSignificand[partCount]) + partCount++; + } while (p <= D.lastSigDigit); + + category = fcNormal; + fs = roundSignificandWithExponent(decSignificand, partCount, + D.exponent, rounding_mode); + + delete [] decSignificand; + } + + return fs; +} + APFloat::opStatus APFloat::convertFromString(const char *p, roundingMode rounding_mode) { + assertArithmeticOK(*semantics); + /* Handle a leading minus sign. */ if(*p == '-') sign = 1, p++; @@ -1687,9 +2231,8 @@ APFloat::convertFromString(const char *p, roundingMode rounding_mode) if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X')) return convertFromHexadecimalString(p + 2, rounding_mode); - - assert(0 && "Decimal to binary conversions not yet implemented"); - abort(); + else + return convertFromDecimalString(p, rounding_mode); } /* Write out a hexadecimal representation of the floating point value @@ -1722,6 +2265,8 @@ APFloat::convertToHexString(char *dst, unsigned int hexDigits, { char *p; + assertArithmeticOK(*semantics); + p = dst; if (sign) *dst++ = '-'; @@ -1897,7 +2442,7 @@ APFloat::getHashValue() const APInt APFloat::convertF80LongDoubleAPFloatToAPInt() const { - assert(semantics == (const llvm::fltSemantics* const)&x87DoubleExtended); + assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended); assert (partCount()==2); uint64_t myexponent, mysignificand; @@ -1927,6 +2472,51 @@ APFloat::convertF80LongDoubleAPFloatToAPInt() const return APInt(80, 2, words); } +APInt +APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const +{ + assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble); + assert (partCount()==2); + + uint64_t myexponent, mysignificand, myexponent2, mysignificand2; + + if (category==fcNormal) { + myexponent = exponent + 1023; //bias + myexponent2 = exponent2 + 1023; + mysignificand = significandParts()[0]; + mysignificand2 = significandParts()[1]; + if (myexponent==1 && !(mysignificand & 0x10000000000000LL)) + myexponent = 0; // denormal + if (myexponent2==1 && !(mysignificand2 & 0x10000000000000LL)) + myexponent2 = 0; // denormal + } else if (category==fcZero) { + myexponent = 0; + mysignificand = 0; + myexponent2 = 0; + mysignificand2 = 0; + } else if (category==fcInfinity) { + myexponent = 0x7ff; + myexponent2 = 0; + mysignificand = 0; + mysignificand2 = 0; + } else { + assert(category == fcNaN && "Unknown category"); + myexponent = 0x7ff; + mysignificand = significandParts()[0]; + myexponent2 = exponent2; + mysignificand2 = significandParts()[1]; + } + + uint64_t words[2]; + words[0] = (((uint64_t)sign & 1) << 63) | + ((myexponent & 0x7ff) << 52) | + (mysignificand & 0xfffffffffffffLL); + words[1] = (((uint64_t)sign2 & 1) << 63) | + ((myexponent2 & 0x7ff) << 52) | + (mysignificand2 & 0xfffffffffffffLL); + return APInt(128, 2, words); +} + APInt APFloat::convertDoubleAPFloatToAPInt() const { @@ -1968,7 +2558,7 @@ APFloat::convertFloatAPFloatToAPInt() const if (category==fcNormal) { myexponent = exponent+127; //bias mysignificand = *significandParts(); - if (myexponent == 1 && !(mysignificand & 0x400000)) + if (myexponent == 1 && !(mysignificand & 0x800000)) myexponent = 0; // denormal } else if (category==fcZero) { myexponent = 0; @@ -1986,16 +2576,23 @@ APFloat::convertFloatAPFloatToAPInt() const (mysignificand & 0x7fffff))); } +// This function creates an APInt that is just a bit map of the floating +// point constant as it would appear in memory. It is not a conversion, +// and treating the result as a normal integer is unlikely to be useful. + APInt APFloat::convertToAPInt() const { - if (semantics == (const llvm::fltSemantics* const)&IEEEsingle) + if (semantics == (const llvm::fltSemantics*)&IEEEsingle) return convertFloatAPFloatToAPInt(); - if (semantics == (const llvm::fltSemantics* const)&IEEEdouble) + if (semantics == (const llvm::fltSemantics*)&IEEEdouble) return convertDoubleAPFloatToAPInt(); - assert(semantics == (const llvm::fltSemantics* const)&x87DoubleExtended && + if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble) + return convertPPCDoubleDoubleAPFloatToAPInt(); + + assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended && "unknown format!"); return convertF80LongDoubleAPFloatToAPInt(); } @@ -2003,7 +2600,7 @@ APFloat::convertToAPInt() const float APFloat::convertToFloat() const { - assert(semantics == (const llvm::fltSemantics* const)&IEEEsingle); + assert(semantics == (const llvm::fltSemantics*)&IEEEsingle); APInt api = convertToAPInt(); return api.bitsToFloat(); } @@ -2011,7 +2608,7 @@ APFloat::convertToFloat() const double APFloat::convertToDouble() const { - assert(semantics == (const llvm::fltSemantics* const)&IEEEdouble); + assert(semantics == (const llvm::fltSemantics*)&IEEEdouble); APInt api = convertToAPInt(); return api.bitsToDouble(); } @@ -2056,6 +2653,56 @@ APFloat::initFromF80LongDoubleAPInt(const APInt &api) } } +void +APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api) +{ + assert(api.getBitWidth()==128); + uint64_t i1 = api.getRawData()[0]; + uint64_t i2 = api.getRawData()[1]; + uint64_t myexponent = (i1 >> 52) & 0x7ff; + uint64_t mysignificand = i1 & 0xfffffffffffffLL; + uint64_t myexponent2 = (i2 >> 52) & 0x7ff; + uint64_t mysignificand2 = i2 & 0xfffffffffffffLL; + + initialize(&APFloat::PPCDoubleDouble); + assert(partCount()==2); + + sign = i1>>63; + sign2 = i2>>63; + if (myexponent==0 && mysignificand==0) { + // exponent, significand meaningless + // exponent2 and significand2 are required to be 0; we don't check + category = fcZero; + } else if (myexponent==0x7ff && mysignificand==0) { + // exponent, significand meaningless + // exponent2 and significand2 are required to be 0; we don't check + category = fcInfinity; + } else if (myexponent==0x7ff && mysignificand!=0) { + // exponent meaningless. So is the whole second word, but keep it + // for determinism. + category = fcNaN; + exponent2 = myexponent2; + significandParts()[0] = mysignificand; + significandParts()[1] = mysignificand2; + } else { + category = fcNormal; + // Note there is no category2; the second word is treated as if it is + // fcNormal, although it might be something else considered by itself. + exponent = myexponent - 1023; + exponent2 = myexponent2 - 1023; + significandParts()[0] = mysignificand; + significandParts()[1] = mysignificand2; + if (myexponent==0) // denormal + exponent = -1022; + else + significandParts()[0] |= 0x10000000000000LL; // integer bit + if (myexponent2==0) + exponent2 = -1022; + else + significandParts()[1] |= 0x10000000000000LL; // integer bit + } +} + void APFloat::initFromDoubleAPInt(const APInt &api) { @@ -2123,11 +2770,11 @@ APFloat::initFromFloatAPInt(const APInt & api) } /// Treat api as containing the bits of a floating point number. Currently -/// we infer the floating point type from the size of the APInt. FIXME: This -/// breaks when we get to PPC128 and IEEE128 (but both cannot exist in the -/// same compile...) +/// we infer the floating point type from the size of the APInt. The +/// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful +/// when the size is anything else). void -APFloat::initFromAPInt(const APInt& api) +APFloat::initFromAPInt(const APInt& api, bool isIEEE) { if (api.getBitWidth() == 32) return initFromFloatAPInt(api); @@ -2135,13 +2782,15 @@ APFloat::initFromAPInt(const APInt& api) return initFromDoubleAPInt(api); else if (api.getBitWidth()==80) return initFromF80LongDoubleAPInt(api); + else if (api.getBitWidth()==128 && !isIEEE) + return initFromPPCDoubleDoubleAPInt(api); else assert(0); } -APFloat::APFloat(const APInt& api) +APFloat::APFloat(const APInt& api, bool isIEEE) { - initFromAPInt(api); + initFromAPInt(api, isIEEE); } APFloat::APFloat(float f)