1//===-- lib/Evaluate/real.cpp ---------------------------------------------===//
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#include "flang/Evaluate/real.h"
10#include "int-power.h"
11#include "flang/Common/idioms.h"
12#include "flang/Decimal/decimal.h"
13#include "flang/Parser/characters.h"
14#include "llvm/Support/raw_ostream.h"
15#include <limits>
16
17namespace Fortran::evaluate::value {
18
19template <typename W, int P> Relation Real<W, P>::Compare(const Real &y) const {
20 if (IsNotANumber() || y.IsNotANumber()) { // NaN vs x, x vs NaN
21 return Relation::Unordered;
22 } else if (IsInfinite()) {
23 if (y.IsInfinite()) {
24 if (IsNegative()) { // -Inf vs +/-Inf
25 return y.IsNegative() ? Relation::Equal : Relation::Less;
26 } else { // +Inf vs +/-Inf
27 return y.IsNegative() ? Relation::Greater : Relation::Equal;
28 }
29 } else { // +/-Inf vs finite
30 return IsNegative() ? Relation::Less : Relation::Greater;
31 }
32 } else if (y.IsInfinite()) { // finite vs +/-Inf
33 return y.IsNegative() ? Relation::Greater : Relation::Less;
34 } else { // two finite numbers
35 bool isNegative{IsNegative()};
36 if (isNegative != y.IsNegative()) {
37 if (word_.IOR(y.word_).IBCLR(bits - 1).IsZero()) {
38 return Relation::Equal; // +/-0.0 == -/+0.0
39 } else {
40 return isNegative ? Relation::Less : Relation::Greater;
41 }
42 } else {
43 // same sign
44 Ordering order{evaluate::Compare(Exponent(), y.Exponent())};
45 if (order == Ordering::Equal) {
46 order = GetSignificand().CompareUnsigned(y.GetSignificand());
47 }
48 if (isNegative) {
49 order = Reverse(order);
50 }
51 return RelationFromOrdering(order);
52 }
53 }
54}
55
56template <typename W, int P>
57ValueWithRealFlags<Real<W, P>> Real<W, P>::Add(
58 const Real &y, Rounding rounding) const {
59 ValueWithRealFlags<Real> result;
60 if (IsNotANumber() || y.IsNotANumber()) {
61 result.value = NotANumber(); // NaN + x -> NaN
62 if (IsSignalingNaN() || y.IsSignalingNaN()) {
63 result.flags.set(RealFlag::InvalidArgument);
64 }
65 return result;
66 }
67 bool isNegative{IsNegative()};
68 bool yIsNegative{y.IsNegative()};
69 if (IsInfinite()) {
70 if (y.IsInfinite()) {
71 if (isNegative == yIsNegative) {
72 result.value = *this; // +/-Inf + +/-Inf -> +/-Inf
73 } else {
74 result.value = NotANumber(); // +/-Inf + -/+Inf -> NaN
75 result.flags.set(RealFlag::InvalidArgument);
76 }
77 } else {
78 result.value = *this; // +/-Inf + x -> +/-Inf
79 }
80 return result;
81 }
82 if (y.IsInfinite()) {
83 result.value = y; // x + +/-Inf -> +/-Inf
84 return result;
85 }
86 int exponent{Exponent()};
87 int yExponent{y.Exponent()};
88 if (exponent < yExponent) {
89 // y is larger in magnitude; simplify by reversing operands
90 return y.Add(*this, rounding);
91 }
92 if (exponent == yExponent && isNegative != yIsNegative) {
93 Ordering order{GetSignificand().CompareUnsigned(y.GetSignificand())};
94 if (order == Ordering::Less) {
95 // Same exponent, opposite signs, and y is larger in magnitude
96 return y.Add(*this, rounding);
97 }
98 if (order == Ordering::Equal) {
99 // x + (-x) -> +0.0 unless rounding is directed downwards
100 if (rounding.mode == common::RoundingMode::Down) {
101 result.value = NegativeZero();
102 }
103 return result;
104 }
105 }
106 // Our exponent is greater than y's, or the exponents match and y is not
107 // of the opposite sign and greater magnitude. So (x+y) will have the
108 // same sign as x.
109 Fraction fraction{GetFraction()};
110 Fraction yFraction{y.GetFraction()};
111 int rshift = exponent - yExponent;
112 if (exponent > 0 && yExponent == 0) {
113 --rshift; // correct overshift when only y is subnormal
114 }
115 RoundingBits roundingBits{yFraction, rshift};
116 yFraction = yFraction.SHIFTR(rshift);
117 bool carry{false};
118 if (isNegative != yIsNegative) {
119 // Opposite signs: subtract via addition of two's complement of y and
120 // the rounding bits.
121 yFraction = yFraction.NOT();
122 carry = roundingBits.Negate();
123 }
124 auto sum{fraction.AddUnsigned(yFraction, carry)};
125 fraction = sum.value;
126 if (isNegative == yIsNegative && sum.carry) {
127 roundingBits.ShiftRight(sum.value.BTEST(0));
128 fraction = fraction.SHIFTR(1).IBSET(fraction.bits - 1);
129 ++exponent;
130 }
131 NormalizeAndRound(
132 result, isNegative, exponent, fraction, rounding, roundingBits);
133 return result;
134}
135
136template <typename W, int P>
137ValueWithRealFlags<Real<W, P>> Real<W, P>::Multiply(
138 const Real &y, Rounding rounding) const {
139 ValueWithRealFlags<Real> result;
140 if (IsNotANumber() || y.IsNotANumber()) {
141 result.value = NotANumber(); // NaN * x -> NaN
142 if (IsSignalingNaN() || y.IsSignalingNaN()) {
143 result.flags.set(RealFlag::InvalidArgument);
144 }
145 } else {
146 bool isNegative{IsNegative() != y.IsNegative()};
147 if (IsInfinite() || y.IsInfinite()) {
148 if (IsZero() || y.IsZero()) {
149 result.value = NotANumber(); // 0 * Inf -> NaN
150 result.flags.set(RealFlag::InvalidArgument);
151 } else {
152 result.value = Infinity(isNegative);
153 }
154 } else {
155 auto product{GetFraction().MultiplyUnsigned(y.GetFraction())};
156 std::int64_t exponent{CombineExponents(y, false)};
157 if (exponent < 1) {
158 int rshift = 1 - exponent;
159 exponent = 1;
160 bool sticky{false};
161 if (rshift >= product.upper.bits + product.lower.bits) {
162 sticky = !product.lower.IsZero() || !product.upper.IsZero();
163 } else if (rshift >= product.lower.bits) {
164 sticky = !product.lower.IsZero() ||
165 !product.upper
166 .IAND(product.upper.MASKR(rshift - product.lower.bits))
167 .IsZero();
168 } else {
169 sticky = !product.lower.IAND(product.lower.MASKR(rshift)).IsZero();
170 }
171 product.lower = product.lower.SHIFTRWithFill(product.upper, rshift);
172 product.upper = product.upper.SHIFTR(rshift);
173 if (sticky) {
174 product.lower = product.lower.IBSET(0);
175 }
176 }
177 int leadz{product.upper.LEADZ()};
178 if (leadz >= product.upper.bits) {
179 leadz += product.lower.LEADZ();
180 }
181 int lshift{leadz};
182 if (lshift > exponent - 1) {
183 lshift = exponent - 1;
184 }
185 exponent -= lshift;
186 product.upper = product.upper.SHIFTLWithFill(product.lower, lshift);
187 product.lower = product.lower.SHIFTL(lshift);
188 RoundingBits roundingBits{product.lower, product.lower.bits};
189 NormalizeAndRound(result, isNegative, exponent, product.upper, rounding,
190 roundingBits, true /*multiply*/);
191 }
192 }
193 return result;
194}
195
196template <typename W, int P>
197ValueWithRealFlags<Real<W, P>> Real<W, P>::Divide(
198 const Real &y, Rounding rounding) const {
199 ValueWithRealFlags<Real> result;
200 if (IsNotANumber() || y.IsNotANumber()) {
201 result.value = NotANumber(); // NaN / x -> NaN, x / NaN -> NaN
202 if (IsSignalingNaN() || y.IsSignalingNaN()) {
203 result.flags.set(RealFlag::InvalidArgument);
204 }
205 } else {
206 bool isNegative{IsNegative() != y.IsNegative()};
207 if (IsInfinite()) {
208 if (y.IsInfinite()) {
209 result.value = NotANumber(); // Inf/Inf -> NaN
210 result.flags.set(RealFlag::InvalidArgument);
211 } else { // Inf/x -> Inf, Inf/0 -> Inf
212 result.value = Infinity(isNegative);
213 }
214 } else if (y.IsZero()) {
215 if (IsZero()) { // 0/0 -> NaN
216 result.value = NotANumber();
217 result.flags.set(RealFlag::InvalidArgument);
218 } else { // x/0 -> Inf, Inf/0 -> Inf
219 result.value = Infinity(isNegative);
220 result.flags.set(RealFlag::DivideByZero);
221 }
222 } else if (IsZero() || y.IsInfinite()) { // 0/x, x/Inf -> 0
223 if (isNegative) {
224 result.value = NegativeZero();
225 }
226 } else {
227 // dividend and divisor are both finite and nonzero numbers
228 Fraction top{GetFraction()}, divisor{y.GetFraction()};
229 std::int64_t exponent{CombineExponents(y, true)};
230 Fraction quotient;
231 bool msb{false};
232 if (!top.BTEST(top.bits - 1) || !divisor.BTEST(divisor.bits - 1)) {
233 // One or two subnormals
234 int topLshift{top.LEADZ()};
235 top = top.SHIFTL(topLshift);
236 int divisorLshift{divisor.LEADZ()};
237 divisor = divisor.SHIFTL(divisorLshift);
238 exponent += divisorLshift - topLshift;
239 }
240 for (int j{1}; j <= quotient.bits; ++j) {
241 if (NextQuotientBit(top, msb, divisor)) {
242 quotient = quotient.IBSET(quotient.bits - j);
243 }
244 }
245 bool guard{NextQuotientBit(top, msb, divisor)};
246 bool round{NextQuotientBit(top, msb, divisor)};
247 bool sticky{msb || !top.IsZero()};
248 RoundingBits roundingBits{guard, round, sticky};
249 if (exponent < 1) {
250 std::int64_t rshift{1 - exponent};
251 for (; rshift > 0; --rshift) {
252 roundingBits.ShiftRight(quotient.BTEST(0));
253 quotient = quotient.SHIFTR(1);
254 }
255 exponent = 1;
256 }
257 NormalizeAndRound(
258 result, isNegative, exponent, quotient, rounding, roundingBits);
259 }
260 }
261 return result;
262}
263
264template <typename W, int P>
265ValueWithRealFlags<Real<W, P>> Real<W, P>::SQRT(Rounding rounding) const {
266 ValueWithRealFlags<Real> result;
267 if (IsNotANumber()) {
268 result.value = NotANumber();
269 if (IsSignalingNaN()) {
270 result.flags.set(RealFlag::InvalidArgument);
271 }
272 } else if (IsNegative()) {
273 if (IsZero()) {
274 // SQRT(-0) == -0 in IEEE-754.
275 result.value = NegativeZero();
276 } else {
277 result.flags.set(RealFlag::InvalidArgument);
278 result.value = NotANumber();
279 }
280 } else if (IsInfinite()) {
281 // SQRT(+Inf) == +Inf
282 result.value = Infinity(false);
283 } else if (IsZero()) {
284 result.value = PositiveZero();
285 } else {
286 int expo{UnbiasedExponent()};
287 if (expo < -1 || expo > 1) {
288 // Reduce the range to [0.5 .. 4.0) by dividing by an integral power
289 // of four to avoid trouble with very large and very small values
290 // (esp. truncation of subnormals).
291 // SQRT(2**(2a) * x) = SQRT(2**(2a)) * SQRT(x) = 2**a * SQRT(x)
292 Real scaled;
293 int adjust{expo / 2};
294 scaled.Normalize(false, expo - 2 * adjust + exponentBias, GetFraction());
295 result = scaled.SQRT(rounding);
296 result.value.Normalize(false,
297 result.value.UnbiasedExponent() + adjust + exponentBias,
298 result.value.GetFraction());
299 return result;
300 }
301 // (-1) <= expo <= 1; use it as a shift to set the desired square.
302 using Extended = typename value::Integer<(binaryPrecision + 2)>;
303 Extended goal{
304 Extended::ConvertUnsigned(GetFraction()).value.SHIFTL(expo + 1)};
305 // Calculate the exact square root by maximizing a value whose square
306 // does not exceed the goal. Use two extra bits of precision for
307 // rounding.
308 bool sticky{true};
309 Extended extFrac{};
310 for (int bit{Extended::bits - 1}; bit >= 0; --bit) {
311 Extended next{extFrac.IBSET(bit)};
312 auto squared{next.MultiplyUnsigned(next)};
313 auto cmp{squared.upper.CompareUnsigned(goal)};
314 if (cmp == Ordering::Less) {
315 extFrac = next;
316 } else if (cmp == Ordering::Equal && squared.lower.IsZero()) {
317 extFrac = next;
318 sticky = false;
319 break; // exact result
320 }
321 }
322 RoundingBits roundingBits{extFrac.BTEST(1), extFrac.BTEST(0), sticky};
323 NormalizeAndRound(result, false, exponentBias,
324 Fraction::ConvertUnsigned(extFrac.SHIFTR(2)).value, rounding,
325 roundingBits);
326 }
327 return result;
328}
329
330template <typename W, int P>
331ValueWithRealFlags<Real<W, P>> Real<W, P>::NEAREST(bool upward) const {
332 ValueWithRealFlags<Real> result;
333 bool isNegative{IsNegative()};
334 if (IsFinite()) {
335 Fraction fraction{GetFraction()};
336 int expo{Exponent()};
337 Fraction one{1};
338 Fraction nearest;
339 if (upward != isNegative) { // upward in magnitude
340 auto next{fraction.AddUnsigned(one)};
341 if (next.carry) {
342 ++expo;
343 nearest = Fraction::Least(); // MSB only
344 } else {
345 nearest = next.value;
346 }
347 } else { // downward in magnitude
348 if (IsZero()) {
349 nearest = 1; // smallest magnitude negative subnormal
350 isNegative = !isNegative;
351 } else {
352 auto sub1{fraction.SubtractSigned(one)};
353 if (sub1.overflow && expo > 1) {
354 nearest = Fraction{0}.NOT();
355 --expo;
356 } else {
357 nearest = sub1.value;
358 }
359 }
360 }
361 result.value.Normalize(isNegative, expo, nearest);
362 } else if (IsInfinite()) {
363 if (upward == isNegative) {
364 result.value =
365 isNegative ? HUGE().Negate() : HUGE(); // largest mag finite
366 } else {
367 result.value = *this;
368 }
369 } else { // NaN
370 result.flags.set(RealFlag::InvalidArgument);
371 result.value = *this;
372 }
373 return result;
374}
375
376// HYPOT(x,y) = SQRT(x**2 + y**2) by definition, but those squared intermediate
377// values are susceptible to over/underflow when computed naively.
378// Assuming that x>=y, calculate instead:
379// HYPOT(x,y) = SQRT(x**2 * (1+(y/x)**2))
380// = ABS(x) * SQRT(1+(y/x)**2)
381template <typename W, int P>
382ValueWithRealFlags<Real<W, P>> Real<W, P>::HYPOT(
383 const Real &y, Rounding rounding) const {
384 ValueWithRealFlags<Real> result;
385 if (IsNotANumber() || y.IsNotANumber()) {
386 result.flags.set(RealFlag::InvalidArgument);
387 result.value = NotANumber();
388 } else if (ABS().Compare(y.ABS()) == Relation::Less) {
389 return y.HYPOT(*this);
390 } else if (IsZero()) {
391 return result; // x==y==0
392 } else {
393 auto yOverX{y.Divide(*this, rounding)}; // y/x
394 bool inexact{yOverX.flags.test(RealFlag::Inexact)};
395 auto squared{yOverX.value.Multiply(yOverX.value, rounding)}; // (y/x)**2
396 inexact |= squared.flags.test(RealFlag::Inexact);
397 Real one;
398 one.Normalize(false, exponentBias, Fraction::MASKL(1)); // 1.0
399 auto sum{squared.value.Add(one, rounding)}; // 1.0 + (y/x)**2
400 inexact |= sum.flags.test(RealFlag::Inexact);
401 auto sqrt{sum.value.SQRT()};
402 inexact |= sqrt.flags.test(RealFlag::Inexact);
403 result = sqrt.value.Multiply(ABS(), rounding);
404 if (inexact) {
405 result.flags.set(RealFlag::Inexact);
406 }
407 }
408 return result;
409}
410
411// MOD(x,y) = x - AINT(x/y)*y in the standard; unfortunately, this definition
412// can be pretty inaccurate when x is much larger than y in magnitude due to
413// cancellation. Implement instead with (essentially) arbitrary precision
414// long division, discarding the quotient and returning the remainder.
415// See runtime/numeric.cpp for more details.
416template <typename W, int P>
417ValueWithRealFlags<Real<W, P>> Real<W, P>::MOD(
418 const Real &p, Rounding rounding) const {
419 ValueWithRealFlags<Real> result;
420 if (IsNotANumber() || p.IsNotANumber() || IsInfinite()) {
421 result.flags.set(RealFlag::InvalidArgument);
422 result.value = NotANumber();
423 } else if (p.IsZero()) {
424 result.flags.set(RealFlag::DivideByZero);
425 result.value = NotANumber();
426 } else if (p.IsInfinite()) {
427 result.value = *this;
428 } else {
429 result.value = ABS();
430 auto pAbs{p.ABS()};
431 Real half, adj;
432 half.Normalize(false, exponentBias - 1, Fraction::MASKL(1)); // 0.5
433 for (adj.Normalize(false, Exponent(), pAbs.GetFraction());
434 result.value.Compare(pAbs) != Relation::Less;
435 adj = adj.Multiply(half).value) {
436 if (result.value.Compare(adj) != Relation::Less) {
437 result.value =
438 result.value.Subtract(adj, rounding).AccumulateFlags(result.flags);
439 if (result.value.IsZero()) {
440 break;
441 }
442 }
443 }
444 if (IsNegative()) {
445 result.value = result.value.Negate();
446 }
447 }
448 return result;
449}
450
451// MODULO(x,y) = x - FLOOR(x/y)*y in the standard; here, it is defined
452// in terms of MOD() with adjustment of the result.
453template <typename W, int P>
454ValueWithRealFlags<Real<W, P>> Real<W, P>::MODULO(
455 const Real &p, Rounding rounding) const {
456 ValueWithRealFlags<Real> result{MOD(p, rounding)};
457 if (IsNegative() != p.IsNegative()) {
458 if (result.value.IsZero()) {
459 result.value = result.value.Negate();
460 } else {
461 result.value =
462 result.value.Add(p, rounding).AccumulateFlags(result.flags);
463 }
464 }
465 return result;
466}
467
468template <typename W, int P>
469ValueWithRealFlags<Real<W, P>> Real<W, P>::DIM(
470 const Real &y, Rounding rounding) const {
471 ValueWithRealFlags<Real> result;
472 if (IsNotANumber() || y.IsNotANumber()) {
473 result.flags.set(RealFlag::InvalidArgument);
474 result.value = NotANumber();
475 } else if (Compare(y) == Relation::Greater) {
476 result = Subtract(y, rounding);
477 } else {
478 // result is already zero
479 }
480 return result;
481}
482
483template <typename W, int P>
484ValueWithRealFlags<Real<W, P>> Real<W, P>::ToWholeNumber(
485 common::RoundingMode mode) const {
486 ValueWithRealFlags<Real> result{*this};
487 if (IsNotANumber()) {
488 result.flags.set(RealFlag::InvalidArgument);
489 result.value = NotANumber();
490 } else if (IsInfinite()) {
491 result.flags.set(RealFlag::Overflow);
492 } else {
493 constexpr int noClipExponent{exponentBias + binaryPrecision - 1};
494 if (Exponent() < noClipExponent) {
495 Real adjust; // ABS(EPSILON(adjust)) == 0.5
496 adjust.Normalize(IsSignBitSet(), noClipExponent, Fraction::MASKL(1));
497 // Compute ival=(*this + adjust), losing any fractional bits; keep flags
498 result = Add(adjust, Rounding{mode});
499 result.flags.reset(RealFlag::Inexact); // result *is* exact
500 // Return (ival-adjust) with original sign in case we've generated a zero.
501 result.value =
502 result.value.Subtract(adjust, Rounding{common::RoundingMode::ToZero})
503 .value.SIGN(*this);
504 }
505 }
506 return result;
507}
508
509template <typename W, int P>
510RealFlags Real<W, P>::Normalize(bool negative, int exponent,
511 const Fraction &fraction, Rounding rounding, RoundingBits *roundingBits) {
512 int lshift{fraction.LEADZ()};
513 if (lshift == fraction.bits /* fraction is zero */ &&
514 (!roundingBits || roundingBits->empty())) {
515 // No fraction, no rounding bits -> +/-0.0
516 exponent = lshift = 0;
517 } else if (lshift < exponent) {
518 exponent -= lshift;
519 } else if (exponent > 0) {
520 lshift = exponent - 1;
521 exponent = 0;
522 } else if (lshift == 0) {
523 exponent = 1;
524 } else {
525 lshift = 0;
526 }
527 if (exponent >= maxExponent) {
528 // Infinity or overflow
529 if (rounding.mode == common::RoundingMode::TiesToEven ||
530 rounding.mode == common::RoundingMode::TiesAwayFromZero ||
531 (rounding.mode == common::RoundingMode::Up && !negative) ||
532 (rounding.mode == common::RoundingMode::Down && negative)) {
533 word_ = Word{maxExponent}.SHIFTL(significandBits); // Inf
534 if constexpr (!isImplicitMSB) {
535 word_ = word_.IBSET(significandBits - 1);
536 }
537 } else {
538 // directed rounding: round to largest finite value rather than infinity
539 // (x86 does this, not sure whether it's standard behavior)
540 word_ = Word{word_.MASKR(word_.bits - 1)};
541 if constexpr (isImplicitMSB) {
542 word_ = word_.IBCLR(significandBits);
543 }
544 }
545 if (negative) {
546 word_ = word_.IBSET(bits - 1);
547 }
548 RealFlags flags{RealFlag::Overflow};
549 if (!fraction.IsZero()) {
550 flags.set(RealFlag::Inexact);
551 }
552 return flags;
553 }
554 word_ = Word::ConvertUnsigned(fraction).value;
555 if (lshift > 0) {
556 word_ = word_.SHIFTL(lshift);
557 if (roundingBits) {
558 for (; lshift > 0; --lshift) {
559 if (roundingBits->ShiftLeft()) {
560 word_ = word_.IBSET(lshift - 1);
561 }
562 }
563 }
564 }
565 if constexpr (isImplicitMSB) {
566 word_ = word_.IBCLR(significandBits);
567 }
568 word_ = word_.IOR(Word{exponent}.SHIFTL(significandBits));
569 if (negative) {
570 word_ = word_.IBSET(bits - 1);
571 }
572 return {};
573}
574
575template <typename W, int P>
576RealFlags Real<W, P>::Round(
577 Rounding rounding, const RoundingBits &bits, bool multiply) {
578 int origExponent{Exponent()};
579 RealFlags flags;
580 bool inexact{!bits.empty()};
581 if (inexact) {
582 flags.set(RealFlag::Inexact);
583 }
584 if (origExponent < maxExponent &&
585 bits.MustRound(rounding, IsNegative(), word_.BTEST(0) /* is odd */)) {
586 typename Fraction::ValueWithCarry sum{
587 GetFraction().AddUnsigned(Fraction{}, true)};
588 int newExponent{origExponent};
589 if (sum.carry) {
590 // The fraction was all ones before rounding; sum.value is now zero
591 sum.value = sum.value.IBSET(binaryPrecision - 1);
592 if (++newExponent >= maxExponent) {
593 flags.set(RealFlag::Overflow); // rounded away to an infinity
594 }
595 }
596 flags |= Normalize(IsNegative(), newExponent, sum.value);
597 }
598 if (inexact && origExponent == 0) {
599 // inexact subnormal input: signal Underflow unless in an x86-specific
600 // edge case
601 if (rounding.x86CompatibleBehavior && Exponent() != 0 && multiply &&
602 bits.sticky() &&
603 (bits.guard() ||
604 (rounding.mode != common::RoundingMode::Up &&
605 rounding.mode != common::RoundingMode::Down))) {
606 // x86 edge case in which Underflow fails to signal when a subnormal
607 // inexact multiplication product rounds to a normal result when
608 // the guard bit is set or we're not using directed rounding
609 } else {
610 flags.set(RealFlag::Underflow);
611 }
612 }
613 return flags;
614}
615
616template <typename W, int P>
617void Real<W, P>::NormalizeAndRound(ValueWithRealFlags<Real> &result,
618 bool isNegative, int exponent, const Fraction &fraction, Rounding rounding,
619 RoundingBits roundingBits, bool multiply) {
620 result.flags |= result.value.Normalize(
621 isNegative, exponent, fraction, rounding, &roundingBits);
622 result.flags |= result.value.Round(rounding, roundingBits, multiply);
623}
624
625inline enum decimal::FortranRounding MapRoundingMode(
626 common::RoundingMode rounding) {
627 switch (rounding) {
628 case common::RoundingMode::TiesToEven:
629 break;
630 case common::RoundingMode::ToZero:
631 return decimal::RoundToZero;
632 case common::RoundingMode::Down:
633 return decimal::RoundDown;
634 case common::RoundingMode::Up:
635 return decimal::RoundUp;
636 case common::RoundingMode::TiesAwayFromZero:
637 return decimal::RoundCompatible;
638 }
639 return decimal::RoundNearest; // dodge gcc warning about lack of result
640}
641
642inline RealFlags MapFlags(decimal::ConversionResultFlags flags) {
643 RealFlags result;
644 if (flags & decimal::Overflow) {
645 result.set(RealFlag::Overflow);
646 }
647 if (flags & decimal::Inexact) {
648 result.set(RealFlag::Inexact);
649 }
650 if (flags & decimal::Invalid) {
651 result.set(RealFlag::InvalidArgument);
652 }
653 return result;
654}
655
656template <typename W, int P>
657ValueWithRealFlags<Real<W, P>> Real<W, P>::Read(
658 const char *&p, Rounding rounding) {
659 auto converted{
660 decimal::ConvertToBinary<P>(p, MapRoundingMode(rounding.mode))};
661 const auto *value{reinterpret_cast<Real<W, P> *>(&converted.binary)};
662 return {*value, MapFlags(converted.flags)};
663}
664
665template <typename W, int P> std::string Real<W, P>::DumpHexadecimal() const {
666 if (IsNotANumber()) {
667 return "NaN0x"s + word_.Hexadecimal();
668 } else if (IsNegative()) {
669 return "-"s + Negate().DumpHexadecimal();
670 } else if (IsInfinite()) {
671 return "Inf"s;
672 } else if (IsZero()) {
673 return "0.0"s;
674 } else {
675 Fraction frac{GetFraction()};
676 std::string result{"0x"};
677 char intPart = '0' + frac.BTEST(frac.bits - 1);
678 result += intPart;
679 result += '.';
680 int trailz{frac.TRAILZ()};
681 if (trailz >= frac.bits - 1) {
682 result += '0';
683 } else {
684 int remainingBits{frac.bits - 1 - trailz};
685 int wholeNybbles{remainingBits / 4};
686 int lostBits{remainingBits - 4 * wholeNybbles};
687 if (wholeNybbles > 0) {
688 std::string fracHex{frac.SHIFTR(trailz + lostBits)
689 .IAND(frac.MASKR(4 * wholeNybbles))
690 .Hexadecimal()};
691 std::size_t field = wholeNybbles;
692 if (fracHex.size() < field) {
693 result += std::string(field - fracHex.size(), '0');
694 }
695 result += fracHex;
696 }
697 if (lostBits > 0) {
698 result += frac.SHIFTR(trailz)
699 .IAND(frac.MASKR(lostBits))
700 .SHIFTL(4 - lostBits)
701 .Hexadecimal();
702 }
703 }
704 result += 'p';
705 int exponent = Exponent() - exponentBias;
706 if (intPart == '0') {
707 exponent += 1;
708 }
709 result += Integer<32>{exponent}.SignedDecimal();
710 return result;
711 }
712}
713
714template <typename W, int P>
715llvm::raw_ostream &Real<W, P>::AsFortran(
716 llvm::raw_ostream &o, int kind, bool minimal) const {
717 if (IsNotANumber()) {
718 o << "(0._" << kind << "/0.)";
719 } else if (IsInfinite()) {
720 if (IsNegative()) {
721 o << "(-1._" << kind << "/0.)";
722 } else {
723 o << "(1._" << kind << "/0.)";
724 }
725 } else {
726 using B = decimal::BinaryFloatingPointNumber<P>;
727 B value{word_.template ToUInt<typename B::RawType>()};
728 char buffer[common::MaxDecimalConversionDigits(P) +
729 EXTRA_DECIMAL_CONVERSION_SPACE];
730 decimal::DecimalConversionFlags flags{}; // default: exact representation
731 if (minimal) {
732 flags = decimal::Minimize;
733 }
734 auto result{decimal::ConvertToDecimal<P>(buffer, sizeof buffer, flags,
735 static_cast<int>(sizeof buffer), decimal::RoundNearest, value)};
736 const char *p{result.str};
737 if (DEREF(p) == '-' || *p == '+') {
738 o << *p++;
739 }
740 int expo{result.decimalExponent};
741 if (*p != '0') {
742 --expo;
743 }
744 o << *p << '.' << (p + 1);
745 if (expo != 0) {
746 o << 'e' << expo;
747 }
748 o << '_' << kind;
749 }
750 return o;
751}
752
753// 16.9.180
754template <typename W, int P> Real<W, P> Real<W, P>::RRSPACING() const {
755 if (IsNotANumber()) {
756 return *this;
757 } else if (IsInfinite()) {
758 return NotANumber();
759 } else {
760 Real result;
761 result.Normalize(false, binaryPrecision + exponentBias - 1, GetFraction());
762 return result;
763 }
764}
765
766// 16.9.180
767template <typename W, int P> Real<W, P> Real<W, P>::SPACING() const {
768 if (IsNotANumber()) {
769 return *this;
770 } else if (IsInfinite()) {
771 return NotANumber();
772 } else if (IsZero() || IsSubnormal()) {
773 return TINY(); // standard & 100% portable
774 } else {
775 Real result;
776 result.Normalize(false, Exponent(), Fraction::MASKR(1));
777 // Can the result be less than TINY()? No, with five commonly
778 // used compilers; yes, with two less commonly used ones.
779 return result.IsZero() || result.IsSubnormal() ? TINY() : result;
780 }
781}
782
783// 16.9.171
784template <typename W, int P>
785Real<W, P> Real<W, P>::SET_EXPONENT(std::int64_t expo) const {
786 if (IsNotANumber()) {
787 return *this;
788 } else if (IsInfinite()) {
789 return NotANumber();
790 } else if (IsZero()) {
791 return *this;
792 } else {
793 return SCALE(Integer<64>(expo - UnbiasedExponent() - 1)).value;
794 }
795}
796
797// 16.9.171
798template <typename W, int P> Real<W, P> Real<W, P>::FRACTION() const {
799 return SET_EXPONENT(0);
800}
801
802template class Real<Integer<16>, 11>;
803template class Real<Integer<16>, 8>;
804template class Real<Integer<32>, 24>;
805template class Real<Integer<64>, 53>;
806template class Real<X87IntegerContainer, 64>;
807template class Real<Integer<128>, 113>;
808} // namespace Fortran::evaluate::value
809

source code of flang/lib/Evaluate/real.cpp