diff --git a/Project.toml b/Project.toml
index 1e55d15..4c96297 100644
--- a/Project.toml
+++ b/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"
diff --git a/benchmark/Project.toml b/benchmark/Project.toml
index 69acf5c..c28f708 100644
--- a/benchmark/Project.toml
+++ b/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"
diff --git a/benchmark/benchmarks.jl b/benchmark/benchmarks.jl
index 7f4820a..4507197 100644
--- a/benchmark/benchmarks.jl
+++ b/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()
 
diff --git a/ext/TransitionMatricesDoubleFloatsExt.jl b/ext/TransitionMatricesDoubleFloatsExt.jl
new file mode 100644
index 0000000..7a8f439
--- /dev/null
+++ b/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
diff --git a/ext/TransitionMatricesGenericFFTExt.jl b/ext/TransitionMatricesGenericFFTExt.jl
new file mode 100644
index 0000000..97c734c
--- /dev/null
+++ b/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
diff --git a/ext/TransitionMatricesQuadmathExt.jl b/ext/TransitionMatricesQuadmathExt.jl
new file mode 100644
index 0000000..7a3bcad
--- /dev/null
+++ b/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
diff --git a/packages/EBCMPrecisionLossEstimators/Project.toml b/packages/EBCMPrecisionLossEstimators/Project.toml
index be74b49..dab21f0 100644
--- a/packages/EBCMPrecisionLossEstimators/Project.toml
+++ b/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"
diff --git a/src/IITM/fourier.jl b/src/IITM/fourier.jl
index a8667f3..62c8d11 100644
--- a/src/IITM/fourier.jl
+++ b/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} ----------
 
diff --git a/src/IITM/nfold.jl b/src/IITM/nfold.jl
index 07d9aff..4bfc2ce 100644
--- a/src/IITM/nfold.jl
+++ b/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)
diff --git a/src/TransitionMatrices.jl b/src/TransitionMatrices.jl
index 182ca4a..b12db00 100644
--- a/src/TransitionMatrices.jl
+++ b/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
diff --git a/src/compat/index.jl b/src/compat/index.jl
index 2ea7c36..dc8418b 100644
--- a/src/compat/index.jl
+++ b/src/compat/index.jl
@@ -1,34 +1,12 @@
-Base.convert(::Type{Float128}, x::ArbLike) = Float128(BigFloat(x))
-Quadmath.Float128(x::ArbLike) = Float128(BigFloat(x))
+# Float128 (Quadmath) ⇄ Arb conversions live in `ext/TransitionMatricesQuadmathExt.jl`
+# (Quadmath is a weak dependency).
 Base.round(x::Arb, ::RoundingMode{:Up}) = ceil(BigFloat(x))
 Base.abs2(x::AcbLike) = abs2(real(x)) + abs2(imag(x))
 Base.complex(::Type{Arb}) = Acb
 
-# 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. (A package upgrade recompiles this module, so the
-# guard is re-evaluated against the new DoubleFloats version.)
-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
+# Double64 (DoubleFloats) shims — the `cbrt` workaround, `precision(Complex{Double64})`,
+# and `Arblib.set!(::Double64)` — live in `ext/TransitionMatricesDoubleFloatsExt.jl`
+# (DoubleFloats is a weak dependency).
 
 function Base.inv(x::Matrix{Arb})
     a = ArbMatrix(x)
@@ -44,8 +22,6 @@ end
 
 Base.precision(::Type{ComplexF32}) = 24
 Base.precision(::Type{ComplexF64}) = 53
-Base.precision(::Type{Complex{Double64}}) = 106
-Base.precision(::Type{ComplexF128}) = 113
 Base.precision(::Type{Complex{Arb}}) = precision(Arb)
 
 function Base.precision(::Type{
@@ -58,13 +34,6 @@ function Arblib.set!(arb::Arblib.ArbLike, dual::ForwardDiff.Dual)
     Arblib.set!(arb, dual.value)
 end
 
-function Arblib.set!(arb::Arblib.ArbLike, val::Float128)
-    Arblib.set!(arb, BigFloat(val))
-end
-
-function Arblib.set!(arb::Arblib.ArbLike, val::Double64)
-    Arblib.set!(arb, BigFloat(val))
-end
 
 Base.Int64(x::ArbLike) = round(Int64, BigFloat(x))
 
diff --git a/src/shapes/superspheroid.jl b/src/shapes/superspheroid.jl
index badd02b..5fd1d0b 100644
--- a/src/shapes/superspheroid.jl
+++ b/src/shapes/superspheroid.jl
@@ -40,16 +40,8 @@ SuperSpheroid(a, c, n, m) = SuperSpheroid{typeof(a), typeof(m)}(a, c, n, m)
 # needless detour and, because `BigFloat(::Arb)` uses the *global* precision, silently
 # truncate high-precision `T` (Arb, large BigFloat). `precision(x)` reads the instance's
 # precision, so Float64/Double64/Float128/BigFloat/Arb shapes each compute at their own.
-# Working precision (bits) from the value's own precision, so high-precision shapes
-# (Double64, BigFloat, Arb instances) are not truncated. `precision(::Float128)` (the
-# instance method) is unimplemented in the current Quadmath/Julia (it routes to a
-# missing `_precision_with_base_2`), so Float128 is pinned to its 113-bit mantissa via
-# the working type method `precision(Float128)`.
-_beta_prec(::Float128) = precision(Float128)
-_beta_prec(x::Real) = precision(x)
-
 function _beta_lgamma(x::T, y::T) where {T <: Real}
-    prec = max(_beta_prec(x), _beta_prec(y))
+    prec = max(precision(x), precision(y))
     lΓx = Arblib.lgamma!(Arb(; prec), Arb(x; prec); prec)
     lΓy = Arblib.lgamma!(Arb(; prec), Arb(y; prec); prec)
     lΓxy = Arblib.lgamma!(Arb(; prec), Arb(x + y; prec); prec)
@@ -79,7 +71,9 @@ has_symmetric_plane(::SuperSpheroid) = true
 
 @testitem "SuperSpheroid utility functions" begin
     using TransitionMatrices: SuperSpheroid, volume, volume_equivalent_radius,
-                              has_symmetric_plane, Double64, Float128, ComplexF128
+                              has_symmetric_plane
+    using DoubleFloats: Double64
+    using Quadmath: Float128, ComplexF128
 
     # n=1 reduces to spheroid: V = 4π/3 * a^2 * c
     @testset "n=1 reduction to spheroid" begin
diff --git a/test/Project.toml b/test/Project.toml
index dd6f894..acd4140 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,8 +1,11 @@
 [deps]
 Arblib = "fb37089c-8514-4489-9461-98f9c8763369"
+DoubleFloats = "497a8b3b-efae-58df-a0af-a86822472b78"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GSL = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
+GenericFFT = "a8297547-1b15-4a5a-a998-a2ac5f1cef28"
+Quadmath = "be4d8f0f-7fa4-5f49-b795-2f01399ab2dd"
 Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
diff --git a/test/compat.jl b/test/compat.jl
index 58f5d04..6c3c138 100644
--- a/test/compat.jl
+++ b/test/compat.jl
@@ -2,6 +2,8 @@
     using TransitionMatrices
     using Arblib
     using ForwardDiff
+    using Quadmath: Quadmath, Float128, ComplexF128   # weak deps: load explicitly
+    using DoubleFloats: Double64
 
     @test precision(ComplexF32) == 24
     @test precision(ComplexF64) == 53
@@ -15,7 +17,7 @@
 
     arb = Arb(1.25)
     @test convert(Float128, arb) == Float128(1.25)
-    @test TransitionMatrices.Quadmath.Float128(arb) == Float128(1.25)
+    @test Quadmath.Float128(arb) == Float128(1.25)
     @test Float32(arb) == Float32(1.25)
     @test Int64(Arb(2.0)) == 2
     @test round(Arb(1.2), RoundUp) == ceil(BigFloat(Arb(1.2)))
@@ -63,6 +65,7 @@ end
 
 @testitem "cbrt on Double64 returns a Double64 (DoubleFloats v1.9 workaround)" begin
     using TransitionMatrices
+    using DoubleFloats: Double64
 
     # Upstream DoubleFloats returns a (hi, lo) tuple here; ensure our shim wraps it.
     @test cbrt(Double64(8.0)) isa Double64