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