Julia 1.8.0: 10.0^-Inf = NaN

General Usage
1

Just switched to Julia 1.8 and was surprised by this:

julia> 10.0^-Inf
NaN

However:

julia> exp10(-Inf)
0.0

I am pretty sure that 10.0^-Inf returned 0.0 in 1.6 and 1.7

5 Likes
2

My fault. Good catch!.

9 Likes
3

The really annoying thing is that I wrote a really good test set for this because I knew there were lots of potential problems here.

for T in (Float16, Float32, Float64)
    for x in (0.0, -0.0, 1.0, 10.0, 2.0, Inf, NaN, -Inf, -NaN)
        for y in (0.0, -0.0, 1.0, -3.0,-10.0 , Inf, NaN, -Inf, -NaN)
            got, expected = T(x)^T(y), T(big(x))^T(y)
            if isnan(expected)
                @test isnan_type(T, got) || (x,y)
            else
                @test got === expected || (x,y)
            end
        end
    end
end

Unfortunately, I wrote got, expected = T(x)^T(y), T(big(x))^T(y) instead of got, expected = T(x)^T(y), T(big(x)^T(y)) so the tests was testing absolutely nothing.

13 Likes
4

PR to fix it on master Fix some pow edge cases by oscardssmith · Pull Request #46412 · JuliaLang/julia · GitHub (I’ll need to manually backport to 1.8)

6 Likes
5

Ouch, we need tests for the tests (and tests for the tests-tests…).

FYI: For a hotfix if you can’t wait for 1.8.1 (and don’t want to downgrade):

julia> import Base:^, @constprop

julia> @constprop :aggressive function ^(x::Float64, y::Float64)
           xu = reinterpret(UInt64, x)
           xu == reinterpret(UInt64, 1.0) && return 1.0
           # Exponents greater than this will always overflow or underflow.
           # Note that NaN can pass through this, but that will end up fine.
           if !(abs(y)<0x1.8p62)
               isnan(y) && return y
               y = sign(y)*0x1.8p62
           end
           yint = unsafe_trunc(Int, y) # This is actually safe since julia freezes the result
           y == yint && return @noinline x^yint
           2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x==0
           x<0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer
           !isfinite(x) && return x*(y>0 || isnan(x))           # x is inf or NaN
           if xu < (UInt64(1)<<52) # x is subnormal
               xu = reinterpret(UInt64, x * 0x1p52) # normalize x
               xu &= ~sign_mask(Float64)
               xu -= UInt64(52) << 52 # mess with the exponent
           end
           return pow_body(xu, y)
       end
^ (generic function with 68 methods)

julia> -10.0^(-Inf)
-0.0

# didn't even need from the PR, for Float64:

julia> @inline function pow_body(x::T, y::T) where T <: Union{Float16, Float32}
           return T(exp2(log2(abs(widen(x))) * y))
       end

# would need to copy too from PR, for the other, less used, types:
@constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32}
[..]
    return pow_body(x, y)
end

Now, it’s not a very serious bug (or at least “poisoning” with NaN so you will know, not get wrpong non-NaN answers), but it got me thinking, could we have Bugfix.jl with such fixes, that people could opt into? For those who don’t (know of), we still need 1.8.1, and it would mean extra work to fix in two places. Could then Bugfix.jl be a stdlib of Julia, meaning 1.8.1 would just have to reference a certain version. I’m not sure if users of an older version of Julia can upgrade stdlibs.

6

This isn’t a valid fix for 1.8 (the _log_ext implimentation changed). The correct version is

import Base: ^
import Base.Math: pow_body, throw_exp_domainerror
@constprop :aggressive function ^(x::Float64, y::Float64)
    x == 1 && return 1.0
    # Exponents greater than this will always overflow or underflow.
    # Note that NaN can pass through this, but that will end up fine.
    if !(abs(y)<0x1.8p62)
        isnan(y) && return y
        y = sign(y)*0x1.8p62
    end
    yint = unsafe_trunc(Int, y) # This is actually safe since julia freezes the result
    y == yint && return x^yint
    x<0 && throw_exp_domainerror(x)
    !isfinite(x) && return x*(y>0 || isnan(x))
    x==0 && return abs(y)*Inf*(!(y>0))
    return pow_body(x, y)
end
4 Likes
7

Nah, that’s where you want an alert reviewer.

2 Likes