faster isfinite

[oscardssmith]

julia> function mytest!(f)
           s=0
           for x in typemin(UInt32):typemax(UInt32)
               s += f(reinterpret(Float32, x)) 
           end
           s
       end
mytest! (generic function with 1 method)

#old
julia> @btime mytest!(isfinite)
  1.971 s (0 allocations: 0 bytes)
4278190080

#new
julia> @btime mytest!(isfinite)
  1.770 s (0 allocations: 0 bytes)
4278190080

This has been tested exhaustively for Float32 and Float16

[mikmoore]

Nice!

It seems brave to assume that x-x === zero(x) for any finite AbstractFloat, since we really have no documented interface or conditions for AbstractFloat. This method would be unsuitable for any mutable (including Base's very own BigFloat, which thankfully has its own method) or some not-completely-unreasonable type where x-x === -zero(x) for some values.

I think this method should apply only to IEEEFloat. The AbstractFloat method should remain intact, become even more generic (e.g., isfinite(x::AbstractFloat) = !(isnan(x) | isinf(x))), or be removed entirely and fall back to Real.

[giordano]

As I said in #46163 (comment), I feel like this should be iszero(x - x), and also iszero can be optimised for IEEFloat with x === zero(x), instead of x == zero(x)

[oscardssmith]

julia> iszero(-0.0)
true

[oscardssmith]

so that would break the optimization. The reason for this change is that === on floating points can turn into a bit-compare rather than a more expensive floating point compare.

[giordano]

Right, sad 😕

[vtjnash]

Instead then:

julia> Base.iszero(x::Float64) = reinterpret(UInt64, x) & 0x7fff_ffff_ffff_ffff === UInt(0)

(btime cannot distinguish these three approaches in performance, and returns 3.5 ns for all of them)

[mikmoore]

I'll define the initially-proposed solution
isfinite_1(x::IEEEFloat) = x-x===zero(x).
Another contender, related to the iszero(x-x) variant but using the other half of the fact that x-x can only take the values +0.0 or NaN, is
isfinite_2(x::IEEEFloat) = !isnan(x-x).
Of course, now that we're limited to IEEEFloat we can instead reach for bit-twiddling
isfinite_3(x::IEEEFloat) = (reinterpret(Unsigned,x) & Base.exponent_mask(typeof(x))) != Base.exponent_mask(typeof(x))

julia> code_native(isfinite_1,(Float64,);debuginfo=:none)
        vsubsd  %xmm0, %xmm0, %xmm0
        vmovq   %xmm0, %rax
        testq   %rax, %rax
        sete    %al
        retq

julia> code_native(isfinite_2,(Float64,);debuginfo=:none)
        vsubsd  %xmm0, %xmm0, %xmm0
        vucomisd        %xmm0, %xmm0
        setnp   %al
        retq

julia> code_native(isfinite_3,(Float64,);debuginfo=:none)
        vmovq   %xmm0, %rax
        movabsq $9218868437227405312, %rcx      # imm = 0x7FF0000000000000
        andnq   %rcx, %rax, %rax
        setne   %al
        retq

Compared to isfinite_1, isfinite_2 is one instruction shorter in isolation (on my x86). isfinite_3 is the same number of instructions as isfinite_1 but one of those is a hoistable movabsq. Another advantage is that it requires no floating point operations. Note that, after inlining, this may not be what these functions look like in the wild.

Keep in mind that mytest! is not a very canonical use of isfinite. Applying a @code_native shows this fact. That said, I see the _1 and _2 variants benchmarking identically and _3 about 15% faster.

EDIT:
Actually, I'm increasingly a fan of the !isnan(x-x) variant. While it's not quite as fast as bit twiddling in this nanobenchmark, I think that it has the benefit of clarity. Further, I think that it /would/ be a suitable ::AbstractFloat definition.

[nlw0]

How about this?

isfinite_4(x::IEEEFloat) = x + one(x) > x

[oscardssmith]

I would think that wouldn't be better since > on floating point isn't completely trivial. It's also wrong since 2.0^53 +1.0 == 2.0^53

[nlw0]

Well observed... The point was just that it this handles the -0.0 case, and naturally fails for NaN or Inf.... Maybe there could be a simple way to generate a valid successor, but that would still fail for whatever is the highest valid number? oh well