diff --git a/src/methods.jl b/src/methods.jl
index 20b56ea9..511ec95e 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -26,6 +26,31 @@ function central_grid(p::Int)
     end
 end
 
+
+# Check the method and derivative orders for consistency
+function _check_p_q(p::Integer, q::Integer)
+    q >= 0 || throw(ArgumentError("order of derivative must be nonnegative"))
+    q < p || throw(ArgumentError("order of the method must be strictly greater than that " *
+                                 "of the derivative"))
+    # Check whether the method can be computed. We require the factorial of the
+    # method order to be computable with regular `Int`s, but `factorial` will overflow
+    # after 20, so 20 is the largest we can allow.
+    p > 20 && throw(ArgumentError("order of the method is too large to be computed"))
+    return
+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)
+    x[q + 1] = factorial(q)
+    return C \ x
+end
+
+# Estimate the bound on the function value and its derivatives at a point
+_estimate_bound(x, cond) = cond * maximum(abs, x) + TINY
+
+
 """
     History
 
@@ -97,9 +122,11 @@ end
 # The below does not apply to Nonstandard, as it has its own constructor
 for D in (:Forward, :Backward, :Central)
     lcname = lowercase(String(D))
-    gridf = Symbol(lcname, "_grid")
-    fdmf = Symbol(lcname, "_fdm")
+    gridf = Symbol(lcname, :_grid)
+    fdmf = Symbol(lcname, :_fdm)
+    precomputes = Symbol(lcname, :_precomputed)
 
+    num_precomputed = 8
     @eval begin
         # Compatibility layer over the "old" API
         function $fdmf(p::Integer, q::Integer; adapt=1, kwargs...)
@@ -107,10 +134,23 @@ for D in (:Forward, :Backward, :Central)
             return $D(p, q; adapt=adapt, kwargs...)
         end
 
+        # precompute a bunch of grids and coefs
+        const $precomputes = ntuple($num_precomputed) do p
+            p == 1 && return nothing  # no grid define for p==1
+            grid = $gridf(p)
+            all_coefs = ntuple(q->_coefs(grid, p, q), p - 1)
+            (grid, all_coefs)
+        end
+
         function $D(p::Integer, q::Integer; adapt=1, kwargs...)
             _check_p_q(p, q)
-            grid = $gridf(p)
-            coefs = _coefs(grid, p, q)
+            if p <= $num_precomputed
+                grid, all_coefs = $precomputes[p]
+                coefs = all_coefs[q]
+            else
+                grid = $gridf(p)
+                coefs = _coefs(grid, p, q)
+            end
             hist = History(; adapt=adapt, kwargs...)
             return $D{typeof(grid), typeof(coefs)}(Int(p), Int(q), grid, coefs, hist)
         end
@@ -143,29 +183,6 @@ function Nonstandard(grid::AbstractVector{<:Real}, q::Integer; adapt=0, kwargs..
     return Nonstandard{typeof(grid), typeof(coefs)}(Int(p), Int(q), grid, coefs, hist)
 end
 
-# Check the method and derivative orders for consistency
-function _check_p_q(p::Integer, q::Integer)
-    q >= 0 || throw(ArgumentError("order of derivative must be nonnegative"))
-    q < p || throw(ArgumentError("order of the method must be strictly greater than that " *
-                                 "of the derivative"))
-    # Check whether the method can be computed. We require the factorial of the
-    # method order to be computable with regular `Int`s, but `factorial` will overflow
-    # after 20, so 20 is the largest we can allow.
-    p > 20 && throw(ArgumentError("order of the method is too large to be computed"))
-    return
-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)
-    x[q + 1] = factorial(q)
-    return C \ x
-end
-
-# Estimate the bound on the function value and its derivatives at a point
-_estimate_bound(x, cond) = cond * maximum(abs, x) + TINY
-
 """
     fdm(m::FiniteDifferenceMethod, f, x[, Val(false)]; kwargs...) -> Real
     fdm(m::FiniteDifferenceMethod, f, x, Val(true); kwargs...) -> Tuple{FiniteDifferenceMethod, Real}
@@ -181,6 +198,7 @@ The recognized keywords are:
 * `bound`: Bound on the value of the function and its derivatives at `x`.
 * `condition`: The condition number. See [`DEFAULT_CONDITION`](@ref).
 * `eps`: The assumed roundoff error. Defaults to `eps()` plus [`TINY`](@ref).
+* `track_history`: wether to update the history of the method `m` with e.g. accuracy stats.
 
 !!! warning
     Bounds can't be adaptively computed over nonstandard grids; passing a value for
@@ -211,14 +229,15 @@ julia> fdm(central_fdm(2, 1), exp, 0, Val(true))
 function fdm(
     m::M,
     f,
-    x,
+    x::T,
     ::Val{true};
     condition=DEFAULT_CONDITION,
-    bound=_estimate_bound(f(x), condition),
-    eps=(Base.eps(float(bound)) + TINY),
+    bound::T=_estimate_bound(f(x)::T, condition)::T,
+    eps::T=(Base.eps(float(bound)) + TINY)::T,
     adapt=m.history.adapt,
     max_step=0.1,
-) where M<:FiniteDifferenceMethod
+    track_history=true  # will be set to false if `Val{false}()` used
+)::Tuple{M, T} where {M<:FiniteDifferenceMethod, T}
     if M <: Nonstandard && adapt > 0
         throw(ArgumentError("can't adaptively compute bounds over Nonstandard grids"))
     end
@@ -237,7 +256,7 @@ function fdm(
 
     # Adaptively compute the bound on the function and derivative values, if applicable.
     if adapt > 0
-        newm = (M.name.wrapper)(p + 1, p)
+        newm = adapted_FDM(m)
         dfdx = fdm(
             newm,
             f,
@@ -256,25 +275,34 @@ function fdm(
     C₂ = bound * sum(n->abs(coefs[n] * grid[n]^p), eachindex(coefs)) / factorial(p)
     ĥ = min((q / (p - q) * C₁ / C₂)^(1 / p), max_step)
 
-    # Estimate the accuracy of the method.
-    accuracy = ĥ^(-q) * C₁ + ĥ^(p - q) * C₂
-
     # Estimate the value of the derivative.
-    dfdx = sum(i->coefs[i] * f(x + ĥ * grid[i]), eachindex(grid)) / ĥ^q
-
-    m.history.eps = eps
-    m.history.bound = bound
-    m.history.step = ĥ
-    m.history.accuracy = accuracy
-
+    dfdx_s = sum(eachindex(grid)) do i
+        (@inbounds coefs[i] * f(x + ĥ * grid[i]))::T
+    end
+    dfdx = dfdx_s / ĥ^q
+
+    if track_history
+        # Estimate the accuracy of the method.
+        accuracy = ĥ^(-q) * C₁ + ĥ^(p - q) * C₂
+        m.history.eps = eps
+        m.history.bound = bound
+        m.history.step = ĥ
+        m.history.accuracy = accuracy
+    end
     return m, dfdx
 end
 
 function fdm(m::FiniteDifferenceMethod, f, x, ::Val{false}=Val(false); kwargs...)
-    _, dfdx = fdm(m, f, x, Val(true); kwargs...)
+    _, dfdx = fdm(m, f, x, Val{true}(); track_history=false, kwargs...)
     return dfdx
 end
 
+_adapted_FDM(m::M) where M = (M.name.wrapper)(m.p + 1, m.p)
+adapted_FDM(m::M) where M = _adapted_FDM(m)
+# Deal with type-instability in Central:
+adapted_FDM(m::Central{Vector{Int}}) where M = _adapted_FDM(m)::Central{UnitRange{Int}}
+adapted_FDM(m::Central{UnitRange{Int}}) where M = _adapted_FDM(m)::Central{Vector{Int}}
+
 
 ## Deprecations