[flang] save work in progress
authorpeter klausler <pklausler@nvidia.com>
Thu, 8 Nov 2018 17:32:28 +0000 (09:32 -0800)
committerpeter klausler <pklausler@nvidia.com>
Thu, 8 Nov 2018 17:38:04 +0000 (09:38 -0800)
Original-commit: flang-compiler/f18@98bac3d297af460ae80467f257f724cd7b041dbd
Reviewed-on: https://github.com/flang-compiler/f18/pull/225
Tree-same-pre-rewrite: false

flang/lib/evaluate/integer.h
flang/lib/evaluate/real.cc
flang/lib/evaluate/real.h
flang/test/evaluate/real.cc

index 6f0dfb6..0c97be3 100644 (file)
@@ -87,6 +87,7 @@ private:
   static constexpr Part topPartMask{static_cast<Part>(~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>;
index 9eef531..3c5778e 100644 (file)
@@ -481,9 +481,160 @@ std::string Real<W, P, IM>::DumpHexadecimal() const {
   }
 }
 
+template<typename W, int P, bool IM>
+std::ostream &Real<W, P, IM>::AsFortran(std::ostream &o, int kind) const {
+  ValueWithRealFlags<ScaledDecimal> 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<typename W, int P, bool IM>
+auto Real<W, P, IM>::AsScaledDecimal() const
+    -> ValueWithRealFlags<ScaledDecimal> {
+
+  // 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<ScaledDecimal> 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<Real> 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<Word>()};
+  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<Integer<16>, 11>;
 template class Real<Integer<32>, 24>;
 template class Real<Integer<64>, 53>;
 template class Real<Integer<80>, 64, false>;
 template class Real<Integer<128>, 112>;
-}
+}
\ No newline at end of file
index 89fa4ef..14ea521 100644 (file)
@@ -20,6 +20,7 @@
 #include "rounding-bits.h"
 #include <cinttypes>
 #include <limits>
+#include <ostream>
 #include <string>
 
 // 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<typename W, int P, bool I> 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<Real> AINT() const {
+    ValueWithRealFlags<Real> 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<typename INT> constexpr ValueWithRealFlags<INT> ToInteger() const {
-    bool isNegative{IsNegative()};
-    std::uint64_t exponent{Exponent()};
-    Fraction fraction{GetFraction()};
     ValueWithRealFlags<INT> result;
-    if (exponent == maxExponent && !fraction.IsZero()) {  // NaN
+    if (IsNotANumber()) {
       result.flags.set(RealFlag::InvalidArgument);
       result.value = result.value.HUGE();
-    } else if (exponent >= maxExponent ||  // +/-Inf
+      return result;
+    }
+    bool isNegative{IsNegative()};
+    std::uint64_t exponent{Exponent()};
+    Fraction fraction{GetFraction()};
+    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<ScaledDecimal> 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<precision>;  // all bits made explicit
   using Significand = Integer<significandBits>;  // no implicit bit
index e14c3c2..9064f79 100644 (file)
@@ -15,6 +15,7 @@
 #include "fp-testing.h"
 #include "testing.h"
 #include "../../lib/evaluate/type.h"
+#include <cmath>
 #include <cstdio>
 #include <cstdlib>
 
@@ -158,6 +159,7 @@ template<typename R> 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<typename R> 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<REAL> 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<long long>(rj));
+      MATCH(actualFlags, FlagsToBits(aint.flags))
+      ("%d AINT(0x%llx)", pass, static_cast<long long>(rj));
+    }
+
+    {
+      MATCH(IsNaN(rj), x.IsNotANumber())
+      ("%d IsNaN(0x%llx)", pass, static_cast<long long>(rj));
+      MATCH(IsInfinite(rj), x.IsInfinite())
+      ("%d IsInfinite(0x%llx)", pass, static_cast<long long>(rj));
+
+      auto scaled{x.AsScaledDecimal()};
+      if (IsNaN(rj)) {
+        TEST(scaled.flags.test(RealFlag::InvalidArgument))
+        ("%d invalid(0x%llx)", pass, static_cast<long long>(rj));
+      } else if (IsInfinite(rj)) {
+        TEST(scaled.flags.test(RealFlag::Overflow))
+        ("%d overflow(0x%llx)", pass, static_cast<long long>(rj));
+      } else {
+        auto integer{scaled.value.integer.ToInt64()};
+        MATCH(x.IsNegative(), scaled.value.negative)
+        ("%d IsNegative(0x%llx)", pass, static_cast<long long>(rj));
+        char buffer[128], *dummy;
+        snprintf(buffer, sizeof buffer, "%c%lld.0E%d",
+            scaled.value.negative ? '-' : ' ', static_cast<long long>(integer),
+            scaled.value.decimalExponent);
+        u.f = std::strtold(buffer, &dummy);
+        MATCH(rj, u.ui)
+        ("%d scaled decimal 0x%llx %s %Lg", pass, static_cast<long long>(rj),
+            buffer, static_cast<long double>(u.f));
+      }
+    }
+
+    // dyadic operations
     for (UINT k{0}; k < opds; ++k) {
       UINT rk{MakeReal(k)};
       u.ui = rk;