• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 //===-- include/flang/Evaluate/real.h ---------------------------*- C++ -*-===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 
9 #ifndef FORTRAN_EVALUATE_REAL_H_
10 #define FORTRAN_EVALUATE_REAL_H_
11 
12 #include "formatting.h"
13 #include "integer.h"
14 #include "rounding-bits.h"
15 #include "flang/Common/real.h"
16 #include "flang/Evaluate/common.h"
17 #include <cinttypes>
18 #include <limits>
19 #include <string>
20 
21 // Some environments, viz. clang on Darwin, allow the macro HUGE
22 // to leak out of <math.h> even when it is never directly included.
23 #undef HUGE
24 
25 namespace llvm {
26 class raw_ostream;
27 }
28 namespace Fortran::evaluate::value {
29 
30 // LOG10(2.)*1E12
31 static constexpr std::int64_t ScaledLogBaseTenOfTwo{301029995664};
32 
33 // Models IEEE binary floating-point numbers (IEEE 754-2008,
34 // ISO/IEC/IEEE 60559.2011).  The first argument to this
35 // class template must be (or look like) an instance of Integer<>;
36 // the second specifies the number of effective bits (binary precision)
37 // in the fraction.
38 template <typename WORD, int PREC>
39 class Real : public common::RealDetails<PREC> {
40 public:
41   using Word = WORD;
42   static constexpr int binaryPrecision{PREC};
43   using Details = common::RealDetails<PREC>;
44   using Details::exponentBias;
45   using Details::exponentBits;
46   using Details::isImplicitMSB;
47   using Details::maxExponent;
48   using Details::significandBits;
49 
50   static constexpr int bits{Word::bits};
51   static_assert(bits >= Details::bits);
52   using Fraction = Integer<binaryPrecision>; // all bits made explicit
53 
54   template <typename W, int P> friend class Real;
55 
Real()56   constexpr Real() {} // +0.0
57   constexpr Real(const Real &) = default;
Real(const Word & bits)58   constexpr Real(const Word &bits) : word_{bits} {}
59   constexpr Real &operator=(const Real &) = default;
60   constexpr Real &operator=(Real &&) = default;
61 
62   constexpr bool operator==(const Real &that) const {
63     return word_ == that.word_;
64   }
65 
66   // TODO: DIM, MAX, MIN, DPROD, FRACTION,
67   // INT/NINT, NEAREST, OUT_OF_RANGE,
68   // RRSPACING/SPACING, SCALE, SET_EXPONENT
69 
IsSignBitSet()70   constexpr bool IsSignBitSet() const { return word_.BTEST(bits - 1); }
IsNegative()71   constexpr bool IsNegative() const {
72     return !IsNotANumber() && IsSignBitSet();
73   }
IsNotANumber()74   constexpr bool IsNotANumber() const {
75     return Exponent() == maxExponent && !GetSignificand().IsZero();
76   }
IsQuietNaN()77   constexpr bool IsQuietNaN() const {
78     return Exponent() == maxExponent &&
79         GetSignificand().BTEST(significandBits - 1);
80   }
IsSignalingNaN()81   constexpr bool IsSignalingNaN() const {
82     return IsNotANumber() && !GetSignificand().BTEST(significandBits - 1);
83   }
IsInfinite()84   constexpr bool IsInfinite() const {
85     return Exponent() == maxExponent && GetSignificand().IsZero();
86   }
IsZero()87   constexpr bool IsZero() const {
88     return Exponent() == 0 && GetSignificand().IsZero();
89   }
IsSubnormal()90   constexpr bool IsSubnormal() const {
91     return Exponent() == 0 && !GetSignificand().IsZero();
92   }
93 
ABS()94   constexpr Real ABS() const { // non-arithmetic, no flags returned
95     return {word_.IBCLR(bits - 1)};
96   }
SetSign(bool toNegative)97   constexpr Real SetSign(bool toNegative) const { // non-arithmetic
98     if (toNegative) {
99       return {word_.IBSET(bits - 1)};
100     } else {
101       return ABS();
102     }
103   }
SIGN(const Real & x)104   constexpr Real SIGN(const Real &x) const { return SetSign(x.IsSignBitSet()); }
105 
Negate()106   constexpr Real Negate() const { return {word_.IEOR(word_.MASKL(1))}; }
107 
108   Relation Compare(const Real &) const;
109   ValueWithRealFlags<Real> Add(
110       const Real &, Rounding rounding = defaultRounding) const;
111   ValueWithRealFlags<Real> Subtract(
112       const Real &y, Rounding rounding = defaultRounding) const {
113     return Add(y.Negate(), rounding);
114   }
115   ValueWithRealFlags<Real> Multiply(
116       const Real &, Rounding rounding = defaultRounding) const;
117   ValueWithRealFlags<Real> Divide(
118       const Real &, Rounding rounding = defaultRounding) const;
119 
120   // SQRT(x**2 + y**2) but computed so as to avoid spurious overflow
121   // TODO: needed for CABS
122   ValueWithRealFlags<Real> HYPOT(
123       const Real &, Rounding rounding = defaultRounding) const;
124 
EXPONENT()125   template <typename INT> constexpr INT EXPONENT() const {
126     if (Exponent() == maxExponent) {
127       return INT::HUGE();
128     } else {
129       return {UnbiasedExponent()};
130     }
131   }
132 
EPSILON()133   static constexpr Real EPSILON() {
134     Real epsilon;
135     epsilon.Normalize(
136         false, exponentBias - binaryPrecision, Fraction::MASKL(1));
137     return epsilon;
138   }
HUGE()139   static constexpr Real HUGE() {
140     Real huge;
141     huge.Normalize(false, maxExponent - 1, Fraction::MASKR(binaryPrecision));
142     return huge;
143   }
TINY()144   static constexpr Real TINY() {
145     Real tiny;
146     tiny.Normalize(false, 1, Fraction::MASKL(1)); // minimum *normal* number
147     return tiny;
148   }
149 
150   static constexpr int DIGITS{binaryPrecision};
151   static constexpr int PRECISION{Details::decimalPrecision};
152   static constexpr int RANGE{Details::decimalRange};
153   static constexpr int MAXEXPONENT{maxExponent - 1 - exponentBias};
154   static constexpr int MINEXPONENT{1 - exponentBias};
155 
FlushSubnormalToZero()156   constexpr Real FlushSubnormalToZero() const {
157     if (IsSubnormal()) {
158       return Real{};
159     }
160     return *this;
161   }
162 
163   // TODO: Configurable NotANumber representations
NotANumber()164   static constexpr Real NotANumber() {
165     return {Word{maxExponent}
166                 .SHIFTL(significandBits)
167                 .IBSET(significandBits - 1)
168                 .IBSET(significandBits - 2)};
169   }
170 
Infinity(bool negative)171   static constexpr Real Infinity(bool negative) {
172     Word infinity{maxExponent};
173     infinity = infinity.SHIFTL(significandBits);
174     if (negative) {
175       infinity = infinity.IBSET(infinity.bits - 1);
176     }
177     return {infinity};
178   }
179 
180   template <typename INT>
181   static ValueWithRealFlags<Real> FromInteger(
182       const INT &n, Rounding rounding = defaultRounding) {
183     bool isNegative{n.IsNegative()};
184     INT absN{n};
185     if (isNegative) {
186       absN = n.Negate().value; // overflow is safe to ignore
187     }
188     int leadz{absN.LEADZ()};
189     if (leadz >= absN.bits) {
190       return {}; // all bits zero -> +0.0
191     }
192     ValueWithRealFlags<Real> result;
193     int exponent{exponentBias + absN.bits - leadz - 1};
194     int bitsNeeded{absN.bits - (leadz + isImplicitMSB)};
195     int bitsLost{bitsNeeded - significandBits};
196     if (bitsLost <= 0) {
197       Fraction fraction{Fraction::ConvertUnsigned(absN).value};
198       result.flags |= result.value.Normalize(
199           isNegative, exponent, fraction.SHIFTL(-bitsLost));
200     } else {
201       Fraction fraction{Fraction::ConvertUnsigned(absN.SHIFTR(bitsLost)).value};
202       result.flags |= result.value.Normalize(isNegative, exponent, fraction);
203       RoundingBits roundingBits{absN, bitsLost};
204       result.flags |= result.value.Round(rounding, roundingBits);
205     }
206     return result;
207   }
208 
209   // Conversion to integer in the same real format (AINT(), ANINT())
210   ValueWithRealFlags<Real> ToWholeNumber(
211       common::RoundingMode = common::RoundingMode::ToZero) const;
212 
213   // Conversion to an integer (INT(), NINT(), FLOOR(), CEILING())
214   template <typename INT>
215   constexpr ValueWithRealFlags<INT> ToInteger(
216       common::RoundingMode mode = common::RoundingMode::ToZero) const {
217     ValueWithRealFlags<INT> result;
218     if (IsNotANumber()) {
219       result.flags.set(RealFlag::InvalidArgument);
220       result.value = result.value.HUGE();
221       return result;
222     }
223     ValueWithRealFlags<Real> intPart{ToWholeNumber(mode)};
224     int exponent{intPart.value.Exponent()};
225     result.flags.set(
226         RealFlag::Overflow, exponent >= exponentBias + result.value.bits);
227     result.flags |= intPart.flags;
228     int shift{
229         exponent - exponentBias - binaryPrecision + 1}; // positive -> left
230     result.value =
231         result.value.ConvertUnsigned(intPart.value.GetFraction().SHIFTR(-shift))
232             .value.SHIFTL(shift);
233     if (IsSignBitSet()) {
234       auto negated{result.value.Negate()};
235       result.value = negated.value;
236       if (negated.overflow) {
237         result.flags.set(RealFlag::Overflow);
238       }
239     }
240     if (result.flags.test(RealFlag::Overflow)) {
241       result.value =
242           IsSignBitSet() ? result.value.MASKL(1) : result.value.HUGE();
243     }
244     return result;
245   }
246 
247   template <typename A>
248   static ValueWithRealFlags<Real> Convert(
249       const A &x, Rounding rounding = defaultRounding) {
250     ValueWithRealFlags<Real> result;
251     if (x.IsNotANumber()) {
252       result.flags.set(RealFlag::InvalidArgument);
253       result.value = NotANumber();
254       return result;
255     }
256     bool isNegative{x.IsNegative()};
257     A absX{x};
258     if (isNegative) {
259       absX = x.Negate();
260     }
261     int exponent{exponentBias + x.UnbiasedExponent()};
262     int bitsLost{A::binaryPrecision - binaryPrecision};
263     if (exponent < 1) {
264       bitsLost += 1 - exponent;
265       exponent = 1;
266     }
267     typename A::Fraction xFraction{x.GetFraction()};
268     if (bitsLost <= 0) {
269       Fraction fraction{
270           Fraction::ConvertUnsigned(xFraction).value.SHIFTL(-bitsLost)};
271       result.flags |= result.value.Normalize(isNegative, exponent, fraction);
272     } else {
273       Fraction fraction{
274           Fraction::ConvertUnsigned(xFraction.SHIFTR(bitsLost)).value};
275       result.flags |= result.value.Normalize(isNegative, exponent, fraction);
276       RoundingBits roundingBits{xFraction, bitsLost};
277       result.flags |= result.value.Round(rounding, roundingBits);
278     }
279     return result;
280   }
281 
RawBits()282   constexpr Word RawBits() const { return word_; }
283 
284   // Extracts "raw" biased exponent field.
Exponent()285   constexpr int Exponent() const {
286     return word_.IBITS(significandBits, exponentBits).ToUInt64();
287   }
288 
289   // Extracts the fraction; any implied bit is made explicit.
GetFraction()290   constexpr Fraction GetFraction() const {
291     Fraction result{Fraction::ConvertUnsigned(word_).value};
292     if constexpr (!isImplicitMSB) {
293       return result;
294     } else {
295       int exponent{Exponent()};
296       if (exponent > 0 && exponent < maxExponent) {
297         return result.IBSET(significandBits);
298       } else {
299         return result.IBCLR(significandBits);
300       }
301     }
302   }
303 
304   // Extracts unbiased exponent value.
305   // Corrects the exponent value of a subnormal number.
UnbiasedExponent()306   constexpr int UnbiasedExponent() const {
307     int exponent{Exponent() - exponentBias};
308     if (IsSubnormal()) {
309       ++exponent;
310     }
311     return exponent;
312   }
313 
314   static ValueWithRealFlags<Real> Read(
315       const char *&, Rounding rounding = defaultRounding);
316   std::string DumpHexadecimal() const;
317 
318   // Emits a character representation for an equivalent Fortran constant
319   // or parenthesized constant expression that produces this value.
320   llvm::raw_ostream &AsFortran(
321       llvm::raw_ostream &, int kind, bool minimal = false) const;
322 
323 private:
324   using Significand = Integer<significandBits>; // no implicit bit
325 
GetSignificand()326   constexpr Significand GetSignificand() const {
327     return Significand::ConvertUnsigned(word_).value;
328   }
329 
CombineExponents(const Real & y,bool forDivide)330   constexpr int CombineExponents(const Real &y, bool forDivide) const {
331     int exponent = Exponent(), yExponent = y.Exponent();
332     // A zero exponent field value has the same weight as 1.
333     exponent += !exponent;
334     yExponent += !yExponent;
335     if (forDivide) {
336       exponent += exponentBias - yExponent;
337     } else {
338       exponent += yExponent - exponentBias + 1;
339     }
340     return exponent;
341   }
342 
NextQuotientBit(Fraction & top,bool & msb,const Fraction & divisor)343   static constexpr bool NextQuotientBit(
344       Fraction &top, bool &msb, const Fraction &divisor) {
345     bool greaterOrEqual{msb || top.CompareUnsigned(divisor) != Ordering::Less};
346     if (greaterOrEqual) {
347       top = top.SubtractSigned(divisor).value;
348     }
349     auto doubled{top.AddUnsigned(top)};
350     top = doubled.value;
351     msb = doubled.carry;
352     return greaterOrEqual;
353   }
354 
355   // Normalizes and marshals the fields of a floating-point number in place.
356   // The value is a number, and a zero fraction means a zero value (i.e.,
357   // a maximal exponent and zero fraction doesn't signify infinity, although
358   // this member function will detect overflow and encode infinities).
359   RealFlags Normalize(bool negative, int exponent, const Fraction &fraction,
360       Rounding rounding = defaultRounding,
361       RoundingBits *roundingBits = nullptr);
362 
363   // Rounds a result, if necessary, in place.
364   RealFlags Round(Rounding, const RoundingBits &, bool multiply = false);
365 
366   static void NormalizeAndRound(ValueWithRealFlags<Real> &result,
367       bool isNegative, int exponent, const Fraction &, Rounding, RoundingBits,
368       bool multiply = false);
369 
370   Word word_{}; // an Integer<>
371 };
372 
373 extern template class Real<Integer<16>, 11>; // IEEE half format
374 extern template class Real<Integer<16>, 8>; // the "other" half format
375 extern template class Real<Integer<32>, 24>; // IEEE single
376 extern template class Real<Integer<64>, 53>; // IEEE double
377 extern template class Real<Integer<80>, 64>; // 80387 extended precision
378 extern template class Real<Integer<128>, 113>; // IEEE quad
379 // N.B. No "double-double" support.
380 } // namespace Fortran::evaluate::value
381 #endif // FORTRAN_EVALUATE_REAL_H_
382