diff --git a/src/grad.jl b/src/grad.jl
index 815a10a5..aef1811c 100644
--- a/src/grad.jl
+++ b/src/grad.jl
@@ -8,33 +8,40 @@ Approximate the gradient of `f` at `xs...` using `fdm`. Assumes that `f(xs...)`
 """
 function grad end
 
-function grad(fdm, f, x::AbstractArray{T}) where T <: Number
+function _grad(fdm, f, x::AbstractArray{T}) where T <: Number
+    # x must be mutable, we will mutate it and then mutate it back.
     dx = similar(x)
-    tmp = similar(x)
     for k in eachindex(x)
         dx[k] = fdm(zero(T)) do ϵ
-            tmp .= x
-            tmp[k] += ϵ
-            return f(tmp)
+            xk = x[k]
+            x[k] = xk +  ϵ
+            ret = f(x)
+            x[k] = xk  # Can't do `x[k] -= ϵ` as floating-point math is not associative
+            return ret::T
         end
     end
     return (dx, )
 end
 
+grad(fdm, f, x::Array{<:Number}) = _grad(fdm, f, x)
+# Fallback for when we don't know `x` will be mutable:
+grad(fdm, f, x::AbstractArray{<:Number}) = _grad(fdm, f, similar(x).=x)
+
 grad(fdm, f, x::Real) = (fdm(f, x), )
 grad(fdm, f, x::Tuple) = (grad(fdm, (xs...)->f(xs), x...), )
 
 function grad(fdm, f, d::Dict{K, V}) where {K, V}
-    dd = Dict{K, V}()
+    ∇d = Dict{K, V}()
     for (k, v) in d
+        dk = d[k]
         function f′(x)
-            tmp = copy(d)
-            tmp[k] = x
-            return f(tmp)
+            d[k] = x
+            return f(d)
         end
-        dd[k] = grad(fdm, f′, v)[1]
+        ∇d[k] = grad(fdm, f′, v)[1]
+        d[k] = dk
     end
-    return (dd, )
+    return (∇d, )
 end
 
 function grad(fdm, f, x)
diff --git a/src/methods.jl b/src/methods.jl
index 20b56ea9..e7911e32 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -181,6 +181,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`: whether 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
@@ -215,10 +216,11 @@ function fdm(
     ::Val{true};
     condition=DEFAULT_CONDITION,
     bound=_estimate_bound(f(x), condition),
-    eps=(Base.eps(float(bound)) + TINY),
+    eps::T=(Base.eps(float(bound)) + TINY),
     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
@@ -256,22 +258,25 @@ 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