diff --git a/Project.toml b/Project.toml
index 17fdd131..ee75d23d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -25,6 +25,7 @@ DataInterpolationsSymbolicsExt = "Symbolics"
 
 [compat]
 Aqua = "0.8"
+BenchmarkTools = "1"
 ChainRulesCore = "1.24"
 FindFirstFunctions = "1.1"
 FiniteDifferences = "0.12.31"
@@ -45,6 +46,7 @@ julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -58,4 +60,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Zygote"]
+test = ["Aqua", "BenchmarkTools", "SafeTestsets", "ChainRulesCore", "Optim", "RegularizationTools", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Zygote"]
diff --git a/docs/src/manual.md b/docs/src/manual.md
index 077cf29f..3818a9e3 100644
--- a/docs/src/manual.md
+++ b/docs/src/manual.md
@@ -14,3 +14,9 @@ CubicHermiteSpline
 PCHIPInterpolation
 QuinticHermiteSpline
 ```
+
+# Utility Functions
+
+```@docs
+DataInterpolations.looks_linear
+```
diff --git a/ext/DataInterpolationsChainRulesCoreExt.jl b/ext/DataInterpolationsChainRulesCoreExt.jl
index 9c33b09c..312e67e3 100644
--- a/ext/DataInterpolationsChainRulesCoreExt.jl
+++ b/ext/DataInterpolationsChainRulesCoreExt.jl
@@ -52,7 +52,7 @@ end
 
 function u_tangent(A::LinearInterpolation, t, Δ)
     out = zero(A.u)
-    idx = get_idx(A.t, t, A.idx_prev[])
+    idx = get_idx(A, t, A.idx_prev[])
     t_factor = (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx])
     out[idx] = Δ * (one(eltype(out)) - t_factor)
     out[idx + 1] = Δ * t_factor
diff --git a/src/derivatives.jl b/src/derivatives.jl
index 75872095..e7cb97fd 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -18,7 +18,7 @@ function derivative(A, t, order = 1)
 end
 
 function _derivative(A::LinearInterpolation, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; idx_shift = -1, ub_shift = -1, side = :first)
+    idx = get_idx(A, t, iguess; idx_shift = -1, ub_shift = -1, side = :first)
     slope = get_parameters(A, idx)
     slope, idx
 end
@@ -108,7 +108,7 @@ function _derivative(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, idx)
 end
 
 function _derivative(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; idx_shift = -1, side = :first)
+    idx = get_idx(A, t, iguess; idx_shift = -1, side = :first)
     j = min(idx, length(A.c))  # for smooth derivative at A.t[end]
     wj = t - A.t[idx]
     (@evalpoly wj A.b[idx] 2A.c[j] 3A.d[j]), idx
@@ -130,14 +130,14 @@ end
 
 # QuadraticSpline Interpolation
 function _derivative(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; lb = 2, ub_shift = 0, side = :first)
+    idx = get_idx(A, t, iguess; lb = 2, ub_shift = 0, side = :first)
     σ = get_parameters(A, idx - 1)
     A.z[idx - 1] + 2σ * (t - A.t[idx - 1]), idx
 end
 
 # CubicSpline Interpolation
 function _derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₁ = t - A.t[idx]
     Δt₂ = A.t[idx + 1] - t
     dI = (-A.z[idx] * Δt₂^2 + A.z[idx + 1] * Δt₁^2) / (2A.h[idx + 1])
@@ -151,7 +151,7 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num
     # change t into param [0 1]
     t < A.t[1] && return zero(A.u[1]), 1
     t > A.t[end] && return zero(A.u[end]), lastindex(t)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     n = length(A.t)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
@@ -173,7 +173,7 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
     # change t into param [0 1]
     t < A.t[1] && return zero(A.u[1]), 1
     t > A.t[end] && return zero(A.u[end]), lastindex(t)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     scale = (A.p[idx + 1] - A.p[idx]) / (A.t[idx + 1] - A.t[idx])
     t_ = A.p[idx] + (t - A.t[idx]) * scale
     N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N
@@ -192,7 +192,7 @@ end
 # Cubic Hermite Spline
 function _derivative(
         A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₀ = t - A.t[idx]
     Δt₁ = t - A.t[idx + 1]
     out = A.du[idx]
@@ -204,7 +204,7 @@ end
 # Quintic Hermite Spline
 function _derivative(
         A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₀ = t - A.t[idx]
     Δt₁ = t - A.t[idx + 1]
     out = A.du[idx] + A.ddu[idx] * Δt₀
diff --git a/src/integral_inverses.jl b/src/integral_inverses.jl
index 31d0853c..aebee9d5 100644
--- a/src/integral_inverses.jl
+++ b/src/integral_inverses.jl
@@ -59,7 +59,7 @@ end
 
 function _interpolate(
         A::LinearInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt = t - A.t[idx]
     x = A.itp.u[idx]
     slope = get_parameters(A.itp, idx)
@@ -104,13 +104,13 @@ end
 
 function _interpolate(
         A::ConstantInterpolationIntInv{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; ub_shift = 0)
+    idx = get_idx(A, t, iguess; ub_shift = 0)
     if A.itp.dir === :left
         # :left means that value to the left is used for interpolation
-        idx_ = get_idx(A.t, t, idx; lb = 1, ub_shift = 0)
+        idx_ = get_idx(A, t, idx; lb = 1, ub_shift = 0)
     else
         # :right means that value to the right is used for interpolation
-        idx_ = get_idx(A.t, t, idx; side = :first, lb = 1, ub_shift = 0)
+        idx_ = get_idx(A, t, idx; side = :first, lb = 1, ub_shift = 0)
     end
     A.u[idx] + (t - A.t[idx]) / A.itp.u[idx_], idx
 end
diff --git a/src/integrals.jl b/src/integrals.jl
index 03ea26c5..3ae893f7 100644
--- a/src/integrals.jl
+++ b/src/integrals.jl
@@ -8,9 +8,9 @@ function integral(A::AbstractInterpolation, t1::Number, t2::Number)
     ((t2 < A.t[1] || t2 > A.t[end]) && !A.extrapolate) && throw(ExtrapolationError())
     !hasfield(typeof(A), :I) && throw(IntegralNotFoundError())
     # the index less than or equal to t1
-    idx1 = get_idx(A.t, t1, 0)
+    idx1 = get_idx(A, t1, 0)
     # the index less than t2
-    idx2 = get_idx(A.t, t2, 0; idx_shift = -1, side = :first)
+    idx2 = get_idx(A, t2, 0; idx_shift = -1, side = :first)
 
     if A.cache_parameters
         total = A.I[idx2] - A.I[idx1]
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index 7f6991b8..71b2b042 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -12,7 +12,12 @@ Extrapolation extends the last linear polynomial on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
-  - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `cache_parameters`: precompute parameters at initialization for faster interpolation
+    computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
     u::uType
@@ -22,18 +27,22 @@ struct LinearInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolati
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    function LinearInterpolation(u, t, I, p, extrapolate, cache_parameters)
+    linear_lookup::Bool
+    function LinearInterpolation(u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.slope), eltype(u)}(
-            u, t, I, p, extrapolate, Ref(1), cache_parameters)
+            u, t, I, p, extrapolate, Ref(1), cache_parameters, linear_lookup)
     end
 end
 
-function LinearInterpolation(u, t; extrapolate = false, cache_parameters = false)
+function LinearInterpolation(
+        u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
     p = LinearParameterCache(u, t, cache_parameters)
-    A = LinearInterpolation(u, t, nothing, p, extrapolate, cache_parameters)
+    A = LinearInterpolation(
+        u, t, nothing, p, extrapolate, cache_parameters, assume_linear_t)
     I = cumulative_integral(A, cache_parameters)
-    LinearInterpolation(u, t, I, p, extrapolate, cache_parameters)
+    LinearInterpolation(u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
 end
 
 """
@@ -52,6 +61,10 @@ Extrapolation extends the last quadratic polynomial on each side.
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpolation{T}
     u::uType
@@ -62,24 +75,30 @@ struct QuadraticInterpolation{uType, tType, IType, pType, T} <: AbstractInterpol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    function QuadraticInterpolation(u, t, I, p, mode, extrapolate, cache_parameters)
+    linear_lookup::Bool
+    function QuadraticInterpolation(
+            u, t, I, p, mode, extrapolate, cache_parameters, assume_linear_t)
         mode ∈ (:Forward, :Backward) ||
             error("mode should be :Forward or :Backward for QuadraticInterpolation")
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.l₀), eltype(u)}(
-            u, t, I, p, mode, extrapolate, Ref(1), cache_parameters)
+            u, t, I, p, mode, extrapolate, Ref(1), cache_parameters, linear_lookup)
     end
 end
 
-function QuadraticInterpolation(u, t, mode; extrapolate = false, cache_parameters = false)
+function QuadraticInterpolation(
+        u, t, mode; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     p = QuadraticParameterCache(u, t, cache_parameters)
-    A = QuadraticInterpolation(u, t, nothing, p, mode, extrapolate, cache_parameters)
+    A = QuadraticInterpolation(
+        u, t, nothing, p, mode, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    QuadraticInterpolation(u, t, I, p, mode, extrapolate, cache_parameters)
+    QuadraticInterpolation(u, t, I, p, mode, extrapolate, cache_parameters, linear_lookup)
 end
 
-function QuadraticInterpolation(u, t; extrapolate = false, cache_parameters = false)
-    QuadraticInterpolation(u, t, :Forward; extrapolate, cache_parameters)
+function QuadraticInterpolation(u, t; kwargs...)
+    QuadraticInterpolation(u, t, :Forward; kwargs...)
 end
 
 """
@@ -145,6 +164,10 @@ Extrapolation extends the last cubic polynomial on each side.
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
        AbstractInterpolation{T}
@@ -157,7 +180,10 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    function AkimaInterpolation(u, t, I, b, c, d, extrapolate, cache_parameters)
+    linear_lookup::Bool
+    function AkimaInterpolation(
+            u, t, I, b, c, d, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(b), typeof(c),
             typeof(d), eltype(u)}(u,
             t,
@@ -167,13 +193,16 @@ struct AkimaInterpolation{uType, tType, IType, bType, cType, dType, T} <:
             d,
             extrapolate,
             Ref(1),
-            cache_parameters
+            cache_parameters,
+            linear_lookup
         )
     end
 end
 
-function AkimaInterpolation(u, t; extrapolate = false, cache_parameters = false)
+function AkimaInterpolation(
+        u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     n = length(t)
     dt = diff(t)
     m = Array{eltype(u)}(undef, n + 3)
@@ -194,9 +223,10 @@ function AkimaInterpolation(u, t; extrapolate = false, cache_parameters = false)
     c = (3.0 .* m[3:(end - 2)] .- 2.0 .* b[1:(end - 1)] .- b[2:end]) ./ dt
     d = (b[1:(end - 1)] .+ b[2:end] .- 2.0 .* m[3:(end - 2)]) ./ dt .^ 2
 
-    A = AkimaInterpolation(u, t, nothing, b, c, d, extrapolate, cache_parameters)
+    A = AkimaInterpolation(
+        u, t, nothing, b, c, d, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    AkimaInterpolation(u, t, I, b, c, d, extrapolate, cache_parameters)
+    AkimaInterpolation(u, t, I, b, c, d, extrapolate, cache_parameters, linear_lookup)
 end
 
 """
@@ -216,6 +246,10 @@ Extrapolation extends the last constant polynomial at the end points on each sid
   - `dir`: indicates which value should be used for interpolation (`:left` or `:right`).
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     u::uType
@@ -226,18 +260,23 @@ struct ConstantInterpolation{uType, tType, IType, T} <: AbstractInterpolation{T}
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    function ConstantInterpolation(u, t, I, dir, extrapolate, cache_parameters)
+    linear_lookup::Bool
+    function ConstantInterpolation(
+            u, t, I, dir, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), eltype(u)}(
-            u, t, I, nothing, dir, extrapolate, Ref(1), cache_parameters)
+            u, t, I, nothing, dir, extrapolate, Ref(1), cache_parameters, linear_lookup)
     end
 end
 
 function ConstantInterpolation(
-        u, t; dir = :left, extrapolate = false, cache_parameters = false)
+        u, t; dir = :left, extrapolate = false,
+        cache_parameters = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
-    A = ConstantInterpolation(u, t, nothing, dir, extrapolate, cache_parameters)
+    A = ConstantInterpolation(
+        u, t, nothing, dir, extrapolate, cache_parameters, assume_linear_t)
     I = cumulative_integral(A, cache_parameters)
-    ConstantInterpolation(u, t, I, dir, extrapolate, cache_parameters)
+    ConstantInterpolation(u, t, I, dir, extrapolate, cache_parameters, assume_linear_t)
 end
 
 """
@@ -255,6 +294,10 @@ Extrapolation extends the last quadratic polynomial on each side.
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
        AbstractInterpolation{T}
@@ -268,7 +311,10 @@ struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    function QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters)
+    linear_lookup::Bool
+    function QuadraticSpline(
+            u, t, I, p, tA, d, z, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.σ), typeof(tA),
             typeof(d), typeof(z), eltype(u)}(u,
             t,
@@ -279,15 +325,18 @@ struct QuadraticSpline{uType, tType, IType, pType, tAType, dType, zType, T} <:
             z,
             extrapolate,
             Ref(1),
-            cache_parameters
+            cache_parameters,
+            linear_lookup
         )
     end
 end
 
 function QuadraticSpline(
         u::uType, t; extrapolate = false,
-        cache_parameters = false) where {uType <: AbstractVector{<:Number}}
+        cache_parameters = false, assume_linear_t = 1e-2) where {uType <:
+                                                                 AbstractVector{<:Number}}
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     s = length(t)
     dl = ones(eltype(t), s - 1)
     d_tmp = ones(eltype(t), s)
@@ -301,15 +350,18 @@ function QuadraticSpline(
     z = tA \ d
 
     p = QuadraticSplineParameterCache(z, t, cache_parameters)
-    A = QuadraticSpline(u, t, nothing, p, tA, d, z, extrapolate, cache_parameters)
+    A = QuadraticSpline(
+        u, t, nothing, p, tA, d, z, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters)
+    QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters, linear_lookup)
 end
 
 function QuadraticSpline(
-        u::uType, t; extrapolate = false, cache_parameters = false) where {uType <:
-                                                                           AbstractVector}
+        u::uType, t; extrapolate = false, cache_parameters = false,
+        assume_linear_t = 1e-2) where {uType <:
+                                       AbstractVector}
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     s = length(t)
     dl = ones(eltype(t), s - 1)
     d_tmp = ones(eltype(t), s)
@@ -324,9 +376,10 @@ function QuadraticSpline(
     z = [z_s for z_s in eachslice(z_, dims = ndims(z_))]
 
     p = QuadraticSplineParameterCache(z, t, cache_parameters)
-    A = QuadraticSpline(u, t, nothing, p, tA, d, z, extrapolate, cache_parameters)
+    A = QuadraticSpline(
+        u, t, nothing, p, tA, d, z, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters)
+    QuadraticSpline(u, t, I, p, tA, d, z, extrapolate, cache_parameters, linear_lookup)
 end
 
 """
@@ -344,6 +397,10 @@ Second derivative on both ends are zero, which are also called "natural" boundar
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInterpolation{T}
     u::uType
@@ -355,7 +412,9 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInter
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    function CubicSpline(u, t, I, p, h, z, extrapolate, cache_parameters)
+    linear_lookup::Bool
+    function CubicSpline(u, t, I, p, h, z, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(p.c₁), typeof(h), typeof(z), eltype(u)}(
             u,
             t,
@@ -365,15 +424,17 @@ struct CubicSpline{uType, tType, IType, pType, hType, zType, T} <: AbstractInter
             z,
             extrapolate,
             Ref(1),
-            cache_parameters
+            cache_parameters,
+            linear_lookup
         )
     end
 end
 
 function CubicSpline(u::uType,
         t;
-        extrapolate = false, cache_parameters = false) where {uType <:
-                                                              AbstractVector{<:Number}}
+        extrapolate = false, cache_parameters = false,
+        assume_linear_t = 1e-2) where {uType <:
+                                       AbstractVector{<:Number}}
     u, t = munge_data(u, t)
     n = length(t) - 1
     h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0)
@@ -391,16 +452,18 @@ function CubicSpline(u::uType,
              6(u[i + 1] - u[i]) / h[i + 1] - 6(u[i] - u[i - 1]) / h[i],
         1:(n + 1))
     z = tA \ d
-
+    linear_lookup = seems_linear(assume_linear_t, t)
     p = CubicSplineParameterCache(u, h, z, cache_parameters)
-    A = CubicSpline(u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters)
+    A = CubicSpline(
+        u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters)
+    CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, linear_lookup)
 end
 
 function CubicSpline(
-        u::uType, t; extrapolate = false, cache_parameters = false) where {uType <:
-                                                                           AbstractVector}
+        u::uType, t; extrapolate = false, cache_parameters = false,
+        assume_linear_t = 1e-2) where {uType <:
+                                       AbstractVector}
     u, t = munge_data(u, t)
     n = length(t) - 1
     h = vcat(0, map(k -> t[k + 1] - t[k], 1:(length(t) - 1)), 0)
@@ -417,9 +480,10 @@ function CubicSpline(
     z = [z_s for z_s in eachslice(z_, dims = ndims(z_))]
 
     p = CubicSplineParameterCache(u, h, z, cache_parameters)
-    A = CubicSpline(u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters)
+    A = CubicSpline(
+        u, t, nothing, p, h[1:(n + 1)], z, extrapolate, cache_parameters, assume_linear_t)
     I = cumulative_integral(A, cache_parameters)
-    CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters)
+    CubicSpline(u, t, I, p, h[1:(n + 1)], z, extrapolate, cache_parameters, assume_linear_t)
 end
 
 """
@@ -439,6 +503,10 @@ Extrapolation is a constant polynomial of the end points on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
        AbstractInterpolation{T}
@@ -453,6 +521,7 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
     knotVecType::Symbol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
+    linear_lookup::Bool
     function BSplineInterpolation(u,
             t,
             d,
@@ -462,7 +531,9 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
             N,
             pVecType,
             knotVecType,
-            extrapolate)
+            extrapolate,
+            assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u,
             t,
             d,
@@ -473,13 +544,14 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, NType, T} <:
             pVecType,
             knotVecType,
             extrapolate,
-            Ref(1)
+            Ref(1),
+            linear_lookup
         )
     end
 end
 
 function BSplineInterpolation(
-        u, t, d, pVecType, knotVecType; extrapolate = false)
+        u, t, d, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
     n = length(t)
     n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d+1) points.")
@@ -544,7 +616,7 @@ function BSplineInterpolation(
     c = vec(N \ u[:, :])
     N = zeros(eltype(t), n)
     BSplineInterpolation(
-        u, t, d, p, k, c, N, pVecType, knotVecType, extrapolate)
+        u, t, d, p, k, c, N, pVecType, knotVecType, extrapolate, assume_linear_t)
 end
 
 """
@@ -566,6 +638,10 @@ Extrapolation is a constant polynomial of the end points on each side.
 ## Keyword Arguments
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
        AbstractInterpolation{T}
@@ -581,6 +657,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
     knotVecType::Symbol
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
+    linear_lookup::Bool
     function BSplineApprox(u,
             t,
             d,
@@ -591,8 +668,10 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
             N,
             pVecType,
             knotVecType,
-            extrapolate
+            extrapolate,
+            assume_linear_t
     )
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(N), eltype(u)}(u,
             t,
             d,
@@ -604,13 +683,14 @@ struct BSplineApprox{uType, tType, pType, kType, cType, NType, T} <:
             pVecType,
             knotVecType,
             extrapolate,
-            Ref(1)
+            Ref(1),
+            linear_lookup
         )
     end
 end
 
 function BSplineApprox(
-        u, t, d, h, pVecType, knotVecType; extrapolate = false)
+        u, t, d, h, pVecType, knotVecType; extrapolate = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
     n = length(t)
     h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d+1) control points.")
@@ -696,7 +776,7 @@ function BSplineApprox(
     c[2:(end - 1)] .= vec(P)
     N = zeros(eltype(t), h)
     BSplineApprox(
-        u, t, d, h, p, k, c, N, pVecType, knotVecType, extrapolate)
+        u, t, d, h, p, k, c, N, pVecType, knotVecType, extrapolate, assume_linear_t)
 end
 
 """
@@ -714,6 +794,10 @@ It is a Cubic Hermite interpolation, which is a piece-wise third degree polynomi
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInterpolation{T}
     du::duType
@@ -724,19 +808,25 @@ struct CubicHermiteSpline{uType, tType, IType, duType, pType, T} <: AbstractInte
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    function CubicHermiteSpline(du, u, t, I, p, extrapolate, cache_parameters)
+    linear_lookup::Bool
+    function CubicHermiteSpline(
+            du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(du), typeof(p.c₁), eltype(u)}(
-            du, u, t, I, p, extrapolate, Ref(1), cache_parameters)
+            du, u, t, I, p, extrapolate, Ref(1), cache_parameters, linear_lookup)
     end
 end
 
-function CubicHermiteSpline(du, u, t; extrapolate = false, cache_parameters = false)
+function CubicHermiteSpline(
+        du, u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
     @assert length(u)==length(du) "Length of `u` is not equal to length of `du`."
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     p = CubicHermiteParameterCache(du, u, t, cache_parameters)
-    A = CubicHermiteSpline(du, u, t, nothing, p, extrapolate, cache_parameters)
+    A = CubicHermiteSpline(
+        du, u, t, nothing, p, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    CubicHermiteSpline(du, u, t, I, p, extrapolate, cache_parameters)
+    CubicHermiteSpline(du, u, t, I, p, extrapolate, cache_parameters, linear_lookup)
 end
 
 """
@@ -755,11 +845,16 @@ section 3.4 for more details.
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
-function PCHIPInterpolation(u, t; extrapolate = false, cache_parameters = false)
+function PCHIPInterpolation(
+        u, t; extrapolate = false, cache_parameters = false, assume_linear_t = 1e-2)
     u, t = munge_data(u, t)
     du = du_PCHIP(u, t)
-    CubicHermiteSpline(du, u, t; extrapolate, cache_parameters)
+    CubicHermiteSpline(du, u, t; extrapolate, cache_parameters, assume_linear_t)
 end
 
 """
@@ -778,6 +873,10 @@ It is a Quintic Hermite interpolation, which is a piece-wise fifth degree polyno
 
   - `extrapolate`: boolean value to allow extrapolation. Defaults to `false`.
   - `cache_parameters`: precompute parameters at initialization for faster interpolation computations. Note: if activated, `u` and `t` should not be modified. Defaults to `false`.
+  - `assume_linear_t`: boolean value to specify a faster index lookup behaviour for
+    evenly-distributed abscissae. Alternatively, a numerical threshold may be specified
+    for a test based on the normalized standard deviation of the difference with respect
+    to the straight line (see [`looks_linear`](@ref)). Defaults to 1e-2.
 """
 struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
        AbstractInterpolation{T}
@@ -790,18 +889,25 @@ struct QuinticHermiteSpline{uType, tType, IType, duType, dduType, pType, T} <:
     extrapolate::Bool
     idx_prev::Base.RefValue{Int}
     cache_parameters::Bool
-    function QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate, cache_parameters)
+    linear_lookup::Bool
+    function QuinticHermiteSpline(
+            ddu, du, u, t, I, p, extrapolate, cache_parameters, assume_linear_t)
+        linear_lookup = seems_linear(assume_linear_t, t)
         new{typeof(u), typeof(t), typeof(I), typeof(du),
             typeof(ddu), typeof(p.c₁), eltype(u)}(
-            ddu, du, u, t, I, p, extrapolate, Ref(1), cache_parameters)
+            ddu, du, u, t, I, p, extrapolate, Ref(1), cache_parameters, linear_lookup)
     end
 end
 
-function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false, cache_parameters = false)
+function QuinticHermiteSpline(ddu, du, u, t; extrapolate = false,
+        cache_parameters = false, assume_linear_t = 1e-2)
     @assert length(u)==length(du)==length(ddu) "Length of `u` is not equal to length of `du` or `ddu`."
     u, t = munge_data(u, t)
+    linear_lookup = seems_linear(assume_linear_t, t)
     p = QuinticHermiteParameterCache(ddu, du, u, t, cache_parameters)
-    A = QuinticHermiteSpline(ddu, du, u, t, nothing, p, extrapolate, cache_parameters)
+    A = QuinticHermiteSpline(
+        ddu, du, u, t, nothing, p, extrapolate, cache_parameters, linear_lookup)
     I = cumulative_integral(A, cache_parameters)
-    QuinticHermiteSpline(ddu, du, u, t, I, p, extrapolate, cache_parameters)
+    QuinticHermiteSpline(
+        ddu, du, u, t, I, p, extrapolate, cache_parameters, linear_lookup)
 end
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index c409d5fc..69a5cf94 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -1,7 +1,8 @@
 function _interpolate(A, t)
     ((t < A.t[1] || t > A.t[end]) && !A.extrapolate) &&
         throw(ExtrapolationError())
-    val, idx_prev = _interpolate(A, t, A.idx_prev[])
+    idx_guess = A.idx_prev[]
+    val, idx_prev = _interpolate(A, t, idx_guess)
     A.idx_prev[] = idx_prev
     return val
 end
@@ -15,7 +16,7 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues
         u1 = u2 = one(eltype(A.u))
         slope = t * get_parameters(A, idx)
     else
-        idx = get_idx(A.t, t, iguess)
+        idx = get_idx(A, t, iguess)
         t1, t2 = A.t[idx], A.t[idx + 1]
         u1, u2 = A.u[idx], A.u[idx + 1]
         slope = get_parameters(A, idx)
@@ -36,7 +37,7 @@ function _interpolate(A::LinearInterpolation{<:AbstractVector}, t::Number, igues
 end
 
 function _interpolate(A::LinearInterpolation{<:AbstractMatrix}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt = t - A.t[idx]
     slope = get_parameters(A, idx)
     return A.u[:, idx] + slope * Δt, idx
@@ -45,7 +46,7 @@ end
 # Quadratic Interpolation
 _quad_interp_indices(A, t) = _quad_interp_indices(A, t, firstindex(A.t) - 1)
 function _quad_interp_indices(A::QuadraticInterpolation, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess; idx_shift = A.mode == :Backward ? -1 : 0, ub_shift = -2)
+    idx = get_idx(A, t, iguess; idx_shift = A.mode == :Backward ? -1 : 0, ub_shift = -2)
     idx, idx + 1, idx + 2
 end
 
@@ -60,7 +61,7 @@ end
 
 # Lagrange Interpolation
 function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     findRequiredIdxs!(A, t, idx)
     if A.t[A.idxs[1]] == t
         return A.u[A.idxs[1]], idx
@@ -89,7 +90,7 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractVector}, t::Number, igu
 end
 
 function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     findRequiredIdxs!(A, t, idx)
     if A.t[A.idxs[1]] == t
         return A.u[:, A.idxs[1]], idx
@@ -118,7 +119,7 @@ function _interpolate(A::LagrangeInterpolation{<:AbstractMatrix}, t::Number, igu
 end
 
 function _interpolate(A::AkimaInterpolation{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     wj = t - A.t[idx]
     (@evalpoly wj A.u[idx] A.b[idx] A.c[idx] A.d[idx]), idx
 end
@@ -127,10 +128,10 @@ end
 function _interpolate(A::ConstantInterpolation{<:AbstractVector}, t::Number, iguess)
     if A.dir === :left
         # :left means that value to the left is used for interpolation
-        idx = get_idx(A.t, t, iguess; lb = 1, ub_shift = 0)
+        idx = get_idx(A, t, iguess; lb = 1, ub_shift = 0)
     else
         # :right means that value to the right is used for interpolation
-        idx = get_idx(A.t, t, iguess; side = :first, lb = 1, ub_shift = 0)
+        idx = get_idx(A, t, iguess; side = :first, lb = 1, ub_shift = 0)
     end
     A.u[idx], idx
 end
@@ -138,17 +139,17 @@ end
 function _interpolate(A::ConstantInterpolation{<:AbstractMatrix}, t::Number, iguess)
     if A.dir === :left
         # :left means that value to the left is used for interpolation
-        idx = get_idx(A.t, t, iguess; lb = 1, ub_shift = 0)
+        idx = get_idx(A, t, iguess; lb = 1, ub_shift = 0)
     else
         # :right means that value to the right is used for interpolation
-        idx = get_idx(A.t, t, iguess; side = :first, lb = 1, ub_shift = 0)
+        idx = get_idx(A, t, iguess; side = :first, lb = 1, ub_shift = 0)
     end
     A.u[:, idx], idx
 end
 
 # QuadraticSpline Interpolation
 function _interpolate(A::QuadraticSpline{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Cᵢ = A.u[idx]
     Δt = t - A.t[idx]
     σ = get_parameters(A, idx)
@@ -157,7 +158,7 @@ end
 
 # CubicSpline Interpolation
 function _interpolate(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₁ = t - A.t[idx]
     Δt₂ = A.t[idx + 1] - t
     I = (A.z[idx] * Δt₂^3 + A.z[idx + 1] * Δt₁^3) / (6A.h[idx + 1])
@@ -174,7 +175,7 @@ function _interpolate(A::BSplineInterpolation{<:AbstractVector{<:Number}},
     t < A.t[1] && return A.u[1], 1
     t > A.t[end] && return A.u[end], lastindex(t)
     # change t into param [0 1]
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
     n = length(A.t)
     N = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.N
@@ -191,7 +192,7 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
     t < A.t[1] && return A.u[1], 1
     t > A.t[end] && return A.u[end], lastindex(t)
     # change t into param [0 1]
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     t = A.p[idx] + (t - A.t[idx]) / (A.t[idx + 1] - A.t[idx]) * (A.p[idx + 1] - A.p[idx])
     N = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.N
     nonzero_coefficient_idxs = spline_coefficients!(N, A.d, A.k, t)
@@ -205,7 +206,7 @@ end
 # Cubic Hermite Spline
 function _interpolate(
         A::CubicHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₀ = t - A.t[idx]
     Δt₁ = t - A.t[idx + 1]
     out = A.u[idx] + Δt₀ * A.du[idx]
@@ -217,7 +218,7 @@ end
 # Quintic Hermite Spline
 function _interpolate(
         A::QuinticHermiteSpline{<:AbstractVector{<:Number}}, t::Number, iguess)
-    idx = get_idx(A.t, t, iguess)
+    idx = get_idx(A, t, iguess)
     Δt₀ = t - A.t[idx]
     Δt₁ = t - A.t[idx + 1]
     out = A.u[idx] + Δt₀ * (A.du[idx] + A.ddu[idx] * Δt₀ / 2)
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index 2a2392bd..ccd3db63 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -94,7 +94,40 @@ function munge_data(U::StridedMatrix, t::AbstractVector)
     return U, t
 end
 
-function get_idx(tvec, t, iguess; lb = 1, ub_shift = -1, idx_shift = 0, side = :last)
+seems_linear(assume_linear_t::Bool, _) = assume_linear_t
+seems_linear(assume_linear_t::Number, t) = looks_linear(t; threshold = assume_linear_t)
+
+"""
+    looks_linear(t; threshold = 1e-2)
+
+Determine if the abscissae `t` are regularly distributed, taking the standard deviation of
+the difference between the array of abscissae with respect to the straight line linking
+its first and last elements, normalized by the range of `t`. If this standard deviation is
+below the given `threshold`, the vector looks linear (return true). Internal function -
+interface may change.
+"""
+function looks_linear(t; threshold = 1e-2)
+    length(t) <= 2 && return true
+    t_0, t_f = first(t), last(t)
+    t_span = t_f - t_0
+    tspan_over_N = t_span * length(t)^(-1)
+    norm_var = sum(
+        (t_i - t_0 - i * tspan_over_N)^2 for (i, t_i) in enumerate(t)
+    ) / (length(t) * t_span^2)
+    norm_var < threshold^2
+end
+
+function get_idx(A::AbstractInterpolation, t, iguess; lb = 1,
+        ub_shift = -1, idx_shift = 0, side = :last)
+    iguess = if hasfield(typeof(A), :linear_lookup) &&
+                A.linear_lookup
+        f = (t - first(A.t)) / (last(A.t) - first(A.t))
+        i_0, i_f = firstindex(A.t), lastindex(A.t)
+        round(typeof(firstindex(A.t)), f * (i_f - i_0) + i_0)
+    else
+        iguess
+    end
+    tvec = A.t
     ub = length(tvec) + ub_shift
     return if side == :last
         clamp(searchsortedlastcorrelated(tvec, t, iguess) + idx_shift, lb, ub)
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index 53fff25e..1d4e2e6b 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -2,6 +2,7 @@ using DataInterpolations
 using FindFirstFunctions: searchsortedfirstcorrelated
 using StableRNGs
 using Optim, ForwardDiff
+using BenchmarkTools
 
 function test_interpolation_type(T)
     @test T <: DataInterpolations.AbstractInterpolation
@@ -800,3 +801,20 @@ f_cubic_spline = c -> square(CubicSpline, c)
 @test ForwardDiff.derivative(f_quadratic_spline, 4.0) ≈ 8.0
 @test ForwardDiff.derivative(f_cubic_spline, 2.0) ≈ 4.0
 @test ForwardDiff.derivative(f_cubic_spline, 4.0) ≈ 8.0
+
+@testset "Linear lookup" begin
+    for N in (100, 1_000, 10_000, 100_000, 1_000_000) # Interpolant size
+        seed = 1234
+        t = collect(LinRange(0, 1, N)) # collect to avoid fast LinRange dispatch
+        y = rand(N)
+        A = LinearInterpolation(y, t)
+        A_fallback = LinearInterpolation(copy(y), copy(t); assume_linear_t = false)
+        @test A.linear_lookup
+        @test !(A_fallback.linear_lookup)
+
+        n_samples = 1_000
+        b_linear = @benchmark $A(rand()) samples=n_samples
+        b_fallback = @benchmark $A_fallback(rand()) samples=n_samples
+        @test mean(b_linear.times) < mean(b_fallback.times)
+    end
+end