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

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

My fault. Good catch!.

9 Likes

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

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

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.

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

Nah, thatâ€™s where you want an alert reviewer.

2 Likes