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.word_ = result.value.word_.IBSET(bits - 1); // -0.0
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.word_ = result.value.word_.IBSET(bits - 1);
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>
ToWholeNumber(common::RoundingMode mode) const265 ValueWithRealFlags<Real<W, P>> Real<W, P>::ToWholeNumber(
266 common::RoundingMode mode) const {
267 ValueWithRealFlags<Real> result{*this};
268 if (IsNotANumber()) {
269 result.flags.set(RealFlag::InvalidArgument);
270 result.value = NotANumber();
271 } else if (IsInfinite()) {
272 result.flags.set(RealFlag::Overflow);
273 } else {
274 constexpr int noClipExponent{exponentBias + binaryPrecision - 1};
275 if (Exponent() < noClipExponent) {
276 Real adjust; // ABS(EPSILON(adjust)) == 0.5
277 adjust.Normalize(IsSignBitSet(), noClipExponent, Fraction::MASKL(1));
278 // Compute ival=(*this + adjust), losing any fractional bits; keep flags
279 result = Add(adjust, Rounding{mode});
280 result.flags.reset(RealFlag::Inexact); // result *is* exact
281 // Return (ival-adjust) with original sign in case we've generated a zero.
282 result.value =
283 result.value.Subtract(adjust, Rounding{common::RoundingMode::ToZero})
284 .value.SIGN(*this);
285 }
286 }
287 return result;
288 }
289
290 template <typename W, int P>
Normalize(bool negative,int exponent,const Fraction & fraction,Rounding rounding,RoundingBits * roundingBits)291 RealFlags Real<W, P>::Normalize(bool negative, int exponent,
292 const Fraction &fraction, Rounding rounding, RoundingBits *roundingBits) {
293 int lshift{fraction.LEADZ()};
294 if (lshift == fraction.bits /* fraction is zero */ &&
295 (!roundingBits || roundingBits->empty())) {
296 // No fraction, no rounding bits -> +/-0.0
297 exponent = lshift = 0;
298 } else if (lshift < exponent) {
299 exponent -= lshift;
300 } else if (exponent > 0) {
301 lshift = exponent - 1;
302 exponent = 0;
303 } else if (lshift == 0) {
304 exponent = 1;
305 } else {
306 lshift = 0;
307 }
308 if (exponent >= maxExponent) {
309 // Infinity or overflow
310 if (rounding.mode == common::RoundingMode::TiesToEven ||
311 rounding.mode == common::RoundingMode::TiesAwayFromZero ||
312 (rounding.mode == common::RoundingMode::Up && !negative) ||
313 (rounding.mode == common::RoundingMode::Down && negative)) {
314 word_ = Word{maxExponent}.SHIFTL(significandBits); // Inf
315 } else {
316 // directed rounding: round to largest finite value rather than infinity
317 // (x86 does this, not sure whether it's standard behavior)
318 word_ = Word{word_.MASKR(word_.bits - 1)}.IBCLR(significandBits);
319 }
320 if (negative) {
321 word_ = word_.IBSET(bits - 1);
322 }
323 RealFlags flags{RealFlag::Overflow};
324 if (!fraction.IsZero()) {
325 flags.set(RealFlag::Inexact);
326 }
327 return flags;
328 }
329 word_ = Word::ConvertUnsigned(fraction).value;
330 if (lshift > 0) {
331 word_ = word_.SHIFTL(lshift);
332 if (roundingBits) {
333 for (; lshift > 0; --lshift) {
334 if (roundingBits->ShiftLeft()) {
335 word_ = word_.IBSET(lshift - 1);
336 }
337 }
338 }
339 }
340 if constexpr (isImplicitMSB) {
341 word_ = word_.IBCLR(significandBits);
342 }
343 word_ = word_.IOR(Word{exponent}.SHIFTL(significandBits));
344 if (negative) {
345 word_ = word_.IBSET(bits - 1);
346 }
347 return {};
348 }
349
350 template <typename W, int P>
Round(Rounding rounding,const RoundingBits & bits,bool multiply)351 RealFlags Real<W, P>::Round(
352 Rounding rounding, const RoundingBits &bits, bool multiply) {
353 int origExponent{Exponent()};
354 RealFlags flags;
355 bool inexact{!bits.empty()};
356 if (inexact) {
357 flags.set(RealFlag::Inexact);
358 }
359 if (origExponent < maxExponent &&
360 bits.MustRound(rounding, IsNegative(), word_.BTEST(0) /* is odd */)) {
361 typename Fraction::ValueWithCarry sum{
362 GetFraction().AddUnsigned(Fraction{}, true)};
363 int newExponent{origExponent};
364 if (sum.carry) {
365 // The fraction was all ones before rounding; sum.value is now zero
366 sum.value = sum.value.IBSET(binaryPrecision - 1);
367 if (++newExponent >= maxExponent) {
368 flags.set(RealFlag::Overflow); // rounded away to an infinity
369 }
370 }
371 flags |= Normalize(IsNegative(), newExponent, sum.value);
372 }
373 if (inexact && origExponent == 0) {
374 // inexact subnormal input: signal Underflow unless in an x86-specific
375 // edge case
376 if (rounding.x86CompatibleBehavior && Exponent() != 0 && multiply &&
377 bits.sticky() &&
378 (bits.guard() ||
379 (rounding.mode != common::RoundingMode::Up &&
380 rounding.mode != common::RoundingMode::Down))) {
381 // x86 edge case in which Underflow fails to signal when a subnormal
382 // inexact multiplication product rounds to a normal result when
383 // the guard bit is set or we're not using directed rounding
384 } else {
385 flags.set(RealFlag::Underflow);
386 }
387 }
388 return flags;
389 }
390
391 template <typename W, int P>
NormalizeAndRound(ValueWithRealFlags<Real> & result,bool isNegative,int exponent,const Fraction & fraction,Rounding rounding,RoundingBits roundingBits,bool multiply)392 void Real<W, P>::NormalizeAndRound(ValueWithRealFlags<Real> &result,
393 bool isNegative, int exponent, const Fraction &fraction, Rounding rounding,
394 RoundingBits roundingBits, bool multiply) {
395 result.flags |= result.value.Normalize(
396 isNegative, exponent, fraction, rounding, &roundingBits);
397 result.flags |= result.value.Round(rounding, roundingBits, multiply);
398 }
399
MapRoundingMode(common::RoundingMode rounding)400 inline enum decimal::FortranRounding MapRoundingMode(
401 common::RoundingMode rounding) {
402 switch (rounding) {
403 case common::RoundingMode::TiesToEven:
404 break;
405 case common::RoundingMode::ToZero:
406 return decimal::RoundToZero;
407 case common::RoundingMode::Down:
408 return decimal::RoundDown;
409 case common::RoundingMode::Up:
410 return decimal::RoundUp;
411 case common::RoundingMode::TiesAwayFromZero:
412 return decimal::RoundCompatible;
413 }
414 return decimal::RoundNearest; // dodge gcc warning about lack of result
415 }
416
MapFlags(decimal::ConversionResultFlags flags)417 inline RealFlags MapFlags(decimal::ConversionResultFlags flags) {
418 RealFlags result;
419 if (flags & decimal::Overflow) {
420 result.set(RealFlag::Overflow);
421 }
422 if (flags & decimal::Inexact) {
423 result.set(RealFlag::Inexact);
424 }
425 if (flags & decimal::Invalid) {
426 result.set(RealFlag::InvalidArgument);
427 }
428 return result;
429 }
430
431 template <typename W, int P>
Read(const char * & p,Rounding rounding)432 ValueWithRealFlags<Real<W, P>> Real<W, P>::Read(
433 const char *&p, Rounding rounding) {
434 auto converted{
435 decimal::ConvertToBinary<P>(p, MapRoundingMode(rounding.mode))};
436 const auto *value{reinterpret_cast<Real<W, P> *>(&converted.binary)};
437 return {*value, MapFlags(converted.flags)};
438 }
439
DumpHexadecimal() const440 template <typename W, int P> std::string Real<W, P>::DumpHexadecimal() const {
441 if (IsNotANumber()) {
442 return "NaN 0x"s + word_.Hexadecimal();
443 } else if (IsNegative()) {
444 return "-"s + Negate().DumpHexadecimal();
445 } else if (IsInfinite()) {
446 return "Inf"s;
447 } else if (IsZero()) {
448 return "0.0"s;
449 } else {
450 Fraction frac{GetFraction()};
451 std::string result{"0x"};
452 char intPart = '0' + frac.BTEST(frac.bits - 1);
453 result += intPart;
454 result += '.';
455 int trailz{frac.TRAILZ()};
456 if (trailz >= frac.bits - 1) {
457 result += '0';
458 } else {
459 int remainingBits{frac.bits - 1 - trailz};
460 int wholeNybbles{remainingBits / 4};
461 int lostBits{remainingBits - 4 * wholeNybbles};
462 if (wholeNybbles > 0) {
463 std::string fracHex{frac.SHIFTR(trailz + lostBits)
464 .IAND(frac.MASKR(4 * wholeNybbles))
465 .Hexadecimal()};
466 std::size_t field = wholeNybbles;
467 if (fracHex.size() < field) {
468 result += std::string(field - fracHex.size(), '0');
469 }
470 result += fracHex;
471 }
472 if (lostBits > 0) {
473 result += frac.SHIFTR(trailz)
474 .IAND(frac.MASKR(lostBits))
475 .SHIFTL(4 - lostBits)
476 .Hexadecimal();
477 }
478 }
479 result += 'p';
480 int exponent = Exponent() - exponentBias;
481 result += Integer<32>{exponent}.SignedDecimal();
482 return result;
483 }
484 }
485
486 template <typename W, int P>
AsFortran(llvm::raw_ostream & o,int kind,bool minimal) const487 llvm::raw_ostream &Real<W, P>::AsFortran(
488 llvm::raw_ostream &o, int kind, bool minimal) const {
489 if (IsNotANumber()) {
490 o << "(0._" << kind << "/0.)";
491 } else if (IsInfinite()) {
492 if (IsNegative()) {
493 o << "(-1._" << kind << "/0.)";
494 } else {
495 o << "(1._" << kind << "/0.)";
496 }
497 } else {
498 using B = decimal::BinaryFloatingPointNumber<P>;
499 B value{word_.template ToUInt<typename B::RawType>()};
500 char buffer[common::MaxDecimalConversionDigits(P) +
501 EXTRA_DECIMAL_CONVERSION_SPACE];
502 decimal::DecimalConversionFlags flags{}; // default: exact representation
503 if (minimal) {
504 flags = decimal::Minimize;
505 }
506 auto result{decimal::ConvertToDecimal<P>(buffer, sizeof buffer, flags,
507 static_cast<int>(sizeof buffer), decimal::RoundNearest, value)};
508 const char *p{result.str};
509 if (DEREF(p) == '-' || *p == '+') {
510 o << *p++;
511 }
512 int expo{result.decimalExponent};
513 if (*p != '0') {
514 --expo;
515 }
516 o << *p << '.' << (p + 1);
517 if (expo != 0) {
518 o << 'e' << expo;
519 }
520 o << '_' << kind;
521 }
522 return o;
523 }
524
525 template class Real<Integer<16>, 11>;
526 template class Real<Integer<16>, 8>;
527 template class Real<Integer<32>, 24>;
528 template class Real<Integer<64>, 53>;
529 template class Real<Integer<80>, 64>;
530 template class Real<Integer<128>, 113>;
531 } // namespace Fortran::evaluate::value
532