Make DoubleFloats, Quadmath & GenericFFT weak dependencies (package extensions) (#14)

Float64-only users no longer load DoubleFloats, Quadmath (and thus the libquadmath C
library), or GenericFFT: all three are `[weakdeps]` activated via package extensions,
pulled in only when the user does `using DoubleFloats` / `using Quadmath` / `using
GenericFFT`. Verified: `using TransitionMatrices` alone loads none of the three
extensions, and a Float64 EBCM computation runs unchanged.

These are all opt-in extras — the core (Float64, FFTW, Arblib) path never needs them:
- ext/TransitionMatricesDoubleFloatsExt.jl: Double64 — precision(Complex{Double64}),
  Arblib.set!(::Double64), and the Double64 `cbrt`/`∛` workaround.
- ext/TransitionMatricesQuadmathExt.jl: Float128 — Float128⇄Arb conversions,
  precision(ComplexF128), the instance precision(::Float128) fix (Quadmath leaves it
  unimplemented), and Arblib.set!(::Float128).
- ext/TransitionMatricesGenericFFTExt.jl: the generic-type azimuthal FFT for the n-fold
  IITM — `_iitm_fft_capable(Complex{<:AbstractFloat})`, the GenericFFT plan, and the
  generic `_apply_forward_dft!`. Without GenericFFT, Complex{Double64}/Complex{BigFloat}
  n-fold IITM gracefully falls back to the direct azimuthal sum (correct, just no FFT
  acceleration); ComplexF64 stays on FFTW and Acb on Arblib.dft! (both hard deps).

Note: weakening Quadmath alone would have no effect, because DoubleFloats lists Quadmath
in its own deps and so transitively loads it — they had to be weakened together for
Float64-only loads to drop libquadmath.

Arblib stays a hard dependency: the arbitrary-precision + special-function backbone
(gausslegendre for non-Float64, factorial(n>20) via gamma! in the Wigner-d functions —
reached even on the Float64 path — and Arb/Acb as first-class numeric types in the EBCM
Arb-matrix inversion and the IITM Acb FFT backend).

Re-exports of Double64 / Float128 / ComplexF128 are dropped (use the respective package);
Arb/Acb stay exported. Tests and the benchmark load the weak deps explicitly;
test/Project.toml and benchmark/Project.toml gain them. Full suite: 48,290 tests pass.

Author: lucifer1004

Project.toml

@@ -1,25 +1,32 @@
 name = "TransitionMatrices"
 uuid = "057c4241-e127-4181-840e-6b4b92e6eef5"
-version = "0.5.0"
+version = "0.6.0"
 authors = ["Gabriel Wu <wuzihua@pku.edu.cn> and contributors"]
 
 [deps]
 Arblib = "fb37089c-8514-4489-9461-98f9c8763369"
-DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
-GenericFFT = "a8297547-1b15-4a5a-a998-a2ac5f1cef28"
 GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
-Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 StableTasks = "91464d47-22a1-43fe-8b7f-2d57ee82463f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TestItems = "1c621080-faea-4a02-84b6-bbd5e436b8fe"
 Wigxjpf = "af901252-fd8a-4391-8647-10b4fde07a1e"
 
+[weakdeps]
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
+GenericFFT = "a8297547-1b15-4a5a-a998-a2ac5f1cef28"
+Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
+
+[extensions]
+TransitionMatricesDoubleFloatsExt = "DoubleFloats"
+TransitionMatricesGenericFFTExt = "GenericFFT"
+TransitionMatricesQuadmathExt = "Quadmath"
+
 [compat]
 Arblib = "1"
 DoubleFloats = "1"

benchmark/Project.toml

@@ -1,7 +1,9 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 TransitionMatrices = "057c4241-e127-4181-840e-6b4b92e6eef5"
 
 [compat]
 BenchmarkTools = "1"
+DoubleFloats = "1"
 julia = "1.10"

benchmark/benchmarks.jl

@@ -19,6 +19,7 @@
 
 using BenchmarkTools
 using TransitionMatrices
+using DoubleFloats: Double64   # weak dep of TransitionMatrices; load explicitly for the precision benchmark
 
 const SUITE = BenchmarkGroup()
 

ext/TransitionMatricesDoubleFloatsExt.jl

@@ -0,0 +1,42 @@
+module TransitionMatricesDoubleFloatsExt
+
+# Double64 (DoubleFloats) support, loaded on demand. DoubleFloats is a *weak* dependency:
+# the core package never needs Double64 (it is a purely user-opt-in numeric type), so the
+# Double64-specific shims live here and activate only when the user has `using DoubleFloats`.
+# (Arblib is a hard dependency, so it is available to this extension.)
+
+using TransitionMatrices
+using DoubleFloats: Double64
+using Arblib: Arblib
+
+Base.precision(::Type{Complex{Double64}}) = 106
+
+function Arblib.set!(arb::Arblib.ArbLike, val::Double64)
+    Arblib.set!(arb, BigFloat(val))
+end
+
+# Workaround for DoubleFloats v1.9.x: `cbrt` (and thus its alias `∛`) on a `Double64`
+# returns the raw `(hi, lo)` component tuple instead of a `Double64` (the upstream
+# `cbrt_db_db` forgets to wrap its result; `sqrt` is unaffected). This breaks every `∛`
+# call reached with a `Double64` argument, e.g. `volume_equivalent_radius` and the
+# Riccati-Bessel term estimators.
+#
+# The fix is behavior-gated: the guard below calls the upstream `cbrt` before our method
+# exists, so it installs the workaround only while upstream is still broken. Once a fixed
+# DoubleFloats returns a `Double64` here, this block is skipped at load time and the
+# upstream method is used unchanged — no sticky shadowing, no dead code.
+if !(cbrt(Double64(8.0)) isa Double64)
+    # Specialized to `Double64` (more specific than upstream's parametric method);
+    # recovers full precision with two Newton steps from a Float64 seed. `cbrt(Float64(a))`
+    # dispatches to Base, so there is no recursion.
+    function Base.cbrt(x::Double64)
+        iszero(x) && return x
+        a = abs(x)
+        y = Double64(cbrt(Float64(a)))
+        y = (2y + a / (y * y)) / 3
+        y = (2y + a / (y * y)) / 3
+        return x < 0 ? -y : y
+    end
+end
+
+end

ext/TransitionMatricesGenericFFTExt.jl

@@ -0,0 +1,27 @@
+module TransitionMatricesGenericFFTExt
+
+# Generic-type azimuthal FFT for the n-fold IITM, loaded on demand. GenericFFT is a *weak*
+# dependency: ComplexF64 uses FFTW and Acb uses Arblib.dft! (both hard deps), so only the
+# Complex{Double64} / Complex{BigFloat} FFT path needs GenericFFT. Without it, those types
+# fall back to the direct azimuthal sum (slower but correct) — see `_iitm_fft_capable`.
+
+using TransitionMatrices
+using GenericFFT: plan_fft
+
+# Make non-Float64 complex-float types FFT-capable. ComplexF64 keeps its more-specific core
+# method (FFTW); Acb keeps Arblib.dft!.
+TransitionMatrices._iitm_fft_capable(::Type{<:Complex{<:AbstractFloat}}) = true
+
+# Fresh GenericFFT plan each call: a BigFloat plan bakes in precision-specific twiddles and
+# must not be cached/reused across precision changes (no FFTW flags here).
+function TransitionMatrices._azimuthal_fft_plan(::Type{CT}, Nφ, Nϑ) where {CT <: Complex{<:AbstractFloat}}
+    return plan_fft(zeros(CT, Nφ, Nϑ), 1)
+end
+
+# GenericFFT's plan supports `plan * A` (allocating), not in-place mul!.
+function TransitionMatrices._apply_forward_dft!(spectrum::Matrix{CT}, contrast::Matrix{CT},
+        plan) where {CT <: Complex{<:AbstractFloat}}
+    spectrum .= plan * contrast
+end
+
+end

ext/TransitionMatricesQuadmathExt.jl

@@ -0,0 +1,25 @@
+module TransitionMatricesQuadmathExt
+
+# Float128 (Quadmath) support, loaded on demand. Quadmath is a *weak* dependency: the
+# core package never needs Float128 (it is a purely user-opt-in numeric type), so the
+# Float128 ⇄ Arb glue lives here and activates only when the user has `using Quadmath`.
+# (Arblib is a hard dependency, so it is available to this extension.)
+
+using TransitionMatrices
+using Quadmath: Quadmath, Float128, ComplexF128
+using Arblib: Arblib, ArbLike
+
+Base.convert(::Type{Float128}, x::ArbLike) = Float128(BigFloat(x))
+Quadmath.Float128(x::ArbLike) = Float128(BigFloat(x))
+Base.precision(::Type{ComplexF128}) = 113
+
+# Quadmath implements `precision(::Type{Float128})` but not the *instance* method
+# `precision(::Float128)` (it routes to a missing `_precision_with_base_2`). Define it so
+# generic `precision(x)` works on Float128 values (Float128 has a 113-bit significand).
+Base.precision(::Float128) = 113
+
+function Arblib.set!(arb::Arblib.ArbLike, val::Float128)
+    Arblib.set!(arb, BigFloat(val))
+end
+
+end

packages/EBCMPrecisionLossEstimators/Project.toml

@@ -11,5 +11,5 @@ TransitionMatrices = "057c4241-e127-4181-840e-6b4b92e6eef5"
 [compat]
 MLJ = "0.23"
 MLJXGBoostInterface = "0.3"
-TransitionMatrices = "0.5"
+TransitionMatrices = "0.6"
 julia = "1.8"

src/IITM/fourier.jl

@@ -6,9 +6,12 @@ _iitm_ldiv(𝐌::AbstractMatrix{<:Union{Arb, Acb}}, 𝐗) = inv(𝐌) * 𝐗
 _iitm_ldiv(𝐌, 𝐗) = 𝐌 \ 𝐗
 
 # Capability predicate: true when the FFT path is available for this complex type.
-_iitm_fft_capable(::Type{<:Complex{<:AbstractFloat}}) = true  # FFTW / GenericFFT
-_iitm_fft_capable(::Type{Acb}) = true                         # Arblib.dft!
-_iitm_fft_capable(::Type) = false                             # direct fallback
+# ComplexF64 → FFTW (hard dep); Acb → Arblib.dft! (hard dep). Other Complex{<:AbstractFloat}
+# (Double64/BigFloat) become capable only when GenericFFT is loaded — see
+# `ext/TransitionMatricesGenericFFTExt.jl`; without it they use the direct azimuthal sum.
+_iitm_fft_capable(::Type{ComplexF64}) = true  # FFTW
+_iitm_fft_capable(::Type{Acb}) = true         # Arblib.dft!
+_iitm_fft_capable(::Type) = false             # direct fallback (incl. generic float w/o GenericFFT)
 
 struct _AzimuthalFourierWorkspace{CT, P}
     contrast::Matrix{CT}
@@ -49,11 +52,8 @@ function _azimuthal_fft_plan(::Type{ComplexF64}, Nφ, Nϑ)
     end
 end
 
-# GenericFFT path: fresh plan each time (bakes in precision-specific twiddles).
-# Do NOT cache — a BigFloat plan would go stale if precision changes.
-function _azimuthal_fft_plan(::Type{CT}, Nφ, Nϑ) where {CT <: Complex{<:AbstractFloat}}
-    return plan_fft(zeros(CT, Nφ, Nϑ), 1)   # GenericFFT; no flags kwarg
-end
+# The GenericFFT plan for non-Float64 complex-float types lives in
+# `ext/TransitionMatricesGenericFFTExt.jl` (GenericFFT is a weak dependency).
 
 # ---------- Workspace constructors ----------
 
@@ -102,13 +102,8 @@ function _apply_forward_dft!(spectrum::Matrix{ComplexF64},
     mul!(spectrum, plan, contrast)
 end
 
-# ---------- Forward column DFT: generic Complex{<:AbstractFloat} (GenericFFT) ----------
-# GenericFFT's DummyFFTPlan only supports `plan * A` (allocating), not mul!.
-function _apply_forward_dft!(spectrum::Matrix{CT},
-        contrast::Matrix{CT},
-        plan) where {CT <: Complex{<:AbstractFloat}}
-    spectrum .= plan * contrast
-end
+# The generic (non-ComplexF64) forward column DFT — used by the GenericFFT plan — lives in
+# `ext/TransitionMatricesGenericFFTExt.jl` (GenericFFT is a weak dependency).
 
 # ---------- Fourier coefficient extraction: generic Complex{<:AbstractFloat} ----------
 

src/IITM/nfold.jl

@@ -216,7 +216,9 @@ function transition_matrix_iitm(s::AbstractNFoldShape{N, T, CT}, λ, nₘₐₓ,
 end
 
 @testitem "Generic FFT: Complex{Double64} Prism matches ComplexF64 reference" begin
-    using TransitionMatrices: Prism, calc_T_iitm, calc_Csca, calc_Cext, Double64
+    using TransitionMatrices: Prism, calc_T_iitm, calc_Csca, calc_Cext
+    using DoubleFloats: Double64
+    using GenericFFT   # weak dep: load it so the Complex{Double64} azimuthal FFT path is exercised
 
     # Reference: ComplexF64 (FFTW path)
     m_f64 = complex(1.5)

src/TransitionMatrices.jl

@@ -2,15 +2,12 @@ module TransitionMatrices
 
 using Arblib
 using Arblib: ArbLike, AcbLike, ArbVectorLike, AcbVectorLike, ArbMatrixLike, AcbMatrixLike
-using DoubleFloats: Double64
 using FastGaussQuadrature: FastGaussQuadrature
 import FFTW
-using GenericFFT
 using ForwardDiff: ForwardDiff
 using GenericLinearAlgebra: Diagonal, GenericLinearAlgebra, cond, inv
 using LinearAlgebra: lu, mul!
 using OffsetArrays: OffsetArray
-using Quadmath: Quadmath, Float128, ComplexF128
 using Rotations: Angle2d, Rotation, RotMatrix2, RotZYZ
 using StableTasks: StableTasks
 using StaticArrays: SVector, SMatrix, SArray, @SVector, @SMatrix, @SArray
@@ -78,7 +75,8 @@ export AbstractShape, AbstractAxisymmetricShape, AbstractNFoldShape, volume,
        rmin, rmax, Spheroid, Cylinder, Chebyshev, Prism,
        SuperSpheroid, SuperEllipsoid, SuperSpheroidRevolved
 
-# Re-exports
-export RotZYZ, Double64, Float128, ComplexF128, Arb, Acb
+# Re-exports. Double64 / Float128 / ComplexF128 are NOT re-exported — DoubleFloats and
+# Quadmath are weak dependencies; `using DoubleFloats` / `using Quadmath` to access them.
+export RotZYZ, Arb, Acb
 
 end