[lkdvos]

Braid and transpose a FusionTreeBlock of a ProductSector through its factors: the fusion trees of a product are pairs of fusion trees of two parts of the product, so the recoupling matrix of a block is the Kronecker product of the recoupling matrices of the two parts.

• A product with both UniqueFusion and other factors is split into its other factors and its UniqueFusion factors. The latter have a single tree pair, such that a block is a block of the other factors with a fixed tree pair and a scalar. All blocks that only differ in their abelian charges therefore share the cached manipulations of the non-abelian factors: for fℤ₂ × U(1) × SU(2) × U(1) below, 26004 blocks map onto 142 SU(2) blocks.
• Any other product (e.g. several non-abelian factors) is split into its first and its remaining factors, recursively.
• The canonical order of a product block, i.e. of FusionTreeBlock{I}(uncoupled, isdual) and therefore of fusionblocks, is now the Kronecker order of the two parts, such that the factorized manipulations produce canonical blocks without reordering. This only affects the order of the trees within a block: the order of fusiontrees(::HomSpace), and thus the data layout of tensors, is unchanged.
• fusionblocks also benefits: the abelian part of a block is constructed first, which skips the non-abelian construction for the (many) combinations of uncoupled sectors with mismatching abelian charges.
• Blocks of UniqueFusion sectors, which hold a single tree pair, are braided and transposed through the tree pair.

The code lives in src/fusiontrees/producttrees.jl; braid and transpose dispatch on ProductSectorBlock.

Cold construction of the recoupling data for one permutation of a tensor from a PEPS CTMRG step (the reproducers of #580 and #581), main → this PR, minimum of 5:

	main	this PR
fℤ₂ × U(1) × SU(3), 1 thread	6.2–6.7 s	3.7 s (1.7×)
fℤ₂ × U(1) × SU(3), 8 threads	1.7–2.0 s	1.0 s (1.7×)
fℤ₂ × U(1) × SU(2) × U(1), 1 thread	4.6–5.0 s	0.85 s (5.5×)
fℤ₂ × U(1) × SU(2) × U(1), 8 threads	2.6–3.1 s	0.75 s (3.8×)

Ranges are over two interleaved rounds; Julia 1.13.2, shared 32-core machine. For SU(2), fusionblocks (~0.4 s, single-threaded) is now the largest part of the remaining time, which is what #581 addresses; for SU(3), the remaining time is mostly braiding genuinely distinct SU(3) blocks, which is what #580 addresses. Both should compose with this PR.

Tested: a new test compares the factorized braid/transpose with the unfactorized fsbraid/fstranspose for all product sectors in the test list. Additionally, for every block of a 3 ← 3 space of every product sector with a simple unit (including A4Irrep ⊠ A4Irrep with multiplicities, FermionParity ⊠ SU2Irrep ⊠ SU2Irrep and the anyonic products), and for the blocks of both reproducers, the results agree with the unfactorized ones to ~1e-15.

🤖 Generated with Claude Code

[leburgel]

This looks great, but I think there might be some additional gain for the specific case where there's only one symmetry factor which doesn't actually have UniqueFusion. This just happens to be the case I'm most interested in so maybe this would make things too artificial, but I think it could be worth it:

For products with a single factor without UniqueFusion (e.g. fℤ₂ × U(1) × SU(2) × U(1)), most of the time currently goes into the recursion around the braids rather than the braids themselves: splitting off the first factor gives three levels here, each copying all trees and building a block, and the result is then checked against the canonical block. All of this can be avoided by choosing the first factor without UniqueFusion as the head instead of the first factor, so that with a single such factor the tail is a product of UniqueFusion factors with a single tree per block. That needs three changes that build on each other:

1. Choose the head. Take the first factor without UniqueFusion. The factors before it have a single vertex label, so the head is still the slowest label of the product, as the splitting of the vertex labels in headtail assumes. With several such factors, the UniqueFusion factors then no longer get a level of their own either:

# blocks of products with `UniqueFusion` have a single tree and don't reach `product_transform` (see 3.)
@generated headfactor(::Type{ProductSector{T}}) where {T} =
    something(findfirst(F -> !(FusionStyle(F) isa UniqueFusion), collect(fieldtypes(T))), 1)
headsector(s::ProductSector) = s.sectors[headfactor(typeof(s))]
function tailsector(s::ProductSector)
    tail = _deletefactor(s.sectors, Val(headfactor(typeof(s))))
    return length(tail) == 1 ? only(tail) : ProductSector(tail)
end
@generated _deletefactor(t::Tuple, ::Val{k}) where {k} =
    Expr(:tuple, (:(t[$i]) for i in 1:fieldcount(t) if i != k)...)

2. Join head and tail at that position. The tree ⊠ puts the head first, so product_transform needs its inverse for the product type I instead (⊠ on trees isn't used elsewhere):

joinsectors(::Type{I}, h, t) where {I} =
    I(_insertfactor(t isa ProductSector ? t.sectors : (t,), h, Val(headfactor(I))))
@generated _insertfactor(t::Tuple, h, ::Val{k}) where {k} =
    Expr(:tuple, (i < k ? :(t[$i]) : i == k ? :h : :(t[$(i - 1)]) for i in 1:(fieldcount(t) + 1))...)

function jointrees(::Type{I}, h::FusionTree, t::FusionTree) where {I}
    vertices = ... # as in `⊠`, see 3.
    join(a, b) = joinsectors(I, a, b)
    return FusionTree(
        map(join, h.uncoupled, t.uncoupled), join(h.coupled, t.coupled), h.isdual,
        map(join, h.innerlines, t.innerlines), vertices
    )
end

with jointrees(I, h₁, t₁), jointrees(I, h₂, t₂) in place of h₁ ⊠ t₁, h₂ ⊠ t₂ in product_transform.

3. Use that the tail has unique fusion. With a single factor without UniqueFusion, the tail has UniqueFusion, which removes the remaining overhead in three places:

• A block of the tail has a single tree, so it doesn't need to recurse either: in braid(::ProductSectorBlock, ...) and transpose(::ProductSectorBlock, ...),

  (isempty(src) || FusionStyle(sectortype(src)) isa UniqueFusion) && return fsbraid((src, p, levels)) # resp. fstranspose

• The vertex labels of the tail are all 1, so the head keeps the labels of the product as they are, without vertexmultiplicities (and the same the other way around, in headtail and jointrees):

  # in `headtail`
  if FusionStyle(headsector(f.coupled)) isa UniqueFusion
      vertices_h = map(Returns(1), f.vertices)
  elseif FusionStyle(tailsector(f.coupled)) isa UniqueFusion
      vertices_t = map(Returns(1), f.vertices)
  else # as now
  ...
  # in `jointrees`
  if FusionStyle(I) isa GenericFusion && length(h) > 1
      vertices = FusionStyle(h) isa UniqueFusion ? t.vertices :
          FusionStyle(t) isa UniqueFusion ? h.vertices :
          map((μ, ν, n) -> ν + (μ - 1) * n, h.vertices, t.vertices, vertexmultiplicities(t))
  end

• With a single tree on one side, the product trees come out in the canonical order of the other side's block, which transform returns in canonical order, so the canonical block doesn't need to be rebuilt; in product_transform, before canonical = ...:

  (FusionStyle(first(first(src_h))) isa UniqueFusion || FusionStyle(first(first(src_t))) isa UniqueFusion) &&
      return dst => U

For products with several factors without UniqueFusion, the general head/tail recursion remains, now only over those factors. With these changes the braid transformers are identical to the current ones for fℤ₂ × U(1) × SU(3) and fℤ₂ × U(1) × SU(2) × U(1).

[leburgel]

One more thing for products with a single factor without UniqueFusion: the tail is now a product of UniqueFusion factors, whose blocks have a single tree, but those still go through product_transform factor by factor. They can go to fsbraid and fstranspose directly, which handle single trees already:

function braid(src::ProductSectorBlock, p::Index2Tuple, levels::Index2Tuple)
    # blocks of products without multiplicities have a single tree
    (isempty(src) || FusionStyle(sectortype(src)) isa UniqueFusion) && return fsbraid((src, p, levels))
    return product_transform(sectorscalartype(sectortype(src)), src) do x
        return braid(x, p, levels)
    end
end

and in the same way with fstranspose((src, p)) in transpose(src::ProductSectorBlock, p).