blob: 80e677780db40e63631c8daccf4254f404a8e22a [file] [log] [blame]
//===-- include/flang/Evaluate/real.h ---------------------------*- C++ -*-===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
#ifndef FORTRAN_EVALUATE_REAL_H_
#define FORTRAN_EVALUATE_REAL_H_
#include "formatting.h"
#include "integer.h"
#include "rounding-bits.h"
#include "flang/Common/real.h"
#include "flang/Evaluate/common.h"
#include <cinttypes>
#include <limits>
#include <string>
// 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
namespace llvm {
class raw_ostream;
}
namespace Fortran::evaluate::value {
// LOG10(2.)*1E12
static constexpr std::int64_t ScaledLogBaseTenOfTwo{301029995664};
// Models IEEE binary floating-point numbers (IEEE 754-2008,
// ISO/IEC/IEEE 60559.2011). The first argument to this
// class template must be (or look like) an instance of Integer<>;
// the second specifies the number of effective bits (binary precision)
// in the fraction.
template <typename WORD, int PREC>
class Real : public common::RealDetails<PREC> {
public:
using Word = WORD;
static constexpr int binaryPrecision{PREC};
using Details = common::RealDetails<PREC>;
using Details::exponentBias;
using Details::exponentBits;
using Details::isImplicitMSB;
using Details::maxExponent;
using Details::significandBits;
static constexpr int bits{Word::bits};
static_assert(bits >= Details::bits);
using Fraction = Integer<binaryPrecision>; // all bits made explicit
template <typename W, int P> friend class Real;
constexpr Real() {} // +0.0
constexpr Real(const Real &) = default;
constexpr Real(Real &&) = default;
constexpr Real(const Word &bits) : word_{bits} {}
constexpr Real &operator=(const Real &) = default;
constexpr Real &operator=(Real &&) = default;
constexpr bool operator==(const Real &that) const {
return word_ == that.word_;
}
constexpr bool IsSignBitSet() const { return word_.BTEST(bits - 1); }
constexpr bool IsNegative() const {
return !IsNotANumber() && IsSignBitSet();
}
constexpr bool IsNotANumber() const {
return Exponent() == maxExponent && !GetSignificand().IsZero();
}
constexpr bool IsQuietNaN() const {
return Exponent() == maxExponent &&
GetSignificand().BTEST(significandBits - 1);
}
constexpr bool IsSignalingNaN() const {
return IsNotANumber() && !GetSignificand().BTEST(significandBits - 1);
}
constexpr bool IsInfinite() const {
return Exponent() == maxExponent && GetSignificand().IsZero();
}
constexpr bool IsFinite() const { return Exponent() != maxExponent; }
constexpr bool IsZero() const {
return Exponent() == 0 && GetSignificand().IsZero();
}
constexpr bool IsSubnormal() const {
return Exponent() == 0 && !GetSignificand().IsZero();
}
constexpr Real ABS() const { // non-arithmetic, no flags returned
return {word_.IBCLR(bits - 1)};
}
constexpr Real SetSign(bool toNegative) const { // non-arithmetic
if (toNegative) {
return {word_.IBSET(bits - 1)};
} else {
return ABS();
}
}
constexpr Real SIGN(const Real &x) const { return SetSign(x.IsSignBitSet()); }
constexpr Real Negate() const { return {word_.IEOR(word_.MASKL(1))}; }
Relation Compare(const Real &) const;
ValueWithRealFlags<Real> Add(
const Real &, Rounding rounding = defaultRounding) const;
ValueWithRealFlags<Real> Subtract(
const Real &y, Rounding rounding = defaultRounding) const {
return Add(y.Negate(), rounding);
}
ValueWithRealFlags<Real> Multiply(
const Real &, Rounding rounding = defaultRounding) const;
ValueWithRealFlags<Real> Divide(
const Real &, Rounding rounding = defaultRounding) const;
ValueWithRealFlags<Real> SQRT(Rounding rounding = defaultRounding) const;
// HYPOT(x,y)=SQRT(x**2 + y**2) computed so as to avoid spurious
// intermediate overflows.
ValueWithRealFlags<Real> HYPOT(
const Real &, Rounding rounding = defaultRounding) const;
template <typename INT> constexpr INT EXPONENT() const {
if (Exponent() == maxExponent) {
return INT::HUGE();
} else {
return {UnbiasedExponent()};
}
}
static constexpr Real EPSILON() {
Real epsilon;
epsilon.Normalize(
false, exponentBias + 1 - binaryPrecision, Fraction::MASKL(1));
return epsilon;
}
static constexpr Real HUGE() {
Real huge;
huge.Normalize(false, maxExponent - 1, Fraction::MASKR(binaryPrecision));
return huge;
}
static constexpr Real TINY() {
Real tiny;
tiny.Normalize(false, 1, Fraction::MASKL(1)); // minimum *normal* number
return tiny;
}
static constexpr int DIGITS{binaryPrecision};
static constexpr int PRECISION{Details::decimalPrecision};
static constexpr int RANGE{Details::decimalRange};
static constexpr int MAXEXPONENT{maxExponent - exponentBias};
static constexpr int MINEXPONENT{2 - exponentBias};
constexpr Real FlushSubnormalToZero() const {
if (IsSubnormal()) {
return Real{};
}
return *this;
}
// TODO: Configurable NotANumber representations
static constexpr Real NotANumber() {
return {Word{maxExponent}
.SHIFTL(significandBits)
.IBSET(significandBits - 1)
.IBSET(significandBits - 2)};
}
static constexpr Real Infinity(bool negative) {
Word infinity{maxExponent};
infinity = infinity.SHIFTL(significandBits);
if (negative) {
infinity = infinity.IBSET(infinity.bits - 1);
}
return {infinity};
}
template <typename INT>
static ValueWithRealFlags<Real> FromInteger(
const INT &n, Rounding rounding = defaultRounding) {
bool isNegative{n.IsNegative()};
INT absN{n};
if (isNegative) {
absN = n.Negate().value; // overflow is safe to ignore
}
int leadz{absN.LEADZ()};
if (leadz >= absN.bits) {
return {}; // all bits zero -> +0.0
}
ValueWithRealFlags<Real> result;
int exponent{exponentBias + absN.bits - leadz - 1};
int bitsNeeded{absN.bits - (leadz + isImplicitMSB)};
int bitsLost{bitsNeeded - significandBits};
if (bitsLost <= 0) {
Fraction fraction{Fraction::ConvertUnsigned(absN).value};
result.flags |= result.value.Normalize(
isNegative, exponent, fraction.SHIFTL(-bitsLost));
} else {
Fraction fraction{Fraction::ConvertUnsigned(absN.SHIFTR(bitsLost)).value};
result.flags |= result.value.Normalize(isNegative, exponent, fraction);
RoundingBits roundingBits{absN, bitsLost};
result.flags |= result.value.Round(rounding, roundingBits);
}
return result;
}
// Conversion to integer in the same real format (AINT(), ANINT())
ValueWithRealFlags<Real> ToWholeNumber(
common::RoundingMode = common::RoundingMode::ToZero) const;
// Conversion to an integer (INT(), NINT(), FLOOR(), CEILING())
template <typename INT>
constexpr ValueWithRealFlags<INT> ToInteger(
common::RoundingMode mode = common::RoundingMode::ToZero) const {
ValueWithRealFlags<INT> result;
if (IsNotANumber()) {
result.flags.set(RealFlag::InvalidArgument);
result.value = result.value.HUGE();
return result;
}
ValueWithRealFlags<Real> intPart{ToWholeNumber(mode)};
result.flags |= intPart.flags;
int exponent{intPart.value.Exponent()};
// shift positive -> left shift, negative -> right shift
int shift{exponent - exponentBias - binaryPrecision + 1};
// Apply any right shift before moving to the result type
auto rshifted{intPart.value.GetFraction().SHIFTR(-shift)};
auto converted{result.value.ConvertUnsigned(rshifted)};
if (converted.overflow) {
result.flags.set(RealFlag::Overflow);
}
result.value = converted.value.SHIFTL(shift);
if (converted.value.CompareUnsigned(result.value.SHIFTR(shift)) !=
Ordering::Equal) {
result.flags.set(RealFlag::Overflow);
}
if (IsSignBitSet()) {
result.value = result.value.Negate().value;
}
if (IsSignBitSet() != result.value.IsNegative()) {
result.flags.set(RealFlag::Overflow);
}
if (result.flags.test(RealFlag::Overflow)) {
result.value =
IsSignBitSet() ? result.value.MASKL(1) : result.value.HUGE();
}
return result;
}
template <typename A>
static ValueWithRealFlags<Real> Convert(
const A &x, Rounding rounding = defaultRounding) {
ValueWithRealFlags<Real> result;
if (x.IsNotANumber()) {
result.flags.set(RealFlag::InvalidArgument);
result.value = NotANumber();
return result;
}
bool isNegative{x.IsNegative()};
A absX{x};
if (isNegative) {
absX = x.Negate();
}
int exponent{exponentBias + x.UnbiasedExponent()};
int bitsLost{A::binaryPrecision - binaryPrecision};
if (exponent < 1) {
bitsLost += 1 - exponent;
exponent = 1;
}
typename A::Fraction xFraction{x.GetFraction()};
if (bitsLost <= 0) {
Fraction fraction{
Fraction::ConvertUnsigned(xFraction).value.SHIFTL(-bitsLost)};
result.flags |= result.value.Normalize(isNegative, exponent, fraction);
} else {
Fraction fraction{
Fraction::ConvertUnsigned(xFraction.SHIFTR(bitsLost)).value};
result.flags |= result.value.Normalize(isNegative, exponent, fraction);
RoundingBits roundingBits{xFraction, bitsLost};
result.flags |= result.value.Round(rounding, roundingBits);
}
return result;
}
constexpr Word RawBits() const { return word_; }
// Extracts "raw" biased exponent field.
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 (!isImplicitMSB) {
return result;
} else {
int exponent{Exponent()};
if (exponent > 0 && exponent < maxExponent) {
return result.IBSET(significandBits);
} else {
return result.IBCLR(significandBits);
}
}
}
// Extracts unbiased exponent value.
// Corrects the exponent value of a subnormal number.
constexpr int UnbiasedExponent() const {
int exponent{Exponent() - exponentBias};
if (IsSubnormal()) {
++exponent;
}
return exponent;
}
static ValueWithRealFlags<Real> Read(
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.
llvm::raw_ostream &AsFortran(
llvm::raw_ostream &, int kind, bool minimal = false) const;
private:
using Significand = Integer<significandBits>; // no implicit bit
constexpr Significand GetSignificand() const {
return Significand::ConvertUnsigned(word_).value;
}
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;
if (forDivide) {
exponent += exponentBias - yExponent;
} else {
exponent += yExponent - exponentBias + 1;
}
return exponent;
}
static constexpr bool NextQuotientBit(
Fraction &top, bool &msb, const Fraction &divisor) {
bool greaterOrEqual{msb || top.CompareUnsigned(divisor) != Ordering::Less};
if (greaterOrEqual) {
top = top.SubtractSigned(divisor).value;
}
auto doubled{top.AddUnsigned(top)};
top = doubled.value;
msb = doubled.carry;
return greaterOrEqual;
}
// Normalizes and marshals the fields of a floating-point number in place.
// 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, 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, int exponent, const Fraction &, Rounding, RoundingBits,
bool multiply = false);
Word word_{}; // an Integer<>
};
extern template class Real<Integer<16>, 11>; // IEEE half format
extern template class Real<Integer<16>, 8>; // the "other" half format
extern template class Real<Integer<32>, 24>; // IEEE single
extern template class Real<Integer<64>, 53>; // IEEE double
extern template class Real<Integer<80>, 64>; // 80387 extended precision
extern template class Real<Integer<128>, 113>; // IEEE quad
// N.B. No "double-double" support.
} // namespace Fortran::evaluate::value
#endif // FORTRAN_EVALUATE_REAL_H_