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 809