xref: /llvm-project/flang/lib/Evaluate/real.cpp (revision 64ab3302d5a130c00b66a6957b2e7f0c9b9c537d) 
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 <limits>
15 
16 namespace Fortran::evaluate::value {
17 
18 template<typename W, int P> Relation Real<W, P>::Compare(const Real &y) const {
19   if (IsNotANumber() || y.IsNotANumber()) {  // NaN vs x, x vs NaN
20     return Relation::Unordered;
21   } else if (IsInfinite()) {
22     if (y.IsInfinite()) {
23       if (IsNegative()) {  // -Inf vs +/-Inf
24         return y.IsNegative() ? Relation::Equal : Relation::Less;
25       } else {  // +Inf vs +/-Inf
26         return y.IsNegative() ? Relation::Greater : Relation::Equal;
27       }
28     } else {  // +/-Inf vs finite
29       return IsNegative() ? Relation::Less : Relation::Greater;
30     }
31   } else if (y.IsInfinite()) {  // finite vs +/-Inf
32     return y.IsNegative() ? Relation::Greater : Relation::Less;
33   } else {  // two finite numbers
34     bool isNegative{IsNegative()};
35     if (isNegative != y.IsNegative()) {
36       if (word_.IOR(y.word_).IBCLR(bits - 1).IsZero()) {
37         return Relation::Equal;  // +/-0.0 == -/+0.0
38       } else {
39         return isNegative ? Relation::Less : Relation::Greater;
40       }
41     } else {
42       // same sign
43       Ordering order{evaluate::Compare(Exponent(), y.Exponent())};
44       if (order == Ordering::Equal) {
45         order = GetSignificand().CompareUnsigned(y.GetSignificand());
46       }
47       if (isNegative) {
48         order = Reverse(order);
49       }
50       return RelationFromOrdering(order);
51     }
52   }
53 }
54 
55 template<typename W, int P>
56 ValueWithRealFlags<Real<W, P>> Real<W, P>::Add(
57     const Real &y, Rounding rounding) const {
58   ValueWithRealFlags<Real> result;
59   if (IsNotANumber() || y.IsNotANumber()) {
60     result.value = NotANumber();  // NaN + x -> NaN
61     if (IsSignalingNaN() || y.IsSignalingNaN()) {
62       result.flags.set(RealFlag::InvalidArgument);
63     }
64     return result;
65   }
66   bool isNegative{IsNegative()};
67   bool yIsNegative{y.IsNegative()};
68   if (IsInfinite()) {
69     if (y.IsInfinite()) {
70       if (isNegative == yIsNegative) {
71         result.value = *this;  // +/-Inf + +/-Inf -> +/-Inf
72       } else {
73         result.value = NotANumber();  // +/-Inf + -/+Inf -> NaN
74         result.flags.set(RealFlag::InvalidArgument);
75       }
76     } else {
77       result.value = *this;  // +/-Inf + x -> +/-Inf
78     }
79     return result;
80   }
81   if (y.IsInfinite()) {
82     result.value = y;  // x + +/-Inf -> +/-Inf
83     return result;
84   }
85   int exponent{Exponent()};
86   int yExponent{y.Exponent()};
87   if (exponent < yExponent) {
88     // y is larger in magnitude; simplify by reversing operands
89     return y.Add(*this, rounding);
90   }
91   if (exponent == yExponent && isNegative != yIsNegative) {
92     Ordering order{GetSignificand().CompareUnsigned(y.GetSignificand())};
93     if (order == Ordering::Less) {
94       // Same exponent, opposite signs, and y is larger in magnitude
95       return y.Add(*this, rounding);
96     }
97     if (order == Ordering::Equal) {
98       // x + (-x) -> +0.0 unless rounding is directed downwards
99       if (rounding.mode == common::RoundingMode::Down) {
100         result.value.word_ = result.value.word_.IBSET(bits - 1);  // -0.0
101       }
102       return result;
103     }
104   }
105   // Our exponent is greater than y's, or the exponents match and y is not
106   // of the opposite sign and greater magnitude.  So (x+y) will have the
107   // same sign as x.
108   Fraction fraction{GetFraction()};
109   Fraction yFraction{y.GetFraction()};
110   int rshift = exponent - yExponent;
111   if (exponent > 0 && yExponent == 0) {
112     --rshift;  // correct overshift when only y is subnormal
113   }
114   RoundingBits roundingBits{yFraction, rshift};
115   yFraction = yFraction.SHIFTR(rshift);
116   bool carry{false};
117   if (isNegative != yIsNegative) {
118     // Opposite signs: subtract via addition of two's complement of y and
119     // the rounding bits.
120     yFraction = yFraction.NOT();
121     carry = roundingBits.Negate();
122   }
123   auto sum{fraction.AddUnsigned(yFraction, carry)};
124   fraction = sum.value;
125   if (isNegative == yIsNegative && sum.carry) {
126     roundingBits.ShiftRight(sum.value.BTEST(0));
127     fraction = fraction.SHIFTR(1).IBSET(fraction.bits - 1);
128     ++exponent;
129   }
130   NormalizeAndRound(
131       result, isNegative, exponent, fraction, rounding, roundingBits);
132   return result;
133 }
134 
135 template<typename W, int P>
136 ValueWithRealFlags<Real<W, P>> Real<W, P>::Multiply(
137     const Real &y, Rounding rounding) const {
138   ValueWithRealFlags<Real> result;
139   if (IsNotANumber() || y.IsNotANumber()) {
140     result.value = NotANumber();  // NaN * x -> NaN
141     if (IsSignalingNaN() || y.IsSignalingNaN()) {
142       result.flags.set(RealFlag::InvalidArgument);
143     }
144   } else {
145     bool isNegative{IsNegative() != y.IsNegative()};
146     if (IsInfinite() || y.IsInfinite()) {
147       if (IsZero() || y.IsZero()) {
148         result.value = NotANumber();  // 0 * Inf -> NaN
149         result.flags.set(RealFlag::InvalidArgument);
150       } else {
151         result.value = Infinity(isNegative);
152       }
153     } else {
154       auto product{GetFraction().MultiplyUnsigned(y.GetFraction())};
155       std::int64_t exponent{CombineExponents(y, false)};
156       if (exponent < 1) {
157         int rshift = 1 - exponent;
158         exponent = 1;
159         bool sticky{false};
160         if (rshift >= product.upper.bits + product.lower.bits) {
161           sticky = !product.lower.IsZero() || !product.upper.IsZero();
162         } else if (rshift >= product.lower.bits) {
163           sticky = !product.lower.IsZero() ||
164               !product.upper
165                    .IAND(product.upper.MASKR(rshift - product.lower.bits))
166                    .IsZero();
167         } else {
168           sticky = !product.lower.IAND(product.lower.MASKR(rshift)).IsZero();
169         }
170         product.lower = product.lower.SHIFTRWithFill(product.upper, rshift);
171         product.upper = product.upper.SHIFTR(rshift);
172         if (sticky) {
173           product.lower = product.lower.IBSET(0);
174         }
175       }
176       int leadz{product.upper.LEADZ()};
177       if (leadz >= product.upper.bits) {
178         leadz += product.lower.LEADZ();
179       }
180       int lshift{leadz};
181       if (lshift > exponent - 1) {
182         lshift = exponent - 1;
183       }
184       exponent -= lshift;
185       product.upper = product.upper.SHIFTLWithFill(product.lower, lshift);
186       product.lower = product.lower.SHIFTL(lshift);
187       RoundingBits roundingBits{product.lower, product.lower.bits};
188       NormalizeAndRound(result, isNegative, exponent, product.upper, rounding,
189           roundingBits, true /*multiply*/);
190     }
191   }
192   return result;
193 }
194 
195 template<typename W, int P>
196 ValueWithRealFlags<Real<W, P>> Real<W, P>::Divide(
197     const Real &y, Rounding rounding) const {
198   ValueWithRealFlags<Real> result;
199   if (IsNotANumber() || y.IsNotANumber()) {
200     result.value = NotANumber();  // NaN / x -> NaN, x / NaN -> NaN
201     if (IsSignalingNaN() || y.IsSignalingNaN()) {
202       result.flags.set(RealFlag::InvalidArgument);
203     }
204   } else {
205     bool isNegative{IsNegative() != y.IsNegative()};
206     if (IsInfinite()) {
207       if (y.IsInfinite()) {
208         result.value = NotANumber();  // Inf/Inf -> NaN
209         result.flags.set(RealFlag::InvalidArgument);
210       } else {  // Inf/x -> Inf,  Inf/0 -> Inf
211         result.value = Infinity(isNegative);
212       }
213     } else if (y.IsZero()) {
214       if (IsZero()) {  // 0/0 -> NaN
215         result.value = NotANumber();
216         result.flags.set(RealFlag::InvalidArgument);
217       } else {  // x/0 -> Inf, Inf/0 -> Inf
218         result.value = Infinity(isNegative);
219         result.flags.set(RealFlag::DivideByZero);
220       }
221     } else if (IsZero() || y.IsInfinite()) {  // 0/x, x/Inf -> 0
222       if (isNegative) {
223         result.value.word_ = result.value.word_.IBSET(bits - 1);
224       }
225     } else {
226       // dividend and divisor are both finite and nonzero numbers
227       Fraction top{GetFraction()}, divisor{y.GetFraction()};
228       std::int64_t exponent{CombineExponents(y, true)};
229       Fraction quotient;
230       bool msb{false};
231       if (!top.BTEST(top.bits - 1) || !divisor.BTEST(divisor.bits - 1)) {
232         // One or two subnormals
233         int topLshift{top.LEADZ()};
234         top = top.SHIFTL(topLshift);
235         int divisorLshift{divisor.LEADZ()};
236         divisor = divisor.SHIFTL(divisorLshift);
237         exponent += divisorLshift - topLshift;
238       }
239       for (int j{1}; j <= quotient.bits; ++j) {
240         if (NextQuotientBit(top, msb, divisor)) {
241           quotient = quotient.IBSET(quotient.bits - j);
242         }
243       }
244       bool guard{NextQuotientBit(top, msb, divisor)};
245       bool round{NextQuotientBit(top, msb, divisor)};
246       bool sticky{msb || !top.IsZero()};
247       RoundingBits roundingBits{guard, round, sticky};
248       if (exponent < 1) {
249         std::int64_t rshift{1 - exponent};
250         for (; rshift > 0; --rshift) {
251           roundingBits.ShiftRight(quotient.BTEST(0));
252           quotient = quotient.SHIFTR(1);
253         }
254         exponent = 1;
255       }
256       NormalizeAndRound(
257           result, isNegative, exponent, quotient, rounding, roundingBits);
258     }
259   }
260   return result;
261 }
262 
263 template<typename W, int P>
264 ValueWithRealFlags<Real<W, P>> Real<W, P>::ToWholeNumber(
265     common::RoundingMode mode) const {
266   ValueWithRealFlags<Real> result{*this};
267   if (IsNotANumber()) {
268     result.flags.set(RealFlag::InvalidArgument);
269     result.value = NotANumber();
270   } else if (IsInfinite()) {
271     result.flags.set(RealFlag::Overflow);
272   } else {
273     constexpr int noClipExponent{exponentBias + binaryPrecision - 1};
274     if (Exponent() < noClipExponent) {
275       Real adjust;  // ABS(EPSILON(adjust)) == 0.5
276       adjust.Normalize(IsSignBitSet(), noClipExponent, Fraction::MASKL(1));
277       // Compute ival=(*this + adjust), losing any fractional bits; keep flags
278       result = Add(adjust, Rounding{mode});
279       result.flags.reset(RealFlag::Inexact);  // result *is* exact
280       // Return (ival-adjust) with original sign in case we've generated a zero.
281       result.value =
282           result.value.Subtract(adjust, Rounding{common::RoundingMode::ToZero})
283               .value.SIGN(*this);
284     }
285   }
286   return result;
287 }
288 
289 template<typename W, int P>
290 RealFlags Real<W, P>::Normalize(bool negative, int exponent,
291     const Fraction &fraction, Rounding rounding, RoundingBits *roundingBits) {
292   int lshift{fraction.LEADZ()};
293   if (lshift == fraction.bits /* fraction is zero */ &&
294       (!roundingBits || roundingBits->empty())) {
295     // No fraction, no rounding bits -> +/-0.0
296     exponent = lshift = 0;
297   } else if (lshift < exponent) {
298     exponent -= lshift;
299   } else if (exponent > 0) {
300     lshift = exponent - 1;
301     exponent = 0;
302   } else if (lshift == 0) {
303     exponent = 1;
304   } else {
305     lshift = 0;
306   }
307   if (exponent >= maxExponent) {
308     // Infinity or overflow
309     if (rounding.mode == common::RoundingMode::TiesToEven ||
310         rounding.mode == common::RoundingMode::TiesAwayFromZero ||
311         (rounding.mode == common::RoundingMode::Up && !negative) ||
312         (rounding.mode == common::RoundingMode::Down && negative)) {
313       word_ = Word{maxExponent}.SHIFTL(significandBits);  // Inf
314     } else {
315       // directed rounding: round to largest finite value rather than infinity
316       // (x86 does this, not sure whether it's standard behavior)
317       word_ = Word{word_.MASKR(word_.bits - 1)}.IBCLR(significandBits);
318     }
319     if (negative) {
320       word_ = word_.IBSET(bits - 1);
321     }
322     RealFlags flags{RealFlag::Overflow};
323     if (!fraction.IsZero()) {
324       flags.set(RealFlag::Inexact);
325     }
326     return flags;
327   }
328   word_ = Word::ConvertUnsigned(fraction).value;
329   if (lshift > 0) {
330     word_ = word_.SHIFTL(lshift);
331     if (roundingBits) {
332       for (; lshift > 0; --lshift) {
333         if (roundingBits->ShiftLeft()) {
334           word_ = word_.IBSET(lshift - 1);
335         }
336       }
337     }
338   }
339   if constexpr (isImplicitMSB) {
340     word_ = word_.IBCLR(significandBits);
341   }
342   word_ = word_.IOR(Word{exponent}.SHIFTL(significandBits));
343   if (negative) {
344     word_ = word_.IBSET(bits - 1);
345   }
346   return {};
347 }
348 
349 template<typename W, int P>
350 RealFlags Real<W, P>::Round(
351     Rounding rounding, const RoundingBits &bits, bool multiply) {
352   int origExponent{Exponent()};
353   RealFlags flags;
354   bool inexact{!bits.empty()};
355   if (inexact) {
356     flags.set(RealFlag::Inexact);
357   }
358   if (origExponent < maxExponent &&
359       bits.MustRound(rounding, IsNegative(), word_.BTEST(0) /* is odd */)) {
360     typename Fraction::ValueWithCarry sum{
361         GetFraction().AddUnsigned(Fraction{}, true)};
362     int newExponent{origExponent};
363     if (sum.carry) {
364       // The fraction was all ones before rounding; sum.value is now zero
365       sum.value = sum.value.IBSET(binaryPrecision - 1);
366       if (++newExponent >= maxExponent) {
367         flags.set(RealFlag::Overflow);  // rounded away to an infinity
368       }
369     }
370     flags |= Normalize(IsNegative(), newExponent, sum.value);
371   }
372   if (inexact && origExponent == 0) {
373     // inexact subnormal input: signal Underflow unless in an x86-specific
374     // edge case
375     if (rounding.x86CompatibleBehavior && Exponent() != 0 && multiply &&
376         bits.sticky() &&
377         (bits.guard() ||
378             (rounding.mode != common::RoundingMode::Up &&
379                 rounding.mode != common::RoundingMode::Down))) {
380       // x86 edge case in which Underflow fails to signal when a subnormal
381       // inexact multiplication product rounds to a normal result when
382       // the guard bit is set or we're not using directed rounding
383     } else {
384       flags.set(RealFlag::Underflow);
385     }
386   }
387   return flags;
388 }
389 
390 template<typename W, int P>
391 void Real<W, P>::NormalizeAndRound(ValueWithRealFlags<Real> &result,
392     bool isNegative, int exponent, const Fraction &fraction, Rounding rounding,
393     RoundingBits roundingBits, bool multiply) {
394   result.flags |= result.value.Normalize(
395       isNegative, exponent, fraction, rounding, &roundingBits);
396   result.flags |= result.value.Round(rounding, roundingBits, multiply);
397 }
398 
399 inline enum decimal::FortranRounding MapRoundingMode(
400     common::RoundingMode rounding) {
401   switch (rounding) {
402   case common::RoundingMode::TiesToEven: break;
403   case common::RoundingMode::ToZero: return decimal::RoundToZero;
404   case common::RoundingMode::Down: return decimal::RoundDown;
405   case common::RoundingMode::Up: return decimal::RoundUp;
406   case common::RoundingMode::TiesAwayFromZero: return decimal::RoundCompatible;
407   }
408   return decimal::RoundNearest;  // dodge gcc warning about lack of result
409 }
410 
411 inline RealFlags MapFlags(decimal::ConversionResultFlags flags) {
412   RealFlags result;
413   if (flags & decimal::Overflow) {
414     result.set(RealFlag::Overflow);
415   }
416   if (flags & decimal::Inexact) {
417     result.set(RealFlag::Inexact);
418   }
419   if (flags & decimal::Invalid) {
420     result.set(RealFlag::InvalidArgument);
421   }
422   return result;
423 }
424 
425 template<typename W, int P>
426 ValueWithRealFlags<Real<W, P>> Real<W, P>::Read(
427     const char *&p, Rounding rounding) {
428   auto converted{
429       decimal::ConvertToBinary<P>(p, MapRoundingMode(rounding.mode))};
430   const auto *value{reinterpret_cast<Real<W, P> *>(&converted.binary)};
431   return {*value, MapFlags(converted.flags)};
432 }
433 
434 template<typename W, int P> std::string Real<W, P>::DumpHexadecimal() const {
435   if (IsNotANumber()) {
436     return "NaN 0x"s + word_.Hexadecimal();
437   } else if (IsNegative()) {
438     return "-"s + Negate().DumpHexadecimal();
439   } else if (IsInfinite()) {
440     return "Inf"s;
441   } else if (IsZero()) {
442     return "0.0"s;
443   } else {
444     Fraction frac{GetFraction()};
445     std::string result{"0x"};
446     char intPart = '0' + frac.BTEST(frac.bits - 1);
447     result += intPart;
448     result += '.';
449     int trailz{frac.TRAILZ()};
450     if (trailz >= frac.bits - 1) {
451       result += '0';
452     } else {
453       int remainingBits{frac.bits - 1 - trailz};
454       int wholeNybbles{remainingBits / 4};
455       int lostBits{remainingBits - 4 * wholeNybbles};
456       if (wholeNybbles > 0) {
457         std::string fracHex{frac.SHIFTR(trailz + lostBits)
458                                 .IAND(frac.MASKR(4 * wholeNybbles))
459                                 .Hexadecimal()};
460         std::size_t field = wholeNybbles;
461         if (fracHex.size() < field) {
462           result += std::string(field - fracHex.size(), '0');
463         }
464         result += fracHex;
465       }
466       if (lostBits > 0) {
467         result += frac.SHIFTR(trailz)
468                       .IAND(frac.MASKR(lostBits))
469                       .SHIFTL(4 - lostBits)
470                       .Hexadecimal();
471       }
472     }
473     result += 'p';
474     int exponent = Exponent() - exponentBias;
475     result += Integer<32>{exponent}.SignedDecimal();
476     return result;
477   }
478 }
479 
480 template<typename W, int P>
481 std::ostream &Real<W, P>::AsFortran(
482     std::ostream &o, int kind, bool minimal) const {
483   if (IsNotANumber()) {
484     o << "(0._" << kind << "/0.)";
485   } else if (IsInfinite()) {
486     if (IsNegative()) {
487       o << "(-1._" << kind << "/0.)";
488     } else {
489       o << "(1._" << kind << "/0.)";
490     }
491   } else {
492     using B = decimal::BinaryFloatingPointNumber<P>;
493     const auto *value{reinterpret_cast<const B *>(this)};
494     char buffer[24000];  // accommodate real*16
495     decimal::DecimalConversionFlags flags{};  // default: exact representation
496     if (minimal) {
497       flags = decimal::Minimize;
498     }
499     auto result{decimal::ConvertToDecimal<P>(buffer, sizeof buffer, flags,
500         static_cast<int>(sizeof buffer), decimal::RoundNearest, *value)};
501     const char *p{result.str};
502     if (DEREF(p) == '-' || *p == '+') {
503       o << *p++;
504     }
505     int expo{result.decimalExponent};
506     if (*p != '0') {
507       --expo;
508     }
509     o << *p << '.' << (p + 1);
510     if (expo != 0) {
511       o << 'e' << expo;
512     }
513     o << '_' << kind;
514   }
515   return o;
516 }
517 
518 template class Real<Integer<16>, 11>;
519 template class Real<Integer<16>, 8>;
520 template class Real<Integer<32>, 24>;
521 template class Real<Integer<64>, 53>;
522 template class Real<Integer<80>, 64>;
523 template class Real<Integer<128>, 112>;
524 }
525