flang/lib/Evaluate/real.cpp

   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
  19 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>
  57 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>
 137 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>
 197 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>
 265 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>
 331 ValueWithRealFlags<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)
 381 template <typename W, int P>
 382 ValueWithRealFlags<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.
 416 template <typename W, int P>
 417 ValueWithRealFlags<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.
 453 template <typename W, int P>
 454 ValueWithRealFlags<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
 468 template <typename W, int P>
 469 ValueWithRealFlags<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
 483 template <typename W, int P>
 484 ValueWithRealFlags<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
 509 template <typename W, int P>
 510 RealFlags 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
 575 template <typename W, int P>
 576 RealFlags 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
 616 template <typename W, int P>
 617 void 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
 625 inline 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
 642 inline 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
 656 template <typename W, int P>
 657 ValueWithRealFlags<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
 665 template <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
 714 template <typename W, int P>
 715 llvm::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
 754 template <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
 767 template <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
 784 template <typename W, int P>
 785 Real<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
 798 template <typename W, int P> Real<W, P> Real<W, P>::FRACTION() const {
 799   return SET_EXPONENT(0);
 800 }
 801
 802 template class Real<Integer<16>, 11>;
 803 template class Real<Integer<16>, 8>;
 804 template class Real<Integer<32>, 24>;
 805 template class Real<Integer<64>, 53>;
 806 template class Real<X87IntegerContainer, 64>;
 807 template class Real<Integer<128>, 113>;
 808 } // namespace Fortran::evaluate::value