[flang] debugging
Original-commit: flang-compiler/f18@9f30eac130 Reviewed-on: https://github.com/flang-compiler/f18/pull/231 Tree-same-pre-rewrite: false
This commit is contained in:
parent
1df60f3ceb
commit
2fe5b128bd
|
@ -14,10 +14,14 @@
|
|||
|
||||
#include "decimal.h"
|
||||
#include "integer.h"
|
||||
#include "leading-zero-bit-count.h"
|
||||
#include "../common/idioms.h"
|
||||
#include <iostream> // TODO pmk rm
|
||||
bool pmk{false};
|
||||
|
||||
namespace Fortran::evaluate::value {
|
||||
|
||||
static std::ostream *debug{nullptr};
|
||||
static std::ostream *debug{nullptr}; // pmk constexpr
|
||||
|
||||
template<typename REAL>
|
||||
std::ostream &Decimal<REAL>::Dump(std::ostream &o) const {
|
||||
|
@ -30,14 +34,15 @@ std::ostream &Decimal<REAL>::Dump(std::ostream &o) const {
|
|||
return o << " e" << exponent_ << '\n';
|
||||
}
|
||||
|
||||
template<typename REAL> void Decimal<REAL>::FromReal(const Real &x) {
|
||||
template<typename REAL> void Decimal<REAL>::FromReal(const REAL &x) {
|
||||
if (x.IsNegative()) {
|
||||
FromReal(x.Negate());
|
||||
isNegative_ = true;
|
||||
return;
|
||||
}
|
||||
if (debug) {
|
||||
*debug << "FromReal(" << x.DumpHexadecimal() << ")\n";
|
||||
*debug << "FromReal(" << x.DumpHexadecimal() << ") bits " << Real::bits
|
||||
<< '\n';
|
||||
}
|
||||
if (x.IsZero()) {
|
||||
return;
|
||||
|
@ -58,7 +63,12 @@ template<typename REAL> void Decimal<REAL>::FromReal(const Real &x) {
|
|||
if (debug) {
|
||||
*debug << "second twoPow " << twoPow << ", lshift " << lshift << '\n';
|
||||
}
|
||||
SetTo(x.GetFraction().SHIFTL(lshift));
|
||||
using Word = typename Real::Word;
|
||||
Word word{Word::ConvertUnsigned(x.GetFraction()).value};
|
||||
SetTo(word.SHIFTL(lshift));
|
||||
if (debug) {
|
||||
Dump(*debug);
|
||||
}
|
||||
|
||||
for (; twoPow > 0 && IsDivisibleBy<5>(); --twoPow) {
|
||||
DivideBy<5>();
|
||||
|
@ -115,9 +125,23 @@ public:
|
|||
using Real = REAL;
|
||||
using Word = typename Real::Word;
|
||||
|
||||
void SetTo(Word n) {
|
||||
word_ = n;
|
||||
void SetTo(std::uint64_t n) {
|
||||
exponent_ = 0;
|
||||
if constexpr (Word::bits >= 8 * sizeof n) {
|
||||
word_ = n;
|
||||
} else {
|
||||
int shift{64 - LeadingZeroBitCount(n) - Word::bits};
|
||||
if (shift <= 0) {
|
||||
word_ = n;
|
||||
} else {
|
||||
word_ = n >> shift;
|
||||
exponent_ = shift;
|
||||
bool sticky{n << (64 - shift) != 0};
|
||||
if (sticky) {
|
||||
word_ = word_.IOR(Word{1});
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
void MultiplyAndAdd(std::uint32_t n, std::uint32_t plus = 0) {
|
||||
|
@ -134,10 +158,11 @@ public:
|
|||
sticky |= product.lower.BTEST(0);
|
||||
product.lower = product.lower.DSHIFTR(product.upper, 1);
|
||||
product.upper = product.upper.SHIFTR(1);
|
||||
++exponent_;
|
||||
}
|
||||
word_ = product.lower;
|
||||
if (sticky) {
|
||||
word_ = word_.IOR(word_.MASKR(1));
|
||||
word_ = word_.IOR(Word{1});
|
||||
}
|
||||
}
|
||||
|
||||
|
@ -149,7 +174,8 @@ public:
|
|||
|
||||
void AdjustExponent(int by = -1) { exponent_ += by; }
|
||||
|
||||
Real ToReal(bool isNegative = false) const;
|
||||
ValueWithRealFlags<Real> ToReal(
|
||||
bool isNegative = false, Rounding rounding = Rounding::TiesToEven) const;
|
||||
|
||||
private:
|
||||
Word word_{0};
|
||||
|
@ -162,34 +188,70 @@ std::ostream &IntermediateFloat<REAL>::Dump(std::ostream &o) const {
|
|||
}
|
||||
|
||||
template<typename REAL>
|
||||
REAL IntermediateFloat<REAL>::ToReal(bool isNegative) const {
|
||||
ValueWithRealFlags<REAL> IntermediateFloat<REAL>::ToReal(
|
||||
bool isNegative, Rounding rounding) const {
|
||||
if (word_.IsNegative()) {
|
||||
// word_ represents an unsigned quantity, so shift it down if the MSB is set
|
||||
IntermediateFloat shifted;
|
||||
shifted.word_ =
|
||||
word_.SHIFTR(1).IOR(word_.IAND(word_.MASKR(1))); // sticky bit
|
||||
Word sticky{word_.IAND(Word{1})};
|
||||
shifted.word_ = word_.SHIFTR(1).IOR(sticky);
|
||||
shifted.exponent_ = exponent_ + 1;
|
||||
return shifted.ToReal(false);
|
||||
if (debug) {
|
||||
shifted.Dump(*debug << "IntermediateFloat::ToReal: shifted: ") << '\n';
|
||||
}
|
||||
return shifted.ToReal(isNegative, rounding);
|
||||
}
|
||||
ValueWithRealFlags<Real> result;
|
||||
if (isNegative) {
|
||||
result = Real::FromInteger(word_.Negate().value, rounding);
|
||||
} else {
|
||||
result = Real::FromInteger(word_, rounding);
|
||||
}
|
||||
if (debug) {
|
||||
*debug << "IntermediateFloat::ToReal: after FromInteger: "
|
||||
<< result.value.DumpHexadecimal() << " * 2**" << exponent_ << '\n';
|
||||
}
|
||||
Real result = Real::FromInteger(word_).value;
|
||||
int expo{exponent_};
|
||||
while (expo + Real::exponentBias < 1) {
|
||||
Real twoPow{Word{1}.SHIFTL(Real::significandBits)};
|
||||
result = result.Multiply(twoPow).value;
|
||||
Real twoPow{Word{1}.SHIFTL(Real::significandBits)}; // min normal value
|
||||
result.value = result.value.Multiply(twoPow).AccumulateFlags(result.flags);
|
||||
expo += Real::exponentBias - 1;
|
||||
if (debug) {
|
||||
*debug << "IntermediateFloat::ToReal: reduced: "
|
||||
<< result.value.DumpHexadecimal() << " * 2**" << expo << '\n';
|
||||
}
|
||||
}
|
||||
while (expo + Real::exponentBias >= Real::maxExponent) {
|
||||
Real twoPow{Word{Real::maxExponent - 1}.SHIFTL(Real::significandBits)};
|
||||
result = result.Multiply(twoPow).value;
|
||||
result.value = result.value.Multiply(twoPow).AccumulateFlags(result.flags);
|
||||
expo += Real::maxExponent - 1 - Real::exponentBias;
|
||||
if (debug) {
|
||||
*debug << "IntermediateFloat::ToReal: magnified: "
|
||||
<< result.value.DumpHexadecimal() << " * 2**" << expo << '\n';
|
||||
}
|
||||
}
|
||||
Real twoPow{Word{expo + Real::exponentBias}.SHIFTL(Real::significandBits)};
|
||||
if (isNegative) {
|
||||
twoPow = twoPow.Negate();
|
||||
if (debug) {
|
||||
*debug << "IntermediateFloat::ToReal: twoPow: " << twoPow.DumpHexadecimal()
|
||||
<< '\n';
|
||||
}
|
||||
return result.Multiply(twoPow).value;
|
||||
result.value = result.value.Multiply(twoPow).AccumulateFlags(result.flags);
|
||||
return result;
|
||||
}
|
||||
|
||||
template<typename REAL> REAL Decimal<REAL>::ToReal(const char *&p) {
|
||||
template<typename REAL>
|
||||
ValueWithRealFlags<REAL> Decimal<REAL>::ToReal(
|
||||
const char *&p, Rounding rounding) {
|
||||
if (std::string{p} == "3.4028234663852885981170418348451692544e38_4") {
|
||||
debug = &std::cerr;
|
||||
pmk = true;
|
||||
} else {
|
||||
debug = nullptr;
|
||||
pmk = false;
|
||||
} // pmk
|
||||
if (debug) {
|
||||
*debug << "ToReal('" << p << "')\n";
|
||||
}
|
||||
while (*p == ' ') {
|
||||
++p;
|
||||
}
|
||||
|
@ -218,11 +280,13 @@ template<typename REAL> REAL Decimal<REAL>::ToReal(const char *&p) {
|
|||
++exponent_;
|
||||
}
|
||||
} else {
|
||||
MultiplyBy<10>(c - '0');
|
||||
int carry{MultiplyBy<10>(c - '0')};
|
||||
CHECK(carry == 0);
|
||||
if (decimalPoint) {
|
||||
--exponent_;
|
||||
}
|
||||
}
|
||||
if (debug) Dump(*debug << "ToReal in loop, p at '" << p << "'\n'");
|
||||
}
|
||||
|
||||
switch (*p) {
|
||||
|
@ -232,9 +296,8 @@ template<typename REAL> REAL Decimal<REAL>::ToReal(const char *&p) {
|
|||
case 'D':
|
||||
case 'q':
|
||||
case 'Q':
|
||||
++p;
|
||||
bool negExpo{*p == '-'};
|
||||
if (*p == '+' || *p == '-') {
|
||||
bool negExpo{*++p == '-'};
|
||||
if (negExpo || *p == '+') {
|
||||
++p;
|
||||
}
|
||||
char *q;
|
||||
|
@ -248,10 +311,14 @@ template<typename REAL> REAL Decimal<REAL>::ToReal(const char *&p) {
|
|||
}
|
||||
|
||||
if (debug) {
|
||||
Dump(*debug << "ToReal start: ");
|
||||
Dump(*debug << "ToReal start, p at '" << p << "'\n");
|
||||
}
|
||||
if (IsZero()) {
|
||||
return Real{}.Negate(); // -0.0
|
||||
ValueWithRealFlags<Real> result;
|
||||
if (isNegative_) {
|
||||
result.value = Real{}.Negate(); // -0.0
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
if (exponent_ < 0) {
|
||||
|
@ -276,12 +343,23 @@ template<typename REAL> REAL Decimal<REAL>::ToReal(const char *&p) {
|
|||
// result. The most significant digit can be moved directly;
|
||||
// lesser-order digits require transfer of carries.
|
||||
if (exponent_ >= -(digits_ - 1) * log10Quintillion) {
|
||||
if (debug) {
|
||||
Dump(*debug << "converting integer part:");
|
||||
}
|
||||
f.SetTo(digit_[--digits_]);
|
||||
if (debug) {
|
||||
f.Dump(*debug << "after converting top digit: ") << '\n';
|
||||
Dump(*debug);
|
||||
}
|
||||
while (exponent_ > -digits_ * log10Quintillion) {
|
||||
digitLimit_ = digits_;
|
||||
int carry{MultiplyBy<10>()};
|
||||
f.MultiplyAndAdd(10, carry);
|
||||
--exponent_;
|
||||
if (debug) {
|
||||
f.Dump(*debug << "foot of loop after carry " << carry << ": ") << '\n';
|
||||
Dump(*debug);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
@ -324,7 +402,8 @@ template<typename REAL> REAL Decimal<REAL>::ToReal(const char *&p) {
|
|||
Dump(*debug);
|
||||
}
|
||||
|
||||
return f.ToReal(isNegative_);
|
||||
auto result{f.ToReal(isNegative_, rounding)};
|
||||
return result;
|
||||
}
|
||||
|
||||
template<typename REAL>
|
||||
|
|
|
@ -15,6 +15,8 @@
|
|||
#ifndef FORTRAN_EVALUATE_DECIMAL_H_
|
||||
#define FORTRAN_EVALUATE_DECIMAL_H_
|
||||
|
||||
#include "common.h"
|
||||
#include "integer.h"
|
||||
#include "real.h"
|
||||
#include <cinttypes>
|
||||
#include <limits>
|
||||
|
@ -62,8 +64,15 @@ public:
|
|||
first_ = 0;
|
||||
exponent_ = 0;
|
||||
}
|
||||
|
||||
void FromReal(const Real &);
|
||||
Real ToReal(const char *&); // arg left pointing to first unparsed char
|
||||
|
||||
// Convert a character representation of a floating-point value to
|
||||
// the underlying Real type. The reference argument is a pointer that
|
||||
// is left pointing to the first character that wasn't included.
|
||||
ValueWithRealFlags<Real> ToReal(
|
||||
const char *&, Rounding rounding = Rounding::TiesToEven);
|
||||
|
||||
std::string ToString(int maxDigits = 1000000) const;
|
||||
|
||||
private:
|
||||
|
@ -92,15 +101,23 @@ private:
|
|||
++exponent_;
|
||||
n = qr.quotient;
|
||||
}
|
||||
while (!n.IsZero() && digits_ < digitLimit_) {
|
||||
auto qr{n.DivideUnsigned(quintillion)};
|
||||
digit_[digits_++] = qr.remainder.ToUInt64();
|
||||
if (digits_ == first_ + 1 && digit_[first_] == 0) {
|
||||
++first_;
|
||||
if constexpr (INT::bits < 60) {
|
||||
// n is necessarily less than a quintillion
|
||||
if (!n.IsZero()) {
|
||||
digit_[digits_++] = n.ToUInt64();
|
||||
}
|
||||
n = qr.quotient;
|
||||
return 0;
|
||||
} else {
|
||||
while (!n.IsZero() && digits_ < digitLimit_) {
|
||||
auto qr{n.DivideUnsigned(quintillion)};
|
||||
digit_[digits_++] = qr.remainder.ToUInt64();
|
||||
if (digits_ == first_ + 1 && digit_[first_] == 0) {
|
||||
++first_;
|
||||
}
|
||||
n = qr.quotient;
|
||||
}
|
||||
return n;
|
||||
}
|
||||
return n;
|
||||
}
|
||||
|
||||
int RemoveLeastOrderZeroDigits() {
|
||||
|
|
|
@ -13,6 +13,7 @@
|
|||
// limitations under the License.
|
||||
|
||||
#include "real.h"
|
||||
#include "decimal.h"
|
||||
#include "int-power.h"
|
||||
#include "../common/idioms.h"
|
||||
#include "../parser/characters.h"
|
||||
|
@ -88,8 +89,8 @@ ValueWithRealFlags<Real<W, P, IM>> Real<W, P, IM>::Add(
|
|||
result.value = y; // x + +/-Inf -> +/-Inf
|
||||
return result;
|
||||
}
|
||||
std::uint64_t exponent{Exponent()};
|
||||
std::uint64_t yExponent{y.Exponent()};
|
||||
int exponent{Exponent()};
|
||||
int yExponent{y.Exponent()};
|
||||
if (exponent < yExponent) {
|
||||
// y is larger in magnitude; simplify by reversing operands
|
||||
return y.Add(*this, rounding);
|
||||
|
@ -267,9 +268,9 @@ ValueWithRealFlags<Real<W, P, IM>> Real<W, P, IM>::Divide(
|
|||
}
|
||||
|
||||
template<typename W, int P, bool IM>
|
||||
RealFlags Real<W, P, IM>::Normalize(bool negative, std::uint64_t exponent,
|
||||
RealFlags Real<W, P, IM>::Normalize(bool negative, int exponent,
|
||||
const Fraction &fraction, Rounding rounding, RoundingBits *roundingBits) {
|
||||
std::uint64_t lshift = fraction.LEADZ();
|
||||
int lshift{fraction.LEADZ()};
|
||||
if (lshift == fraction.bits /* fraction is zero */ &&
|
||||
(roundingBits == nullptr || roundingBits->empty())) {
|
||||
// No fraction, no rounding bits -> +/-0.0
|
||||
|
@ -329,7 +330,7 @@ RealFlags Real<W, P, IM>::Normalize(bool negative, std::uint64_t exponent,
|
|||
template<typename W, int P, bool IM>
|
||||
RealFlags Real<W, P, IM>::Round(
|
||||
Rounding rounding, const RoundingBits &bits, bool multiply) {
|
||||
std::uint64_t origExponent{Exponent()};
|
||||
int origExponent{Exponent()};
|
||||
RealFlags flags;
|
||||
bool inexact{!bits.empty()};
|
||||
if (inexact) {
|
||||
|
@ -339,7 +340,7 @@ RealFlags Real<W, P, IM>::Round(
|
|||
bits.MustRound(rounding, IsNegative(), word_.BTEST(0) /* is odd */)) {
|
||||
typename Fraction::ValueWithCarry sum{
|
||||
GetFraction().AddUnsigned(Fraction{}, true)};
|
||||
std::uint64_t newExponent{origExponent};
|
||||
int newExponent{origExponent};
|
||||
if (sum.carry) {
|
||||
// The fraction was all ones before rounding; sum.value is now zero
|
||||
sum.value = sum.value.IBSET(precision - 1);
|
||||
|
@ -372,8 +373,8 @@ RealFlags Real<W, P, IM>::Round(
|
|||
|
||||
template<typename W, int P, bool IM>
|
||||
void Real<W, P, IM>::NormalizeAndRound(ValueWithRealFlags<Real> &result,
|
||||
bool isNegative, std::uint64_t exponent, const Fraction &fraction,
|
||||
Rounding rounding, RoundingBits roundingBits, bool multiply) {
|
||||
bool isNegative, int exponent, const Fraction &fraction, Rounding rounding,
|
||||
RoundingBits roundingBits, bool multiply) {
|
||||
result.flags |= result.value.Normalize(
|
||||
isNegative, exponent, fraction, rounding, &roundingBits);
|
||||
result.flags |= result.value.Round(rounding, roundingBits, multiply);
|
||||
|
@ -382,126 +383,8 @@ void Real<W, P, IM>::NormalizeAndRound(ValueWithRealFlags<Real> &result,
|
|||
template<typename W, int P, bool IM>
|
||||
ValueWithRealFlags<Real<W, P, IM>> Real<W, P, IM>::Read(
|
||||
const char *&p, Rounding rounding) {
|
||||
ValueWithRealFlags<Real> result;
|
||||
while (*p == ' ') {
|
||||
++p;
|
||||
}
|
||||
bool isNegative{*p == '-'};
|
||||
if (*p == '-' || *p == '+') {
|
||||
++p;
|
||||
}
|
||||
Word integer{0};
|
||||
int decimalExponent{0};
|
||||
bool full{false};
|
||||
bool inFraction{false};
|
||||
for (;; ++p) {
|
||||
if (*p == '.') {
|
||||
if (inFraction) {
|
||||
break;
|
||||
}
|
||||
inFraction = true;
|
||||
} else if (!parser::IsDecimalDigit(*p)) {
|
||||
break;
|
||||
} else if (full) {
|
||||
if (!inFraction) {
|
||||
++decimalExponent;
|
||||
}
|
||||
} else {
|
||||
auto times10{integer.MultiplyUnsigned(Word{10})};
|
||||
if (!times10.upper.IsZero()) {
|
||||
full = true;
|
||||
if (!inFraction) {
|
||||
++decimalExponent;
|
||||
}
|
||||
} else {
|
||||
auto augmented{times10.lower.AddUnsigned(Word{*p - '0'})};
|
||||
if (augmented.carry) {
|
||||
full = true;
|
||||
if (!inFraction) {
|
||||
++decimalExponent;
|
||||
}
|
||||
} else {
|
||||
integer = augmented.value;
|
||||
if (inFraction) {
|
||||
--decimalExponent;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (parser::IsLetter(*p)) {
|
||||
++p;
|
||||
bool negExpo{*p == '-'};
|
||||
if (*p == '+' || *p == '-') {
|
||||
++p;
|
||||
}
|
||||
auto expo{Integer<32>::ReadUnsigned(p)};
|
||||
std::int64_t expoVal{expo.value.ToInt64()};
|
||||
if (expo.overflow) {
|
||||
expoVal = std::numeric_limits<std::int32_t>::max();
|
||||
} else if (negExpo) {
|
||||
expoVal *= -1;
|
||||
}
|
||||
decimalExponent += expoVal;
|
||||
}
|
||||
|
||||
int binaryExponent{exponentBias + bits - 1};
|
||||
if (integer.IsZero()) {
|
||||
decimalExponent = 0;
|
||||
} else {
|
||||
int leadz{integer.LEADZ()};
|
||||
binaryExponent -= leadz;
|
||||
integer = integer.SHIFTL(leadz);
|
||||
}
|
||||
for (; decimalExponent > 0; --decimalExponent) {
|
||||
auto times5{integer.MultiplyUnsigned(Word{5})};
|
||||
++binaryExponent;
|
||||
integer = times5.lower;
|
||||
for (; !times5.upper.IsZero(); times5.upper = times5.upper.SHIFTR(1)) {
|
||||
++binaryExponent;
|
||||
integer = integer.SHIFTR(1);
|
||||
if (times5.upper.BTEST(0)) {
|
||||
integer = integer.IBSET(bits - 1);
|
||||
}
|
||||
}
|
||||
}
|
||||
for (; decimalExponent < 0; ++decimalExponent) {
|
||||
auto div5{integer.DivideUnsigned(Word{5})};
|
||||
--binaryExponent;
|
||||
integer = div5.quotient;
|
||||
std::uint8_t lost = div5.remainder.ToUInt64() * 0x33;
|
||||
while (!integer.BTEST(bits - 1)) {
|
||||
integer = integer.SHIFTL(1);
|
||||
if (lost & 0x80) {
|
||||
integer = integer.IBSET(0);
|
||||
}
|
||||
lost <<= 1;
|
||||
--binaryExponent;
|
||||
}
|
||||
}
|
||||
|
||||
RoundingBits roundingBits;
|
||||
for (int j{0}; bits - j > precision; ++j) {
|
||||
roundingBits.ShiftRight(integer.BTEST(0));
|
||||
integer = integer.SHIFTR(1);
|
||||
}
|
||||
|
||||
Fraction fraction{Fraction::ConvertUnsigned(integer).value};
|
||||
while (binaryExponent < 1) {
|
||||
if (fraction.IsZero()) {
|
||||
binaryExponent = 0;
|
||||
break;
|
||||
} else {
|
||||
++binaryExponent;
|
||||
roundingBits.ShiftRight(fraction.BTEST(0));
|
||||
fraction = fraction.SHIFTR(1);
|
||||
}
|
||||
}
|
||||
|
||||
NormalizeAndRound(
|
||||
result, isNegative, binaryExponent, fraction, rounding, roundingBits);
|
||||
return result;
|
||||
Decimal<Real> decimal;
|
||||
return decimal.ToReal(p, rounding);
|
||||
}
|
||||
|
||||
template<typename W, int P, bool IM>
|
||||
|
@ -553,151 +436,30 @@ 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)) {
|
||||
if (IsNotANumber()) {
|
||||
o << "(0._" << kind << "/0.)";
|
||||
} else if (scaled.flags.test(RealFlag::Overflow)) {
|
||||
} else if (IsInfinite()) {
|
||||
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;
|
||||
bool parenthesize{IsNegative()};
|
||||
if (parenthesize) {
|
||||
o << "(";
|
||||
}
|
||||
Decimal<Real> decimal;
|
||||
decimal.FromReal(*this);
|
||||
o << decimal.ToString();
|
||||
o << '_' << kind;
|
||||
if (scaled.value.negative) {
|
||||
if (parenthesize) {
|
||||
o << ')';
|
||||
}
|
||||
}
|
||||
return o;
|
||||
}
|
||||
|
||||
template<typename W, int P, bool IM>
|
||||
auto Real<W, P, IM>::AsScaledDecimal(Rounding rounding) const
|
||||
-> ValueWithRealFlags<ScaledDecimal> {
|
||||
|
||||
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()) {
|
||||
if (rounding == Rounding::Up) {
|
||||
rounding = Rounding::Down;
|
||||
} else if (rounding == Rounding::Down) {
|
||||
rounding = Rounding::Up;
|
||||
}
|
||||
result = Negate().AsScaledDecimal(rounding);
|
||||
result.value.negative = true;
|
||||
return result;
|
||||
}
|
||||
if (IsZero()) {
|
||||
return result;
|
||||
}
|
||||
|
||||
Word fraction{Word::ConvertUnsigned(GetFraction()).value};
|
||||
Word asInt{fraction.SHIFTL(bits - precision)};
|
||||
int twoPower{UnbiasedExponent() - bits + 1};
|
||||
|
||||
// The original Real, a finite positive number, is now represented as a
|
||||
// left-justified integer and a binary exponent. In other words, asInt
|
||||
// has its sign bit set and the value of "asInt * 2.0**twoPower" equals
|
||||
// the original Real exactly; both asInt and twoPower are integers.
|
||||
//
|
||||
// To generate the scaled decimal result, we multiply or divide asInt by
|
||||
// 2**twoPower iteratively. Whenever the result of a multiplication or
|
||||
// division would overflow (or lose precision, respectively), we scale
|
||||
// asInt by dividing it by ten (or multiplying, respectively). Since
|
||||
// 10==2*5, the scaling is actually by a factor of five since the factor
|
||||
// of two can be accommodated by shifting. These scalings are recorded
|
||||
// in a decimal exponent. Once all of the "2.0**twoPower" scaling is
|
||||
// done, we have a value that is represented as
|
||||
// "asInt * 10.0**decimalExponent", and that is the penultimate result.
|
||||
|
||||
static constexpr Word five{5};
|
||||
if (twoPower > 0) {
|
||||
// Multiply asInt by 2**twoPower.
|
||||
static constexpr Word two{2};
|
||||
for (; twoPower > 0; --twoPower) {
|
||||
if (asInt.BTEST(bits - 1)) {
|
||||
// asInt * 2 would overflow: divide by five and increment the
|
||||
// decimal exponent (effectively dividing by ten and then multiplying
|
||||
// by two).
|
||||
auto qr{asInt.DivideUnsigned(five)};
|
||||
asInt = qr.quotient;
|
||||
++result.value.decimalExponent;
|
||||
if (twoPower > 1 &&
|
||||
qr.remainder.CompareUnsigned(two) == Ordering::Greater) {
|
||||
--twoPower;
|
||||
asInt = asInt.SHIFTL(1).IBSET(0);
|
||||
}
|
||||
} else {
|
||||
asInt = asInt.SHIFTL(1);
|
||||
}
|
||||
}
|
||||
} else {
|
||||
// Divide asInt by 2**(-twoPower).
|
||||
std::uint32_t lower{0};
|
||||
for (; twoPower < 0; ++twoPower) {
|
||||
auto times5{asInt.MultiplyUnsigned(five)};
|
||||
if (!times5.upper.IsZero()) {
|
||||
// asInt is too big to need scaling, just shift it down.
|
||||
lower >>= 1;
|
||||
if (asInt.BTEST(0)) {
|
||||
lower |= 1 << 31;
|
||||
}
|
||||
asInt = asInt.SHIFTR(1);
|
||||
} else {
|
||||
std::uint64_t lowerTimes5{lower * static_cast<std::uint64_t>(5)};
|
||||
std::uint32_t round = lowerTimes5 >> 32;
|
||||
auto rounded{times5.lower.AddUnsigned(Word{round})};
|
||||
if (rounded.carry) {
|
||||
// asInt is still too big to need scaling (rounding would overflow)
|
||||
lower >>= 1;
|
||||
if (asInt.BTEST(0)) {
|
||||
lower |= 1 << 31;
|
||||
}
|
||||
asInt = asInt.SHIFTR(1);
|
||||
} else {
|
||||
// asInt is small enough to be scaled; do so.
|
||||
--result.value.decimalExponent;
|
||||
lower = lowerTimes5;
|
||||
asInt = rounded.value;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Canonicalize to the minimum integer value: the final result will not
|
||||
// be a multiple of ten.
|
||||
result.value.integer = asInt;
|
||||
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>;
|
||||
|
|
|
@ -23,6 +23,9 @@
|
|||
#include <ostream>
|
||||
#include <string>
|
||||
|
||||
#include <iostream> // TODO pmk rm
|
||||
extern bool pmk;
|
||||
|
||||
// Some environments, viz. clang on Darwin, allow the macro HUGE
|
||||
// to leak out of <math.h> even when it is never directly included.
|
||||
#undef HUGE
|
||||
|
@ -39,17 +42,17 @@ namespace Fortran::evaluate::value {
|
|||
template<typename WORD, int PRECISION, bool IMPLICIT_MSB = true> class Real {
|
||||
public:
|
||||
using Word = WORD;
|
||||
static constexpr int bits{Word::bits};
|
||||
static constexpr int precision{PRECISION};
|
||||
using Fraction = Integer<precision>; // all bits made explicit
|
||||
|
||||
static constexpr int bits{Word::bits};
|
||||
static constexpr bool implicitMSB{IMPLICIT_MSB};
|
||||
static constexpr int significandBits{precision - implicitMSB};
|
||||
static constexpr int exponentBits{bits - significandBits - 1 /*sign*/};
|
||||
static_assert(precision > 0);
|
||||
static_assert(exponentBits > 1);
|
||||
static constexpr std::uint64_t maxExponent{(1 << exponentBits) - 1};
|
||||
static constexpr std::uint64_t exponentBias{maxExponent / 2};
|
||||
static_assert(exponentBits <= 16);
|
||||
static constexpr int maxExponent{(1 << exponentBits) - 1};
|
||||
static constexpr int exponentBias{maxExponent / 2};
|
||||
|
||||
// Decimal precision of a binary floating-point representation is
|
||||
// actually the same as the base-5 precision, as factors of two
|
||||
|
@ -58,14 +61,6 @@ public:
|
|||
// Calculate "precision*0.43" with integer arithmetic so as to be constexpr.
|
||||
static constexpr int decimalDigits{(precision * 43) / 100};
|
||||
|
||||
// Associates a decimal exponent with an integral value for formatting.
|
||||
// TODO pmk remove
|
||||
struct ScaledDecimal {
|
||||
bool negative{false};
|
||||
Word integer; // unsigned
|
||||
int decimalExponent{0}; // Exxx
|
||||
};
|
||||
|
||||
template<typename W, int P, bool I> friend class Real;
|
||||
|
||||
constexpr Real() {} // +0.0
|
||||
|
@ -125,11 +120,11 @@ public:
|
|||
const Real &, Rounding rounding = defaultRounding) const;
|
||||
|
||||
template<typename INT> constexpr INT EXPONENT() const {
|
||||
std::uint64_t exponent{Exponent()};
|
||||
int exponent{Exponent()};
|
||||
if (exponent == maxExponent) {
|
||||
return INT::HUGE();
|
||||
} else {
|
||||
int result = exponent - exponentBias;
|
||||
int result{exponent - exponentBias};
|
||||
if (IsDenormal()) {
|
||||
++result;
|
||||
}
|
||||
|
@ -180,9 +175,13 @@ public:
|
|||
return {}; // all bits zero -> +0.0
|
||||
}
|
||||
ValueWithRealFlags<Real> result;
|
||||
std::uint64_t exponent{exponentBias + absN.bits - leadz - 1};
|
||||
int exponent{exponentBias + absN.bits - leadz - 1};
|
||||
int bitsNeeded{absN.bits - (leadz + implicitMSB)};
|
||||
int bitsLost{bitsNeeded - significandBits};
|
||||
if (pmk)
|
||||
std::cerr << "pmk real.h exponent " << exponent << " bitsLost "
|
||||
<< bitsLost << " rounding " << static_cast<int>(rounding)
|
||||
<< '\n';
|
||||
if (bitsLost <= 0) {
|
||||
Fraction fraction{Fraction::ConvertUnsigned(absN).value};
|
||||
result.flags |= result.value.Normalize(
|
||||
|
@ -190,8 +189,13 @@ public:
|
|||
} else {
|
||||
Fraction fraction{Fraction::ConvertUnsigned(absN.SHIFTR(bitsLost)).value};
|
||||
result.flags |= result.value.Normalize(isNegative, exponent, fraction);
|
||||
if (pmk)
|
||||
std::cerr << "pmk Normalized " << result.value.DumpHexadecimal()
|
||||
<< '\n';
|
||||
RoundingBits roundingBits{absN, bitsLost};
|
||||
result.flags |= result.value.Round(rounding, roundingBits);
|
||||
if (pmk)
|
||||
std::cerr << "pmk Rounded " << result.value.DumpHexadecimal() << '\n';
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
@ -205,11 +209,11 @@ public:
|
|||
} else if (IsInfinite()) {
|
||||
result.flags.set(RealFlag::Overflow);
|
||||
} else {
|
||||
std::uint64_t exponent{Exponent()};
|
||||
int 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};
|
||||
constexpr int noClipExponent{exponentBias + precision - 1};
|
||||
if (int clip = noClipExponent - exponent; clip > 0) {
|
||||
result.value.word_ = result.value.word_.IAND(Word::MASKR(clip).NOT());
|
||||
}
|
||||
|
@ -226,7 +230,7 @@ public:
|
|||
return result;
|
||||
}
|
||||
bool isNegative{IsNegative()};
|
||||
std::uint64_t exponent{Exponent()};
|
||||
int exponent{Exponent()};
|
||||
Fraction fraction{GetFraction()};
|
||||
if (exponent >= maxExponent || // +/-Inf
|
||||
exponent >= exponentBias + result.value.bits) { // too big
|
||||
|
@ -243,7 +247,7 @@ public:
|
|||
}
|
||||
} else {
|
||||
// finite number |x| >= 1.0
|
||||
constexpr std::uint64_t noShiftExponent{exponentBias + precision - 1};
|
||||
constexpr int noShiftExponent{exponentBias + precision - 1};
|
||||
if (exponent < noShiftExponent) {
|
||||
int rshift = noShiftExponent - exponent;
|
||||
if (!fraction.IBITS(0, rshift).IsZero()) {
|
||||
|
@ -288,7 +292,7 @@ public:
|
|||
absX = x.Negate();
|
||||
}
|
||||
ValueWithRealFlags<Real> result;
|
||||
std::uint64_t exponent{exponentBias + x.Exponent() - A::exponentBias};
|
||||
int exponent{exponentBias + x.Exponent() - A::exponentBias};
|
||||
int bitsLost{A::precision - precision};
|
||||
typename A::Fraction xFraction{x.GetFraction()};
|
||||
if (bitsLost <= 0) {
|
||||
|
@ -304,17 +308,28 @@ public:
|
|||
return result;
|
||||
}
|
||||
|
||||
// Represents the number as "J*(10**K)" where J and K are integers.
|
||||
ValueWithRealFlags<ScaledDecimal> AsScaledDecimal(
|
||||
Rounding rounding = defaultRounding) const;
|
||||
|
||||
constexpr Word RawBits() const { return word_; }
|
||||
|
||||
// Extracts "raw" biased exponent field.
|
||||
constexpr std::uint64_t Exponent() const {
|
||||
constexpr int Exponent() const {
|
||||
return word_.IBITS(significandBits, exponentBits).ToUInt64();
|
||||
}
|
||||
|
||||
// Extracts the fraction; any implied bit is made explicit.
|
||||
constexpr Fraction GetFraction() const {
|
||||
Fraction result{Fraction::ConvertUnsigned(word_).value};
|
||||
if constexpr (!implicitMSB) {
|
||||
return result;
|
||||
} else {
|
||||
int exponent{Exponent()};
|
||||
if (exponent > 0 && exponent < maxExponent) {
|
||||
return result.IBSET(significandBits);
|
||||
} else {
|
||||
return result.IBCLR(significandBits);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// Extracts unbiased exponent value.
|
||||
constexpr int UnbiasedExponent() const {
|
||||
int exponent =
|
||||
|
@ -325,21 +340,6 @@ public:
|
|||
return exponent;
|
||||
}
|
||||
|
||||
// Extracts the fraction; any implied bit is made explicit.
|
||||
constexpr Fraction GetFraction() const {
|
||||
Fraction result{Fraction::ConvertUnsigned(word_).value};
|
||||
if constexpr (!implicitMSB) {
|
||||
return result;
|
||||
} else {
|
||||
std::uint64_t exponent{Exponent()};
|
||||
if (exponent > 0 && exponent < maxExponent) {
|
||||
return result.IBSET(significandBits);
|
||||
} else {
|
||||
return result.IBCLR(significandBits);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
static ValueWithRealFlags<Real> Read(
|
||||
const char *&, Rounding rounding = defaultRounding);
|
||||
std::string DumpHexadecimal() const;
|
||||
|
@ -355,8 +355,8 @@ private:
|
|||
return Significand::ConvertUnsigned(word_).value;
|
||||
}
|
||||
|
||||
constexpr std::int64_t CombineExponents(const Real &y, bool forDivide) const {
|
||||
std::int64_t exponent = Exponent(), yExponent = y.Exponent();
|
||||
constexpr int CombineExponents(const Real &y, bool forDivide) const {
|
||||
int exponent = Exponent(), yExponent = y.Exponent();
|
||||
// A zero exponent field value has the same weight as 1.
|
||||
exponent += !exponent;
|
||||
yExponent += !yExponent;
|
||||
|
@ -384,16 +384,16 @@ private:
|
|||
// The value is a number, and a zero fraction means a zero value (i.e.,
|
||||
// a maximal exponent and zero fraction doesn't signify infinity, although
|
||||
// this member function will detect overflow and encode infinities).
|
||||
RealFlags Normalize(bool negative, std::uint64_t exponent,
|
||||
const Fraction &fraction, Rounding rounding = defaultRounding,
|
||||
RealFlags Normalize(bool negative, int exponent, const Fraction &fraction,
|
||||
Rounding rounding = defaultRounding,
|
||||
RoundingBits *roundingBits = nullptr);
|
||||
|
||||
// Rounds a result, if necessary, in place.
|
||||
RealFlags Round(Rounding, const RoundingBits &, bool multiply = false);
|
||||
|
||||
static void NormalizeAndRound(ValueWithRealFlags<Real> &result,
|
||||
bool isNegative, std::uint64_t exponent, const Fraction &, Rounding,
|
||||
RoundingBits, bool multiply = false);
|
||||
bool isNegative, int exponent, const Fraction &, Rounding, RoundingBits,
|
||||
bool multiply = false);
|
||||
|
||||
Word word_{}; // an Integer<>
|
||||
};
|
||||
|
|
|
@ -18,6 +18,7 @@
|
|||
#include <cmath>
|
||||
#include <cstdio>
|
||||
#include <cstdlib>
|
||||
#include <sstream>
|
||||
#include <type_traits>
|
||||
|
||||
using namespace Fortran::evaluate;
|
||||
|
@ -228,6 +229,12 @@ inline bool IsInfinite(std::uint64_t x) {
|
|||
return (x & 0x7fffffffffffffff) == 0x7ff0000000000000;
|
||||
}
|
||||
|
||||
inline bool IsNegative(std::uint32_t x) { return (x & 0x80000000) != 0; }
|
||||
|
||||
inline bool IsNegative(std::uint64_t x) {
|
||||
return (x & 0x8000000000000000) != 0;
|
||||
}
|
||||
|
||||
inline std::uint32_t NormalizeNaN(std::uint32_t x) {
|
||||
if (IsNaN(x)) {
|
||||
x = 0x7fe00000;
|
||||
|
@ -386,28 +393,37 @@ void subsetTests(int pass, Rounding rounding, std::uint32_t opds) {
|
|||
MATCH(IsInfinite(rj), x.IsInfinite())
|
||||
("%d IsInfinite(0x%llx)", pass, static_cast<long long>(rj));
|
||||
|
||||
if (rounding == Rounding::TiesToEven) {
|
||||
auto scaled{x.AsScaledDecimal()};
|
||||
if (rounding == Rounding::TiesToEven ||
|
||||
rounding == Rounding::ToZero) { // pmk
|
||||
int kind{REAL::bits / 8};
|
||||
std::stringstream ss, css;
|
||||
x.AsFortran(ss, kind);
|
||||
std::string s{ss.str()};
|
||||
if (IsNaN(rj)) {
|
||||
TEST(scaled.flags.test(RealFlag::InvalidArgument))
|
||||
css << "(0._" << kind << "/0.)";
|
||||
MATCH(css.str(), s)
|
||||
("%d invalid(0x%llx)", pass, static_cast<long long>(rj));
|
||||
} else if (IsInfinite(rj)) {
|
||||
TEST(scaled.flags.test(RealFlag::Overflow))
|
||||
css << '(';
|
||||
if (IsNegative(rj)) {
|
||||
css << '-';
|
||||
}
|
||||
css << "1._" << kind << "/0.)";
|
||||
MATCH(css.str(), s)
|
||||
("%d overflow(0x%llx)", pass, static_cast<long long>(rj));
|
||||
} else {
|
||||
auto integer{scaled.value.integer.ToUInt64()};
|
||||
MATCH(x.IsNegative(), scaled.value.negative)
|
||||
("%d IsNegative(0x%llx)", pass, static_cast<long long>(rj));
|
||||
char buffer[128];
|
||||
const char *p = buffer;
|
||||
snprintf(buffer, sizeof buffer, "%c%llu.0E%d",
|
||||
"+-"[scaled.value.negative],
|
||||
static_cast<unsigned long long>(integer),
|
||||
scaled.value.decimalExponent);
|
||||
const char *p = s.data();
|
||||
if (*p == '(') {
|
||||
++p;
|
||||
}
|
||||
auto readBack{REAL::Read(p, rounding)};
|
||||
MATCH(rj, readBack.value.RawBits().ToUInt64())
|
||||
("%d scaled decimal 0x%llx %s", pass, static_cast<long long>(rj),
|
||||
buffer);
|
||||
("%d Read(AsFortran()) 0x%llx %s %g", pass,
|
||||
static_cast<long long>(rj), s.data(), static_cast<double>(fj));
|
||||
MATCH('_', *p)
|
||||
("%d Read(AsFortran()) 0x%llx %s %d", pass,
|
||||
static_cast<long long>(rj), s.data(),
|
||||
static_cast<int>(p - s.data()));
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
|
@ -74,6 +74,11 @@ FailureDetailPrinter Match(const char *file, int line, const char *want,
|
|||
}
|
||||
}
|
||||
|
||||
FailureDetailPrinter Match(const char *file, int line, const std::string &want,
|
||||
const char *gots, const std::string &got) {
|
||||
return Match(file, line, want.data(), gots, got);
|
||||
}
|
||||
|
||||
FailureDetailPrinter Compare(const char *file, int line, const char *xs,
|
||||
const char *rel, const char *ys, unsigned long long x,
|
||||
unsigned long long y) {
|
||||
|
|
|
@ -41,6 +41,8 @@ FailureDetailPrinter Match(const char *file, int line, unsigned long long want,
|
|||
const char *gots, unsigned long long got);
|
||||
FailureDetailPrinter Match(const char *file, int line, const char *want,
|
||||
const char *gots, const std::string &got);
|
||||
FailureDetailPrinter Match(const char *file, int line, const std::string &want,
|
||||
const char *gots, const std::string &got);
|
||||
FailureDetailPrinter Compare(const char *file, int line, const char *xs,
|
||||
const char *rel, const char *ys, unsigned long long x,
|
||||
unsigned long long y);
|
||||
|
|
Loading…
Reference in New Issue