817 lines · 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-exception6//7//===----------------------------------------------------------------------===//8 9#include "flang/Evaluate/real.h"10#include "int-power.h"11#include "flang/Common/idioms.h"12#include "flang/Decimal/decimal.h"13#include "flang/Parser/characters.h"14#include "llvm/Support/raw_ostream.h"15#include <limits>16 17namespace Fortran::evaluate::value {18 19template <typename W, int P> Relation Real<W, P>::Compare(const Real &y) const {20 if (IsNotANumber() || y.IsNotANumber()) { // NaN vs x, x vs NaN21 return Relation::Unordered;22 } else if (IsInfinite()) {23 if (y.IsInfinite()) {24 if (IsNegative()) { // -Inf vs +/-Inf25 return y.IsNegative() ? Relation::Equal : Relation::Less;26 } else { // +Inf vs +/-Inf27 return y.IsNegative() ? Relation::Greater : Relation::Equal;28 }29 } else { // +/-Inf vs finite30 return IsNegative() ? Relation::Less : Relation::Greater;31 }32 } else if (y.IsInfinite()) { // finite vs +/-Inf33 return y.IsNegative() ? Relation::Greater : Relation::Less;34 } else { // two finite numbers35 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.039 } else {40 return isNegative ? Relation::Less : Relation::Greater;41 }42 } else {43 // same sign44 Ordering order{evaluate::Compare(Exponent(), y.Exponent())};45 if (order == Ordering::Equal) {46 order = GetSignificand().CompareUnsigned(y.GetSignificand());47 }48 if (isNegative) {49 order = Reverse(order);50 }51 return RelationFromOrdering(order);52 }53 }54}55 56template <typename W, int P>57ValueWithRealFlags<Real<W, P>> Real<W, P>::Add(58 const Real &y, Rounding rounding) const {59 ValueWithRealFlags<Real> result;60 if (IsNotANumber() || y.IsNotANumber()) {61 result.value = NotANumber(); // NaN + x -> NaN62 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 -> +/-Inf73 } else {74 result.value = NotANumber(); // +/-Inf + -/+Inf -> NaN75 result.flags.set(RealFlag::InvalidArgument);76 }77 } else {78 result.value = *this; // +/-Inf + x -> +/-Inf79 }80 return result;81 }82 if (y.IsInfinite()) {83 result.value = y; // x + +/-Inf -> +/-Inf84 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 operands90 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 magnitude96 return y.Add(*this, rounding);97 }98 if (order == Ordering::Equal) {99 // x + (-x) -> +0.0 unless rounding is directed downwards100 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 not107 // of the opposite sign and greater magnitude. So (x+y) will have the108 // 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 subnormal114 }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 and120 // the rounding bits.121 yFraction = yFraction.NOT();122 carry = roundingBits.Negate();123 }124 auto sum{fraction.AddUnsigned(yFraction, carry)};125 fraction = sum.value;126 if (isNegative == yIsNegative && sum.carry) {127 roundingBits.ShiftRight(sum.value.BTEST(0));128 fraction = fraction.SHIFTR(1).IBSET(fraction.bits - 1);129 ++exponent;130 }131 NormalizeAndRound(132 result, isNegative, exponent, fraction, rounding, roundingBits);133 return result;134}135 136template <typename W, int P>137ValueWithRealFlags<Real<W, P>> Real<W, P>::Multiply(138 const Real &y, Rounding rounding) const {139 ValueWithRealFlags<Real> result;140 if (IsNotANumber() || y.IsNotANumber()) {141 result.value = NotANumber(); // NaN * x -> NaN142 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 -> NaN150 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.upper166 .IAND(product.upper.MASKR(rshift - product.lower.bits))167 .IsZero();168 } else {169 sticky = !product.lower.IAND(product.lower.MASKR(rshift)).IsZero();170 }171 product.lower = product.lower.SHIFTRWithFill(product.upper, rshift);172 product.upper = product.upper.SHIFTR(rshift);173 if (sticky) {174 product.lower = product.lower.IBSET(0);175 }176 }177 int leadz{product.upper.LEADZ()};178 if (leadz >= product.upper.bits) {179 leadz += product.lower.LEADZ();180 }181 int lshift{leadz};182 if (lshift > exponent - 1) {183 lshift = exponent - 1;184 }185 exponent -= lshift;186 product.upper = product.upper.SHIFTLWithFill(product.lower, lshift);187 product.lower = product.lower.SHIFTL(lshift);188 RoundingBits roundingBits{product.lower, product.lower.bits};189 NormalizeAndRound(result, isNegative, exponent, product.upper, rounding,190 roundingBits, true /*multiply*/);191 }192 }193 return result;194}195 196template <typename W, int P>197ValueWithRealFlags<Real<W, P>> Real<W, P>::Divide(198 const Real &y, Rounding rounding) const {199 ValueWithRealFlags<Real> result;200 if (IsNotANumber() || y.IsNotANumber()) {201 result.value = NotANumber(); // NaN / x -> NaN, x / NaN -> NaN202 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 -> NaN210 result.flags.set(RealFlag::InvalidArgument);211 } else { // Inf/x -> Inf, Inf/0 -> Inf212 result.value = Infinity(isNegative);213 }214 } else if (y.IsZero()) {215 if (IsZero()) { // 0/0 -> NaN216 result.value = NotANumber();217 result.flags.set(RealFlag::InvalidArgument);218 } else { // x/0 -> Inf, Inf/0 -> Inf219 result.value = Infinity(isNegative);220 result.flags.set(RealFlag::DivideByZero);221 }222 } else if (IsZero() || y.IsInfinite()) { // 0/x, x/Inf -> 0223 if (isNegative) {224 result.value = NegativeZero();225 }226 } else {227 // dividend and divisor are both finite and nonzero numbers228 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 subnormals234 int topLshift{top.LEADZ()};235 top = top.SHIFTL(topLshift);236 int divisorLshift{divisor.LEADZ()};237 divisor = divisor.SHIFTL(divisorLshift);238 exponent += divisorLshift - topLshift;239 }240 for (int j{1}; j <= quotient.bits; ++j) {241 if (NextQuotientBit(top, msb, divisor)) {242 quotient = quotient.IBSET(quotient.bits - j);243 }244 }245 bool guard{NextQuotientBit(top, msb, divisor)};246 bool round{NextQuotientBit(top, msb, divisor)};247 bool sticky{msb || !top.IsZero()};248 RoundingBits roundingBits{guard, round, sticky};249 if (exponent < 1) {250 std::int64_t rshift{1 - exponent};251 for (; rshift > 0; --rshift) {252 roundingBits.ShiftRight(quotient.BTEST(0));253 quotient = quotient.SHIFTR(1);254 }255 exponent = 1;256 }257 NormalizeAndRound(258 result, isNegative, exponent, quotient, rounding, roundingBits);259 }260 }261 return result;262}263 264template <typename W, int P>265ValueWithRealFlags<Real<W, P>> Real<W, P>::SQRT(Rounding rounding) const {266 ValueWithRealFlags<Real> result;267 if (IsNotANumber()) {268 result.value = NotANumber();269 if (IsSignalingNaN()) {270 result.flags.set(RealFlag::InvalidArgument);271 }272 } else if (IsNegative()) {273 if (IsZero()) {274 // SQRT(-0) == -0 in IEEE-754.275 result.value = NegativeZero();276 } else {277 result.flags.set(RealFlag::InvalidArgument);278 result.value = NotANumber();279 }280 } else if (IsInfinite()) {281 // SQRT(+Inf) == +Inf282 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 power289 // of four to avoid trouble with very large and very small values290 // (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 square306 // does not exceed the goal. Use two extra bits of precision for307 // 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 result320 }321 }322 RoundingBits roundingBits{extFrac.BTEST(1), extFrac.BTEST(0), sticky};323 NormalizeAndRound(result, false, exponentBias,324 Fraction::ConvertUnsigned(extFrac.SHIFTR(2)).value, rounding,325 roundingBits);326 }327 return result;328}329 330template <typename W, int P>331ValueWithRealFlags<Real<W, P>> Real<W, P>::NEAREST(bool upward) const {332 ValueWithRealFlags<Real> result;333 bool isNegative{IsNegative()};334 if (IsFinite()) {335 Fraction fraction{GetFraction()};336 int expo{Exponent()};337 Fraction one{1};338 Fraction nearest;339 if (upward != isNegative) { // upward in magnitude340 auto next{fraction.AddUnsigned(one)};341 if (next.carry) {342 ++expo;343 nearest = Fraction::Least(); // MSB only344 } else {345 nearest = next.value;346 }347 } else { // downward in magnitude348 if (IsZero()) {349 nearest = 1; // smallest magnitude negative subnormal350 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 finite366 } else {367 result.value = *this;368 }369 } else { // NaN370 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 intermediate377// values are susceptible to over/underflow when computed naively.378// Assuming that x>=y, calculate instead:379// HYPOT(x,y) = SQRT(x**2 * (1+(y/x)**2))380// = ABS(x) * SQRT(1+(y/x)**2)381template <typename W, int P>382ValueWithRealFlags<Real<W, P>> Real<W, P>::HYPOT(383 const Real &y, Rounding rounding) const {384 ValueWithRealFlags<Real> result;385 if (IsNotANumber() || y.IsNotANumber()) {386 result.flags.set(RealFlag::InvalidArgument);387 result.value = NotANumber();388 } else if (ABS().Compare(y.ABS()) == Relation::Less) {389 return y.HYPOT(*this);390 } else if (IsZero()) {391 return result; // x==y==0392 } else {393 auto yOverX{y.Divide(*this, rounding)}; // y/x394 bool inexact{yOverX.flags.test(RealFlag::Inexact)};395 auto squared{yOverX.value.Multiply(yOverX.value, rounding)}; // (y/x)**2396 inexact |= squared.flags.test(RealFlag::Inexact);397 Real one;398 one.Normalize(false, exponentBias, Fraction::MASKL(1)); // 1.0399 auto sum{squared.value.Add(one, rounding)}; // 1.0 + (y/x)**2400 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 definition412// can be pretty inaccurate when x is much larger than y in magnitude due to413// cancellation. Implement instead with (essentially) arbitrary precision414// long division, discarding the quotient and returning the remainder.415// See runtime/numeric.cpp for more details.416template <typename W, int P>417ValueWithRealFlags<Real<W, P>> Real<W, P>::MOD(418 const Real &p, Rounding rounding) const {419 ValueWithRealFlags<Real> result;420 if (IsNotANumber() || p.IsNotANumber() || IsInfinite()) {421 result.flags.set(RealFlag::InvalidArgument);422 result.value = NotANumber();423 } else if (p.IsZero()) {424 result.flags.set(RealFlag::DivideByZero);425 result.value = NotANumber();426 } else if (p.IsInfinite()) {427 result.value = *this;428 } else {429 result.value = ABS();430 auto pAbs{p.ABS()};431 Real half, adj;432 half.Normalize(false, exponentBias - 1, Fraction::MASKL(1)); // 0.5433 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 defined452// in terms of MOD() with adjustment of the result.453template <typename W, int P>454ValueWithRealFlags<Real<W, P>> Real<W, P>::MODULO(455 const Real &p, Rounding rounding) const {456 ValueWithRealFlags<Real> result{MOD(p, rounding)};457 if (IsNegative() != p.IsNegative()) {458 if (result.value.IsZero()) {459 result.value = result.value.Negate();460 } else {461 result.value =462 result.value.Add(p, rounding).AccumulateFlags(result.flags);463 }464 }465 return result;466}467 468template <typename W, int P>469ValueWithRealFlags<Real<W, P>> Real<W, P>::DIM(470 const Real &y, Rounding rounding) const {471 ValueWithRealFlags<Real> result;472 if (IsNotANumber() || y.IsNotANumber()) {473 result.flags.set(RealFlag::InvalidArgument);474 result.value = NotANumber();475 } else if (Compare(y) == Relation::Greater) {476 result = Subtract(y, rounding);477 } else {478 // result is already zero479 }480 return result;481}482 483template <typename W, int P>484ValueWithRealFlags<Real<W, P>> Real<W, P>::ToWholeNumber(485 common::RoundingMode mode) const {486 ValueWithRealFlags<Real> result{*this};487 if (IsNotANumber()) {488 result.flags.set(RealFlag::InvalidArgument);489 result.value = NotANumber();490 } else if (IsInfinite()) {491 result.flags.set(RealFlag::Overflow);492 } else {493 constexpr int noClipExponent{exponentBias + binaryPrecision - 1};494 if (Exponent() < noClipExponent) {495 Real adjust; // ABS(EPSILON(adjust)) == 0.5496 adjust.Normalize(IsSignBitSet(), noClipExponent, Fraction::MASKL(1));497 // Compute ival=(*this + adjust), losing any fractional bits; keep flags498 result = Add(adjust, Rounding{mode});499 result.flags.reset(RealFlag::Inexact); // result *is* exact500 // Return (ival-adjust) with original sign in case we've generated a zero.501 result.value =502 result.value.Subtract(adjust, Rounding{common::RoundingMode::ToZero})503 .value.SIGN(*this);504 }505 }506 return result;507}508 509template <typename W, int P>510RealFlags Real<W, P>::Normalize(bool negative, int exponent,511 const Fraction &fraction, Rounding rounding, RoundingBits *roundingBits) {512 int lshift{fraction.LEADZ()};513 if (lshift == fraction.bits /* fraction is zero */ &&514 (!roundingBits || roundingBits->empty())) {515 // No fraction, no rounding bits -> +/-0.0516 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 overflow529 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); // Inf534 if constexpr (!isImplicitMSB) {535 word_ = word_.IBSET(significandBits - 1);536 }537 } else {538 // directed rounding: round to largest finite value rather than infinity539 // (x86 does this, not sure whether it's standard behavior)540 word_ = Word{word_.MASKR(word_.bits - 1)};541 if constexpr (isImplicitMSB) {542 word_ = word_.IBCLR(significandBits);543 }544 }545 if (negative) {546 word_ = word_.IBSET(bits - 1);547 }548 RealFlags flags{RealFlag::Overflow};549 if (!fraction.IsZero()) {550 flags.set(RealFlag::Inexact);551 }552 return flags;553 }554 word_ = Word::ConvertUnsigned(fraction).value;555 if (lshift > 0) {556 word_ = word_.SHIFTL(lshift);557 if (roundingBits) {558 for (; lshift > 0; --lshift) {559 if (roundingBits->ShiftLeft()) {560 word_ = word_.IBSET(lshift - 1);561 }562 }563 }564 }565 if constexpr (isImplicitMSB) {566 word_ = word_.IBCLR(significandBits);567 }568 word_ = word_.IOR(Word{exponent}.SHIFTL(significandBits));569 if (negative) {570 word_ = word_.IBSET(bits - 1);571 }572 return {};573}574 575template <typename W, int P>576RealFlags Real<W, P>::Round(577 Rounding rounding, const RoundingBits &bits, bool multiply) {578 int origExponent{Exponent()};579 RealFlags flags;580 bool inexact{!bits.empty()};581 if (inexact) {582 flags.set(RealFlag::Inexact);583 }584 if (origExponent < maxExponent &&585 bits.MustRound(rounding, IsNegative(), word_.BTEST(0) /* is odd */)) {586 typename Fraction::ValueWithCarry sum{587 GetFraction().AddUnsigned(Fraction{}, true)};588 int newExponent{origExponent};589 if (sum.carry) {590 // The fraction was all ones before rounding; sum.value is now zero591 sum.value = sum.value.IBSET(binaryPrecision - 1);592 if (++newExponent >= maxExponent) {593 flags.set(RealFlag::Overflow); // rounded away to an infinity594 }595 }596 flags |= Normalize(IsNegative(), newExponent, sum.value);597 }598 if (inexact && origExponent == 0) {599 // inexact subnormal input: signal Underflow unless in an x86-specific600 // edge case601 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 subnormal607 // inexact multiplication product rounds to a normal result when608 // the guard bit is set or we're not using directed rounding609 } else {610 flags.set(RealFlag::Underflow);611 }612 }613 return flags;614}615 616template <typename W, int P>617void Real<W, P>::NormalizeAndRound(ValueWithRealFlags<Real> &result,618 bool isNegative, int exponent, const Fraction &fraction, Rounding rounding,619 RoundingBits roundingBits, bool multiply) {620 result.flags |= result.value.Normalize(621 isNegative, exponent, fraction, rounding, &roundingBits);622 result.flags |= result.value.Round(rounding, roundingBits, multiply);623}624 625inline enum decimal::FortranRounding MapRoundingMode(626 common::RoundingMode rounding) {627 switch (rounding) {628 case common::RoundingMode::TiesToEven:629 break;630 case common::RoundingMode::ToZero:631 return decimal::RoundToZero;632 case common::RoundingMode::Down:633 return decimal::RoundDown;634 case common::RoundingMode::Up:635 return decimal::RoundUp;636 case common::RoundingMode::TiesAwayFromZero:637 return decimal::RoundCompatible;638 }639 return decimal::RoundNearest; // dodge gcc warning about lack of result640}641 642inline RealFlags MapFlags(decimal::ConversionResultFlags flags) {643 RealFlags result;644 if (flags & decimal::Overflow) {645 result.set(RealFlag::Overflow);646 }647 if (flags & decimal::Inexact) {648 result.set(RealFlag::Inexact);649 }650 if (flags & decimal::Invalid) {651 result.set(RealFlag::InvalidArgument);652 }653 return result;654}655 656template <typename W, int P>657ValueWithRealFlags<Real<W, P>> Real<W, P>::Read(658 const char *&p, Rounding rounding) {659 auto converted{660 decimal::ConvertToBinary<P>(p, MapRoundingMode(rounding.mode))};661 const auto *value{reinterpret_cast<Real<W, P> *>(&converted.binary)};662 return {*value, MapFlags(converted.flags)};663}664 665template <typename W, int P> std::string Real<W, P>::DumpHexadecimal() const {666 if (IsNotANumber()) {667 return "NaN0x"s + word_.Hexadecimal();668 } else if (IsNegative()) {669 return "-"s + Negate().DumpHexadecimal();670 } else if (IsInfinite()) {671 return "Inf"s;672 } else if (IsZero()) {673 return "0.0"s;674 } else {675 Fraction frac{GetFraction()};676 std::string result{"0x"};677 char intPart = '0' + frac.BTEST(frac.bits - 1);678 result += intPart;679 result += '.';680 int trailz{frac.TRAILZ()};681 if (trailz >= frac.bits - 1) {682 result += '0';683 } else {684 int remainingBits{frac.bits - 1 - trailz};685 int wholeNybbles{remainingBits / 4};686 int lostBits{remainingBits - 4 * wholeNybbles};687 if (wholeNybbles > 0) {688 std::string fracHex{frac.SHIFTR(trailz + lostBits)689 .IAND(frac.MASKR(4 * wholeNybbles))690 .Hexadecimal()};691 std::size_t field = wholeNybbles;692 if (fracHex.size() < field) {693 result += std::string(field - fracHex.size(), '0');694 }695 result += fracHex;696 }697 if (lostBits > 0) {698 result += frac.SHIFTR(trailz)699 .IAND(frac.MASKR(lostBits))700 .SHIFTL(4 - lostBits)701 .Hexadecimal();702 }703 }704 result += 'p';705 int exponent = Exponent() - exponentBias;706 if (intPart == '0') {707 exponent += 1;708 }709 result += Integer<32>{exponent}.SignedDecimal();710 return result;711 }712}713 714template <typename W, int P>715llvm::raw_ostream &Real<W, P>::AsFortran(716 llvm::raw_ostream &o, int kind, bool minimal) const {717 if (IsNotANumber()) {718 o << "(0._" << kind << "/0.)";719 } else if (IsInfinite()) {720 if (IsNegative()) {721 o << "(-1._" << kind << "/0.)";722 } else {723 o << "(1._" << kind << "/0.)";724 }725 } else {726 using B = decimal::BinaryFloatingPointNumber<P>;727 B value{word_.template ToUInt<typename B::RawType>()};728 char buffer[common::MaxDecimalConversionDigits(P) +729 EXTRA_DECIMAL_CONVERSION_SPACE];730 decimal::DecimalConversionFlags flags{}; // default: exact representation731 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 753template <typename W, int P>754std::string Real<W, P>::AsFortran(int kind, bool minimal) const {755 std::string result;756 llvm::raw_string_ostream sstream(result);757 AsFortran(sstream, kind, minimal);758 return result;759}760 761// 16.9.180762template <typename W, int P> Real<W, P> Real<W, P>::RRSPACING() const {763 if (IsNotANumber()) {764 return *this;765 } else if (IsInfinite()) {766 return NotANumber();767 } else {768 Real result;769 result.Normalize(false, binaryPrecision + exponentBias - 1, GetFraction());770 return result;771 }772}773 774// 16.9.180775template <typename W, int P> Real<W, P> Real<W, P>::SPACING() const {776 if (IsNotANumber()) {777 return *this;778 } else if (IsInfinite()) {779 return NotANumber();780 } else if (IsZero() || IsSubnormal()) {781 return TINY(); // standard & 100% portable782 } else {783 Real result;784 result.Normalize(false, Exponent(), Fraction::MASKR(1));785 // Can the result be less than TINY()? No, with five commonly786 // used compilers; yes, with two less commonly used ones.787 return result.IsZero() || result.IsSubnormal() ? TINY() : result;788 }789}790 791// 16.9.171792template <typename W, int P>793Real<W, P> Real<W, P>::SET_EXPONENT(std::int64_t expo) const {794 if (IsNotANumber()) {795 return *this;796 } else if (IsInfinite()) {797 return NotANumber();798 } else if (IsZero()) {799 return *this;800 } else {801 return SCALE(Integer<64>(expo - UnbiasedExponent() - 1)).value;802 }803}804 805// 16.9.171806template <typename W, int P> Real<W, P> Real<W, P>::FRACTION() const {807 return SET_EXPONENT(0);808}809 810template class Real<Integer<16>, 11>;811template class Real<Integer<16>, 8>;812template class Real<Integer<32>, 24>;813template class Real<Integer<64>, 53>;814template class Real<X87IntegerContainer, 64>;815template class Real<Integer<128>, 113>;816} // namespace Fortran::evaluate::value817