> Do you know a way to get the IEEE 754 rounding behavior without invoking C undefined behavior?

One option is to hard-code the actual boundary, which is 2**128 * (1 - 2**-25) (as opposed to FLT_MAX, which is 2**128 * (1 - 2**-24)): values equal to or larger than 2**128 * (1 - 2**-25) in absolute value should raise. But that means assuming IEEE 754 and round-ties-to-even, which isn't an outrageous assumption but does make the solution feel a bit fragile.

An alternative would be to scale values in the range (FLT_MAX, 2.0 * FLT_MAX] by 0.5 before doing the conversion, something like this:

        if (fabs(x) > FLT_MAX && !Py_IS_INFINITY(x)) {
          double half_x = 0.5 * x;
          if (half_x > FLT_MAX) {
            goto Overflow;
          }
          float half_y = (float)half_x;
          if (half_y > 0.5 * FLT_MAX) {
            goto Overflow;
          }
        }