[dgrinko]

I got severe numerical instabilities with large number of systems, so I asked Claude to fix it.

Numerically stable computation of lower-weight CGCs

lower_weight_CGC! obtains the CGCs of every non-highest weight from those of its parent
weights by solving the lowering equations

J⁻₃ |m₃⟩ = (J⁻₁ ⊗ 𝕀 + 𝕀 ⊗ J⁻₂) |m₁, m₂⟩

weight by weight. This recursion is numerically unstable: rounding errors along the other irreps
s3′ ⊂ s1 ⊗ s2 are amplified at every step. For SU(2), the factor per step is
c_{J′}(M) / c_J(M) > 1 for J′ > J, which compounds to binomial(J′ + J, 2J) over a full
weight chain. For CGC(SUNIrrep{2}(2j, 0), SUNIrrep{2}(2j, 0), s3), the largest error over all s3
compared with the exact Racah formula is:

2j	10	20	30	40	50	60	70	80
max error, before	1.5e-13	4.7e-9	2.0e-5	0.23	7.3e3	4.3e7	1.0e12	1.2e16
max error, after	8.3e-16	1.1e-15	1.4e-15	2.4e-15	3.1e-15	9.2e-15	9.3e-15	1.0e-14

After the fix the error stays at 1e-14 up to 2j = 100; that floor comes from purge! dropping
entries below TOL_PURGE = 1e-14.

Larger SU(3) and SU(4) irreps are affected too. The table gives the largest violation of the defining
properties of the CGCs (isometry, J± C = C J±₃ for all simple roots):

product	before	after
SU(3) (6,3,0) ⊗ (6,3,0) → all 25	5.1e-12	3.9e-14
SU(3) (14,0,0) ⊗ (14,14,0) → (14,7,0)	5.5e-6	8.2e-15
SU(3) (16,0,0) ⊗ (16,16,0) → (16,8,0)	2.7e-4	1.0e-14
SU(4) (10,0,0,0) ⊗ (10,10,10,0) → (10,5,5,0)	1.8e-6	3.3e-14

Fix

At every weight, the raising equations (J⁺₁ ⊗ 𝕀 + 𝕀 ⊗ J⁺₂) X = Y J⁺₃ (with J⁺ = J⁻ᵀ) are now
imposed as well, and both sets are solved in the least-squares sense. Let Lₗ and Kₗ be the
blocks of J⁻ₗ on s1 ⊗ s2 and on s3 from the parent weight w3 + αₗ to w3, and let Yₗ be
the CGC block of that parent. Then

min_X ∑ₗ ‖X Kₗ - Lₗ Yₗ‖² + ‖Lₗᵀ X - Yₗ Kₗᵀ‖²

has the normal equations X A + B X = R, where:

• A = ∑ₗ Kₗ Kₗᵀ is small and positive definite;
• B = ∑ₗ Lₗ Lₗᵀ is positive semidefinite. The sum runs over all l, including those where
  w3 + αₗ is not a weight of s3; for those, the raising equations read Lₗᵀ X = 0;
• R = 2 ∑ₗ Lₗ Yₗ Kₗᵀ.

The exact CGCs satisfy both sets of equations, and the solution is unique, so nothing changes in
exact arithmetic, including the gauge. On all 1323 CGCs between small irreps of SU(2) to SU(5), the
old and new code agree to 1e-12, which is the old code's own error.

For SU(2), each step now multiplies the error along J′ by 2ab / (a² + b²) ≤ 1 (with
a = c_{J′}, b = c_J) instead of a / b. For SU(N ≥ 3) this is not an exact statement per irrep,
because ∑ₗ J⁻ₗ J⁺ₗ over the simple roots is not diagonal in the irrep decomposition. The tables
above show the same damping in practice.

The system is solved by diagonalizing A and doing one dense Cholesky factorization of
B + μₖ 𝕀 per eigenvalue μₖ. B is accumulated from the sparse columns of the raising operators.
In lower_weight_CGC! itself, only the qr! solve is replaced (by lowering_solve). The other
additions are the helpers lowering_solve, raising_columns and raising_gram.

Cost

Computing all CGCs between irreps of dimension ≤ 31 (SU(2)), ≤ 64 (SU(3), SU(4)) and ≤ 50 (SU(5))
takes 1.3–2.0× as long as before. Single larger products, like those in the tables above, take
2–2.7× as long. The dense Cholesky factorization of the weight-space block is the main extra cost.
CGCs are cached in RAM and on disk, so this cost is paid once per CGC.

Disk cache

Cached CGCs written by the old code can be wrong, so the cache file names now carry a version
suffix (CGC_CACHE_VERSION = "v2"), and old entries are never read back. As a result, existing disk
caches get recomputed. The old files stay on disk until clear_disk_cache!() is called.

Tests

• test/clebschgordan.jl existed but was empty and was not included in runtests.jl. It now
  contains:

  • a comparison of SU(2) CGCs up to spin 40 with the exact ones of SU2Irrep (from
    WignerSymbols);
  • isometry and intertwiner (J⁻ C = C J⁻₃) checks for SU(3) and SU(4) products, including
    outer multiplicities 3 and 4.

  On main, 12 of these 19 tests fail.

• The full test suite passes locally (fusion trees, pentagon/hexagon, caching, Aqua), and Runic
  reports no formatting changes.