Some fastpow limits differ from ^

[andreasnoack]

julia> for x in (0.0, Inf), y in (0.0, Inf)
         @eval @show $x^$y
         @eval @show DiffEqBase.fastpow($x, $y)
         println()
       end
0.0 ^ 0.0 = 1.0
DiffEqBase.fastpow(0.0, 0.0) = 1.0f0

0.0 ^ Inf = 0.0
DiffEqBase.fastpow(0.0, Inf) = NaN32

Inf ^ 0.0 = 1.0
DiffEqBase.fastpow(Inf, 0.0) = 1.0f0

Inf ^ Inf = Inf
DiffEqBase.fastpow(Inf, Inf) = NaN32

I'm wondering if this could be causing https://github.com/PumasAI/Pumas.jl/issues/761

[ChrisRackauckas]

0.0 ^ Inf = 0.0
DiffEqBase.fastpow(0.0, Inf) = NaN32

This one would. It would make almost no error in the integration be a NaN. The linear ODE be a case where methods like Rosenbrock could solve it with zero error and throw a NaN, which would then cause a NaN dt and exit. @YingboMa

[andreasnoack]

Below is a more complete list and there are a couple of other places things could go wrong. I guess x and y might be non-negative in this application so maybe many of the cases are irrelevant but based on your comment I'd suspect that

julia> DiffEqBase.fastpow(1.0, Inf)
NaN32

could be a problem.

julia> vals = (-Inf, -1.0, -0.0, 0.0, 1.0, Inf, NaN)

julia> for x in vals, y in vals
         @eval @show $x^$y
         @eval @show DiffEqBase.fastpow($x, $y)
         println()
       end
(-Inf) ^ -Inf = 0.0
DiffEqBase.fastpow(-Inf, -Inf) = NaN32

(-Inf) ^ -1.0 = -0.0
DiffEqBase.fastpow(-Inf, -1.0) = 2.938736f-39

(-Inf) ^ -0.0 = 1.0
DiffEqBase.fastpow(-Inf, -0.0) = 1.0f0

(-Inf) ^ 0.0 = 1.0
DiffEqBase.fastpow(-Inf, 0.0) = 1.0f0

(-Inf) ^ 1.0 = -Inf
DiffEqBase.fastpow(-Inf, 1.0) = Inf32

(-Inf) ^ Inf = Inf
DiffEqBase.fastpow(-Inf, Inf) = NaN32

(-Inf) ^ NaN = NaN
DiffEqBase.fastpow(-Inf, NaN) = NaN32

(-1.0) ^ -Inf = 1.0
DiffEqBase.fastpow(-1.0, -Inf) = NaN32

(-1.0) ^ -1.0 = -1.0
DiffEqBase.fastpow(-1.0, -1.0) = 1.0f0

(-1.0) ^ -0.0 = 1.0
DiffEqBase.fastpow(-1.0, -0.0) = 1.0f0

(-1.0) ^ 0.0 = 1.0
DiffEqBase.fastpow(-1.0, 0.0) = 1.0f0

(-1.0) ^ 1.0 = -1.0
DiffEqBase.fastpow(-1.0, 1.0) = 1.0f0

(-1.0) ^ Inf = 1.0
DiffEqBase.fastpow(-1.0, Inf) = NaN32

(-1.0) ^ NaN = NaN
DiffEqBase.fastpow(-1.0, NaN) = NaN32

-0.0 ^ -Inf = Inf
DiffEqBase.fastpow(-0.0, -Inf) = NaN32

-0.0 ^ -1.0 = -Inf
DiffEqBase.fastpow(-0.0, -1.0) = 1.7014118f38

-0.0 ^ -0.0 = 1.0
DiffEqBase.fastpow(-0.0, -0.0) = 1.0f0

-0.0 ^ 0.0 = 1.0
DiffEqBase.fastpow(-0.0, 0.0) = 1.0f0

-0.0 ^ 1.0 = -0.0
DiffEqBase.fastpow(-0.0, 1.0) = 5.877472f-39

-0.0 ^ Inf = 0.0
DiffEqBase.fastpow(-0.0, Inf) = NaN32

-0.0 ^ NaN = NaN
DiffEqBase.fastpow(-0.0, NaN) = NaN32

0.0 ^ -Inf = Inf
DiffEqBase.fastpow(0.0, -Inf) = NaN32

0.0 ^ -1.0 = Inf
DiffEqBase.fastpow(0.0, -1.0) = 1.7014118f38

0.0 ^ -0.0 = 1.0
DiffEqBase.fastpow(0.0, -0.0) = 1.0f0

0.0 ^ 0.0 = 1.0
DiffEqBase.fastpow(0.0, 0.0) = 1.0f0

0.0 ^ 1.0 = 0.0
DiffEqBase.fastpow(0.0, 1.0) = 5.877472f-39

0.0 ^ Inf = 0.0
DiffEqBase.fastpow(0.0, Inf) = NaN32

0.0 ^ NaN = NaN
DiffEqBase.fastpow(0.0, NaN) = NaN32

1.0 ^ -Inf = 1.0
DiffEqBase.fastpow(1.0, -Inf) = NaN32

1.0 ^ -1.0 = 1.0
DiffEqBase.fastpow(1.0, -1.0) = 1.0f0

1.0 ^ -0.0 = 1.0
DiffEqBase.fastpow(1.0, -0.0) = 1.0f0

1.0 ^ 0.0 = 1.0
DiffEqBase.fastpow(1.0, 0.0) = 1.0f0

1.0 ^ 1.0 = 1.0
DiffEqBase.fastpow(1.0, 1.0) = 1.0f0

1.0 ^ Inf = 1.0
DiffEqBase.fastpow(1.0, Inf) = NaN32

1.0 ^ NaN = 1.0
DiffEqBase.fastpow(1.0, NaN) = NaN32

Inf ^ -Inf = 0.0
DiffEqBase.fastpow(Inf, -Inf) = NaN32

Inf ^ -1.0 = 0.0
DiffEqBase.fastpow(Inf, -1.0) = 2.938736f-39

Inf ^ -0.0 = 1.0
DiffEqBase.fastpow(Inf, -0.0) = 1.0f0

Inf ^ 0.0 = 1.0
DiffEqBase.fastpow(Inf, 0.0) = 1.0f0

Inf ^ 1.0 = Inf
DiffEqBase.fastpow(Inf, 1.0) = Inf32

Inf ^ Inf = Inf
DiffEqBase.fastpow(Inf, Inf) = NaN32

Inf ^ NaN = NaN
DiffEqBase.fastpow(Inf, NaN) = NaN32

NaN ^ -Inf = NaN
DiffEqBase.fastpow(NaN, -Inf) = NaN32

NaN ^ -1.0 = NaN
DiffEqBase.fastpow(NaN, -1.0) = 1.958973f-39

NaN ^ -0.0 = 1.0
DiffEqBase.fastpow(NaN, -0.0) = 1.0f0

NaN ^ 0.0 = 1.0
DiffEqBase.fastpow(NaN, 0.0) = 1.0f0

NaN ^ 1.0 = NaN
DiffEqBase.fastpow(NaN, 1.0) = Inf32

NaN ^ Inf = NaN
DiffEqBase.fastpow(NaN, Inf) = NaN32

NaN ^ NaN = NaN
DiffEqBase.fastpow(NaN, NaN) = NaN32

[ChrisRackauckas]

The power part is a fixed (positive) finite exponent and the other part is positive as well. Given that, the only thing that sticks out to me is

0.0 ^ 1.0 = 0.0
DiffEqBase.fastpow(0.0, 1.0) = 5.877472f-39

but that shouldn't be a major issue. I'm not 100% convinced the problem is fastpow anymore, but rather something to do with stiffness detection and switching.

[oscardssmith]

This can probably be closed at this point.