Adding faster BigFloat overload of cl2

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)))

src/Cl2.jl

@@ -56,3 +56,136 @@ 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)
+
+    # +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
+
+    # 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)
+
+    # Horner: S = v(d₁ + v(d₂ + ... + v dₖ))
+    s = c[K]
+    for k in (K - 1):-1:1
+         s = muladd(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] + s
+            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

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

test/Cl2.jl

@@ -24,3 +24,51 @@
     @test !signbit(ClausenFunctions.cl2(0.0))
     @test signbit(ClausenFunctions.cl2(-0.0))
 end
+
+
+@testset "cl2 BigFloat" begin
+    decimal_digits = 100
+    binary_digits = ceil(Int, decimal_digits*log(10)/log(2))
+
+    setprecision(BigFloat, binary_digits) do
+        eps = 10.0^(1 - decimal_digits)
+
+        @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(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
+    decimal_digits = 200
+    binary_digits = ceil(Int, decimal_digits*log(10)/log(2))
+
+    setprecision(BigFloat, binary_digits) do
+        eps = 10.0^(1 - decimal_digits)
+
+        @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

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(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