[flang] Fix ICE for sqrt(0.0) evaluation

[FruitClover]

Summary

During real range reduction to [0.5, 4) with

       SQRT(2**(2a) * x) = SQRT(2**(2a)) * SQRT(x) = 2**a * SQRT(x)

we fall into inf. recursion if IsZero() == true.

Explicitly handle SQRT(0.0) instead of additional checks during folding. Also
add helpers for +0.0/-0.0 generation to clean up a bit.

[FruitClover]

Fix typo

[klausler]

What was the "ICE"?

flang/include/flang/Evaluate/real.h

201

MASKL(1) is clearer, since you're moving things around and can clean them up a little.

[FruitClover]

Replaced IBSET -> MASKL, updated summary details.

flang/include/flang/Evaluate/real.h

201

Updated, thanks.

flang/include/flang/Evaluate/real.h

public:
   // TODO: Configurable NotANumber representations
   static constexpr Real NotANumber() {
     return {Word{maxExponent}
                 .SHIFTL(significandBits)
                 .IBSET(significandBits - 1)
                 .IBSET(significandBits - 2)};
   }
 
+  static constexpr Real PositiveZero() { return Real{}; }
+
+  static constexpr Real NegativeZero() { return {Word{}.MASKL(1)}; }
+
   static constexpr Real Infinity(bool negative) {
     Word infinity{maxExponent};
     infinity = infinity.SHIFTL(significandBits);
     if (negative) {
       infinity = infinity.IBSET(infinity.bits - 1);
     }
     return {infinity};
   }

flang/lib/Evaluate/real.cpp

if (exponent == yExponent && isNegative != yIsNegative) {
     Ordering order{GetSignificand().CompareUnsigned(y.GetSignificand())};
     if (order == Ordering::Less) {
       // Same exponent, opposite signs, and y is larger in magnitude
       return y.Add(*this, rounding);
     }
     if (order == Ordering::Equal) {
       // x + (-x) -> +0.0 unless rounding is directed downwards
       if (rounding.mode == common::RoundingMode::Down) {
-        result.value.word_ = result.value.word_.IBSET(bits - 1); // -0.0
+        result.value = NegativeZero();
       }
       return result;
     }
   }
   // Our exponent is greater than y's, or the exponents match and y is not
   // of the opposite sign and greater magnitude.  So (x+y) will have the
   // same sign as x.
   Fraction fraction{GetFraction()};
if (IsInfinite()) {
         result.value = NotANumber();
         result.flags.set(RealFlag::InvalidArgument);
       } else { // x/0 -> Inf, Inf/0 -> Inf
         result.value = Infinity(isNegative);
         result.flags.set(RealFlag::DivideByZero);
       }
     } else if (IsZero() || y.IsInfinite()) { // 0/x, x/Inf -> 0
       if (isNegative) {
-        result.value.word_ = result.value.word_.IBSET(bits - 1);
+        result.value = NegativeZero();
       }
     } else {
       // dividend and divisor are both finite and nonzero numbers
       Fraction top{GetFraction()}, divisor{y.GetFraction()};
       std::int64_t exponent{CombineExponents(y, true)};
       Fraction quotient;
       bool msb{false};
       if (!top.BTEST(top.bits - 1) || !divisor.BTEST(divisor.bits - 1)) {
ValueWithRealFlags<Real<W, P>> Real<W, P>::SQRT(Rounding rounding) const {
   if (IsNotANumber()) {
     result.value = NotANumber();
     if (IsSignalingNaN()) {
       result.flags.set(RealFlag::InvalidArgument);
     }
   } else if (IsNegative()) {
     if (IsZero()) {
       // SQRT(-0) == -0 in IEEE-754.
-      result.value.word_ = result.value.word_.IBSET(bits - 1);
+      result.value = NegativeZero();
     } else {
       result.value = NotANumber();
     }
   } else if (IsInfinite()) {
     // SQRT(+Inf) == +Inf
     result.value = Infinity(false);
+  } else if (IsZero()) {
+    result.value = PositiveZero();
   } else {
     int expo{UnbiasedExponent()};
     if (expo < -1 || expo > 1) {
       // Reduce the range to [0.5 .. 4.0) by dividing by an integral power
       // of four to avoid trouble with very large and very small values
       // (esp. truncation of subnormals).
       // SQRT(2**(2a) * x) = SQRT(2**(2a)) * SQRT(x) = 2**a * SQRT(x)
       Real scaled;

flang/test/Evaluate/folding28.f90

module m
   real(4), parameter :: sqrt_under1 = sqrt(.96875)
   logical, parameter :: test_sqrt_under1 = sqrt_under1 == .984250962734222412109375
   ! oddball case: the value before 1. is also its own sqrt, but not its own square
   real(4), parameter :: before_1 = z'3f7fffff' ! .999999940395355224609375
   real(4), parameter :: sqrt_before_1 = sqrt(before_1)
   logical, parameter :: test_before_1 = sqrt_before_1 == before_1
   real(4), parameter :: sq_sqrt_before_1 = sqrt_before_1 * sqrt_before_1
   logical, parameter :: test_sq_before_1 = sq_sqrt_before_1 < before_1
+  ! ICE at 0.0
+  real(4), parameter :: sqrt_zero_4 = sqrt(0.0)
+  logical, parameter :: test_sqrt_zero_4 = sqrt_zero_4 == 0.0
+  real(8), parameter :: sqrt_zero_8 = sqrt(0.0)
+  logical, parameter :: test_sqrt_zero_8 = sqrt_zero_8 == 0.0
 end module