diff --git a/Project.toml b/Project.toml
index 894332f7..cb763806 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -73,11 +74,10 @@ SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "AllocCheck", "BenchmarkTools", "SafeTestsets", "SciMLTesting", "ChainRulesCore", "Mooncake", "Optim", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "StaticArrays", "Symbolics", "Unitful", "Zygote", "SparseConnectivityTracer"]
+test = ["Aqua", "AllocCheck", "BenchmarkTools", "SafeTestsets", "SciMLTesting", "ChainRulesCore", "Mooncake", "Optim", "Test", "StableRNGs", "FiniteDifferences", "QuadGK", "ForwardDiff", "Symbolics", "Unitful", "Zygote", "SparseConnectivityTracer"]
diff --git a/src/DataInterpolations.jl b/src/DataInterpolations.jl
index ffaf8157..73f0c34b 100644
--- a/src/DataInterpolations.jl
+++ b/src/DataInterpolations.jl
@@ -9,6 +9,7 @@ using RecipesBase: RecipesBase, @recipe, @series
 using PrettyTables: PrettyTables, pretty_table
 using ForwardDiff: ForwardDiff
 using EnumX: EnumX, @enumx
+using StaticArrays: MVector, SVector
 import FindFirstFunctions
 import FindFirstFunctions: Guesser
 
diff --git a/src/derivatives.jl b/src/derivatives.jl
index 268e3ec0..68bb015e 100644
--- a/src/derivatives.jl
+++ b/src/derivatives.jl
@@ -256,14 +256,16 @@ function _derivative(A::BSplineInterpolation{<:AbstractVector{<:Number}}, t::Num
     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_)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
         ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2])
     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])
+        # Stack-allocated basis window: differentiation must be reentrant (#532)
+        vals, offset, m = bspline_nonzero_coefficients(A.d - 1, A.k, t_, n)
+        @inbounds for l in 1:m
+            i = offset + l - 1
+            (1 <= i <= n - 1) || continue
+            ducum += vals[l] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
         end
     end
     return ducum * A.d * scale
@@ -280,15 +282,17 @@ function _derivative(
     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_)
     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])
     else
-        for i in 1:(n - 1)
+        # Stack-allocated basis window: differentiation must be reentrant (#532)
+        vals, offset, m = bspline_nonzero_coefficients(A.d - 1, A.k, t_, n)
+        @inbounds for l in 1:m
+            i = offset + l - 1
+            (1 <= i <= n - 1) || continue
             ducum = ducum +
-                sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
+                vals[l] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
                 (A.k[i + A.d + 1] - A.k[i + 1])
         end
     end
@@ -302,14 +306,16 @@ function _derivative(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, ig
     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_)
     ducum = zero(eltype(A.u))
     if t == A.t[1]
         ducum = (A.c[2] - A.c[1]) / (A.k[A.d + 2])
     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])
+        # Stack-allocated basis window: differentiation must be reentrant (#532)
+        vals, offset, m = bspline_nonzero_coefficients(A.d - 1, A.k, t_, A.h)
+        @inbounds for l in 1:m
+            i = offset + l - 1
+            (1 <= i <= A.h - 1) || continue
+            ducum += vals[l] * (A.c[i + 1] - A.c[i]) / (A.k[i + A.d + 1] - A.k[i + 1])
         end
     end
     return ducum * A.d * scale
@@ -325,14 +331,16 @@ function _derivative(
     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_)
     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])
     else
-        for i in 1:(A.h - 1)
-            ducum += sc[i + 1] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
+        # Stack-allocated basis window: differentiation must be reentrant (#532)
+        vals, offset, m = bspline_nonzero_coefficients(A.d - 1, A.k, t_, A.h)
+        @inbounds for l in 1:m
+            i = offset + l - 1
+            (1 <= i <= A.h - 1) || continue
+            ducum = ducum + vals[l] * (A.c[ax_u..., i + 1] - A.c[ax_u..., i]) /
                 (A.k[i + A.d + 1] - A.k[i + 1])
         end
     end
diff --git a/src/interpolation_caches.jl b/src/interpolation_caches.jl
index cdcd6995..ea59cb85 100644
--- a/src/interpolation_caches.jl
+++ b/src/interpolation_caches.jl
@@ -668,7 +668,7 @@ function QuadraticSpline(
     k, A = quadratic_spline_params(t, sc)
     c = A \ u
 
-    p = QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters)
+    p = QuadraticSplineParameterCache(u, t, k, c, cache_parameters)
     A = QuadraticSpline(
         u, t, nothing, p, k, c, sc, extrapolation_left,
         extrapolation_right, cache_parameters, t_props
@@ -710,7 +710,7 @@ function QuadraticSpline(
         end
     end
 
-    p = QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters)
+    p = QuadraticSplineParameterCache(u, t, k, c, cache_parameters)
     A = QuadraticSpline(
         u, t, nothing, p, k, c, sc, extrapolation_left,
         extrapolation_right, cache_parameters, t_props
@@ -1020,6 +1020,8 @@ function BSplineInterpolation(
     t_props = something(search_properties, FindFirstFunctions.SearchProperties(t))
     n = length(t)
     n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d + 1) points.")
+    d ≥ BSPLINE_STACK_MAXLEN &&
+        error("BSplineInterpolation supports degree d < $(BSPLINE_STACK_MAXLEN); got d = $d.")
     s = zero(eltype(u))
     p = zero(t)
     k = zeros(eltype(t), n + d + 1)
@@ -1101,6 +1103,8 @@ function BSplineInterpolation(
     t_props = something(search_properties, FindFirstFunctions.SearchProperties(t))
     n = length(t)
     n < d + 1 && error("BSplineInterpolation needs at least d + 1, i.e. $(d + 1) points.")
+    d ≥ BSPLINE_STACK_MAXLEN &&
+        error("BSplineInterpolation supports degree d < $(BSPLINE_STACK_MAXLEN); got d = $d.")
     s = zero(eltype(u))
     p = zero(t)
     k = zeros(eltype(t), n + d + 1)
@@ -1274,6 +1278,8 @@ function BSplineApprox(
     t_props = something(search_properties, FindFirstFunctions.SearchProperties(t))
     n = length(t)
     h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d + 1) control points.")
+    d ≥ BSPLINE_STACK_MAXLEN &&
+        error("BSplineApprox supports degree d < $(BSPLINE_STACK_MAXLEN); got d = $d.")
     s = zero(eltype(u))
     p = zero(t)
     k = zeros(eltype(t), h + d + 1)
@@ -1376,6 +1382,8 @@ function BSplineApprox(
     t_props = something(search_properties, FindFirstFunctions.SearchProperties(t))
     n = length(t)
     h < d + 1 && error("BSplineApprox needs at least d + 1, i.e. $(d + 1) control points.")
+    d ≥ BSPLINE_STACK_MAXLEN &&
+        error("BSplineApprox supports degree d < $(BSPLINE_STACK_MAXLEN); got d = $d.")
     s = zero(eltype(u))
     p = zero(t)
     k = zeros(eltype(t), h + d + 1)
diff --git a/src/interpolation_methods.jl b/src/interpolation_methods.jl
index aecdeec9..1929a0cf 100644
--- a/src/interpolation_methods.jl
+++ b/src/interpolation_methods.jl
@@ -1264,11 +1264,11 @@ function _interpolate(
     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)
+    # Stack-allocated basis window: evaluation must be reentrant for thread safety (#532)
+    vals, offset, m = bspline_nonzero_coefficients(A.d, A.k, t, n)
     ucum = zero(eltype(A.u))
-    for i in nonzero_coefficient_idxs
-        ucum += sc[i] * A.c[i]
+    @inbounds for l in 1:m
+        ucum += vals[l] * A.c[offset + l]
     end
     return ucum
 end
@@ -1285,11 +1285,11 @@ function _interpolate(
     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)
+    # Stack-allocated basis window: evaluation must be reentrant for thread safety (#532)
+    vals, offset, m = bspline_nonzero_coefficients(A.d, A.k, t, n)
     ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...)
-    for i in nonzero_coefficient_idxs
-        ucum = ucum + (sc[i] * A.c[ax_u..., i])
+    @inbounds for l in 1:m
+        ucum = ucum + (vals[l] * A.c[ax_u..., offset + l])
     end
     return ucum
 end
@@ -1301,11 +1301,11 @@ function _interpolate(A::BSplineApprox{<:AbstractVector{<:Number}}, t::Number, i
     # 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)
+    # Stack-allocated basis window: evaluation must be reentrant for thread safety (#532)
+    vals, offset, m = bspline_nonzero_coefficients(A.d, A.k, t, A.h)
     ucum = zero(eltype(A.u))
-    for i in nonzero_coefficient_idxs
-        ucum += sc[i] * A.c[i]
+    @inbounds for l in 1:m
+        ucum += vals[l] * A.c[offset + l]
     end
     return ucum
 end
@@ -1319,11 +1319,11 @@ function _interpolate(
     # 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)
+    # Stack-allocated basis window: evaluation must be reentrant for thread safety (#532)
+    vals, offset, m = bspline_nonzero_coefficients(A.d, A.k, t, A.h)
     ucum = zeros(eltype(A.u), size(A.u)[1:(end - 1)]...)
-    for i in nonzero_coefficient_idxs
-        ucum = ucum + (sc[i] * A.c[ax_u..., i])
+    @inbounds for l in 1:m
+        ucum = ucum + (vals[l] * A.c[ax_u..., offset + l])
     end
     return ucum
 end
diff --git a/src/interpolation_utils.jl b/src/interpolation_utils.jl
index ceea6aa6..72bf7c6f 100644
--- a/src/interpolation_utils.jl
+++ b/src/interpolation_utils.jl
@@ -67,6 +67,46 @@ function spline_coefficients!(N, d, k, u::AbstractVector)
     return nothing
 end
 
+# A degree-`d` B-spline has only `d + 1` nonzero basis functions at any point, so
+# evaluation needs a scratch buffer of just `d + 1` entries — not one sized to the
+# full knot/control-point vector. `bspline_nonzero_coefficients` computes those
+# values into a stack-allocated `MVector` (no heap allocation, hence reentrant /
+# thread-safe, #532) for degrees up to `BSPLINE_STACK_MAXLEN - 1`; callers fall back
+# to the heap `spline_coefficients!` for larger degrees, which essentially never occur.
+const BSPLINE_STACK_MAXLEN = 16
+
+# Returns `(vals, offset, m)`: `vals[l]` is the basis weight of control point
+# `offset + l` for `l in 1:m`. `ncp` is the number of control points (needed only to
+# place the single nonzero weight at the right boundary). Mirrors the locator and
+# recurrence of `spline_coefficients!`, but writes into a local window indexed `1:m`
+# instead of absolute knot indices.
+@inline function bspline_nonzero_coefficients(d::Integer, k, u::T, ncp::Integer) where {T}
+    N = zero(MVector{BSPLINE_STACK_MAXLEN, T})
+    if u == k[1]
+        @inbounds N[1] = one(u)
+        return SVector(N), 0, 1
+    elseif u == k[end]
+        @inbounds N[1] = one(u)
+        return SVector(N), ncp - 1, 1
+    end
+    i = searchsortedlast(k, u)
+    # Local index of global knot index `g` is `g + off` (so `i - d → 1`, `i → d + 1`).
+    off = d + 1 - i
+    @inbounds begin
+        N[i + off] = one(u)
+        for deg in 1:d
+            N[i - deg + off] = (k[i + 1] - u) / (k[i + 1] - k[i - deg + 1]) *
+                N[i - deg + 1 + off]
+            for j in (i - deg + 1):(i - 1)
+                N[j + off] = (u - k[j]) / (k[j + deg] - k[j]) * N[j + off] +
+                    (k[j + deg + 1] - u) / (k[j + deg + 1] - k[j + 1]) * N[j + 1 + off]
+            end
+            N[i + off] = (u - k[i]) / (k[i + deg] - k[i]) * N[i + off]
+        end
+    end
+    return SVector(N), i - d - 1, d + 1
+end
+
 function quadratic_spline_params(t::AbstractVector, sc::AbstractVector)
     # Duplicate time points make the collocation system singular
     if any(i -> t[i] == t[i + 1], 1:(length(t) - 1))
@@ -301,7 +341,7 @@ function get_parameters(A::QuadraticSpline, idx)
     return if A.cache_parameters
         A.p.α[idx], A.p.β[idx]
     else
-        quadratic_spline_parameters(A.u, A.t, A.k, A.c, A.sc, idx)
+        quadratic_spline_parameters(A.u, A.t, A.k, A.c, idx)
     end
 end
 
diff --git a/src/parameter_caches.jl b/src/parameter_caches.jl
index d30930a4..cccc4e7d 100644
--- a/src/parameter_caches.jl
+++ b/src/parameter_caches.jl
@@ -131,21 +131,21 @@ struct QuadraticSplineParameterCache{pType}
     β::pType
 end
 
-function QuadraticSplineParameterCache(u, t, k, c, sc, cache_parameters)
+function QuadraticSplineParameterCache(u, t, k, c, cache_parameters)
     return if cache_parameters
         parameters = quadratic_spline_parameters.(
-            Ref(u), Ref(t), Ref(k), Ref(c), Ref(sc), 1:(length(t) - 1)
+            Ref(u), Ref(t), Ref(k), Ref(c), 1:(length(t) - 1)
         )
         α, β = collect.(eachrow(stack(collect.(parameters))))
         QuadraticSplineParameterCache(α, β)
     else
         # Compute parameters once to infer types
-        α, β = quadratic_spline_parameters(u, t, k, c, sc, 1)
+        α, β = quadratic_spline_parameters(u, t, k, c, 1)
         QuadraticSplineParameterCache(typeof(α)[], typeof(β)[])
     end
 end
 
-function quadratic_spline_parameters(u, t, k, c, sc, idx)
+function quadratic_spline_parameters(u, t, k, c, idx)
     uᵢ₊ = if length(t) == 2
         # For 2 data points the knot vector has boundary multiplicity 2 < degree + 1,
         # so the B-spline basis cannot be evaluated at interior points; the spline
@@ -153,11 +153,24 @@ function quadratic_spline_parameters(u, t, k, c, sc, idx)
         (u[1] + u[2]) / 2
     else
         tᵢ₊ = (t[idx] + t[idx + 1]) / 2
-        nonzero_coefficient_idxs = spline_coefficients!(sc, 2, k, tᵢ₊)
+        # Value of the spline at the segment midpoint via the (degree 2) B-spline basis.
+        # The three nonzero basis values are evaluated into local variables rather than a
+        # shared scratch buffer so that evaluation is reentrant / thread-safe (#532).
+        # `tᵢ₊` is always interior, so only the non-boundary branch of the Cox-de Boor
+        # recursion is needed (cf. `spline_coefficients!`).
+        i = findfirst(x -> x > tᵢ₊, k)::Int - 1
+        w₁ = (k[i + 1] - tᵢ₊) / (k[i + 1] - k[i])
+        w₂ = (tᵢ₊ - k[i]) / (k[i + 1] - k[i])
+        N₁ = (k[i + 1] - tᵢ₊) / (k[i + 1] - k[i - 1]) * w₁
+        N₂ = (tᵢ₊ - k[i - 1]) / (k[i + 1] - k[i - 1]) * w₁ +
+            (k[i + 2] - tᵢ₊) / (k[i + 2] - k[i]) * w₂
+        N₃ = (tᵢ₊ - k[i]) / (k[i + 2] - k[i]) * w₂
+        # Seed the accumulator with `zero(first(u))` (as the buffer-based version did) so
+        # `uᵢ₊` keeps the element type of `u` for vector-valued data.
         uᵢ₊ = zero(first(u))
-        for j in nonzero_coefficient_idxs
-            uᵢ₊ += sc[j] * c[j]
-        end
+        uᵢ₊ += N₁ * c[i - 2]
+        uᵢ₊ += N₂ * c[i - 1]
+        uᵢ₊ += N₃ * c[i]
         uᵢ₊
     end
     α = 2 * (u[idx + 1] + u[idx]) - 4uᵢ₊
diff --git a/test/Methods/thread_safety_tests.jl b/test/Methods/thread_safety_tests.jl
new file mode 100644
index 00000000..f7a327ec
--- /dev/null
+++ b/test/Methods/thread_safety_tests.jl
@@ -0,0 +1,70 @@
+using DataInterpolations
+using DataInterpolations: derivative
+using Test
+using Base.Threads
+
+# Regression tests for issue #532: a `BSplineInterpolation`, `BSplineApprox` or
+# (default, uncached) `QuadraticSpline` kept a scratch buffer in its `sc` field
+# and overwrote it in place on every evaluation. Because a "read" (evaluation)
+# mutated shared state, evaluating one shared interpolant from several threads
+# silently returned wrong values. Evaluation (and differentiation) must be
+# reentrant: it must not write to any state stored in the interpolant.
+
+function thread_safety_interpolants()
+    u = [0.0, 1.0, 0.0, 1.0, 0.0, 1.0]
+    t = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
+    umat = [
+        0.0 1.0 0.0 1.0 0.0 1.0
+        1.0 0.0 1.0 0.0 1.0 0.0
+    ]            # 2 × length(t): last axis indexes t
+    return [
+        ("BSplineInterpolation", BSplineInterpolation(u, t, 3, :Uniform, :Average)),
+        # degree 1 == "linear" B-spline: same `_interpolate` method, so equally affected
+        ("BSplineInterpolation (degree 1)", BSplineInterpolation(u, t, 1, :Uniform, :Average)),
+        ("BSplineApprox", BSplineApprox(u, t, 3, 4, :Uniform, :Average)),
+        ("BSplineInterpolation (array)", BSplineInterpolation(umat, t, 3, :Uniform, :Average)),
+        ("BSplineApprox (array)", BSplineApprox(umat, t, 3, 4, :Uniform, :Average)),
+        ("QuadraticSpline", QuadraticSpline(u, t)),                       # cache_parameters = false
+        ("QuadraticSpline (cache_parameters)", QuadraticSpline(u, t; cache_parameters = true)),
+    ]
+end
+
+@testset "Thread safety / reentrancy (issue #532)" begin
+    pts = collect(range(1.0, 6.0; length = 41))   # includes the knots themselves
+
+    # Deterministic guard (catches the bug even on a single thread): repeated and
+    # interleaved evaluation must leave the interpolant's scratch buffer — and its
+    # results — unchanged.
+    @testset "evaluation does not mutate the interpolant: $name" for (name, A) in
+        thread_safety_interpolants()
+        sc_before = copy(A.sc)
+        ref = [A(x) for x in pts]
+        dref = [derivative(A, x) for x in pts]
+        for _ in 1:200, x in pts
+            A(x)
+            derivative(A, x)
+        end
+        @test A.sc == sc_before
+        @test [A(x) for x in pts] == ref
+        @test [derivative(A, x) for x in pts] == dref
+    end
+
+    # End-to-end check: many threads hammering one shared interpolant must all
+    # agree with the single-threaded reference. Only meaningful with >1 thread.
+    if Threads.nthreads() > 1
+        @testset "concurrent evaluation is race-free: $name" for (name, A) in
+            thread_safety_interpolants()
+            ref = [A(x) for x in pts]
+            wrong = Threads.Atomic{Int}(0)
+            Threads.@threads for _ in 1:50_000
+                for (j, x) in enumerate(pts)
+                    A(x) == ref[j] || Threads.atomic_add!(wrong, 1)
+                end
+            end
+            @test wrong[] == 0
+        end
+    else
+        @info "issue #532: concurrent stress test skipped; start Julia with more " *
+            "than one thread (e.g. `julia -t auto`) to exercise it."
+    end
+end
diff --git a/test/qa/alloc_tests.jl b/test/qa/alloc_tests.jl
index 4235f843..a5e187b4 100644
--- a/test/qa/alloc_tests.jl
+++ b/test/qa/alloc_tests.jl
@@ -1,8 +1,10 @@
 using DataInterpolations
+using DataInterpolations: derivative
 using AllocCheck: @check_allocs
 using StaticArrays: SVector
 
 @check_allocs(test_allocs(itp, x) = itp(x)) # Reuse function definition to save on compilation time
+@check_allocs(test_derivative_allocs(itp, x) = derivative(itp, x))
 
 @testset "Allocation-free interpolation tests" begin
     @testset "LinearInterpolation" begin
@@ -53,6 +55,32 @@ using StaticArrays: SVector
         @test_nowarn test_allocs(A_s, 0.7)
     end
 
+    # B-spline evaluation must be allocation-free *and* reentrant: it computes the
+    # d+1 nonzero basis functions into a stack buffer instead of a shared scratch
+    # field, which is what makes it thread-safe (#532). Scalar data keeps the return
+    # value itself off the heap.
+    @testset "BSplineInterpolation" begin
+        t = collect(range(0.0, 2π; length = 12))
+        u = sin.(t)
+        for d in (1, 2, 3)
+            A = BSplineInterpolation(
+                u, t, d, :Uniform, :Average; extrapolation = ExtrapolationType.Extension
+            )
+            @test_nowarn test_allocs(A, 3.1)
+            @test_nowarn test_derivative_allocs(A, 3.1)
+        end
+    end
+
+    @testset "BSplineApprox" begin
+        t = collect(range(0.0, 2π; length = 12))
+        u = sin.(t)
+        A = BSplineApprox(
+            u, t, 3, 8, :Uniform, :Average; extrapolation = ExtrapolationType.Extension
+        )
+        @test_nowarn test_allocs(A, 3.1)
+        @test_nowarn test_derivative_allocs(A, 3.1)
+    end
+
     @testset "CubicSpline" begin
         t = [-1.0, 0.0, 1.0]
         u_ = [0.0, 1.0, 3.0]' .* ones(4)