Improved algorithm for U1 implementation of Complex phi4

src/models/phi4_complex.jl

@@ -253,58 +253,75 @@ function phi4_complex(::Type{Z2Irrep ⊠ Z2Irrep}, K::Integer, μ0::Float64, λ:
 end
 function phi4_complex(::Type{U1Irrep}, K::Integer, μ0::Float64, λ::Float64; T::Type{<:Number} = Float64)
     if K % 2 != 0
-        error("K must be even to split into even/odd groups")
+        error("K must be even")
     end
 
-    # precompute
     moments = precompute_moments_complex(K, μ0, λ)
-    # log factorials 0..K-1
     logfact = log.(factorial.(0:(K - 1)))
 
-    T_arr = zeros(T, K, K, K, K, K, K, K, K)
+    V1 = U1Space([U1Irrep(q) => 1 for q in 0:(K - 1)]...)
+    V2 = U1Space([U1Irrep(q) => 1 for q in 0:-1:(-K + 1)]...)
+    W = fuse(V1, V2)
+
+    # Build multiplicity index lookup:
+    # For fused charge q = a - B, the pairs (a, B) are ordered by increasing a.
+    # mult_index[q][a] = the 1-based index into the W-block for that pair.
+    mult_index = Dict{Int, Dict{Int, Int}}()
+    for q in (-(K - 1)):(K - 1)
+        pairs = Int[]
+        for a in max(0, q):(min(K - 1, q + K - 1))
+            B = a - q
+            (0 <= B <= K - 1) || continue
+            push!(pairs, a)
+        end
+        mult_index[q] = Dict(a => i for (i, a) in enumerate(pairs))
+    end
 
-    @threads for a in 0:(K - 1)
-        for b in 0:(K - 1), c in 0:(K - 1), d in 0:(K - 1), e in 0:(K - 1), f in 0:(K - 1), g in 0:(K - 1)
-            # solve delta for l4:
-            # b + d + e + g = a + c + f + h
-            h = e + g + b + d - a - c - f
+    T_fused = TensorMap(zeros(T, K^2, K^2, K^2, K^2), W ⊗ W ← W ⊗ W)
 
-            if h < 0 || h > K - 1
-                continue
-            end
+    for (split_tree, fuse_tree) in fusiontrees(T_fused)
+        q1 = Int(split_tree.uncoupled[1].charge)
+        q2 = Int(split_tree.uncoupled[2].charge)
+        q3 = Int(fuse_tree.uncoupled[1].charge)
+        q4 = Int(fuse_tree.uncoupled[2].charge)
 
-            # total power
-            sum_power = a + b + c + d + e + f + g + h
-            n = 1 + sum_power
-            # quick skip if moment is zero
-            M = moments[n + 1]
-            if M == 0.0
-                continue
-            end
+        block = T_fused[split_tree, fuse_tree]
+        # block has shape (mult(q1) × mult(q2)) × (mult(q3) × mult(q4))
+        # index as block[i1, i2, i3, i4] where each i is the multiplicity index
 
-            # denomenator via logfacts
-            logdenom = 0.5 * (
-                log(2) * sum_power +
-                    logfact[a + 1] + logfact[b + 1] + logfact[c + 1] + logfact[d + 1] + logfact[e + 1] + logfact[f + 1] + logfact[g + 1] + logfact[h + 1]
-            )
-            denom = exp(logdenom)
+        for a in max(0, q1):min(K - 1, q1 + K - 1)
+            B = a - q1;  (0 <= B <= K - 1) || continue
+            i1 = mult_index[q1][a]
 
-            val = 2π * M / denom
+            for c in max(0, q2):min(K - 1, q2 + K - 1)
+                D = c - q2;  (0 <= D <= K - 1) || continue
+                i2 = mult_index[q2][c]
 
-            # store into array (indices +1)
-            T_arr[a + 1, b + 1, c + 1, d + 1, e + 1, f + 1, g + 1, h + 1] = val
-        end
-    end
+                for e in max(0, q3):min(K - 1, q3 + K - 1)
+                    F = e - q3;  (0 <= F <= K - 1) || continue
+                    i3 = mult_index[q3][e]
 
-    # Build U1 spaces
-    V1 = U1Space([U1Irrep(q) => 1 for q in 0:(K - 1)]...)
-    V2 = U1Space([U1Irrep(q) => 1 for q in 0:-1:(-K + 1)]...)
-    T_unfused = TensorMap(T_arr, V1 ⊗ V2 ⊗ V1 ⊗ V2 ← V1 ⊗ V2 ⊗ V1 ⊗ V2)
+                    for g in max(0, q4):min(K - 1, q4 + K - 1)
+                        H = g - q4;  (0 <= H <= K - 1) || continue
+                        i4 = mult_index[q4][g]
 
-    U = isometry(fuse(V1, V2), V1 ⊗ V2)
-    Udg = adjoint(U)
+                        sum_power = a + B + c + D + e + F + g + H
+                        M = moments[sum_power + 2]
+                        (M == 0.0) && continue
+
+                        logdenom = 0.5 * (
+                            log(2) * sum_power +
+                                logfact[a + 1] + logfact[B + 1] + logfact[c + 1] + logfact[D + 1] +
+                                logfact[e + 1] + logfact[F + 1] + logfact[g + 1] + logfact[H + 1]
+                        )
+
+                        block[i1, i2, i3, i4] += 2π * M / exp(logdenom)
+                    end
+                end
+            end
+        end
+    end
 
-    @tensor T_fused[-1 -2; -3 -4] := T_unfused[1 2 3 4; 5 6 7 8] * U[-1; 1 2] * U[-2; 3 4] * Udg[5 6; -3] * Udg[7 8; -4]
     return T_fused
 end