From 3aaee003b5fab028eb0df5dcffb45e9a175d4d6e Mon Sep 17 00:00:00 2001
From: Alexander Voigt <ahgvoigt@gmail.com>
Date: Sun, 29 Mar 2026 19:52:45 +0200
Subject: [PATCH 1/2] Make cache for cl2 thread-safe

---
 src/Cl2.jl  | 53 ++++++++++++++++++++++++++++-------------------------
 test/Cl2.jl | 28 ++++++++++++++++++++++++++++
 2 files changed, 56 insertions(+), 25 deletions(-)

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..92ca738 100644
--- a/test/Cl2.jl
+++ b/test/Cl2.jl
@@ -72,3 +72,31 @@ 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.
+    ref_256 = setprecision(BigFloat, 256) do
+        ClausenFunctions.cl2(BigFloat(pi)/4)
+    end
+    ref_512 = setprecision(BigFloat, 512) do
+        ClausenFunctions.cl2(BigFloat(pi)/4)
+    end
+
+    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
+
+    for (r256, r512) in results
+        @test BigFloat(r256, precision=256) == BigFloat(ref_256, precision=256)
+        @test BigFloat(r512, precision=512) == BigFloat(ref_512, precision=512)
+    end
+end

From 09c3e20eccd5c59f4c9c26f980f404cc8922347e Mon Sep 17 00:00:00 2001
From: Alexander Voigt <ahgvoigt@gmail.com>
Date: Sun, 29 Mar 2026 20:05:04 +0200
Subject: [PATCH 2/2] use julia 1.0 friendly syntax

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

diff --git a/test/Cl2.jl b/test/Cl2.jl
index 92ca738..4831693 100644
--- a/test/Cl2.jl
+++ b/test/Cl2.jl
@@ -77,6 +77,8 @@ 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
@@ -84,6 +86,7 @@ end
         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
@@ -95,8 +98,13 @@ end
         results[i] = (r256, r512)
     end
 
+    # test equality
     for (r256, r512) in results
-        @test BigFloat(r256, precision=256) == BigFloat(ref_256, precision=256)
-        @test BigFloat(r512, precision=512) == BigFloat(ref_512, precision=512)
+        setprecision(BigFloat, 256) do
+            @test BigFloat(r256) == BigFloat(ref_256)
+        end
+        setprecision(BigFloat, 512) do
+            @test BigFloat(r512) == BigFloat(ref_512)
+        end
     end
 end