1 //===-- APFloat.cpp - Implement APFloat class -----------------------------===//
3 // The LLVM Compiler Infrastructure
5 // This file is distributed under the University of Illinois Open Source
6 // License. See LICENSE.TXT for details.
8 //===----------------------------------------------------------------------===//
10 // This file implements a class to represent arbitrary precision floating
11 // point values and provide a variety of arithmetic operations on them.
13 //===----------------------------------------------------------------------===//
15 #include "llvm/ADT/APFloat.h"
16 #include "llvm/ADT/StringRef.h"
17 #include "llvm/ADT/FoldingSet.h"
18 #include "llvm/Support/ErrorHandling.h"
19 #include "llvm/Support/MathExtras.h"
25 #define convolve(lhs, rhs) ((lhs) * 4 + (rhs))
27 /* Assumed in hexadecimal significand parsing, and conversion to
28 hexadecimal strings. */
29 #define COMPILE_TIME_ASSERT(cond) extern int CTAssert[(cond) ? 1 : -1]
30 COMPILE_TIME_ASSERT(integerPartWidth % 4 == 0);
34 /* Represents floating point arithmetic semantics. */
36 /* The largest E such that 2^E is representable; this matches the
37 definition of IEEE 754. */
38 exponent_t maxExponent;
40 /* The smallest E such that 2^E is a normalized number; this
41 matches the definition of IEEE 754. */
42 exponent_t minExponent;
44 /* Number of bits in the significand. This includes the integer
46 unsigned int precision;
48 /* True if arithmetic is supported. */
49 unsigned int arithmeticOK;
52 const fltSemantics APFloat::IEEEhalf = { 15, -14, 11, true };
53 const fltSemantics APFloat::IEEEsingle = { 127, -126, 24, true };
54 const fltSemantics APFloat::IEEEdouble = { 1023, -1022, 53, true };
55 const fltSemantics APFloat::IEEEquad = { 16383, -16382, 113, true };
56 const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64, true };
57 const fltSemantics APFloat::Bogus = { 0, 0, 0, true };
59 // The PowerPC format consists of two doubles. It does not map cleanly
60 // onto the usual format above. For now only storage of constants of
61 // this type is supported, no arithmetic.
62 const fltSemantics APFloat::PPCDoubleDouble = { 1023, -1022, 106, false };
64 /* A tight upper bound on number of parts required to hold the value
67 power * 815 / (351 * integerPartWidth) + 1
69 However, whilst the result may require only this many parts,
70 because we are multiplying two values to get it, the
71 multiplication may require an extra part with the excess part
72 being zero (consider the trivial case of 1 * 1, tcFullMultiply
73 requires two parts to hold the single-part result). So we add an
74 extra one to guarantee enough space whilst multiplying. */
75 const unsigned int maxExponent = 16383;
76 const unsigned int maxPrecision = 113;
77 const unsigned int maxPowerOfFiveExponent = maxExponent + maxPrecision - 1;
78 const unsigned int maxPowerOfFiveParts = 2 + ((maxPowerOfFiveExponent * 815)
79 / (351 * integerPartWidth));
82 /* A bunch of private, handy routines. */
84 static inline unsigned int
85 partCountForBits(unsigned int bits)
87 return ((bits) + integerPartWidth - 1) / integerPartWidth;
90 /* Returns 0U-9U. Return values >= 10U are not digits. */
91 static inline unsigned int
92 decDigitValue(unsigned int c)
98 hexDigitValue(unsigned int c)
118 assertArithmeticOK(const llvm::fltSemantics &semantics) {
119 assert(semantics.arithmeticOK &&
120 "Compile-time arithmetic does not support these semantics");
123 /* Return the value of a decimal exponent of the form
126 If the exponent overflows, returns a large exponent with the
129 readExponent(StringRef::iterator begin, StringRef::iterator end)
132 unsigned int absExponent;
133 const unsigned int overlargeExponent = 24000; /* FIXME. */
134 StringRef::iterator p = begin;
136 assert(p != end && "Exponent has no digits");
138 isNegative = (*p == '-');
139 if (*p == '-' || *p == '+') {
141 assert(p != end && "Exponent has no digits");
144 absExponent = decDigitValue(*p++);
145 assert(absExponent < 10U && "Invalid character in exponent");
147 for (; p != end; ++p) {
150 value = decDigitValue(*p);
151 assert(value < 10U && "Invalid character in exponent");
153 value += absExponent * 10;
154 if (absExponent >= overlargeExponent) {
155 absExponent = overlargeExponent;
161 assert(p == end && "Invalid exponent in exponent");
164 return -(int) absExponent;
166 return (int) absExponent;
169 /* This is ugly and needs cleaning up, but I don't immediately see
170 how whilst remaining safe. */
172 totalExponent(StringRef::iterator p, StringRef::iterator end,
173 int exponentAdjustment)
175 int unsignedExponent;
176 bool negative, overflow;
179 assert(p != end && "Exponent has no digits");
181 negative = *p == '-';
182 if (*p == '-' || *p == '+') {
184 assert(p != end && "Exponent has no digits");
187 unsignedExponent = 0;
189 for (; p != end; ++p) {
192 value = decDigitValue(*p);
193 assert(value < 10U && "Invalid character in exponent");
195 unsignedExponent = unsignedExponent * 10 + value;
196 if (unsignedExponent > 65535)
200 if (exponentAdjustment > 65535 || exponentAdjustment < -65536)
204 exponent = unsignedExponent;
206 exponent = -exponent;
207 exponent += exponentAdjustment;
208 if (exponent > 65535 || exponent < -65536)
213 exponent = negative ? -65536: 65535;
218 static StringRef::iterator
219 skipLeadingZeroesAndAnyDot(StringRef::iterator begin, StringRef::iterator end,
220 StringRef::iterator *dot)
222 StringRef::iterator p = begin;
224 while (*p == '0' && p != end)
230 assert(end - begin != 1 && "Significand has no digits");
232 while (*p == '0' && p != end)
239 /* Given a normal decimal floating point number of the form
243 where the decimal point and exponent are optional, fill out the
244 structure D. Exponent is appropriate if the significand is
245 treated as an integer, and normalizedExponent if the significand
246 is taken to have the decimal point after a single leading
249 If the value is zero, V->firstSigDigit points to a non-digit, and
250 the return exponent is zero.
253 const char *firstSigDigit;
254 const char *lastSigDigit;
256 int normalizedExponent;
260 interpretDecimal(StringRef::iterator begin, StringRef::iterator end,
263 StringRef::iterator dot = end;
264 StringRef::iterator p = skipLeadingZeroesAndAnyDot (begin, end, &dot);
266 D->firstSigDigit = p;
268 D->normalizedExponent = 0;
270 for (; p != end; ++p) {
272 assert(dot == end && "String contains multiple dots");
277 if (decDigitValue(*p) >= 10U)
282 assert((*p == 'e' || *p == 'E') && "Invalid character in significand");
283 assert(p != begin && "Significand has no digits");
284 assert((dot == end || p - begin != 1) && "Significand has no digits");
286 /* p points to the first non-digit in the string */
287 D->exponent = readExponent(p + 1, end);
289 /* Implied decimal point? */
294 /* If number is all zeroes accept any exponent. */
295 if (p != D->firstSigDigit) {
296 /* Drop insignificant trailing zeroes. */
301 while (p != begin && *p == '0');
302 while (p != begin && *p == '.');
305 /* Adjust the exponents for any decimal point. */
306 D->exponent += static_cast<exponent_t>((dot - p) - (dot > p));
307 D->normalizedExponent = (D->exponent +
308 static_cast<exponent_t>((p - D->firstSigDigit)
309 - (dot > D->firstSigDigit && dot < p)));
315 /* Return the trailing fraction of a hexadecimal number.
316 DIGITVALUE is the first hex digit of the fraction, P points to
319 trailingHexadecimalFraction(StringRef::iterator p, StringRef::iterator end,
320 unsigned int digitValue)
322 unsigned int hexDigit;
324 /* If the first trailing digit isn't 0 or 8 we can work out the
325 fraction immediately. */
327 return lfMoreThanHalf;
328 else if (digitValue < 8 && digitValue > 0)
329 return lfLessThanHalf;
331 /* Otherwise we need to find the first non-zero digit. */
335 assert(p != end && "Invalid trailing hexadecimal fraction!");
337 hexDigit = hexDigitValue(*p);
339 /* If we ran off the end it is exactly zero or one-half, otherwise
342 return digitValue == 0 ? lfExactlyZero: lfExactlyHalf;
344 return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf;
347 /* Return the fraction lost were a bignum truncated losing the least
348 significant BITS bits. */
350 lostFractionThroughTruncation(const integerPart *parts,
351 unsigned int partCount,
356 lsb = APInt::tcLSB(parts, partCount);
358 /* Note this is guaranteed true if bits == 0, or LSB == -1U. */
360 return lfExactlyZero;
362 return lfExactlyHalf;
363 if (bits <= partCount * integerPartWidth &&
364 APInt::tcExtractBit(parts, bits - 1))
365 return lfMoreThanHalf;
367 return lfLessThanHalf;
370 /* Shift DST right BITS bits noting lost fraction. */
372 shiftRight(integerPart *dst, unsigned int parts, unsigned int bits)
374 lostFraction lost_fraction;
376 lost_fraction = lostFractionThroughTruncation(dst, parts, bits);
378 APInt::tcShiftRight(dst, parts, bits);
380 return lost_fraction;
383 /* Combine the effect of two lost fractions. */
385 combineLostFractions(lostFraction moreSignificant,
386 lostFraction lessSignificant)
388 if (lessSignificant != lfExactlyZero) {
389 if (moreSignificant == lfExactlyZero)
390 moreSignificant = lfLessThanHalf;
391 else if (moreSignificant == lfExactlyHalf)
392 moreSignificant = lfMoreThanHalf;
395 return moreSignificant;
398 /* The error from the true value, in half-ulps, on multiplying two
399 floating point numbers, which differ from the value they
400 approximate by at most HUE1 and HUE2 half-ulps, is strictly less
401 than the returned value.
403 See "How to Read Floating Point Numbers Accurately" by William D
406 HUerrBound(bool inexactMultiply, unsigned int HUerr1, unsigned int HUerr2)
408 assert(HUerr1 < 2 || HUerr2 < 2 || (HUerr1 + HUerr2 < 8));
410 if (HUerr1 + HUerr2 == 0)
411 return inexactMultiply * 2; /* <= inexactMultiply half-ulps. */
413 return inexactMultiply + 2 * (HUerr1 + HUerr2);
416 /* The number of ulps from the boundary (zero, or half if ISNEAREST)
417 when the least significant BITS are truncated. BITS cannot be
420 ulpsFromBoundary(const integerPart *parts, unsigned int bits, bool isNearest)
422 unsigned int count, partBits;
423 integerPart part, boundary;
428 count = bits / integerPartWidth;
429 partBits = bits % integerPartWidth + 1;
431 part = parts[count] & (~(integerPart) 0 >> (integerPartWidth - partBits));
434 boundary = (integerPart) 1 << (partBits - 1);
439 if (part - boundary <= boundary - part)
440 return part - boundary;
442 return boundary - part;
445 if (part == boundary) {
448 return ~(integerPart) 0; /* A lot. */
451 } else if (part == boundary - 1) {
454 return ~(integerPart) 0; /* A lot. */
459 return ~(integerPart) 0; /* A lot. */
462 /* Place pow(5, power) in DST, and return the number of parts used.
463 DST must be at least one part larger than size of the answer. */
465 powerOf5(integerPart *dst, unsigned int power)
467 static const integerPart firstEightPowers[] = { 1, 5, 25, 125, 625, 3125,
469 integerPart pow5s[maxPowerOfFiveParts * 2 + 5];
470 pow5s[0] = 78125 * 5;
472 unsigned int partsCount[16] = { 1 };
473 integerPart scratch[maxPowerOfFiveParts], *p1, *p2, *pow5;
475 assert(power <= maxExponent);
480 *p1 = firstEightPowers[power & 7];
486 for (unsigned int n = 0; power; power >>= 1, n++) {
491 /* Calculate pow(5,pow(2,n+3)) if we haven't yet. */
493 pc = partsCount[n - 1];
494 APInt::tcFullMultiply(pow5, pow5 - pc, pow5 - pc, pc, pc);
496 if (pow5[pc - 1] == 0)
504 APInt::tcFullMultiply(p2, p1, pow5, result, pc);
506 if (p2[result - 1] == 0)
509 /* Now result is in p1 with partsCount parts and p2 is scratch
511 tmp = p1, p1 = p2, p2 = tmp;
518 APInt::tcAssign(dst, p1, result);
523 /* Zero at the end to avoid modular arithmetic when adding one; used
524 when rounding up during hexadecimal output. */
525 static const char hexDigitsLower[] = "0123456789abcdef0";
526 static const char hexDigitsUpper[] = "0123456789ABCDEF0";
527 static const char infinityL[] = "infinity";
528 static const char infinityU[] = "INFINITY";
529 static const char NaNL[] = "nan";
530 static const char NaNU[] = "NAN";
532 /* Write out an integerPart in hexadecimal, starting with the most
533 significant nibble. Write out exactly COUNT hexdigits, return
536 partAsHex (char *dst, integerPart part, unsigned int count,
537 const char *hexDigitChars)
539 unsigned int result = count;
541 assert(count != 0 && count <= integerPartWidth / 4);
543 part >>= (integerPartWidth - 4 * count);
545 dst[count] = hexDigitChars[part & 0xf];
552 /* Write out an unsigned decimal integer. */
554 writeUnsignedDecimal (char *dst, unsigned int n)
570 /* Write out a signed decimal integer. */
572 writeSignedDecimal (char *dst, int value)
576 dst = writeUnsignedDecimal(dst, -(unsigned) value);
578 dst = writeUnsignedDecimal(dst, value);
585 APFloat::initialize(const fltSemantics *ourSemantics)
589 semantics = ourSemantics;
592 significand.parts = new integerPart[count];
596 APFloat::freeSignificand()
599 delete [] significand.parts;
603 APFloat::assign(const APFloat &rhs)
605 assert(semantics == rhs.semantics);
608 category = rhs.category;
609 exponent = rhs.exponent;
611 exponent2 = rhs.exponent2;
612 if (category == fcNormal || category == fcNaN)
613 copySignificand(rhs);
617 APFloat::copySignificand(const APFloat &rhs)
619 assert(category == fcNormal || category == fcNaN);
620 assert(rhs.partCount() >= partCount());
622 APInt::tcAssign(significandParts(), rhs.significandParts(),
626 /* Make this number a NaN, with an arbitrary but deterministic value
627 for the significand. If double or longer, this is a signalling NaN,
628 which may not be ideal. If float, this is QNaN(0). */
629 void APFloat::makeNaN(bool SNaN, bool Negative, const APInt *fill)
634 integerPart *significand = significandParts();
635 unsigned numParts = partCount();
637 // Set the significand bits to the fill.
638 if (!fill || fill->getNumWords() < numParts)
639 APInt::tcSet(significand, 0, numParts);
641 APInt::tcAssign(significand, fill->getRawData(),
642 std::min(fill->getNumWords(), numParts));
644 // Zero out the excess bits of the significand.
645 unsigned bitsToPreserve = semantics->precision - 1;
646 unsigned part = bitsToPreserve / 64;
647 bitsToPreserve %= 64;
648 significand[part] &= ((1ULL << bitsToPreserve) - 1);
649 for (part++; part != numParts; ++part)
650 significand[part] = 0;
653 unsigned QNaNBit = semantics->precision - 2;
656 // We always have to clear the QNaN bit to make it an SNaN.
657 APInt::tcClearBit(significand, QNaNBit);
659 // If there are no bits set in the payload, we have to set
660 // *something* to make it a NaN instead of an infinity;
661 // conventionally, this is the next bit down from the QNaN bit.
662 if (APInt::tcIsZero(significand, numParts))
663 APInt::tcSetBit(significand, QNaNBit - 1);
665 // We always have to set the QNaN bit to make it a QNaN.
666 APInt::tcSetBit(significand, QNaNBit);
669 // For x87 extended precision, we want to make a NaN, not a
670 // pseudo-NaN. Maybe we should expose the ability to make
672 if (semantics == &APFloat::x87DoubleExtended)
673 APInt::tcSetBit(significand, QNaNBit + 1);
676 APFloat APFloat::makeNaN(const fltSemantics &Sem, bool SNaN, bool Negative,
678 APFloat value(Sem, uninitialized);
679 value.makeNaN(SNaN, Negative, fill);
684 APFloat::operator=(const APFloat &rhs)
687 if (semantics != rhs.semantics) {
689 initialize(rhs.semantics);
698 APFloat::bitwiseIsEqual(const APFloat &rhs) const {
701 if (semantics != rhs.semantics ||
702 category != rhs.category ||
705 if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble &&
708 if (category==fcZero || category==fcInfinity)
710 else if (category==fcNormal && exponent!=rhs.exponent)
712 else if (semantics==(const llvm::fltSemantics*)&PPCDoubleDouble &&
713 exponent2!=rhs.exponent2)
717 const integerPart* p=significandParts();
718 const integerPart* q=rhs.significandParts();
719 for (; i>0; i--, p++, q++) {
727 APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value)
729 assertArithmeticOK(ourSemantics);
730 initialize(&ourSemantics);
733 exponent = ourSemantics.precision - 1;
734 significandParts()[0] = value;
735 normalize(rmNearestTiesToEven, lfExactlyZero);
738 APFloat::APFloat(const fltSemantics &ourSemantics) {
739 assertArithmeticOK(ourSemantics);
740 initialize(&ourSemantics);
745 APFloat::APFloat(const fltSemantics &ourSemantics, uninitializedTag tag) {
746 assertArithmeticOK(ourSemantics);
747 // Allocates storage if necessary but does not initialize it.
748 initialize(&ourSemantics);
751 APFloat::APFloat(const fltSemantics &ourSemantics,
752 fltCategory ourCategory, bool negative)
754 assertArithmeticOK(ourSemantics);
755 initialize(&ourSemantics);
756 category = ourCategory;
758 if (category == fcNormal)
760 else if (ourCategory == fcNaN)
764 APFloat::APFloat(const fltSemantics &ourSemantics, const StringRef& text)
766 assertArithmeticOK(ourSemantics);
767 initialize(&ourSemantics);
768 convertFromString(text, rmNearestTiesToEven);
771 APFloat::APFloat(const APFloat &rhs)
773 initialize(rhs.semantics);
782 // Profile - This method 'profiles' an APFloat for use with FoldingSet.
783 void APFloat::Profile(FoldingSetNodeID& ID) const {
784 ID.Add(bitcastToAPInt());
788 APFloat::partCount() const
790 return partCountForBits(semantics->precision + 1);
794 APFloat::semanticsPrecision(const fltSemantics &semantics)
796 return semantics.precision;
800 APFloat::significandParts() const
802 return const_cast<APFloat *>(this)->significandParts();
806 APFloat::significandParts()
808 assert(category == fcNormal || category == fcNaN);
811 return significand.parts;
813 return &significand.part;
817 APFloat::zeroSignificand()
820 APInt::tcSet(significandParts(), 0, partCount());
823 /* Increment an fcNormal floating point number's significand. */
825 APFloat::incrementSignificand()
829 carry = APInt::tcIncrement(significandParts(), partCount());
831 /* Our callers should never cause us to overflow. */
835 /* Add the significand of the RHS. Returns the carry flag. */
837 APFloat::addSignificand(const APFloat &rhs)
841 parts = significandParts();
843 assert(semantics == rhs.semantics);
844 assert(exponent == rhs.exponent);
846 return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
849 /* Subtract the significand of the RHS with a borrow flag. Returns
852 APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
856 parts = significandParts();
858 assert(semantics == rhs.semantics);
859 assert(exponent == rhs.exponent);
861 return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
865 /* Multiply the significand of the RHS. If ADDEND is non-NULL, add it
866 on to the full-precision result of the multiplication. Returns the
869 APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
871 unsigned int omsb; // One, not zero, based MSB.
872 unsigned int partsCount, newPartsCount, precision;
873 integerPart *lhsSignificand;
874 integerPart scratch[4];
875 integerPart *fullSignificand;
876 lostFraction lost_fraction;
879 assert(semantics == rhs.semantics);
881 precision = semantics->precision;
882 newPartsCount = partCountForBits(precision * 2);
884 if (newPartsCount > 4)
885 fullSignificand = new integerPart[newPartsCount];
887 fullSignificand = scratch;
889 lhsSignificand = significandParts();
890 partsCount = partCount();
892 APInt::tcFullMultiply(fullSignificand, lhsSignificand,
893 rhs.significandParts(), partsCount, partsCount);
895 lost_fraction = lfExactlyZero;
896 omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
897 exponent += rhs.exponent;
900 Significand savedSignificand = significand;
901 const fltSemantics *savedSemantics = semantics;
902 fltSemantics extendedSemantics;
904 unsigned int extendedPrecision;
906 /* Normalize our MSB. */
907 extendedPrecision = precision + precision - 1;
908 if (omsb != extendedPrecision) {
909 APInt::tcShiftLeft(fullSignificand, newPartsCount,
910 extendedPrecision - omsb);
911 exponent -= extendedPrecision - omsb;
914 /* Create new semantics. */
915 extendedSemantics = *semantics;
916 extendedSemantics.precision = extendedPrecision;
918 if (newPartsCount == 1)
919 significand.part = fullSignificand[0];
921 significand.parts = fullSignificand;
922 semantics = &extendedSemantics;
924 APFloat extendedAddend(*addend);
925 status = extendedAddend.convert(extendedSemantics, rmTowardZero, &ignored);
926 assert(status == opOK);
927 lost_fraction = addOrSubtractSignificand(extendedAddend, false);
929 /* Restore our state. */
930 if (newPartsCount == 1)
931 fullSignificand[0] = significand.part;
932 significand = savedSignificand;
933 semantics = savedSemantics;
935 omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
938 exponent -= (precision - 1);
940 if (omsb > precision) {
941 unsigned int bits, significantParts;
944 bits = omsb - precision;
945 significantParts = partCountForBits(omsb);
946 lf = shiftRight(fullSignificand, significantParts, bits);
947 lost_fraction = combineLostFractions(lf, lost_fraction);
951 APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
953 if (newPartsCount > 4)
954 delete [] fullSignificand;
956 return lost_fraction;
959 /* Multiply the significands of LHS and RHS to DST. */
961 APFloat::divideSignificand(const APFloat &rhs)
963 unsigned int bit, i, partsCount;
964 const integerPart *rhsSignificand;
965 integerPart *lhsSignificand, *dividend, *divisor;
966 integerPart scratch[4];
967 lostFraction lost_fraction;
969 assert(semantics == rhs.semantics);
971 lhsSignificand = significandParts();
972 rhsSignificand = rhs.significandParts();
973 partsCount = partCount();
976 dividend = new integerPart[partsCount * 2];
980 divisor = dividend + partsCount;
982 /* Copy the dividend and divisor as they will be modified in-place. */
983 for (i = 0; i < partsCount; i++) {
984 dividend[i] = lhsSignificand[i];
985 divisor[i] = rhsSignificand[i];
986 lhsSignificand[i] = 0;
989 exponent -= rhs.exponent;
991 unsigned int precision = semantics->precision;
993 /* Normalize the divisor. */
994 bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
997 APInt::tcShiftLeft(divisor, partsCount, bit);
1000 /* Normalize the dividend. */
1001 bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
1004 APInt::tcShiftLeft(dividend, partsCount, bit);
1007 /* Ensure the dividend >= divisor initially for the loop below.
1008 Incidentally, this means that the division loop below is
1009 guaranteed to set the integer bit to one. */
1010 if (APInt::tcCompare(dividend, divisor, partsCount) < 0) {
1012 APInt::tcShiftLeft(dividend, partsCount, 1);
1013 assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
1016 /* Long division. */
1017 for (bit = precision; bit; bit -= 1) {
1018 if (APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
1019 APInt::tcSubtract(dividend, divisor, 0, partsCount);
1020 APInt::tcSetBit(lhsSignificand, bit - 1);
1023 APInt::tcShiftLeft(dividend, partsCount, 1);
1026 /* Figure out the lost fraction. */
1027 int cmp = APInt::tcCompare(dividend, divisor, partsCount);
1030 lost_fraction = lfMoreThanHalf;
1032 lost_fraction = lfExactlyHalf;
1033 else if (APInt::tcIsZero(dividend, partsCount))
1034 lost_fraction = lfExactlyZero;
1036 lost_fraction = lfLessThanHalf;
1041 return lost_fraction;
1045 APFloat::significandMSB() const
1047 return APInt::tcMSB(significandParts(), partCount());
1051 APFloat::significandLSB() const
1053 return APInt::tcLSB(significandParts(), partCount());
1056 /* Note that a zero result is NOT normalized to fcZero. */
1058 APFloat::shiftSignificandRight(unsigned int bits)
1060 /* Our exponent should not overflow. */
1061 assert((exponent_t) (exponent + bits) >= exponent);
1065 return shiftRight(significandParts(), partCount(), bits);
1068 /* Shift the significand left BITS bits, subtract BITS from its exponent. */
1070 APFloat::shiftSignificandLeft(unsigned int bits)
1072 assert(bits < semantics->precision);
1075 unsigned int partsCount = partCount();
1077 APInt::tcShiftLeft(significandParts(), partsCount, bits);
1080 assert(!APInt::tcIsZero(significandParts(), partsCount));
1085 APFloat::compareAbsoluteValue(const APFloat &rhs) const
1089 assert(semantics == rhs.semantics);
1090 assert(category == fcNormal);
1091 assert(rhs.category == fcNormal);
1093 compare = exponent - rhs.exponent;
1095 /* If exponents are equal, do an unsigned bignum comparison of the
1098 compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
1102 return cmpGreaterThan;
1103 else if (compare < 0)
1109 /* Handle overflow. Sign is preserved. We either become infinity or
1110 the largest finite number. */
1112 APFloat::handleOverflow(roundingMode rounding_mode)
1115 if (rounding_mode == rmNearestTiesToEven ||
1116 rounding_mode == rmNearestTiesToAway ||
1117 (rounding_mode == rmTowardPositive && !sign) ||
1118 (rounding_mode == rmTowardNegative && sign)) {
1119 category = fcInfinity;
1120 return (opStatus) (opOverflow | opInexact);
1123 /* Otherwise we become the largest finite number. */
1124 category = fcNormal;
1125 exponent = semantics->maxExponent;
1126 APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
1127 semantics->precision);
1132 /* Returns TRUE if, when truncating the current number, with BIT the
1133 new LSB, with the given lost fraction and rounding mode, the result
1134 would need to be rounded away from zero (i.e., by increasing the
1135 signficand). This routine must work for fcZero of both signs, and
1136 fcNormal numbers. */
1138 APFloat::roundAwayFromZero(roundingMode rounding_mode,
1139 lostFraction lost_fraction,
1140 unsigned int bit) const
1142 /* NaNs and infinities should not have lost fractions. */
1143 assert(category == fcNormal || category == fcZero);
1145 /* Current callers never pass this so we don't handle it. */
1146 assert(lost_fraction != lfExactlyZero);
1148 switch (rounding_mode) {
1150 llvm_unreachable(0);
1152 case rmNearestTiesToAway:
1153 return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
1155 case rmNearestTiesToEven:
1156 if (lost_fraction == lfMoreThanHalf)
1159 /* Our zeroes don't have a significand to test. */
1160 if (lost_fraction == lfExactlyHalf && category != fcZero)
1161 return APInt::tcExtractBit(significandParts(), bit);
1168 case rmTowardPositive:
1169 return sign == false;
1171 case rmTowardNegative:
1172 return sign == true;
1177 APFloat::normalize(roundingMode rounding_mode,
1178 lostFraction lost_fraction)
1180 unsigned int omsb; /* One, not zero, based MSB. */
1183 if (category != fcNormal)
1186 /* Before rounding normalize the exponent of fcNormal numbers. */
1187 omsb = significandMSB() + 1;
1190 /* OMSB is numbered from 1. We want to place it in the integer
1191 bit numbered PRECISON if possible, with a compensating change in
1193 exponentChange = omsb - semantics->precision;
1195 /* If the resulting exponent is too high, overflow according to
1196 the rounding mode. */
1197 if (exponent + exponentChange > semantics->maxExponent)
1198 return handleOverflow(rounding_mode);
1200 /* Subnormal numbers have exponent minExponent, and their MSB
1201 is forced based on that. */
1202 if (exponent + exponentChange < semantics->minExponent)
1203 exponentChange = semantics->minExponent - exponent;
1205 /* Shifting left is easy as we don't lose precision. */
1206 if (exponentChange < 0) {
1207 assert(lost_fraction == lfExactlyZero);
1209 shiftSignificandLeft(-exponentChange);
1214 if (exponentChange > 0) {
1217 /* Shift right and capture any new lost fraction. */
1218 lf = shiftSignificandRight(exponentChange);
1220 lost_fraction = combineLostFractions(lf, lost_fraction);
1222 /* Keep OMSB up-to-date. */
1223 if (omsb > (unsigned) exponentChange)
1224 omsb -= exponentChange;
1230 /* Now round the number according to rounding_mode given the lost
1233 /* As specified in IEEE 754, since we do not trap we do not report
1234 underflow for exact results. */
1235 if (lost_fraction == lfExactlyZero) {
1236 /* Canonicalize zeroes. */
1243 /* Increment the significand if we're rounding away from zero. */
1244 if (roundAwayFromZero(rounding_mode, lost_fraction, 0)) {
1246 exponent = semantics->minExponent;
1248 incrementSignificand();
1249 omsb = significandMSB() + 1;
1251 /* Did the significand increment overflow? */
1252 if (omsb == (unsigned) semantics->precision + 1) {
1253 /* Renormalize by incrementing the exponent and shifting our
1254 significand right one. However if we already have the
1255 maximum exponent we overflow to infinity. */
1256 if (exponent == semantics->maxExponent) {
1257 category = fcInfinity;
1259 return (opStatus) (opOverflow | opInexact);
1262 shiftSignificandRight(1);
1268 /* The normal case - we were and are not denormal, and any
1269 significand increment above didn't overflow. */
1270 if (omsb == semantics->precision)
1273 /* We have a non-zero denormal. */
1274 assert(omsb < semantics->precision);
1276 /* Canonicalize zeroes. */
1280 /* The fcZero case is a denormal that underflowed to zero. */
1281 return (opStatus) (opUnderflow | opInexact);
1285 APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
1287 switch (convolve(category, rhs.category)) {
1289 llvm_unreachable(0);
1291 case convolve(fcNaN, fcZero):
1292 case convolve(fcNaN, fcNormal):
1293 case convolve(fcNaN, fcInfinity):
1294 case convolve(fcNaN, fcNaN):
1295 case convolve(fcNormal, fcZero):
1296 case convolve(fcInfinity, fcNormal):
1297 case convolve(fcInfinity, fcZero):
1300 case convolve(fcZero, fcNaN):
1301 case convolve(fcNormal, fcNaN):
1302 case convolve(fcInfinity, fcNaN):
1304 copySignificand(rhs);
1307 case convolve(fcNormal, fcInfinity):
1308 case convolve(fcZero, fcInfinity):
1309 category = fcInfinity;
1310 sign = rhs.sign ^ subtract;
1313 case convolve(fcZero, fcNormal):
1315 sign = rhs.sign ^ subtract;
1318 case convolve(fcZero, fcZero):
1319 /* Sign depends on rounding mode; handled by caller. */
1322 case convolve(fcInfinity, fcInfinity):
1323 /* Differently signed infinities can only be validly
1325 if (((sign ^ rhs.sign)!=0) != subtract) {
1332 case convolve(fcNormal, fcNormal):
1337 /* Add or subtract two normal numbers. */
1339 APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
1342 lostFraction lost_fraction;
1345 /* Determine if the operation on the absolute values is effectively
1346 an addition or subtraction. */
1347 subtract ^= (sign ^ rhs.sign) ? true : false;
1349 /* Are we bigger exponent-wise than the RHS? */
1350 bits = exponent - rhs.exponent;
1352 /* Subtraction is more subtle than one might naively expect. */
1354 APFloat temp_rhs(rhs);
1358 reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
1359 lost_fraction = lfExactlyZero;
1360 } else if (bits > 0) {
1361 lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
1362 shiftSignificandLeft(1);
1365 lost_fraction = shiftSignificandRight(-bits - 1);
1366 temp_rhs.shiftSignificandLeft(1);
1371 carry = temp_rhs.subtractSignificand
1372 (*this, lost_fraction != lfExactlyZero);
1373 copySignificand(temp_rhs);
1376 carry = subtractSignificand
1377 (temp_rhs, lost_fraction != lfExactlyZero);
1380 /* Invert the lost fraction - it was on the RHS and
1382 if (lost_fraction == lfLessThanHalf)
1383 lost_fraction = lfMoreThanHalf;
1384 else if (lost_fraction == lfMoreThanHalf)
1385 lost_fraction = lfLessThanHalf;
1387 /* The code above is intended to ensure that no borrow is
1392 APFloat temp_rhs(rhs);
1394 lost_fraction = temp_rhs.shiftSignificandRight(bits);
1395 carry = addSignificand(temp_rhs);
1397 lost_fraction = shiftSignificandRight(-bits);
1398 carry = addSignificand(rhs);
1401 /* We have a guard bit; generating a carry cannot happen. */
1405 return lost_fraction;
1409 APFloat::multiplySpecials(const APFloat &rhs)
1411 switch (convolve(category, rhs.category)) {
1413 llvm_unreachable(0);
1415 case convolve(fcNaN, fcZero):
1416 case convolve(fcNaN, fcNormal):
1417 case convolve(fcNaN, fcInfinity):
1418 case convolve(fcNaN, fcNaN):
1421 case convolve(fcZero, fcNaN):
1422 case convolve(fcNormal, fcNaN):
1423 case convolve(fcInfinity, fcNaN):
1425 copySignificand(rhs);
1428 case convolve(fcNormal, fcInfinity):
1429 case convolve(fcInfinity, fcNormal):
1430 case convolve(fcInfinity, fcInfinity):
1431 category = fcInfinity;
1434 case convolve(fcZero, fcNormal):
1435 case convolve(fcNormal, fcZero):
1436 case convolve(fcZero, fcZero):
1440 case convolve(fcZero, fcInfinity):
1441 case convolve(fcInfinity, fcZero):
1445 case convolve(fcNormal, fcNormal):
1451 APFloat::divideSpecials(const APFloat &rhs)
1453 switch (convolve(category, rhs.category)) {
1455 llvm_unreachable(0);
1457 case convolve(fcNaN, fcZero):
1458 case convolve(fcNaN, fcNormal):
1459 case convolve(fcNaN, fcInfinity):
1460 case convolve(fcNaN, fcNaN):
1461 case convolve(fcInfinity, fcZero):
1462 case convolve(fcInfinity, fcNormal):
1463 case convolve(fcZero, fcInfinity):
1464 case convolve(fcZero, fcNormal):
1467 case convolve(fcZero, fcNaN):
1468 case convolve(fcNormal, fcNaN):
1469 case convolve(fcInfinity, fcNaN):
1471 copySignificand(rhs);
1474 case convolve(fcNormal, fcInfinity):
1478 case convolve(fcNormal, fcZero):
1479 category = fcInfinity;
1482 case convolve(fcInfinity, fcInfinity):
1483 case convolve(fcZero, fcZero):
1487 case convolve(fcNormal, fcNormal):
1493 APFloat::modSpecials(const APFloat &rhs)
1495 switch (convolve(category, rhs.category)) {
1497 llvm_unreachable(0);
1499 case convolve(fcNaN, fcZero):
1500 case convolve(fcNaN, fcNormal):
1501 case convolve(fcNaN, fcInfinity):
1502 case convolve(fcNaN, fcNaN):
1503 case convolve(fcZero, fcInfinity):
1504 case convolve(fcZero, fcNormal):
1505 case convolve(fcNormal, fcInfinity):
1508 case convolve(fcZero, fcNaN):
1509 case convolve(fcNormal, fcNaN):
1510 case convolve(fcInfinity, fcNaN):
1512 copySignificand(rhs);
1515 case convolve(fcNormal, fcZero):
1516 case convolve(fcInfinity, fcZero):
1517 case convolve(fcInfinity, fcNormal):
1518 case convolve(fcInfinity, fcInfinity):
1519 case convolve(fcZero, fcZero):
1523 case convolve(fcNormal, fcNormal):
1530 APFloat::changeSign()
1532 /* Look mummy, this one's easy. */
1537 APFloat::clearSign()
1539 /* So is this one. */
1544 APFloat::copySign(const APFloat &rhs)
1550 /* Normalized addition or subtraction. */
1552 APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1557 assertArithmeticOK(*semantics);
1559 fs = addOrSubtractSpecials(rhs, subtract);
1561 /* This return code means it was not a simple case. */
1562 if (fs == opDivByZero) {
1563 lostFraction lost_fraction;
1565 lost_fraction = addOrSubtractSignificand(rhs, subtract);
1566 fs = normalize(rounding_mode, lost_fraction);
1568 /* Can only be zero if we lost no fraction. */
1569 assert(category != fcZero || lost_fraction == lfExactlyZero);
1572 /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1573 positive zero unless rounding to minus infinity, except that
1574 adding two like-signed zeroes gives that zero. */
1575 if (category == fcZero) {
1576 if (rhs.category != fcZero || (sign == rhs.sign) == subtract)
1577 sign = (rounding_mode == rmTowardNegative);
1583 /* Normalized addition. */
1585 APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1587 return addOrSubtract(rhs, rounding_mode, false);
1590 /* Normalized subtraction. */
1592 APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1594 return addOrSubtract(rhs, rounding_mode, true);
1597 /* Normalized multiply. */
1599 APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1603 assertArithmeticOK(*semantics);
1605 fs = multiplySpecials(rhs);
1607 if (category == fcNormal) {
1608 lostFraction lost_fraction = multiplySignificand(rhs, 0);
1609 fs = normalize(rounding_mode, lost_fraction);
1610 if (lost_fraction != lfExactlyZero)
1611 fs = (opStatus) (fs | opInexact);
1617 /* Normalized divide. */
1619 APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1623 assertArithmeticOK(*semantics);
1625 fs = divideSpecials(rhs);
1627 if (category == fcNormal) {
1628 lostFraction lost_fraction = divideSignificand(rhs);
1629 fs = normalize(rounding_mode, lost_fraction);
1630 if (lost_fraction != lfExactlyZero)
1631 fs = (opStatus) (fs | opInexact);
1637 /* Normalized remainder. This is not currently correct in all cases. */
1639 APFloat::remainder(const APFloat &rhs)
1643 unsigned int origSign = sign;
1645 assertArithmeticOK(*semantics);
1646 fs = V.divide(rhs, rmNearestTiesToEven);
1647 if (fs == opDivByZero)
1650 int parts = partCount();
1651 integerPart *x = new integerPart[parts];
1653 fs = V.convertToInteger(x, parts * integerPartWidth, true,
1654 rmNearestTiesToEven, &ignored);
1655 if (fs==opInvalidOp)
1658 fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1659 rmNearestTiesToEven);
1660 assert(fs==opOK); // should always work
1662 fs = V.multiply(rhs, rmNearestTiesToEven);
1663 assert(fs==opOK || fs==opInexact); // should not overflow or underflow
1665 fs = subtract(V, rmNearestTiesToEven);
1666 assert(fs==opOK || fs==opInexact); // likewise
1669 sign = origSign; // IEEE754 requires this
1674 /* Normalized llvm frem (C fmod).
1675 This is not currently correct in all cases. */
1677 APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1680 assertArithmeticOK(*semantics);
1681 fs = modSpecials(rhs);
1683 if (category == fcNormal && rhs.category == fcNormal) {
1685 unsigned int origSign = sign;
1687 fs = V.divide(rhs, rmNearestTiesToEven);
1688 if (fs == opDivByZero)
1691 int parts = partCount();
1692 integerPart *x = new integerPart[parts];
1694 fs = V.convertToInteger(x, parts * integerPartWidth, true,
1695 rmTowardZero, &ignored);
1696 if (fs==opInvalidOp)
1699 fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1700 rmNearestTiesToEven);
1701 assert(fs==opOK); // should always work
1703 fs = V.multiply(rhs, rounding_mode);
1704 assert(fs==opOK || fs==opInexact); // should not overflow or underflow
1706 fs = subtract(V, rounding_mode);
1707 assert(fs==opOK || fs==opInexact); // likewise
1710 sign = origSign; // IEEE754 requires this
1716 /* Normalized fused-multiply-add. */
1718 APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1719 const APFloat &addend,
1720 roundingMode rounding_mode)
1724 assertArithmeticOK(*semantics);
1726 /* Post-multiplication sign, before addition. */
1727 sign ^= multiplicand.sign;
1729 /* If and only if all arguments are normal do we need to do an
1730 extended-precision calculation. */
1731 if (category == fcNormal &&
1732 multiplicand.category == fcNormal &&
1733 addend.category == fcNormal) {
1734 lostFraction lost_fraction;
1736 lost_fraction = multiplySignificand(multiplicand, &addend);
1737 fs = normalize(rounding_mode, lost_fraction);
1738 if (lost_fraction != lfExactlyZero)
1739 fs = (opStatus) (fs | opInexact);
1741 /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1742 positive zero unless rounding to minus infinity, except that
1743 adding two like-signed zeroes gives that zero. */
1744 if (category == fcZero && sign != addend.sign)
1745 sign = (rounding_mode == rmTowardNegative);
1747 fs = multiplySpecials(multiplicand);
1749 /* FS can only be opOK or opInvalidOp. There is no more work
1750 to do in the latter case. The IEEE-754R standard says it is
1751 implementation-defined in this case whether, if ADDEND is a
1752 quiet NaN, we raise invalid op; this implementation does so.
1754 If we need to do the addition we can do so with normal
1757 fs = addOrSubtract(addend, rounding_mode, false);
1763 /* Comparison requires normalized numbers. */
1765 APFloat::compare(const APFloat &rhs) const
1769 assertArithmeticOK(*semantics);
1770 assert(semantics == rhs.semantics);
1772 switch (convolve(category, rhs.category)) {
1774 llvm_unreachable(0);
1776 case convolve(fcNaN, fcZero):
1777 case convolve(fcNaN, fcNormal):
1778 case convolve(fcNaN, fcInfinity):
1779 case convolve(fcNaN, fcNaN):
1780 case convolve(fcZero, fcNaN):
1781 case convolve(fcNormal, fcNaN):
1782 case convolve(fcInfinity, fcNaN):
1783 return cmpUnordered;
1785 case convolve(fcInfinity, fcNormal):
1786 case convolve(fcInfinity, fcZero):
1787 case convolve(fcNormal, fcZero):
1791 return cmpGreaterThan;
1793 case convolve(fcNormal, fcInfinity):
1794 case convolve(fcZero, fcInfinity):
1795 case convolve(fcZero, fcNormal):
1797 return cmpGreaterThan;
1801 case convolve(fcInfinity, fcInfinity):
1802 if (sign == rhs.sign)
1807 return cmpGreaterThan;
1809 case convolve(fcZero, fcZero):
1812 case convolve(fcNormal, fcNormal):
1816 /* Two normal numbers. Do they have the same sign? */
1817 if (sign != rhs.sign) {
1819 result = cmpLessThan;
1821 result = cmpGreaterThan;
1823 /* Compare absolute values; invert result if negative. */
1824 result = compareAbsoluteValue(rhs);
1827 if (result == cmpLessThan)
1828 result = cmpGreaterThan;
1829 else if (result == cmpGreaterThan)
1830 result = cmpLessThan;
1837 /// APFloat::convert - convert a value of one floating point type to another.
1838 /// The return value corresponds to the IEEE754 exceptions. *losesInfo
1839 /// records whether the transformation lost information, i.e. whether
1840 /// converting the result back to the original type will produce the
1841 /// original value (this is almost the same as return value==fsOK, but there
1842 /// are edge cases where this is not so).
1845 APFloat::convert(const fltSemantics &toSemantics,
1846 roundingMode rounding_mode, bool *losesInfo)
1848 lostFraction lostFraction;
1849 unsigned int newPartCount, oldPartCount;
1852 assertArithmeticOK(*semantics);
1853 assertArithmeticOK(toSemantics);
1854 lostFraction = lfExactlyZero;
1855 newPartCount = partCountForBits(toSemantics.precision + 1);
1856 oldPartCount = partCount();
1858 /* Handle storage complications. If our new form is wider,
1859 re-allocate our bit pattern into wider storage. If it is
1860 narrower, we ignore the excess parts, but if narrowing to a
1861 single part we need to free the old storage.
1862 Be careful not to reference significandParts for zeroes
1863 and infinities, since it aborts. */
1864 if (newPartCount > oldPartCount) {
1865 integerPart *newParts;
1866 newParts = new integerPart[newPartCount];
1867 APInt::tcSet(newParts, 0, newPartCount);
1868 if (category==fcNormal || category==fcNaN)
1869 APInt::tcAssign(newParts, significandParts(), oldPartCount);
1871 significand.parts = newParts;
1872 } else if (newPartCount < oldPartCount) {
1873 /* Capture any lost fraction through truncation of parts so we get
1874 correct rounding whilst normalizing. */
1875 if (category==fcNormal)
1876 lostFraction = lostFractionThroughTruncation
1877 (significandParts(), oldPartCount, toSemantics.precision);
1878 if (newPartCount == 1) {
1879 integerPart newPart = 0;
1880 if (category==fcNormal || category==fcNaN)
1881 newPart = significandParts()[0];
1883 significand.part = newPart;
1887 if (category == fcNormal) {
1888 /* Re-interpret our bit-pattern. */
1889 exponent += toSemantics.precision - semantics->precision;
1890 semantics = &toSemantics;
1891 fs = normalize(rounding_mode, lostFraction);
1892 *losesInfo = (fs != opOK);
1893 } else if (category == fcNaN) {
1894 int shift = toSemantics.precision - semantics->precision;
1895 // Do this now so significandParts gets the right answer
1896 const fltSemantics *oldSemantics = semantics;
1897 semantics = &toSemantics;
1899 // No normalization here, just truncate
1901 APInt::tcShiftLeft(significandParts(), newPartCount, shift);
1902 else if (shift < 0) {
1903 unsigned ushift = -shift;
1904 // Figure out if we are losing information. This happens
1905 // if are shifting out something other than 0s, or if the x87 long
1906 // double input did not have its integer bit set (pseudo-NaN), or if the
1907 // x87 long double input did not have its QNan bit set (because the x87
1908 // hardware sets this bit when converting a lower-precision NaN to
1909 // x87 long double).
1910 if (APInt::tcLSB(significandParts(), newPartCount) < ushift)
1912 if (oldSemantics == &APFloat::x87DoubleExtended &&
1913 (!(*significandParts() & 0x8000000000000000ULL) ||
1914 !(*significandParts() & 0x4000000000000000ULL)))
1916 APInt::tcShiftRight(significandParts(), newPartCount, ushift);
1918 // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
1919 // does not give you back the same bits. This is dubious, and we
1920 // don't currently do it. You're really supposed to get
1921 // an invalid operation signal at runtime, but nobody does that.
1924 semantics = &toSemantics;
1932 /* Convert a floating point number to an integer according to the
1933 rounding mode. If the rounded integer value is out of range this
1934 returns an invalid operation exception and the contents of the
1935 destination parts are unspecified. If the rounded value is in
1936 range but the floating point number is not the exact integer, the C
1937 standard doesn't require an inexact exception to be raised. IEEE
1938 854 does require it so we do that.
1940 Note that for conversions to integer type the C standard requires
1941 round-to-zero to always be used. */
1943 APFloat::convertToSignExtendedInteger(integerPart *parts, unsigned int width,
1945 roundingMode rounding_mode,
1946 bool *isExact) const
1948 lostFraction lost_fraction;
1949 const integerPart *src;
1950 unsigned int dstPartsCount, truncatedBits;
1952 assertArithmeticOK(*semantics);
1956 /* Handle the three special cases first. */
1957 if (category == fcInfinity || category == fcNaN)
1960 dstPartsCount = partCountForBits(width);
1962 if (category == fcZero) {
1963 APInt::tcSet(parts, 0, dstPartsCount);
1964 // Negative zero can't be represented as an int.
1969 src = significandParts();
1971 /* Step 1: place our absolute value, with any fraction truncated, in
1974 /* Our absolute value is less than one; truncate everything. */
1975 APInt::tcSet(parts, 0, dstPartsCount);
1976 /* For exponent -1 the integer bit represents .5, look at that.
1977 For smaller exponents leftmost truncated bit is 0. */
1978 truncatedBits = semantics->precision -1U - exponent;
1980 /* We want the most significant (exponent + 1) bits; the rest are
1982 unsigned int bits = exponent + 1U;
1984 /* Hopelessly large in magnitude? */
1988 if (bits < semantics->precision) {
1989 /* We truncate (semantics->precision - bits) bits. */
1990 truncatedBits = semantics->precision - bits;
1991 APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits);
1993 /* We want at least as many bits as are available. */
1994 APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0);
1995 APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision);
2000 /* Step 2: work out any lost fraction, and increment the absolute
2001 value if we would round away from zero. */
2002 if (truncatedBits) {
2003 lost_fraction = lostFractionThroughTruncation(src, partCount(),
2005 if (lost_fraction != lfExactlyZero &&
2006 roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) {
2007 if (APInt::tcIncrement(parts, dstPartsCount))
2008 return opInvalidOp; /* Overflow. */
2011 lost_fraction = lfExactlyZero;
2014 /* Step 3: check if we fit in the destination. */
2015 unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1;
2019 /* Negative numbers cannot be represented as unsigned. */
2023 /* It takes omsb bits to represent the unsigned integer value.
2024 We lose a bit for the sign, but care is needed as the
2025 maximally negative integer is a special case. */
2026 if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb)
2029 /* This case can happen because of rounding. */
2034 APInt::tcNegate (parts, dstPartsCount);
2036 if (omsb >= width + !isSigned)
2040 if (lost_fraction == lfExactlyZero) {
2047 /* Same as convertToSignExtendedInteger, except we provide
2048 deterministic values in case of an invalid operation exception,
2049 namely zero for NaNs and the minimal or maximal value respectively
2050 for underflow or overflow.
2051 The *isExact output tells whether the result is exact, in the sense
2052 that converting it back to the original floating point type produces
2053 the original value. This is almost equivalent to result==opOK,
2054 except for negative zeroes.
2057 APFloat::convertToInteger(integerPart *parts, unsigned int width,
2059 roundingMode rounding_mode, bool *isExact) const
2063 fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode,
2066 if (fs == opInvalidOp) {
2067 unsigned int bits, dstPartsCount;
2069 dstPartsCount = partCountForBits(width);
2071 if (category == fcNaN)
2076 bits = width - isSigned;
2078 APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits);
2079 if (sign && isSigned)
2080 APInt::tcShiftLeft(parts, dstPartsCount, width - 1);
2086 /* Convert an unsigned integer SRC to a floating point number,
2087 rounding according to ROUNDING_MODE. The sign of the floating
2088 point number is not modified. */
2090 APFloat::convertFromUnsignedParts(const integerPart *src,
2091 unsigned int srcCount,
2092 roundingMode rounding_mode)
2094 unsigned int omsb, precision, dstCount;
2096 lostFraction lost_fraction;
2098 assertArithmeticOK(*semantics);
2099 category = fcNormal;
2100 omsb = APInt::tcMSB(src, srcCount) + 1;
2101 dst = significandParts();
2102 dstCount = partCount();
2103 precision = semantics->precision;
2105 /* We want the most significant PRECISON bits of SRC. There may not
2106 be that many; extract what we can. */
2107 if (precision <= omsb) {
2108 exponent = omsb - 1;
2109 lost_fraction = lostFractionThroughTruncation(src, srcCount,
2111 APInt::tcExtract(dst, dstCount, src, precision, omsb - precision);
2113 exponent = precision - 1;
2114 lost_fraction = lfExactlyZero;
2115 APInt::tcExtract(dst, dstCount, src, omsb, 0);
2118 return normalize(rounding_mode, lost_fraction);
2122 APFloat::convertFromAPInt(const APInt &Val,
2124 roundingMode rounding_mode)
2126 unsigned int partCount = Val.getNumWords();
2130 if (isSigned && api.isNegative()) {
2135 return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2138 /* Convert a two's complement integer SRC to a floating point number,
2139 rounding according to ROUNDING_MODE. ISSIGNED is true if the
2140 integer is signed, in which case it must be sign-extended. */
2142 APFloat::convertFromSignExtendedInteger(const integerPart *src,
2143 unsigned int srcCount,
2145 roundingMode rounding_mode)
2149 assertArithmeticOK(*semantics);
2151 APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) {
2154 /* If we're signed and negative negate a copy. */
2156 copy = new integerPart[srcCount];
2157 APInt::tcAssign(copy, src, srcCount);
2158 APInt::tcNegate(copy, srcCount);
2159 status = convertFromUnsignedParts(copy, srcCount, rounding_mode);
2163 status = convertFromUnsignedParts(src, srcCount, rounding_mode);
2169 /* FIXME: should this just take a const APInt reference? */
2171 APFloat::convertFromZeroExtendedInteger(const integerPart *parts,
2172 unsigned int width, bool isSigned,
2173 roundingMode rounding_mode)
2175 unsigned int partCount = partCountForBits(width);
2176 APInt api = APInt(width, partCount, parts);
2179 if (isSigned && APInt::tcExtractBit(parts, width - 1)) {
2184 return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2188 APFloat::convertFromHexadecimalString(const StringRef &s,
2189 roundingMode rounding_mode)
2191 lostFraction lost_fraction = lfExactlyZero;
2192 integerPart *significand;
2193 unsigned int bitPos, partsCount;
2194 StringRef::iterator dot, firstSignificantDigit;
2198 category = fcNormal;
2200 significand = significandParts();
2201 partsCount = partCount();
2202 bitPos = partsCount * integerPartWidth;
2204 /* Skip leading zeroes and any (hexa)decimal point. */
2205 StringRef::iterator begin = s.begin();
2206 StringRef::iterator end = s.end();
2207 StringRef::iterator p = skipLeadingZeroesAndAnyDot(begin, end, &dot);
2208 firstSignificantDigit = p;
2211 integerPart hex_value;
2214 assert(dot == end && "String contains multiple dots");
2221 hex_value = hexDigitValue(*p);
2222 if (hex_value == -1U) {
2231 /* Store the number whilst 4-bit nibbles remain. */
2234 hex_value <<= bitPos % integerPartWidth;
2235 significand[bitPos / integerPartWidth] |= hex_value;
2237 lost_fraction = trailingHexadecimalFraction(p, end, hex_value);
2238 while (p != end && hexDigitValue(*p) != -1U)
2245 /* Hex floats require an exponent but not a hexadecimal point. */
2246 assert(p != end && "Hex strings require an exponent");
2247 assert((*p == 'p' || *p == 'P') && "Invalid character in significand");
2248 assert(p != begin && "Significand has no digits");
2249 assert((dot == end || p - begin != 1) && "Significand has no digits");
2251 /* Ignore the exponent if we are zero. */
2252 if (p != firstSignificantDigit) {
2255 /* Implicit hexadecimal point? */
2259 /* Calculate the exponent adjustment implicit in the number of
2260 significant digits. */
2261 expAdjustment = static_cast<int>(dot - firstSignificantDigit);
2262 if (expAdjustment < 0)
2264 expAdjustment = expAdjustment * 4 - 1;
2266 /* Adjust for writing the significand starting at the most
2267 significant nibble. */
2268 expAdjustment += semantics->precision;
2269 expAdjustment -= partsCount * integerPartWidth;
2271 /* Adjust for the given exponent. */
2272 exponent = totalExponent(p + 1, end, expAdjustment);
2275 return normalize(rounding_mode, lost_fraction);
2279 APFloat::roundSignificandWithExponent(const integerPart *decSigParts,
2280 unsigned sigPartCount, int exp,
2281 roundingMode rounding_mode)
2283 unsigned int parts, pow5PartCount;
2284 fltSemantics calcSemantics = { 32767, -32767, 0, true };
2285 integerPart pow5Parts[maxPowerOfFiveParts];
2288 isNearest = (rounding_mode == rmNearestTiesToEven ||
2289 rounding_mode == rmNearestTiesToAway);
2291 parts = partCountForBits(semantics->precision + 11);
2293 /* Calculate pow(5, abs(exp)). */
2294 pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp);
2296 for (;; parts *= 2) {
2297 opStatus sigStatus, powStatus;
2298 unsigned int excessPrecision, truncatedBits;
2300 calcSemantics.precision = parts * integerPartWidth - 1;
2301 excessPrecision = calcSemantics.precision - semantics->precision;
2302 truncatedBits = excessPrecision;
2304 APFloat decSig(calcSemantics, fcZero, sign);
2305 APFloat pow5(calcSemantics, fcZero, false);
2307 sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount,
2308 rmNearestTiesToEven);
2309 powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount,
2310 rmNearestTiesToEven);
2311 /* Add exp, as 10^n = 5^n * 2^n. */
2312 decSig.exponent += exp;
2314 lostFraction calcLostFraction;
2315 integerPart HUerr, HUdistance;
2316 unsigned int powHUerr;
2319 /* multiplySignificand leaves the precision-th bit set to 1. */
2320 calcLostFraction = decSig.multiplySignificand(pow5, NULL);
2321 powHUerr = powStatus != opOK;
2323 calcLostFraction = decSig.divideSignificand(pow5);
2324 /* Denormal numbers have less precision. */
2325 if (decSig.exponent < semantics->minExponent) {
2326 excessPrecision += (semantics->minExponent - decSig.exponent);
2327 truncatedBits = excessPrecision;
2328 if (excessPrecision > calcSemantics.precision)
2329 excessPrecision = calcSemantics.precision;
2331 /* Extra half-ulp lost in reciprocal of exponent. */
2332 powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2;
2335 /* Both multiplySignificand and divideSignificand return the
2336 result with the integer bit set. */
2337 assert(APInt::tcExtractBit
2338 (decSig.significandParts(), calcSemantics.precision - 1) == 1);
2340 HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK,
2342 HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(),
2343 excessPrecision, isNearest);
2345 /* Are we guaranteed to round correctly if we truncate? */
2346 if (HUdistance >= HUerr) {
2347 APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(),
2348 calcSemantics.precision - excessPrecision,
2350 /* Take the exponent of decSig. If we tcExtract-ed less bits
2351 above we must adjust our exponent to compensate for the
2352 implicit right shift. */
2353 exponent = (decSig.exponent + semantics->precision
2354 - (calcSemantics.precision - excessPrecision));
2355 calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(),
2358 return normalize(rounding_mode, calcLostFraction);
2364 APFloat::convertFromDecimalString(const StringRef &str, roundingMode rounding_mode)
2369 /* Scan the text. */
2370 StringRef::iterator p = str.begin();
2371 interpretDecimal(p, str.end(), &D);
2373 /* Handle the quick cases. First the case of no significant digits,
2374 i.e. zero, and then exponents that are obviously too large or too
2375 small. Writing L for log 10 / log 2, a number d.ddddd*10^exp
2376 definitely overflows if
2378 (exp - 1) * L >= maxExponent
2380 and definitely underflows to zero where
2382 (exp + 1) * L <= minExponent - precision
2384 With integer arithmetic the tightest bounds for L are
2386 93/28 < L < 196/59 [ numerator <= 256 ]
2387 42039/12655 < L < 28738/8651 [ numerator <= 65536 ]
2390 if (decDigitValue(*D.firstSigDigit) >= 10U) {
2394 /* Check whether the normalized exponent is high enough to overflow
2395 max during the log-rebasing in the max-exponent check below. */
2396 } else if (D.normalizedExponent - 1 > INT_MAX / 42039) {
2397 fs = handleOverflow(rounding_mode);
2399 /* If it wasn't, then it also wasn't high enough to overflow max
2400 during the log-rebasing in the min-exponent check. Check that it
2401 won't overflow min in either check, then perform the min-exponent
2403 } else if (D.normalizedExponent - 1 < INT_MIN / 42039 ||
2404 (D.normalizedExponent + 1) * 28738 <=
2405 8651 * (semantics->minExponent - (int) semantics->precision)) {
2406 /* Underflow to zero and round. */
2408 fs = normalize(rounding_mode, lfLessThanHalf);
2410 /* We can finally safely perform the max-exponent check. */
2411 } else if ((D.normalizedExponent - 1) * 42039
2412 >= 12655 * semantics->maxExponent) {
2413 /* Overflow and round. */
2414 fs = handleOverflow(rounding_mode);
2416 integerPart *decSignificand;
2417 unsigned int partCount;
2419 /* A tight upper bound on number of bits required to hold an
2420 N-digit decimal integer is N * 196 / 59. Allocate enough space
2421 to hold the full significand, and an extra part required by
2423 partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1;
2424 partCount = partCountForBits(1 + 196 * partCount / 59);
2425 decSignificand = new integerPart[partCount + 1];
2428 /* Convert to binary efficiently - we do almost all multiplication
2429 in an integerPart. When this would overflow do we do a single
2430 bignum multiplication, and then revert again to multiplication
2431 in an integerPart. */
2433 integerPart decValue, val, multiplier;
2441 if (p == str.end()) {
2445 decValue = decDigitValue(*p++);
2446 assert(decValue < 10U && "Invalid character in significand");
2448 val = val * 10 + decValue;
2449 /* The maximum number that can be multiplied by ten with any
2450 digit added without overflowing an integerPart. */
2451 } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10);
2453 /* Multiply out the current part. */
2454 APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val,
2455 partCount, partCount + 1, false);
2457 /* If we used another part (likely but not guaranteed), increase
2459 if (decSignificand[partCount])
2461 } while (p <= D.lastSigDigit);
2463 category = fcNormal;
2464 fs = roundSignificandWithExponent(decSignificand, partCount,
2465 D.exponent, rounding_mode);
2467 delete [] decSignificand;
2474 APFloat::convertFromString(const StringRef &str, roundingMode rounding_mode)
2476 assertArithmeticOK(*semantics);
2477 assert(!str.empty() && "Invalid string length");
2479 /* Handle a leading minus sign. */
2480 StringRef::iterator p = str.begin();
2481 size_t slen = str.size();
2482 sign = *p == '-' ? 1 : 0;
2483 if (*p == '-' || *p == '+') {
2486 assert(slen && "String has no digits");
2489 if (slen >= 2 && p[0] == '0' && (p[1] == 'x' || p[1] == 'X')) {
2490 assert(slen - 2 && "Invalid string");
2491 return convertFromHexadecimalString(StringRef(p + 2, slen - 2),
2495 return convertFromDecimalString(StringRef(p, slen), rounding_mode);
2498 /* Write out a hexadecimal representation of the floating point value
2499 to DST, which must be of sufficient size, in the C99 form
2500 [-]0xh.hhhhp[+-]d. Return the number of characters written,
2501 excluding the terminating NUL.
2503 If UPPERCASE, the output is in upper case, otherwise in lower case.
2505 HEXDIGITS digits appear altogether, rounding the value if
2506 necessary. If HEXDIGITS is 0, the minimal precision to display the
2507 number precisely is used instead. If nothing would appear after
2508 the decimal point it is suppressed.
2510 The decimal exponent is always printed and has at least one digit.
2511 Zero values display an exponent of zero. Infinities and NaNs
2512 appear as "infinity" or "nan" respectively.
2514 The above rules are as specified by C99. There is ambiguity about
2515 what the leading hexadecimal digit should be. This implementation
2516 uses whatever is necessary so that the exponent is displayed as
2517 stored. This implies the exponent will fall within the IEEE format
2518 range, and the leading hexadecimal digit will be 0 (for denormals),
2519 1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with
2520 any other digits zero).
2523 APFloat::convertToHexString(char *dst, unsigned int hexDigits,
2524 bool upperCase, roundingMode rounding_mode) const
2528 assertArithmeticOK(*semantics);
2536 memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1);
2537 dst += sizeof infinityL - 1;
2541 memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1);
2542 dst += sizeof NaNU - 1;
2547 *dst++ = upperCase ? 'X': 'x';
2549 if (hexDigits > 1) {
2551 memset (dst, '0', hexDigits - 1);
2552 dst += hexDigits - 1;
2554 *dst++ = upperCase ? 'P': 'p';
2559 dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode);
2565 return static_cast<unsigned int>(dst - p);
2568 /* Does the hard work of outputting the correctly rounded hexadecimal
2569 form of a normal floating point number with the specified number of
2570 hexadecimal digits. If HEXDIGITS is zero the minimum number of
2571 digits necessary to print the value precisely is output. */
2573 APFloat::convertNormalToHexString(char *dst, unsigned int hexDigits,
2575 roundingMode rounding_mode) const
2577 unsigned int count, valueBits, shift, partsCount, outputDigits;
2578 const char *hexDigitChars;
2579 const integerPart *significand;
2584 *dst++ = upperCase ? 'X': 'x';
2587 hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower;
2589 significand = significandParts();
2590 partsCount = partCount();
2592 /* +3 because the first digit only uses the single integer bit, so
2593 we have 3 virtual zero most-significant-bits. */
2594 valueBits = semantics->precision + 3;
2595 shift = integerPartWidth - valueBits % integerPartWidth;
2597 /* The natural number of digits required ignoring trailing
2598 insignificant zeroes. */
2599 outputDigits = (valueBits - significandLSB () + 3) / 4;
2601 /* hexDigits of zero means use the required number for the
2602 precision. Otherwise, see if we are truncating. If we are,
2603 find out if we need to round away from zero. */
2605 if (hexDigits < outputDigits) {
2606 /* We are dropping non-zero bits, so need to check how to round.
2607 "bits" is the number of dropped bits. */
2609 lostFraction fraction;
2611 bits = valueBits - hexDigits * 4;
2612 fraction = lostFractionThroughTruncation (significand, partsCount, bits);
2613 roundUp = roundAwayFromZero(rounding_mode, fraction, bits);
2615 outputDigits = hexDigits;
2618 /* Write the digits consecutively, and start writing in the location
2619 of the hexadecimal point. We move the most significant digit
2620 left and add the hexadecimal point later. */
2623 count = (valueBits + integerPartWidth - 1) / integerPartWidth;
2625 while (outputDigits && count) {
2628 /* Put the most significant integerPartWidth bits in "part". */
2629 if (--count == partsCount)
2630 part = 0; /* An imaginary higher zero part. */
2632 part = significand[count] << shift;
2635 part |= significand[count - 1] >> (integerPartWidth - shift);
2637 /* Convert as much of "part" to hexdigits as we can. */
2638 unsigned int curDigits = integerPartWidth / 4;
2640 if (curDigits > outputDigits)
2641 curDigits = outputDigits;
2642 dst += partAsHex (dst, part, curDigits, hexDigitChars);
2643 outputDigits -= curDigits;
2649 /* Note that hexDigitChars has a trailing '0'. */
2652 *q = hexDigitChars[hexDigitValue (*q) + 1];
2653 } while (*q == '0');
2656 /* Add trailing zeroes. */
2657 memset (dst, '0', outputDigits);
2658 dst += outputDigits;
2661 /* Move the most significant digit to before the point, and if there
2662 is something after the decimal point add it. This must come
2663 after rounding above. */
2670 /* Finally output the exponent. */
2671 *dst++ = upperCase ? 'P': 'p';
2673 return writeSignedDecimal (dst, exponent);
2676 // For good performance it is desirable for different APFloats
2677 // to produce different integers.
2679 APFloat::getHashValue() const
2681 if (category==fcZero) return sign<<8 | semantics->precision ;
2682 else if (category==fcInfinity) return sign<<9 | semantics->precision;
2683 else if (category==fcNaN) return 1<<10 | semantics->precision;
2685 uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
2686 const integerPart* p = significandParts();
2687 for (int i=partCount(); i>0; i--, p++)
2688 hash ^= ((uint32_t)*p) ^ (uint32_t)((*p)>>32);
2693 // Conversion from APFloat to/from host float/double. It may eventually be
2694 // possible to eliminate these and have everybody deal with APFloats, but that
2695 // will take a while. This approach will not easily extend to long double.
2696 // Current implementation requires integerPartWidth==64, which is correct at
2697 // the moment but could be made more general.
2699 // Denormals have exponent minExponent in APFloat, but minExponent-1 in
2700 // the actual IEEE respresentations. We compensate for that here.
2703 APFloat::convertF80LongDoubleAPFloatToAPInt() const
2705 assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended);
2706 assert(partCount()==2);
2708 uint64_t myexponent, mysignificand;
2710 if (category==fcNormal) {
2711 myexponent = exponent+16383; //bias
2712 mysignificand = significandParts()[0];
2713 if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
2714 myexponent = 0; // denormal
2715 } else if (category==fcZero) {
2718 } else if (category==fcInfinity) {
2719 myexponent = 0x7fff;
2720 mysignificand = 0x8000000000000000ULL;
2722 assert(category == fcNaN && "Unknown category");
2723 myexponent = 0x7fff;
2724 mysignificand = significandParts()[0];
2728 words[0] = mysignificand;
2729 words[1] = ((uint64_t)(sign & 1) << 15) |
2730 (myexponent & 0x7fffLL);
2731 return APInt(80, 2, words);
2735 APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const
2737 assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble);
2738 assert(partCount()==2);
2740 uint64_t myexponent, mysignificand, myexponent2, mysignificand2;
2742 if (category==fcNormal) {
2743 myexponent = exponent + 1023; //bias
2744 myexponent2 = exponent2 + 1023;
2745 mysignificand = significandParts()[0];
2746 mysignificand2 = significandParts()[1];
2747 if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2748 myexponent = 0; // denormal
2749 if (myexponent2==1 && !(mysignificand2 & 0x10000000000000LL))
2750 myexponent2 = 0; // denormal
2751 } else if (category==fcZero) {
2756 } else if (category==fcInfinity) {
2762 assert(category == fcNaN && "Unknown category");
2764 mysignificand = significandParts()[0];
2765 myexponent2 = exponent2;
2766 mysignificand2 = significandParts()[1];
2770 words[0] = ((uint64_t)(sign & 1) << 63) |
2771 ((myexponent & 0x7ff) << 52) |
2772 (mysignificand & 0xfffffffffffffLL);
2773 words[1] = ((uint64_t)(sign2 & 1) << 63) |
2774 ((myexponent2 & 0x7ff) << 52) |
2775 (mysignificand2 & 0xfffffffffffffLL);
2776 return APInt(128, 2, words);
2780 APFloat::convertQuadrupleAPFloatToAPInt() const
2782 assert(semantics == (const llvm::fltSemantics*)&IEEEquad);
2783 assert(partCount()==2);
2785 uint64_t myexponent, mysignificand, mysignificand2;
2787 if (category==fcNormal) {
2788 myexponent = exponent+16383; //bias
2789 mysignificand = significandParts()[0];
2790 mysignificand2 = significandParts()[1];
2791 if (myexponent==1 && !(mysignificand2 & 0x1000000000000LL))
2792 myexponent = 0; // denormal
2793 } else if (category==fcZero) {
2795 mysignificand = mysignificand2 = 0;
2796 } else if (category==fcInfinity) {
2797 myexponent = 0x7fff;
2798 mysignificand = mysignificand2 = 0;
2800 assert(category == fcNaN && "Unknown category!");
2801 myexponent = 0x7fff;
2802 mysignificand = significandParts()[0];
2803 mysignificand2 = significandParts()[1];
2807 words[0] = mysignificand;
2808 words[1] = ((uint64_t)(sign & 1) << 63) |
2809 ((myexponent & 0x7fff) << 48) |
2810 (mysignificand2 & 0xffffffffffffLL);
2812 return APInt(128, 2, words);
2816 APFloat::convertDoubleAPFloatToAPInt() const
2818 assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2819 assert(partCount()==1);
2821 uint64_t myexponent, mysignificand;
2823 if (category==fcNormal) {
2824 myexponent = exponent+1023; //bias
2825 mysignificand = *significandParts();
2826 if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2827 myexponent = 0; // denormal
2828 } else if (category==fcZero) {
2831 } else if (category==fcInfinity) {
2835 assert(category == fcNaN && "Unknown category!");
2837 mysignificand = *significandParts();
2840 return APInt(64, ((((uint64_t)(sign & 1) << 63) |
2841 ((myexponent & 0x7ff) << 52) |
2842 (mysignificand & 0xfffffffffffffLL))));
2846 APFloat::convertFloatAPFloatToAPInt() const
2848 assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2849 assert(partCount()==1);
2851 uint32_t myexponent, mysignificand;
2853 if (category==fcNormal) {
2854 myexponent = exponent+127; //bias
2855 mysignificand = (uint32_t)*significandParts();
2856 if (myexponent == 1 && !(mysignificand & 0x800000))
2857 myexponent = 0; // denormal
2858 } else if (category==fcZero) {
2861 } else if (category==fcInfinity) {
2865 assert(category == fcNaN && "Unknown category!");
2867 mysignificand = (uint32_t)*significandParts();
2870 return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
2871 (mysignificand & 0x7fffff)));
2875 APFloat::convertHalfAPFloatToAPInt() const
2877 assert(semantics == (const llvm::fltSemantics*)&IEEEhalf);
2878 assert(partCount()==1);
2880 uint32_t myexponent, mysignificand;
2882 if (category==fcNormal) {
2883 myexponent = exponent+15; //bias
2884 mysignificand = (uint32_t)*significandParts();
2885 if (myexponent == 1 && !(mysignificand & 0x400))
2886 myexponent = 0; // denormal
2887 } else if (category==fcZero) {
2890 } else if (category==fcInfinity) {
2894 assert(category == fcNaN && "Unknown category!");
2896 mysignificand = (uint32_t)*significandParts();
2899 return APInt(16, (((sign&1) << 15) | ((myexponent&0x1f) << 10) |
2900 (mysignificand & 0x3ff)));
2903 // This function creates an APInt that is just a bit map of the floating
2904 // point constant as it would appear in memory. It is not a conversion,
2905 // and treating the result as a normal integer is unlikely to be useful.
2908 APFloat::bitcastToAPInt() const
2910 if (semantics == (const llvm::fltSemantics*)&IEEEhalf)
2911 return convertHalfAPFloatToAPInt();
2913 if (semantics == (const llvm::fltSemantics*)&IEEEsingle)
2914 return convertFloatAPFloatToAPInt();
2916 if (semantics == (const llvm::fltSemantics*)&IEEEdouble)
2917 return convertDoubleAPFloatToAPInt();
2919 if (semantics == (const llvm::fltSemantics*)&IEEEquad)
2920 return convertQuadrupleAPFloatToAPInt();
2922 if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble)
2923 return convertPPCDoubleDoubleAPFloatToAPInt();
2925 assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended &&
2927 return convertF80LongDoubleAPFloatToAPInt();
2931 APFloat::convertToFloat() const
2933 assert(semantics == (const llvm::fltSemantics*)&IEEEsingle &&
2934 "Float semantics are not IEEEsingle");
2935 APInt api = bitcastToAPInt();
2936 return api.bitsToFloat();
2940 APFloat::convertToDouble() const
2942 assert(semantics == (const llvm::fltSemantics*)&IEEEdouble &&
2943 "Float semantics are not IEEEdouble");
2944 APInt api = bitcastToAPInt();
2945 return api.bitsToDouble();
2948 /// Integer bit is explicit in this format. Intel hardware (387 and later)
2949 /// does not support these bit patterns:
2950 /// exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
2951 /// exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
2952 /// exponent = 0, integer bit 1 ("pseudodenormal")
2953 /// exponent!=0 nor all 1's, integer bit 0 ("unnormal")
2954 /// At the moment, the first two are treated as NaNs, the second two as Normal.
2956 APFloat::initFromF80LongDoubleAPInt(const APInt &api)
2958 assert(api.getBitWidth()==80);
2959 uint64_t i1 = api.getRawData()[0];
2960 uint64_t i2 = api.getRawData()[1];
2961 uint64_t myexponent = (i2 & 0x7fff);
2962 uint64_t mysignificand = i1;
2964 initialize(&APFloat::x87DoubleExtended);
2965 assert(partCount()==2);
2967 sign = static_cast<unsigned int>(i2>>15);
2968 if (myexponent==0 && mysignificand==0) {
2969 // exponent, significand meaningless
2971 } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
2972 // exponent, significand meaningless
2973 category = fcInfinity;
2974 } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
2975 // exponent meaningless
2977 significandParts()[0] = mysignificand;
2978 significandParts()[1] = 0;
2980 category = fcNormal;
2981 exponent = myexponent - 16383;
2982 significandParts()[0] = mysignificand;
2983 significandParts()[1] = 0;
2984 if (myexponent==0) // denormal
2990 APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api)
2992 assert(api.getBitWidth()==128);
2993 uint64_t i1 = api.getRawData()[0];
2994 uint64_t i2 = api.getRawData()[1];
2995 uint64_t myexponent = (i1 >> 52) & 0x7ff;
2996 uint64_t mysignificand = i1 & 0xfffffffffffffLL;
2997 uint64_t myexponent2 = (i2 >> 52) & 0x7ff;
2998 uint64_t mysignificand2 = i2 & 0xfffffffffffffLL;
3000 initialize(&APFloat::PPCDoubleDouble);
3001 assert(partCount()==2);
3003 sign = static_cast<unsigned int>(i1>>63);
3004 sign2 = static_cast<unsigned int>(i2>>63);
3005 if (myexponent==0 && mysignificand==0) {
3006 // exponent, significand meaningless
3007 // exponent2 and significand2 are required to be 0; we don't check
3009 } else if (myexponent==0x7ff && mysignificand==0) {
3010 // exponent, significand meaningless
3011 // exponent2 and significand2 are required to be 0; we don't check
3012 category = fcInfinity;
3013 } else if (myexponent==0x7ff && mysignificand!=0) {
3014 // exponent meaningless. So is the whole second word, but keep it
3017 exponent2 = myexponent2;
3018 significandParts()[0] = mysignificand;
3019 significandParts()[1] = mysignificand2;
3021 category = fcNormal;
3022 // Note there is no category2; the second word is treated as if it is
3023 // fcNormal, although it might be something else considered by itself.
3024 exponent = myexponent - 1023;
3025 exponent2 = myexponent2 - 1023;
3026 significandParts()[0] = mysignificand;
3027 significandParts()[1] = mysignificand2;
3028 if (myexponent==0) // denormal
3031 significandParts()[0] |= 0x10000000000000LL; // integer bit
3035 significandParts()[1] |= 0x10000000000000LL; // integer bit
3040 APFloat::initFromQuadrupleAPInt(const APInt &api)
3042 assert(api.getBitWidth()==128);
3043 uint64_t i1 = api.getRawData()[0];
3044 uint64_t i2 = api.getRawData()[1];
3045 uint64_t myexponent = (i2 >> 48) & 0x7fff;
3046 uint64_t mysignificand = i1;
3047 uint64_t mysignificand2 = i2 & 0xffffffffffffLL;
3049 initialize(&APFloat::IEEEquad);
3050 assert(partCount()==2);
3052 sign = static_cast<unsigned int>(i2>>63);
3053 if (myexponent==0 &&
3054 (mysignificand==0 && mysignificand2==0)) {
3055 // exponent, significand meaningless
3057 } else if (myexponent==0x7fff &&
3058 (mysignificand==0 && mysignificand2==0)) {
3059 // exponent, significand meaningless
3060 category = fcInfinity;
3061 } else if (myexponent==0x7fff &&
3062 (mysignificand!=0 || mysignificand2 !=0)) {
3063 // exponent meaningless
3065 significandParts()[0] = mysignificand;
3066 significandParts()[1] = mysignificand2;
3068 category = fcNormal;
3069 exponent = myexponent - 16383;
3070 significandParts()[0] = mysignificand;
3071 significandParts()[1] = mysignificand2;
3072 if (myexponent==0) // denormal
3075 significandParts()[1] |= 0x1000000000000LL; // integer bit
3080 APFloat::initFromDoubleAPInt(const APInt &api)
3082 assert(api.getBitWidth()==64);
3083 uint64_t i = *api.getRawData();
3084 uint64_t myexponent = (i >> 52) & 0x7ff;
3085 uint64_t mysignificand = i & 0xfffffffffffffLL;
3087 initialize(&APFloat::IEEEdouble);
3088 assert(partCount()==1);
3090 sign = static_cast<unsigned int>(i>>63);
3091 if (myexponent==0 && mysignificand==0) {
3092 // exponent, significand meaningless
3094 } else if (myexponent==0x7ff && mysignificand==0) {
3095 // exponent, significand meaningless
3096 category = fcInfinity;
3097 } else if (myexponent==0x7ff && mysignificand!=0) {
3098 // exponent meaningless
3100 *significandParts() = mysignificand;
3102 category = fcNormal;
3103 exponent = myexponent - 1023;
3104 *significandParts() = mysignificand;
3105 if (myexponent==0) // denormal
3108 *significandParts() |= 0x10000000000000LL; // integer bit
3113 APFloat::initFromFloatAPInt(const APInt & api)
3115 assert(api.getBitWidth()==32);
3116 uint32_t i = (uint32_t)*api.getRawData();
3117 uint32_t myexponent = (i >> 23) & 0xff;
3118 uint32_t mysignificand = i & 0x7fffff;
3120 initialize(&APFloat::IEEEsingle);
3121 assert(partCount()==1);
3124 if (myexponent==0 && mysignificand==0) {
3125 // exponent, significand meaningless
3127 } else if (myexponent==0xff && mysignificand==0) {
3128 // exponent, significand meaningless
3129 category = fcInfinity;
3130 } else if (myexponent==0xff && mysignificand!=0) {
3131 // sign, exponent, significand meaningless
3133 *significandParts() = mysignificand;
3135 category = fcNormal;
3136 exponent = myexponent - 127; //bias
3137 *significandParts() = mysignificand;
3138 if (myexponent==0) // denormal
3141 *significandParts() |= 0x800000; // integer bit
3146 APFloat::initFromHalfAPInt(const APInt & api)
3148 assert(api.getBitWidth()==16);
3149 uint32_t i = (uint32_t)*api.getRawData();
3150 uint32_t myexponent = (i >> 10) & 0x1f;
3151 uint32_t mysignificand = i & 0x3ff;
3153 initialize(&APFloat::IEEEhalf);
3154 assert(partCount()==1);
3157 if (myexponent==0 && mysignificand==0) {
3158 // exponent, significand meaningless
3160 } else if (myexponent==0x1f && mysignificand==0) {
3161 // exponent, significand meaningless
3162 category = fcInfinity;
3163 } else if (myexponent==0x1f && mysignificand!=0) {
3164 // sign, exponent, significand meaningless
3166 *significandParts() = mysignificand;
3168 category = fcNormal;
3169 exponent = myexponent - 15; //bias
3170 *significandParts() = mysignificand;
3171 if (myexponent==0) // denormal
3174 *significandParts() |= 0x400; // integer bit
3178 /// Treat api as containing the bits of a floating point number. Currently
3179 /// we infer the floating point type from the size of the APInt. The
3180 /// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful
3181 /// when the size is anything else).
3183 APFloat::initFromAPInt(const APInt& api, bool isIEEE)
3185 if (api.getBitWidth() == 16)
3186 return initFromHalfAPInt(api);
3187 else if (api.getBitWidth() == 32)
3188 return initFromFloatAPInt(api);
3189 else if (api.getBitWidth()==64)
3190 return initFromDoubleAPInt(api);
3191 else if (api.getBitWidth()==80)
3192 return initFromF80LongDoubleAPInt(api);
3193 else if (api.getBitWidth()==128)
3195 initFromQuadrupleAPInt(api) : initFromPPCDoubleDoubleAPInt(api));
3197 llvm_unreachable(0);
3200 APFloat APFloat::getLargest(const fltSemantics &Sem, bool Negative) {
3201 APFloat Val(Sem, fcNormal, Negative);
3203 // We want (in interchange format):
3204 // sign = {Negative}
3206 // significand = 1..1
3208 Val.exponent = Sem.maxExponent; // unbiased
3210 // 1-initialize all bits....
3211 Val.zeroSignificand();
3212 integerPart *significand = Val.significandParts();
3213 unsigned N = partCountForBits(Sem.precision);
3214 for (unsigned i = 0; i != N; ++i)
3215 significand[i] = ~((integerPart) 0);
3217 // ...and then clear the top bits for internal consistency.
3219 (((integerPart) 1) << ((Sem.precision % integerPartWidth) - 1)) - 1;
3224 APFloat APFloat::getSmallest(const fltSemantics &Sem, bool Negative) {
3225 APFloat Val(Sem, fcNormal, Negative);
3227 // We want (in interchange format):
3228 // sign = {Negative}
3230 // significand = 0..01
3232 Val.exponent = Sem.minExponent; // unbiased
3233 Val.zeroSignificand();
3234 Val.significandParts()[0] = 1;
3238 APFloat APFloat::getSmallestNormalized(const fltSemantics &Sem, bool Negative) {
3239 APFloat Val(Sem, fcNormal, Negative);
3241 // We want (in interchange format):
3242 // sign = {Negative}
3244 // significand = 10..0
3246 Val.exponent = Sem.minExponent;
3247 Val.zeroSignificand();
3248 Val.significandParts()[partCountForBits(Sem.precision)-1] |=
3249 (((integerPart) 1) << ((Sem.precision % integerPartWidth) - 1));
3254 APFloat::APFloat(const APInt& api, bool isIEEE)
3256 initFromAPInt(api, isIEEE);
3259 APFloat::APFloat(float f)
3261 APInt api = APInt(32, 0);
3262 initFromAPInt(api.floatToBits(f));
3265 APFloat::APFloat(double d)
3267 APInt api = APInt(64, 0);
3268 initFromAPInt(api.doubleToBits(d));
3272 static void append(SmallVectorImpl<char> &Buffer,
3273 unsigned N, const char *Str) {
3274 unsigned Start = Buffer.size();
3275 Buffer.set_size(Start + N);
3276 memcpy(&Buffer[Start], Str, N);
3279 template <unsigned N>
3280 void append(SmallVectorImpl<char> &Buffer, const char (&Str)[N]) {
3281 append(Buffer, N, Str);
3284 /// Removes data from the given significand until it is no more
3285 /// precise than is required for the desired precision.
3286 void AdjustToPrecision(APInt &significand,
3287 int &exp, unsigned FormatPrecision) {
3288 unsigned bits = significand.getActiveBits();
3290 // 196/59 is a very slight overestimate of lg_2(10).
3291 unsigned bitsRequired = (FormatPrecision * 196 + 58) / 59;
3293 if (bits <= bitsRequired) return;
3295 unsigned tensRemovable = (bits - bitsRequired) * 59 / 196;
3296 if (!tensRemovable) return;
3298 exp += tensRemovable;
3300 APInt divisor(significand.getBitWidth(), 1);
3301 APInt powten(significand.getBitWidth(), 10);
3303 if (tensRemovable & 1)
3305 tensRemovable >>= 1;
3306 if (!tensRemovable) break;
3310 significand = significand.udiv(divisor);
3312 // Truncate the significand down to its active bit count, but
3313 // don't try to drop below 32.
3314 unsigned newPrecision = std::max(32U, significand.getActiveBits());
3315 significand.trunc(newPrecision);
3319 void AdjustToPrecision(SmallVectorImpl<char> &buffer,
3320 int &exp, unsigned FormatPrecision) {
3321 unsigned N = buffer.size();
3322 if (N <= FormatPrecision) return;
3324 // The most significant figures are the last ones in the buffer.
3325 unsigned FirstSignificant = N - FormatPrecision;
3328 // FIXME: this probably shouldn't use 'round half up'.
3330 // Rounding down is just a truncation, except we also want to drop
3331 // trailing zeros from the new result.
3332 if (buffer[FirstSignificant - 1] < '5') {
3333 while (buffer[FirstSignificant] == '0')
3336 exp += FirstSignificant;
3337 buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3341 // Rounding up requires a decimal add-with-carry. If we continue
3342 // the carry, the newly-introduced zeros will just be truncated.
3343 for (unsigned I = FirstSignificant; I != N; ++I) {
3344 if (buffer[I] == '9') {
3352 // If we carried through, we have exactly one digit of precision.
3353 if (FirstSignificant == N) {
3354 exp += FirstSignificant;
3356 buffer.push_back('1');
3360 exp += FirstSignificant;
3361 buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3365 void APFloat::toString(SmallVectorImpl<char> &Str,
3366 unsigned FormatPrecision,
3367 unsigned FormatMaxPadding) const {
3371 return append(Str, "-Inf");
3373 return append(Str, "+Inf");
3375 case fcNaN: return append(Str, "NaN");
3381 if (!FormatMaxPadding)
3382 append(Str, "0.0E+0");
3394 // Decompose the number into an APInt and an exponent.
3395 int exp = exponent - ((int) semantics->precision - 1);
3396 APInt significand(semantics->precision,
3397 partCountForBits(semantics->precision),
3398 significandParts());
3400 // Set FormatPrecision if zero. We want to do this before we
3401 // truncate trailing zeros, as those are part of the precision.
3402 if (!FormatPrecision) {
3403 // It's an interesting question whether to use the nominal
3404 // precision or the active precision here for denormals.
3406 // FormatPrecision = ceil(significandBits / lg_2(10))
3407 FormatPrecision = (semantics->precision * 59 + 195) / 196;
3410 // Ignore trailing binary zeros.
3411 int trailingZeros = significand.countTrailingZeros();
3412 exp += trailingZeros;
3413 significand = significand.lshr(trailingZeros);
3415 // Change the exponent from 2^e to 10^e.
3418 } else if (exp > 0) {
3420 significand.zext(semantics->precision + exp);
3421 significand <<= exp;
3423 } else { /* exp < 0 */
3426 // We transform this using the identity:
3427 // (N)(2^-e) == (N)(5^e)(10^-e)
3428 // This means we have to multiply N (the significand) by 5^e.
3429 // To avoid overflow, we have to operate on numbers large
3430 // enough to store N * 5^e:
3431 // log2(N * 5^e) == log2(N) + e * log2(5)
3432 // <= semantics->precision + e * 137 / 59
3433 // (log_2(5) ~ 2.321928 < 2.322034 ~ 137/59)
3435 unsigned precision = semantics->precision + 137 * texp / 59;
3437 // Multiply significand by 5^e.
3438 // N * 5^0101 == N * 5^(1*1) * 5^(0*2) * 5^(1*4) * 5^(0*8)
3439 significand.zext(precision);
3440 APInt five_to_the_i(precision, 5);
3442 if (texp & 1) significand *= five_to_the_i;
3446 five_to_the_i *= five_to_the_i;
3450 AdjustToPrecision(significand, exp, FormatPrecision);
3452 llvm::SmallVector<char, 256> buffer;
3455 unsigned precision = significand.getBitWidth();
3456 APInt ten(precision, 10);
3457 APInt digit(precision, 0);
3459 bool inTrail = true;
3460 while (significand != 0) {
3461 // digit <- significand % 10
3462 // significand <- significand / 10
3463 APInt::udivrem(significand, ten, significand, digit);
3465 unsigned d = digit.getZExtValue();
3467 // Drop trailing zeros.
3468 if (inTrail && !d) exp++;
3470 buffer.push_back((char) ('0' + d));
3475 assert(!buffer.empty() && "no characters in buffer!");
3477 // Drop down to FormatPrecision.
3478 // TODO: don't do more precise calculations above than are required.
3479 AdjustToPrecision(buffer, exp, FormatPrecision);
3481 unsigned NDigits = buffer.size();
3483 // Check whether we should use scientific notation.
3484 bool FormatScientific;
3485 if (!FormatMaxPadding)
3486 FormatScientific = true;
3491 // But we shouldn't make the number look more precise than it is.
3492 FormatScientific = ((unsigned) exp > FormatMaxPadding ||
3493 NDigits + (unsigned) exp > FormatPrecision);
3495 // Power of the most significant digit.
3496 int MSD = exp + (int) (NDigits - 1);
3499 FormatScientific = false;
3501 // 765e-5 == 0.00765
3503 FormatScientific = ((unsigned) -MSD) > FormatMaxPadding;
3508 // Scientific formatting is pretty straightforward.
3509 if (FormatScientific) {
3510 exp += (NDigits - 1);
3512 Str.push_back(buffer[NDigits-1]);
3517 for (unsigned I = 1; I != NDigits; ++I)
3518 Str.push_back(buffer[NDigits-1-I]);
3521 Str.push_back(exp >= 0 ? '+' : '-');
3522 if (exp < 0) exp = -exp;
3523 SmallVector<char, 6> expbuf;
3525 expbuf.push_back((char) ('0' + (exp % 10)));
3528 for (unsigned I = 0, E = expbuf.size(); I != E; ++I)
3529 Str.push_back(expbuf[E-1-I]);
3533 // Non-scientific, positive exponents.
3535 for (unsigned I = 0; I != NDigits; ++I)
3536 Str.push_back(buffer[NDigits-1-I]);
3537 for (unsigned I = 0; I != (unsigned) exp; ++I)
3542 // Non-scientific, negative exponents.
3544 // The number of digits to the left of the decimal point.
3545 int NWholeDigits = exp + (int) NDigits;
3548 if (NWholeDigits > 0) {
3549 for (; I != (unsigned) NWholeDigits; ++I)
3550 Str.push_back(buffer[NDigits-I-1]);
3553 unsigned NZeros = 1 + (unsigned) -NWholeDigits;
3557 for (unsigned Z = 1; Z != NZeros; ++Z)
3561 for (; I != NDigits; ++I)
3562 Str.push_back(buffer[NDigits-I-1]);