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)
910 APInt::tcShiftLeft(fullSignificand, newPartsCount,
911 extendedPrecision - omsb);
912 exponent -= extendedPrecision - omsb;
915 /* Create new semantics. */
916 extendedSemantics = *semantics;
917 extendedSemantics.precision = extendedPrecision;
919 if(newPartsCount == 1)
920 significand.part = fullSignificand[0];
922 significand.parts = fullSignificand;
923 semantics = &extendedSemantics;
925 APFloat extendedAddend(*addend);
926 status = extendedAddend.convert(extendedSemantics, rmTowardZero, &ignored);
927 assert(status == opOK);
928 lost_fraction = addOrSubtractSignificand(extendedAddend, false);
930 /* Restore our state. */
931 if(newPartsCount == 1)
932 fullSignificand[0] = significand.part;
933 significand = savedSignificand;
934 semantics = savedSemantics;
936 omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
939 exponent -= (precision - 1);
941 if(omsb > precision) {
942 unsigned int bits, significantParts;
945 bits = omsb - precision;
946 significantParts = partCountForBits(omsb);
947 lf = shiftRight(fullSignificand, significantParts, bits);
948 lost_fraction = combineLostFractions(lf, lost_fraction);
952 APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
954 if(newPartsCount > 4)
955 delete [] fullSignificand;
957 return lost_fraction;
960 /* Multiply the significands of LHS and RHS to DST. */
962 APFloat::divideSignificand(const APFloat &rhs)
964 unsigned int bit, i, partsCount;
965 const integerPart *rhsSignificand;
966 integerPart *lhsSignificand, *dividend, *divisor;
967 integerPart scratch[4];
968 lostFraction lost_fraction;
970 assert(semantics == rhs.semantics);
972 lhsSignificand = significandParts();
973 rhsSignificand = rhs.significandParts();
974 partsCount = partCount();
977 dividend = new integerPart[partsCount * 2];
981 divisor = dividend + partsCount;
983 /* Copy the dividend and divisor as they will be modified in-place. */
984 for(i = 0; i < partsCount; i++) {
985 dividend[i] = lhsSignificand[i];
986 divisor[i] = rhsSignificand[i];
987 lhsSignificand[i] = 0;
990 exponent -= rhs.exponent;
992 unsigned int precision = semantics->precision;
994 /* Normalize the divisor. */
995 bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
998 APInt::tcShiftLeft(divisor, partsCount, bit);
1001 /* Normalize the dividend. */
1002 bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
1005 APInt::tcShiftLeft(dividend, partsCount, bit);
1008 /* Ensure the dividend >= divisor initially for the loop below.
1009 Incidentally, this means that the division loop below is
1010 guaranteed to set the integer bit to one. */
1011 if(APInt::tcCompare(dividend, divisor, partsCount) < 0) {
1013 APInt::tcShiftLeft(dividend, partsCount, 1);
1014 assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
1017 /* Long division. */
1018 for(bit = precision; bit; bit -= 1) {
1019 if(APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
1020 APInt::tcSubtract(dividend, divisor, 0, partsCount);
1021 APInt::tcSetBit(lhsSignificand, bit - 1);
1024 APInt::tcShiftLeft(dividend, partsCount, 1);
1027 /* Figure out the lost fraction. */
1028 int cmp = APInt::tcCompare(dividend, divisor, partsCount);
1031 lost_fraction = lfMoreThanHalf;
1033 lost_fraction = lfExactlyHalf;
1034 else if(APInt::tcIsZero(dividend, partsCount))
1035 lost_fraction = lfExactlyZero;
1037 lost_fraction = lfLessThanHalf;
1042 return lost_fraction;
1046 APFloat::significandMSB() const
1048 return APInt::tcMSB(significandParts(), partCount());
1052 APFloat::significandLSB() const
1054 return APInt::tcLSB(significandParts(), partCount());
1057 /* Note that a zero result is NOT normalized to fcZero. */
1059 APFloat::shiftSignificandRight(unsigned int bits)
1061 /* Our exponent should not overflow. */
1062 assert((exponent_t) (exponent + bits) >= exponent);
1066 return shiftRight(significandParts(), partCount(), bits);
1069 /* Shift the significand left BITS bits, subtract BITS from its exponent. */
1071 APFloat::shiftSignificandLeft(unsigned int bits)
1073 assert(bits < semantics->precision);
1076 unsigned int partsCount = partCount();
1078 APInt::tcShiftLeft(significandParts(), partsCount, bits);
1081 assert(!APInt::tcIsZero(significandParts(), partsCount));
1086 APFloat::compareAbsoluteValue(const APFloat &rhs) const
1090 assert(semantics == rhs.semantics);
1091 assert(category == fcNormal);
1092 assert(rhs.category == fcNormal);
1094 compare = exponent - rhs.exponent;
1096 /* If exponents are equal, do an unsigned bignum comparison of the
1099 compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
1103 return cmpGreaterThan;
1104 else if(compare < 0)
1110 /* Handle overflow. Sign is preserved. We either become infinity or
1111 the largest finite number. */
1113 APFloat::handleOverflow(roundingMode rounding_mode)
1116 if(rounding_mode == rmNearestTiesToEven
1117 || rounding_mode == rmNearestTiesToAway
1118 || (rounding_mode == rmTowardPositive && !sign)
1119 || (rounding_mode == rmTowardNegative && sign))
1121 category = fcInfinity;
1122 return (opStatus) (opOverflow | opInexact);
1125 /* Otherwise we become the largest finite number. */
1126 category = fcNormal;
1127 exponent = semantics->maxExponent;
1128 APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
1129 semantics->precision);
1134 /* Returns TRUE if, when truncating the current number, with BIT the
1135 new LSB, with the given lost fraction and rounding mode, the result
1136 would need to be rounded away from zero (i.e., by increasing the
1137 signficand). This routine must work for fcZero of both signs, and
1138 fcNormal numbers. */
1140 APFloat::roundAwayFromZero(roundingMode rounding_mode,
1141 lostFraction lost_fraction,
1142 unsigned int bit) const
1144 /* NaNs and infinities should not have lost fractions. */
1145 assert(category == fcNormal || category == fcZero);
1147 /* Current callers never pass this so we don't handle it. */
1148 assert(lost_fraction != lfExactlyZero);
1150 switch (rounding_mode) {
1152 llvm_unreachable(0);
1154 case rmNearestTiesToAway:
1155 return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
1157 case rmNearestTiesToEven:
1158 if(lost_fraction == lfMoreThanHalf)
1161 /* Our zeroes don't have a significand to test. */
1162 if(lost_fraction == lfExactlyHalf && category != fcZero)
1163 return APInt::tcExtractBit(significandParts(), bit);
1170 case rmTowardPositive:
1171 return sign == false;
1173 case rmTowardNegative:
1174 return sign == true;
1179 APFloat::normalize(roundingMode rounding_mode,
1180 lostFraction lost_fraction)
1182 unsigned int omsb; /* One, not zero, based MSB. */
1185 if(category != fcNormal)
1188 /* Before rounding normalize the exponent of fcNormal numbers. */
1189 omsb = significandMSB() + 1;
1192 /* OMSB is numbered from 1. We want to place it in the integer
1193 bit numbered PRECISON if possible, with a compensating change in
1195 exponentChange = omsb - semantics->precision;
1197 /* If the resulting exponent is too high, overflow according to
1198 the rounding mode. */
1199 if(exponent + exponentChange > semantics->maxExponent)
1200 return handleOverflow(rounding_mode);
1202 /* Subnormal numbers have exponent minExponent, and their MSB
1203 is forced based on that. */
1204 if(exponent + exponentChange < semantics->minExponent)
1205 exponentChange = semantics->minExponent - exponent;
1207 /* Shifting left is easy as we don't lose precision. */
1208 if(exponentChange < 0) {
1209 assert(lost_fraction == lfExactlyZero);
1211 shiftSignificandLeft(-exponentChange);
1216 if(exponentChange > 0) {
1219 /* Shift right and capture any new lost fraction. */
1220 lf = shiftSignificandRight(exponentChange);
1222 lost_fraction = combineLostFractions(lf, lost_fraction);
1224 /* Keep OMSB up-to-date. */
1225 if(omsb > (unsigned) exponentChange)
1226 omsb -= exponentChange;
1232 /* Now round the number according to rounding_mode given the lost
1235 /* As specified in IEEE 754, since we do not trap we do not report
1236 underflow for exact results. */
1237 if(lost_fraction == lfExactlyZero) {
1238 /* Canonicalize zeroes. */
1245 /* Increment the significand if we're rounding away from zero. */
1246 if(roundAwayFromZero(rounding_mode, lost_fraction, 0)) {
1248 exponent = semantics->minExponent;
1250 incrementSignificand();
1251 omsb = significandMSB() + 1;
1253 /* Did the significand increment overflow? */
1254 if(omsb == (unsigned) semantics->precision + 1) {
1255 /* Renormalize by incrementing the exponent and shifting our
1256 significand right one. However if we already have the
1257 maximum exponent we overflow to infinity. */
1258 if(exponent == semantics->maxExponent) {
1259 category = fcInfinity;
1261 return (opStatus) (opOverflow | opInexact);
1264 shiftSignificandRight(1);
1270 /* The normal case - we were and are not denormal, and any
1271 significand increment above didn't overflow. */
1272 if(omsb == semantics->precision)
1275 /* We have a non-zero denormal. */
1276 assert(omsb < semantics->precision);
1278 /* Canonicalize zeroes. */
1282 /* The fcZero case is a denormal that underflowed to zero. */
1283 return (opStatus) (opUnderflow | opInexact);
1287 APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
1289 switch (convolve(category, rhs.category)) {
1291 llvm_unreachable(0);
1293 case convolve(fcNaN, fcZero):
1294 case convolve(fcNaN, fcNormal):
1295 case convolve(fcNaN, fcInfinity):
1296 case convolve(fcNaN, fcNaN):
1297 case convolve(fcNormal, fcZero):
1298 case convolve(fcInfinity, fcNormal):
1299 case convolve(fcInfinity, fcZero):
1302 case convolve(fcZero, fcNaN):
1303 case convolve(fcNormal, fcNaN):
1304 case convolve(fcInfinity, fcNaN):
1306 copySignificand(rhs);
1309 case convolve(fcNormal, fcInfinity):
1310 case convolve(fcZero, fcInfinity):
1311 category = fcInfinity;
1312 sign = rhs.sign ^ subtract;
1315 case convolve(fcZero, fcNormal):
1317 sign = rhs.sign ^ subtract;
1320 case convolve(fcZero, fcZero):
1321 /* Sign depends on rounding mode; handled by caller. */
1324 case convolve(fcInfinity, fcInfinity):
1325 /* Differently signed infinities can only be validly
1327 if(((sign ^ rhs.sign)!=0) != subtract) {
1334 case convolve(fcNormal, fcNormal):
1339 /* Add or subtract two normal numbers. */
1341 APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
1344 lostFraction lost_fraction;
1347 /* Determine if the operation on the absolute values is effectively
1348 an addition or subtraction. */
1349 subtract ^= (sign ^ rhs.sign) ? true : false;
1351 /* Are we bigger exponent-wise than the RHS? */
1352 bits = exponent - rhs.exponent;
1354 /* Subtraction is more subtle than one might naively expect. */
1356 APFloat temp_rhs(rhs);
1360 reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
1361 lost_fraction = lfExactlyZero;
1362 } else if (bits > 0) {
1363 lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
1364 shiftSignificandLeft(1);
1367 lost_fraction = shiftSignificandRight(-bits - 1);
1368 temp_rhs.shiftSignificandLeft(1);
1373 carry = temp_rhs.subtractSignificand
1374 (*this, lost_fraction != lfExactlyZero);
1375 copySignificand(temp_rhs);
1378 carry = subtractSignificand
1379 (temp_rhs, lost_fraction != lfExactlyZero);
1382 /* Invert the lost fraction - it was on the RHS and
1384 if(lost_fraction == lfLessThanHalf)
1385 lost_fraction = lfMoreThanHalf;
1386 else if(lost_fraction == lfMoreThanHalf)
1387 lost_fraction = lfLessThanHalf;
1389 /* The code above is intended to ensure that no borrow is
1394 APFloat temp_rhs(rhs);
1396 lost_fraction = temp_rhs.shiftSignificandRight(bits);
1397 carry = addSignificand(temp_rhs);
1399 lost_fraction = shiftSignificandRight(-bits);
1400 carry = addSignificand(rhs);
1403 /* We have a guard bit; generating a carry cannot happen. */
1407 return lost_fraction;
1411 APFloat::multiplySpecials(const APFloat &rhs)
1413 switch (convolve(category, rhs.category)) {
1415 llvm_unreachable(0);
1417 case convolve(fcNaN, fcZero):
1418 case convolve(fcNaN, fcNormal):
1419 case convolve(fcNaN, fcInfinity):
1420 case convolve(fcNaN, fcNaN):
1423 case convolve(fcZero, fcNaN):
1424 case convolve(fcNormal, fcNaN):
1425 case convolve(fcInfinity, fcNaN):
1427 copySignificand(rhs);
1430 case convolve(fcNormal, fcInfinity):
1431 case convolve(fcInfinity, fcNormal):
1432 case convolve(fcInfinity, fcInfinity):
1433 category = fcInfinity;
1436 case convolve(fcZero, fcNormal):
1437 case convolve(fcNormal, fcZero):
1438 case convolve(fcZero, fcZero):
1442 case convolve(fcZero, fcInfinity):
1443 case convolve(fcInfinity, fcZero):
1447 case convolve(fcNormal, fcNormal):
1453 APFloat::divideSpecials(const APFloat &rhs)
1455 switch (convolve(category, rhs.category)) {
1457 llvm_unreachable(0);
1459 case convolve(fcNaN, fcZero):
1460 case convolve(fcNaN, fcNormal):
1461 case convolve(fcNaN, fcInfinity):
1462 case convolve(fcNaN, fcNaN):
1463 case convolve(fcInfinity, fcZero):
1464 case convolve(fcInfinity, fcNormal):
1465 case convolve(fcZero, fcInfinity):
1466 case convolve(fcZero, fcNormal):
1469 case convolve(fcZero, fcNaN):
1470 case convolve(fcNormal, fcNaN):
1471 case convolve(fcInfinity, fcNaN):
1473 copySignificand(rhs);
1476 case convolve(fcNormal, fcInfinity):
1480 case convolve(fcNormal, fcZero):
1481 category = fcInfinity;
1484 case convolve(fcInfinity, fcInfinity):
1485 case convolve(fcZero, fcZero):
1489 case convolve(fcNormal, fcNormal):
1495 APFloat::modSpecials(const APFloat &rhs)
1497 switch (convolve(category, rhs.category)) {
1499 llvm_unreachable(0);
1501 case convolve(fcNaN, fcZero):
1502 case convolve(fcNaN, fcNormal):
1503 case convolve(fcNaN, fcInfinity):
1504 case convolve(fcNaN, fcNaN):
1505 case convolve(fcZero, fcInfinity):
1506 case convolve(fcZero, fcNormal):
1507 case convolve(fcNormal, fcInfinity):
1510 case convolve(fcZero, fcNaN):
1511 case convolve(fcNormal, fcNaN):
1512 case convolve(fcInfinity, fcNaN):
1514 copySignificand(rhs);
1517 case convolve(fcNormal, fcZero):
1518 case convolve(fcInfinity, fcZero):
1519 case convolve(fcInfinity, fcNormal):
1520 case convolve(fcInfinity, fcInfinity):
1521 case convolve(fcZero, fcZero):
1525 case convolve(fcNormal, fcNormal):
1532 APFloat::changeSign()
1534 /* Look mummy, this one's easy. */
1539 APFloat::clearSign()
1541 /* So is this one. */
1546 APFloat::copySign(const APFloat &rhs)
1552 /* Normalized addition or subtraction. */
1554 APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1559 assertArithmeticOK(*semantics);
1561 fs = addOrSubtractSpecials(rhs, subtract);
1563 /* This return code means it was not a simple case. */
1564 if(fs == opDivByZero) {
1565 lostFraction lost_fraction;
1567 lost_fraction = addOrSubtractSignificand(rhs, subtract);
1568 fs = normalize(rounding_mode, lost_fraction);
1570 /* Can only be zero if we lost no fraction. */
1571 assert(category != fcZero || lost_fraction == lfExactlyZero);
1574 /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1575 positive zero unless rounding to minus infinity, except that
1576 adding two like-signed zeroes gives that zero. */
1577 if(category == fcZero) {
1578 if(rhs.category != fcZero || (sign == rhs.sign) == subtract)
1579 sign = (rounding_mode == rmTowardNegative);
1585 /* Normalized addition. */
1587 APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1589 return addOrSubtract(rhs, rounding_mode, false);
1592 /* Normalized subtraction. */
1594 APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1596 return addOrSubtract(rhs, rounding_mode, true);
1599 /* Normalized multiply. */
1601 APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1605 assertArithmeticOK(*semantics);
1607 fs = multiplySpecials(rhs);
1609 if(category == fcNormal) {
1610 lostFraction lost_fraction = multiplySignificand(rhs, 0);
1611 fs = normalize(rounding_mode, lost_fraction);
1612 if(lost_fraction != lfExactlyZero)
1613 fs = (opStatus) (fs | opInexact);
1619 /* Normalized divide. */
1621 APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1625 assertArithmeticOK(*semantics);
1627 fs = divideSpecials(rhs);
1629 if(category == fcNormal) {
1630 lostFraction lost_fraction = divideSignificand(rhs);
1631 fs = normalize(rounding_mode, lost_fraction);
1632 if(lost_fraction != lfExactlyZero)
1633 fs = (opStatus) (fs | opInexact);
1639 /* Normalized remainder. This is not currently correct in all cases. */
1641 APFloat::remainder(const APFloat &rhs)
1645 unsigned int origSign = sign;
1647 assertArithmeticOK(*semantics);
1648 fs = V.divide(rhs, rmNearestTiesToEven);
1649 if (fs == opDivByZero)
1652 int parts = partCount();
1653 integerPart *x = new integerPart[parts];
1655 fs = V.convertToInteger(x, parts * integerPartWidth, true,
1656 rmNearestTiesToEven, &ignored);
1657 if (fs==opInvalidOp)
1660 fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1661 rmNearestTiesToEven);
1662 assert(fs==opOK); // should always work
1664 fs = V.multiply(rhs, rmNearestTiesToEven);
1665 assert(fs==opOK || fs==opInexact); // should not overflow or underflow
1667 fs = subtract(V, rmNearestTiesToEven);
1668 assert(fs==opOK || fs==opInexact); // likewise
1671 sign = origSign; // IEEE754 requires this
1676 /* Normalized llvm frem (C fmod).
1677 This is not currently correct in all cases. */
1679 APFloat::mod(const APFloat &rhs, roundingMode rounding_mode)
1682 assertArithmeticOK(*semantics);
1683 fs = modSpecials(rhs);
1685 if (category == fcNormal && rhs.category == fcNormal) {
1687 unsigned int origSign = sign;
1689 fs = V.divide(rhs, rmNearestTiesToEven);
1690 if (fs == opDivByZero)
1693 int parts = partCount();
1694 integerPart *x = new integerPart[parts];
1696 fs = V.convertToInteger(x, parts * integerPartWidth, true,
1697 rmTowardZero, &ignored);
1698 if (fs==opInvalidOp)
1701 fs = V.convertFromZeroExtendedInteger(x, parts * integerPartWidth, true,
1702 rmNearestTiesToEven);
1703 assert(fs==opOK); // should always work
1705 fs = V.multiply(rhs, rounding_mode);
1706 assert(fs==opOK || fs==opInexact); // should not overflow or underflow
1708 fs = subtract(V, rounding_mode);
1709 assert(fs==opOK || fs==opInexact); // likewise
1712 sign = origSign; // IEEE754 requires this
1718 /* Normalized fused-multiply-add. */
1720 APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1721 const APFloat &addend,
1722 roundingMode rounding_mode)
1726 assertArithmeticOK(*semantics);
1728 /* Post-multiplication sign, before addition. */
1729 sign ^= multiplicand.sign;
1731 /* If and only if all arguments are normal do we need to do an
1732 extended-precision calculation. */
1733 if(category == fcNormal
1734 && multiplicand.category == fcNormal
1735 && addend.category == fcNormal) {
1736 lostFraction lost_fraction;
1738 lost_fraction = multiplySignificand(multiplicand, &addend);
1739 fs = normalize(rounding_mode, lost_fraction);
1740 if(lost_fraction != lfExactlyZero)
1741 fs = (opStatus) (fs | opInexact);
1743 /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1744 positive zero unless rounding to minus infinity, except that
1745 adding two like-signed zeroes gives that zero. */
1746 if(category == fcZero && sign != addend.sign)
1747 sign = (rounding_mode == rmTowardNegative);
1749 fs = multiplySpecials(multiplicand);
1751 /* FS can only be opOK or opInvalidOp. There is no more work
1752 to do in the latter case. The IEEE-754R standard says it is
1753 implementation-defined in this case whether, if ADDEND is a
1754 quiet NaN, we raise invalid op; this implementation does so.
1756 If we need to do the addition we can do so with normal
1759 fs = addOrSubtract(addend, rounding_mode, false);
1765 /* Comparison requires normalized numbers. */
1767 APFloat::compare(const APFloat &rhs) const
1771 assertArithmeticOK(*semantics);
1772 assert(semantics == rhs.semantics);
1774 switch (convolve(category, rhs.category)) {
1776 llvm_unreachable(0);
1778 case convolve(fcNaN, fcZero):
1779 case convolve(fcNaN, fcNormal):
1780 case convolve(fcNaN, fcInfinity):
1781 case convolve(fcNaN, fcNaN):
1782 case convolve(fcZero, fcNaN):
1783 case convolve(fcNormal, fcNaN):
1784 case convolve(fcInfinity, fcNaN):
1785 return cmpUnordered;
1787 case convolve(fcInfinity, fcNormal):
1788 case convolve(fcInfinity, fcZero):
1789 case convolve(fcNormal, fcZero):
1793 return cmpGreaterThan;
1795 case convolve(fcNormal, fcInfinity):
1796 case convolve(fcZero, fcInfinity):
1797 case convolve(fcZero, fcNormal):
1799 return cmpGreaterThan;
1803 case convolve(fcInfinity, fcInfinity):
1804 if(sign == rhs.sign)
1809 return cmpGreaterThan;
1811 case convolve(fcZero, fcZero):
1814 case convolve(fcNormal, fcNormal):
1818 /* Two normal numbers. Do they have the same sign? */
1819 if(sign != rhs.sign) {
1821 result = cmpLessThan;
1823 result = cmpGreaterThan;
1825 /* Compare absolute values; invert result if negative. */
1826 result = compareAbsoluteValue(rhs);
1829 if(result == cmpLessThan)
1830 result = cmpGreaterThan;
1831 else if(result == cmpGreaterThan)
1832 result = cmpLessThan;
1839 /// APFloat::convert - convert a value of one floating point type to another.
1840 /// The return value corresponds to the IEEE754 exceptions. *losesInfo
1841 /// records whether the transformation lost information, i.e. whether
1842 /// converting the result back to the original type will produce the
1843 /// original value (this is almost the same as return value==fsOK, but there
1844 /// are edge cases where this is not so).
1847 APFloat::convert(const fltSemantics &toSemantics,
1848 roundingMode rounding_mode, bool *losesInfo)
1850 lostFraction lostFraction;
1851 unsigned int newPartCount, oldPartCount;
1854 assertArithmeticOK(*semantics);
1855 assertArithmeticOK(toSemantics);
1856 lostFraction = lfExactlyZero;
1857 newPartCount = partCountForBits(toSemantics.precision + 1);
1858 oldPartCount = partCount();
1860 /* Handle storage complications. If our new form is wider,
1861 re-allocate our bit pattern into wider storage. If it is
1862 narrower, we ignore the excess parts, but if narrowing to a
1863 single part we need to free the old storage.
1864 Be careful not to reference significandParts for zeroes
1865 and infinities, since it aborts. */
1866 if (newPartCount > oldPartCount) {
1867 integerPart *newParts;
1868 newParts = new integerPart[newPartCount];
1869 APInt::tcSet(newParts, 0, newPartCount);
1870 if (category==fcNormal || category==fcNaN)
1871 APInt::tcAssign(newParts, significandParts(), oldPartCount);
1873 significand.parts = newParts;
1874 } else if (newPartCount < oldPartCount) {
1875 /* Capture any lost fraction through truncation of parts so we get
1876 correct rounding whilst normalizing. */
1877 if (category==fcNormal)
1878 lostFraction = lostFractionThroughTruncation
1879 (significandParts(), oldPartCount, toSemantics.precision);
1880 if (newPartCount == 1) {
1881 integerPart newPart = 0;
1882 if (category==fcNormal || category==fcNaN)
1883 newPart = significandParts()[0];
1885 significand.part = newPart;
1889 if(category == fcNormal) {
1890 /* Re-interpret our bit-pattern. */
1891 exponent += toSemantics.precision - semantics->precision;
1892 semantics = &toSemantics;
1893 fs = normalize(rounding_mode, lostFraction);
1894 *losesInfo = (fs != opOK);
1895 } else if (category == fcNaN) {
1896 int shift = toSemantics.precision - semantics->precision;
1897 // Do this now so significandParts gets the right answer
1898 const fltSemantics *oldSemantics = semantics;
1899 semantics = &toSemantics;
1901 // No normalization here, just truncate
1903 APInt::tcShiftLeft(significandParts(), newPartCount, shift);
1904 else if (shift < 0) {
1905 unsigned ushift = -shift;
1906 // Figure out if we are losing information. This happens
1907 // if are shifting out something other than 0s, or if the x87 long
1908 // double input did not have its integer bit set (pseudo-NaN), or if the
1909 // x87 long double input did not have its QNan bit set (because the x87
1910 // hardware sets this bit when converting a lower-precision NaN to
1911 // x87 long double).
1912 if (APInt::tcLSB(significandParts(), newPartCount) < ushift)
1914 if (oldSemantics == &APFloat::x87DoubleExtended &&
1915 (!(*significandParts() & 0x8000000000000000ULL) ||
1916 !(*significandParts() & 0x4000000000000000ULL)))
1918 APInt::tcShiftRight(significandParts(), newPartCount, ushift);
1920 // gcc forces the Quiet bit on, which means (float)(double)(float_sNan)
1921 // does not give you back the same bits. This is dubious, and we
1922 // don't currently do it. You're really supposed to get
1923 // an invalid operation signal at runtime, but nobody does that.
1926 semantics = &toSemantics;
1934 /* Convert a floating point number to an integer according to the
1935 rounding mode. If the rounded integer value is out of range this
1936 returns an invalid operation exception and the contents of the
1937 destination parts are unspecified. If the rounded value is in
1938 range but the floating point number is not the exact integer, the C
1939 standard doesn't require an inexact exception to be raised. IEEE
1940 854 does require it so we do that.
1942 Note that for conversions to integer type the C standard requires
1943 round-to-zero to always be used. */
1945 APFloat::convertToSignExtendedInteger(integerPart *parts, unsigned int width,
1947 roundingMode rounding_mode,
1948 bool *isExact) const
1950 lostFraction lost_fraction;
1951 const integerPart *src;
1952 unsigned int dstPartsCount, truncatedBits;
1954 assertArithmeticOK(*semantics);
1958 /* Handle the three special cases first. */
1959 if(category == fcInfinity || category == fcNaN)
1962 dstPartsCount = partCountForBits(width);
1964 if(category == fcZero) {
1965 APInt::tcSet(parts, 0, dstPartsCount);
1966 // Negative zero can't be represented as an int.
1971 src = significandParts();
1973 /* Step 1: place our absolute value, with any fraction truncated, in
1976 /* Our absolute value is less than one; truncate everything. */
1977 APInt::tcSet(parts, 0, dstPartsCount);
1978 /* For exponent -1 the integer bit represents .5, look at that.
1979 For smaller exponents leftmost truncated bit is 0. */
1980 truncatedBits = semantics->precision -1U - exponent;
1982 /* We want the most significant (exponent + 1) bits; the rest are
1984 unsigned int bits = exponent + 1U;
1986 /* Hopelessly large in magnitude? */
1990 if (bits < semantics->precision) {
1991 /* We truncate (semantics->precision - bits) bits. */
1992 truncatedBits = semantics->precision - bits;
1993 APInt::tcExtract(parts, dstPartsCount, src, bits, truncatedBits);
1995 /* We want at least as many bits as are available. */
1996 APInt::tcExtract(parts, dstPartsCount, src, semantics->precision, 0);
1997 APInt::tcShiftLeft(parts, dstPartsCount, bits - semantics->precision);
2002 /* Step 2: work out any lost fraction, and increment the absolute
2003 value if we would round away from zero. */
2004 if (truncatedBits) {
2005 lost_fraction = lostFractionThroughTruncation(src, partCount(),
2007 if (lost_fraction != lfExactlyZero
2008 && roundAwayFromZero(rounding_mode, lost_fraction, truncatedBits)) {
2009 if (APInt::tcIncrement(parts, dstPartsCount))
2010 return opInvalidOp; /* Overflow. */
2013 lost_fraction = lfExactlyZero;
2016 /* Step 3: check if we fit in the destination. */
2017 unsigned int omsb = APInt::tcMSB(parts, dstPartsCount) + 1;
2021 /* Negative numbers cannot be represented as unsigned. */
2025 /* It takes omsb bits to represent the unsigned integer value.
2026 We lose a bit for the sign, but care is needed as the
2027 maximally negative integer is a special case. */
2028 if (omsb == width && APInt::tcLSB(parts, dstPartsCount) + 1 != omsb)
2031 /* This case can happen because of rounding. */
2036 APInt::tcNegate (parts, dstPartsCount);
2038 if (omsb >= width + !isSigned)
2042 if (lost_fraction == lfExactlyZero) {
2049 /* Same as convertToSignExtendedInteger, except we provide
2050 deterministic values in case of an invalid operation exception,
2051 namely zero for NaNs and the minimal or maximal value respectively
2052 for underflow or overflow.
2053 The *isExact output tells whether the result is exact, in the sense
2054 that converting it back to the original floating point type produces
2055 the original value. This is almost equivalent to result==opOK,
2056 except for negative zeroes.
2059 APFloat::convertToInteger(integerPart *parts, unsigned int width,
2061 roundingMode rounding_mode, bool *isExact) const
2065 fs = convertToSignExtendedInteger(parts, width, isSigned, rounding_mode,
2068 if (fs == opInvalidOp) {
2069 unsigned int bits, dstPartsCount;
2071 dstPartsCount = partCountForBits(width);
2073 if (category == fcNaN)
2078 bits = width - isSigned;
2080 APInt::tcSetLeastSignificantBits(parts, dstPartsCount, bits);
2081 if (sign && isSigned)
2082 APInt::tcShiftLeft(parts, dstPartsCount, width - 1);
2088 /* Convert an unsigned integer SRC to a floating point number,
2089 rounding according to ROUNDING_MODE. The sign of the floating
2090 point number is not modified. */
2092 APFloat::convertFromUnsignedParts(const integerPart *src,
2093 unsigned int srcCount,
2094 roundingMode rounding_mode)
2096 unsigned int omsb, precision, dstCount;
2098 lostFraction lost_fraction;
2100 assertArithmeticOK(*semantics);
2101 category = fcNormal;
2102 omsb = APInt::tcMSB(src, srcCount) + 1;
2103 dst = significandParts();
2104 dstCount = partCount();
2105 precision = semantics->precision;
2107 /* We want the most significant PRECISON bits of SRC. There may not
2108 be that many; extract what we can. */
2109 if (precision <= omsb) {
2110 exponent = omsb - 1;
2111 lost_fraction = lostFractionThroughTruncation(src, srcCount,
2113 APInt::tcExtract(dst, dstCount, src, precision, omsb - precision);
2115 exponent = precision - 1;
2116 lost_fraction = lfExactlyZero;
2117 APInt::tcExtract(dst, dstCount, src, omsb, 0);
2120 return normalize(rounding_mode, lost_fraction);
2124 APFloat::convertFromAPInt(const APInt &Val,
2126 roundingMode rounding_mode)
2128 unsigned int partCount = Val.getNumWords();
2132 if (isSigned && api.isNegative()) {
2137 return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2140 /* Convert a two's complement integer SRC to a floating point number,
2141 rounding according to ROUNDING_MODE. ISSIGNED is true if the
2142 integer is signed, in which case it must be sign-extended. */
2144 APFloat::convertFromSignExtendedInteger(const integerPart *src,
2145 unsigned int srcCount,
2147 roundingMode rounding_mode)
2151 assertArithmeticOK(*semantics);
2153 && APInt::tcExtractBit(src, srcCount * integerPartWidth - 1)) {
2156 /* If we're signed and negative negate a copy. */
2158 copy = new integerPart[srcCount];
2159 APInt::tcAssign(copy, src, srcCount);
2160 APInt::tcNegate(copy, srcCount);
2161 status = convertFromUnsignedParts(copy, srcCount, rounding_mode);
2165 status = convertFromUnsignedParts(src, srcCount, rounding_mode);
2171 /* FIXME: should this just take a const APInt reference? */
2173 APFloat::convertFromZeroExtendedInteger(const integerPart *parts,
2174 unsigned int width, bool isSigned,
2175 roundingMode rounding_mode)
2177 unsigned int partCount = partCountForBits(width);
2178 APInt api = APInt(width, partCount, parts);
2181 if(isSigned && APInt::tcExtractBit(parts, width - 1)) {
2186 return convertFromUnsignedParts(api.getRawData(), partCount, rounding_mode);
2190 APFloat::convertFromHexadecimalString(const StringRef &s,
2191 roundingMode rounding_mode)
2193 lostFraction lost_fraction = lfExactlyZero;
2194 integerPart *significand;
2195 unsigned int bitPos, partsCount;
2196 StringRef::iterator dot, firstSignificantDigit;
2200 category = fcNormal;
2202 significand = significandParts();
2203 partsCount = partCount();
2204 bitPos = partsCount * integerPartWidth;
2206 /* Skip leading zeroes and any (hexa)decimal point. */
2207 StringRef::iterator begin = s.begin();
2208 StringRef::iterator end = s.end();
2209 StringRef::iterator p = skipLeadingZeroesAndAnyDot(begin, end, &dot);
2210 firstSignificantDigit = p;
2213 integerPart hex_value;
2216 assert(dot == end && "String contains multiple dots");
2223 hex_value = hexDigitValue(*p);
2224 if(hex_value == -1U) {
2233 /* Store the number whilst 4-bit nibbles remain. */
2236 hex_value <<= bitPos % integerPartWidth;
2237 significand[bitPos / integerPartWidth] |= hex_value;
2239 lost_fraction = trailingHexadecimalFraction(p, end, hex_value);
2240 while(p != end && hexDigitValue(*p) != -1U)
2247 /* Hex floats require an exponent but not a hexadecimal point. */
2248 assert(p != end && "Hex strings require an exponent");
2249 assert((*p == 'p' || *p == 'P') && "Invalid character in significand");
2250 assert(p != begin && "Significand has no digits");
2251 assert((dot == end || p - begin != 1) && "Significand has no digits");
2253 /* Ignore the exponent if we are zero. */
2254 if(p != firstSignificantDigit) {
2257 /* Implicit hexadecimal point? */
2261 /* Calculate the exponent adjustment implicit in the number of
2262 significant digits. */
2263 expAdjustment = static_cast<int>(dot - firstSignificantDigit);
2264 if(expAdjustment < 0)
2266 expAdjustment = expAdjustment * 4 - 1;
2268 /* Adjust for writing the significand starting at the most
2269 significant nibble. */
2270 expAdjustment += semantics->precision;
2271 expAdjustment -= partsCount * integerPartWidth;
2273 /* Adjust for the given exponent. */
2274 exponent = totalExponent(p + 1, end, expAdjustment);
2277 return normalize(rounding_mode, lost_fraction);
2281 APFloat::roundSignificandWithExponent(const integerPart *decSigParts,
2282 unsigned sigPartCount, int exp,
2283 roundingMode rounding_mode)
2285 unsigned int parts, pow5PartCount;
2286 fltSemantics calcSemantics = { 32767, -32767, 0, true };
2287 integerPart pow5Parts[maxPowerOfFiveParts];
2290 isNearest = (rounding_mode == rmNearestTiesToEven
2291 || rounding_mode == rmNearestTiesToAway);
2293 parts = partCountForBits(semantics->precision + 11);
2295 /* Calculate pow(5, abs(exp)). */
2296 pow5PartCount = powerOf5(pow5Parts, exp >= 0 ? exp: -exp);
2298 for (;; parts *= 2) {
2299 opStatus sigStatus, powStatus;
2300 unsigned int excessPrecision, truncatedBits;
2302 calcSemantics.precision = parts * integerPartWidth - 1;
2303 excessPrecision = calcSemantics.precision - semantics->precision;
2304 truncatedBits = excessPrecision;
2306 APFloat decSig(calcSemantics, fcZero, sign);
2307 APFloat pow5(calcSemantics, fcZero, false);
2309 sigStatus = decSig.convertFromUnsignedParts(decSigParts, sigPartCount,
2310 rmNearestTiesToEven);
2311 powStatus = pow5.convertFromUnsignedParts(pow5Parts, pow5PartCount,
2312 rmNearestTiesToEven);
2313 /* Add exp, as 10^n = 5^n * 2^n. */
2314 decSig.exponent += exp;
2316 lostFraction calcLostFraction;
2317 integerPart HUerr, HUdistance;
2318 unsigned int powHUerr;
2321 /* multiplySignificand leaves the precision-th bit set to 1. */
2322 calcLostFraction = decSig.multiplySignificand(pow5, NULL);
2323 powHUerr = powStatus != opOK;
2325 calcLostFraction = decSig.divideSignificand(pow5);
2326 /* Denormal numbers have less precision. */
2327 if (decSig.exponent < semantics->minExponent) {
2328 excessPrecision += (semantics->minExponent - decSig.exponent);
2329 truncatedBits = excessPrecision;
2330 if (excessPrecision > calcSemantics.precision)
2331 excessPrecision = calcSemantics.precision;
2333 /* Extra half-ulp lost in reciprocal of exponent. */
2334 powHUerr = (powStatus == opOK && calcLostFraction == lfExactlyZero) ? 0:2;
2337 /* Both multiplySignificand and divideSignificand return the
2338 result with the integer bit set. */
2339 assert(APInt::tcExtractBit
2340 (decSig.significandParts(), calcSemantics.precision - 1) == 1);
2342 HUerr = HUerrBound(calcLostFraction != lfExactlyZero, sigStatus != opOK,
2344 HUdistance = 2 * ulpsFromBoundary(decSig.significandParts(),
2345 excessPrecision, isNearest);
2347 /* Are we guaranteed to round correctly if we truncate? */
2348 if (HUdistance >= HUerr) {
2349 APInt::tcExtract(significandParts(), partCount(), decSig.significandParts(),
2350 calcSemantics.precision - excessPrecision,
2352 /* Take the exponent of decSig. If we tcExtract-ed less bits
2353 above we must adjust our exponent to compensate for the
2354 implicit right shift. */
2355 exponent = (decSig.exponent + semantics->precision
2356 - (calcSemantics.precision - excessPrecision));
2357 calcLostFraction = lostFractionThroughTruncation(decSig.significandParts(),
2360 return normalize(rounding_mode, calcLostFraction);
2366 APFloat::convertFromDecimalString(const StringRef &str, roundingMode rounding_mode)
2371 /* Scan the text. */
2372 StringRef::iterator p = str.begin();
2373 interpretDecimal(p, str.end(), &D);
2375 /* Handle the quick cases. First the case of no significant digits,
2376 i.e. zero, and then exponents that are obviously too large or too
2377 small. Writing L for log 10 / log 2, a number d.ddddd*10^exp
2378 definitely overflows if
2380 (exp - 1) * L >= maxExponent
2382 and definitely underflows to zero where
2384 (exp + 1) * L <= minExponent - precision
2386 With integer arithmetic the tightest bounds for L are
2388 93/28 < L < 196/59 [ numerator <= 256 ]
2389 42039/12655 < L < 28738/8651 [ numerator <= 65536 ]
2392 if (decDigitValue(*D.firstSigDigit) >= 10U) {
2396 /* Check whether the normalized exponent is high enough to overflow
2397 max during the log-rebasing in the max-exponent check below. */
2398 } else if (D.normalizedExponent - 1 > INT_MAX / 42039) {
2399 fs = handleOverflow(rounding_mode);
2401 /* If it wasn't, then it also wasn't high enough to overflow max
2402 during the log-rebasing in the min-exponent check. Check that it
2403 won't overflow min in either check, then perform the min-exponent
2405 } else if (D.normalizedExponent - 1 < INT_MIN / 42039 ||
2406 (D.normalizedExponent + 1) * 28738 <=
2407 8651 * (semantics->minExponent - (int) semantics->precision)) {
2408 /* Underflow to zero and round. */
2410 fs = normalize(rounding_mode, lfLessThanHalf);
2412 /* We can finally safely perform the max-exponent check. */
2413 } else if ((D.normalizedExponent - 1) * 42039
2414 >= 12655 * semantics->maxExponent) {
2415 /* Overflow and round. */
2416 fs = handleOverflow(rounding_mode);
2418 integerPart *decSignificand;
2419 unsigned int partCount;
2421 /* A tight upper bound on number of bits required to hold an
2422 N-digit decimal integer is N * 196 / 59. Allocate enough space
2423 to hold the full significand, and an extra part required by
2425 partCount = static_cast<unsigned int>(D.lastSigDigit - D.firstSigDigit) + 1;
2426 partCount = partCountForBits(1 + 196 * partCount / 59);
2427 decSignificand = new integerPart[partCount + 1];
2430 /* Convert to binary efficiently - we do almost all multiplication
2431 in an integerPart. When this would overflow do we do a single
2432 bignum multiplication, and then revert again to multiplication
2433 in an integerPart. */
2435 integerPart decValue, val, multiplier;
2443 if (p == str.end()) {
2447 decValue = decDigitValue(*p++);
2448 assert(decValue < 10U && "Invalid character in significand");
2450 val = val * 10 + decValue;
2451 /* The maximum number that can be multiplied by ten with any
2452 digit added without overflowing an integerPart. */
2453 } while (p <= D.lastSigDigit && multiplier <= (~ (integerPart) 0 - 9) / 10);
2455 /* Multiply out the current part. */
2456 APInt::tcMultiplyPart(decSignificand, decSignificand, multiplier, val,
2457 partCount, partCount + 1, false);
2459 /* If we used another part (likely but not guaranteed), increase
2461 if (decSignificand[partCount])
2463 } while (p <= D.lastSigDigit);
2465 category = fcNormal;
2466 fs = roundSignificandWithExponent(decSignificand, partCount,
2467 D.exponent, rounding_mode);
2469 delete [] decSignificand;
2476 APFloat::convertFromString(const StringRef &str, roundingMode rounding_mode)
2478 assertArithmeticOK(*semantics);
2479 assert(!str.empty() && "Invalid string length");
2481 /* Handle a leading minus sign. */
2482 StringRef::iterator p = str.begin();
2483 size_t slen = str.size();
2484 sign = *p == '-' ? 1 : 0;
2485 if(*p == '-' || *p == '+') {
2488 assert(slen && "String has no digits");
2491 if(slen >= 2 && p[0] == '0' && (p[1] == 'x' || p[1] == 'X')) {
2492 assert(slen - 2 && "Invalid string");
2493 return convertFromHexadecimalString(StringRef(p + 2, slen - 2),
2497 return convertFromDecimalString(StringRef(p, slen), rounding_mode);
2500 /* Write out a hexadecimal representation of the floating point value
2501 to DST, which must be of sufficient size, in the C99 form
2502 [-]0xh.hhhhp[+-]d. Return the number of characters written,
2503 excluding the terminating NUL.
2505 If UPPERCASE, the output is in upper case, otherwise in lower case.
2507 HEXDIGITS digits appear altogether, rounding the value if
2508 necessary. If HEXDIGITS is 0, the minimal precision to display the
2509 number precisely is used instead. If nothing would appear after
2510 the decimal point it is suppressed.
2512 The decimal exponent is always printed and has at least one digit.
2513 Zero values display an exponent of zero. Infinities and NaNs
2514 appear as "infinity" or "nan" respectively.
2516 The above rules are as specified by C99. There is ambiguity about
2517 what the leading hexadecimal digit should be. This implementation
2518 uses whatever is necessary so that the exponent is displayed as
2519 stored. This implies the exponent will fall within the IEEE format
2520 range, and the leading hexadecimal digit will be 0 (for denormals),
2521 1 (normal numbers) or 2 (normal numbers rounded-away-from-zero with
2522 any other digits zero).
2525 APFloat::convertToHexString(char *dst, unsigned int hexDigits,
2526 bool upperCase, roundingMode rounding_mode) const
2530 assertArithmeticOK(*semantics);
2538 memcpy (dst, upperCase ? infinityU: infinityL, sizeof infinityU - 1);
2539 dst += sizeof infinityL - 1;
2543 memcpy (dst, upperCase ? NaNU: NaNL, sizeof NaNU - 1);
2544 dst += sizeof NaNU - 1;
2549 *dst++ = upperCase ? 'X': 'x';
2551 if (hexDigits > 1) {
2553 memset (dst, '0', hexDigits - 1);
2554 dst += hexDigits - 1;
2556 *dst++ = upperCase ? 'P': 'p';
2561 dst = convertNormalToHexString (dst, hexDigits, upperCase, rounding_mode);
2567 return static_cast<unsigned int>(dst - p);
2570 /* Does the hard work of outputting the correctly rounded hexadecimal
2571 form of a normal floating point number with the specified number of
2572 hexadecimal digits. If HEXDIGITS is zero the minimum number of
2573 digits necessary to print the value precisely is output. */
2575 APFloat::convertNormalToHexString(char *dst, unsigned int hexDigits,
2577 roundingMode rounding_mode) const
2579 unsigned int count, valueBits, shift, partsCount, outputDigits;
2580 const char *hexDigitChars;
2581 const integerPart *significand;
2586 *dst++ = upperCase ? 'X': 'x';
2589 hexDigitChars = upperCase ? hexDigitsUpper: hexDigitsLower;
2591 significand = significandParts();
2592 partsCount = partCount();
2594 /* +3 because the first digit only uses the single integer bit, so
2595 we have 3 virtual zero most-significant-bits. */
2596 valueBits = semantics->precision + 3;
2597 shift = integerPartWidth - valueBits % integerPartWidth;
2599 /* The natural number of digits required ignoring trailing
2600 insignificant zeroes. */
2601 outputDigits = (valueBits - significandLSB () + 3) / 4;
2603 /* hexDigits of zero means use the required number for the
2604 precision. Otherwise, see if we are truncating. If we are,
2605 find out if we need to round away from zero. */
2607 if (hexDigits < outputDigits) {
2608 /* We are dropping non-zero bits, so need to check how to round.
2609 "bits" is the number of dropped bits. */
2611 lostFraction fraction;
2613 bits = valueBits - hexDigits * 4;
2614 fraction = lostFractionThroughTruncation (significand, partsCount, bits);
2615 roundUp = roundAwayFromZero(rounding_mode, fraction, bits);
2617 outputDigits = hexDigits;
2620 /* Write the digits consecutively, and start writing in the location
2621 of the hexadecimal point. We move the most significant digit
2622 left and add the hexadecimal point later. */
2625 count = (valueBits + integerPartWidth - 1) / integerPartWidth;
2627 while (outputDigits && count) {
2630 /* Put the most significant integerPartWidth bits in "part". */
2631 if (--count == partsCount)
2632 part = 0; /* An imaginary higher zero part. */
2634 part = significand[count] << shift;
2637 part |= significand[count - 1] >> (integerPartWidth - shift);
2639 /* Convert as much of "part" to hexdigits as we can. */
2640 unsigned int curDigits = integerPartWidth / 4;
2642 if (curDigits > outputDigits)
2643 curDigits = outputDigits;
2644 dst += partAsHex (dst, part, curDigits, hexDigitChars);
2645 outputDigits -= curDigits;
2651 /* Note that hexDigitChars has a trailing '0'. */
2654 *q = hexDigitChars[hexDigitValue (*q) + 1];
2655 } while (*q == '0');
2658 /* Add trailing zeroes. */
2659 memset (dst, '0', outputDigits);
2660 dst += outputDigits;
2663 /* Move the most significant digit to before the point, and if there
2664 is something after the decimal point add it. This must come
2665 after rounding above. */
2672 /* Finally output the exponent. */
2673 *dst++ = upperCase ? 'P': 'p';
2675 return writeSignedDecimal (dst, exponent);
2678 // For good performance it is desirable for different APFloats
2679 // to produce different integers.
2681 APFloat::getHashValue() const
2683 if (category==fcZero) return sign<<8 | semantics->precision ;
2684 else if (category==fcInfinity) return sign<<9 | semantics->precision;
2685 else if (category==fcNaN) return 1<<10 | semantics->precision;
2687 uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
2688 const integerPart* p = significandParts();
2689 for (int i=partCount(); i>0; i--, p++)
2690 hash ^= ((uint32_t)*p) ^ (uint32_t)((*p)>>32);
2695 // Conversion from APFloat to/from host float/double. It may eventually be
2696 // possible to eliminate these and have everybody deal with APFloats, but that
2697 // will take a while. This approach will not easily extend to long double.
2698 // Current implementation requires integerPartWidth==64, which is correct at
2699 // the moment but could be made more general.
2701 // Denormals have exponent minExponent in APFloat, but minExponent-1 in
2702 // the actual IEEE respresentations. We compensate for that here.
2705 APFloat::convertF80LongDoubleAPFloatToAPInt() const
2707 assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended);
2708 assert(partCount()==2);
2710 uint64_t myexponent, mysignificand;
2712 if (category==fcNormal) {
2713 myexponent = exponent+16383; //bias
2714 mysignificand = significandParts()[0];
2715 if (myexponent==1 && !(mysignificand & 0x8000000000000000ULL))
2716 myexponent = 0; // denormal
2717 } else if (category==fcZero) {
2720 } else if (category==fcInfinity) {
2721 myexponent = 0x7fff;
2722 mysignificand = 0x8000000000000000ULL;
2724 assert(category == fcNaN && "Unknown category");
2725 myexponent = 0x7fff;
2726 mysignificand = significandParts()[0];
2730 words[0] = mysignificand;
2731 words[1] = ((uint64_t)(sign & 1) << 15) |
2732 (myexponent & 0x7fffLL);
2733 return APInt(80, 2, words);
2737 APFloat::convertPPCDoubleDoubleAPFloatToAPInt() const
2739 assert(semantics == (const llvm::fltSemantics*)&PPCDoubleDouble);
2740 assert(partCount()==2);
2742 uint64_t myexponent, mysignificand, myexponent2, mysignificand2;
2744 if (category==fcNormal) {
2745 myexponent = exponent + 1023; //bias
2746 myexponent2 = exponent2 + 1023;
2747 mysignificand = significandParts()[0];
2748 mysignificand2 = significandParts()[1];
2749 if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2750 myexponent = 0; // denormal
2751 if (myexponent2==1 && !(mysignificand2 & 0x10000000000000LL))
2752 myexponent2 = 0; // denormal
2753 } else if (category==fcZero) {
2758 } else if (category==fcInfinity) {
2764 assert(category == fcNaN && "Unknown category");
2766 mysignificand = significandParts()[0];
2767 myexponent2 = exponent2;
2768 mysignificand2 = significandParts()[1];
2772 words[0] = ((uint64_t)(sign & 1) << 63) |
2773 ((myexponent & 0x7ff) << 52) |
2774 (mysignificand & 0xfffffffffffffLL);
2775 words[1] = ((uint64_t)(sign2 & 1) << 63) |
2776 ((myexponent2 & 0x7ff) << 52) |
2777 (mysignificand2 & 0xfffffffffffffLL);
2778 return APInt(128, 2, words);
2782 APFloat::convertQuadrupleAPFloatToAPInt() const
2784 assert(semantics == (const llvm::fltSemantics*)&IEEEquad);
2785 assert(partCount()==2);
2787 uint64_t myexponent, mysignificand, mysignificand2;
2789 if (category==fcNormal) {
2790 myexponent = exponent+16383; //bias
2791 mysignificand = significandParts()[0];
2792 mysignificand2 = significandParts()[1];
2793 if (myexponent==1 && !(mysignificand2 & 0x1000000000000LL))
2794 myexponent = 0; // denormal
2795 } else if (category==fcZero) {
2797 mysignificand = mysignificand2 = 0;
2798 } else if (category==fcInfinity) {
2799 myexponent = 0x7fff;
2800 mysignificand = mysignificand2 = 0;
2802 assert(category == fcNaN && "Unknown category!");
2803 myexponent = 0x7fff;
2804 mysignificand = significandParts()[0];
2805 mysignificand2 = significandParts()[1];
2809 words[0] = mysignificand;
2810 words[1] = ((uint64_t)(sign & 1) << 63) |
2811 ((myexponent & 0x7fff) << 48) |
2812 (mysignificand2 & 0xffffffffffffLL);
2814 return APInt(128, 2, words);
2818 APFloat::convertDoubleAPFloatToAPInt() const
2820 assert(semantics == (const llvm::fltSemantics*)&IEEEdouble);
2821 assert(partCount()==1);
2823 uint64_t myexponent, mysignificand;
2825 if (category==fcNormal) {
2826 myexponent = exponent+1023; //bias
2827 mysignificand = *significandParts();
2828 if (myexponent==1 && !(mysignificand & 0x10000000000000LL))
2829 myexponent = 0; // denormal
2830 } else if (category==fcZero) {
2833 } else if (category==fcInfinity) {
2837 assert(category == fcNaN && "Unknown category!");
2839 mysignificand = *significandParts();
2842 return APInt(64, ((((uint64_t)(sign & 1) << 63) |
2843 ((myexponent & 0x7ff) << 52) |
2844 (mysignificand & 0xfffffffffffffLL))));
2848 APFloat::convertFloatAPFloatToAPInt() const
2850 assert(semantics == (const llvm::fltSemantics*)&IEEEsingle);
2851 assert(partCount()==1);
2853 uint32_t myexponent, mysignificand;
2855 if (category==fcNormal) {
2856 myexponent = exponent+127; //bias
2857 mysignificand = (uint32_t)*significandParts();
2858 if (myexponent == 1 && !(mysignificand & 0x800000))
2859 myexponent = 0; // denormal
2860 } else if (category==fcZero) {
2863 } else if (category==fcInfinity) {
2867 assert(category == fcNaN && "Unknown category!");
2869 mysignificand = (uint32_t)*significandParts();
2872 return APInt(32, (((sign&1) << 31) | ((myexponent&0xff) << 23) |
2873 (mysignificand & 0x7fffff)));
2877 APFloat::convertHalfAPFloatToAPInt() const
2879 assert(semantics == (const llvm::fltSemantics*)&IEEEhalf);
2880 assert(partCount()==1);
2882 uint32_t myexponent, mysignificand;
2884 if (category==fcNormal) {
2885 myexponent = exponent+15; //bias
2886 mysignificand = (uint32_t)*significandParts();
2887 if (myexponent == 1 && !(mysignificand & 0x400))
2888 myexponent = 0; // denormal
2889 } else if (category==fcZero) {
2892 } else if (category==fcInfinity) {
2896 assert(category == fcNaN && "Unknown category!");
2898 mysignificand = (uint32_t)*significandParts();
2901 return APInt(16, (((sign&1) << 15) | ((myexponent&0x1f) << 10) |
2902 (mysignificand & 0x3ff)));
2905 // This function creates an APInt that is just a bit map of the floating
2906 // point constant as it would appear in memory. It is not a conversion,
2907 // and treating the result as a normal integer is unlikely to be useful.
2910 APFloat::bitcastToAPInt() const
2912 if (semantics == (const llvm::fltSemantics*)&IEEEhalf)
2913 return convertHalfAPFloatToAPInt();
2915 if (semantics == (const llvm::fltSemantics*)&IEEEsingle)
2916 return convertFloatAPFloatToAPInt();
2918 if (semantics == (const llvm::fltSemantics*)&IEEEdouble)
2919 return convertDoubleAPFloatToAPInt();
2921 if (semantics == (const llvm::fltSemantics*)&IEEEquad)
2922 return convertQuadrupleAPFloatToAPInt();
2924 if (semantics == (const llvm::fltSemantics*)&PPCDoubleDouble)
2925 return convertPPCDoubleDoubleAPFloatToAPInt();
2927 assert(semantics == (const llvm::fltSemantics*)&x87DoubleExtended &&
2929 return convertF80LongDoubleAPFloatToAPInt();
2933 APFloat::convertToFloat() const
2935 assert(semantics == (const llvm::fltSemantics*)&IEEEsingle &&
2936 "Float semantics are not IEEEsingle");
2937 APInt api = bitcastToAPInt();
2938 return api.bitsToFloat();
2942 APFloat::convertToDouble() const
2944 assert(semantics == (const llvm::fltSemantics*)&IEEEdouble &&
2945 "Float semantics are not IEEEdouble");
2946 APInt api = bitcastToAPInt();
2947 return api.bitsToDouble();
2950 /// Integer bit is explicit in this format. Intel hardware (387 and later)
2951 /// does not support these bit patterns:
2952 /// exponent = all 1's, integer bit 0, significand 0 ("pseudoinfinity")
2953 /// exponent = all 1's, integer bit 0, significand nonzero ("pseudoNaN")
2954 /// exponent = 0, integer bit 1 ("pseudodenormal")
2955 /// exponent!=0 nor all 1's, integer bit 0 ("unnormal")
2956 /// At the moment, the first two are treated as NaNs, the second two as Normal.
2958 APFloat::initFromF80LongDoubleAPInt(const APInt &api)
2960 assert(api.getBitWidth()==80);
2961 uint64_t i1 = api.getRawData()[0];
2962 uint64_t i2 = api.getRawData()[1];
2963 uint64_t myexponent = (i2 & 0x7fff);
2964 uint64_t mysignificand = i1;
2966 initialize(&APFloat::x87DoubleExtended);
2967 assert(partCount()==2);
2969 sign = static_cast<unsigned int>(i2>>15);
2970 if (myexponent==0 && mysignificand==0) {
2971 // exponent, significand meaningless
2973 } else if (myexponent==0x7fff && mysignificand==0x8000000000000000ULL) {
2974 // exponent, significand meaningless
2975 category = fcInfinity;
2976 } else if (myexponent==0x7fff && mysignificand!=0x8000000000000000ULL) {
2977 // exponent meaningless
2979 significandParts()[0] = mysignificand;
2980 significandParts()[1] = 0;
2982 category = fcNormal;
2983 exponent = myexponent - 16383;
2984 significandParts()[0] = mysignificand;
2985 significandParts()[1] = 0;
2986 if (myexponent==0) // denormal
2992 APFloat::initFromPPCDoubleDoubleAPInt(const APInt &api)
2994 assert(api.getBitWidth()==128);
2995 uint64_t i1 = api.getRawData()[0];
2996 uint64_t i2 = api.getRawData()[1];
2997 uint64_t myexponent = (i1 >> 52) & 0x7ff;
2998 uint64_t mysignificand = i1 & 0xfffffffffffffLL;
2999 uint64_t myexponent2 = (i2 >> 52) & 0x7ff;
3000 uint64_t mysignificand2 = i2 & 0xfffffffffffffLL;
3002 initialize(&APFloat::PPCDoubleDouble);
3003 assert(partCount()==2);
3005 sign = static_cast<unsigned int>(i1>>63);
3006 sign2 = static_cast<unsigned int>(i2>>63);
3007 if (myexponent==0 && mysignificand==0) {
3008 // exponent, significand meaningless
3009 // exponent2 and significand2 are required to be 0; we don't check
3011 } else if (myexponent==0x7ff && mysignificand==0) {
3012 // exponent, significand meaningless
3013 // exponent2 and significand2 are required to be 0; we don't check
3014 category = fcInfinity;
3015 } else if (myexponent==0x7ff && mysignificand!=0) {
3016 // exponent meaningless. So is the whole second word, but keep it
3019 exponent2 = myexponent2;
3020 significandParts()[0] = mysignificand;
3021 significandParts()[1] = mysignificand2;
3023 category = fcNormal;
3024 // Note there is no category2; the second word is treated as if it is
3025 // fcNormal, although it might be something else considered by itself.
3026 exponent = myexponent - 1023;
3027 exponent2 = myexponent2 - 1023;
3028 significandParts()[0] = mysignificand;
3029 significandParts()[1] = mysignificand2;
3030 if (myexponent==0) // denormal
3033 significandParts()[0] |= 0x10000000000000LL; // integer bit
3037 significandParts()[1] |= 0x10000000000000LL; // integer bit
3042 APFloat::initFromQuadrupleAPInt(const APInt &api)
3044 assert(api.getBitWidth()==128);
3045 uint64_t i1 = api.getRawData()[0];
3046 uint64_t i2 = api.getRawData()[1];
3047 uint64_t myexponent = (i2 >> 48) & 0x7fff;
3048 uint64_t mysignificand = i1;
3049 uint64_t mysignificand2 = i2 & 0xffffffffffffLL;
3051 initialize(&APFloat::IEEEquad);
3052 assert(partCount()==2);
3054 sign = static_cast<unsigned int>(i2>>63);
3055 if (myexponent==0 &&
3056 (mysignificand==0 && mysignificand2==0)) {
3057 // exponent, significand meaningless
3059 } else if (myexponent==0x7fff &&
3060 (mysignificand==0 && mysignificand2==0)) {
3061 // exponent, significand meaningless
3062 category = fcInfinity;
3063 } else if (myexponent==0x7fff &&
3064 (mysignificand!=0 || mysignificand2 !=0)) {
3065 // exponent meaningless
3067 significandParts()[0] = mysignificand;
3068 significandParts()[1] = mysignificand2;
3070 category = fcNormal;
3071 exponent = myexponent - 16383;
3072 significandParts()[0] = mysignificand;
3073 significandParts()[1] = mysignificand2;
3074 if (myexponent==0) // denormal
3077 significandParts()[1] |= 0x1000000000000LL; // integer bit
3082 APFloat::initFromDoubleAPInt(const APInt &api)
3084 assert(api.getBitWidth()==64);
3085 uint64_t i = *api.getRawData();
3086 uint64_t myexponent = (i >> 52) & 0x7ff;
3087 uint64_t mysignificand = i & 0xfffffffffffffLL;
3089 initialize(&APFloat::IEEEdouble);
3090 assert(partCount()==1);
3092 sign = static_cast<unsigned int>(i>>63);
3093 if (myexponent==0 && mysignificand==0) {
3094 // exponent, significand meaningless
3096 } else if (myexponent==0x7ff && mysignificand==0) {
3097 // exponent, significand meaningless
3098 category = fcInfinity;
3099 } else if (myexponent==0x7ff && mysignificand!=0) {
3100 // exponent meaningless
3102 *significandParts() = mysignificand;
3104 category = fcNormal;
3105 exponent = myexponent - 1023;
3106 *significandParts() = mysignificand;
3107 if (myexponent==0) // denormal
3110 *significandParts() |= 0x10000000000000LL; // integer bit
3115 APFloat::initFromFloatAPInt(const APInt & api)
3117 assert(api.getBitWidth()==32);
3118 uint32_t i = (uint32_t)*api.getRawData();
3119 uint32_t myexponent = (i >> 23) & 0xff;
3120 uint32_t mysignificand = i & 0x7fffff;
3122 initialize(&APFloat::IEEEsingle);
3123 assert(partCount()==1);
3126 if (myexponent==0 && mysignificand==0) {
3127 // exponent, significand meaningless
3129 } else if (myexponent==0xff && mysignificand==0) {
3130 // exponent, significand meaningless
3131 category = fcInfinity;
3132 } else if (myexponent==0xff && mysignificand!=0) {
3133 // sign, exponent, significand meaningless
3135 *significandParts() = mysignificand;
3137 category = fcNormal;
3138 exponent = myexponent - 127; //bias
3139 *significandParts() = mysignificand;
3140 if (myexponent==0) // denormal
3143 *significandParts() |= 0x800000; // integer bit
3148 APFloat::initFromHalfAPInt(const APInt & api)
3150 assert(api.getBitWidth()==16);
3151 uint32_t i = (uint32_t)*api.getRawData();
3152 uint32_t myexponent = (i >> 10) & 0x1f;
3153 uint32_t mysignificand = i & 0x3ff;
3155 initialize(&APFloat::IEEEhalf);
3156 assert(partCount()==1);
3159 if (myexponent==0 && mysignificand==0) {
3160 // exponent, significand meaningless
3162 } else if (myexponent==0x1f && mysignificand==0) {
3163 // exponent, significand meaningless
3164 category = fcInfinity;
3165 } else if (myexponent==0x1f && mysignificand!=0) {
3166 // sign, exponent, significand meaningless
3168 *significandParts() = mysignificand;
3170 category = fcNormal;
3171 exponent = myexponent - 15; //bias
3172 *significandParts() = mysignificand;
3173 if (myexponent==0) // denormal
3176 *significandParts() |= 0x400; // integer bit
3180 /// Treat api as containing the bits of a floating point number. Currently
3181 /// we infer the floating point type from the size of the APInt. The
3182 /// isIEEE argument distinguishes between PPC128 and IEEE128 (not meaningful
3183 /// when the size is anything else).
3185 APFloat::initFromAPInt(const APInt& api, bool isIEEE)
3187 if (api.getBitWidth() == 16)
3188 return initFromHalfAPInt(api);
3189 else if (api.getBitWidth() == 32)
3190 return initFromFloatAPInt(api);
3191 else if (api.getBitWidth()==64)
3192 return initFromDoubleAPInt(api);
3193 else if (api.getBitWidth()==80)
3194 return initFromF80LongDoubleAPInt(api);
3195 else if (api.getBitWidth()==128)
3197 initFromQuadrupleAPInt(api) : initFromPPCDoubleDoubleAPInt(api));
3199 llvm_unreachable(0);
3202 APFloat APFloat::getLargest(const fltSemantics &Sem, bool Negative) {
3203 APFloat Val(Sem, fcNormal, Negative);
3205 // We want (in interchange format):
3206 // sign = {Negative}
3208 // significand = 1..1
3210 Val.exponent = Sem.maxExponent; // unbiased
3212 // 1-initialize all bits....
3213 Val.zeroSignificand();
3214 integerPart *significand = Val.significandParts();
3215 unsigned N = partCountForBits(Sem.precision);
3216 for (unsigned i = 0; i != N; ++i)
3217 significand[i] = ~((integerPart) 0);
3219 // ...and then clear the top bits for internal consistency.
3221 &= (((integerPart) 1) << ((Sem.precision % integerPartWidth) - 1)) - 1;
3226 APFloat APFloat::getSmallest(const fltSemantics &Sem, bool Negative) {
3227 APFloat Val(Sem, fcNormal, Negative);
3229 // We want (in interchange format):
3230 // sign = {Negative}
3232 // significand = 0..01
3234 Val.exponent = Sem.minExponent; // unbiased
3235 Val.zeroSignificand();
3236 Val.significandParts()[0] = 1;
3240 APFloat APFloat::getSmallestNormalized(const fltSemantics &Sem, bool Negative) {
3241 APFloat Val(Sem, fcNormal, Negative);
3243 // We want (in interchange format):
3244 // sign = {Negative}
3246 // significand = 10..0
3248 Val.exponent = Sem.minExponent;
3249 Val.zeroSignificand();
3250 Val.significandParts()[partCountForBits(Sem.precision)-1]
3251 |= (((integerPart) 1) << ((Sem.precision % integerPartWidth) - 1));
3256 APFloat::APFloat(const APInt& api, bool isIEEE)
3258 initFromAPInt(api, isIEEE);
3261 APFloat::APFloat(float f)
3263 APInt api = APInt(32, 0);
3264 initFromAPInt(api.floatToBits(f));
3267 APFloat::APFloat(double d)
3269 APInt api = APInt(64, 0);
3270 initFromAPInt(api.doubleToBits(d));
3274 static void append(SmallVectorImpl<char> &Buffer,
3275 unsigned N, const char *Str) {
3276 unsigned Start = Buffer.size();
3277 Buffer.set_size(Start + N);
3278 memcpy(&Buffer[Start], Str, N);
3281 template <unsigned N>
3282 void append(SmallVectorImpl<char> &Buffer, const char (&Str)[N]) {
3283 append(Buffer, N, Str);
3286 /// Removes data from the given significand until it is no more
3287 /// precise than is required for the desired precision.
3288 void AdjustToPrecision(APInt &significand,
3289 int &exp, unsigned FormatPrecision) {
3290 unsigned bits = significand.getActiveBits();
3292 // 196/59 is a very slight overestimate of lg_2(10).
3293 unsigned bitsRequired = (FormatPrecision * 196 + 58) / 59;
3295 if (bits <= bitsRequired) return;
3297 unsigned tensRemovable = (bits - bitsRequired) * 59 / 196;
3298 if (!tensRemovable) return;
3300 exp += tensRemovable;
3302 APInt divisor(significand.getBitWidth(), 1);
3303 APInt powten(significand.getBitWidth(), 10);
3305 if (tensRemovable & 1)
3307 tensRemovable >>= 1;
3308 if (!tensRemovable) break;
3312 significand = significand.udiv(divisor);
3314 // Truncate the significand down to its active bit count, but
3315 // don't try to drop below 32.
3316 unsigned newPrecision = std::max(32U, significand.getActiveBits());
3317 significand.trunc(newPrecision);
3321 void AdjustToPrecision(SmallVectorImpl<char> &buffer,
3322 int &exp, unsigned FormatPrecision) {
3323 unsigned N = buffer.size();
3324 if (N <= FormatPrecision) return;
3326 // The most significant figures are the last ones in the buffer.
3327 unsigned FirstSignificant = N - FormatPrecision;
3330 // FIXME: this probably shouldn't use 'round half up'.
3332 // Rounding down is just a truncation, except we also want to drop
3333 // trailing zeros from the new result.
3334 if (buffer[FirstSignificant - 1] < '5') {
3335 while (buffer[FirstSignificant] == '0')
3338 exp += FirstSignificant;
3339 buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3343 // Rounding up requires a decimal add-with-carry. If we continue
3344 // the carry, the newly-introduced zeros will just be truncated.
3345 for (unsigned I = FirstSignificant; I != N; ++I) {
3346 if (buffer[I] == '9') {
3354 // If we carried through, we have exactly one digit of precision.
3355 if (FirstSignificant == N) {
3356 exp += FirstSignificant;
3358 buffer.push_back('1');
3362 exp += FirstSignificant;
3363 buffer.erase(&buffer[0], &buffer[FirstSignificant]);
3367 void APFloat::toString(SmallVectorImpl<char> &Str,
3368 unsigned FormatPrecision,
3369 unsigned FormatMaxPadding) const {
3373 return append(Str, "-Inf");
3375 return append(Str, "+Inf");
3377 case fcNaN: return append(Str, "NaN");
3383 if (!FormatMaxPadding)
3384 append(Str, "0.0E+0");
3396 // Decompose the number into an APInt and an exponent.
3397 int exp = exponent - ((int) semantics->precision - 1);
3398 APInt significand(semantics->precision,
3399 partCountForBits(semantics->precision),
3400 significandParts());
3402 // Set FormatPrecision if zero. We want to do this before we
3403 // truncate trailing zeros, as those are part of the precision.
3404 if (!FormatPrecision) {
3405 // It's an interesting question whether to use the nominal
3406 // precision or the active precision here for denormals.
3408 // FormatPrecision = ceil(significandBits / lg_2(10))
3409 FormatPrecision = (semantics->precision * 59 + 195) / 196;
3412 // Ignore trailing binary zeros.
3413 int trailingZeros = significand.countTrailingZeros();
3414 exp += trailingZeros;
3415 significand = significand.lshr(trailingZeros);
3417 // Change the exponent from 2^e to 10^e.
3420 } else if (exp > 0) {
3422 significand.zext(semantics->precision + exp);
3423 significand <<= exp;
3425 } else { /* exp < 0 */
3428 // We transform this using the identity:
3429 // (N)(2^-e) == (N)(5^e)(10^-e)
3430 // This means we have to multiply N (the significand) by 5^e.
3431 // To avoid overflow, we have to operate on numbers large
3432 // enough to store N * 5^e:
3433 // log2(N * 5^e) == log2(N) + e * log2(5)
3434 // <= semantics->precision + e * 137 / 59
3435 // (log_2(5) ~ 2.321928 < 2.322034 ~ 137/59)
3437 unsigned precision = semantics->precision + 137 * texp / 59;
3439 // Multiply significand by 5^e.
3440 // N * 5^0101 == N * 5^(1*1) * 5^(0*2) * 5^(1*4) * 5^(0*8)
3441 significand.zext(precision);
3442 APInt five_to_the_i(precision, 5);
3444 if (texp & 1) significand *= five_to_the_i;
3448 five_to_the_i *= five_to_the_i;
3452 AdjustToPrecision(significand, exp, FormatPrecision);
3454 llvm::SmallVector<char, 256> buffer;
3457 unsigned precision = significand.getBitWidth();
3458 APInt ten(precision, 10);
3459 APInt digit(precision, 0);
3461 bool inTrail = true;
3462 while (significand != 0) {
3463 // digit <- significand % 10
3464 // significand <- significand / 10
3465 APInt::udivrem(significand, ten, significand, digit);
3467 unsigned d = digit.getZExtValue();
3469 // Drop trailing zeros.
3470 if (inTrail && !d) exp++;
3472 buffer.push_back((char) ('0' + d));
3477 assert(!buffer.empty() && "no characters in buffer!");
3479 // Drop down to FormatPrecision.
3480 // TODO: don't do more precise calculations above than are required.
3481 AdjustToPrecision(buffer, exp, FormatPrecision);
3483 unsigned NDigits = buffer.size();
3485 // Check whether we should use scientific notation.
3486 bool FormatScientific;
3487 if (!FormatMaxPadding)
3488 FormatScientific = true;
3493 // But we shouldn't make the number look more precise than it is.
3494 FormatScientific = ((unsigned) exp > FormatMaxPadding ||
3495 NDigits + (unsigned) exp > FormatPrecision);
3497 // Power of the most significant digit.
3498 int MSD = exp + (int) (NDigits - 1);
3501 FormatScientific = false;
3503 // 765e-5 == 0.00765
3505 FormatScientific = ((unsigned) -MSD) > FormatMaxPadding;
3510 // Scientific formatting is pretty straightforward.
3511 if (FormatScientific) {
3512 exp += (NDigits - 1);
3514 Str.push_back(buffer[NDigits-1]);
3519 for (unsigned I = 1; I != NDigits; ++I)
3520 Str.push_back(buffer[NDigits-1-I]);
3523 Str.push_back(exp >= 0 ? '+' : '-');
3524 if (exp < 0) exp = -exp;
3525 SmallVector<char, 6> expbuf;
3527 expbuf.push_back((char) ('0' + (exp % 10)));
3530 for (unsigned I = 0, E = expbuf.size(); I != E; ++I)
3531 Str.push_back(expbuf[E-1-I]);
3535 // Non-scientific, positive exponents.
3537 for (unsigned I = 0; I != NDigits; ++I)
3538 Str.push_back(buffer[NDigits-1-I]);
3539 for (unsigned I = 0; I != (unsigned) exp; ++I)
3544 // Non-scientific, negative exponents.
3546 // The number of digits to the left of the decimal point.
3547 int NWholeDigits = exp + (int) NDigits;
3550 if (NWholeDigits > 0) {
3551 for (; I != (unsigned) NWholeDigits; ++I)
3552 Str.push_back(buffer[NDigits-I-1]);
3555 unsigned NZeros = 1 + (unsigned) -NWholeDigits;
3559 for (unsigned Z = 1; Z != NZeros; ++Z)
3563 for (; I != NDigits; ++I)
3564 Str.push_back(buffer[NDigits-I-1]);