Compute coeffs using Rationals

Project.toml

@@ -1,6 +1,6 @@
 name = "FiniteDifferences"
 uuid = "26cc04aa-876d-5657-8c51-4c34ba976000"
-version = "0.9.3"
+version = "0.9.4"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/methods.jl

@@ -157,10 +157,12 @@ end
 
 # Compute coefficients for the method
 function _coefs(grid::AbstractVector{<:Real}, p::Integer, q::Integer)
-    C = [g^i for i in 0:(p - 1), g in grid]
-    x = zeros(Int, p)
+    # For high precision on the \ we use Rational, and to prevent overflows we use Int128
+    # At the end we go to Float64 for fast floating point math (rather than rational math)
+    C = [Rational{Int128}(g^i) for i in 0:(p - 1), g in grid]
+    x = zeros(Rational{Int128}, p)
     x[q + 1] = factorial(q)
-    return C \ x
+    return Float64.(C \ x)
 end
 
 # Estimate the bound on the function value and its derivatives at a point

test/methods.jl

@@ -37,6 +37,17 @@ using FiniteDifferences: Forward, Backward, Central, Nonstandard
         @test central_fdm(5, 1)(abs, 0.001) ≈ 1.0
     end
 
+    @testset "Accuracy at high orders, with high adapt" begin
+        # Regression test against issues with precision during computation of _coeffs
+        # see https://github.com/JuliaDiff/FiniteDifferences.jl/issues/64
+
+        @test fdm(central_fdm(9, 5), exp, 1.0, adapt=4) ≈ exp(1) atol=1e-7
+
+        poly(x) = 4x^3 + 3x^2 + 2x + 1
+        @test fdm(central_fdm(9, 3), poly, 1.0, adapt=4) ≈ 24 atol=1e-11
+    end
+
+
     @testset "Printing FiniteDifferenceMethods" begin
         @test sprint(show, central_fdm(2, 1)) == """
             FiniteDifferenceMethod: