diff --git a/src/Cl2.jl b/src/Cl2.jl
index 9bd43e8..0985234 100644
--- a/src/Cl2.jl
+++ b/src/Cl2.jl
@@ -85,6 +85,7 @@ end
 const _cl2_coeff_cache = Ref{Vector{BigFloat}}(BigFloat[])
 const _cl2_coeff_prec = Ref{Int}(0)
 const _cl2_coeff_K = Ref{Int}(0)
+const _cl2_coeff_lock = ReentrantLock()
 
 
 function _cl2(x::BigFloat)
@@ -161,31 +162,33 @@ end
 function _cl2_ensure_coeffs(K::Int)
     prec = precision(BigFloat)
 
-    if _cl2_coeff_prec[] >= prec && _cl2_coeff_K[] >= K
-        _cl2_coeff_cache[]
-    else
-        Beven = _bernoulli_even(K)
-
-        inv_twopi_sq = 1/(2*BigFloat(pi))^2
-        c = Vector{BigFloat}(undef, K)
-        fac = one(BigFloat)
-        inv_twopi_pow = one(BigFloat)
-
-        for k in 1:K
-            m = 2k
-            fac *= (m - 1)*m # (2k)!
-            inv_twopi_pow *= inv_twopi_sq # 1/(2π)^{2k}
-            # Standard Bernoulli coefficient: |B_{2k}|/(2k(2k+1)(2k)!)
-            orig = abs(Beven[k + 1])/(m*(m + 1)*fac)
-            # Kummer correction: 1/(k(2k+1)(2π)^{2k})
-            corr = inv_twopi_pow/(k*(m + 1))
-            c[k] = orig - corr
+    lock(_cl2_coeff_lock) do
+        if _cl2_coeff_prec[] >= prec && _cl2_coeff_K[] >= K
+            _cl2_coeff_cache[]
+        else
+            Beven = _bernoulli_even(K)
+
+            inv_twopi_sq = 1/(2*BigFloat(pi))^2
+            c = Vector{BigFloat}(undef, K)
+            fac = one(BigFloat)
+            inv_twopi_pow = one(BigFloat)
+
+            for k in 1:K
+                m = 2k
+                fac *= (m - 1)*m # (2k)!
+                inv_twopi_pow *= inv_twopi_sq # 1/(2π)^{2k}
+                # Standard Bernoulli coefficient: |B_{2k}|/(2k(2k+1)(2k)!)
+                orig = abs(Beven[k + 1])/(m*(m + 1)*fac)
+                # Kummer correction: 1/(k(2k+1)(2π)^{2k})
+                corr = inv_twopi_pow/(k*(m + 1))
+                c[k] = orig - corr
+            end
+
+            _cl2_coeff_cache[] = c
+            _cl2_coeff_prec[] = prec
+            _cl2_coeff_K[] = K
+
+            c
         end
-
-        _cl2_coeff_cache[] = c
-        _cl2_coeff_prec[] = prec
-        _cl2_coeff_K[] = K
-
-        c
     end
 end
diff --git a/test/Cl2.jl b/test/Cl2.jl
index 61cfbb0..4831693 100644
--- a/test/Cl2.jl
+++ b/test/Cl2.jl
@@ -72,3 +72,39 @@ end
         @test ClausenFunctions.cl(2,BigFloat("1"))  ≈ BigFloat("1.01395913236076850429457433888591468756117928007771731687704851226813781234607955733638821865477122042157440086434150311308425232269856980238591034316884447292797795707328765622668664434914352893878283") rtol=eps
     end
 end
+
+
+@testset "cl2 BigFloat thread safety" begin
+    # Hammer _cl2_ensure_coeffs from multiple threads with varying
+    # precisions to exercise the cache lock.
+
+    # reference value, calculated sequentially in the main thread
+    ref_256 = setprecision(BigFloat, 256) do
+        ClausenFunctions.cl2(BigFloat(pi)/4)
+    end
+    ref_512 = setprecision(BigFloat, 512) do
+        ClausenFunctions.cl2(BigFloat(pi)/4)
+    end
+
+    # values calculated in parallel in multiple threads
+    results = Vector{Tuple{BigFloat,BigFloat}}(undef, 64)
+    Threads.@threads for i in 1:64
+        r256 = setprecision(BigFloat, 256) do
+            ClausenFunctions.cl2(BigFloat(pi)/4)
+        end
+        r512 = setprecision(BigFloat, 512) do
+            ClausenFunctions.cl2(BigFloat(pi)/4)
+        end
+        results[i] = (r256, r512)
+    end
+
+    # test equality
+    for (r256, r512) in results
+        setprecision(BigFloat, 256) do
+            @test BigFloat(r256) == BigFloat(ref_256)
+        end
+        setprecision(BigFloat, 512) do
+            @test BigFloat(r512) == BigFloat(ref_512)
+        end
+    end
+end