xref: /llvm-project/flang/include/flang/Evaluate/real.h (revision fc97d2e68b03bc2979395e84b645e5b3ba35aecd) 
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/target.h"
17 #include <cinttypes>
18 #include <limits>
19 #include <string>
20 
21 // Some environments, viz. glibc 2.17 and *BSD, allow the macro HUGE
22 // to leak out of <math.h>.
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> class Real {
39 public:
40   using Word = WORD;
41   static constexpr int binaryPrecision{PREC};
42   static constexpr common::RealCharacteristics realChars{PREC};
43   static constexpr int exponentBias{realChars.exponentBias};
44   static constexpr int exponentBits{realChars.exponentBits};
45   static constexpr int isImplicitMSB{realChars.isImplicitMSB};
46   static constexpr int maxExponent{realChars.maxExponent};
47   static constexpr int significandBits{realChars.significandBits};
48 
49   static constexpr int bits{Word::bits};
50   static_assert(bits >= realChars.bits);
51   using Fraction = Integer<binaryPrecision>; // all bits made explicit
52 
53   template <typename W, int P> friend class Real;
54 
55   constexpr Real() {} // +0.0
56   constexpr Real(const Real &) = default;
57   constexpr Real(Real &&) = default;
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   constexpr bool IsSignBitSet() const { return word_.BTEST(bits - 1); }
67   constexpr bool IsNegative() const {
68     return !IsNotANumber() && IsSignBitSet();
69   }
70   constexpr bool IsNotANumber() const {
71     auto expo{Exponent()};
72     auto sig{GetSignificand()};
73     if constexpr (bits == 80) { // x87
74       // 7FFF8000000000000000 is Infinity, not NaN, on 80387 & later.
75       if (expo == maxExponent) {
76         return sig != Significand{}.IBSET(63);
77       } else {
78         return expo != 0 && !sig.BTEST(63);
79       }
80     } else {
81       return expo == maxExponent && !sig.IsZero();
82     }
83   }
84   constexpr bool IsQuietNaN() const {
85     auto expo{Exponent()};
86     auto sig{GetSignificand()};
87     if constexpr (bits == 80) { // x87
88       if (expo == maxExponent) {
89         return sig.IBITS(62, 2) == 3;
90       } else {
91         return expo != 0 && !sig.BTEST(63);
92       }
93     } else {
94       return expo == maxExponent && sig.BTEST(significandBits - 1);
95     }
96   }
97   constexpr bool IsSignalingNaN() const {
98     auto expo{Exponent()};
99     auto sig{GetSignificand()};
100     if constexpr (bits == 80) { // x87
101       return expo == maxExponent && sig != Significand{}.IBSET(63) &&
102           sig.IBITS(62, 2) != 3;
103     } else {
104       return expo == maxExponent && !sig.IsZero() &&
105           !sig.BTEST(significandBits - 1);
106     }
107   }
108   constexpr bool IsInfinite() const {
109     if constexpr (bits == 80) { // x87
110       // 7FFF8000000000000000 is Infinity, not NaN, on 80387 & later.
111       return Exponent() == maxExponent &&
112           GetSignificand() == Significand{}.IBSET(63);
113     } else {
114       return Exponent() == maxExponent && GetSignificand().IsZero();
115     }
116   }
117   constexpr bool IsFinite() const {
118     auto expo{Exponent()};
119     if constexpr (bits == 80) { // x87
120       return expo != maxExponent && (expo == 0 || GetSignificand().BTEST(63));
121     } else {
122       return expo != maxExponent;
123     }
124   }
125   constexpr bool IsZero() const {
126     return Exponent() == 0 && GetSignificand().IsZero();
127   }
128   constexpr bool IsSubnormal() const {
129     return Exponent() == 0 && !GetSignificand().IsZero();
130   }
131   constexpr bool IsNormal() const {
132     return !(IsInfinite() || IsNotANumber() || IsSubnormal());
133   }
134 
135   constexpr Real ABS() const { // non-arithmetic, no flags returned
136     return {word_.IBCLR(bits - 1)};
137   }
138   constexpr Real SetSign(bool toNegative) const { // non-arithmetic
139     if (toNegative) {
140       return {word_.IBSET(bits - 1)};
141     } else {
142       return ABS();
143     }
144   }
145   constexpr Real SIGN(const Real &x) const { return SetSign(x.IsSignBitSet()); }
146 
147   constexpr Real Negate() const { return {word_.IEOR(word_.MASKL(1))}; }
148 
149   Relation Compare(const Real &) const;
150   ValueWithRealFlags<Real> Add(const Real &,
151       Rounding rounding = TargetCharacteristics::defaultRounding) const;
152   ValueWithRealFlags<Real> Subtract(const Real &y,
153       Rounding rounding = TargetCharacteristics::defaultRounding) const {
154     return Add(y.Negate(), rounding);
155   }
156   ValueWithRealFlags<Real> Multiply(const Real &,
157       Rounding rounding = TargetCharacteristics::defaultRounding) const;
158   ValueWithRealFlags<Real> Divide(const Real &,
159       Rounding rounding = TargetCharacteristics::defaultRounding) const;
160 
161   ValueWithRealFlags<Real> SQRT(
162       Rounding rounding = TargetCharacteristics::defaultRounding) const;
163   // NEAREST(), IEEE_NEXT_AFTER(), IEEE_NEXT_UP(), and IEEE_NEXT_DOWN()
164   ValueWithRealFlags<Real> NEAREST(bool upward) const;
165   // HYPOT(x,y)=SQRT(x**2 + y**2) computed so as to avoid spurious
166   // intermediate overflows.
167   ValueWithRealFlags<Real> HYPOT(const Real &,
168       Rounding rounding = TargetCharacteristics::defaultRounding) const;
169   // DIM(X,Y) = MAX(X-Y, 0)
170   ValueWithRealFlags<Real> DIM(const Real &,
171       Rounding rounding = TargetCharacteristics::defaultRounding) const;
172   // MOD(x,y) = x - AINT(x/y)*y (in the standard)
173   // MODULO(x,y) = x - FLOOR(x/y)*y (in the standard)
174   ValueWithRealFlags<Real> MOD(const Real &,
175       Rounding rounding = TargetCharacteristics::defaultRounding) const;
176   ValueWithRealFlags<Real> MODULO(const Real &,
177       Rounding rounding = TargetCharacteristics::defaultRounding) const;
178 
179   template <typename INT> constexpr INT EXPONENT() const {
180     if (Exponent() == maxExponent) {
181       return INT::HUGE();
182     } else if (IsZero()) {
183       return {0};
184     } else {
185       return {UnbiasedExponent() + 1};
186     }
187   }
188 
189   static constexpr Real EPSILON() {
190     Real epsilon;
191     epsilon.Normalize(
192         false, exponentBias + 1 - binaryPrecision, Fraction::MASKL(1));
193     return epsilon;
194   }
195   static constexpr Real HUGE() {
196     Real huge;
197     huge.Normalize(false, maxExponent - 1, Fraction::MASKR(binaryPrecision));
198     return huge;
199   }
200   static constexpr Real TINY() {
201     Real tiny;
202     tiny.Normalize(false, 1, Fraction::MASKL(1)); // minimum *normal* number
203     return tiny;
204   }
205 
206   static constexpr int DIGITS{binaryPrecision};
207   static constexpr int PRECISION{realChars.decimalPrecision};
208   static constexpr int RANGE{realChars.decimalRange};
209   static constexpr int MAXEXPONENT{maxExponent - exponentBias};
210   static constexpr int MINEXPONENT{2 - exponentBias};
211   Real RRSPACING() const;
212   Real SPACING() const;
213   Real SET_EXPONENT(std::int64_t) const;
214   Real FRACTION() const;
215 
216   // SCALE(); also known as IEEE_SCALB and (in IEEE-754 '08) ScaleB.
217   template <typename INT>
218   ValueWithRealFlags<Real> SCALE(const INT &by,
219       Rounding rounding = TargetCharacteristics::defaultRounding) const {
220     // Normalize a fraction with just its LSB set and then multiply.
221     // (Set the LSB, not the MSB, in case the scale factor needs to
222     //  be subnormal.)
223     constexpr auto adjust{exponentBias + binaryPrecision - 1};
224     constexpr auto maxCoeffExpo{maxExponent + binaryPrecision - 1};
225     auto expo{adjust + by.ToInt64()};
226     RealFlags flags;
227     int rMask{1};
228     if (IsZero()) {
229       expo = exponentBias; // ignore by, don't overflow
230     } else if (expo > maxCoeffExpo) {
231       if (Exponent() < exponentBias) {
232         // Must implement with two multiplications
233         return SCALE(INT{exponentBias})
234             .value.SCALE(by.SubtractSigned(INT{exponentBias}).value, rounding);
235       } else { // overflow
236         expo = maxCoeffExpo;
237       }
238     } else if (expo < 0) {
239       if (Exponent() > exponentBias) {
240         // Must implement with two multiplications
241         return SCALE(INT{-exponentBias})
242             .value.SCALE(by.AddSigned(INT{exponentBias}).value, rounding);
243       } else { // underflow to zero
244         expo = 0;
245         rMask = 0;
246         flags.set(RealFlag::Underflow);
247       }
248     }
249     Real twoPow;
250     flags |=
251         twoPow.Normalize(false, static_cast<int>(expo), Fraction::MASKR(rMask));
252     ValueWithRealFlags<Real> result{Multiply(twoPow, rounding)};
253     result.flags |= flags;
254     return result;
255   }
256 
257   constexpr Real FlushSubnormalToZero() const {
258     if (IsSubnormal()) {
259       return Real{};
260     }
261     return *this;
262   }
263 
264   // TODO: Configurable NotANumber representations
265   static constexpr Real NotANumber() {
266     return {Word{maxExponent}
267                 .SHIFTL(significandBits)
268                 .IBSET(significandBits - 1)
269                 .IBSET(significandBits - 2)};
270   }
271 
272   static constexpr Real PositiveZero() { return Real{}; }
273 
274   static constexpr Real NegativeZero() { return {Word{}.MASKL(1)}; }
275 
276   static constexpr Real Infinity(bool negative) {
277     Word infinity{maxExponent};
278     infinity = infinity.SHIFTL(significandBits);
279     if (negative) {
280       infinity = infinity.IBSET(infinity.bits - 1);
281     }
282     if constexpr (bits == 80) { // x87
283       // 7FFF8000000000000000 is Infinity, not NaN, on 80387 & later.
284       infinity = infinity.IBSET(63);
285     }
286     return {infinity};
287   }
288 
289   template <typename INT>
290   static ValueWithRealFlags<Real> FromInteger(const INT &n,
291       bool isUnsigned = false,
292       Rounding rounding = TargetCharacteristics::defaultRounding) {
293     bool isNegative{!isUnsigned && n.IsNegative()};
294     INT absN{n};
295     if (isNegative) {
296       absN = n.Negate().value; // overflow is safe to ignore
297     }
298     int leadz{absN.LEADZ()};
299     if (leadz >= absN.bits) {
300       return {}; // all bits zero -> +0.0
301     }
302     ValueWithRealFlags<Real> result;
303     int exponent{exponentBias + absN.bits - leadz - 1};
304     int bitsNeeded{absN.bits - (leadz + isImplicitMSB)};
305     int bitsLost{bitsNeeded - significandBits};
306     if (bitsLost <= 0) {
307       Fraction fraction{Fraction::ConvertUnsigned(absN).value};
308       result.flags |= result.value.Normalize(
309           isNegative, exponent, fraction.SHIFTL(-bitsLost));
310     } else {
311       Fraction fraction{Fraction::ConvertUnsigned(absN.SHIFTR(bitsLost)).value};
312       result.flags |= result.value.Normalize(isNegative, exponent, fraction);
313       RoundingBits roundingBits{absN, bitsLost};
314       result.flags |= result.value.Round(rounding, roundingBits);
315     }
316     return result;
317   }
318 
319   // Conversion to integer in the same real format (AINT(), ANINT())
320   ValueWithRealFlags<Real> ToWholeNumber(
321       common::RoundingMode = common::RoundingMode::ToZero) const;
322 
323   // Conversion to an integer (INT(), NINT(), FLOOR(), CEILING())
324   template <typename INT>
325   constexpr ValueWithRealFlags<INT> ToInteger(
326       common::RoundingMode mode = common::RoundingMode::ToZero) const {
327     ValueWithRealFlags<INT> result;
328     if (IsNotANumber()) {
329       result.flags.set(RealFlag::InvalidArgument);
330       result.value = result.value.HUGE();
331       return result;
332     }
333     ValueWithRealFlags<Real> intPart{ToWholeNumber(mode)};
334     result.flags |= intPart.flags;
335     int exponent{intPart.value.Exponent()};
336     // shift positive -> left shift, negative -> right shift
337     int shift{exponent - exponentBias - binaryPrecision + 1};
338     // Apply any right shift before moving to the result type
339     auto rshifted{intPart.value.GetFraction().SHIFTR(-shift)};
340     auto converted{result.value.ConvertUnsigned(rshifted)};
341     if (converted.overflow) {
342       result.flags.set(RealFlag::Overflow);
343     }
344     result.value = converted.value.SHIFTL(shift);
345     if (converted.value.CompareUnsigned(result.value.SHIFTR(shift)) !=
346         Ordering::Equal) {
347       result.flags.set(RealFlag::Overflow);
348     }
349     if (IsSignBitSet()) {
350       result.value = result.value.Negate().value;
351     }
352     if (!result.value.IsZero()) {
353       if (IsSignBitSet() != result.value.IsNegative()) {
354         result.flags.set(RealFlag::Overflow);
355       }
356     }
357     if (result.flags.test(RealFlag::Overflow)) {
358       result.value =
359           IsSignBitSet() ? result.value.MASKL(1) : result.value.HUGE();
360     }
361     return result;
362   }
363 
364   template <typename A>
365   static ValueWithRealFlags<Real> Convert(
366       const A &x, Rounding rounding = TargetCharacteristics::defaultRounding) {
367     ValueWithRealFlags<Real> result;
368     if (x.IsNotANumber()) {
369       result.flags.set(RealFlag::InvalidArgument);
370       result.value = NotANumber();
371       return result;
372     }
373     bool isNegative{x.IsNegative()};
374     if (x.IsInfinite()) {
375       result.value = Infinity(isNegative);
376       return result;
377     }
378     A absX{x};
379     if (isNegative) {
380       absX = x.Negate();
381     }
382     int exponent{exponentBias + x.UnbiasedExponent()};
383     int bitsLost{A::binaryPrecision - binaryPrecision};
384     if (exponent < 1) {
385       bitsLost += 1 - exponent;
386       exponent = 1;
387     }
388     typename A::Fraction xFraction{x.GetFraction()};
389     if (bitsLost <= 0) {
390       Fraction fraction{
391           Fraction::ConvertUnsigned(xFraction).value.SHIFTL(-bitsLost)};
392       result.flags |= result.value.Normalize(isNegative, exponent, fraction);
393     } else {
394       Fraction fraction{
395           Fraction::ConvertUnsigned(xFraction.SHIFTR(bitsLost)).value};
396       result.flags |= result.value.Normalize(isNegative, exponent, fraction);
397       RoundingBits roundingBits{xFraction, bitsLost};
398       result.flags |= result.value.Round(rounding, roundingBits);
399     }
400     return result;
401   }
402 
403   constexpr Word RawBits() const { return word_; }
404 
405   // Extracts "raw" biased exponent field.
406   constexpr int Exponent() const {
407     return word_.IBITS(significandBits, exponentBits).ToUInt64();
408   }
409 
410   // Extracts the fraction; any implied bit is made explicit.
411   constexpr Fraction GetFraction() const {
412     Fraction result{Fraction::ConvertUnsigned(word_).value};
413     if constexpr (!isImplicitMSB) {
414       return result;
415     } else {
416       int exponent{Exponent()};
417       if (exponent > 0 && exponent < maxExponent) {
418         return result.IBSET(significandBits);
419       } else {
420         return result.IBCLR(significandBits);
421       }
422     }
423   }
424 
425   // Extracts unbiased exponent value.
426   // Corrects the exponent value of a subnormal number.
427   // Note that the result is one less than the EXPONENT intrinsic;
428   // UnbiasedExponent(1.0) is 0, not 1.
429   constexpr int UnbiasedExponent() const {
430     int exponent{Exponent() - exponentBias};
431     if (IsSubnormal()) {
432       ++exponent;
433     }
434     return exponent;
435   }
436 
437   static ValueWithRealFlags<Real> Read(const char *&,
438       Rounding rounding = TargetCharacteristics::defaultRounding);
439   std::string DumpHexadecimal() const;
440 
441   // Emits a character representation for an equivalent Fortran constant
442   // or parenthesized constant expression that produces this value.
443   llvm::raw_ostream &AsFortran(
444       llvm::raw_ostream &, int kind, bool minimal = false) const;
445 
446 private:
447   using Significand = Integer<significandBits>; // no implicit bit
448 
449   constexpr Significand GetSignificand() const {
450     return Significand::ConvertUnsigned(word_).value;
451   }
452 
453   constexpr int CombineExponents(const Real &y, bool forDivide) const {
454     int exponent = Exponent(), yExponent = y.Exponent();
455     // A zero exponent field value has the same weight as 1.
456     exponent += !exponent;
457     yExponent += !yExponent;
458     if (forDivide) {
459       exponent += exponentBias - yExponent;
460     } else {
461       exponent += yExponent - exponentBias + 1;
462     }
463     return exponent;
464   }
465 
466   static constexpr bool NextQuotientBit(
467       Fraction &top, bool &msb, const Fraction &divisor) {
468     bool greaterOrEqual{msb || top.CompareUnsigned(divisor) != Ordering::Less};
469     if (greaterOrEqual) {
470       top = top.SubtractSigned(divisor).value;
471     }
472     auto doubled{top.AddUnsigned(top)};
473     top = doubled.value;
474     msb = doubled.carry;
475     return greaterOrEqual;
476   }
477 
478   // Normalizes and marshals the fields of a floating-point number in place.
479   // The value is a number, and a zero fraction means a zero value (i.e.,
480   // a maximal exponent and zero fraction doesn't signify infinity, although
481   // this member function will detect overflow and encode infinities).
482   RealFlags Normalize(bool negative, int exponent, const Fraction &fraction,
483       Rounding rounding = TargetCharacteristics::defaultRounding,
484       RoundingBits *roundingBits = nullptr);
485 
486   // Rounds a result, if necessary, in place.
487   RealFlags Round(Rounding, const RoundingBits &, bool multiply = false);
488 
489   static void NormalizeAndRound(ValueWithRealFlags<Real> &result,
490       bool isNegative, int exponent, const Fraction &, Rounding, RoundingBits,
491       bool multiply = false);
492 
493   Word word_{}; // an Integer<>
494 };
495 
496 extern template class Real<Integer<16>, 11>; // IEEE half format
497 extern template class Real<Integer<16>, 8>; // the "other" half format
498 extern template class Real<Integer<32>, 24>; // IEEE single
499 extern template class Real<Integer<64>, 53>; // IEEE double
500 extern template class Real<X87IntegerContainer, 64>; // 80387 extended precision
501 extern template class Real<Integer<128>, 113>; // IEEE quad
502 // N.B. No "double-double" support.
503 } // namespace Fortran::evaluate::value
504 #endif // FORTRAN_EVALUATE_REAL_H_
505