static constexpr Part topPartMask{static_cast<Part>(~0) >> extraTopPartBits};
public:
+ // Some types used for member function results
struct ValueWithOverflow {
Integer value;
bool overflow;
}
}
- Part part_[parts];
+ Part part_[parts]{};
};
extern template class Integer<8>;
}
}
+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
#include "rounding-bits.h"
#include <cinttypes>
#include <limits>
+#include <ostream>
#include <string>
// Some environments, viz. clang on Darwin, allow the macro HUGE
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
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
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
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();
}
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
#include "fp-testing.h"
#include "testing.h"
#include "../../lib/evaluate/type.h"
+#include <cmath>
#include <cstdio>
#include <cstdlib>
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;
MATCH(x, ivf.value.ToUInt64())(ldesc);
MATCH(nx, ivf.value.ToInt64())(ldesc);
}
+ TEST(vr.value.AINT().value.Compare(vr.value) == Relation::Equal)(ldesc);
}
}
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;
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;