BigFloats give Float64 accuracy — doesn't use the precision of the arguments

[stevengj]

This is related to #98 — Float64 precision seems to be baked into the package, whereas it would be more flexible (and more Julian) to use the precision of the arguments. For example, using BigFloat in the following example only gives 16 digits of accuracy:

julia> using FiniteDifferences

julia> extrapolate_fdm(central_fdm(2, 1), sin, big"1.0")[1] - cos(big"1.0")
-6.71174699531887290713204926350530055924229547805016487900793424335727248883486e-17

In contrast, "manually" calling Richardson extrapolation with a 2nd-order finite-difference approximation gives 70 digits (in 10 iterations):

julia> using Richardson

julia> dsin,_ = extrapolate(big"0.1", rtol=0, power=2) do h
           @show Float64(h)
           (sin(1+h) - sin(1-h)) / 2h
       end
Float64(h) = 0.1
Float64(h) = 0.0125
Float64(h) = 0.0015625
Float64(h) = 0.0001953125
Float64(h) = 2.44140625e-5
Float64(h) = 3.0517578125e-6
Float64(h) = 3.814697265625e-7
Float64(h) = 4.76837158203125e-8
Float64(h) = 5.960464477539063e-9
Float64(h) = 7.450580596923829e-10
(0.5403023058681397174009366074429766037323104206179222276700972553811007395485809, 3.557003874037244232691736616015521320518299017661431723628543167235436096570656e-70)

julia> dsin - cos(big"1.0")
3.44774105579695699101361233110829910923644073773066771188994182638266121185316e-70

So, someplace in FiniteDifferences is either hard-coding a tolerance (rather than using eps(float(x))), or "contaminating" the calculation with an inexact Float64 literal.

cc @hammy4815

[oxinabox]

This is true.
It is hard coded to use Float64 in a bunch of places, like when storing the precomputed grids.
It would take a fair bit to fix, but would be good if someone was interested in doing so.
I don't have the need or the time right now

[sinhtrung]

I wrote the following code to explore this issue,

using FiniteDifferences, Cthulhu
fdm = central_fdm(5, 1);
@descend fdm(sin, big"1.0")

Delving a bit more into the code with Cthulhu, I obtained this,

_limit_step(m::FiniteDifferenceMethod, x::T, step::Real, acc::Real) where T<:AbstractFloat @ FiniteDifferences ~/.julia/dev/FiniteDifferences/src/methods.jl:408
408 function (_limit_step(
409     m::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}}::FiniteDifferenceMethod, x::BigFloat::T, step::BigFloat::Real, acc::BigFloat::Real,
410 ) where {T<:AbstractFloat})::Tuple{Union{Float64, BigFloat}, Union{Float64, BigFloat}}
411     # First, limit the step size based on the maximum range.
412     step_max::Float64 = (m::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}}.max_range::Float64 / maximum(abs::Core.Const(abs).(m::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}}.grid::StaticArraysCore.SVector{5, Int64})::StaticArraysCore.SVector{5, Int64})::Int64)::Float64
413     if (step::BigFloat > step_max::Float64)::Bool
414         step::Float64 = step_max::Float64
415         acc = NaN
416     end
417     # Second, prevent very large step sizes, which can occur for high-order methods or
418     # slowly-varying functions.
419     step_default::BigFloat, _ = _compute_step_acc_default(m::FiniteDifferences.AdaptedFiniteDifferenceMethod{5, 1, FiniteDifferences.UnadaptedFiniteDifferenceMethod{7, 5}}, x::BigFloat)::Core.Const(2)
420     step_max_default::BigFloat = 1000step_default::BigFloat::BigFloat
421     if (step::Union{Float64, BigFloat} > step_max_default::BigFloat)::Bool
422         step::BigFloat = step_max_default::BigFloat
423         acc = NaN
424     end
425     return step::Union{Float64, BigFloat}, acc::Union{Float64, BigFloat}::Tuple{Union{Float64, BigFloat}, Union{Float64, BigFloat}}
426 end

From what I can understand in this Cthulhu output, on the function return of _limit_step both step and acc had unions with Float64 and BigFloat step::Union{Float64, BigFloat}, acc::Union{Float64, BigFloat}::Tuple{Union{Float64, BigFloat}. However, on the original call step and acc had type signatures without unions step::BigFloat::Real, acc::BigFloat::Real.