diff --git a/README.md b/README.md
index d31ce669..e168f60b 100644
--- a/README.md
+++ b/README.md
@@ -55,12 +55,11 @@ corresponding to `(u,t)` pairs.
   - `QuadraticSpline(u,t)` - A quadratic spline interpolation.
   - `CubicSpline(u,t)` - A cubic spline interpolation.
   - `AkimaInterpolation(u, t)` - Akima spline interpolation provides a smoothing effect and is computationally efficient.
-  - `BSplineInterpolation(u,t,d,pVec,knotVec)` - An interpolation B-spline. This is a B-spline which hits each of the data points. The argument choices are:
+  - `BSplineInterpolation(u,t,d,knotVec)` - An interpolation B-spline. This is a B-spline which hits each of the data points. The argument choices are:
     
       + `d` - degree of B-spline
-      + `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniform spaced parameters and `pVec = :ArcLen` for parameters generated by chord length method.
       + `knotVec` - Symbol to Knot Vector, `knotVec = :Uniform` for uniform knot vector, `knotVec = :Average` for average spaced knot vector.
-  - `BSplineApprox(u,t,d,h,pVec,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
+  - `BSplineApprox(u,t,d,h,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
   - `CubicHermiteSpline(du, u, t)` - A third order Hermite interpolation, which matches the values and first (`du`) order derivatives in the data points exactly.
   - `PCHIPInterpolation(u, t)` - a type of `CubicHermiteSpline` where the derivative values `du` are derived from the input data in such a way that the interpolation never overshoots the data.
   - `QuinticHermiteSpline(ddu, du, u, t)` - A fifth order Hermite interpolation, which matches the values and first (`du`) and second (`ddu`) order derivatives in the data points exactly.
diff --git a/docs/src/clarification.md b/docs/src/clarification.md
index d0e8501d..5d8c0e24 100644
--- a/docs/src/clarification.md
+++ b/docs/src/clarification.md
@@ -9,7 +9,7 @@ using DataInterpolations, Plots
 using DataInterpolations: derivative
 u = [2.0, 1.0, 5.0, 4.0, 5.0, 4.0]
 t = [0.0, 2.0, 3.5, 4.0, 5.0, 6.5]
-bspline = BSplineInterpolation(u, t, 0, :Uniform, :Uniform; extrapolation_left=ExtrapolationType.Extension,  extrapolation_right=ExtrapolationType.Extension)
+bspline = BSplineInterpolation(u, t, 0, :Uniform; extrapolation_left=ExtrapolationType.Extension,  extrapolation_right=ExtrapolationType.Extension)
 plot(bspline)
 ```
 
@@ -19,7 +19,7 @@ Thus, the plot for B-Spline interpolation does not appear the same as the plot f
 # Derivative behavior of quadratic B-Spline 
 
 ```@example interpclarity
-bspline = BSplineInterpolation(u, t, 2, :Uniform, :Uniform; extrapolation_left=ExtrapolationType.Extension,  extrapolation_right=ExtrapolationType.Extension)
+bspline = BSplineInterpolation(u, t, 2, :Uniform; extrapolation_left=ExtrapolationType.Extension,  extrapolation_right=ExtrapolationType.Extension)
 plot(t->derivative(bspline, t))
 ```
 
diff --git a/docs/src/index.md b/docs/src/index.md
index a6ce9058..9447bc83 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -25,12 +25,11 @@ corresponding to `(u,t)` pairs.
   - `QuadraticSpline(u,t)` - A quadratic spline interpolation.
   - `CubicSpline(u,t)` - A cubic spline interpolation.
   - `AkimaInterpolation(u, t)` - Akima spline interpolation provides a smoothing effect and is computationally efficient.
-  - `BSplineInterpolation(u,t,d,pVec,knotVec)` - An interpolation B-spline. This is a B-spline that hits each of the data points. The argument choices are:
+  - `BSplineInterpolation(u,t,d,knotVec)` - An interpolation B-spline. This is a B-spline that hits each of the data points. The argument choices are:
     
       + `d` - degree of B-spline
-      + `pVec` - Symbol to Parameters Vector, `pVec = :Uniform` for uniformly spaced parameters, and `pVec = :ArcLen` for parameters generated by the chord length method.
       + `knotVec` - Symbol to Knot Vector, `knotVec = :Uniform` for uniform knot vector, `knotVec = :Average` for average spaced knot vector.
-  - `BSplineApprox(u,t,d,h,pVec,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
+  - `BSplineApprox(u,t,d,h,knotVec)` - A regression B-spline which smooths the fitting curve. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h<length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
   - `CubicHermiteSpline(du, u, t)` - A third order Hermite interpolation, which matches the values and first (`du`) order derivatives in the data points exactly.
   - `PCHIPInterpolation(u, t)` - a type of `CubicHermiteSpline` where the derivative values `du` are derived from the input data in such a way that the interpolation never overshoots the data.
   - `QuinticHermiteSpline(ddu, du, u, t)` - a fifth order Hermite interpolation, which matches the values and first (`du`) and second (`ddu`) order derivatives in the data points exactly.
diff --git a/docs/src/methods.md b/docs/src/methods.md
index 6a9a8432..b6575d3b 100644
--- a/docs/src/methods.md
+++ b/docs/src/methods.md
@@ -124,10 +124,10 @@ that every data point is taken into account for each point of the curve.
 The interpolating B-spline is the version which hits each of the points. This
 method is described in more detail [here](https://pages.mtu.edu/%7Eshene/COURSES/cs3621/NOTES/INT-APP/CURVE-INT-global.html).
 Let's plot a cubic B-spline (3rd order). Since the data points are not close to
-uniformly spaced, we will use the `:ArcLen` and `:Average` choices:
+uniformly spaced, we will use the `:Average` knot vector:
 
 ```@example tutorial
-A = BSplineInterpolation(u, t, 3, :ArcLen, :Average)
+A = BSplineInterpolation(u, t, 3, :Average)
 plot(A)
 ```
 
@@ -138,7 +138,7 @@ is a least square approximation. This has a natural effect of smoothing the
 data. For example, if we use 4 control points, we get the result:
 
 ```@example tutorial
-A = BSplineApprox(u, t, 3, 4, :ArcLen, :Average)
+A = BSplineApprox(u, t, 3, 4, :Average)
 plot(A)
 ```
 
diff --git a/src/derivatives.jl b/src/derivatives.jl
index 268e3ec0..c8c9b70b 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -249,94 +249,102 @@ function _derivative(A::CubicSpline{<:AbstractVector}, t::Number, iguess)
 end
 
 function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Number, iguess)
-    # change t into param [0 1]
-    t < A.t[1] && return zero(A.u[1])
-    t > A.t[end] && return zero(A.u[end])
-    idx = get_idx(A, t, iguess)
+    if A.d == 0
+        return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : typed_nan(A.u)
+    end
     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
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
-    spline_coefficients!(sc, A.d - 1, A.k, t_)
+    spline_coefficients!(sc, A.d - 1, A.k, t)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
-        ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2])
+        denom = A.k[A.d + 2] - A.k[2]
+        ducum = denom != 0 ? (A.c[2] - A.c[1]) / denom : zero(eltype(A.u))
     else
         for i in 1:(n - 1)
-            ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
+            denom = A.k[i + A.d + 1] - A.k[i + 1]
+            if denom != 0
+                ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / denom
+            end
         end
     end
-    return ducum * A.d * scale
+    return ducum * A.d
 end
 
 function _derivative(
         A::BSplineInterpolation{<:AbstractArray{<:Number}}, t::Number, iguess
     )
-    # change t into param [0 1]
+    if A.d == 0
+        return isempty(searchsorted(A.t, t)) ? zero(A.u[:, 1]) :
+            typed_nan(A.u) .* A.u[:, 1]
+    end
     ax_u = axes(A.u)[1:(end - 1)]
-    t < A.t[1] && return zeros(size(A.u)[1:(end - 1)]...)
-    t > A.t[end] && return zeros(size(A.u)[1:(end - 1)]...)
-    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
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
-    spline_coefficients!(sc, A.d - 1, A.k, t_)
+    spline_coefficients!(sc, A.d - 1, A.k, t)
     ducum = zeros(size(A.u)[1:(end - 1)]...)
     if t == A.t[1]
-        ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2])
+        denom = A.k[A.d + 2] - A.k[2]
+        if denom != 0
+            ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / denom
+        end
     else
         for i in 1:(n - 1)
-            ducum = ducum +
-                sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
-                (A.k[i + A.d + 1] - A.k[i + 1])
+            denom = A.k[i + A.d + 1] - A.k[i + 1]
+            if denom != 0
+                ducum = ducum +
+                    sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) / denom
+            end
         end
     end
-    return ducum * A.d * scale
+    return ducum * A.d
 end
 # BSpline Curve Approx
 function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, iguess)
-    # change t into param [0 1]
-    t < A.t[1] && return zero(A.u[1])
-    t > A.t[end] && return zero(A.u[end])
-    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
+    if A.d == 0
+        return isempty(searchsorted(A.t, t)) ? zero(A.u[1]) : typed_nan(A.u)
+    end
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
-    spline_coefficients!(sc, A.d - 1, A.k, t_)
+    spline_coefficients!(sc, A.d - 1, A.k, t)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
-        ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2])
+        denom = A.k[A.d + 2] - A.k[2]
+        ducum = denom != 0 ? (A.c[2] - A.c[1]) / denom : zero(eltype(A.u))
     else
         for i in 1:(A.h - 1)
-            ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
+            denom = A.k[i + A.d + 1] - A.k[i + 1]
+            if denom != 0
+                ducum += sc[i + 1] * (A.c[i + 1] - A.c[i]) / denom
+            end
         end
     end
-    return ducum * A.d * scale
+    return ducum * A.d
 end
 
 function _derivative(
         A::BSplineApprox{<:AbstractArray{<:Number}}, t::Number, iguess
     )
-    # change t into param [0 1]
+    if A.d == 0
+        return isempty(searchsorted(A.t, t)) ? zero(A.u[:, 1]) :
+            typed_nan(A.u) .* A.u[:, 1]
+    end
     ax_u = axes(A.u)[1:(end - 1)]
-    t < A.t[1] && return zeros(size(A.u)[1:(end - 1)]...)
-    t > A.t[end] && return zeros(size(A.u)[1:(end - 1)]...)
-    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
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
-    spline_coefficients!(sc, A.d - 1, A.k, t_)
+    spline_coefficients!(sc, A.d - 1, A.k, t)
     ducum = zeros(size(A.u)[1:(end - 1)]...)
     if t == A.t[1]
-        ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / (A.k[A.d + 2])
+        denom = A.k[A.d + 2] - A.k[2]
+        if denom != 0
+            ducum = (A.c[ax_u..., 2] - A.c[ax_u..., 1]) / denom
+        end
     else
         for i in 1:(A.h - 1)
-            ducum += sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
-                (A.k[i + A.d + 1] - A.k[i + 1])
+            denom = A.k[i + A.d + 1] - A.k[i + 1]
+            if denom != 0
+                ducum += sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) / denom
+            end
         end
     end
-    return ducum * A.d * scale
+    return ducum * A.d
 end
 # Cubic Hermite Spline
 function _derivative(
diff --git a/src/integrals.jl b/src/integrals.jl
index 92255892..28734339 100644
--- a/src/integrals.jl
+++ b/src/integrals.jl
@@ -287,11 +287,86 @@ end
 function _integral(A::LagrangeInterpolation, idx::Number, t1::Number, t2::Number)
     throw(IntegralNotFoundError())
 end
-function _integral(A::BSplineInterpolation, idx::Number, t1::Number, t2::Number)
-    throw(IntegralNotFoundError())
+# Evaluate the antiderivative of a B-spline at a point.
+# The antiderivative of a degree-d B-spline is a degree-(d+1) B-spline with
+# extended knot vector and coefficients derived from the original.
+function _bspline_antiderivative_val(c, d, k, t_eval)
+    nc = length(c)
+    dp1 = d + 1
+    # Antiderivative coefficients: C[1] = 0, C[i+1] = C[i] + c[i]*(k[i+d+1]-k[i])/(d+1)
+    T = promote_type(eltype(c), eltype(k), typeof(t_eval))
+    C = zeros(T, nc + 1)
+    for i in 1:nc
+        C[i + 1] = C[i] + c[i] * (k[i + dp1] - k[i]) / dp1
+    end
+    # Extended knot vector: prepend k[1], append k[end]
+    nk = length(k)
+    k_ext = zeros(eltype(k), nk + 2)
+    k_ext[1] = k[1]
+    for i in 1:nk
+        k_ext[i + 1] = k[i]
+    end
+    k_ext[end] = k[end]
+    # Evaluate degree-(d+1) B-spline with coefficients C on k_ext
+    sc = zeros(T, nc + 1)
+    nonzero = spline_coefficients!(sc, dp1, k_ext, t_eval)
+    result = zero(T)
+    for i in nonzero
+        result += sc[i] * C[i]
+    end
+    return result
 end
-function _integral(A::BSplineApprox, idx::Number, t1::Number, t2::Number)
-    throw(IntegralNotFoundError())
+
+function _integral(
+        A::BSplineInterpolation{<:AbstractVector{<:Number}}, idx::Number, t1::Number, t2::Number
+    )
+    return _bspline_antiderivative_val(A.c, A.d, A.k, t2) -
+        _bspline_antiderivative_val(A.c, A.d, A.k, t1)
+end
+function _integral(A::BSplineApprox{<:AbstractVector{<:Number}}, idx::Number, t1::Number, t2::Number)
+    return _bspline_antiderivative_val(A.c, A.d, A.k, t2) -
+        _bspline_antiderivative_val(A.c, A.d, A.k, t1)
+end
+
+# Override integral to bypass the hasfield(:I) check in the generic method.
+# The antiderivative is computed on the fly, so no cached I field is needed.
+const _BSplineTypes = Union{
+    BSplineInterpolation{<:AbstractVector{<:Number}},
+    BSplineApprox{<:AbstractVector{<:Number}},
+}
+
+function integral(A::_BSplineTypes, t::Number)
+    return integral(A, first(A.t), t)
+end
+
+function integral(A::_BSplineTypes, t1::Number, t2::Number)
+    t1 == t2 && return zero(eltype(A.u))
+    t1 > t2 && return -integral(A, t2, t1)
+
+    total = zero(eltype(A.u))
+    lo = t1
+    hi = t2
+
+    if lo < first(A.t)
+        if hi <= first(A.t)
+            return _extrapolate_integral_left(A, lo) - _extrapolate_integral_left(A, hi)
+        end
+        total += _extrapolate_integral_left(A, lo)
+        lo = first(A.t)
+    end
+
+    if hi > last(A.t)
+        if lo >= last(A.t)
+            return _extrapolate_integral_right(A, hi) - _extrapolate_integral_right(A, lo)
+        end
+        total += _extrapolate_integral_right(A, hi)
+        hi = last(A.t)
+    end
+
+    total += _bspline_antiderivative_val(A.c, A.d, A.k, hi) -
+        _bspline_antiderivative_val(A.c, A.d, A.k, lo)
+
+    return total
 end
 
 # Cubic Hermite Spline
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index 5c4a9f4a..b67bb361 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -814,10 +814,12 @@ function CubicSpline(
 end
 
 """
-    BSplineInterpolation(u, t, d, pVecType, knotVecType; extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
+    BSplineInterpolation(u, t, d, knotVecType; extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
 
-It is a curve defined by the linear combination of `n` basis functions of degree `d` where `n` is the number of data points. For more information, refer [https://pages.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve.html](https://pages.mtu.edu/%7Eshene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve.html).
+It is a B-spline interpolation of degree `d` that passes through all data points. The knot vector
+is constructed directly from the data sites `t`, and basis functions are evaluated directly at `t`
+values (no reparameterization). For more information, refer to de Boor's "A Practical Guide to Splines".
 Extrapolation is a constant polynomial of the end points on each side.
 
 ## Arguments
@@ -825,7 +827,6 @@ Extrapolation is a constant polynomial of the end points on each side.
   - `u`: data points.
   - `t`: time points.
   - `d`: degree of the piecewise polynomial.
-  - `pVecType`: symbol to parameters vector, `:Uniform` for uniform spaced parameters and `:ArcLen` for parameters generated by chord length method.
   - `knotVecType`: symbol to knot vector, `:Uniform` for uniform knot vector, `:Average` for average spaced knot vector.
 
 ## Keyword Arguments
@@ -842,16 +843,14 @@ Extrapolation is a constant polynomial of the end points on each side.
     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, scType, T} <:
+struct BSplineInterpolation{uType, tType, kType, cType, scType, T} <:
     AbstractInterpolation{T}
     u::uType
     t::tType
     d::Int    # degree
-    p::pType  # params vector
     k::kType  # knot vector
     c::cType  # control points
     sc::scType  # Spline coefficients (preallocated memory)
-    pVecType::Symbol
     knotVecType::Symbol
     extrapolation_left::ExtrapolationType.T
     extrapolation_right::ExtrapolationType.T
@@ -861,26 +860,22 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T} <:
             u,
             t,
             d,
-            p,
             k,
             c,
             sc,
-            pVecType,
             knotVecType,
             extrapolation_left,
             extrapolation_right,
             assume_linear_t
         )
         linear_lookup = seems_linear(assume_linear_t, t)
-        return new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u)}(
+        return new{typeof(u), typeof(t), typeof(k), typeof(c), typeof(sc), eltype(u)}(
             u,
             t,
             d,
-            p,
             k,
             c,
             sc,
-            pVecType,
             knotVecType,
             extrapolation_left,
             extrapolation_right,
@@ -891,7 +886,7 @@ struct BSplineInterpolation{uType, tType, pType, kType, cType, scType, T} <:
 end
 
 function BSplineInterpolation(
-        u::AbstractVector, t, d, pVecType, knotVecType;
+        u::AbstractVector, t, d, knotVecType;
         extrapolation::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, assume_linear_t = 1.0e-2
@@ -903,74 +898,38 @@ function BSplineInterpolation(
     u, t = munge_data(u, t)
     n = length(t)
     n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d + 1) points.")
-    s = zero(eltype(u))
-    p = zero(t)
     k = zeros(eltype(t), n + d + 1)
-    l = zeros(eltype(u), n - 1)
-    p[1] = zero(eltype(t))
-    p[end] = one(eltype(t))
-
-    for i in 2:n
-        s += hypot(t[i] - t[i - 1], u[i] - u[i - 1])
-        l[i - 1] = s
+    # Clamped knot vector endpoints
+    for i in 1:(d + 1)
+        k[i] = t[1]
+        k[n + i] = t[end]
     end
-    if pVecType == :Uniform
-        for i in 2:(n - 1)
-            p[i] = p[1] + (i - 1) * (p[end] - p[1]) / (n - 1)
-        end
-    elseif pVecType == :ArcLen
-        for i in 2:(n - 1)
-            p[i] = p[1] + l[i - 1] / s * (p[end] - p[1])
-        end
-    end
-
-    lidx = 1
-    ridx = length(k)
-    while lidx <= (d + 1) && ridx >= (length(k) - d)
-        k[lidx] = p[1]
-        k[ridx] = p[end]
-        lidx += 1
-        ridx -= 1
-    end
-
-    ps = zeros(eltype(t), n - 2)
-    s = zero(eltype(t))
-    for i in 2:(n - 1)
-        s += p[i]
-        ps[i - 1] = s
-    end
-
     if knotVecType == :Uniform
-        # uniformly spaced knot vector
-        # this method is not recommended because, if it is used with the chord length method for global interpolation,
-        # the system of linear equations would be singular.
+        # Uniformly spaced interior knots
         for i in (d + 2):n
-            k[i] = k[1] + (i - d - 1) // (n - d) * (k[end] - k[1])
+            k[i] = t[1] + (i - d - 1) // (n - d) * (t[end] - t[1])
         end
     elseif knotVecType == :Average
-        # average spaced knot vector
-        idx = 1
-        if d + 2 <= n
-            k[d + 2] = 1 // d * ps[d]
-        end
-        for i in (d + 3):n
-            k[i] = 1 // d * (ps[idx + d] - ps[idx])
-            idx += 1
+        # Average (de Boor) interior knots, generalized for all d >= 0
+        # using max(d, 1) to handle d = 0 (places knots at data sites)
+        denom = max(d, 1)
+        for j in 1:(n - d - 1)
+            k[d + 1 + j] = sum(t[(j + 1):(j + denom)]) / denom
         end
     end
     # control points
     sc = zeros(eltype(t), n, n)
-    spline_coefficients!(sc, d, k, p)
+    spline_coefficients!(sc, d, k, t)
     c = vec(sc \ u[:, :])
     sc = zeros(eltype(t), n)
     return BSplineInterpolation(
-        u, t, d, p, k, c, sc, pVecType, knotVecType,
+        u, t, d, k, c, sc, knotVecType,
         extrapolation_left, extrapolation_right, assume_linear_t
     )
 end
 
 function BSplineInterpolation(
-        u::AbstractArray, t, d, pVecType, knotVecType;
+        u::AbstractArray, t, d, knotVecType;
         extrapolation::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
@@ -983,81 +942,41 @@ function BSplineInterpolation(
     u, t = munge_data(u, t)
     n = length(t)
     n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d + 1) points.")
-    s = zero(eltype(u))
-    p = zero(t)
     k = zeros(eltype(t), n + d + 1)
-    l = zeros(eltype(u), n - 1)
-    p[1] = zero(eltype(t))
-    p[end] = one(eltype(t))
-
-    ax_u = axes(u)[1:(end - 1)]
-
-    for i in 2:n
-        s += √((t[i] - t[i - 1])^2 + sum((u[ax_u..., i] - u[ax_u..., i - 1]) .^ 2))
-        l[i - 1] = s
-    end
-    if pVecType == :Uniform
-        for i in 2:(n - 1)
-            p[i] = p[1] + (i - 1) * (p[end] - p[1]) / (n - 1)
-        end
-    elseif pVecType == :ArcLen
-        for i in 2:(n - 1)
-            p[i] = p[1] + l[i - 1] / s * (p[end] - p[1])
-        end
-    end
-
-    lidx = 1
-    ridx = length(k)
-    while lidx <= (d + 1) && ridx >= (length(k) - d)
-        k[lidx] = p[1]
-        k[ridx] = p[end]
-        lidx += 1
-        ridx -= 1
-    end
-
-    ps = zeros(eltype(t), n - 2)
-    s = zero(eltype(t))
-    for i in 2:(n - 1)
-        s += p[i]
-        ps[i - 1] = s
+    # Clamped knot vector endpoints
+    for i in 1:(d + 1)
+        k[i] = t[1]
+        k[n + i] = t[end]
     end
-
     if knotVecType == :Uniform
-        # uniformly spaced knot vector
-        # this method is not recommended because, if it is used with the chord length method for global interpolation,
-        # the system of linear equations would be singular.
+        # Uniformly spaced interior knots
         for i in (d + 2):n
-            k[i] = k[1] + (i - d - 1) // (n - d) * (k[end] - k[1])
+            k[i] = t[1] + (i - d - 1) // (n - d) * (t[end] - t[1])
         end
     elseif knotVecType == :Average
-        # average spaced knot vector
-        idx = 1
-        if d + 2 <= n
-            k[d + 2] = 1 // d * ps[d]
-        end
-        for i in (d + 3):n
-            k[i] = 1 // d * (ps[idx + d] - ps[idx])
-            idx += 1
+        # Average (de Boor) interior knots, generalized for all d >= 0
+        denom = max(d, 1)
+        for j in 1:(n - d - 1)
+            k[d + 1 + j] = sum(t[(j + 1):(j + denom)]) / denom
         end
     end
     # control points
     sc = zeros(eltype(t), n, n)
-    spline_coefficients!(sc, d, k, p)
+    spline_coefficients!(sc, d, k, t)
     c = (sc \ reshape(u, prod(size(u)[1:(end - 1)]), :)')'
     c = reshape(c, size(u)...)
     sc = zeros(eltype(t), n)
     return BSplineInterpolation(
-        u, t, d, p, k, c, sc, pVecType, knotVecType,
+        u, t, d, k, c, sc, knotVecType,
         extrapolation_left, extrapolation_right, assume_linear_t
     )
 end
 
 """
-    BSplineApprox(u, t, d, h, pVecType, knotVecType; extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
+    BSplineApprox(u, t, d, h, knotVecType; extrapolation::ExtrapolationType.T = ExtrapolationType.None, extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None)
 
 It is a regression based B-spline. The argument choices are the same as the `BSplineInterpolation`, with the additional parameter `h < length(t)` which is the number of control points to use, with smaller `h` indicating more smoothing.
-For more information, refer [http://www.cad.zju.edu.cn/home/zhx/GM/009/00-bsia.pdf](http://www.cad.zju.edu.cn/home/zhx/GM/009/00-bsia.pdf).
 Extrapolation is a constant polynomial of the end points on each side.
 
 ## Arguments
@@ -1066,7 +985,6 @@ Extrapolation is a constant polynomial of the end points on each side.
   - `t`: time points.
   - `d`: degree of the piecewise polynomial.
   - `h`: number of control points to use.
-  - `pVecType`: symbol to parameters vector, `:Uniform` for uniform spaced parameters and `:ArcLen` for parameters generated by chord length method.
   - `knotVecType`: symbol to knot vector, `:Uniform` for uniform knot vector, `:Average` for average spaced knot vector.
 
 ## Keyword Arguments
@@ -1083,17 +1001,15 @@ Extrapolation is a constant polynomial of the end points on each side.
     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, scType, T} <:
+struct BSplineApprox{uType, tType, kType, cType, scType, T} <:
     AbstractInterpolation{T}
     u::uType
     t::tType
     d::Int    # degree
     h::Int    # number of control points (n => h >= d >= 1)
-    p::pType  # params vector
     k::kType  # knot vector
     c::cType  # control points
     sc::scType  # Spline coefficients (preallocated memory)
-    pVecType::Symbol
     knotVecType::Symbol
     extrapolation_left::ExtrapolationType.T
     extrapolation_right::ExtrapolationType.T
@@ -1104,27 +1020,23 @@ struct BSplineApprox{uType, tType, pType, kType, cType, scType, T} <:
             t,
             d,
             h,
-            p,
             k,
             c,
             sc,
-            pVecType,
             knotVecType,
             extrapolation_left,
             extrapolation_right,
             assume_linear_t
         )
         linear_lookup = seems_linear(assume_linear_t, t)
-        return new{typeof(u), typeof(t), typeof(p), typeof(k), typeof(c), typeof(sc), eltype(u)}(
+        return new{typeof(u), typeof(t), typeof(k), typeof(c), typeof(sc), eltype(u)}(
             u,
             t,
             d,
             h,
-            p,
             k,
             c,
             sc,
-            pVecType,
             knotVecType,
             extrapolation_left,
             extrapolation_right,
@@ -1135,7 +1047,7 @@ struct BSplineApprox{uType, tType, pType, kType, cType, scType, T} <:
 end
 
 function BSplineApprox(
-        u::AbstractVector, t, d, h, pVecType, knotVecType;
+        u::AbstractVector, t, d, h, knotVecType;
         extrapolation::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None, assume_linear_t = 1.0e-2
@@ -1147,60 +1059,26 @@ function BSplineApprox(
     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.")
-    s = zero(eltype(u))
-    p = zero(t)
     k = zeros(eltype(t), h + d + 1)
-    l = zeros(eltype(u), n - 1)
-    p[1] = zero(eltype(t))
-    p[end] = one(eltype(t))
-
-    for i in 2:n
-        s += hypot(t[i] - t[i - 1], u[i] - u[i - 1])
-        l[i - 1] = s
-    end
-    if pVecType == :Uniform
-        for i in 2:(n - 1)
-            p[i] = p[1] + (i - 1) * (p[end] - p[1]) / (n - 1)
-        end
-    elseif pVecType == :ArcLen
-        for i in 2:(n - 1)
-            p[i] = p[1] + l[i - 1] / s * (p[end] - p[1])
-        end
-    end
-
-    lidx = 1
-    ridx = length(k)
-    while lidx <= (d + 1) && ridx >= (length(k) - d)
-        k[lidx] = p[1]
-        k[ridx] = p[end]
-        lidx += 1
-        ridx -= 1
-    end
-
-    ps = zeros(eltype(t), n - 2)
-    s = zero(eltype(t))
-    for i in 2:(n - 1)
-        s += p[i]
-        ps[i - 1] = s
+    # Clamped knot vector endpoints
+    for i in 1:(d + 1)
+        k[i] = t[1]
+        k[h + i] = t[end]
     end
-
     if knotVecType == :Uniform
-        # uniformly spaced knot vector
-        # this method is not recommended because, if it is used with the chord length method for global interpolation,
-        # the system of linear equations would be singular.
+        # Uniformly spaced interior knots
         for i in (d + 2):h
-            k[i] = k[1] + (i - d - 1) // (h - d) * (k[end] - k[1])
+            k[i] = t[1] + (i - d - 1) // (h - d) * (t[end] - t[1])
         end
     elseif knotVecType == :Average
-        # NOTE: verify that average method can be applied when size of k is less than size of p
-        # average spaced knot vector
-        idx = 1
-        if d + 2 <= h
-            k[d + 2] = 1 // d * ps[d]
-        end
-        for i in (d + 3):h
-            k[i] = 1 // d * (ps[idx + d] - ps[idx])
-            idx += 1
+        # Knot placement using Piegl-Tiller method for approximation
+        denom = max(d, 1)
+        delta = n / (h - denom + 1)
+        for j in 1:(h - d - 1)
+            frac = j * delta
+            i = min(floor(Int, frac), n - 1)
+            alpha = frac - i
+            k[d + 1 + j] = (1 - alpha) * t[i] + alpha * t[i + 1]
         end
     end
     # control points
@@ -1210,7 +1088,7 @@ function BSplineApprox(
     q = zeros(eltype(u), n)
     sc = zeros(eltype(t), n, h)
     for i in 1:n
-        spline_coefficients!(view(sc, i, :), d, k, p[i])
+        spline_coefficients!(view(sc, i, :), d, k, t[i])
     end
     for k in 2:(n - 1)
         q[k] = u[k] - sc[k, 1] * u[1] - sc[k, h] * u[end]
@@ -1229,13 +1107,13 @@ function BSplineApprox(
     c[2:(end - 1)] .= vec(P)
     sc = zeros(eltype(t), h)
     return BSplineApprox(
-        u, t, d, h, p, k, c, sc, pVecType, knotVecType,
+        u, t, d, h, k, c, sc, knotVecType,
         extrapolation_left, extrapolation_right, assume_linear_t
     )
 end
 
 function BSplineApprox(
-        u::AbstractArray{T, N}, t, d, h, pVecType, knotVecType;
+        u::AbstractArray{T, N}, t, d, h, knotVecType;
         extrapolation::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_left::ExtrapolationType.T = ExtrapolationType.None,
         extrapolation_right::ExtrapolationType.T = ExtrapolationType.None,
@@ -1248,62 +1126,27 @@ function BSplineApprox(
     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.")
-    s = zero(eltype(u))
-    p = zero(t)
-    k = zeros(eltype(t), h + d + 1)
-    l = zeros(eltype(u), n - 1)
-    p[1] = zero(eltype(t))
-    p[end] = one(eltype(t))
-
     ax_u = axes(u)[1:(end - 1)]
-
-    for i in 2:n
-        s += √((t[i] - t[i - 1])^2 + sum((u[ax_u..., i] - u[ax_u..., i - 1]) .^ 2))
-        l[i - 1] = s
-    end
-    if pVecType == :Uniform
-        for i in 2:(n - 1)
-            p[i] = p[1] + (i - 1) * (p[end] - p[1]) / (n - 1)
-        end
-    elseif pVecType == :ArcLen
-        for i in 2:(n - 1)
-            p[i] = p[1] + l[i - 1] / s * (p[end] - p[1])
-        end
-    end
-
-    lidx = 1
-    ridx = length(k)
-    while lidx <= (d + 1) && ridx >= (length(k) - d)
-        k[lidx] = p[1]
-        k[ridx] = p[end]
-        lidx += 1
-        ridx -= 1
-    end
-
-    ps = zeros(eltype(t), n - 2)
-    s = zero(eltype(t))
-    for i in 2:(n - 1)
-        s += p[i]
-        ps[i - 1] = s
+    k = zeros(eltype(t), h + d + 1)
+    # Clamped knot vector endpoints
+    for i in 1:(d + 1)
+        k[i] = t[1]
+        k[h + i] = t[end]
     end
-
     if knotVecType == :Uniform
-        # uniformly spaced knot vector
-        # this method is not recommended because, if it is used with the chord length method for global interpolation,
-        # the system of linear equations would be singular.
+        # Uniformly spaced interior knots
         for i in (d + 2):h
-            k[i] = k[1] + (i - d - 1) // (h - d) * (k[end] - k[1])
+            k[i] = t[1] + (i - d - 1) // (h - d) * (t[end] - t[1])
         end
     elseif knotVecType == :Average
-        # NOTE: verify that average method can be applied when size of k is less than size of p
-        # average spaced knot vector
-        idx = 1
-        if d + 2 <= h
-            k[d + 2] = 1 // d * ps[d]
-        end
-        for i in (d + 3):h
-            k[i] = 1 // d * (ps[idx + d] - ps[idx])
-            idx += 1
+        # Knot placement using Piegl-Tiller method for approximation
+        denom = max(d, 1)
+        delta = n / (h - denom + 1)
+        for j in 1:(h - d - 1)
+            frac = j * delta
+            i = min(floor(Int, frac), n - 1)
+            alpha = frac - i
+            k[d + 1 + j] = (1 - alpha) * t[i] + alpha * t[i + 1]
         end
     end
     # control points
@@ -1313,7 +1156,7 @@ function BSplineApprox(
     q = zeros(eltype(u), size(u)[1:(end - 1)]..., n)
     sc = zeros(eltype(t), n, h)
     for i in 1:n
-        spline_coefficients!(view(sc, i, :), d, k, p[i])
+        spline_coefficients!(view(sc, i, :), d, k, t[i])
     end
     for k in 2:(n - 1)
         q[
@@ -1338,7 +1181,7 @@ function BSplineApprox(
     c[ax_u..., 2:(end - 1)] = P
     sc = zeros(eltype(t), h)
     return BSplineApprox(
-        u, t, d, h, p, k, c, sc, pVecType, knotVecType,
+        u, t, d, h, k, c, sc, knotVecType,
         extrapolation_left, extrapolation_right, assume_linear_t
     )
 end
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index 47e18f26..478a0ec7 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -283,11 +283,6 @@ function _interpolate(
         t::Number,
         iguess
     )
-    t < A.t[1] && return A.u[1]
-    t > A.t[end] && return A.u[end]
-    # change t into param [0 1]
-    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)
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
     nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
@@ -304,11 +299,6 @@ function _interpolate(
         iguess
     )
     ax_u = axes(A.u)[1:(end - 1)]
-    t < A.t[1] && return A.u[ax_u..., 1]
-    t > A.t[end] && return A.u[ax_u..., end]
-    # change t into param [0 1]
-    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)
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), n) : A.sc
     nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
@@ -321,11 +311,6 @@ end
 
 # BSpline Curve Approx
 function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, iguess)
-    t < A.t[1] && return A.u[1]
-    t > A.t[end] && return A.u[end]
-    # change t into param [0 1]
-    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])
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
     nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
     ucum = zero(eltype(A.u))
@@ -339,11 +324,6 @@ function _interpolate(
         A::BSplineApprox{<:AbstractArray{<:Number}}, t::Number, iguess
     )
     ax_u = axes(A.u)[1:(end - 1)]
-    t < A.t[1] && return A.u[ax_u..., 1]
-    t > A.t[end] && return A.u[ax_u..., end]
-    # change t into param [0 1]
-    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])
     sc = t isa ForwardDiff.Dual ? zeros(eltype(t), A.h) : A.sc
     nonzero_coefficient_idxs = spline_coefficients!(sc, A.d, A.k, t)
     ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...)
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index c5ff26ef..6621e1c3 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -31,24 +31,45 @@ end
 
 function spline_coefficients!(N, d, k, u::Number)
     N .= zero(u)
+    n = length(N)
     if u == k[1]
         N[1] = one(u)
         return 1:1
     elseif u == k[end]
         N[end] = one(u)
-        return length(N):length(N)
+        return n:n
     else
-        i = findfirst(x -> x > u, k)::Int - 1
+        idx = findfirst(x -> x > u, k)
+        # For out-of-range points, extend the boundary polynomial span
+        i = if idx === nothing
+            # u > k[end]: use last span
+            findlast(j -> k[j] < k[end], 1:length(k))::Int
+        elseif idx == 1
+            # u < k[1]: use first span
+            findfirst(j -> k[j + 1] > k[1], 1:(length(k) - 1))::Int
+        else
+            idx - 1
+        end
         N[i] = one(u)
         for deg in 1:d
-            N[i - deg] = (k[i + 1] - u) / (k[i + 1] - k[i - deg + 1]) * N[i - deg + 1]
-            for j in (i - deg + 1):(i - 1)
-                N[j] = (u - k[j]) / (k[j + deg] - k[j]) * N[j] +
-                    (k[j + deg + 1] - u) / (k[j + deg + 1] - k[j + 1]) * N[j + 1]
+            ii = i - deg
+            if ii >= 1
+                denom = k[i + 1] - k[ii + 1]
+                N[ii] = denom != 0 ? (k[i + 1] - u) / denom * N[ii + 1] : zero(u)
+            end
+            for j in max(ii + 1, 1):(i - 1)
+                denom1 = k[j + deg] - k[j]
+                denom2 = k[j + deg + 1] - k[j + 1]
+                left = denom1 != 0 ? (u - k[j]) / denom1 * N[j] : zero(u)
+                right = denom2 != 0 ? (k[j + deg + 1] - u) / denom2 * N[j + 1] : zero(u)
+                N[j] = left + right
             end
-            N[i] = (u - k[i]) / (k[i + deg] - k[i]) * N[i]
+            denom = k[i + deg] - k[i]
+            N[i] = denom != 0 ? (u - k[i]) / denom * N[i] : zero(u)
         end
-        return (i - d):i
+        lo = max(i - d, 1)
+        hi = min(i, n)
+        return lo:hi
     end
 end
 
diff --git a/test/derivative_tests.jl b/test/derivative_tests.jl
index 6abab98f..9785e275 100644
--- a/test/derivative_tests.jl
+++ b/test/derivative_tests.jl
@@ -241,7 +241,6 @@ end
         args = [
             u, t, 2,
             :Uniform,
-            :Uniform,
         ],
         name = "BSpline Interpolation (Uniform, Uniform)"
     )
@@ -249,7 +248,6 @@ end
         BSplineInterpolation;
         args = [
             u, t, 2,
-            :ArcLen,
             :Average,
         ],
         name = "BSpline Interpolation (Arclen, Average)"
@@ -261,7 +259,6 @@ end
             3,
             4,
             :Uniform,
-            :Uniform,
         ],
         name = "BSpline Approx (Uniform, Uniform)"
     )
@@ -279,7 +276,6 @@ end
             u3d, t3d,
             2,
             :Uniform,
-            :Uniform,
         ],
         name = "BSpline Interpolation (Uniform, Uniform): AbstractArray"
     )
@@ -289,7 +285,6 @@ end
         args = [
             u3d, t3d,
             2,
-            :ArcLen,
             :Average,
         ],
         name = "BSpline Interpolation (Arclen, Average): AbstractArray"
@@ -302,7 +297,6 @@ end
             3,
             4,
             :Uniform,
-            :Uniform,
         ],
         name = "BSpline Approx (Uniform, Uniform): AbstractArray"
     )
@@ -313,7 +307,6 @@ end
             u3d, t3d,
             3,
             4,
-            :ArcLen,
             :Average,
         ],
         name = "BSpline Approx (Arclen, Average): AbstractArray"
diff --git a/test/integral_tests.jl b/test/integral_tests.jl
index 7f388df4..260cec01 100644
--- a/test/integral_tests.jl
+++ b/test/integral_tests.jl
@@ -6,7 +6,10 @@ using RegularizationTools
 using StableRNGs
 using Unitful
 
-function test_integral(method; args = [], kwargs = [], name::String)
+function test_integral(
+        method; args = [], kwargs = [], name::String,
+        test_cache_parameters::Bool = true
+    )
     func = method(args...; kwargs..., extrapolation = ExtrapolationType.Extension)
     (; t) = func
     t1 = minimum(t)
@@ -65,14 +68,16 @@ function test_integral(method; args = [], kwargs = [], name::String)
         func, t[1], t[end] + 1.0
     )
 
-    # Test integration with cached parameters
-    func = method(
-        args...; kwargs..., cache_parameters = true,
-        extrapolation = ExtrapolationType.Extension
-    )
-    qint, err = quadgk(func, t1 - 1, t1; atol = 1.0e-12, rtol = 1.0e-12)
-    aint = integral(func, t1 - 1, t1)
-    return @test isapprox(qint, aint, atol = 1.0e-6, rtol = 1.0e-8)
+    return if test_cache_parameters
+        # Test integration with cached parameters
+        func = method(
+            args...; kwargs..., cache_parameters = true,
+            extrapolation = ExtrapolationType.Extension
+        )
+        qint, err = quadgk(func, t1 - 1, t1; atol = 1.0e-12, rtol = 1.0e-12)
+        aint = integral(func, t1 - 1, t1)
+        @test isapprox(qint, aint, atol = 1.0e-6, rtol = 1.0e-8)
+    end
 end
 
 @testset "LinearInterpolation" begin
@@ -253,19 +258,37 @@ end
 end
 
 @testset "BSplineInterpolation" begin
-    t = [0, 62.25, 109.66, 162.66, 205.8, 252.3]
+    t = Float64[0, 62.25, 109.66, 162.66, 205.8, 252.3]
     u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
-    A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform)
-    @test_throws DataInterpolations.IntegralNotFoundError integral(A, 1.0, 100.0)
-    @test_throws DataInterpolations.IntegralNotFoundError integral(A, 50.0)
+    test_integral(
+        BSplineInterpolation;
+        args = [u, t, 2, :Uniform],
+        name = "BSpline Interpolation (d=2, Uniform)",
+        test_cache_parameters = false
+    )
+    test_integral(
+        BSplineInterpolation;
+        args = [u, t, 3, :Average],
+        name = "BSpline Interpolation (d=3, Average)",
+        test_cache_parameters = false
+    )
 end
 
 @testset "BSplineApprox" begin
-    t = [0, 62.25, 109.66, 162.66, 205.8, 252.3]
+    t = Float64[0, 62.25, 109.66, 162.66, 205.8, 252.3]
     u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
-    A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform)
-    @test_throws DataInterpolations.IntegralNotFoundError integral(A, 1.0, 100.0)
-    @test_throws DataInterpolations.IntegralNotFoundError integral(A, 50.0)
+    test_integral(
+        BSplineApprox;
+        args = [u, t, 2, 4, :Uniform],
+        name = "BSpline Approx (d=2, Uniform)",
+        test_cache_parameters = false
+    )
+    test_integral(
+        BSplineApprox;
+        args = [u, t, 3, 4, :Average],
+        name = "BSpline Approx (d=3, Average)",
+        test_cache_parameters = false
+    )
 end
 
 # issue #385
diff --git a/test/interface.jl b/test/interface.jl
index 5be8527d..125f51f3 100644
--- a/test/interface.jl
+++ b/test/interface.jl
@@ -40,12 +40,12 @@ end
         @inferred method(u, t)
     end
     @testset "BSplineInterpolation" begin
-        @inferred BSplineInterpolation(u, t, 3, :Uniform, :Uniform)
-        @inferred BSplineInterpolation(u, t, 3, :ArcLen, :Average)
+        @inferred BSplineInterpolation(u, t, 3, :Uniform)
+        @inferred BSplineInterpolation(u, t, 3, :Average)
     end
     @testset "BSplineApprox" begin
-        @inferred BSplineApprox(u, t, 3, 5, :Uniform, :Uniform)
-        @inferred BSplineApprox(u, t, 3, 5, :ArcLen, :Average)
+        @inferred BSplineApprox(u, t, 3, 5, :Uniform)
+        @inferred BSplineApprox(u, t, 3, 5, :Average)
     end
     du = ones(10)
     ddu = zeros(10)
diff --git a/test/interpolation_tests.jl b/test/interpolation_tests.jl
index 8a1c7f6f..47c2b871 100644
--- a/test/interpolation_tests.jl
+++ b/test/interpolation_tests.jl
@@ -935,10 +935,10 @@ end
         t = [0, 62.25, 109.66, 162.66, 205.8, 252.3]
         u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
         test_interpolation_type(BSplineInterpolation)
-        A = @inferred(BSplineInterpolation(u, t, 2, :Uniform, :Uniform))
+        A = @inferred(BSplineInterpolation(u, t, 2, :Uniform))
 
-        @test [A(25.0), A(80.0)] == [13.454197730061425, 10.305633616059845]
-        @test [A(190.0), A(225.0)] == [14.07428439395079, 11.057784141519251]
+        @test [A(25.0), A(80.0)] == [14.411908462307684, 9.99346697254525]
+        @test [A(190.0), A(225.0)] == [13.56561617594697, 11.297503333875742]
         @test [A(t[1]), A(t[end])] == [u[1], u[end]]
         test_cached_index(A)
         @test @inferred(output_dim(A)) == 0
@@ -947,45 +947,45 @@ end
         # Test extrapolation
         A = @inferred(
             BSplineInterpolation(
-                u, t, 2, :Uniform, :Uniform; extrapolation = ExtrapolationType.Extension
+                u, t, 2, :Uniform; extrapolation = ExtrapolationType.Constant
             )
         )
         @test A(-1.0) == u[1]
         @test A(300.0) == u[end]
-        A = @inferred(BSplineInterpolation(u, t, 2, :Uniform, :Uniform))
+        A = @inferred(BSplineInterpolation(u, t, 2, :Uniform))
         @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
         @test_throws DataInterpolations.RightExtrapolationError A(300.0)
 
-        A = @inferred(BSplineInterpolation(u, t, 2, :ArcLen, :Average))
+        A = @inferred(BSplineInterpolation(u, t, 2, :Average))
 
-        @test [A(25.0), A(80.0)] ≈ [13.363814458968484, 10.685201117692609]
-        @test [A(190.0), A(225.0)] ≈ [13.437481084762863, 11.367034741256461]
+        @test [A(25.0), A(80.0)] ≈ [13.364285794535945, 10.683641750973738]
+        @test [A(190.0), A(225.0)] ≈ [13.438136405352909, 11.365386175733823]
         @test [A(t[1]), A(t[end])] ≈ [u[1], u[end]]
 
         @test_throws ErrorException("BSplineInterpolation needs at least d + 1, i.e. 4 points.") BSplineInterpolation(
-            u[1:3], t[1:3], 3, :Uniform, :Uniform
+            u[1:3], t[1:3], 3, :Uniform
         )
         @test_throws ErrorException("BSplineInterpolation needs at least d + 1, i.e. 5 points.") BSplineInterpolation(
-            u[1:4], t[1:4], 4, :ArcLen, :Average
+            u[1:4], t[1:4], 4, :Average
         )
-        @test_nowarn BSplineInterpolation(u[1:3], t[1:3], 2, :Uniform, :Uniform)
+        @test_nowarn BSplineInterpolation(u[1:3], t[1:3], 2, :Uniform)
 
         # Test extrapolation
         A = @inferred(
             BSplineInterpolation(
-                u, t, 2, :ArcLen, :Average; extrapolation = ExtrapolationType.Extension
+                u, t, 2, :Average; extrapolation = ExtrapolationType.Constant
             )
         )
         @test A(-1.0) == u[1]
         @test A(300.0) == u[end]
-        A = @inferred(BSplineInterpolation(u, t, 2, :ArcLen, :Average))
+        A = @inferred(BSplineInterpolation(u, t, 2, :Average))
         @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
         @test_throws DataInterpolations.RightExtrapolationError A(300.0)
 
         @testset "AbstractMatrix" begin
             t = 0.1:0.1:1.0
             u2d = [sin.(t) cos.(t)]' |> collect
-            A = @inferred(BSplineInterpolation(u2d, t, 2, :Uniform, :Uniform))
+            A = @inferred(BSplineInterpolation(u2d, t, 2, :Uniform))
             t_test = 0.1:0.05:1.0
             u_test = reduce(hcat, A.(t_test))
             @test isapprox(u_test[1, :], sin.(t_test), atol = 1.0e-3)
@@ -993,7 +993,7 @@ end
             @test @inferred(output_dim(A)) == 1
             @test @inferred(output_size(A)) == (2,)
 
-            A = @inferred(BSplineInterpolation(u2d, t, 2, :ArcLen, :Average))
+            A = @inferred(BSplineInterpolation(u2d, t, 2, :Average))
             u_test = reduce(hcat, A.(t_test))
             @test isapprox(u_test[1, :], sin.(t_test), atol = 1.0e-3)
             @test isapprox(u_test[2, :], cos.(t_test), atol = 1.0e-3)
@@ -1007,7 +1007,7 @@ end
             ]
             t = 0.1:0.1:1.0
             u3d = cat(f3d.(t)..., dims = 3)
-            A = @inferred(BSplineInterpolation(u3d, t, 2, :Uniform, :Uniform))
+            A = @inferred(BSplineInterpolation(u3d, t, 2, :Uniform))
             t_test = 0.1:0.05:1.0
             u_test = reduce(hcat, A.(t_test))
             f_test = reduce(hcat, f3d.(t_test))
@@ -1015,7 +1015,7 @@ end
             @test @inferred(output_dim(A)) == 2
             @test @inferred(output_size(A)) == (2, 2)
 
-            A = @inferred(BSplineInterpolation(u3d, t, 2, :ArcLen, :Average))
+            A = @inferred(BSplineInterpolation(u3d, t, 2, :Average))
             t_test = 0.1:0.05:1.0
             u_test = reduce(hcat, A.(t_test))
             @test isapprox(u_test, f_test, atol = 1.0e-2)
@@ -1028,41 +1028,41 @@ end
         test_interpolation_type(BSplineApprox)
         t = [0, 62.25, 109.66, 162.66, 205.8, 252.3]
         u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
-        A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform)
+        A = BSplineApprox(u, t, 2, 4, :Uniform)
 
-        @test [A(25.0), A(80.0)] ≈ [12.979802931218234, 10.914310609953178]
-        @test [A(190.0), A(225.0)] ≈ [13.851245975109263, 12.963685868886575]
+        @test [A(25.0), A(80.0)] ≈ [12.653438633006644, 10.829963801404205]
+        @test [A(190.0), A(225.0)] ≈ [13.688412211160633, 12.785204994978452]
         @test [A(t[1]), A(t[end])] ≈ [u[1], u[end]]
         test_cached_index(A)
 
         @test_throws ErrorException("BSplineApprox needs at least d + 1, i.e. 3 control points.") BSplineApprox(
-            u, t, 2, 2, :Uniform, :Uniform
+            u, t, 2, 2, :Uniform
         )
         @test_throws ErrorException("BSplineApprox needs at least d + 1, i.e. 4 control points.") BSplineApprox(
-            u, t, 3, 3, :ArcLen, :Average
+            u, t, 3, 3, :Average
         )
-        @test_nowarn BSplineApprox(u, t, 2, 3, :Uniform, :Uniform)
+        @test_nowarn BSplineApprox(u, t, 2, 3, :Uniform)
 
         # Test extrapolation
         A = BSplineApprox(
-            u, t, 2, 4, :Uniform, :Uniform; extrapolation = ExtrapolationType.Extension
+            u, t, 2, 4, :Uniform; extrapolation = ExtrapolationType.Constant
         )
         @test A(-1.0) == u[1]
         @test A(300.0) == u[end]
-        A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform)
+        A = BSplineApprox(u, t, 2, 4, :Uniform)
         @test_throws DataInterpolations.LeftExtrapolationError A(-1.0)
         @test_throws DataInterpolations.RightExtrapolationError A(300.0)
 
         @testset "AbstractMatrix" begin
             t = 0.1:0.1:1.0
             u2d = [sin.(t) cos.(t)]' |> collect
-            A = BSplineApprox(u2d, t, 2, 5, :Uniform, :Uniform)
+            A = BSplineApprox(u2d, t, 2, 5, :Uniform)
             t_test = 0.1:0.05:1.0
             u_test = reduce(hcat, A.(t_test))
             @test isapprox(u_test[1, :], sin.(t_test), atol = 1.0e-3)
             @test isapprox(u_test[2, :], cos.(t_test), atol = 1.0e-3)
 
-            A = BSplineApprox(u2d, t, 2, 5, :ArcLen, :Average)
+            A = BSplineApprox(u2d, t, 2, 5, :Average)
             u_test = reduce(hcat, A.(t_test))
             @test isapprox(u_test[1, :], sin.(t_test), atol = 1.0e-2)
             @test isapprox(u_test[2, :], cos.(t_test), atol = 1.0e-2)
@@ -1074,13 +1074,13 @@ end
             ]
             t = 0.1:0.1:1.0
             u3d = cat(f3d.(t)..., dims = 3)
-            A = BSplineApprox(u3d, t, 2, 6, :Uniform, :Uniform)
+            A = BSplineApprox(u3d, t, 2, 6, :Uniform)
             t_test = 0.1:0.05:1.0
             u_test = reduce(hcat, A.(t_test))
             f_test = reduce(hcat, f3d.(t_test))
             @test isapprox(u_test, f_test, atol = 1.0e-2)
 
-            A = BSplineApprox(u3d, t, 2, 7, :ArcLen, :Average)
+            A = BSplineApprox(u3d, t, 2, 7, :Average)
             t_test = 0.1:0.05:1.0
             u_test = reduce(hcat, A.(t_test))
             @test isapprox(u_test, f_test, atol = 1.0e-2)
@@ -1385,6 +1385,6 @@ end
     @test_nowarn cs = CubicSpline(x, t)
 
     @test_throws Exception ai = AkimaInterpolation(x, t)
-    @test_throws Exception bsi = BSplineInterpolation(x, t, 3, :ArcLen, :Average)
+    @test_throws Exception bsi = BSplineInterpolation(x, t, 3, :Average)
     @test_throws Exception pc = PCHIPInterpolation(x, t)
 end
diff --git a/test/mooncake_tests.jl b/test/mooncake_tests.jl
index 13db757a..05182119 100644
--- a/test/mooncake_tests.jl
+++ b/test/mooncake_tests.jl
@@ -130,9 +130,9 @@ end
     t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
     test_mooncake_ugrad(
         BSplineInterpolation, u, t;
-        args_after = [2, :Uniform, :Uniform], name = "BSpline Interpolation"
+        args_after = [2, :Uniform], name = "BSpline Interpolation"
     )
-    A = BSplineInterpolation(u, t, 2, :Uniform, :Uniform; extrapolation = ExtrapolationType.Extension)
+    A = BSplineInterpolation(u, t, 2, :Uniform; extrapolation = ExtrapolationType.Extension)
     test_mooncake_tgrad(A, t; name = "BSpline Interpolation")
 end
 
@@ -141,8 +141,8 @@ end
     t = [0.0, 62.25, 109.66, 162.66, 205.8, 252.3]
     test_mooncake_ugrad(
         BSplineApprox, u, t;
-        args_after = [2, 4, :Uniform, :Uniform], name = "BSpline Approximation"
+        args_after = [2, 4, :Uniform], name = "BSpline Approximation"
     )
-    A = BSplineApprox(u, t, 2, 4, :Uniform, :Uniform; extrapolation = ExtrapolationType.Extension)
+    A = BSplineApprox(u, t, 2, 4, :Uniform; extrapolation = ExtrapolationType.Extension)
     test_mooncake_tgrad(A, t; name = "BSpline Approximation")
 end
diff --git a/test/show.jl b/test/show.jl
index aa188d48..feb43d9e 100644
--- a/test/show.jl
+++ b/test/show.jl
@@ -46,14 +46,14 @@ end
         )
     end
     @testset "BSplineInterpolation" begin
-        A = BSplineInterpolation(x, t, 3, :Uniform, :Uniform)
+        A = BSplineInterpolation(x, t, 3, :Uniform)
         @test startswith(
             sprint(io -> show(io, MIME"text/plain"(), A)),
             "BSplineInterpolation with 5 points, with degree 3\n"
         )
     end
     @testset "BSplineApprox" begin
-        A = BSplineApprox(x, t, 2, 4, :Uniform, :Uniform)
+        A = BSplineApprox(x, t, 2, 4, :Uniform)
         @test startswith(
             sprint(io -> show(io, MIME"text/plain"(), A)),
             "BSplineApprox with 5 points, with degree 2, number of control points 4\n"
diff --git a/test/sparseconnectivitytracer_tests.jl b/test/sparseconnectivitytracer_tests.jl
index fc6a7307..bc2245d0 100644
--- a/test/sparseconnectivitytracer_tests.jl
+++ b/test/sparseconnectivitytracer_tests.jl
@@ -241,8 +241,8 @@ end
             InterpolationTest(AkimaInterpolation(u, t)),
             InterpolationTest(QuadraticSpline(u, t)),
             InterpolationTest(CubicSpline(u, t)),
-            InterpolationTest(BSplineInterpolation(u, t, 3, :ArcLen, :Average)),
-            InterpolationTest(BSplineApprox(u, t, 3, 4, :ArcLen, :Average)),
+            InterpolationTest(BSplineInterpolation(u, t, 3, :Average)),
+            InterpolationTest(BSplineApprox(u, t, 3, 4, :Average)),
             InterpolationTest(PCHIPInterpolation(u, t)),
             InterpolationTest(CubicHermiteSpline(du, u, t)),
             InterpolationTest(QuinticHermiteSpline(ddu, du, u, t)),
@@ -267,10 +267,10 @@ for N in (2, 5)
                 InterpolationTest(LagrangeInterpolation(um, t)),
                 ## The following interpolations appear to not be supported on N dimensions as of DataInterpolations v6.2.0:
                 # InterpolationTest(AkimaInterpolation(um, t)),
-                # InterpolationTest(BSplineApprox(um, t, 3, 4, :ArcLen, :Average)),
+                # InterpolationTest(BSplineApprox(um, t, 3, 4, :Average)),
                 # InterpolationTest(QuadraticSpline(um, t)),
                 # InterpolationTest(CubicSpline(um, t)),
-                # InterpolationTest(BSplineInterpolation(um, t, 3, :ArcLen, :Average)),
+                # InterpolationTest(BSplineInterpolation(um, t, 3, :Average)),
                 # InterpolationTest(PCHIPInterpolation(um, t)),
             )
             test_jacobian(t)
diff --git a/test/zygote_tests.jl b/test/zygote_tests.jl
index 7fcf98f8..5e9ea4c8 100644
--- a/test/zygote_tests.jl
+++ b/test/zygote_tests.jl
@@ -126,11 +126,11 @@ end
     t = [0, 62.25, 109.66, 162.66, 205.8, 252.3]
     u = [14.7, 11.51, 10.41, 14.95, 12.24, 11.22]
     test_zygote(
-        BSplineInterpolation, u, t; args_after = [2, :Uniform, :Uniform],
+        BSplineInterpolation, u, t; args_after = [2, :Uniform],
         name = "BSpline Interpolation"
     )
     test_zygote(
-        BSplineApprox, u, t; args_after = [2, 4, :Uniform, :Uniform],
+        BSplineApprox, u, t; args_after = [2, 4, :Uniform],
         name = "BSpline approximation"
     )
 end