Portable way to check whether a floating point division would end in +-inf

We could search for the limit using trial and error:

#include <iostream>
#include <limits>

#include <cmath>

int main() {
    float limit = 0.0f;
    float result = 1.0f / limit;
    while (
        result == std::numeric_limits<float>::infinity()
        or std::isnan(result)
    ) {
        limit = std::nextafter(limit, 1.0f);
        result = 1.0f / limit;
    }
    std::cout << "Limit = " << limit << std::endl;
    std::cout << "1.0f / Limit = " << 1.0f / limit << std::endl;
}

This outputs on my system:

Limit = 2.93874e-39
1.0f / Limit = 3.40282e+38

This is however, not a very efficient solution. If we could make this algorithm constexpr, this would mitigate the issue but unfortunately std::nextafter() isn't constexpr.

If you know your environment is using IEEE-754 airthmetic, these limits are probably constant, but as you ask for portability, we can't always assume that.

Prerequisites

C++ does not mandate IEEE-754 or a particular rounding method. For this answer, I will assume IEEE-754 is used with a binary format and round-to-nearest, ties-to-even.

Conclusion

1/x overflows iff fabs(x) <= std::ldexp(1, -std::numeric_limits<float>::max_exponent). For a constant expression, you can use std::numeric_limits<float>::min()/4.

Discussion

Rounding at the end of the finite range is performed as if the exponent kept going. E.g., using decimal for illustration, if the highest representable finite number were 9.99•10¹⁷, the next representable number, if the exponent were not limited, would be 1.00•10¹⁸. The midpoint between those two is 9.995•10¹⁷, so numbers under that are rounded down and numbers above that are rounded up. With ties-to-even, 9.995•10¹⁷ is rounded up.

For a binary format, the greatest representable value is (2−ε)•2^q, where ε is the “machine epsilon” (the ULP of 1, so 2-ε is the greatest representable significand) and q is the maximum exponent. Then the point where rounding occurs is (2−½ε)•2^q.

If 1/x < (2−½ε)•2^q, the result is rounded downward. Otherwise, it is rounded upward, to ∞. Thus, the result is less than ∞ iff x > 1/((2−½ε)•2^q) = 2^−q/(2-½ε).

1/(2-½ε) is slightly greater than ½, by less than ½ε, so the greatest representable value less than or equal to it is ½. Thus, the result of 1/x is less than ∞ iff x > 2^−q/2 = 2^−q−1.

C++ tells us the maximum exponent with std::numeric_limits<double>::max_exponent (defined in the header <limits>). However, C++ calibrates this exponent for a significand range of [½, 1) instead of IEEE-754’s [1, 2), so it is one greater than q. Thus the −q−1 we want is simply -std::numeric_limits<double>::max_exponent.

We can calculate 2^−q−1 with the ldexp function (declared in <cmath>): std::ldexp(1, -std::numeric_limits<float>::max_exponent).

With Apple Clang 11, this program:

#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>


int main(void)
{
    float x = std::ldexp(1, -std::numeric_limits<float>::max_exponent);

    std::cout << std::setprecision(20) << x << " is too small, result will overflow:\n";
    std::cout << "\t" << 1/x << ".\n";

    x = std::nexttoward(x, INFINITY);

    std::cout << std::setprecision(20) << x << " is just big enough, result will not overflow:\n";
    std::cout << "\t" << 1/x << ".\n";
}

produces:

2.9387358770557187699e-39 is too small, result will overflow:
    inf.
2.9387372783541830947e-39 is just big enough, result will not overflow:
    3.4028220466166163425e+38.

Accounting for negative numbers as well, 1/x overflows iff fabs(x) <= std::ldexp(1, -std::numeric_limits<float>::max_exponent).

Because of the way IEEE-754 specifies the exponent range, std::ldexp(1, -std::numeric_limits<float>::max_exponent) equals std::numeric_limits<float>::min()/4. (IEEE-754 specifies the minimum normal exponent to be 1−q, so the −q−1 we desire is (1-q)-2.)

Since you're looking for constant values, we can actually use an SMT solver to find the minimum/maximum values of x where the division 1/x will produce infinity. I did this using Microsoft's z3 SMT solver (https://github.com/Z3Prover/z3), scripting it from Haskell SBV bindings (http://leventerkok.github.io/sbv/). I get:

Prelude Data.SBV> optimize Lexicographic $ do x <- sFloat "x"; constrain (fpIsInfinite (1/x)); minimize "x" x
Optimal model:
  x   = -2.938736e-39 :: Float
  x_0 =    2145386495 :: Word32
Prelude Data.SBV> optimize Lexicographic $ do x <- sFloat "x"; constrain (fpIsInfinite (1/x)); maximize "x" x
Optimal model:
  x   = 2.938736e-39 :: Float
  x_0 =   2149580800 :: Word32

If you squint at this, you'll find that the values it's suggesting are between -2.938736e-39 and 2.938736e-39 where 1/x becomes infinity.

If you want to write this "portably" without any rounding issues, you should use the hexadecimal notation, i.e., -0x1p-128 and 0x1p-128.

I believe these numbers match @Eric Postpischill's values; and his analysis is very useful of course, but this is another way of finding such values using automated theorem-proving techniques.