Warning: This file is not a C or C++ file. It does not have highlighting.

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
25namespace llvm {
26class raw_ostream;
27}
28namespace Fortran::evaluate::value {
29
30// LOG10(2.)*1E12
31static 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.
38template <typename WORD, int PREC> class Real {
39public:
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.IBSET(63);
285 }
286 return {infinity};
287 }
288
289 template <typename INT>
290 static ValueWithRealFlags<Real> FromInteger(const INT &n,
291 Rounding rounding = TargetCharacteristics::defaultRounding) {
292 bool isNegative{n.IsNegative()};
293 INT absN{n};
294 if (isNegative) {
295 absN = n.Negate().value; // overflow is safe to ignore
296 }
297 int leadz{absN.LEADZ()};
298 if (leadz >= absN.bits) {
299 return {}; // all bits zero -> +0.0
300 }
301 ValueWithRealFlags<Real> result;
302 int exponent{exponentBias + absN.bits - leadz - 1};
303 int bitsNeeded{absN.bits - (leadz + isImplicitMSB)};
304 int bitsLost{bitsNeeded - significandBits};
305 if (bitsLost <= 0) {
306 Fraction fraction{Fraction::ConvertUnsigned(absN).value};
307 result.flags |= result.value.Normalize(
308 isNegative, exponent, fraction.SHIFTL(-bitsLost));
309 } else {
310 Fraction fraction{Fraction::ConvertUnsigned(absN.SHIFTR(bitsLost)).value};
311 result.flags |= result.value.Normalize(isNegative, exponent, fraction);
312 RoundingBits roundingBits{absN, bitsLost};
313 result.flags |= result.value.Round(rounding, roundingBits);
314 }
315 return result;
316 }
317
318 // Conversion to integer in the same real format (AINT(), ANINT())
319 ValueWithRealFlags<Real> ToWholeNumber(
320 common::RoundingMode = common::RoundingMode::ToZero) const;
321
322 // Conversion to an integer (INT(), NINT(), FLOOR(), CEILING())
323 template <typename INT>
324 constexpr ValueWithRealFlags<INT> ToInteger(
325 common::RoundingMode mode = common::RoundingMode::ToZero) const {
326 ValueWithRealFlags<INT> result;
327 if (IsNotANumber()) {
328 result.flags.set(RealFlag::InvalidArgument);
329 result.value = result.value.HUGE();
330 return result;
331 }
332 ValueWithRealFlags<Real> intPart{ToWholeNumber(mode)};
333 result.flags |= intPart.flags;
334 int exponent{intPart.value.Exponent()};
335 // shift positive -> left shift, negative -> right shift
336 int shift{exponent - exponentBias - binaryPrecision + 1};
337 // Apply any right shift before moving to the result type
338 auto rshifted{intPart.value.GetFraction().SHIFTR(-shift)};
339 auto converted{result.value.ConvertUnsigned(rshifted)};
340 if (converted.overflow) {
341 result.flags.set(RealFlag::Overflow);
342 }
343 result.value = converted.value.SHIFTL(shift);
344 if (converted.value.CompareUnsigned(result.value.SHIFTR(shift)) !=
345 Ordering::Equal) {
346 result.flags.set(RealFlag::Overflow);
347 }
348 if (IsSignBitSet()) {
349 result.value = result.value.Negate().value;
350 }
351 if (!result.value.IsZero()) {
352 if (IsSignBitSet() != result.value.IsNegative()) {
353 result.flags.set(RealFlag::Overflow);
354 }
355 }
356 if (result.flags.test(RealFlag::Overflow)) {
357 result.value =
358 IsSignBitSet() ? result.value.MASKL(1) : result.value.HUGE();
359 }
360 return result;
361 }
362
363 template <typename A>
364 static ValueWithRealFlags<Real> Convert(
365 const A &x, Rounding rounding = TargetCharacteristics::defaultRounding) {
366 ValueWithRealFlags<Real> result;
367 if (x.IsNotANumber()) {
368 result.flags.set(RealFlag::InvalidArgument);
369 result.value = NotANumber();
370 return result;
371 }
372 bool isNegative{x.IsNegative()};
373 if (x.IsInfinite()) {
374 result.value = Infinity(isNegative);
375 return result;
376 }
377 A absX{x};
378 if (isNegative) {
379 absX = x.Negate();
380 }
381 int exponent{exponentBias + x.UnbiasedExponent()};
382 int bitsLost{A::binaryPrecision - binaryPrecision};
383 if (exponent < 1) {
384 bitsLost += 1 - exponent;
385 exponent = 1;
386 }
387 typename A::Fraction xFraction{x.GetFraction()};
388 if (bitsLost <= 0) {
389 Fraction fraction{
390 Fraction::ConvertUnsigned(xFraction).value.SHIFTL(-bitsLost)};
391 result.flags |= result.value.Normalize(isNegative, exponent, fraction);
392 } else {
393 Fraction fraction{
394 Fraction::ConvertUnsigned(xFraction.SHIFTR(bitsLost)).value};
395 result.flags |= result.value.Normalize(isNegative, exponent, fraction);
396 RoundingBits roundingBits{xFraction, bitsLost};
397 result.flags |= result.value.Round(rounding, roundingBits);
398 }
399 return result;
400 }
401
402 constexpr Word RawBits() const { return word_; }
403
404 // Extracts "raw" biased exponent field.
405 constexpr int Exponent() const {
406 return word_.IBITS(significandBits, exponentBits).ToUInt64();
407 }
408
409 // Extracts the fraction; any implied bit is made explicit.
410 constexpr Fraction GetFraction() const {
411 Fraction result{Fraction::ConvertUnsigned(word_).value};
412 if constexpr (!isImplicitMSB) {
413 return result;
414 } else {
415 int exponent{Exponent()};
416 if (exponent > 0 && exponent < maxExponent) {
417 return result.IBSET(significandBits);
418 } else {
419 return result.IBCLR(significandBits);
420 }
421 }
422 }
423
424 // Extracts unbiased exponent value.
425 // Corrects the exponent value of a subnormal number.
426 // Note that the result is one less than the EXPONENT intrinsic;
427 // UnbiasedExponent(1.0) is 0, not 1.
428 constexpr int UnbiasedExponent() const {
429 int exponent{Exponent() - exponentBias};
430 if (IsSubnormal()) {
431 ++exponent;
432 }
433 return exponent;
434 }
435
436 static ValueWithRealFlags<Real> Read(const char *&,
437 Rounding rounding = TargetCharacteristics::defaultRounding);
438 std::string DumpHexadecimal() const;
439
440 // Emits a character representation for an equivalent Fortran constant
441 // or parenthesized constant expression that produces this value.
442 llvm::raw_ostream &AsFortran(
443 llvm::raw_ostream &, int kind, bool minimal = false) const;
444
445private:
446 using Significand = Integer<significandBits>; // no implicit bit
447
448 constexpr Significand GetSignificand() const {
449 return Significand::ConvertUnsigned(word_).value;
450 }
451
452 constexpr int CombineExponents(const Real &y, bool forDivide) const {
453 int exponent = Exponent(), yExponent = y.Exponent();
454 // A zero exponent field value has the same weight as 1.
455 exponent += !exponent;
456 yExponent += !yExponent;
457 if (forDivide) {
458 exponent += exponentBias - yExponent;
459 } else {
460 exponent += yExponent - exponentBias + 1;
461 }
462 return exponent;
463 }
464
465 static constexpr bool NextQuotientBit(
466 Fraction &top, bool &msb, const Fraction &divisor) {
467 bool greaterOrEqual{msb || top.CompareUnsigned(divisor) != Ordering::Less};
468 if (greaterOrEqual) {
469 top = top.SubtractSigned(divisor).value;
470 }
471 auto doubled{top.AddUnsigned(top)};
472 top = doubled.value;
473 msb = doubled.carry;
474 return greaterOrEqual;
475 }
476
477 // Normalizes and marshals the fields of a floating-point number in place.
478 // The value is a number, and a zero fraction means a zero value (i.e.,
479 // a maximal exponent and zero fraction doesn't signify infinity, although
480 // this member function will detect overflow and encode infinities).
481 RealFlags Normalize(bool negative, int exponent, const Fraction &fraction,
482 Rounding rounding = TargetCharacteristics::defaultRounding,
483 RoundingBits *roundingBits = nullptr);
484
485 // Rounds a result, if necessary, in place.
486 RealFlags Round(Rounding, const RoundingBits &, bool multiply = false);
487
488 static void NormalizeAndRound(ValueWithRealFlags<Real> &result,
489 bool isNegative, int exponent, const Fraction &, Rounding, RoundingBits,
490 bool multiply = false);
491
492 Word word_{}; // an Integer<>
493};
494
495extern template class Real<Integer<16>, 11>; // IEEE half format
496extern template class Real<Integer<16>, 8>; // the "other" half format
497extern template class Real<Integer<32>, 24>; // IEEE single
498extern template class Real<Integer<64>, 53>; // IEEE double
499extern template class Real<X87IntegerContainer, 64>; // 80387 extended precision
500extern template class Real<Integer<128>, 113>; // IEEE quad
501// N.B. No "double-double" support.
502} // namespace Fortran::evaluate::value
503#endif // FORTRAN_EVALUATE_REAL_H_
504

Warning: This file is not a C or C++ file. It does not have highlighting.

source code of flang/include/flang/Evaluate/real.h