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 | |
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.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 | |
445 | private: |
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 | |
495 | extern template class Real<Integer<16>, 11>; // IEEE half format |
496 | extern template class Real<Integer<16>, 8>; // the "other" half format |
497 | extern template class Real<Integer<32>, 24>; // IEEE single |
498 | extern template class Real<Integer<64>, 53>; // IEEE double |
499 | extern template class Real<X87IntegerContainer, 64>; // 80387 extended precision |
500 | extern 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.