Julia: Fallback math methods for ::Real are likely to call themselves

[tkoolen]

On 0.6.2:

julia> using ForwardDiff

julia> ForwardDiff.derivative(sinpi, 1.)
0.0

on b571283a4c:

julia> using ForwardDiff

julia> ForwardDiff.derivative(sinpi, 1.)
ERROR: StackOverflowError:
Stacktrace:
 [1] sinpi(::ForwardDiff.Dual{ForwardDiff.Tag{typeof(sinpi),Float64},Float64,1}) at ./special/trig.jl:865 (repeats 80000 times)

This is because of https://github.com/JuliaLang/julia/blob/7674acada0c6b645319c54ed2a6606b5408f5693/base/special/trig.jl#L865

which seems unfortunate. Similarly, log1p required a workaround in DiffRules. I wholeheartedly agree with https://github.com/JuliaDiff/ForwardDiff.jl/issues/261#issuecomment-330308814:

Looks like Base.Math.JuliaLibm.log1p(x::Real) just calls Base.Math.JuliaLibm.log1p(float(x)), which basically guarantees infinite recursion in the general case

log and cospi also have this issue. Maybe others as well?

[StefanKarpinski]

Do you have a different way to express this in mind that wouldn't have this issue?

[tkoolen]

How about turning

sinpi(x::T) where T<:AbstractFloat

into

_sinpi(x::T) where T<:AbstractFloat

and have

sinpi(x::Real) = _sinpi(float(x)) 

[mbauman]

This is similarly handled in other places with a generic AbstractFloat specialization.

[tkoolen]

That assumes that float always returns an AbstractFloat. ForwardDiff.Dual doesn't adhere to that, and it probably shouldn't be an AbstractFloat, since the backing type can also be non-AbstractFloat.

[martinholters]

We could make that assumption more explicit by using convert(AbstractFloat, x) instead of float (x).

[tkoolen]

That would certainly be an improvement, MethodError is better than StackOverflowError. But the sinpi/log1p/... changes from 0.6.2 to 0.7 would still be breaking, and it's not like the old implementation returned wrong results (as far as I understand). Restoring the old behavior using something like https://github.com/JuliaLang/julia/issues/26552#issuecomment-374958654 seems to me like it would save a lot of hassle for package authors.

[chriscoey]

any path forward on this? I am struggling with this issue: https://github.com/JuliaDiff/ForwardDiff.jl/issues/324

[tpapp]

I think that the original intention for Base.float was to provide a fallback to already defined special functions. The implicit assumption in various instances of

my_calculation(x::Real) = my_calculation(float(x))

is that float returns something my_calculation can deal with.

With the combination of user-defined types and functions, this is very difficult to enforce or even explicitly specify. Eg one could require that float(x)::AbstractFloat, but as a package can just define a new <:AbstractFloat type, and a my_calculation in Base or another package may not be equipped to deal with that.

In retrospect, I think

1. float should have been named Base.basefloat,
2. with an explicit requirement to return only BigFloat, Float16, Float32, Float64 (ie subtypes(AbstractFloat) in Base,
3. perhaps not exported.

(1) and (3) are breaking suggestions, but (2) can be addressed in the docstring of float for now.

Packages should not define a method for float unless they can comply with (2). It is tempting to expect an automatic fallback from it, but it almost always violates implicit assumptions and leads to bugs like the above.