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 
17 namespace Fortran::evaluate::value {
18 
Compare(const Real & y) const19 template <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 
56 template <typename W, int P>
Add(const Real & y,Rounding rounding) const57 ValueWithRealFlags<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 
136 template <typename W, int P>
Multiply(const Real & y,Rounding rounding) const137 ValueWithRealFlags<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 
196 template <typename W, int P>
Divide(const Real & y,Rounding rounding) const197 ValueWithRealFlags<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 
264 template <typename W, int P>
SQRT(Rounding rounding) const265 ValueWithRealFlags<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 
330 template <typename W, int P>
NEAREST(bool upward) const331 ValueWithRealFlags<Real<W, P>> Real<W, P>::NEAREST(bool upward) const {
332   ValueWithRealFlags<Real> result;
333   if (IsFinite()) {
334     Fraction fraction{GetFraction()};
335     int expo{Exponent()};
336     Fraction one{1};
337     Fraction nearest;
338     bool isNegative{IsNegative()};
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) {
354           nearest = Fraction{0}.NOT();
355           --expo;
356         } else {
357           nearest = sub1.value;
358         }
359       }
360     }
361     result.flags = result.value.Normalize(isNegative, expo, nearest);
362   } else {
363     result.flags.set(RealFlag::InvalidArgument);
364     result.value = *this;
365   }
366   return result;
367 }
368 
369 // HYPOT(x,y) = SQRT(x**2 + y**2) by definition, but those squared intermediate
370 // values are susceptible to over/underflow when computed naively.
371 // Assuming that x>=y, calculate instead:
372 //   HYPOT(x,y) = SQRT(x**2 * (1+(y/x)**2))
373 //              = ABS(x) * SQRT(1+(y/x)**2)
374 template <typename W, int P>
HYPOT(const Real & y,Rounding rounding) const375 ValueWithRealFlags<Real<W, P>> Real<W, P>::HYPOT(
376     const Real &y, Rounding rounding) const {
377   ValueWithRealFlags<Real> result;
378   if (IsNotANumber() || y.IsNotANumber()) {
379     result.flags.set(RealFlag::InvalidArgument);
380     result.value = NotANumber();
381   } else if (ABS().Compare(y.ABS()) == Relation::Less) {
382     return y.HYPOT(*this);
383   } else if (IsZero()) {
384     return result; // x==y==0
385   } else {
386     auto yOverX{y.Divide(*this, rounding)}; // y/x
387     bool inexact{yOverX.flags.test(RealFlag::Inexact)};
388     auto squared{yOverX.value.Multiply(yOverX.value, rounding)}; // (y/x)**2
389     inexact |= squared.flags.test(RealFlag::Inexact);
390     Real one;
391     one.Normalize(false, exponentBias, Fraction::MASKL(1)); // 1.0
392     auto sum{squared.value.Add(one, rounding)}; // 1.0 + (y/x)**2
393     inexact |= sum.flags.test(RealFlag::Inexact);
394     auto sqrt{sum.value.SQRT()};
395     inexact |= sqrt.flags.test(RealFlag::Inexact);
396     result = sqrt.value.Multiply(ABS(), rounding);
397     if (inexact) {
398       result.flags.set(RealFlag::Inexact);
399     }
400   }
401   return result;
402 }
403 
404 // MOD(x,y) = x - AINT(x/y)*y
405 template <typename W, int P>
MOD(const Real & y,Rounding rounding) const406 ValueWithRealFlags<Real<W, P>> Real<W, P>::MOD(
407     const Real &y, Rounding rounding) const {
408   ValueWithRealFlags<Real> result;
409   auto quotient{Divide(y, rounding)};
410   if (quotient.value.IsInfinite() && IsFinite() && y.IsFinite() &&
411       !y.IsZero()) {
412     // x/y overflowed -- so it must be an integer in this representation and
413     // the result must be a zero.
414     if (IsNegative()) {
415       result.value = Real{}.Negate(); // -0.
416     }
417   } else {
418     Real toInt{quotient.AccumulateFlags(result.flags)
419                    .ToWholeNumber(common::RoundingMode::ToZero)
420                    .AccumulateFlags(result.flags)};
421     Real product{toInt.Multiply(y, rounding).AccumulateFlags(result.flags)};
422     result.value = Subtract(product, rounding).AccumulateFlags(result.flags);
423   }
424   return result;
425 }
426 
427 // MODULO(x,y) = x - FLOOR(x/y)*y
428 template <typename W, int P>
MODULO(const Real & y,Rounding rounding) const429 ValueWithRealFlags<Real<W, P>> Real<W, P>::MODULO(
430     const Real &y, Rounding rounding) const {
431   ValueWithRealFlags<Real> result;
432   auto quotient{Divide(y, rounding)};
433   if (quotient.value.IsInfinite() && IsFinite() && y.IsFinite() &&
434       !y.IsZero()) {
435     // x/y overflowed -- so it must be an integer in this representation and
436     // the result must be a zero.
437     if (y.IsNegative()) {
438       result.value = Real{}.Negate(); // -0.
439     }
440   } else {
441     Real toInt{quotient.AccumulateFlags(result.flags)
442                    .ToWholeNumber(common::RoundingMode::Down)
443                    .AccumulateFlags(result.flags)};
444     Real product{toInt.Multiply(y, rounding).AccumulateFlags(result.flags)};
445     result.value = Subtract(product, rounding).AccumulateFlags(result.flags);
446   }
447   return result;
448 }
449 
450 template <typename W, int P>
DIM(const Real & y,Rounding rounding) const451 ValueWithRealFlags<Real<W, P>> Real<W, P>::DIM(
452     const Real &y, Rounding rounding) const {
453   ValueWithRealFlags<Real> result;
454   if (IsNotANumber() || y.IsNotANumber()) {
455     result.flags.set(RealFlag::InvalidArgument);
456     result.value = NotANumber();
457   } else if (Compare(y) == Relation::Greater) {
458     result = Subtract(y, rounding);
459   } else {
460     // result is already zero
461   }
462   return result;
463 }
464 
465 template <typename W, int P>
ToWholeNumber(common::RoundingMode mode) const466 ValueWithRealFlags<Real<W, P>> Real<W, P>::ToWholeNumber(
467     common::RoundingMode mode) const {
468   ValueWithRealFlags<Real> result{*this};
469   if (IsNotANumber()) {
470     result.flags.set(RealFlag::InvalidArgument);
471     result.value = NotANumber();
472   } else if (IsInfinite()) {
473     result.flags.set(RealFlag::Overflow);
474   } else {
475     constexpr int noClipExponent{exponentBias + binaryPrecision - 1};
476     if (Exponent() < noClipExponent) {
477       Real adjust; // ABS(EPSILON(adjust)) == 0.5
478       adjust.Normalize(IsSignBitSet(), noClipExponent, Fraction::MASKL(1));
479       // Compute ival=(*this + adjust), losing any fractional bits; keep flags
480       result = Add(adjust, Rounding{mode});
481       result.flags.reset(RealFlag::Inexact); // result *is* exact
482       // Return (ival-adjust) with original sign in case we've generated a zero.
483       result.value =
484           result.value.Subtract(adjust, Rounding{common::RoundingMode::ToZero})
485               .value.SIGN(*this);
486     }
487   }
488   return result;
489 }
490 
491 template <typename W, int P>
Normalize(bool negative,int exponent,const Fraction & fraction,Rounding rounding,RoundingBits * roundingBits)492 RealFlags Real<W, P>::Normalize(bool negative, int exponent,
493     const Fraction &fraction, Rounding rounding, RoundingBits *roundingBits) {
494   int lshift{fraction.LEADZ()};
495   if (lshift == fraction.bits /* fraction is zero */ &&
496       (!roundingBits || roundingBits->empty())) {
497     // No fraction, no rounding bits -> +/-0.0
498     exponent = lshift = 0;
499   } else if (lshift < exponent) {
500     exponent -= lshift;
501   } else if (exponent > 0) {
502     lshift = exponent - 1;
503     exponent = 0;
504   } else if (lshift == 0) {
505     exponent = 1;
506   } else {
507     lshift = 0;
508   }
509   if (exponent >= maxExponent) {
510     // Infinity or overflow
511     if (rounding.mode == common::RoundingMode::TiesToEven ||
512         rounding.mode == common::RoundingMode::TiesAwayFromZero ||
513         (rounding.mode == common::RoundingMode::Up && !negative) ||
514         (rounding.mode == common::RoundingMode::Down && negative)) {
515       word_ = Word{maxExponent}.SHIFTL(significandBits); // Inf
516     } else {
517       // directed rounding: round to largest finite value rather than infinity
518       // (x86 does this, not sure whether it's standard behavior)
519       word_ = Word{word_.MASKR(word_.bits - 1)}.IBCLR(significandBits);
520     }
521     if (negative) {
522       word_ = word_.IBSET(bits - 1);
523     }
524     RealFlags flags{RealFlag::Overflow};
525     if (!fraction.IsZero()) {
526       flags.set(RealFlag::Inexact);
527     }
528     return flags;
529   }
530   word_ = Word::ConvertUnsigned(fraction).value;
531   if (lshift > 0) {
532     word_ = word_.SHIFTL(lshift);
533     if (roundingBits) {
534       for (; lshift > 0; --lshift) {
535         if (roundingBits->ShiftLeft()) {
536           word_ = word_.IBSET(lshift - 1);
537         }
538       }
539     }
540   }
541   if constexpr (isImplicitMSB) {
542     word_ = word_.IBCLR(significandBits);
543   }
544   word_ = word_.IOR(Word{exponent}.SHIFTL(significandBits));
545   if (negative) {
546     word_ = word_.IBSET(bits - 1);
547   }
548   return {};
549 }
550 
551 template <typename W, int P>
Round(Rounding rounding,const RoundingBits & bits,bool multiply)552 RealFlags Real<W, P>::Round(
553     Rounding rounding, const RoundingBits &bits, bool multiply) {
554   int origExponent{Exponent()};
555   RealFlags flags;
556   bool inexact{!bits.empty()};
557   if (inexact) {
558     flags.set(RealFlag::Inexact);
559   }
560   if (origExponent < maxExponent &&
561       bits.MustRound(rounding, IsNegative(), word_.BTEST(0) /* is odd */)) {
562     typename Fraction::ValueWithCarry sum{
563         GetFraction().AddUnsigned(Fraction{}, true)};
564     int newExponent{origExponent};
565     if (sum.carry) {
566       // The fraction was all ones before rounding; sum.value is now zero
567       sum.value = sum.value.IBSET(binaryPrecision - 1);
568       if (++newExponent >= maxExponent) {
569         flags.set(RealFlag::Overflow); // rounded away to an infinity
570       }
571     }
572     flags |= Normalize(IsNegative(), newExponent, sum.value);
573   }
574   if (inexact && origExponent == 0) {
575     // inexact subnormal input: signal Underflow unless in an x86-specific
576     // edge case
577     if (rounding.x86CompatibleBehavior && Exponent() != 0 && multiply &&
578         bits.sticky() &&
579         (bits.guard() ||
580             (rounding.mode != common::RoundingMode::Up &&
581                 rounding.mode != common::RoundingMode::Down))) {
582       // x86 edge case in which Underflow fails to signal when a subnormal
583       // inexact multiplication product rounds to a normal result when
584       // the guard bit is set or we're not using directed rounding
585     } else {
586       flags.set(RealFlag::Underflow);
587     }
588   }
589   return flags;
590 }
591 
592 template <typename W, int P>
NormalizeAndRound(ValueWithRealFlags<Real> & result,bool isNegative,int exponent,const Fraction & fraction,Rounding rounding,RoundingBits roundingBits,bool multiply)593 void Real<W, P>::NormalizeAndRound(ValueWithRealFlags<Real> &result,
594     bool isNegative, int exponent, const Fraction &fraction, Rounding rounding,
595     RoundingBits roundingBits, bool multiply) {
596   result.flags |= result.value.Normalize(
597       isNegative, exponent, fraction, rounding, &roundingBits);
598   result.flags |= result.value.Round(rounding, roundingBits, multiply);
599 }
600 
MapRoundingMode(common::RoundingMode rounding)601 inline enum decimal::FortranRounding MapRoundingMode(
602     common::RoundingMode rounding) {
603   switch (rounding) {
604   case common::RoundingMode::TiesToEven:
605     break;
606   case common::RoundingMode::ToZero:
607     return decimal::RoundToZero;
608   case common::RoundingMode::Down:
609     return decimal::RoundDown;
610   case common::RoundingMode::Up:
611     return decimal::RoundUp;
612   case common::RoundingMode::TiesAwayFromZero:
613     return decimal::RoundCompatible;
614   }
615   return decimal::RoundNearest; // dodge gcc warning about lack of result
616 }
617 
MapFlags(decimal::ConversionResultFlags flags)618 inline RealFlags MapFlags(decimal::ConversionResultFlags flags) {
619   RealFlags result;
620   if (flags & decimal::Overflow) {
621     result.set(RealFlag::Overflow);
622   }
623   if (flags & decimal::Inexact) {
624     result.set(RealFlag::Inexact);
625   }
626   if (flags & decimal::Invalid) {
627     result.set(RealFlag::InvalidArgument);
628   }
629   return result;
630 }
631 
632 template <typename W, int P>
Read(const char * & p,Rounding rounding)633 ValueWithRealFlags<Real<W, P>> Real<W, P>::Read(
634     const char *&p, Rounding rounding) {
635   auto converted{
636       decimal::ConvertToBinary<P>(p, MapRoundingMode(rounding.mode))};
637   const auto *value{reinterpret_cast<Real<W, P> *>(&converted.binary)};
638   return {*value, MapFlags(converted.flags)};
639 }
640 
DumpHexadecimal() const641 template <typename W, int P> std::string Real<W, P>::DumpHexadecimal() const {
642   if (IsNotANumber()) {
643     return "NaN0x"s + word_.Hexadecimal();
644   } else if (IsNegative()) {
645     return "-"s + Negate().DumpHexadecimal();
646   } else if (IsInfinite()) {
647     return "Inf"s;
648   } else if (IsZero()) {
649     return "0.0"s;
650   } else {
651     Fraction frac{GetFraction()};
652     std::string result{"0x"};
653     char intPart = '0' + frac.BTEST(frac.bits - 1);
654     result += intPart;
655     result += '.';
656     int trailz{frac.TRAILZ()};
657     if (trailz >= frac.bits - 1) {
658       result += '0';
659     } else {
660       int remainingBits{frac.bits - 1 - trailz};
661       int wholeNybbles{remainingBits / 4};
662       int lostBits{remainingBits - 4 * wholeNybbles};
663       if (wholeNybbles > 0) {
664         std::string fracHex{frac.SHIFTR(trailz + lostBits)
665                                 .IAND(frac.MASKR(4 * wholeNybbles))
666                                 .Hexadecimal()};
667         std::size_t field = wholeNybbles;
668         if (fracHex.size() < field) {
669           result += std::string(field - fracHex.size(), '0');
670         }
671         result += fracHex;
672       }
673       if (lostBits > 0) {
674         result += frac.SHIFTR(trailz)
675                       .IAND(frac.MASKR(lostBits))
676                       .SHIFTL(4 - lostBits)
677                       .Hexadecimal();
678       }
679     }
680     result += 'p';
681     int exponent = Exponent() - exponentBias;
682     if (intPart == '0') {
683       exponent += 1;
684     }
685     result += Integer<32>{exponent}.SignedDecimal();
686     return result;
687   }
688 }
689 
690 template <typename W, int P>
AsFortran(llvm::raw_ostream & o,int kind,bool minimal) const691 llvm::raw_ostream &Real<W, P>::AsFortran(
692     llvm::raw_ostream &o, int kind, bool minimal) const {
693   if (IsNotANumber()) {
694     o << "(0._" << kind << "/0.)";
695   } else if (IsInfinite()) {
696     if (IsNegative()) {
697       o << "(-1._" << kind << "/0.)";
698     } else {
699       o << "(1._" << kind << "/0.)";
700     }
701   } else {
702     using B = decimal::BinaryFloatingPointNumber<P>;
703     B value{word_.template ToUInt<typename B::RawType>()};
704     char buffer[common::MaxDecimalConversionDigits(P) +
705         EXTRA_DECIMAL_CONVERSION_SPACE];
706     decimal::DecimalConversionFlags flags{}; // default: exact representation
707     if (minimal) {
708       flags = decimal::Minimize;
709     }
710     auto result{decimal::ConvertToDecimal<P>(buffer, sizeof buffer, flags,
711         static_cast<int>(sizeof buffer), decimal::RoundNearest, value)};
712     const char *p{result.str};
713     if (DEREF(p) == '-' || *p == '+') {
714       o << *p++;
715     }
716     int expo{result.decimalExponent};
717     if (*p != '0') {
718       --expo;
719     }
720     o << *p << '.' << (p + 1);
721     if (expo != 0) {
722       o << 'e' << expo;
723     }
724     o << '_' << kind;
725   }
726   return o;
727 }
728 
729 // 16.9.180
RRSPACING() const730 template <typename W, int P> Real<W, P> Real<W, P>::RRSPACING() const {
731   if (IsNotANumber()) {
732     return *this;
733   } else if (IsInfinite()) {
734     return NotANumber();
735   } else {
736     Real result;
737     result.Normalize(false, binaryPrecision + exponentBias - 1, GetFraction());
738     return result;
739   }
740 }
741 
742 // 16.9.180
SPACING() const743 template <typename W, int P> Real<W, P> Real<W, P>::SPACING() const {
744   if (IsNotANumber()) {
745     return *this;
746   } else if (IsInfinite()) {
747     return NotANumber();
748   } else if (IsZero()) {
749     return TINY();
750   } else {
751     Real result;
752     result.Normalize(
753         false, Exponent() - binaryPrecision + 1, Fraction::MASKL(1));
754     return result;
755   }
756 }
757 
758 // 16.9.171
759 template <typename W, int P>
SET_EXPONENT(int expo) const760 Real<W, P> Real<W, P>::SET_EXPONENT(int expo) const {
761   if (IsNotANumber()) {
762     return *this;
763   } else if (IsInfinite()) {
764     return NotANumber();
765   } else if (IsZero()) {
766     return *this;
767   } else {
768     Real result;
769     result.Normalize(IsNegative(), expo + exponentBias - 1, GetFraction());
770     return result;
771   }
772 }
773 
774 // 16.9.171
FRACTION() const775 template <typename W, int P> Real<W, P> Real<W, P>::FRACTION() const {
776   return SET_EXPONENT(0);
777 }
778 
779 template class Real<Integer<16>, 11>;
780 template class Real<Integer<16>, 8>;
781 template class Real<Integer<32>, 24>;
782 template class Real<Integer<64>, 53>;
783 template class Real<Integer<80>, 64>;
784 template class Real<Integer<128>, 113>;
785 } // namespace Fortran::evaluate::value
786