diff --git a/flang/lib/evaluate/integer.h b/flang/lib/evaluate/integer.h index 6f0dfb6e55a7..0c97be3dc17d 100644 --- a/flang/lib/evaluate/integer.h +++ b/flang/lib/evaluate/integer.h @@ -87,6 +87,7 @@ private: static constexpr Part topPartMask{static_cast(~0) >> extraTopPartBits}; public: + // Some types used for member function results struct ValueWithOverflow { Integer value; bool overflow; @@ -994,7 +995,7 @@ private: } } - Part part_[parts]; + Part part_[parts]{}; }; extern template class Integer<8>; diff --git a/flang/lib/evaluate/real.cc b/flang/lib/evaluate/real.cc index 9eef5312c897..3c5778eeb0ab 100644 --- a/flang/lib/evaluate/real.cc +++ b/flang/lib/evaluate/real.cc @@ -481,9 +481,160 @@ std::string Real::DumpHexadecimal() const { } } +template +std::ostream &Real::AsFortran(std::ostream &o, int kind) const { + ValueWithRealFlags scaled{AsScaledDecimal()}; + if (scaled.flags.test(RealFlag::InvalidArgument)) { + o << "(0._" << kind << "/0.)"; + } else if (scaled.flags.test(RealFlag::Overflow)) { + if (IsNegative()) { + o << "(-1._" << kind << "/0.)"; + } else { + o << "(1._" << kind << "/0.)"; + } + } else { + if (scaled.value.negative) { + o << "(-"; + } + std::string digits{scaled.value.integer.UnsignedDecimal()}; + int exponent = scaled.value.decimalExponent + digits.size() - 1; + o << digits[0] << '.' << digits.substr(1); + if (exponent != 0) { + o << 'e' << exponent; + } + o << '_' << kind; + if (scaled.value.negative) { + o << ')'; + } + } + return o; +} + +template +auto Real::AsScaledDecimal() const + -> ValueWithRealFlags { + + // This constant is max x such that REAL(INT(x)) == x without loss + static constexpr Real maxExactSignedValue{ + Word{exponentBias + bits - 2} + .SHIFTL(significandBits) + .IOR(Word::MASKR(significandBits))}; + // This constant is min x such that x + 1.0 == NEXTAFTER(x) + // (or, equivalently, for all representable y >= x, y == AINT(y)) + static constexpr Real minFractionlessValue{ + Word{exponentBias + precision - 1} + .SHIFTL(significandBits) + .IOR(Word::MASKR(significandBits))}; + + ValueWithRealFlags result; + if (IsNotANumber()) { + result.flags.set(RealFlag::InvalidArgument); + return result; + } + if (IsInfinite()) { + result.flags.set(RealFlag::Overflow); + result.value.integer = Word::HUGE(); + return result; + } + if (IsNegative()) { + result = Negate().AsScaledDecimal(); + result.value.negative = true; + return result; + } + if (IsZero()) { + return result; + } + + // N.B. This code could be made more generic and work with output bases + // other than ten, so long as "decimalDigits" were modified accordingly. + Real ten{FromInteger(Word{10}).value}; // 10.0 + + // The maximum exact power of ten can't be calculated as + // a constexpr, so it's computed on the first call to + // this function. + static Real maxExactPowerOfTen; // 10.0 ** decimalDigits + if (maxExactPowerOfTen.IsZero()) { + maxExactPowerOfTen = ten; + for (int j{1}; j < decimalDigits; ++j) { + auto next{maxExactPowerOfTen.Multiply(ten)}; + CHECK(!next.value.IsZero()); + maxExactPowerOfTen = next.value; + } + } + + Real one{FromInteger(Word{1}).value}; // 1.0 + Real bigStep{one}; + Real smallStep{one}; + ValueWithRealFlags scaled; + Rounding rounding{Rounding::TiesToEven}; // pmk? try ToZero + if (Compare(minFractionlessValue) == Relation::Less) { + // Scale smaller value up by a power of ten so that it loses no bit when + // converted to integer. + Real next{maxExactPowerOfTen}; + while (Multiply(next, rounding).value.Compare(maxExactSignedValue) == + Relation::Less) { + result.value.decimalExponent -= decimalDigits; + bigStep = next; + next = next.Multiply(maxExactPowerOfTen, rounding).value; + } + Real bigger{Multiply(bigStep, rounding).value}; + next = ten; + while (bigger.Multiply(next, rounding).value.Compare(maxExactSignedValue) == + Relation::Less) { + --result.value.decimalExponent; + smallStep = next; + next = next.Multiply(ten, rounding).value; + } + scaled = Multiply(smallStep, rounding); + scaled.value = + scaled.value.Multiply(bigStep, rounding).AccumulateFlags(scaled.flags); + } else { + // Scale larger value down by a power of ten so that it does not overflow + // when converted to integer. + Real last{maxExactSignedValue}; + Real next{last.Multiply(maxExactPowerOfTen, rounding).value}; + while (Compare(next) != Relation::Less) { + result.value.decimalExponent += decimalDigits; + bigStep = bigStep.Multiply(maxExactPowerOfTen, rounding).value; + last = next; + next = next.Multiply(maxExactPowerOfTen, rounding).value; + } + next = last; + while (Compare(next) == Relation::Greater) { + ++result.value.decimalExponent; + smallStep = smallStep.Multiply(ten, rounding).value; + next = last.Multiply(smallStep, rounding).value; + } + scaled = Divide(smallStep, rounding); + scaled.value = + scaled.value.Divide(bigStep, rounding).AccumulateFlags(scaled.flags); + } + + scaled.flags.reset(RealFlag::Inexact); + CHECK(scaled.flags.empty()); + auto asInteger{scaled.value.template ToInteger()}; + asInteger.flags.reset(RealFlag::Inexact); + CHECK(asInteger.flags.empty()); + result.value.integer = asInteger.value; + + // Canonicalize to the minimum integer value + if (!result.value.integer.IsZero()) { + while (true) { + auto qr{result.value.integer.DivideUnsigned(10)}; + if (!qr.remainder.IsZero()) { + break; + } + ++result.value.decimalExponent; + result.value.integer = qr.quotient; + } + } + + return result; +} + template class Real, 11>; template class Real, 24>; template class Real, 53>; template class Real, 64, false>; template class Real, 112>; -} +} \ No newline at end of file diff --git a/flang/lib/evaluate/real.h b/flang/lib/evaluate/real.h index 89fa4ef1b7f0..14ea52110551 100644 --- a/flang/lib/evaluate/real.h +++ b/flang/lib/evaluate/real.h @@ -20,6 +20,7 @@ #include "rounding-bits.h" #include #include +#include #include // Some environments, viz. clang on Darwin, allow the macro HUGE @@ -48,6 +49,17 @@ public: static constexpr std::uint64_t maxExponent{(1 << exponentBits) - 1}; static constexpr std::uint64_t exponentBias{maxExponent / 2}; + // Decimal precision + // log(2)/log(10) = 0.30102999... in any base; avoid floating constexpr + static constexpr int decimalDigits{(precision * 30103) / 100000}; + + // Associates a decimal exponent with an integral value for formmatting. + struct ScaledDecimal { + bool negative{false}; + Word integer; + int decimalExponent{0}; // Exxx + }; + template friend class Real; constexpr Real() {} // +0.0 @@ -56,7 +68,7 @@ public: constexpr Real &operator=(const Real &) = default; constexpr Real &operator=(Real &&) = default; - // TODO AINT/ANINT, CEILING, FLOOR, DIM, MAX, MIN, DPROD, FRACTION + // TODO ANINT, CEILING, FLOOR, DIM, MAX, MIN, DPROD, FRACTION // HUGE, INT/NINT, MAXEXPONENT, MINEXPONENT, NEAREST, OUT_OF_RANGE, // PRECISION, HUGE, TINY, RRSPACING/SPACING, SCALE, SET_EXPONENT, SIGN @@ -174,20 +186,44 @@ public: return result; } + // Truncation to integer in same real format. + constexpr ValueWithRealFlags AINT() const { + ValueWithRealFlags result{*this}; + if (IsNotANumber()) { + result.flags.set(RealFlag::InvalidArgument); + result.value = NotANumber(); + } else if (IsInfinite()) { + result.flags.set(RealFlag::Overflow); + } else { + std::uint64_t exponent{Exponent()}; + if (exponent < exponentBias) { // |x| < 1.0 + result.value.Normalize(IsNegative(), 0, Fraction{}); // +/-0.0 + } else { + constexpr std::uint64_t noClipExponent{exponentBias + precision - 1}; + if (int clip = noClipExponent - exponent; clip > 0) { + result.value.word_ = result.value.word_.IAND(Word::MASKR(clip).NOT()); + } + } + } + return result; + } + template constexpr ValueWithRealFlags ToInteger() const { + ValueWithRealFlags result; + if (IsNotANumber()) { + result.flags.set(RealFlag::InvalidArgument); + result.value = result.value.HUGE(); + return result; + } bool isNegative{IsNegative()}; std::uint64_t exponent{Exponent()}; Fraction fraction{GetFraction()}; - ValueWithRealFlags result; - if (exponent == maxExponent && !fraction.IsZero()) { // NaN - result.flags.set(RealFlag::InvalidArgument); - result.value = result.value.HUGE(); - } else if (exponent >= maxExponent || // +/-Inf + if (exponent >= maxExponent || // +/-Inf exponent >= exponentBias + result.value.bits) { // too big if (isNegative) { - result.value = result.value.MASKL(1); + result.value = result.value.MASKL(1); // most negative integer value } else { - result.value = result.value.HUGE(); + result.value = result.value.HUGE(); // most positive integer value } result.flags.set(RealFlag::Overflow); } else if (exponent < exponentBias) { // |x| < 1.0 -> 0 @@ -258,8 +294,12 @@ public: return result; } + // Represents the number as "J*(10**K)" where J and K are integers. + ValueWithRealFlags AsScaledDecimal() const; + constexpr Word RawBits() const { return word_; } + // Extracts "raw" biased exponent field. constexpr std::uint64_t Exponent() const { return word_.IBITS(significandBits, exponentBits).ToUInt64(); } @@ -268,6 +308,10 @@ public: const char *&, Rounding rounding = defaultRounding); std::string DumpHexadecimal() const; + // Emits a character representation for an equivalent Fortran constant + // or parenthesized constant expression that produces this value. + std::ostream &AsFortran(std::ostream &, int kind) const; + private: using Fraction = Integer; // all bits made explicit using Significand = Integer; // no implicit bit diff --git a/flang/test/evaluate/real.cc b/flang/test/evaluate/real.cc index e14c3c2fc35b..9064f793ae23 100644 --- a/flang/test/evaluate/real.cc +++ b/flang/test/evaluate/real.cc @@ -15,6 +15,7 @@ #include "fp-testing.h" #include "testing.h" #include "../../lib/evaluate/type.h" +#include #include #include @@ -158,6 +159,7 @@ template void basicTests(int rm, Rounding rounding) { TEST(ivf.flags.empty())(ldesc); MATCH(x, ivf.value.ToUInt64())(ldesc); } + TEST(vr.value.AINT().value.Compare(vr.value) == Relation::Equal)(ldesc); ix = ix.Negate().value; TEST(ix.IsNegative())(ldesc); x = -x; @@ -181,6 +183,7 @@ template void basicTests(int rm, Rounding rounding) { MATCH(x, ivf.value.ToUInt64())(ldesc); MATCH(nx, ivf.value.ToInt64())(ldesc); } + TEST(vr.value.AINT().value.Compare(vr.value) == Relation::Equal)(ldesc); } } @@ -207,40 +210,64 @@ std::uint64_t MakeReal(std::uint64_t n) { return x; } -std::uint32_t NormalizeNaN(std::uint32_t x) { - if ((x & 0x7f800000) == 0x7f800000 && (x & 0x007fffff) != 0) { +inline bool IsNaN(std::uint32_t x) { + return (x & 0x7f800000) == 0x7f800000 && (x & 0x007fffff) != 0; +} + +inline bool IsNaN(std::uint64_t x) { + return (x & 0x7ff0000000000000) == 0x7ff0000000000000 && + (x & 0x000fffffffffffff) != 0; +} + +inline bool IsInfinite(std::uint32_t x) { + return (x & 0x7fffffff) == 0x7f800000; +} + +inline bool IsInfinite(std::uint64_t x) { + return (x & 0x7fffffffffffffff) == 0x7ff0000000000000; +} + +inline std::uint32_t NormalizeNaN(std::uint32_t x) { + if (IsNaN(x)) { x = 0x7fe00000; } return x; } -std::uint64_t NormalizeNaN(std::uint64_t x) { - if ((x & 0x7ff0000000000000) == 0x7ff0000000000000 && - (x & 0x000fffffffffffff) != 0) { +inline std::uint64_t NormalizeNaN(std::uint64_t x) { + if (IsNaN(x)) { x = 0x7ffc000000000000; } return x; } +enum FlagBits { + Overflow = 1, + DivideByZero = 2, + InvalidArgument = 4, + Underflow = 8, + Inexact = 16, +}; + std::uint32_t FlagsToBits(const RealFlags &flags) { std::uint32_t bits{0}; #ifndef __clang__ // TODO: clang support for fenv.h is broken, so tests of flag settings // are disabled. if (flags.test(RealFlag::Overflow)) { - bits |= 1; + bits |= Overflow; } if (flags.test(RealFlag::DivideByZero)) { - bits |= 2; + bits |= DivideByZero; } if (flags.test(RealFlag::InvalidArgument)) { - bits |= 4; + bits |= InvalidArgument; } if (flags.test(RealFlag::Underflow)) { - bits |= 8; + bits |= Underflow; } if (flags.test(RealFlag::Inexact)) { - bits |= 0x10; + bits |= Inexact; } #endif // __clang__ return bits; @@ -288,10 +315,62 @@ void subsetTests(int pass, Rounding rounding, std::uint32_t opds) { ScopedHostFloatingPointEnvironment fpenv; for (UINT j{0}; j < opds; ++j) { + UINT rj{MakeReal(j)}; u.ui = rj; FLT fj{u.f}; REAL x{typename REAL::Word{std::uint64_t{rj}}}; + + // unary operations + { + ValueWithRealFlags aint{x.AINT()}; +#ifndef __clang__ // broken and also slow + fpenv.ClearFlags(); +#endif + FLT fcheck{std::trunc(fj)}; + auto actualFlags{FlagsToBits(fpenv.CurrentFlags())}; + actualFlags &= ~Inexact; // x86 std::trunc can set Inexact; AINT ain't + u.f = fcheck; + if (IsNaN(u.ui)) { + actualFlags |= InvalidArgument; // x86 std::trunc(NaN) WAR + } + UINT rcheck{NormalizeNaN(u.ui)}; + UINT check = aint.value.RawBits().ToUInt64(); + MATCH(rcheck, check) + ("%d AINT(0x%llx)", pass, static_cast(rj)); + MATCH(actualFlags, FlagsToBits(aint.flags)) + ("%d AINT(0x%llx)", pass, static_cast(rj)); + } + + { + MATCH(IsNaN(rj), x.IsNotANumber()) + ("%d IsNaN(0x%llx)", pass, static_cast(rj)); + MATCH(IsInfinite(rj), x.IsInfinite()) + ("%d IsInfinite(0x%llx)", pass, static_cast(rj)); + + auto scaled{x.AsScaledDecimal()}; + if (IsNaN(rj)) { + TEST(scaled.flags.test(RealFlag::InvalidArgument)) + ("%d invalid(0x%llx)", pass, static_cast(rj)); + } else if (IsInfinite(rj)) { + TEST(scaled.flags.test(RealFlag::Overflow)) + ("%d overflow(0x%llx)", pass, static_cast(rj)); + } else { + auto integer{scaled.value.integer.ToInt64()}; + MATCH(x.IsNegative(), scaled.value.negative) + ("%d IsNegative(0x%llx)", pass, static_cast(rj)); + char buffer[128], *dummy; + snprintf(buffer, sizeof buffer, "%c%lld.0E%d", + scaled.value.negative ? '-' : ' ', static_cast(integer), + scaled.value.decimalExponent); + u.f = std::strtold(buffer, &dummy); + MATCH(rj, u.ui) + ("%d scaled decimal 0x%llx %s %Lg", pass, static_cast(rj), + buffer, static_cast(u.f)); + } + } + + // dyadic operations for (UINT k{0}; k < opds; ++k) { UINT rk{MakeReal(k)}; u.ui = rk;