1 //===-- APFloat.cpp - Implement APFloat class -----------------------------===//
3 // The LLVM Compiler Infrastructure
5 // This file was developed by Neil Booth and is distributed under the
6 // University of Illinois Open Source 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 //===----------------------------------------------------------------------===//
16 #include "llvm/ADT/APFloat.h"
20 #define convolve(lhs, rhs) ((lhs) * 4 + (rhs))
22 /* Assumed in hexadecimal significand parsing. */
23 COMPILE_TIME_ASSERT(integerPartWidth % 4 == 0);
27 /* Represents floating point arithmetic semantics. */
29 /* The largest E such that 2^E is representable; this matches the
30 definition of IEEE 754. */
31 exponent_t maxExponent;
33 /* The smallest E such that 2^E is a normalized number; this
34 matches the definition of IEEE 754. */
35 exponent_t minExponent;
37 /* Number of bits in the significand. This includes the integer
39 unsigned char precision;
41 /* If the target format has an implicit integer bit. */
42 bool implicitIntegerBit;
45 const fltSemantics APFloat::IEEEsingle = { 127, -126, 24, true };
46 const fltSemantics APFloat::IEEEdouble = { 1023, -1022, 53, true };
47 const fltSemantics APFloat::IEEEquad = { 16383, -16382, 113, true };
48 const fltSemantics APFloat::x87DoubleExtended = { 16383, -16382, 64, false };
49 const fltSemantics APFloat::Bogus = { 0, 0, 0, false };
52 /* Put a bunch of private, handy routines in an anonymous namespace. */
56 partCountForBits(unsigned int bits)
58 return ((bits) + integerPartWidth - 1) / integerPartWidth;
62 digitValue(unsigned int c)
74 hexDigitValue (unsigned int c)
93 /* This is ugly and needs cleaning up, but I don't immediately see
94 how whilst remaining safe. */
96 totalExponent(const char *p, int exponentAdjustment)
98 integerPart unsignedExponent;
99 bool negative, overflow;
102 /* Move past the exponent letter and sign to the digits. */
104 negative = *p == '-';
105 if(*p == '-' || *p == '+')
108 unsignedExponent = 0;
113 value = digitValue(*p);
118 unsignedExponent = unsignedExponent * 10 + value;
119 if(unsignedExponent > 65535)
123 if(exponentAdjustment > 65535 || exponentAdjustment < -65536)
127 exponent = unsignedExponent;
129 exponent = -exponent;
130 exponent += exponentAdjustment;
131 if(exponent > 65535 || exponent < -65536)
136 exponent = negative ? -65536: 65535;
142 skipLeadingZeroesAndAnyDot(const char *p, const char **dot)
157 /* Return the trailing fraction of a hexadecimal number.
158 DIGITVALUE is the first hex digit of the fraction, P points to
161 trailingHexadecimalFraction(const char *p, unsigned int digitValue)
163 unsigned int hexDigit;
165 /* If the first trailing digit isn't 0 or 8 we can work out the
166 fraction immediately. */
168 return lfMoreThanHalf;
169 else if(digitValue < 8 && digitValue > 0)
170 return lfLessThanHalf;
172 /* Otherwise we need to find the first non-zero digit. */
176 hexDigit = hexDigitValue(*p);
178 /* If we ran off the end it is exactly zero or one-half, otherwise
181 return digitValue == 0 ? lfExactlyZero: lfExactlyHalf;
183 return digitValue == 0 ? lfLessThanHalf: lfMoreThanHalf;
186 /* Return the fraction lost were a bignum truncated. */
188 lostFractionThroughTruncation(integerPart *parts,
189 unsigned int partCount,
194 lsb = APInt::tcLSB(parts, partCount);
196 /* Note this is guaranteed true if bits == 0, or LSB == -1U. */
198 return lfExactlyZero;
200 return lfExactlyHalf;
201 if(bits <= partCount * integerPartWidth
202 && APInt::tcExtractBit(parts, bits - 1))
203 return lfMoreThanHalf;
205 return lfLessThanHalf;
208 /* Shift DST right BITS bits noting lost fraction. */
210 shiftRight(integerPart *dst, unsigned int parts, unsigned int bits)
212 lostFraction lost_fraction;
214 lost_fraction = lostFractionThroughTruncation(dst, parts, bits);
216 APInt::tcShiftRight(dst, parts, bits);
218 return lost_fraction;
224 APFloat::initialize(const fltSemantics *ourSemantics)
228 semantics = ourSemantics;
231 significand.parts = new integerPart[count];
235 APFloat::freeSignificand()
238 delete [] significand.parts;
242 APFloat::assign(const APFloat &rhs)
244 assert(semantics == rhs.semantics);
247 category = rhs.category;
248 exponent = rhs.exponent;
249 if(category == fcNormal)
250 copySignificand(rhs);
254 APFloat::copySignificand(const APFloat &rhs)
256 assert(category == fcNormal);
257 assert(rhs.partCount() >= partCount());
259 APInt::tcAssign(significandParts(), rhs.significandParts(),
264 APFloat::operator=(const APFloat &rhs)
267 if(semantics != rhs.semantics) {
269 initialize(rhs.semantics);
278 APFloat::operator==(const APFloat &rhs) const {
281 if (semantics != rhs.semantics ||
282 category != rhs.category)
284 if (category==fcQNaN)
286 else if (category==fcZero || category==fcInfinity)
287 return sign==rhs.sign;
289 if (sign!=rhs.sign || exponent!=rhs.exponent)
292 const integerPart* p=significandParts();
293 const integerPart* q=rhs.significandParts();
294 for (; i>0; i--, p++, q++) {
302 APFloat::APFloat(const fltSemantics &ourSemantics, integerPart value)
304 initialize(&ourSemantics);
307 exponent = ourSemantics.precision - 1;
308 significandParts()[0] = value;
309 normalize(rmNearestTiesToEven, lfExactlyZero);
312 APFloat::APFloat(const fltSemantics &ourSemantics,
313 fltCategory ourCategory, bool negative)
315 initialize(&ourSemantics);
316 category = ourCategory;
318 if(category == fcNormal)
322 APFloat::APFloat(const fltSemantics &ourSemantics, const char *text)
324 initialize(&ourSemantics);
325 convertFromString(text, rmNearestTiesToEven);
328 APFloat::APFloat(const APFloat &rhs)
330 initialize(rhs.semantics);
340 APFloat::partCount() const
342 return partCountForBits(semantics->precision + 1);
346 APFloat::semanticsPrecision(const fltSemantics &semantics)
348 return semantics.precision;
352 APFloat::significandParts() const
354 return const_cast<APFloat *>(this)->significandParts();
358 APFloat::significandParts()
360 assert(category == fcNormal);
363 return significand.parts;
365 return &significand.part;
368 /* Combine the effect of two lost fractions. */
370 APFloat::combineLostFractions(lostFraction moreSignificant,
371 lostFraction lessSignificant)
373 if(lessSignificant != lfExactlyZero) {
374 if(moreSignificant == lfExactlyZero)
375 moreSignificant = lfLessThanHalf;
376 else if(moreSignificant == lfExactlyHalf)
377 moreSignificant = lfMoreThanHalf;
380 return moreSignificant;
384 APFloat::zeroSignificand()
387 APInt::tcSet(significandParts(), 0, partCount());
390 /* Increment an fcNormal floating point number's significand. */
392 APFloat::incrementSignificand()
396 carry = APInt::tcIncrement(significandParts(), partCount());
398 /* Our callers should never cause us to overflow. */
402 /* Add the significand of the RHS. Returns the carry flag. */
404 APFloat::addSignificand(const APFloat &rhs)
408 parts = significandParts();
410 assert(semantics == rhs.semantics);
411 assert(exponent == rhs.exponent);
413 return APInt::tcAdd(parts, rhs.significandParts(), 0, partCount());
416 /* Subtract the significand of the RHS with a borrow flag. Returns
419 APFloat::subtractSignificand(const APFloat &rhs, integerPart borrow)
423 parts = significandParts();
425 assert(semantics == rhs.semantics);
426 assert(exponent == rhs.exponent);
428 return APInt::tcSubtract(parts, rhs.significandParts(), borrow,
432 /* Multiply the significand of the RHS. If ADDEND is non-NULL, add it
433 on to the full-precision result of the multiplication. Returns the
436 APFloat::multiplySignificand(const APFloat &rhs, const APFloat *addend)
438 unsigned int omsb; // One, not zero, based MSB.
439 unsigned int partsCount, newPartsCount, precision;
440 integerPart *lhsSignificand;
441 integerPart scratch[4];
442 integerPart *fullSignificand;
443 lostFraction lost_fraction;
445 assert(semantics == rhs.semantics);
447 precision = semantics->precision;
448 newPartsCount = partCountForBits(precision * 2);
450 if(newPartsCount > 4)
451 fullSignificand = new integerPart[newPartsCount];
453 fullSignificand = scratch;
455 lhsSignificand = significandParts();
456 partsCount = partCount();
458 APInt::tcFullMultiply(fullSignificand, lhsSignificand,
459 rhs.significandParts(), partsCount);
461 lost_fraction = lfExactlyZero;
462 omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
463 exponent += rhs.exponent;
466 Significand savedSignificand = significand;
467 const fltSemantics *savedSemantics = semantics;
468 fltSemantics extendedSemantics;
470 unsigned int extendedPrecision;
472 /* Normalize our MSB. */
473 extendedPrecision = precision + precision - 1;
474 if(omsb != extendedPrecision)
476 APInt::tcShiftLeft(fullSignificand, newPartsCount,
477 extendedPrecision - omsb);
478 exponent -= extendedPrecision - omsb;
481 /* Create new semantics. */
482 extendedSemantics = *semantics;
483 extendedSemantics.precision = extendedPrecision;
485 if(newPartsCount == 1)
486 significand.part = fullSignificand[0];
488 significand.parts = fullSignificand;
489 semantics = &extendedSemantics;
491 APFloat extendedAddend(*addend);
492 status = extendedAddend.convert(extendedSemantics, rmTowardZero);
493 assert(status == opOK);
494 lost_fraction = addOrSubtractSignificand(extendedAddend, false);
496 /* Restore our state. */
497 if(newPartsCount == 1)
498 fullSignificand[0] = significand.part;
499 significand = savedSignificand;
500 semantics = savedSemantics;
502 omsb = APInt::tcMSB(fullSignificand, newPartsCount) + 1;
505 exponent -= (precision - 1);
507 if(omsb > precision) {
508 unsigned int bits, significantParts;
511 bits = omsb - precision;
512 significantParts = partCountForBits(omsb);
513 lf = shiftRight(fullSignificand, significantParts, bits);
514 lost_fraction = combineLostFractions(lf, lost_fraction);
518 APInt::tcAssign(lhsSignificand, fullSignificand, partsCount);
520 if(newPartsCount > 4)
521 delete [] fullSignificand;
523 return lost_fraction;
526 /* Multiply the significands of LHS and RHS to DST. */
528 APFloat::divideSignificand(const APFloat &rhs)
530 unsigned int bit, i, partsCount;
531 const integerPart *rhsSignificand;
532 integerPart *lhsSignificand, *dividend, *divisor;
533 integerPart scratch[4];
534 lostFraction lost_fraction;
536 assert(semantics == rhs.semantics);
538 lhsSignificand = significandParts();
539 rhsSignificand = rhs.significandParts();
540 partsCount = partCount();
543 dividend = new integerPart[partsCount * 2];
547 divisor = dividend + partsCount;
549 /* Copy the dividend and divisor as they will be modified in-place. */
550 for(i = 0; i < partsCount; i++) {
551 dividend[i] = lhsSignificand[i];
552 divisor[i] = rhsSignificand[i];
553 lhsSignificand[i] = 0;
556 exponent -= rhs.exponent;
558 unsigned int precision = semantics->precision;
560 /* Normalize the divisor. */
561 bit = precision - APInt::tcMSB(divisor, partsCount) - 1;
564 APInt::tcShiftLeft(divisor, partsCount, bit);
567 /* Normalize the dividend. */
568 bit = precision - APInt::tcMSB(dividend, partsCount) - 1;
571 APInt::tcShiftLeft(dividend, partsCount, bit);
574 if(APInt::tcCompare(dividend, divisor, partsCount) < 0) {
576 APInt::tcShiftLeft(dividend, partsCount, 1);
577 assert(APInt::tcCompare(dividend, divisor, partsCount) >= 0);
581 for(bit = precision; bit; bit -= 1) {
582 if(APInt::tcCompare(dividend, divisor, partsCount) >= 0) {
583 APInt::tcSubtract(dividend, divisor, 0, partsCount);
584 APInt::tcSetBit(lhsSignificand, bit - 1);
587 APInt::tcShiftLeft(dividend, partsCount, 1);
590 /* Figure out the lost fraction. */
591 int cmp = APInt::tcCompare(dividend, divisor, partsCount);
594 lost_fraction = lfMoreThanHalf;
596 lost_fraction = lfExactlyHalf;
597 else if(APInt::tcIsZero(dividend, partsCount))
598 lost_fraction = lfExactlyZero;
600 lost_fraction = lfLessThanHalf;
605 return lost_fraction;
609 APFloat::significandMSB() const
611 return APInt::tcMSB(significandParts(), partCount());
615 APFloat::significandLSB() const
617 return APInt::tcLSB(significandParts(), partCount());
620 /* Note that a zero result is NOT normalized to fcZero. */
622 APFloat::shiftSignificandRight(unsigned int bits)
624 /* Our exponent should not overflow. */
625 assert((exponent_t) (exponent + bits) >= exponent);
629 return shiftRight(significandParts(), partCount(), bits);
632 /* Shift the significand left BITS bits, subtract BITS from its exponent. */
634 APFloat::shiftSignificandLeft(unsigned int bits)
636 assert(bits < semantics->precision);
639 unsigned int partsCount = partCount();
641 APInt::tcShiftLeft(significandParts(), partsCount, bits);
644 assert(!APInt::tcIsZero(significandParts(), partsCount));
649 APFloat::compareAbsoluteValue(const APFloat &rhs) const
653 assert(semantics == rhs.semantics);
654 assert(category == fcNormal);
655 assert(rhs.category == fcNormal);
657 compare = exponent - rhs.exponent;
659 /* If exponents are equal, do an unsigned bignum comparison of the
662 compare = APInt::tcCompare(significandParts(), rhs.significandParts(),
666 return cmpGreaterThan;
673 /* Handle overflow. Sign is preserved. We either become infinity or
674 the largest finite number. */
676 APFloat::handleOverflow(roundingMode rounding_mode)
679 if(rounding_mode == rmNearestTiesToEven
680 || rounding_mode == rmNearestTiesToAway
681 || (rounding_mode == rmTowardPositive && !sign)
682 || (rounding_mode == rmTowardNegative && sign))
684 category = fcInfinity;
685 return (opStatus) (opOverflow | opInexact);
688 /* Otherwise we become the largest finite number. */
690 exponent = semantics->maxExponent;
691 APInt::tcSetLeastSignificantBits(significandParts(), partCount(),
692 semantics->precision);
697 /* This routine must work for fcZero of both signs, and fcNormal
700 APFloat::roundAwayFromZero(roundingMode rounding_mode,
701 lostFraction lost_fraction)
703 /* QNaNs and infinities should not have lost fractions. */
704 assert(category == fcNormal || category == fcZero);
706 /* Our caller has already handled this case. */
707 assert(lost_fraction != lfExactlyZero);
709 switch(rounding_mode) {
713 case rmNearestTiesToAway:
714 return lost_fraction == lfExactlyHalf || lost_fraction == lfMoreThanHalf;
716 case rmNearestTiesToEven:
717 if(lost_fraction == lfMoreThanHalf)
720 /* Our zeroes don't have a significand to test. */
721 if(lost_fraction == lfExactlyHalf && category != fcZero)
722 return significandParts()[0] & 1;
729 case rmTowardPositive:
730 return sign == false;
732 case rmTowardNegative:
738 APFloat::normalize(roundingMode rounding_mode,
739 lostFraction lost_fraction)
741 unsigned int omsb; /* One, not zero, based MSB. */
744 if(category != fcNormal)
747 /* Before rounding normalize the exponent of fcNormal numbers. */
748 omsb = significandMSB() + 1;
751 /* OMSB is numbered from 1. We want to place it in the integer
752 bit numbered PRECISON if possible, with a compensating change in
754 exponentChange = omsb - semantics->precision;
756 /* If the resulting exponent is too high, overflow according to
757 the rounding mode. */
758 if(exponent + exponentChange > semantics->maxExponent)
759 return handleOverflow(rounding_mode);
761 /* Subnormal numbers have exponent minExponent, and their MSB
762 is forced based on that. */
763 if(exponent + exponentChange < semantics->minExponent)
764 exponentChange = semantics->minExponent - exponent;
766 /* Shifting left is easy as we don't lose precision. */
767 if(exponentChange < 0) {
768 assert(lost_fraction == lfExactlyZero);
770 shiftSignificandLeft(-exponentChange);
775 if(exponentChange > 0) {
778 /* Shift right and capture any new lost fraction. */
779 lf = shiftSignificandRight(exponentChange);
781 lost_fraction = combineLostFractions(lf, lost_fraction);
783 /* Keep OMSB up-to-date. */
784 if(omsb > (unsigned) exponentChange)
785 omsb -= (unsigned) exponentChange;
791 /* Now round the number according to rounding_mode given the lost
794 /* As specified in IEEE 754, since we do not trap we do not report
795 underflow for exact results. */
796 if(lost_fraction == lfExactlyZero) {
797 /* Canonicalize zeroes. */
804 /* Increment the significand if we're rounding away from zero. */
805 if(roundAwayFromZero(rounding_mode, lost_fraction)) {
807 exponent = semantics->minExponent;
809 incrementSignificand();
810 omsb = significandMSB() + 1;
812 /* Did the significand increment overflow? */
813 if(omsb == (unsigned) semantics->precision + 1) {
814 /* Renormalize by incrementing the exponent and shifting our
815 significand right one. However if we already have the
816 maximum exponent we overflow to infinity. */
817 if(exponent == semantics->maxExponent) {
818 category = fcInfinity;
820 return (opStatus) (opOverflow | opInexact);
823 shiftSignificandRight(1);
829 /* The normal case - we were and are not denormal, and any
830 significand increment above didn't overflow. */
831 if(omsb == semantics->precision)
834 /* We have a non-zero denormal. */
835 assert(omsb < semantics->precision);
836 assert(exponent == semantics->minExponent);
838 /* Canonicalize zeroes. */
842 /* The fcZero case is a denormal that underflowed to zero. */
843 return (opStatus) (opUnderflow | opInexact);
847 APFloat::addOrSubtractSpecials(const APFloat &rhs, bool subtract)
849 switch(convolve(category, rhs.category)) {
853 case convolve(fcQNaN, fcZero):
854 case convolve(fcQNaN, fcNormal):
855 case convolve(fcQNaN, fcInfinity):
856 case convolve(fcQNaN, fcQNaN):
857 case convolve(fcNormal, fcZero):
858 case convolve(fcInfinity, fcNormal):
859 case convolve(fcInfinity, fcZero):
862 case convolve(fcZero, fcQNaN):
863 case convolve(fcNormal, fcQNaN):
864 case convolve(fcInfinity, fcQNaN):
868 case convolve(fcNormal, fcInfinity):
869 case convolve(fcZero, fcInfinity):
870 category = fcInfinity;
871 sign = rhs.sign ^ subtract;
874 case convolve(fcZero, fcNormal):
876 sign = rhs.sign ^ subtract;
879 case convolve(fcZero, fcZero):
880 /* Sign depends on rounding mode; handled by caller. */
883 case convolve(fcInfinity, fcInfinity):
884 /* Differently signed infinities can only be validly
886 if(sign ^ rhs.sign != subtract) {
893 case convolve(fcNormal, fcNormal):
898 /* Add or subtract two normal numbers. */
900 APFloat::addOrSubtractSignificand(const APFloat &rhs, bool subtract)
903 lostFraction lost_fraction;
906 /* Determine if the operation on the absolute values is effectively
907 an addition or subtraction. */
908 subtract ^= (sign ^ rhs.sign);
910 /* Are we bigger exponent-wise than the RHS? */
911 bits = exponent - rhs.exponent;
913 /* Subtraction is more subtle than one might naively expect. */
915 APFloat temp_rhs(rhs);
919 reverse = compareAbsoluteValue(temp_rhs) == cmpLessThan;
920 lost_fraction = lfExactlyZero;
921 } else if(bits > 0) {
922 lost_fraction = temp_rhs.shiftSignificandRight(bits - 1);
923 shiftSignificandLeft(1);
925 } else if(bits < 0) {
926 lost_fraction = shiftSignificandRight(-bits - 1);
927 temp_rhs.shiftSignificandLeft(1);
932 carry = temp_rhs.subtractSignificand
933 (*this, lost_fraction != lfExactlyZero);
934 copySignificand(temp_rhs);
937 carry = subtractSignificand
938 (temp_rhs, lost_fraction != lfExactlyZero);
941 /* Invert the lost fraction - it was on the RHS and
943 if(lost_fraction == lfLessThanHalf)
944 lost_fraction = lfMoreThanHalf;
945 else if(lost_fraction == lfMoreThanHalf)
946 lost_fraction = lfLessThanHalf;
948 /* The code above is intended to ensure that no borrow is
953 APFloat temp_rhs(rhs);
955 lost_fraction = temp_rhs.shiftSignificandRight(bits);
956 carry = addSignificand(temp_rhs);
958 lost_fraction = shiftSignificandRight(-bits);
959 carry = addSignificand(rhs);
962 /* We have a guard bit; generating a carry cannot happen. */
966 return lost_fraction;
970 APFloat::multiplySpecials(const APFloat &rhs)
972 switch(convolve(category, rhs.category)) {
976 case convolve(fcQNaN, fcZero):
977 case convolve(fcQNaN, fcNormal):
978 case convolve(fcQNaN, fcInfinity):
979 case convolve(fcQNaN, fcQNaN):
980 case convolve(fcZero, fcQNaN):
981 case convolve(fcNormal, fcQNaN):
982 case convolve(fcInfinity, fcQNaN):
986 case convolve(fcNormal, fcInfinity):
987 case convolve(fcInfinity, fcNormal):
988 case convolve(fcInfinity, fcInfinity):
989 category = fcInfinity;
992 case convolve(fcZero, fcNormal):
993 case convolve(fcNormal, fcZero):
994 case convolve(fcZero, fcZero):
998 case convolve(fcZero, fcInfinity):
999 case convolve(fcInfinity, fcZero):
1003 case convolve(fcNormal, fcNormal):
1009 APFloat::divideSpecials(const APFloat &rhs)
1011 switch(convolve(category, rhs.category)) {
1015 case convolve(fcQNaN, fcZero):
1016 case convolve(fcQNaN, fcNormal):
1017 case convolve(fcQNaN, fcInfinity):
1018 case convolve(fcQNaN, fcQNaN):
1019 case convolve(fcInfinity, fcZero):
1020 case convolve(fcInfinity, fcNormal):
1021 case convolve(fcZero, fcInfinity):
1022 case convolve(fcZero, fcNormal):
1025 case convolve(fcZero, fcQNaN):
1026 case convolve(fcNormal, fcQNaN):
1027 case convolve(fcInfinity, fcQNaN):
1031 case convolve(fcNormal, fcInfinity):
1035 case convolve(fcNormal, fcZero):
1036 category = fcInfinity;
1039 case convolve(fcInfinity, fcInfinity):
1040 case convolve(fcZero, fcZero):
1044 case convolve(fcNormal, fcNormal):
1051 APFloat::changeSign()
1053 /* Look mummy, this one's easy. */
1057 /* Normalized addition or subtraction. */
1059 APFloat::addOrSubtract(const APFloat &rhs, roundingMode rounding_mode,
1064 fs = addOrSubtractSpecials(rhs, subtract);
1066 /* This return code means it was not a simple case. */
1067 if(fs == opDivByZero) {
1068 lostFraction lost_fraction;
1070 lost_fraction = addOrSubtractSignificand(rhs, subtract);
1071 fs = normalize(rounding_mode, lost_fraction);
1073 /* Can only be zero if we lost no fraction. */
1074 assert(category != fcZero || lost_fraction == lfExactlyZero);
1077 /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1078 positive zero unless rounding to minus infinity, except that
1079 adding two like-signed zeroes gives that zero. */
1080 if(category == fcZero) {
1081 if(rhs.category != fcZero || (sign == rhs.sign) == subtract)
1082 sign = (rounding_mode == rmTowardNegative);
1088 /* Normalized addition. */
1090 APFloat::add(const APFloat &rhs, roundingMode rounding_mode)
1092 return addOrSubtract(rhs, rounding_mode, false);
1095 /* Normalized subtraction. */
1097 APFloat::subtract(const APFloat &rhs, roundingMode rounding_mode)
1099 return addOrSubtract(rhs, rounding_mode, true);
1102 /* Normalized multiply. */
1104 APFloat::multiply(const APFloat &rhs, roundingMode rounding_mode)
1109 fs = multiplySpecials(rhs);
1111 if(category == fcNormal) {
1112 lostFraction lost_fraction = multiplySignificand(rhs, 0);
1113 fs = normalize(rounding_mode, lost_fraction);
1114 if(lost_fraction != lfExactlyZero)
1115 fs = (opStatus) (fs | opInexact);
1121 /* Normalized divide. */
1123 APFloat::divide(const APFloat &rhs, roundingMode rounding_mode)
1128 fs = divideSpecials(rhs);
1130 if(category == fcNormal) {
1131 lostFraction lost_fraction = divideSignificand(rhs);
1132 fs = normalize(rounding_mode, lost_fraction);
1133 if(lost_fraction != lfExactlyZero)
1134 fs = (opStatus) (fs | opInexact);
1140 /* Normalized fused-multiply-add. */
1142 APFloat::fusedMultiplyAdd(const APFloat &multiplicand,
1143 const APFloat &addend,
1144 roundingMode rounding_mode)
1148 /* Post-multiplication sign, before addition. */
1149 sign ^= multiplicand.sign;
1151 /* If and only if all arguments are normal do we need to do an
1152 extended-precision calculation. */
1153 if(category == fcNormal
1154 && multiplicand.category == fcNormal
1155 && addend.category == fcNormal) {
1156 lostFraction lost_fraction;
1158 lost_fraction = multiplySignificand(multiplicand, &addend);
1159 fs = normalize(rounding_mode, lost_fraction);
1160 if(lost_fraction != lfExactlyZero)
1161 fs = (opStatus) (fs | opInexact);
1163 /* If two numbers add (exactly) to zero, IEEE 754 decrees it is a
1164 positive zero unless rounding to minus infinity, except that
1165 adding two like-signed zeroes gives that zero. */
1166 if(category == fcZero && sign != addend.sign)
1167 sign = (rounding_mode == rmTowardNegative);
1169 fs = multiplySpecials(multiplicand);
1171 /* FS can only be opOK or opInvalidOp. There is no more work
1172 to do in the latter case. The IEEE-754R standard says it is
1173 implementation-defined in this case whether, if ADDEND is a
1174 quiet QNaN, we raise invalid op; this implementation does so.
1176 If we need to do the addition we can do so with normal
1179 fs = addOrSubtract(addend, rounding_mode, false);
1185 /* Comparison requires normalized numbers. */
1187 APFloat::compare(const APFloat &rhs) const
1191 assert(semantics == rhs.semantics);
1193 switch(convolve(category, rhs.category)) {
1197 case convolve(fcQNaN, fcZero):
1198 case convolve(fcQNaN, fcNormal):
1199 case convolve(fcQNaN, fcInfinity):
1200 case convolve(fcQNaN, fcQNaN):
1201 case convolve(fcZero, fcQNaN):
1202 case convolve(fcNormal, fcQNaN):
1203 case convolve(fcInfinity, fcQNaN):
1204 return cmpUnordered;
1206 case convolve(fcInfinity, fcNormal):
1207 case convolve(fcInfinity, fcZero):
1208 case convolve(fcNormal, fcZero):
1212 return cmpGreaterThan;
1214 case convolve(fcNormal, fcInfinity):
1215 case convolve(fcZero, fcInfinity):
1216 case convolve(fcZero, fcNormal):
1218 return cmpGreaterThan;
1222 case convolve(fcInfinity, fcInfinity):
1223 if(sign == rhs.sign)
1228 return cmpGreaterThan;
1230 case convolve(fcZero, fcZero):
1233 case convolve(fcNormal, fcNormal):
1237 /* Two normal numbers. Do they have the same sign? */
1238 if(sign != rhs.sign) {
1240 result = cmpLessThan;
1242 result = cmpGreaterThan;
1244 /* Compare absolute values; invert result if negative. */
1245 result = compareAbsoluteValue(rhs);
1248 if(result == cmpLessThan)
1249 result = cmpGreaterThan;
1250 else if(result == cmpGreaterThan)
1251 result = cmpLessThan;
1259 APFloat::convert(const fltSemantics &toSemantics,
1260 roundingMode rounding_mode)
1262 unsigned int newPartCount;
1265 newPartCount = partCountForBits(toSemantics.precision + 1);
1267 /* If our new form is wider, re-allocate our bit pattern into wider
1269 if(newPartCount > partCount()) {
1270 integerPart *newParts;
1272 newParts = new integerPart[newPartCount];
1273 APInt::tcSet(newParts, 0, newPartCount);
1274 APInt::tcAssign(newParts, significandParts(), partCount());
1276 significand.parts = newParts;
1279 if(category == fcNormal) {
1280 /* Re-interpret our bit-pattern. */
1281 exponent += toSemantics.precision - semantics->precision;
1282 semantics = &toSemantics;
1283 fs = normalize(rounding_mode, lfExactlyZero);
1285 semantics = &toSemantics;
1292 /* Convert a floating point number to an integer according to the
1293 rounding mode. If the rounded integer value is out of range this
1294 returns an invalid operation exception. If the rounded value is in
1295 range but the floating point number is not the exact integer, the C
1296 standard doesn't require an inexact exception to be raised. IEEE
1297 854 does require it so we do that.
1299 Note that for conversions to integer type the C standard requires
1300 round-to-zero to always be used. */
1302 APFloat::convertToInteger(integerPart *parts, unsigned int width,
1304 roundingMode rounding_mode) const
1306 lostFraction lost_fraction;
1307 unsigned int msb, partsCount;
1310 /* Handle the three special cases first. */
1311 if(category == fcInfinity || category == fcQNaN)
1314 partsCount = partCountForBits(width);
1316 if(category == fcZero) {
1317 APInt::tcSet(parts, 0, partsCount);
1321 /* Shift the bit pattern so the fraction is lost. */
1324 bits = (int) semantics->precision - 1 - exponent;
1327 lost_fraction = tmp.shiftSignificandRight(bits);
1329 tmp.shiftSignificandLeft(-bits);
1330 lost_fraction = lfExactlyZero;
1333 if(lost_fraction != lfExactlyZero
1334 && tmp.roundAwayFromZero(rounding_mode, lost_fraction))
1335 tmp.incrementSignificand();
1337 msb = tmp.significandMSB();
1339 /* Negative numbers cannot be represented as unsigned. */
1340 if(!isSigned && tmp.sign && msb != -1U)
1343 /* It takes exponent + 1 bits to represent the truncated floating
1344 point number without its sign. We lose a bit for the sign, but
1345 the maximally negative integer is a special case. */
1346 if(msb + 1 > width) /* !! Not same as msb >= width !! */
1349 if(isSigned && msb + 1 == width
1350 && (!tmp.sign || tmp.significandLSB() != msb))
1353 APInt::tcAssign(parts, tmp.significandParts(), partsCount);
1356 APInt::tcNegate(parts, partsCount);
1358 if(lost_fraction == lfExactlyZero)
1365 APFloat::convertFromUnsignedInteger(integerPart *parts,
1366 unsigned int partCount,
1367 roundingMode rounding_mode)
1369 unsigned int msb, precision;
1370 lostFraction lost_fraction;
1372 msb = APInt::tcMSB(parts, partCount) + 1;
1373 precision = semantics->precision;
1375 category = fcNormal;
1376 exponent = precision - 1;
1378 if(msb > precision) {
1379 exponent += (msb - precision);
1380 lost_fraction = shiftRight(parts, partCount, msb - precision);
1383 lost_fraction = lfExactlyZero;
1385 /* Copy the bit image. */
1387 APInt::tcAssign(significandParts(), parts, partCountForBits(msb));
1389 return normalize(rounding_mode, lost_fraction);
1393 APFloat::convertFromInteger(const integerPart *parts,
1394 unsigned int partCount, bool isSigned,
1395 roundingMode rounding_mode)
1401 copy = new integerPart[partCount];
1402 APInt::tcAssign(copy, parts, partCount);
1404 width = partCount * integerPartWidth;
1407 if(isSigned && APInt::tcExtractBit(parts, width - 1)) {
1409 APInt::tcNegate(copy, partCount);
1412 status = convertFromUnsignedInteger(copy, partCount, rounding_mode);
1419 APFloat::convertFromHexadecimalString(const char *p,
1420 roundingMode rounding_mode)
1422 lostFraction lost_fraction;
1423 integerPart *significand;
1424 unsigned int bitPos, partsCount;
1425 const char *dot, *firstSignificantDigit;
1429 category = fcNormal;
1431 significand = significandParts();
1432 partsCount = partCount();
1433 bitPos = partsCount * integerPartWidth;
1435 /* Skip leading zeroes and any(hexa)decimal point. */
1436 p = skipLeadingZeroesAndAnyDot(p, &dot);
1437 firstSignificantDigit = p;
1440 integerPart hex_value;
1447 hex_value = hexDigitValue(*p);
1448 if(hex_value == -1U) {
1449 lost_fraction = lfExactlyZero;
1455 /* Store the number whilst 4-bit nibbles remain. */
1458 hex_value <<= bitPos % integerPartWidth;
1459 significand[bitPos / integerPartWidth] |= hex_value;
1461 lost_fraction = trailingHexadecimalFraction(p, hex_value);
1462 while(hexDigitValue(*p) != -1U)
1468 /* Hex floats require an exponent but not a hexadecimal point. */
1469 assert(*p == 'p' || *p == 'P');
1471 /* Ignore the exponent if we are zero. */
1472 if(p != firstSignificantDigit) {
1475 /* Implicit hexadecimal point? */
1479 /* Calculate the exponent adjustment implicit in the number of
1480 significant digits. */
1481 expAdjustment = dot - firstSignificantDigit;
1482 if(expAdjustment < 0)
1484 expAdjustment = expAdjustment * 4 - 1;
1486 /* Adjust for writing the significand starting at the most
1487 significant nibble. */
1488 expAdjustment += semantics->precision;
1489 expAdjustment -= partsCount * integerPartWidth;
1491 /* Adjust for the given exponent. */
1492 exponent = totalExponent(p, expAdjustment);
1495 return normalize(rounding_mode, lost_fraction);
1499 APFloat::convertFromString(const char *p, roundingMode rounding_mode)
1501 /* Handle a leading minus sign. */
1507 if(p[0] == '0' && (p[1] == 'x' || p[1] == 'X'))
1508 return convertFromHexadecimalString(p + 2, rounding_mode);
1511 assert(0 && "Decimal to binary conversions not yet implemented");
1516 // For good performance it is desirable for different APFloats
1517 // to produce different integers.
1519 APFloat::getHashValue() const {
1520 if (category==fcZero) return sign<<8 | semantics->precision ;
1521 else if (category==fcInfinity) return sign<<9 | semantics->precision;
1522 else if (category==fcQNaN) return 1<<10 | semantics->precision;
1524 uint32_t hash = sign<<11 | semantics->precision | exponent<<12;
1525 const integerPart* p = significandParts();
1526 for (int i=partCount(); i>0; i--, p++)
1527 hash ^= ((uint32_t)*p) ^ (*p)>>32;
1532 // Conversion from APFloat to/from host float/double. It may eventually be
1533 // possible to eliminate these and have everybody deal with APFloats, but that
1534 // will take a while. This approach will not easily extend to long double.
1535 // Current implementation requires partCount()==1, which is correct at the
1536 // moment but could be made more general.
1539 APFloat::convertToDouble() const {
1544 assert(semantics == (const llvm::fltSemantics* const)&IEEEdouble);
1545 assert (partCount()==1);
1547 uint64_t myexponent, mysign, mysignificand;
1549 if (category==fcNormal) {
1551 mysignificand = *significandParts();
1552 myexponent = exponent+1023; //bias
1553 } else if (category==fcZero) {
1557 } else if (category==fcInfinity) {
1561 } else if (category==fcQNaN) {
1564 mysignificand = 0xfffffffffffffLL;
1568 u.i = ((mysign & 1) << 63) | ((myexponent & 0x7ff) << 52) |
1569 (mysignificand & 0xfffffffffffffLL);
1574 APFloat::convertToFloat() const {
1579 assert(semantics == (const llvm::fltSemantics* const)&IEEEsingle);
1580 assert (partCount()==1);
1582 uint32_t mysign, myexponent, mysignificand;
1584 if (category==fcNormal) {
1586 myexponent = exponent+127; //bias
1587 mysignificand = *significandParts();
1588 } else if (category==fcZero) {
1592 } else if (category==fcInfinity) {
1596 } else if (category==fcQNaN) {
1599 mysignificand = 0x7fffff;
1603 u.i = ((mysign&1) << 31) | ((myexponent&0xff) << 23) |
1604 ((mysignificand & 0x7fffff));
1608 APFloat::APFloat(double d) {
1609 initialize(&APFloat::IEEEdouble);
1615 assert(partCount()==1);
1617 uint64_t mysign, myexponent, mysignificand;
1620 myexponent = (u.i >> 52) & 0x7ff;
1621 mysignificand = u.i & 0xfffffffffffffLL;
1623 if (myexponent==0 && mysignificand==0) {
1624 // exponent, significand meaningless
1627 } else if (myexponent==0x7ff && mysignificand==0) {
1628 // exponent, significand meaningless
1629 category = fcInfinity;
1631 } else if (myexponent==0x7ff && (mysignificand & 0x8000000000000LL)) {
1632 // sign, exponent, significand meaningless
1636 category = fcNormal;
1637 exponent = myexponent - 1023;
1638 *significandParts() = mysignificand | 0x100000000000000LL;
1642 APFloat::APFloat(float f) {
1643 initialize(&APFloat::IEEEsingle);
1649 assert(partCount()==1);
1651 uint32_t mysign, myexponent, mysignificand;
1654 myexponent = (u.i >> 23) & 0xff;
1655 mysignificand = u.i & 0x7fffff;
1657 if (myexponent==0 && mysignificand==0) {
1658 // exponent, significand meaningless
1661 } else if (myexponent==0xff && mysignificand==0) {
1662 // exponent, significand meaningless
1663 category = fcInfinity;
1665 } else if (myexponent==0xff && (mysignificand & 0x400000)) {
1666 // sign, exponent, significand meaningless
1669 category = fcNormal;
1671 exponent = myexponent - 127; //bias
1672 *significandParts() = mysignificand | 0x800000; // integer bit