API reference

Classical(; exact = false)

The q → 1 limit: ordinary SU(2) recoupling. With exact = true the result is a rational/radical value rather than floating point. Prefer Exact() for this exact classical request; Classical(exact=true) remains supported.

Level(k; T = Float64)

Floating-point values at the root of unity q = exp(iπ/(k+2)). k may be a range, which sweeps levels.

q6j(Level(10), 1, 1, 1, 1, 1, 1)
q6j(Level(2:50), 1, 1, 1, 1, 1, 1)          # one value per level
q6j(Level(10; T = BigFloat), 1, 1, 1, 1, 1, 1)

Exact()
Exact(k; form=:x)

Exact() requests exact classical evaluation at q = 1. It is shorthand for Classical(exact=true), reusing the direct classical evaluator and its rational/radical result types without constructing a level field. Use Symbolic() to retain q.

Exact(k) requests exact level values:

• form=:x (the default) returns an ExactX: the value in the real basis x = q + q⁻¹ = 2cos(π/(k+2)), at half the degree of the cyclotomic field, printed as the polynomial in x it is stored as. radical rewrites it in nested square roots where they exist. It is closed under *, inv and ^, and sums of it are ExactXSum.
• form=:canonical returns canonical cyclotomic coefficients in ℚ(ζ₂ₕ). This was the default before the real basis existed, and it is deprecated: it will be removed later.

Both use deterministic arithmetic; modular screening alone never establishes an exact zero. Exact braiding phases retain their QPhase form in either.

q6j(Exact(), 1, 1, 1, 1, 1, 1)                     # 1//6
q6j(Exact(3), 1, 1, 1, 1, 1, 1)                    # x − 2 at x = 2cos(π/5)
radical(q6j(Exact(3), 1, 1, 1, 1, 1, 1))           # (√5 − 3)/2, the Fibonacci level

At(q)

Numerical evaluation at a given real or complex q, without assuming a root of unity. Use Symbolic() instead for a generic expression.

Symbolic()

Request a parameter-independent representation instead of evaluation. Recoupling symbols return a SymbolicValue retaining the factorial rule without carrying out the sum. Display shows the deferred factorial rule in x for every label size; it never carries out the sum. Use x_form(value) to request full expansion, or phi_form(value) for cyclotomic factorization. qdim and theta_value use the same x-form interface. rmatrix and twist retain QPhase objects, displayed as algebraic phases over x with their q branch. No evaluation keywords are accepted. Use qeval(Symbolic(), rule) to lower a FactorialSum the same way. The weighted coefficients q3j and qcg do not yet support this target.

q6j(j1, j2, j3, j4, j5, j6; k=nothing, q=nothing, exact=false, T=Float64, workspace=nothing)

Evaluate the q6j symbol; the default is its classical value. Use k for a root-of-unity level or q for an analytic parameter (mutually exclusive). exact=true requests an exact classical or level value. q6j(Symbolic(), ...) retains a factorial rule with bounded x-form display; x_form expands it explicitly. Numerical classical/level calls evaluate the rule directly. eager=true is deprecated and uses the same evaluator.

q3j(j1, j2, j3, m1, m2, m3 = -m1-m2; k=nothing, q=nothing, exact=false, T=Float64, workspace=nothing)

The quantum 3j symbol of U_q(sl₂),

(j₁ j₂ j₃; m₁ m₂ m₃)_q = (−1)^{j₁−j₂−m₃} [2j₃+1]^{−1/2} ⟨j₁m₁; j₂m₂|j₃, −m₃⟩_q,

with the Clebsch–Gordan coefficient of qcg. The default is its classical value, the Wigner 3j symbol, which it reproduces bit for bit. Use q for a real or complex parameter and k for the level q = e^{iπ/(k+2)} (mutually exclusive). At q ≠ 1 the symbol is complex in general and is related to its relabellings by q-powers: (j₂ j₁ j₃; m₂ m₁ m₃)q = (−1)^{j₁+j₂+j₃} (j₁ j₂ j₃; m₁ m₂ m₃){1/q}, the same for m → −m, and (j₂ j₃ j₁; m₂ m₃ m₁)q = q^{m₂} (j₁ j₂ j₃; m₁ m₂ m₃)q.

From v0.5, the deformation-dependent weights are included; classical values are unchanged. Exact() is supported; Exact(k) and Symbolic() are not available yet. eager=true remains a deprecated alias for the standard evaluator.

qcg(j1, m1, j2, m2, j, m = m1 + m2; k=nothing, q=nothing, exact=false, T=Float64, workspace=nothing)
qcg((j1, m1), (j2, m2), (j, m); kwargs...)

The quantum Clebsch-Gordan coefficient ⟨j₁m₁; j₂m₂|j m⟩q of Uq(sl₂), with K|m⟩ = q^m|m⟩ and coproduct Δ(E) = E⊗K + K⁻¹⊗E; the default is its classical value, the ordinary Clebsch-Gordan coefficient. The arguments are the (jᵢ, mᵢ) pairs in the order of ⟨j₁m₁; j₂m₂|j m⟩, with m = m₁ + m₂ by default, or as three (j, m) tuples.

Use q for a real or complex parameter and k for the level, q = e^{iπ/(k+2)} (mutually exclusive). At real q > 0 the coefficients form an orthogonal matrix in (j, m₁) at fixed m; at complex q, including a level, orthogonality is bilinear, Cᵀ C = 1 (at a level, for complete sectors). Negative real q is evaluated as complex(q). At a level, labels outside the fusion rule j₁ + j₂ + j ≤ k give zero. exact = true and Exact() give the exact classical value; exact level values and Symbolic() are not available yet. qcg_matrix and qcg_row build whole matrices and rows by recurrence.

Square roots follow the balanced-factor convention: [n] = Π_{d|n,d>1} Ψ_d(q), Ψ_d(q) = q^{−φ(d)} Φ_d(q²), each Ψ_d rooted separately; Δ(F) = F⊗K + K⁻¹⊗F, with F the transpose of E.

The value is the 3j rule with term z weighted by q^{-z(j₁+j₂+j+1)} and the whole by q^{½(j₁+j₂-j)(j₁+j₂+j+1) + j₁m₂ - j₂m₁}, evaluated by the same scaled, escalating kernel. See also q3j.

qcg_matrix(j1, j2, m; k=nothing, q=nothing, T=nothing) -> (C, m1, j)

One sector of the coupling matrix: C[a, b] = qcg(j1, m1[a], j2, m − m1[a], j[b], m), with m1 running down from its largest value and j up over the spins of j1 ⊗ j2 that reach m (inside the fusion rule at a level). At positive real q and classically C is square and orthogonal, transpose(C) * C ≈ I. At a level, the retained columns satisfy the same bilinear identity, but fusion truncation can make C rectangular; then C * transpose(C) need not be the identity.

Built column by column from the Casimir recurrence in double words, O(1) work per entry, at real q > 0, classically and at a level; every entry whose own error estimate misses the package's accuracy promise is recomputed by qcg. Complex or negative q and BigFloat entries are evaluated entry by entry. T sets a precision floor, as in qcg: level and negative-q values retain their complex type, and a higher-precision q is never narrowed. For example, k = 10, T = BigFloat gives complex BigFloat entries.

qcg_matrix(j1, j2; k=nothing, q=nothing, T=nothing) -> (C, product, coupled)

The coupling matrix of j1 ⊗ j2: C[a, b] = qcg(j1, m1, j2, m2, j, m) with product[a] = (m1, m2) and coupled[b] = (j, m). Rows follow kron order — m1 from j1 down to −j1 and, for each, m2 from j2 down — so a column is the coupled vector |j m⟩ in the basis kron(e_{m1}, e_{m2}); columns run over j upwards and, within each j, m from j down. At positive real q and classically C is orthogonal; at a level only the j inside the fusion rule are columns, and transpose(C) * C ≈ I is a bilinear identity, without conjugation.

The matrix is dense, with (2j1+1)(2j2+1) rows and as many columns before fusion truncation (about 110 MB for Float64 at j1 = j2 = 30). qcg_matrix(j1, j2, m) holds one sector.

using LinearAlgebra
C, prod, coup = qcg_matrix(1, 1//2; q = 0.8)
transpose(C) * C ≈ I
b = findfirst(==((3//2, 3//2)), coup)
C[1, b] ≈ qcg(1, 1, 1//2, 1//2, 3//2, 3//2; q = 0.8)

qcg_row(j1, m1, j2, m2; k=nothing, q=nothing, T=nothing) -> (c, j)

The product state |j1 m1⟩⊗|j2 m2⟩ in the coupled basis: c[i] = qcg(j1, m1, j2, m2, j[i]) for every j of j1 ⊗ j2 that reaches m = m1 + m2 (inside the fusion rule at a level), j ascending. It is one row of qcg_matrix, restricted to its allowed coupling channels (some coefficients can vanish). At positive real q and classically sum(abs2, c) ≈ 1.

At real q > 0, classically and at a level the row comes from the three-term recurrence in j that q^{2m1} satisfies in the coupled basis (the q-analogue of J1z, from the q-Hahn polynomials), in double words: O(1) work per entry. Entries whose error estimate misses the package's promise are recomputed by qcg; complex or negative q and BigFloat are evaluated entry by entry. T is a precision floor; k = 10, T = BigFloat, for example, returns complex BigFloat entries.

fsymbol(j1, j2, j3, j4, j5, j6; k=nothing, q=nothing, exact=false, T=Float64, workspace=nothing)

Evaluate (−1)^(j1+j2+j4+j5) √([2j3+1][2j6+1]) {6j}, classically by default. Use k for a root-of-unity level or q for an analytic parameter (mutually exclusive). exact=true requests an exact classical or level value. fsymbol(Symbolic(), ...) retains a factorial rule with bounded x-form display; x_form expands it explicitly. Numerical classical/level calls evaluate the rule directly.

gsymbol(j1, j2, j3, j4, j5, j6; k=nothing, q=nothing, exact=false, T=Float64, workspace=nothing)

Evaluate √(Πᵢ[2ji+1]) {6j}, classically by default. Use k for a root-of-unity level or q for an analytic parameter (mutually exclusive). exact=true requests an exact classical or level value. gsymbol(Symbolic(), ...) retains a factorial rule with bounded x-form display; x_form expands it explicitly. Numerical classical/level calls evaluate the rule directly.

qdim(j; k=nothing, q=nothing, exact=false, T=Float64)

Quantum dimension [2j+1], classical by default. Use qdim(Symbolic(), j) for a rule-backed x-form.

rmatrix(j1::Spin, j2::Spin, j3::Spin; k=nothing, q=nothing, exact::Bool=false, T::Type=ComplexF64)

Returns the R-matrix phase. Formula: R = (-1)^{j1 + j2 - j3} q^{j3(j3+1) - j1(j1+1) - j2(j_2+1)}. The default is the classical value. rmatrix(Symbolic(), ...) returns a QPhase. At a level, exact=true retains the exact QPhase representation. Numerical level phases are complex even when T or Level(k; T=...) names a real type.

qint(n, p = 1; k, q, exact, T)

The q-integer [n]^p, classical by default — [n] = (qⁿ − q⁻ⁿ)/(q − q⁻¹), with the Turaev–Viro convention [0] = 1. n < 0 is zero.

The keywords are the ones every symbol takes: k for a level, q for a parameter, exact = true for the exact classical rational or, with k, the real-basis level value. qint(Symbolic(), n) returns the rule-backed symbolic x-form instead; use x_form to expand it explicitly.

qint(5)                 # 5.0
qint(5; k = 10)         # [5] at q = exp(iπ/12)
qint(5; q = 0.7)
qint(Exact(10), 5)      # exact, as a polynomial in x; `radical` for the surd

qfact(n, p = 1; k, q, exact, T)

The q-factorial ([n]!)^p = ([1][2]⋯[n])^p, classical by default, with the same keywords as qint. n < 0 is zero and [0]! = [1]! = 1.

qbinomial(n, m; k, q, exact, T)

The Gaussian binomial [n]! / ([m]! [n−m]!), classical by default, with the same keywords as qint. Outside 0 ≤ m ≤ n it is zero.

fmatrix(a, b, c, d; k=nothing, q=nothing, T=nothing) -> (F, e, f)

The F-matrix (associativity matrix) for four spins, relating the two fusion bases

((a b)ₑ c)_d   ⟷   (a (b c)_f)_d ,

with the intermediate spins labelling its rows (e) and columns (f). Entries are

F[i, j] = fsymbol(a, b, e[i], c, d, f[j]) = (-1)^{a+b+c+d} √([2e+1][2f+1]) {a b e; c d f} ,

so transpose(F) * F ≈ I. At complex q the matrix is complex orthogonal, not unitary.

Use k for a level, q for a real or complex parameter, or neither for the classical limit (k and q are exclusive). At a level the ranges hold only admissible intermediates, so F can be smaller or empty. Negative real q uses the branch of complex(q). T sets a precision floor: machine precision uses the column recurrence (O(1) per entry); BigFloat requests and inputs use scalar fsymbol calls instead.

using LinearAlgebra
F, e, f = fmatrix(1, 1, 1, 1; k = 6)
F' * F ≈ I                      # orthogonal
F[1, 1] ≈ fsymbol(1, 1, e[1], 1, 1, f[1]; k=6)

G, e, f = fmatrix(1, 1, 1, 1; q = 0.8 + 0.3im)
transpose(G) * G ≈ I            # complex orthogonal off the real axis, not unitary

fmatrix_labels(a, b, c, d; k=nothing) -> (e, f)

The row and column intermediates of fmatrix without computing it — the admissible e ∈ a⊗b with e⊗c ∋ d, and f ∈ b⊗c with a⊗f ∋ d.

bmatrix(a, b, c, d; k=nothing, q=nothing, T=ComplexF64, inverse=false) -> (B, e, e′)

The braid matrix: the change of basis that carries ((a b)_e c)_d to ((a c)_{e′} b)_d, i.e. braids b past c inside the fusion tree.

[B^{abc}_d]_{e e′} = Σ_f [F^{abc}_d]_{ef} · R^{bc}_f · [F^{acb}_d]_{e′f}

— conjugate by the F-matrix into the basis where the braiding is diagonal, apply it there, and come back. The two F-matrices share their column labels (both index f ∈ b⊗c with a⊗f ∋ d), which is what makes the middle factor a diagonal of R-phases and the whole thing one O(n³) pass with no Racah sums. Rows carry incoming channels e; columns carry outgoing channels e′. For a column vector of incoming amplitudes, the outgoing amplitudes are transpose(B) * state.

inverse = true gives the same construction with R⁻¹, so that

bmatrix(a, c, b, d; k, inverse = true) * bmatrix(a, b, c, d; k) == I

at any q. At a level B is also unitary, because the R-eigenvalues are phases there. At generic positive real q, those eigenvalues are real and need not have unit modulus, so B is generally not unitary. Measured over 107 matrices at k = 2…14: unitarity and the inverse relation both to 4.4e-16.

B, e, e′ = bmatrix(1, 1, 1, 1; k = 6)
B * B' ≈ I

smatrix(k; T = Float64) -> (S, labels)

The modular S-matrix of SU(2)k, `S{ab} = √(2/h)·sin((2a+1)(2b+1)π/h)withh = k+2`, together with the spins labelling it.

S is real, symmetric and orthogonal, S² = C is the charge conjugation matrix — the identity here, since every SU(2)k object is self-dual — and `S{0b}/S_{00}is the quantum dimension ofb`.

Built in O(n²) table lookups with no trigonometry per entry: along a row the argument advances by a constant, folded back into range by one comparison.

S, js = smatrix(6)
S * S' ≈ I
S[1, :] ./ S[1, 1] ≈ qdim.(js; k = 6)

tmatrix(k; T = ComplexF64, anomaly = true) -> (Tm, labels)

The modular T-matrix, T_{ab} = δ_{ab}·exp(-2πi c/24)·θ_a. With anomaly = false the framing phase is dropped and the diagonal is the bare twists, which is what most fusion-category conventions print.

Returned as a Diagonal-shaped dense matrix so that (S*T)^3 ≈ S^2 can be written directly.

twist(j; k = nothing, q = nothing, exact = false, T = ComplexF64)

The topological spin θ_j = q^{J(J+2)/2} = exp(2πi j(j+1)/(k+2)), J = 2j.

This is the convention the package's rmatrix already uses: rmatrix(j, j, 0) = (−1)^{2j} θ_j⁻¹, the self-braiding of a self-dual object through the vacuum. With exact = true the answer is the QPhase q^{J(J+2)/2}, exact at any level; the classical limit is 1.

q evaluates the same monomial away from a root of unity, which is what the rest of the package means by a generic parameter. It is the one place the modular layer has a generic-q value at all — smatrix, tmatrix and the rest are statements about a modular category and need the level.

QPhase

Represents the formal monomial sign * q^(q_pow), with sign in (-1, 0, 1). It retains the chosen q branch, not a level. Rational powers act on the formal exponent and the real sign; they need not equal a principal power after evaluation.

The affine factorial [a z + b]!^c.

FactorialSum

sign0 · (√)Π pre · Σ_{z=zlo}^{zhi} (−1)^{z·alternating} Π fac(z), where pre holds pairs n => c meaning [n]!^c and sits under a square root when sqrt_pre is set.

pre and fac are tuples for the fixed-shape symbols (6j, 3j, F, G), so a rule is a plain bits value: building one allocates nothing and the kernels are compiled for its exact shape. Any iterable collection works for other shapes.

qseries(factors, range; prefactor=(), sqrt_prefactor=false, alternating=false, sign=1)

Build a compact FactorialSum from affine factorial triples (a,b,c) representing [a*z+b]!^c. Construction retains the rule without evaluating the sum. Example: qeval(qseries([(2,0,1)], 1:3); q=1) is 2!+4!+6! = 746.

qeval(rule::FactorialSum; k=nothing, q=nothing, exact=false, T=Float64, workspace=nothing)

Evaluate a finite factorial rule, classically when no target is supplied (k and q are exclusive). Classical and level values use ratio kernels, exact classical values integer arithmetic, and real or complex q scaled ratios with adaptive precision; workspace reuses fixed-q tables. Factorial arguments must be nonnegative, and individually polar terms are rejected. qeval(Symbolic(), rule) returns a SymbolicValue retaining the rule; use x_form to expand it explicitly.

Evaluate the retained factorial rule with the same targets and keywords as qeval(::FactorialSum).

EvaluationWorkspace()

Reusable Float64 scratch and fixed-q analytic tables. Pass workspace=work to q6j, q3j, qcg, fsymbol, gsymbol, or qeval(::FactorialSum). Reuse it sequentially, never across concurrent tasks (batches give each worker its own). Tables are keyed by q, type and precision, with at most four tiers; CG scratch keeps one sector.

all_6j(; k, jmax = k, canonical = false, doubled = false)

Every 6j label set admissible at level k, as spin tuples with each spin at most jmax. canonical = true keeps one representative per class of the 144 symmetries (tetrahedral and Regge); F- and G-symbols are not invariant, so sweep them over all labels. doubled = true returns DoubledLabels, which skips the conversion in batches.

labels = all_6j(k = 6)
vals = q6j(labels; k = 6)        # tables built once, work threaded

empty_caches!()

Clear the cached level-k tables (numeric, root-of-unity, q-integer, modular and exact cyclotomic). Useful for freeing numerical-table memory in long sessions or before benchmarking. Reusable symbolic polynomials and number fields are retained, as is the small classical prime-power sieve. The per-level numerical caches also empty themselves when one passes QRecoupling.LEVEL_CACHE_LIMIT[] bytes (512 MiB by default), so a sweep over many levels does not exhaust memory.

RadExpr

A real number written as rat + Σᵢ cᵢ √(eᵢ), with each eᵢ again a RadExpr. This is what the Lagrange descent produces and what radical returns on success; float evaluates it, so a rendered formula can be checked against the value it claims to be.

phi_form(value; maxdeg = 400, basis = :q) -> PhiForm

Expand a SymbolicValue in q as Φ factors and a remaining polynomial. Obtain the input with Symbolic(), for example phi_form(q6j(Symbolic(), 1,1,1,1,1,1)). The expensive expansion: it factors the numerator over ℤ[q], about 38× the cost of x_form on {6 6 6; 6 6 6}. Above maxdeg the numerator is reported by size, not factored; that is a different answer, recorded in truncated, not a cheaper one. basis = :x writes the remainder in x = q + q⁻¹ (the Φ factors stay in q).

prove_identity(lhs, rhs; kmax = 4096) -> NamedTuple

Decide lhs == rhs for XValue or XSum arguments, in ℤ[x]. The verdict is :proved or :refuted, and a :proved verdict is a statement about every q, hence about every level at once: a polynomial identity in x specialises to x = 2cos(π/h) for every h.

Returned fields:

• verdict — :proved or :refuted
• degree — the degree of the compared numerator, i.e. the size of the certificate
• classes — the square classes involved
• exceptional — levels at which a denominator vanishes, where the identity is silent rather than false
• witness — for :refuted, the nonzero difference

No tolerance, prime or bound enters: normalised XValues compare structurally, and for XSum distinct classes are independent over ℚ(x), so an empty difference is a proof.

xpolynomial(v::ExactX) -> Vector{Rational{BigInt}}

Coefficients of P in 1, x, x², …, so that v = P(x)·√(R(x)) with x = 2cos(π/(k+2)).

radicand(v::ExactX) -> Vector{Rational{BigInt}}

Coefficients of R, the polynomial under the square root. Empty radical class gives R = 1.

iszero_at(k, j1, j2, j3, j4, j5, j6) -> Bool
iszero_at(symbol, k, labels...) -> Bool
iszero_at(k, labels::AbstractVector) -> BitVector

Decide whether a symbol vanishes at level k, including cancellation between finite terms. symbol selects q6j, fsymbol or gsymbol (default q6j). Given a collection of label tuples, the test runs over the batch with shared tables and threads.

Structural identities prove zeros directly. Other candidates pass through modular screening and exact confirmation. In a batch, screening can use threads; exact confirmations run serially and can cost more at high degree. Use level_spectrum for screening without exact confirmation.

iszero_at(20, 5, 5, 5, 5, 5, 5)          # true: a cancellation zero
iszero_at(fsymbol, 10, 1, 1, 1, 1, 1, 1)
count(iszero_at(6, all_6j(k = 6); threads=1))       # how many symbols vanish at level 6

issingular_at(k, j1, ..., j6) -> Bool
issingular_at(symbol, k, labels...) -> Bool

Whether these labels are singular at level k: the level-k rule has a contributing term of negative valuation, i.e. a q-factorial in a denominator that vanishes. Decided from the valuations alone, with no arithmetic. symbol selects q6j, fsymbol or gsymbol (default q6j).

Singular labels are exactly the ones a level cannot represent — for the 6j, those with a triangle sum above 2k. Labels that are admissible at the level are never singular (their prefactor and term valuations are both non-negative), so this query is false throughout the admissible set; it answers a question about the labels, not about the value the symbol functions return. Those follow the Turaev–Viro convention and give 0 outside the admissible set, and level_spectrum reports such levels as :inadmissible.

issingular_at(12, 10, 10, 10, 10, 10, 10)     # true: triangle sums 30 > 2k
issingular_at(30, 10, 10, 10, 10, 10, 10)     # false: admissible from k = 30
issingular_at(4, 1, 1, 5, 1, 1, 5)            # false: no triangle, so no representation to be singular

level_spectrum(j1, ..., j6; k = 2:100, cancellation = true, prove = false) -> Vector{Symbol}
level_spectrum(symbol, labels...; k, cancellation, prove) -> Vector{Symbol}

How one symbol behaves as the level varies, one entry per level in k: :inadmissible (the labels are not in the theory at that level, so the symbol is zero by convention), :pole, :zero (every term vanishes), :cancels (the finite terms cancel) or :finite. Cheap because the valuations are closed-form in the level; pass cancellation = false to skip the cancellation test and report :finite wherever terms contribute.

By default :cancels is a modular candidate, not a proof. With prove = true every candidate is settled: by the structural identities (pairwise or reflection cancellation) or, failing those, by exact arithmetic at the level, so :cancels then means a proved zero and an unconfirmed candidate is reported as :finite. Screening runs on threads; exact confirmations run serially afterwards, as in iszero_at, and can cost milliseconds each at high degree.

level_spectrum(5, 5, 5, 5, 5, 5; k = 2:60)                 # screen
level_spectrum(5, 5, 5, 5, 5, 5; k = 2:60, prove = true)   # every :cancels proved