From 666a1eac6695f9bb0592bb48de053193d21734ac Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Thu, 28 Aug 2025 11:41:47 +0300
Subject: [PATCH 01/18] Add Unitful.jl to dependencies

---
 Project.toml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/Project.toml b/Project.toml
index d6c26eb..ae8e87b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Richardson = "708f8203-808e-40c0-ba2d-98a6953ed40d"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
 ChainRulesCore = "0.10.3, 1"
@@ -19,6 +20,7 @@ Random = "<0.0.1, 1"
 Richardson = "1.2"
 SparseArrays = "<0.0.1, 1"
 StaticArrays = "0.12, 1.0"
+Unitful = "1.24.0"
 julia = "1"
 
 [extras]

From 633ec9b0cdbeabffe9513948bfac4c43a7f2b3ca Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Thu, 28 Aug 2025 13:15:58 +0300
Subject: [PATCH 02/18] Use Unitful in package root file

---
 src/FiniteDifferences.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/FiniteDifferences.jl b/src/FiniteDifferences.jl
index 82dd64c..21f8cd3 100644
--- a/src/FiniteDifferences.jl
+++ b/src/FiniteDifferences.jl
@@ -7,6 +7,7 @@ using Random
 using Richardson
 using SparseArrays
 using StaticArrays
+using Unitful
 
 export to_vec, grad, jacobian, jvp, j′vp
 

From 15eeec13273f4d8c33304de16fa9b838d6469761 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Thu, 28 Aug 2025 13:19:11 +0300
Subject: [PATCH 03/18] methods: change all instances of Real and AbstractFloat
 to Number

... in the hopes of supporting Unitful numbers, since they are <: Number.
This alone does not work though.
---
 src/methods.jl | 92 +++++++++++++++++++++++++-------------------------
 1 file changed, 46 insertions(+), 46 deletions(-)

diff --git a/src/methods.jl b/src/methods.jl
index f3d8ea7..0cea1ce 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -107,9 +107,9 @@ end
     FiniteDifferenceMethod(
         grid::AbstractVector{Int},
         q::Int;
-        condition::Real=DEFAULT_CONDITION,
-        factor::Real=DEFAULT_FACTOR,
-        max_range::Real=Inf
+        condition::Number=DEFAULT_CONDITION,
+        factor::Number=DEFAULT_FACTOR,
+        max_range::Number=Inf
     )
 
 Construct a finite difference method.
@@ -120,9 +120,9 @@ Construct a finite difference method.
 - `q::Int`: Order of the derivative to estimate.
 
 # Keywords
-- `condition::Real`: Condition number. See [`DEFAULT_CONDITION`](@ref).
-- `factor::Real`: Factor number. See [`DEFAULT_FACTOR`](@ref).
-- `max_range::Real=Inf`: Maximum distance that a function is evaluated from the input at
+- `condition::Number`: Condition number. See [`DEFAULT_CONDITION`](@ref).
+- `factor::Number`: Factor number. See [`DEFAULT_FACTOR`](@ref).
+- `max_range::Number=Inf`: Maximum distance that a function is evaluated from the input at
     which the derivative is estimated.
 
 # Returns
@@ -131,9 +131,9 @@ Construct a finite difference method.
 function FiniteDifferenceMethod(
     grid::SVector{P,Int},
     q::Int;
-    condition::Real=DEFAULT_CONDITION,
-    factor::Real=DEFAULT_FACTOR,
-    max_range::Real=Inf,
+    condition::Number=DEFAULT_CONDITION,
+    factor::Number=DEFAULT_FACTOR,
+    max_range::Number=Inf,
 ) where P
     _check_p_q(P, q)
     coefs, coefs_neighbourhood, ∇f_magnitude_mult, f_error_mult = _coefs_mults(grid, q)
@@ -153,7 +153,7 @@ function FiniteDifferenceMethod(grid::AbstractVector{Int}, q::Int; kw_args...)
 end
 
 """
-    (m::FiniteDifferenceMethod)(f, x::T) where T<:AbstractFloat
+    (m::FiniteDifferenceMethod)(f, x::T) where T<:Number
 
 Estimate the derivative of `f` at `x` using the finite differencing method `m` and an
 automatically determined step size.
@@ -188,12 +188,12 @@ julia> FiniteDifferences.estimate_step(fdm, sin, 1.0)  # Computes step size and
 # We loop over all concrete subtypes of `FiniteDifferenceMethod` for Julia v1.0 compatibility.
 for T in (UnadaptedFiniteDifferenceMethod, AdaptedFiniteDifferenceMethod)
     @eval begin
-        function (m::$T)(f::TF, x::Real) where TF
+        function (m::$T)(f::TF, x::Number) where TF
             x = float(x)  # Assume that converting to float is desired, if it isn't already.
             step = first(estimate_step(m, f, x))
             return m(f, x, step)
         end
-        function (m::$T{P,0})(f::TF, x::Real) where {P,TF}
+        function (m::$T{P,0})(f::TF, x::Number) where {P,TF}
             # The automatic step size calculation fails if `Q == 0`, so handle that edge
             # case.
             return f(x)
@@ -202,7 +202,7 @@ for T in (UnadaptedFiniteDifferenceMethod, AdaptedFiniteDifferenceMethod)
 end
 
 """
-    (m::FiniteDifferenceMethod)(f, x::T, step::Real) where T<:AbstractFloat
+    (m::FiniteDifferenceMethod)(f, x::T, step::Number) where T<:Number
 
 Estimate the derivative of `f` at `x` using the finite differencing method `m` and a given
 step size.
@@ -210,7 +210,7 @@ step size.
 # Arguments
 - `f`: Function to estimate derivative of.
 - `x::T`: Input to estimate derivative at.
-- `step::Real`: Step size.
+- `step::Number`: Step size.
 
 # Returns
 - Estimate of the derivative.
@@ -235,7 +235,7 @@ julia> fdm(sin, 1, 1e-3) - cos(1)  # Check the error.
 # We loop over all concrete subtypes of `FiniteDifferenceMethod` for 1.0 compatibility.
 for T in (UnadaptedFiniteDifferenceMethod, AdaptedFiniteDifferenceMethod)
     @eval begin
-        function (m::$T{P,Q})(f::TF, x::Real, step::Real) where {P,Q,TF}
+        function (m::$T{P,Q})(f::TF, x::Number, step::Number) where {P,Q,TF}
             x = float(x)  # Assume that converting to float is desired, if it isn't already.
             fs = _eval_function(m, f, x, step)
             return _compute_estimate(m, fs, x, step, m.coefs)
@@ -244,8 +244,8 @@ for T in (UnadaptedFiniteDifferenceMethod, AdaptedFiniteDifferenceMethod)
 end
 
 function _eval_function(
-    m::FiniteDifferenceMethod, f::TF, x::T, step::Real,
-) where {TF,T<:AbstractFloat}
+    m::FiniteDifferenceMethod, f::TF, x::T, step::Number,
+) where {TF,T<:Number}
     return f.(x .+ T(step) .* m.grid)
 end
 
@@ -253,9 +253,9 @@ function _compute_estimate(
     m::FiniteDifferenceMethod{P,Q},
     fs::SVector{P,TF},
     x::T,
-    step::Real,
+    step::Number,
     coefs::SVector{P,Float64},
-) where {P,Q,TF,T<:AbstractFloat}
+) where {P,Q,TF,T<:Number}
     # If we substitute `T.(coefs)` in the expression below, then allocations occur. We
     # therefore perform the broadcasting first. See
     # https://github.com/JuliaLang/julia/issues/39151.
@@ -338,7 +338,7 @@ end
         m::FiniteDifferenceMethod,
         f,
         x::T
-    ) where T<:AbstractFloat
+    ) where T<:Number
 
 Estimate the step size for a finite difference method `m`. Also estimates the error of the
 estimate of the derivative.
@@ -349,19 +349,19 @@ estimate of the derivative.
 - `x::T`: Point to estimate the derivative at.
 
 # Returns
-- `Tuple{<:AbstractFloat, <:AbstractFloat}`: Estimated step size and an estimate of the
+- `Tuple{<:Number, <:Number}`: Estimated step size and an estimate of the
     error of the finite difference estimate. The error will be `NaN` if the method failed
     to estimate the error.
 """
 function estimate_step(
     m::UnadaptedFiniteDifferenceMethod, f::TF, x::T,
-) where {TF,T<:AbstractFloat}
+) where {TF,T<:Number}
     step, acc = _compute_step_acc_default(m, x)
     return _limit_step(m, x, step, acc)
 end
 function estimate_step(
     m::AdaptedFiniteDifferenceMethod{P,Q}, f::TF, x::T,
-) where {P,Q,TF,T<:AbstractFloat}
+) where {P,Q,TF,T<:Number}
     ∇f_magnitude, f_magnitude = _estimate_magnitudes(m.bound_estimator, f, x)
     if ∇f_magnitude == 0.0 || f_magnitude == 0.0
         step, acc = _compute_step_acc_default(m, x)
@@ -373,7 +373,7 @@ end
 
 function _estimate_magnitudes(
     m::FiniteDifferenceMethod{P,Q}, f::TF, x::T,
-) where {P,Q,TF,T<:AbstractFloat}
+) where {P,Q,TF,T<:Number}
     step = first(estimate_step(m, f, x))
     fs = _eval_function(m, f, x, step)
     # Estimate magnitude of `∇f` in a neighbourhood of `x`.
@@ -388,13 +388,13 @@ function _estimate_magnitudes(
     return ∇f_magnitude, f_magnitude
 end
 
-function _compute_step_acc_default(m::FiniteDifferenceMethod, x::T) where {T<:AbstractFloat}
+function _compute_step_acc_default(m::FiniteDifferenceMethod, x::T) where {T<:Number}
     # Compute a default step size using a heuristic and [`DEFAULT_CONDITION`](@ref).
     return _compute_step_acc(m, m.condition, eps(T))
 end
 
 function _compute_step_acc(
-    m::FiniteDifferenceMethod{P,Q}, ∇f_magnitude::Real, f_error::Real,
+    m::FiniteDifferenceMethod{P,Q}, ∇f_magnitude::Number, f_error::Number,
 ) where {P,Q}
     # Set the step size by minimising an upper bound on the error of the estimate.
     C₁ = f_error * m.f_error_mult * m.factor
@@ -406,8 +406,8 @@ function _compute_step_acc(
 end
 
 function _limit_step(
-    m::FiniteDifferenceMethod, x::T, step::Real, acc::Real,
-) where {T<:AbstractFloat}
+    m::FiniteDifferenceMethod, x::T, step::Number, acc::Number,
+) where {T<:Number}
     # First, limit the step size based on the maximum range.
     step_max = m.max_range / maximum(abs.(m.grid))
     if step > step_max
@@ -433,9 +433,9 @@ for direction in [:forward, :central, :backward]
             p::Int,
             q::Int;
             adapt::Int=1,
-            condition::Real=DEFAULT_CONDITION,
-            factor::Real=DEFAULT_FACTOR,
-            max_range::Real=Inf,
+            condition::Number=DEFAULT_CONDITION,
+            factor::Number=DEFAULT_FACTOR,
+            max_range::Number=Inf,
             geom::Bool=false
         )
             _check_p_q(p, q)
@@ -484,9 +484,9 @@ for direction in [:forward, :central, :backward]
         p::Int,
         q::Int;
         adapt::Int=1,
-        condition::Real=DEFAULT_CONDITION,
-        factor::Real=DEFAULT_FACTOR,
-        max_range::Real=Inf,
+        condition::Number=DEFAULT_CONDITION,
+        factor::Number=DEFAULT_FACTOR,
+        max_range::Number=Inf,
         geom::Bool=false
     )
 
@@ -500,9 +500,9 @@ Construct a finite difference method at a $($(Meta.quot(direction))) grid of `p`
 - `adapt::Int=1`: Use another finite difference method to estimate the magnitude of the
     `p`th order derivative, which is important for the step size computation. Recurse
     this procedure `adapt` times.
-- `condition::Real`: Condition number. See [`DEFAULT_CONDITION`](@ref).
-- `factor::Real`: Factor number. See [`DEFAULT_FACTOR`](@ref).
-- `max_range::Real=Inf`: Maximum distance that a function is evaluated from the input at
+- `condition::Number`: Condition number. See [`DEFAULT_CONDITION`](@ref).
+- `factor::Number`: Factor number. See [`DEFAULT_FACTOR`](@ref).
+- `max_range::Number=Inf`: Maximum distance that a function is evaluated from the input at
     which the derivative is estimated.
 - `geom::Bool`: Use geometrically spaced points instead of linearly spaced points.
 
@@ -552,10 +552,10 @@ end
     extrapolate_fdm(
         m::FiniteDifferenceMethod,
         f,
-        x::Real,
-        initial_step::Real=10,
+        x::Number,
+        initial_step::Number=10,
         power::Int=1,
-        breaktol::Real=Inf,
+        breaktol::Number=Inf,
         kw_args...
     )
 
@@ -568,19 +568,19 @@ automatically sets `power = 2` if `m` is symmetric and `power = 1`. Moreover, it
 # Arguments
 - `m::FiniteDifferenceMethod`: Finite difference method to estimate the step size for.
 - `f`: Function to evaluate the derivative of.
-- `x::Real`: Point to estimate the derivative at.
-- `initial_step::Real=10`: Initial step size.
+- `x::Number`: Point to estimate the derivative at.
+- `initial_step::Number=10`: Initial step size.
 
 # Returns
-- `Tuple{<:AbstractFloat, <:AbstractFloat}`: Estimate of the derivative and error.
+- `Tuple{<:Number, <:Number}`: Estimate of the derivative and error.
 """
 function extrapolate_fdm(
     m::FiniteDifferenceMethod,
     f,
-    x::Real,
-    initial_step::Real=10;
+    x::Number,
+    initial_step::Number=10;
     power::Int=1,
-    breaktol::Real=Inf,
+    breaktol::Number=Inf,
     kw_args...
 )
     (power == 1 && _is_symmetric(m)) && (power = 2)

From 311bf8a3e6faf6ca6f51b76051e0f3d695b388db Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Thu, 28 Aug 2025 13:30:16 +0300
Subject: [PATCH 04/18] _compute_step_acc: compute Pth root in step computation
 using Rationals instead of floats

This makes working with units more exact, since the powers of Unitful
units are always rational. Without this, we get powers such as
74981274124 / 1249141204412 for units.
---
 src/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/methods.jl b/src/methods.jl
index 0cea1ce..c3274a6 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -399,7 +399,7 @@ function _compute_step_acc(
     # Set the step size by minimising an upper bound on the error of the estimate.
     C₁ = f_error * m.f_error_mult * m.factor
     C₂ = ∇f_magnitude * m.∇f_magnitude_mult
-    step = (Q / (P - Q) * (C₁ / C₂))^(1 / P)
+    step = (Q / (P - Q) * (C₁ / C₂))^(1 // P)
     # Estimate the accuracy of the method.
     acc = C₁ * step^(-Q) + C₂ * step^(P - Q)
     return step, acc

From a1bfb2bcbed947848784d768363fec6b336a752f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Thu, 28 Aug 2025 16:29:07 +0300
Subject: [PATCH 05/18] _compute_estimate: change order of T() and ^Q in
 denominator

Running a simple test

    using Unitful: @u_str
    import Interpolations as Ilib
    import Interpolations
    import FiniteDifferences as FDlib

    tt = (0.0 : 0.01 : 1.1) .* u"s"

    fn = 3u"J / s" .* tt

    ii = Ilib.interpolate((tt,),fn, ComposedFunction(Ilib.Gridded, Ilib.Linear)())

    cdm = FDlib.central_fdm(5,1)

    @show dd = cdm(ii, 0.6u"s")

without this change in the order of operations causes the code to crash
in _limit_step, because of a dimension mismatch. This is a temporary "fix",
that actually not a fix. The given order of operations should really
work as-is...
---
 src/methods.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/src/methods.jl b/src/methods.jl
index c3274a6..e05dacd 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -260,7 +260,12 @@ function _compute_estimate(
     # therefore perform the broadcasting first. See
     # https://github.com/JuliaLang/julia/issues/39151.
     _coefs = T.(coefs)
-    return sum(fs .* _coefs) ./ T(step)^Q
+    # NOTE the denominator: while doing T(step)^Q should technically be possible and correct,
+    # without flipping the order of ^ Q and T(), the current code crashes in _limit_step
+    # because of dimension mismatch. However, the output now has wrong units.
+    #
+    # TODO: fix the output units. The numerical value is correct in a simple test.
+    return sum(fs .* _coefs) ./ T(step ^ Q)
 end
 
 # Check the method and derivative orders for consistency.

From b3fe4e724e2622009ef038f7635bd4ebc25657ab Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Fri, 29 Aug 2025 10:05:22 +0300
Subject: [PATCH 06/18] methods: define and use function withUnit to enforce
 correct units

There might still be a few places where units need to be enforced,
since comparing quantities and pure numbers does not work,
and my manual tests have just not hit those branches yet.
---
 src/methods.jl | 65 +++++++++++++++++++++++++++++++++++++-------------
 1 file changed, 48 insertions(+), 17 deletions(-)

diff --git a/src/methods.jl b/src/methods.jl
index e05dacd..e6ce6ae 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -259,13 +259,10 @@ function _compute_estimate(
     # If we substitute `T.(coefs)` in the expression below, then allocations occur. We
     # therefore perform the broadcasting first. See
     # https://github.com/JuliaLang/julia/issues/39151.
-    _coefs = T.(coefs)
-    # NOTE the denominator: while doing T(step)^Q should technically be possible and correct,
-    # without flipping the order of ^ Q and T(), the current code crashes in _limit_step
-    # because of dimension mismatch. However, the output now has wrong units.
     #
-    # TODO: fix the output units. The numerical value is correct in a simple test.
-    return sum(fs .* _coefs) ./ T(step ^ Q)
+    # We strip units because the estimate coefficients are just weights for values of f.
+    _coefs = ustrip(T.(coefs))
+    return sum(fs .* _coefs) ./ T(step) ^ Q
 end
 
 # Check the method and derivative orders for consistency.
@@ -361,18 +358,27 @@ estimate of the derivative.
 function estimate_step(
     m::UnadaptedFiniteDifferenceMethod, f::TF, x::T,
 ) where {TF,T<:Number}
-    step, acc = _compute_step_acc_default(m, x)
+    step, acc = withUnit.(
+        unit(x),
+        _compute_step_acc_default(m, x) .* unit(x)
+    )
     return _limit_step(m, x, step, acc)
 end
 function estimate_step(
     m::AdaptedFiniteDifferenceMethod{P,Q}, f::TF, x::T,
 ) where {P,Q,TF,T<:Number}
-    ∇f_magnitude, f_magnitude = _estimate_magnitudes(m.bound_estimator, f, x)
-    if ∇f_magnitude == 0.0 || f_magnitude == 0.0
-        step, acc = _compute_step_acc_default(m, x)
-    else
-        step, acc = _compute_step_acc(m, ∇f_magnitude, eps(f_magnitude))
-    end
+    @show ∇f_magnitude, f_magnitude = _estimate_magnitudes(m.bound_estimator, f, x)
+    step, acc = withUnit.(
+        (
+            unit(x),
+            unit(f |> eltype) / unit(x) ^ Q
+        ),
+        if ∇f_magnitude == 0.0 || f_magnitude == 0.0
+            _compute_step_acc_default(m, x)
+        else
+            _compute_step_acc(m, ∇f_magnitude, eps(f_magnitude))
+        end
+    )
     return _limit_step(m, x, step, acc)
 end
 
@@ -413,19 +419,26 @@ end
 function _limit_step(
     m::FiniteDifferenceMethod, x::T, step::Number, acc::Number,
 ) where {T<:Number}
+    xunit = unit(x)
     # First, limit the step size based on the maximum range.
-    step_max = m.max_range / maximum(abs.(m.grid))
+    step_max = withUnit(
+        xunit,
+        m.max_range / maximum(abs.(m.grid))
+    )
     if step > step_max
         step = step_max
-        acc = NaN
+        acc = withUnit(xunit,NaN)
     end
     # Second, prevent very large step sizes, which can occur for high-order methods or
     # slowly-varying functions.
-    step_default, _ = _compute_step_acc_default(m, x)
+    step_default, _ = withUnit.(
+        xunit,
+        _compute_step_acc_default(m, x)
+    )
     step_max_default = 1000step_default
     if step > step_max_default
         step = step_max_default
-        acc = NaN
+        acc = withUnit(xunit,NaN)
     end
     return step, acc
 end
@@ -597,3 +610,21 @@ function extrapolate_fdm(
         kw_args...
     )
 end
+
+"""
+Attaches a given unit to a given value, if it is dimensionless.
+If the value is not dimensionless, attempts a conversion to the given unit.
+"""
+function withUnit(targetUnit, value)
+
+    if Unitful.dimension(value) == Unitful.NoDims
+
+        value .* targetUnit
+
+    else
+
+        Unitful.uconvert(targetUnit, value)
+
+    end # if
+
+end # function

From c22698606985433b8307c4ae13ebed96a7835023 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Fri, 29 Aug 2025 10:20:35 +0300
Subject: [PATCH 07/18] methods: add missing withUnit calls to 2 numerical
 comparisons

---
 src/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/methods.jl b/src/methods.jl
index e6ce6ae..9d2df65 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -373,7 +373,7 @@ function estimate_step(
             unit(x),
             unit(f |> eltype) / unit(x) ^ Q
         ),
-        if ∇f_magnitude == 0.0 || f_magnitude == 0.0
+        if ∇f_magnitude == withUnit(∇f_magnitude,0.0) || f_magnitude == withUnit(unit(f_magnitude), 0.0)
             _compute_step_acc_default(m, x)
         else
             _compute_step_acc(m, ∇f_magnitude, eps(f_magnitude))

From 9f74bed42ef22d2c1513112abe50078f8eee764f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Fri, 29 Aug 2025 10:21:29 +0300
Subject: [PATCH 08/18] methods: remove a @show invocation

---
 src/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/methods.jl b/src/methods.jl
index 9d2df65..1b4d92a 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -367,7 +367,7 @@ end
 function estimate_step(
     m::AdaptedFiniteDifferenceMethod{P,Q}, f::TF, x::T,
 ) where {P,Q,TF,T<:Number}
-    @show ∇f_magnitude, f_magnitude = _estimate_magnitudes(m.bound_estimator, f, x)
+    ∇f_magnitude, f_magnitude = _estimate_magnitudes(m.bound_estimator, f, x)
     step, acc = withUnit.(
         (
             unit(x),

From b02fadb5c37bd3ffd5c075d6ba9340922068d464 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Fri, 29 Aug 2025 10:29:31 +0300
Subject: [PATCH 09/18] methods: as per a deprecation warning, call
 broadcasting version of ustrip and remove extra .* unit(x)

---
 src/methods.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/methods.jl b/src/methods.jl
index 1b4d92a..32b06b9 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -261,7 +261,7 @@ function _compute_estimate(
     # https://github.com/JuliaLang/julia/issues/39151.
     #
     # We strip units because the estimate coefficients are just weights for values of f.
-    _coefs = ustrip(T.(coefs))
+    _coefs = ustrip.(T.(coefs))
     return sum(fs .* _coefs) ./ T(step) ^ Q
 end
 
@@ -360,7 +360,7 @@ function estimate_step(
 ) where {TF,T<:Number}
     step, acc = withUnit.(
         unit(x),
-        _compute_step_acc_default(m, x) .* unit(x)
+        _compute_step_acc_default(m, x)
     )
     return _limit_step(m, x, step, acc)
 end

From a902f0fedb0fd7916917e6ad0ca939cca0ef2e43 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Mon, 1 Sep 2025 12:45:48 +0300
Subject: [PATCH 10/18] Add test for unitful output (needs fixing, see details)

For some reason, the input value x of the fn being differentiated should
have the inverse units of the derivative, for the units to come out
correctly. This should be lookied into...
---
 test/methods.jl  | 8 ++++++++
 test/runtests.jl | 1 +
 2 files changed, 9 insertions(+)

diff --git a/test/methods.jl b/test/methods.jl
index af8e15b..799c48a 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -89,6 +89,14 @@ struct NotAFunction end # not <: Function on purpose, cf #224
         @test isapprox(central_fdm(5, 1)(x -> 5, 0), 0; rtol=1e-12, atol=1e-12)
     end
 
+    # Integration test to ensure that Unitful-output functions can be tested.
+    @testset "Unitful output" begin
+        fn(x) = 5u"J/s" * x
+        @show derivativeVal = central_fdm(5, 1)(fn, 1u"s/J")
+        @test unit(derivativeVal) == u"J/s"
+        @test isapprox(derivativeVal, 5u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
+    end
+
     @testset "Adaptation improves estimate" begin
         @test abs(forward_fdm(5, 1, adapt=0)(log, 0.001) - 1000) > 1
         @test forward_fdm(5, 1, adapt=1)(log, 0.001) ≈ 1000
diff --git a/test/runtests.jl b/test/runtests.jl
index 63039f7..a2d524f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,6 +7,7 @@ using Random
 using SparseArrays
 using StaticArrays
 using Test
+using Unitful
 
 
 Random.seed!(1)

From 81b51322b8202d7f7334ae6c861da56e379c6e8e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Mon, 1 Sep 2025 12:50:21 +0300
Subject: [PATCH 11/18] test/methods: do the (still broken) Unitful test with
 floats and not integers (although both work)

Again, the input value x to fn needs to have non-sensible units for the
desired units to come out of the finite difference method.
---
 test/methods.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/methods.jl b/test/methods.jl
index 799c48a..d4365aa 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -91,10 +91,10 @@ struct NotAFunction end # not <: Function on purpose, cf #224
 
     # Integration test to ensure that Unitful-output functions can be tested.
     @testset "Unitful output" begin
-        fn(x) = 5u"J/s" * x
-        @show derivativeVal = central_fdm(5, 1)(fn, 1u"s/J")
+        fn(x) = 5.0u"J/s" * x
+        @show derivativeVal = central_fdm(5, 1)(fn, 1.0u"s/J")
         @test unit(derivativeVal) == u"J/s"
-        @test isapprox(derivativeVal, 5u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
+        @test isapprox(derivativeVal, 5.0u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
     end
 
     @testset "Adaptation improves estimate" begin

From 33a4af7d10fcb24c56747c898e2ac6af5dc1a707 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Mon, 1 Sep 2025 13:40:36 +0300
Subject: [PATCH 12/18] test/methods: (in progress) add explicit adapt
 parameter cases to Unitful tests

These seem to at least produce correct units, but the numerical value is
wrong for the adapt=0 case (not surprising, based on the documentation).

The question is, why the heck do things go wrong, when you do not give
an explicit adapt value, but leave it to the compiler-generated method
that inserts the value adapt=1 automatically.
---
 test/methods.jl | 23 ++++++++++++++++++++---
 1 file changed, 20 insertions(+), 3 deletions(-)

diff --git a/test/methods.jl b/test/methods.jl
index d4365aa..c8f31d9 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -91,10 +91,27 @@ struct NotAFunction end # not <: Function on purpose, cf #224
 
     # Integration test to ensure that Unitful-output functions can be tested.
     @testset "Unitful output" begin
+
         fn(x) = 5.0u"J/s" * x
-        @show derivativeVal = central_fdm(5, 1)(fn, 1.0u"s/J")
-        @test unit(derivativeVal) == u"J/s"
-        @test isapprox(derivativeVal, 5.0u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
+
+        derivativeVal = 5.0u"J/s"
+
+        # NOTE: the input value unit is wrong here to make this run at all!
+        # The unit should be u"s", but it has to be the inverse of the desired output unit.
+        # TODO: Find out the faulty code path using a debugger.
+        twoPlaceDerivativeVal = central_fdm(5, 1)(fn, 1.0u"s/J")
+
+        unadaptedDerivativeVal = central_fdm(5, 1)(fn, 1.0u"s",0)
+
+        adaptedDerivativeVal = central_fdm(5, 1)(fn, 1.0u"s",1)
+
+        @test unit(twoPlaceDerivativeVal) == u"J/s"
+        @test unit(unadaptedDerivativeVal) == u"J/s"
+        @test unit(adaptedDerivativeVal) == u"J/s"
+
+        @test isapprox(twoPlaceDerivativeVal, 5.0u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
+        @test isapprox(unadaptedDerivativeVal, 5.0u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
+        @test isapprox(adaptedDerivativeVal, 5.0u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
     end
 
     @testset "Adaptation improves estimate" begin

From 26431a33e155b4b99bd7772bb3edcfd89fe5d797 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Mon, 1 Sep 2025 15:37:58 +0300
Subject: [PATCH 13/18] methods: assume that f argument of estimate_step is a
 function-like object instead of a container

This means not relying on eltype.
---
 src/methods.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/methods.jl b/src/methods.jl
index 32b06b9..6d130cb 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -371,9 +371,9 @@ function estimate_step(
     step, acc = withUnit.(
         (
             unit(x),
-            unit(f |> eltype) / unit(x) ^ Q
+            unit(first(f(x))) / unit(x) ^ Q
         ),
-        if ∇f_magnitude == withUnit(∇f_magnitude,0.0) || f_magnitude == withUnit(unit(f_magnitude), 0.0)
+        if ∇f_magnitude == withUnit(unit(∇f_magnitude),0.0) || f_magnitude == withUnit(unit(f_magnitude), 0.0)
             _compute_step_acc_default(m, x)
         else
             _compute_step_acc(m, ∇f_magnitude, eps(f_magnitude))

From 3bf342b32a8feb835774c32773e75345e159db97 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Mon, 1 Sep 2025 15:47:59 +0300
Subject: [PATCH 14/18] test/methods: assume correct unit in Unitful test and
 remove useless test

Since the commit where the assumption about the function-likeness of
differentiated objects was made in step estimation, the correct unit
actually works.

The step with adaptive=0 was also useless, as a method with no
iterations does not converge, of course.
---
 test/methods.jl | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/test/methods.jl b/test/methods.jl
index c8f31d9..458c5c5 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -99,18 +99,14 @@ struct NotAFunction end # not <: Function on purpose, cf #224
         # NOTE: the input value unit is wrong here to make this run at all!
         # The unit should be u"s", but it has to be the inverse of the desired output unit.
         # TODO: Find out the faulty code path using a debugger.
-        twoPlaceDerivativeVal = central_fdm(5, 1)(fn, 1.0u"s/J")
-
-        unadaptedDerivativeVal = central_fdm(5, 1)(fn, 1.0u"s",0)
+        twoPlaceDerivativeVal = central_fdm(5, 1)(fn, 1.0u"s")
 
         adaptedDerivativeVal = central_fdm(5, 1)(fn, 1.0u"s",1)
 
         @test unit(twoPlaceDerivativeVal) == u"J/s"
-        @test unit(unadaptedDerivativeVal) == u"J/s"
         @test unit(adaptedDerivativeVal) == u"J/s"
 
         @test isapprox(twoPlaceDerivativeVal, 5.0u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
-        @test isapprox(unadaptedDerivativeVal, 5.0u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
         @test isapprox(adaptedDerivativeVal, 5.0u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
     end
 

From 9d049dadf49954677c5369e1f0500caf55fdc0b5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Mon, 1 Sep 2025 16:06:01 +0300
Subject: [PATCH 15/18] withUnit: do not even attempt broadcasted
 multiplication

This function itself should be called in a broadcasted manner, if
broadcasting is desired.
---
 src/methods.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/methods.jl b/src/methods.jl
index 6d130cb..880a3f4 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -619,7 +619,7 @@ function withUnit(targetUnit, value)
 
     if Unitful.dimension(value) == Unitful.NoDims
 
-        value .* targetUnit
+        value * targetUnit
 
     else
 

From 2b447e09b3f1b7473de46a6ac8f2d386c46d3187 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Mon, 1 Sep 2025 16:09:24 +0300
Subject: [PATCH 16/18] Remove a redundant Unitful test and clean up TODO
 comment

---
 test/methods.jl | 10 ++--------
 1 file changed, 2 insertions(+), 8 deletions(-)

diff --git a/test/methods.jl b/test/methods.jl
index 458c5c5..574a0c5 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -96,18 +96,12 @@ struct NotAFunction end # not <: Function on purpose, cf #224
 
         derivativeVal = 5.0u"J/s"
 
-        # NOTE: the input value unit is wrong here to make this run at all!
-        # The unit should be u"s", but it has to be the inverse of the desired output unit.
-        # TODO: Find out the faulty code path using a debugger.
-        twoPlaceDerivativeVal = central_fdm(5, 1)(fn, 1.0u"s")
-
         adaptedDerivativeVal = central_fdm(5, 1)(fn, 1.0u"s",1)
 
-        @test unit(twoPlaceDerivativeVal) == u"J/s"
         @test unit(adaptedDerivativeVal) == u"J/s"
 
-        @test isapprox(twoPlaceDerivativeVal, 5.0u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
-        @test isapprox(adaptedDerivativeVal, 5.0u"J/s"; rtol=1e-12, atol=1e-12u"J/s")
+        @test isapprox(adaptedDerivativeVal, derivativeVal; rtol=1e-12, atol=1e-12u"J/s")
+
     end
 
     @testset "Adaptation improves estimate" begin

From b6467500d61f247531b1f90ea37457204b7d27ac Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Wed, 3 Sep 2025 10:38:07 +0300
Subject: [PATCH 17/18] Reduce allocations by not calling dotted version of
 withUnit or ustrip

Instead, attach units to step and acc with two separate calls to
withUnit, and determine the number type N of method coefs intead of
first adding units and then stripping them.

A fair amount of allocations still remain for Unitful calls, though.
---
 src/methods.jl | 47 ++++++++++++++++++++++++++++++++---------------
 1 file changed, 32 insertions(+), 15 deletions(-)

diff --git a/src/methods.jl b/src/methods.jl
index 880a3f4..63d6314 100644
--- a/src/methods.jl
+++ b/src/methods.jl
@@ -261,7 +261,8 @@ function _compute_estimate(
     # https://github.com/JuliaLang/julia/issues/39151.
     #
     # We strip units because the estimate coefficients are just weights for values of f.
-    _coefs = ustrip.(T.(coefs))
+    N = numType(T)
+    _coefs = N.(coefs)
     return sum(fs .* _coefs) ./ T(step) ^ Q
 end
 
@@ -358,27 +359,26 @@ estimate of the derivative.
 function estimate_step(
     m::UnadaptedFiniteDifferenceMethod, f::TF, x::T,
 ) where {TF,T<:Number}
-    step, acc = withUnit.(
-        unit(x),
-        _compute_step_acc_default(m, x)
-    )
+    step, acc = _compute_step_acc_default(m, x)
+    xunit = unit(x)
+    step = withUnit(xunit, step)
+    acc = withUnit(xunit,acc)
     return _limit_step(m, x, step, acc)
 end
 function estimate_step(
     m::AdaptedFiniteDifferenceMethod{P,Q}, f::TF, x::T,
 ) where {P,Q,TF,T<:Number}
     ∇f_magnitude, f_magnitude = _estimate_magnitudes(m.bound_estimator, f, x)
-    step, acc = withUnit.(
-        (
-            unit(x),
-            unit(first(f(x))) / unit(x) ^ Q
-        ),
+    xunit = unit(x)
+    dfunit = unit(first(f(x))) / unit(x) ^ Q
+    step, acc =
         if ∇f_magnitude == withUnit(unit(∇f_magnitude),0.0) || f_magnitude == withUnit(unit(f_magnitude), 0.0)
             _compute_step_acc_default(m, x)
         else
             _compute_step_acc(m, ∇f_magnitude, eps(f_magnitude))
         end
-    )
+    step = withUnit(xunit, step)
+    acc = withUnit(dfunit, acc)
     return _limit_step(m, x, step, acc)
 end
 
@@ -431,10 +431,8 @@ function _limit_step(
     end
     # Second, prevent very large step sizes, which can occur for high-order methods or
     # slowly-varying functions.
-    step_default, _ = withUnit.(
-        xunit,
-        _compute_step_acc_default(m, x)
-    )
+    step_default, _ = _compute_step_acc_default(m, x)
+    step_default = withUnit(xunit, step_default)
     step_max_default = 1000step_default
     if step > step_max_default
         step = step_max_default
@@ -628,3 +626,22 @@ function withUnit(targetUnit, value)
     end # if
 
 end # function
+
+"""
+Retrieves the number type of a quantity, or returns the type itself in the case of a raw number.
+"""
+function numType(x::Number)
+    typeof(x)
+end # function
+
+function numType(x::Type{<:Number})
+    x
+end
+
+function numType(x::Unitful.AbstractQuantity)
+    Unitful.numtype(typeof(x))
+end # function
+
+function numType(x::Type{<:Unitful.AbstractQuantity})
+    Unitful.numtype(x)
+end

From 51d05ede00450b4ca00abb6143d49cb1dcd14065 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Santtu=20S=C3=B6derholm?= <santtu.soderholm@tuni.fi>
Date: Wed, 3 Sep 2025 10:43:29 +0300
Subject: [PATCH 18/18] Add a test for Unitful allocations

---
 test/methods.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/test/methods.jl b/test/methods.jl
index 574a0c5..37f6f4c 100644
--- a/test/methods.jl
+++ b/test/methods.jl
@@ -82,6 +82,8 @@ struct NotAFunction end # not <: Function on purpose, cf #224
     @testset "Test allocations" begin
         m = central_fdm(5, 2, adapt=2)
         @test @ballocated($m(sin, 1)) == 0
+        fn(t) = sin(t / u"s")
+        @test @ballocated($m($fn, 1u"s")) == 0
     end
 
     # Integration test to ensure that Integer-output functions can be tested.