From 589bbbaff1232a2c4675b01a15ae22f31f6cfae0 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Thu, 26 Mar 2026 21:46:08 +0100 Subject: [PATCH 01/11] adding BigFloat overload of cl2 with series acceleration --- src/Cl2.jl | 125 +++++++++++++++++++++++++++++++++++++++++ src/range_reduction.jl | 21 +++++-- test/Cl2.jl | 44 +++++++++++++++ 3 files changed, 184 insertions(+), 6 deletions(-) diff --git a/src/Cl2.jl b/src/Cl2.jl index 42e276f..feb7f7e 100644 --- a/src/Cl2.jl +++ b/src/Cl2.jl @@ -56,3 +56,128 @@ function _cl2(x::Float64)::Float64 sgn*y*p/q end end + + +# BigFloat Cl₂ via Kummer-accelerated Bernoulli series: +# +# Cl₂(θ)/θ = 3 - log[θ(1 - θ²/(4π²))] - (2π/θ)log((2π+θ)/(2π-θ)) +# + Σ_{n=1}^K (ζ(2n)-1)/(n(2n+1)) (θ/(2π))^{2n} +# +# valid for |θ| < 2π (always true after range reduction to [0, π]). +# +# Rewriting in the θ^{2n} basis avoids large intermediate values: +# +# d_n = (ζ(2n)-1)/(n(2n+1)(2π)^{2n}) +# = |B_{2n}|/(2n(2n+1)(2n)!) - 1/(n(2n+1)(2π)^{2n}) +# +# The first term is exactly the standard Bernoulli coefficient (cheap, +# small values); the second is a trivial geometric correction. +# +# Compared to the unaccelerated series, ζ(2n)-1 ≈ 1/4^n roughly halves +# K at the worst case θ = π, while the coefficient computation stays +# as cheap as the original. +# +# Convergence rate: (θ/(4π))² per term (for large n). +# At θ = π (worst case): ~1.2 digits per term → K ≈ 0.83P for P digits. +# At θ = π/4: ~2.4 digits per term → K ≈ 0.42P. + +# Cache: stores d_n for n = 1, ..., K. +const _cl2_coeff_cache = Ref{Vector{BigFloat}}(BigFloat[]) +const _cl2_coeff_prec = Ref{Int}(0) +const _cl2_coeff_K = Ref{Int}(0) + + +function _cl2(x::BigFloat) + (x, sgn) = range_reduce_even(x) + + iszero(x) && return x + x == BigFloat(pi) && return zero(x) + + prec_bits = precision(BigFloat) + target_digits = ceil(Int, prec_bits*0.30103) + 10 + + twopi = 2*BigFloat(pi) + v = x*x + + # ≈ -2log₁₀(θ/(2π)) + 0.6 digits per term (the 0.6 from ζ(2n)-1 ∼ 1/4ⁿ) + ratio = Float64(x)/(2*Float64(pi)) + digits_per_term = max(-2*log10(ratio) + 0.6, 0.5) + K = ceil(Int, target_digits/digits_per_term) + 30 + + c = _cl2_ensure_coeffs(K) + + # Horner: S = v(d₁ + v(d₂ + ... + v dₖ)) + s = c[K] + for k in (K - 1):-1:1 + s = s*v + c[k] + end + s *= v + + # Prefix: 3 - log[θ(1 - θ²/(4π²))] - (2π/θ)log((2π+θ)/(2π-θ)) + # log1p avoids cancellation for small θ where (2π+θ)/(2π-θ) ≈ 1. + prefix = 3 - log(x*(one(x) - v/twopi^2)) - (twopi/x)*log1p(2*x/(twopi - x)) + + sgn*x*(prefix + s) +end + + +# Compute Bernoulli numbers B_0, B_2, ..., B_{2K} via the standard recurrence. +function _bernoulli_even(K::Int) + n = 2K + B = zeros(BigFloat, n + 1) + B[1] = one(BigFloat) + B[2] = -one(BigFloat)/2 + + for m in 2:n + s = zero(BigFloat) + binom = one(BigFloat) + for k in 0:m-1 + s += binom*B[k + 1] + binom = binom*(m + 1 - k)/(k + 1) + end + B[m + 1] = -s/(m + 1) + end + + Beven = Vector{BigFloat}(undef, K + 1) + + for i in 0:K + Beven[i + 1] = B[2i + 1] + end + + Beven +end + + +# Compute and cache d_n = |B_{2n}|/(2n(2n+1)(2n)!) - 1/(n(2n+1)(2π)^{2n}) +# for n = 1,...,K. +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 + end + + _cl2_coeff_cache[] = c + _cl2_coeff_prec[] = prec + _cl2_coeff_K[] = K + + c + end +end diff --git a/src/range_reduction.jl b/src/range_reduction.jl index 6fa414e..75de300 100644 --- a/src/range_reduction.jl +++ b/src/range_reduction.jl @@ -15,17 +15,24 @@ function two_pi_minus(x::Float64) (p0 - x) + p1 end +# Generic fallback for BigFloat and other arbitrary-precision types +function two_pi_minus(x::Real) + 2*oftype(x, pi) - x +end + # returns range-reduced x in [0,pi] for odd n function range_reduce_odd(x::Real) if x < zero(x) x = -x end - if x >= 2pi - x = mod(x, 2pi) + twopi = 2*oftype(x, pi) + + if x >= twopi + x = mod(x, twopi) end - if x > pi + if x > oftype(x, pi) x = two_pi_minus(x) end @@ -41,11 +48,13 @@ function range_reduce_even(x::Real) sgn = -one(x) end - if x >= 2pi - x = mod(x, 2pi) + twopi = 2*oftype(x, pi) + + if x >= twopi + x = mod(x, twopi) end - if x > pi + if x > oftype(x, pi) x = two_pi_minus(x) sgn = -sgn end diff --git a/test/Cl2.jl b/test/Cl2.jl index aa00339..2dae626 100644 --- a/test/Cl2.jl +++ b/test/Cl2.jl @@ -24,3 +24,47 @@ @test !signbit(ClausenFunctions.cl2(0.0)) @test signbit(ClausenFunctions.cl2(-0.0)) end + + +@testset "cl2 BigFloat" begin + # Reference values computed with PolyLog.li(2, cis(θ)) at 1024-bit precision. + + setprecision(BigFloat, 512) do + # Cl₂(π/6) + @test ClausenFunctions.cl2(BigFloat(π)/6) ≈ big"0.86437913105389274962503631519021947866818857640368970418203768977532471558206" rtol=big"1e-60" + + # Cl₂(π/4) + @test ClausenFunctions.cl2(BigFloat(π)/4) ≈ big"0.98187215105020335671792454306019566713079097166073046157661313465315566504977" rtol=big"1e-60" + + # Cl₂(π/3) + @test ClausenFunctions.cl2(BigFloat(π)/3) ≈ big"1.01494160640965362502120255427452028594168930753029979201748910677659747625824" rtol=big"1e-60" + + # Cl₂(π/2) = Catalan constant G + @test ClausenFunctions.cl2(BigFloat(π)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" + + # Cl₂(1) + @test ClausenFunctions.cl2(big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" + + # Cl₂(2π/3) + @test ClausenFunctions.cl2(2*BigFloat(π)/3) ≈ big"0.67662773760643575001413503618301352396112620502019986134499273785106498417216" rtol=big"1e-60" + + # Special values + @test ClausenFunctions.cl2(big"0.0") == big"0.0" + @test ClausenFunctions.cl2(BigFloat(π)) ≈ big"0.0" atol=big"1e-60" + + # Antisymmetry: Cl₂(-x) = -Cl₂(x) + @test ClausenFunctions.cl2(-BigFloat(π)/4) ≈ -ClausenFunctions.cl2(BigFloat(π)/4) rtol=big"1e-60" + + # Periodicity: Cl₂(x + 2π) = Cl₂(x) + @test ClausenFunctions.cl2(big"1" + 2*BigFloat(π)) ≈ ClausenFunctions.cl2(big"1") rtol=big"1e-40" + end + + # Test at higher precision (1000 bits ≈ 300 digits) + setprecision(BigFloat, 1000) do + # Cl₂(π/2) = Catalan constant G + @test ClausenFunctions.cl2(BigFloat(π)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" + + # Cl₂(1) + @test ClausenFunctions.cl2(big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" + end +end From 712bce53ee1bfd0b9016e7392e79612132578c61 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Thu, 26 Mar 2026 21:52:29 +0100 Subject: [PATCH 02/11] use fast BigFloat overload for cl(2,x) --- src/Cl.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/Cl.jl b/src/Cl.jl index 4fb3356..42100ce 100644 --- a/src/Cl.jl +++ b/src/Cl.jl @@ -242,7 +242,9 @@ function _cl(n::Integer, x::Float64)::Float64 end function _cl(n::Integer, x::Real) - if iseven(n) + if n == 2 && typeof(x) == BigFloat + cl2(x) # use faster BigFloat overload for n = 2 + elseif iseven(n) imag(PolyLog.li(n, cis(x))) else real(PolyLog.li(n, cis(x))) From d27c249f4577904fee0f52841c52107392bfe253 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Thu, 26 Mar 2026 21:55:52 +0100 Subject: [PATCH 03/11] adding test for new overload of two_pi_minus --- test/range_reduction.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/range_reduction.jl b/test/range_reduction.jl index 98db305..95bb1d3 100644 --- a/test/range_reduction.jl +++ b/test/range_reduction.jl @@ -47,4 +47,12 @@ end @test f(nextfloat(2pi, 2)) ≈ -2*eps(2pi) atol=2*eps(Float64) @test f(prevfloat(2pi)) ≈ eps(2pi) atol=2*eps(Float64) @test f(prevfloat(2pi, 2)) ≈ 2*eps(2pi) atol=2*eps(Float64) + + # BigFloat + twopi = 2*BigFloat(pi) + @test f(twopi) ≈ zero(BigFloat) atol=2*eps(BigFloat) + @test f(nextfloat(twopi)) ≈ -eps(twopi) atol=2*eps(BigFloat) + @test f(nextfloat(twopi, 2)) ≈ -2*eps(twopi) atol=2*eps(BigFloat) + @test f(prevfloat(twopi)) ≈ eps(twopi) atol=2*eps(BigFloat) + @test f(prevfloat(twopi, 2)) ≈ 2*eps(twopi) atol=2*eps(BigFloat) end From f7398e98ea4fbf2ff4eb3f49ef316d74c6c58c1c Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Thu, 26 Mar 2026 22:10:43 +0100 Subject: [PATCH 04/11] use nextfloat and prevfloat only if defined --- test/range_reduction.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/test/range_reduction.jl b/test/range_reduction.jl index 95bb1d3..0121301 100644 --- a/test/range_reduction.jl +++ b/test/range_reduction.jl @@ -50,9 +50,13 @@ end # BigFloat twopi = 2*BigFloat(pi) - @test f(twopi) ≈ zero(BigFloat) atol=2*eps(BigFloat) - @test f(nextfloat(twopi)) ≈ -eps(twopi) atol=2*eps(BigFloat) - @test f(nextfloat(twopi, 2)) ≈ -2*eps(twopi) atol=2*eps(BigFloat) - @test f(prevfloat(twopi)) ≈ eps(twopi) atol=2*eps(BigFloat) - @test f(prevfloat(twopi, 2)) ≈ 2*eps(twopi) atol=2*eps(BigFloat) + @test f(twopi) ≈ zero(BigFloat) atol=2*eps(BigFloat) + if isdefined(Base, :nextfloat) + @test f(nextfloat(twopi)) ≈ -eps(twopi) atol=2*eps(BigFloat) + @test f(nextfloat(twopi, 2)) ≈ -2*eps(twopi) atol=2*eps(BigFloat) + end + if isdefined(Base, :prevfloat) + @test f(prevfloat(twopi)) ≈ eps(twopi) atol=2*eps(BigFloat) + @test f(prevfloat(twopi, 2)) ≈ 2*eps(twopi) atol=2*eps(BigFloat) + end end From 143f710e5158c4234ecb059600f136d051f725c8 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Thu, 26 Mar 2026 22:18:08 +0100 Subject: [PATCH 05/11] can't use overload of nextfloat with two arguments for julia 1.0 --- test/range_reduction.jl | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/test/range_reduction.jl b/test/range_reduction.jl index 0121301..dbcc922 100644 --- a/test/range_reduction.jl +++ b/test/range_reduction.jl @@ -50,13 +50,9 @@ end # BigFloat twopi = 2*BigFloat(pi) - @test f(twopi) ≈ zero(BigFloat) atol=2*eps(BigFloat) - if isdefined(Base, :nextfloat) - @test f(nextfloat(twopi)) ≈ -eps(twopi) atol=2*eps(BigFloat) - @test f(nextfloat(twopi, 2)) ≈ -2*eps(twopi) atol=2*eps(BigFloat) - end - if isdefined(Base, :prevfloat) - @test f(prevfloat(twopi)) ≈ eps(twopi) atol=2*eps(BigFloat) - @test f(prevfloat(twopi, 2)) ≈ 2*eps(twopi) atol=2*eps(BigFloat) - end + @test f(twopi) ≈ zero(BigFloat) atol=2*eps(BigFloat) + @test f(nextfloat(twopi)) ≈ -eps(twopi) atol=2*eps(BigFloat) + @test f(nextfloat(nextfloat(twopi))) ≈ -2*eps(twopi) atol=2*eps(BigFloat) + @test f(prevfloat(twopi)) ≈ eps(twopi) atol=2*eps(BigFloat) + @test f(prevfloat(prevfloat(twopi))) ≈ 2*eps(twopi) atol=2*eps(BigFloat) end From 465bc233f95aa6ddea0efb9a66ffd1b2e68a90f7 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Thu, 26 Mar 2026 22:39:49 +0100 Subject: [PATCH 06/11] test cl(2,::BigFloat) --- test/Cl2.jl | 54 +++++++++++++++++++++++++---------------------------- 1 file changed, 25 insertions(+), 29 deletions(-) diff --git a/test/Cl2.jl b/test/Cl2.jl index 2dae626..9a3581c 100644 --- a/test/Cl2.jl +++ b/test/Cl2.jl @@ -30,41 +30,37 @@ end # Reference values computed with PolyLog.li(2, cis(θ)) at 1024-bit precision. setprecision(BigFloat, 512) do - # Cl₂(π/6) - @test ClausenFunctions.cl2(BigFloat(π)/6) ≈ big"0.86437913105389274962503631519021947866818857640368970418203768977532471558206" rtol=big"1e-60" - - # Cl₂(π/4) - @test ClausenFunctions.cl2(BigFloat(π)/4) ≈ big"0.98187215105020335671792454306019566713079097166073046157661313465315566504977" rtol=big"1e-60" - - # Cl₂(π/3) - @test ClausenFunctions.cl2(BigFloat(π)/3) ≈ big"1.01494160640965362502120255427452028594168930753029979201748910677659747625824" rtol=big"1e-60" - - # Cl₂(π/2) = Catalan constant G - @test ClausenFunctions.cl2(BigFloat(π)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" - - # Cl₂(1) - @test ClausenFunctions.cl2(big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" - - # Cl₂(2π/3) - @test ClausenFunctions.cl2(2*BigFloat(π)/3) ≈ big"0.67662773760643575001413503618301352396112620502019986134499273785106498417216" rtol=big"1e-60" - - # Special values - @test ClausenFunctions.cl2(big"0.0") == big"0.0" - @test ClausenFunctions.cl2(BigFloat(π)) ≈ big"0.0" atol=big"1e-60" + @test ClausenFunctions.cl2(BigFloat(pi)/6) ≈ big"0.86437913105389274962503631519021947866818857640368970418203768977532471558206" rtol=big"1e-60" + @test ClausenFunctions.cl(2,BigFloat(pi)/6) ≈ big"0.86437913105389274962503631519021947866818857640368970418203768977532471558206" rtol=big"1e-60" + @test ClausenFunctions.cl2(BigFloat(pi)/4) ≈ big"0.98187215105020335671792454306019566713079097166073046157661313465315566504977" rtol=big"1e-60" + @test ClausenFunctions.cl(2,BigFloat(pi)/4) ≈ big"0.98187215105020335671792454306019566713079097166073046157661313465315566504977" rtol=big"1e-60" + @test ClausenFunctions.cl2(BigFloat(pi)/3) ≈ big"1.01494160640965362502120255427452028594168930753029979201748910677659747625824" rtol=big"1e-60" + @test ClausenFunctions.cl(2,BigFloat(pi)/3) ≈ big"1.01494160640965362502120255427452028594168930753029979201748910677659747625824" rtol=big"1e-60" + @test ClausenFunctions.cl2(BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" + @test ClausenFunctions.cl(2,BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" + @test ClausenFunctions.cl2(big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" + @test ClausenFunctions.cl(2,big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" + @test ClausenFunctions.cl2(2*BigFloat(pi)/3) ≈ big"0.67662773760643575001413503618301352396112620502019986134499273785106498417216" rtol=big"1e-60" + @test ClausenFunctions.cl(2,2*BigFloat(pi)/3) ≈ big"0.67662773760643575001413503618301352396112620502019986134499273785106498417216" rtol=big"1e-60" + @test ClausenFunctions.cl2(big"0.0") == big"0.0" + @test ClausenFunctions.cl(2,big"0.0") == big"0.0" + @test ClausenFunctions.cl2(BigFloat(pi)) ≈ big"0.0" atol=big"1e-60" + @test ClausenFunctions.cl(2,BigFloat(pi)) ≈ big"0.0" atol=big"1e-60" # Antisymmetry: Cl₂(-x) = -Cl₂(x) - @test ClausenFunctions.cl2(-BigFloat(π)/4) ≈ -ClausenFunctions.cl2(BigFloat(π)/4) rtol=big"1e-60" + @test ClausenFunctions.cl2(-BigFloat(pi)/4) ≈ -ClausenFunctions.cl2(BigFloat(pi)/4) rtol=big"1e-60" + @test ClausenFunctions.cl(2,-BigFloat(pi)/4) ≈ -ClausenFunctions.cl2(BigFloat(pi)/4) rtol=big"1e-60" - # Periodicity: Cl₂(x + 2π) = Cl₂(x) - @test ClausenFunctions.cl2(big"1" + 2*BigFloat(π)) ≈ ClausenFunctions.cl2(big"1") rtol=big"1e-40" + # Periodicity: Cl₂(x + 2pi) = Cl₂(x) + @test ClausenFunctions.cl2(big"1" + 2*BigFloat(pi)) ≈ ClausenFunctions.cl2(big"1") rtol=big"1e-40" + @test ClausenFunctions.cl(2,big"1" + 2*BigFloat(pi)) ≈ ClausenFunctions.cl2(big"1") rtol=big"1e-40" end # Test at higher precision (1000 bits ≈ 300 digits) setprecision(BigFloat, 1000) do - # Cl₂(π/2) = Catalan constant G - @test ClausenFunctions.cl2(BigFloat(π)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" - - # Cl₂(1) - @test ClausenFunctions.cl2(big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" + @test ClausenFunctions.cl2(BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" + @test ClausenFunctions.cl(2,BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" + @test ClausenFunctions.cl2(big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" + @test ClausenFunctions.cl(2,big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" end end From 9c141939d06e91cb5c304d09d525c48fa11f6834 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Fri, 27 Mar 2026 00:27:00 +0100 Subject: [PATCH 07/11] adapt test precision for cl2 for BigFloat --- test/Cl2.jl | 60 ++++++++++++++++++++++++++++++----------------------- 1 file changed, 34 insertions(+), 26 deletions(-) diff --git a/test/Cl2.jl b/test/Cl2.jl index 9a3581c..29fde80 100644 --- a/test/Cl2.jl +++ b/test/Cl2.jl @@ -27,40 +27,48 @@ end @testset "cl2 BigFloat" begin - # Reference values computed with PolyLog.li(2, cis(θ)) at 1024-bit precision. + decimal_digits = 100 + binary_digits = ceil(Int, decimal_digits*log(10)/log(2)) - setprecision(BigFloat, 512) do - @test ClausenFunctions.cl2(BigFloat(pi)/6) ≈ big"0.86437913105389274962503631519021947866818857640368970418203768977532471558206" rtol=big"1e-60" - @test ClausenFunctions.cl(2,BigFloat(pi)/6) ≈ big"0.86437913105389274962503631519021947866818857640368970418203768977532471558206" rtol=big"1e-60" - @test ClausenFunctions.cl2(BigFloat(pi)/4) ≈ big"0.98187215105020335671792454306019566713079097166073046157661313465315566504977" rtol=big"1e-60" - @test ClausenFunctions.cl(2,BigFloat(pi)/4) ≈ big"0.98187215105020335671792454306019566713079097166073046157661313465315566504977" rtol=big"1e-60" - @test ClausenFunctions.cl2(BigFloat(pi)/3) ≈ big"1.01494160640965362502120255427452028594168930753029979201748910677659747625824" rtol=big"1e-60" - @test ClausenFunctions.cl(2,BigFloat(pi)/3) ≈ big"1.01494160640965362502120255427452028594168930753029979201748910677659747625824" rtol=big"1e-60" - @test ClausenFunctions.cl2(BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" - @test ClausenFunctions.cl(2,BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" - @test ClausenFunctions.cl2(big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" - @test ClausenFunctions.cl(2,big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" - @test ClausenFunctions.cl2(2*BigFloat(pi)/3) ≈ big"0.67662773760643575001413503618301352396112620502019986134499273785106498417216" rtol=big"1e-60" - @test ClausenFunctions.cl(2,2*BigFloat(pi)/3) ≈ big"0.67662773760643575001413503618301352396112620502019986134499273785106498417216" rtol=big"1e-60" + setprecision(BigFloat, binary_digits) do + eps = 10.0^(23 - decimal_digits) # @todo: why are we loosing so many digits? + + @test ClausenFunctions.cl2(BigFloat(pi)/6) ≈ big"0.8643791310538927496250363151902194786681885764036897041820376897753247155820641851870219305078075779" rtol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)/6) ≈ big"0.8643791310538927496250363151902194786681885764036897041820376897753247155820641851870219305078075779" rtol=eps + @test ClausenFunctions.cl2(BigFloat(pi)/4) ≈ big"0.9818721510502033567179245430601956671307909716607304615766131346531556650497696362249028028843877241" rtol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)/4) ≈ big"0.9818721510502033567179245430601956671307909716607304615766131346531556650497696362249028028843877241" rtol=eps + @test ClausenFunctions.cl2(BigFloat(pi)/3) ≈ big"1.0149416064096536250212025542745202859416893075302997920174891067765974762582440221364703542282566949" rtol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)/3) ≈ big"1.0149416064096536250212025542745202859416893075302997920174891067765974762582440221364703542282566949" rtol=eps + @test ClausenFunctions.cl2(BigFloat(pi)/2) ≈ big"0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062" rtol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)/2) ≈ big"0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062" rtol=eps + @test ClausenFunctions.cl2(big"1") ≈ big"1.0139591323607685042945743388859146875611792800777173168770485122681378123460795573363882186547712204" rtol=eps + @test ClausenFunctions.cl(2,big"1") ≈ big"1.0139591323607685042945743388859146875611792800777173168770485122681378123460795573363882186547712204" rtol=eps + @test ClausenFunctions.cl2(2*BigFloat(pi)/3) ≈ big"0.6766277376064357500141350361830135239611262050201998613449927378510649841721626814243135694855044633" rtol=eps + @test ClausenFunctions.cl(2,2*BigFloat(pi)/3) ≈ big"0.6766277376064357500141350361830135239611262050201998613449927378510649841721626814243135694855044633" rtol=eps @test ClausenFunctions.cl2(big"0.0") == big"0.0" @test ClausenFunctions.cl(2,big"0.0") == big"0.0" - @test ClausenFunctions.cl2(BigFloat(pi)) ≈ big"0.0" atol=big"1e-60" - @test ClausenFunctions.cl(2,BigFloat(pi)) ≈ big"0.0" atol=big"1e-60" + @test ClausenFunctions.cl2(BigFloat(pi)) ≈ big"0.0" atol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)) ≈ big"0.0" atol=eps # Antisymmetry: Cl₂(-x) = -Cl₂(x) - @test ClausenFunctions.cl2(-BigFloat(pi)/4) ≈ -ClausenFunctions.cl2(BigFloat(pi)/4) rtol=big"1e-60" - @test ClausenFunctions.cl(2,-BigFloat(pi)/4) ≈ -ClausenFunctions.cl2(BigFloat(pi)/4) rtol=big"1e-60" + @test ClausenFunctions.cl2(-BigFloat(pi)/4) ≈ -ClausenFunctions.cl2(BigFloat(pi)/4) rtol=eps + @test ClausenFunctions.cl(2,-BigFloat(pi)/4) ≈ -ClausenFunctions.cl2(BigFloat(pi)/4) rtol=eps # Periodicity: Cl₂(x + 2pi) = Cl₂(x) - @test ClausenFunctions.cl2(big"1" + 2*BigFloat(pi)) ≈ ClausenFunctions.cl2(big"1") rtol=big"1e-40" - @test ClausenFunctions.cl(2,big"1" + 2*BigFloat(pi)) ≈ ClausenFunctions.cl2(big"1") rtol=big"1e-40" + @test ClausenFunctions.cl2(big"1" + 2*BigFloat(pi)) ≈ ClausenFunctions.cl2(big"1") rtol=eps + @test ClausenFunctions.cl(2,big"1" + 2*BigFloat(pi)) ≈ ClausenFunctions.cl2(big"1") rtol=eps end - # Test at higher precision (1000 bits ≈ 300 digits) - setprecision(BigFloat, 1000) do - @test ClausenFunctions.cl2(BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" - @test ClausenFunctions.cl(2,BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625" rtol=big"1e-60" - @test ClausenFunctions.cl2(big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" - @test ClausenFunctions.cl(2,big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234608" rtol=big"1e-60" + # Test at higher precision + decimal_digits = 200 + binary_digits = ceil(Int, decimal_digits*log(10)/log(2)) + + setprecision(BigFloat, binary_digits) do + eps = 10.0^(125 - decimal_digits) # @todo: why are we loosing so many digits? + + @test ClausenFunctions.cl2(BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625476947935651292611510624857442261919619957903589880332585905943159473748115840699533202877331946051903872747816408786590902" rtol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625476947935651292611510624857442261919619957903589880332585905943159473748115840699533202877331946051903872747816408786590902" rtol=eps + @test ClausenFunctions.cl2(big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234607955733638821865477122042157440086434150311308425232269856980238591034316884447292797795707328765622668664434914352893878283" rtol=eps + @test ClausenFunctions.cl(2,big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234607955733638821865477122042157440086434150311308425232269856980238591034316884447292797795707328765622668664434914352893878283" rtol=eps end end From ecc9c6397ffd3648fb1f0b49b6cf6c406ddda7f9 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Fri, 27 Mar 2026 08:15:14 +0100 Subject: [PATCH 08/11] use muladd for improved speed and accuarcy --- src/Cl2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cl2.jl b/src/Cl2.jl index feb7f7e..7d726a5 100644 --- a/src/Cl2.jl +++ b/src/Cl2.jl @@ -109,7 +109,7 @@ function _cl2(x::BigFloat) # Horner: S = v(d₁ + v(d₂ + ... + v dₖ)) s = c[K] for k in (K - 1):-1:1 - s = s*v + c[k] + s = muladd(s, v, c[k]) end s *= v From e2a640a44b0e1838351c15948bac8e301382d033 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Fri, 27 Mar 2026 08:23:04 +0100 Subject: [PATCH 09/11] another use of muladd --- src/Cl2.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Cl2.jl b/src/Cl2.jl index 7d726a5..ef1bf4e 100644 --- a/src/Cl2.jl +++ b/src/Cl2.jl @@ -132,7 +132,7 @@ function _bernoulli_even(K::Int) s = zero(BigFloat) binom = one(BigFloat) for k in 0:m-1 - s += binom*B[k + 1] + s = binom*B[k + 1] + s binom = binom*(m + 1 - k)/(k + 1) end B[m + 1] = -s/(m + 1) From 9df7528766d9e8831ec886e56c69182be9baa2c6 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Fri, 27 Mar 2026 18:17:48 +0100 Subject: [PATCH 10/11] don't use big"" notation as it does not respect setprecision(...) --- test/Cl2.jl | 48 ++++++++++++++++++++++++------------------------ 1 file changed, 24 insertions(+), 24 deletions(-) diff --git a/test/Cl2.jl b/test/Cl2.jl index 29fde80..61cfbb0 100644 --- a/test/Cl2.jl +++ b/test/Cl2.jl @@ -31,32 +31,32 @@ end binary_digits = ceil(Int, decimal_digits*log(10)/log(2)) setprecision(BigFloat, binary_digits) do - eps = 10.0^(23 - decimal_digits) # @todo: why are we loosing so many digits? + eps = 10.0^(1 - decimal_digits) - @test ClausenFunctions.cl2(BigFloat(pi)/6) ≈ big"0.8643791310538927496250363151902194786681885764036897041820376897753247155820641851870219305078075779" rtol=eps - @test ClausenFunctions.cl(2,BigFloat(pi)/6) ≈ big"0.8643791310538927496250363151902194786681885764036897041820376897753247155820641851870219305078075779" rtol=eps - @test ClausenFunctions.cl2(BigFloat(pi)/4) ≈ big"0.9818721510502033567179245430601956671307909716607304615766131346531556650497696362249028028843877241" rtol=eps - @test ClausenFunctions.cl(2,BigFloat(pi)/4) ≈ big"0.9818721510502033567179245430601956671307909716607304615766131346531556650497696362249028028843877241" rtol=eps - @test ClausenFunctions.cl2(BigFloat(pi)/3) ≈ big"1.0149416064096536250212025542745202859416893075302997920174891067765974762582440221364703542282566949" rtol=eps - @test ClausenFunctions.cl(2,BigFloat(pi)/3) ≈ big"1.0149416064096536250212025542745202859416893075302997920174891067765974762582440221364703542282566949" rtol=eps - @test ClausenFunctions.cl2(BigFloat(pi)/2) ≈ big"0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062" rtol=eps - @test ClausenFunctions.cl(2,BigFloat(pi)/2) ≈ big"0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062" rtol=eps - @test ClausenFunctions.cl2(big"1") ≈ big"1.0139591323607685042945743388859146875611792800777173168770485122681378123460795573363882186547712204" rtol=eps - @test ClausenFunctions.cl(2,big"1") ≈ big"1.0139591323607685042945743388859146875611792800777173168770485122681378123460795573363882186547712204" rtol=eps - @test ClausenFunctions.cl2(2*BigFloat(pi)/3) ≈ big"0.6766277376064357500141350361830135239611262050201998613449927378510649841721626814243135694855044633" rtol=eps - @test ClausenFunctions.cl(2,2*BigFloat(pi)/3) ≈ big"0.6766277376064357500141350361830135239611262050201998613449927378510649841721626814243135694855044633" rtol=eps - @test ClausenFunctions.cl2(big"0.0") == big"0.0" - @test ClausenFunctions.cl(2,big"0.0") == big"0.0" - @test ClausenFunctions.cl2(BigFloat(pi)) ≈ big"0.0" atol=eps - @test ClausenFunctions.cl(2,BigFloat(pi)) ≈ big"0.0" atol=eps + @test ClausenFunctions.cl2(BigFloat(pi)/6) ≈ BigFloat("0.8643791310538927496250363151902194786681885764036897041820376897753247155820641851870219305078075779") rtol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)/6) ≈ BigFloat("0.8643791310538927496250363151902194786681885764036897041820376897753247155820641851870219305078075779") rtol=eps + @test ClausenFunctions.cl2(BigFloat(pi)/4) ≈ BigFloat("0.9818721510502033567179245430601956671307909716607304615766131346531556650497696362249028028843877241") rtol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)/4) ≈ BigFloat("0.9818721510502033567179245430601956671307909716607304615766131346531556650497696362249028028843877241") rtol=eps + @test ClausenFunctions.cl2(BigFloat(pi)/3) ≈ BigFloat("1.0149416064096536250212025542745202859416893075302997920174891067765974762582440221364703542282566949") rtol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)/3) ≈ BigFloat("1.0149416064096536250212025542745202859416893075302997920174891067765974762582440221364703542282566949") rtol=eps + @test ClausenFunctions.cl2(BigFloat(pi)/2) ≈ BigFloat("0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062") rtol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)/2) ≈ BigFloat("0.9159655941772190150546035149323841107741493742816721342664981196217630197762547694793565129261151062") rtol=eps + @test ClausenFunctions.cl2(BigFloat("1")) ≈ BigFloat("1.0139591323607685042945743388859146875611792800777173168770485122681378123460795573363882186547712204") rtol=eps + @test ClausenFunctions.cl(2,BigFloat("1")) ≈ BigFloat("1.0139591323607685042945743388859146875611792800777173168770485122681378123460795573363882186547712204") rtol=eps + @test ClausenFunctions.cl2(2*BigFloat(pi)/3) ≈ BigFloat("0.6766277376064357500141350361830135239611262050201998613449927378510649841721626814243135694855044633") rtol=eps + @test ClausenFunctions.cl(2,2*BigFloat(pi)/3) ≈ BigFloat("0.6766277376064357500141350361830135239611262050201998613449927378510649841721626814243135694855044633") rtol=eps + @test ClausenFunctions.cl2(BigFloat("0.0")) == BigFloat("0.0") + @test ClausenFunctions.cl(2,BigFloat("0.0")) == BigFloat("0.0") + @test ClausenFunctions.cl2(BigFloat(pi)) ≈ BigFloat("0.0") atol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)) ≈ BigFloat("0.0") atol=eps # Antisymmetry: Cl₂(-x) = -Cl₂(x) @test ClausenFunctions.cl2(-BigFloat(pi)/4) ≈ -ClausenFunctions.cl2(BigFloat(pi)/4) rtol=eps @test ClausenFunctions.cl(2,-BigFloat(pi)/4) ≈ -ClausenFunctions.cl2(BigFloat(pi)/4) rtol=eps # Periodicity: Cl₂(x + 2pi) = Cl₂(x) - @test ClausenFunctions.cl2(big"1" + 2*BigFloat(pi)) ≈ ClausenFunctions.cl2(big"1") rtol=eps - @test ClausenFunctions.cl(2,big"1" + 2*BigFloat(pi)) ≈ ClausenFunctions.cl2(big"1") rtol=eps + @test ClausenFunctions.cl2(BigFloat("1") + 2*BigFloat(pi)) ≈ ClausenFunctions.cl2(BigFloat("1")) rtol=eps + @test ClausenFunctions.cl(2,BigFloat("1") + 2*BigFloat(pi)) ≈ ClausenFunctions.cl2(BigFloat("1")) rtol=eps end # Test at higher precision @@ -64,11 +64,11 @@ end binary_digits = ceil(Int, decimal_digits*log(10)/log(2)) setprecision(BigFloat, binary_digits) do - eps = 10.0^(125 - decimal_digits) # @todo: why are we loosing so many digits? + eps = 10.0^(1 - decimal_digits) - @test ClausenFunctions.cl2(BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625476947935651292611510624857442261919619957903589880332585905943159473748115840699533202877331946051903872747816408786590902" rtol=eps - @test ClausenFunctions.cl(2,BigFloat(pi)/2) ≈ big"0.91596559417721901505460351493238411077414937428167213426649811962176301977625476947935651292611510624857442261919619957903589880332585905943159473748115840699533202877331946051903872747816408786590902" rtol=eps - @test ClausenFunctions.cl2(big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234607955733638821865477122042157440086434150311308425232269856980238591034316884447292797795707328765622668664434914352893878283" rtol=eps - @test ClausenFunctions.cl(2,big"1") ≈ big"1.01395913236076850429457433888591468756117928007771731687704851226813781234607955733638821865477122042157440086434150311308425232269856980238591034316884447292797795707328765622668664434914352893878283" rtol=eps + @test ClausenFunctions.cl2(BigFloat(pi)/2) ≈ BigFloat("0.91596559417721901505460351493238411077414937428167213426649811962176301977625476947935651292611510624857442261919619957903589880332585905943159473748115840699533202877331946051903872747816408786590902") rtol=eps + @test ClausenFunctions.cl(2,BigFloat(pi)/2) ≈ BigFloat("0.91596559417721901505460351493238411077414937428167213426649811962176301977625476947935651292611510624857442261919619957903589880332585905943159473748115840699533202877331946051903872747816408786590902") rtol=eps + @test ClausenFunctions.cl2(BigFloat("1")) ≈ BigFloat("1.01395913236076850429457433888591468756117928007771731687704851226813781234607955733638821865477122042157440086434150311308425232269856980238591034316884447292797795707328765622668664434914352893878283") rtol=eps + @test ClausenFunctions.cl(2,BigFloat("1")) ≈ BigFloat("1.01395913236076850429457433888591468756117928007771731687704851226813781234607955733638821865477122042157440086434150311308425232269856980238591034316884447292797795707328765622668664434914352893878283") rtol=eps end end From 05c499e9b3e47ccb7e5ad1a5eba115fdbae3abc4 Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Sun, 29 Mar 2026 14:20:46 +0200 Subject: [PATCH 11/11] adding comments --- src/Cl2.jl | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/src/Cl2.jl b/src/Cl2.jl index ef1bf4e..9bd43e8 100644 --- a/src/Cl2.jl +++ b/src/Cl2.jl @@ -93,15 +93,23 @@ function _cl2(x::BigFloat) iszero(x) && return x x == BigFloat(pi) && return zero(x) - prec_bits = precision(BigFloat) - target_digits = ceil(Int, prec_bits*0.30103) + 10 + # +10: guard digits to absorb rounding in the Horner loop and prefix. + target_digits = ceil(Int, precision(BigFloat)*log10(2)) + 10 twopi = 2*BigFloat(pi) v = x*x - # ≈ -2log₁₀(θ/(2π)) + 0.6 digits per term (the 0.6 from ζ(2n)-1 ∼ 1/4ⁿ) + # Each series term decays by a factor of u = (θ/(2π))², giving + # -2·log₁₀(θ/(2π)) decimal digits per term. The Kummer acceleration + # replaces ζ(2n) with ζ(2n)-1 ≈ 1/2^{2n} = 1/4ⁿ, contributing an + # extra log₁₀(4) ≈ 0.602 digits per term (rounded down to 0.6). + # The floor of 0.5 prevents division by zero for θ near 2π (which + # cannot occur after range reduction to [0, π], but is defensive). ratio = Float64(x)/(2*Float64(pi)) digits_per_term = max(-2*log10(ratio) + 0.6, 0.5) + # +30: safety margin covering (a) the first ~10 terms where + # ζ(2n)-1 ≫ 1/4ⁿ so convergence is slower than the asymptotic + # estimate, and (b) any rounding in the Float64 K estimate. K = ceil(Int, target_digits/digits_per_term) + 30 c = _cl2_ensure_coeffs(K)