From 70864a3630d99844be4af6c56465870d3cb7dcad Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 14 Oct 2026 12:59:06 +0200 Subject: [PATCH 01/12] Construct `FusionTreeBlock`s from the fusion product of the shorter side For sectors with a simple unit, a superset of the coupled sectors suffices, as the fusion tree iterators yield nothing for the others: the unit if one side is empty, and otherwise the fusion product of the shorter side only, instead of the intersection of the fusion products of both sides. The trees and their order are unchanged. Co-Authored-By: Claude Opus 5.5 --- src/fusiontrees/fusiontrees.jl | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index f20445f36..9c21de89a 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/src/fusiontrees/fusiontrees.jl @@ -140,7 +140,15 @@ function FusionTreeBlock{I}( end # note: cs is unique due to ⊗ returning unique fusion results - cs = if N₁ == 0 + cs = if UnitStyle(I) isa SimpleUnit + # a superset of the coupled sectors suffices, as there are no fusion trees for the others: + # this avoids the fusion product of the longer side + if N₁ == 0 || N₂ == 0 + [unit(I)] + else + sort!(vec(collect(⊗(uncoupled[N₁ <= N₂ ? 1 : 2]...)))) # `vec`: products of `ProductSector`s have a shape + end + elseif N₁ == 0 sort!(collect(filter(isunit, ⊗(uncoupled[2]...)))) elseif N₂ == 0 sort!(collect(filter(isunit, ⊗(uncoupled[1]...)))) @@ -148,10 +156,12 @@ function FusionTreeBlock{I}( sort!(collect(intersect(⊗(uncoupled[1]...), ⊗(uncoupled[2]...)))) end - for c in cs, - f₂ in fusiontrees(uncoupled[2], c, isdual[2]), - f₁ in fusiontrees(uncoupled[1], c, isdual[1]) - push!(trees, (f₁, f₂)) + for c in cs + isempty(fusiontrees(uncoupled[1], c, isdual[1])) && continue + for f₂ in fusiontrees(uncoupled[2], c, isdual[2]), + f₁ in fusiontrees(uncoupled[1], c, isdual[1]) + push!(trees, (f₁, f₂)) + end end return FusionTreeBlock(trees) end From 4f33f8d51d5d890a297cc21ee5b3678882d7a676 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 14 Oct 2026 12:59:06 +0200 Subject: [PATCH 02/12] Group the cached fusion trees in `fusionblocks` Instead of constructing a block for every combination of codomain and domain sectors, `fusionblocks` groups the cached `fusiontrees(W)` by their uncoupled sectors, in the same order as before. Co-Authored-By: Claude Opus 5.5 --- src/spaces/homspace.jl | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index b41bd45be..45631db92 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -210,13 +210,18 @@ grouped by their uncoupled charges. function fusionblocks(W::HomSpace) I = sectortype(W) N₁, N₂ = numout(W), numin(W) - isdual_src = (map(isdual, codomain(W)), map(isdual, domain(W))) - fblocks = Vector{FusionTreeBlock{I, N₁, N₂, fusiontreetype(I, N₁, N₂)}}() - for dom_uncoupled_src in sectors(domain(W)), cod_uncoupled_src in sectors(codomain(W)) - fs_src = FusionTreeBlock{I}((cod_uncoupled_src, dom_uncoupled_src), isdual_src) - isempty(fs_src) || push!(fblocks, fs_src) + F = fusiontreetype(I, N₁, N₂) + # group the (cached) fusion trees of `W` by their uncoupled sectors + groups = Dict{Tuple{NTuple{N₁, I}, NTuple{N₂, I}}, Vector{F}}() + for (f₁, f₂) in fusiontrees(W) + push!(get!(Vector{F}, groups, (f₁.uncoupled, f₂.uncoupled)), (f₁, f₂)) end - return fblocks + # in the order of `sectors(domain(W))` (outer) and `sectors(codomain(W))` (inner), each of + # which runs fastest over its first space + position(V) = Dict{I, Int}(s => i for (i, s) in enumerate(sectors(V))) + codpos, dompos = map(position, codomain(W)), map(position, domain(W)) + rank((u₁, u₂)) = (reverse(map(getindex, dompos, u₂)), reverse(map(getindex, codpos, u₁))) + return map(u -> FusionTreeBlock(groups[u]), sort!(collect(keys(groups)); by = rank)) end function diagonalblockstructure(W::HomSpace) From e2b7e99563a69e8109bda245cfd2d507f8f471b1 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 14 Oct 2026 14:47:29 +0200 Subject: [PATCH 03/12] Test the blocks of a HomSpace; changelog Co-Authored-By: Claude Opus 5.5 --- docs/src/Changelog.md | 1 + test/symmetries/doubletree.jl | 15 +++++++++++++++ 2 files changed, 16 insertions(+) diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md index 2a853124a..b91e0e68b 100644 --- a/docs/src/Changelog.md +++ b/docs/src/Changelog.md @@ -42,6 +42,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec ### Performance - Recoupling matrices are stored such that they can leverage BLAS, also for mixed complex tensor with real recoupling coefficients. +- Constructing a `FusionTreeBlock` only takes the fusion product of the shorter side, and `fusionblocks` groups the cached fusion trees of the space instead of constructing a block for every combination of sectors (internal) ## [0.17.2](https://github.com/QuantumKitHub/TensorKit.jl/compare/v0.17.1...v0.17.2) - 2026-09-20 diff --git a/test/symmetries/doubletree.jl b/test/symmetries/doubletree.jl index 4fbe43e90..908e0497a 100644 --- a/test/symmetries/doubletree.jl +++ b/test/symmetries/doubletree.jl @@ -235,6 +235,21 @@ using TensorKitSectors end end + UnitStyle(I) isa SimpleUnit && @testset "Double fusion tree: blocks of a HomSpace" begin + V = Vect[I](Dict(s => 1 for s in (out..., incoming))) + W = V ⊗ V' ← V ⊗ V' # never empty: the unit couples both sides + trees = collect(fusiontrees(W)) + blocks = TK.fusionblocks(W) + # every block holds the trees of `W` with its uncoupled sectors, in the same order + for block in blocks + @test fusiontrees(block) == filter(f -> (f[1].uncoupled, f[2].uncoupled) == block.uncoupled, trees) + end + @test sum(length, blocks) == length(trees) + # the blocks are ordered as `sectors(domain(W))` (outer) and `sectors(codomain(W))` (inner) + uncoupled = [(u₁, u₂) for u₂ in sectors(domain(W)) for u₁ in sectors(codomain(W))] + @test map(block -> block.uncoupled, blocks) == filter(u -> any(f -> (f[1].uncoupled, f[2].uncoupled) == u, trees), uncoupled) + end + @testset "Double fusion tree: planar trace" begin if FusionStyle(I) isa UniqueFusion f1, f1 = src From 067d0dabcab2ad141d7af375445d965b38b9a920 Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 15 Oct 2026 13:22:59 +0200 Subject: [PATCH 04/12] Compute the rank of each block once when sorting in `fusionblocks` `sort!(...; by = rank)` evaluated `rank`, and with it the `Dict` lookups of (product) sectors, at every comparison. Co-Authored-By: Claude Opus 5.5 --- src/spaces/homspace.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index 45631db92..dc5bd589b 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -221,7 +221,8 @@ function fusionblocks(W::HomSpace) position(V) = Dict{I, Int}(s => i for (i, s) in enumerate(sectors(V))) codpos, dompos = map(position, codomain(W)), map(position, domain(W)) rank((u₁, u₂)) = (reverse(map(getindex, dompos, u₂)), reverse(map(getindex, codpos, u₁))) - return map(u -> FusionTreeBlock(groups[u]), sort!(collect(keys(groups)); by = rank)) + us = collect(keys(groups)) + return map(i -> FusionTreeBlock(groups[us[i]]), sortperm(map(rank, us))) # each rank computed once end function diagonalblockstructure(W::HomSpace) From 485afaf28aca4ef7ac6ad1e135b9b207422c0317 Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 15 Oct 2026 13:25:51 +0200 Subject: [PATCH 05/12] Use the fusion product of the shorter side for every UnitStyle; check only the longer side for trees Co-Authored-By: Claude Opus 5.5 --- src/fusiontrees/fusiontrees.jl | 23 +++++++++-------------- 1 file changed, 9 insertions(+), 14 deletions(-) diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index 9c21de89a..23c466f11 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/src/fusiontrees/fusiontrees.jl @@ -140,24 +140,19 @@ function FusionTreeBlock{I}( end # note: cs is unique due to ⊗ returning unique fusion results - cs = if UnitStyle(I) isa SimpleUnit - # a superset of the coupled sectors suffices, as there are no fusion trees for the others: - # this avoids the fusion product of the longer side - if N₁ == 0 || N₂ == 0 - [unit(I)] - else - sort!(vec(collect(⊗(uncoupled[N₁ <= N₂ ? 1 : 2]...)))) # `vec`: products of `ProductSector`s have a shape - end - elseif N₁ == 0 - sort!(collect(filter(isunit, ⊗(uncoupled[2]...)))) - elseif N₂ == 0 - sort!(collect(filter(isunit, ⊗(uncoupled[1]...)))) + # a superset of the coupled sectors suffices, as there are no fusion trees for the others: + # the units, or the fusion product of the shorter side, which avoids that of the longer side + # (`vec`: the units and products of `ProductSector`s have a shape) + cs = if N₁ == 0 || N₂ == 0 + sort!(vec(collect(allunits(I)))) else - sort!(collect(intersect(⊗(uncoupled[1]...), ⊗(uncoupled[2]...)))) + sort!(vec(collect(⊗(uncoupled[N₁ <= N₂ ? 1 : 2]...)))) end for c in cs - isempty(fusiontrees(uncoupled[1], c, isdual[1])) && continue + # the shorter side has fusion trees for every candidate, the longer side not necessarily + longer = N₁ <= N₂ ? fusiontrees(uncoupled[2], c, isdual[2]) : fusiontrees(uncoupled[1], c, isdual[1]) + isempty(longer) && continue for f₂ in fusiontrees(uncoupled[2], c, isdual[2]), f₁ in fusiontrees(uncoupled[1], c, isdual[1]) push!(trees, (f₁, f₂)) From 994eab41ca30c2f33c7df401b5be03b8ed136542 Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 15 Oct 2026 20:37:14 +0200 Subject: [PATCH 06/12] Use the unit of the other side as the coupled sector of a block with no legs on one side Co-Authored-By: Claude Opus 5.5 --- src/fusiontrees/fusiontrees.jl | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index 23c466f11..e8bb91f40 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/src/fusiontrees/fusiontrees.jl @@ -141,10 +141,12 @@ function FusionTreeBlock{I}( # note: cs is unique due to ⊗ returning unique fusion results # a superset of the coupled sectors suffices, as there are no fusion trees for the others: - # the units, or the fusion product of the shorter side, which avoids that of the longer side - # (`vec`: the units and products of `ProductSector`s have a shape) - cs = if N₁ == 0 || N₂ == 0 - sort!(vec(collect(allunits(I)))) + # the unit of the other side, or the fusion product of the shorter side, which avoids that of + # the longer side (`vec`: products of `ProductSector`s have a shape) + cs = if N₁ == 0 + [leftunit(first(uncoupled[2]))] + elseif N₂ == 0 + [leftunit(first(uncoupled[1]))] else sort!(vec(collect(⊗(uncoupled[N₁ <= N₂ ? 1 : 2]...)))) end From 8e8029aebb24a009bcfe96cacbcfe6788e065bbd Mon Sep 17 00:00:00 2001 From: leburgel Date: Thu, 15 Oct 2026 22:21:32 +0200 Subject: [PATCH 07/12] Only try coupled sectors whose units match the ends of the other side For several units, a coupled sector of a block must be compatible with the units at both ends of either side; otherwise the fusion tree iterator throws instead of being empty. Co-Authored-By: Claude Opus 5.5 --- src/fusiontrees/fusiontrees.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index e8bb91f40..6dcfd2bdf 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/src/fusiontrees/fusiontrees.jl @@ -143,16 +143,20 @@ function FusionTreeBlock{I}( # a superset of the coupled sectors suffices, as there are no fusion trees for the others: # the unit of the other side, or the fusion product of the shorter side, which avoids that of # the longer side (`vec`: products of `ProductSector`s have a shape) - cs = if N₁ == 0 - [leftunit(first(uncoupled[2]))] - elseif N₂ == 0 - [leftunit(first(uncoupled[1]))] + cs = if N₁ == 0 || N₂ == 0 + # the only unit that can couple: on the left of the first and on the right of the last sector + us = N₁ == 0 ? uncoupled[2] : uncoupled[1] + leftunit(first(us)) == rightunit(last(us)) ? [leftunit(first(us))] : I[] else sort!(vec(collect(⊗(uncoupled[N₁ <= N₂ ? 1 : 2]...)))) end for c in cs # the shorter side has fusion trees for every candidate, the longer side not necessarily + # (with several units, `c` must also have the units at the ends of the longer side) + us = N₁ <= N₂ ? uncoupled[2] : uncoupled[1] + UnitStyle(I) isa SimpleUnit || N₁ == 0 || N₂ == 0 || + (leftunit(c) == leftunit(first(us)) && rightunit(c) == rightunit(last(us))) || continue longer = N₁ <= N₂ ? fusiontrees(uncoupled[2], c, isdual[2]) : fusiontrees(uncoupled[1], c, isdual[1]) isempty(longer) && continue for f₂ in fusiontrees(uncoupled[2], c, isdual[2]), From abc00bb4e30bc90ca28ac712c7def73380fa3dfd Mon Sep 17 00:00:00 2001 From: leburgel Date: Tue, 20 Oct 2026 18:05:24 +0200 Subject: [PATCH 08/12] Order the blocks by comparing their sectors Co-Authored-By: Claude Opus 5.5 --- src/spaces/homspace.jl | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index dc5bd589b..bae972f4e 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -217,12 +217,10 @@ function fusionblocks(W::HomSpace) push!(get!(Vector{F}, groups, (f₁.uncoupled, f₂.uncoupled)), (f₁, f₂)) end # in the order of `sectors(domain(W))` (outer) and `sectors(codomain(W))` (inner), each of - # which runs fastest over its first space - position(V) = Dict{I, Int}(s => i for (i, s) in enumerate(sectors(V))) - codpos, dompos = map(position, codomain(W)), map(position, domain(W)) - rank((u₁, u₂)) = (reverse(map(getindex, dompos, u₂)), reverse(map(getindex, codpos, u₁))) - us = collect(keys(groups)) - return map(i -> FusionTreeBlock(groups[us[i]]), sortperm(map(rank, us))) # each rank computed once + # which runs fastest over its first space (the sectors of a dual space are sorted by their dual) + key(V, s) = isdual(V) ? dual(s) : s + rank((u₁, u₂)) = (reverse(map(key, domain(W).spaces, u₂)), reverse(map(key, codomain(W).spaces, u₁))) + return map(u -> FusionTreeBlock(groups[u]), sort!(collect(keys(groups)); by = rank)) end function diagonalblockstructure(W::HomSpace) From ff0af2c158490ba087d9a1f61ddf6c4e29be47f4 Mon Sep 17 00:00:00 2001 From: leburgel Date: Tue, 20 Oct 2026 21:14:36 +0200 Subject: [PATCH 09/12] Compute the sort key of each block once Co-Authored-By: Claude Opus 5.5 --- src/spaces/homspace.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index bae972f4e..30c364fca 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -220,7 +220,8 @@ function fusionblocks(W::HomSpace) # which runs fastest over its first space (the sectors of a dual space are sorted by their dual) key(V, s) = isdual(V) ? dual(s) : s rank((u₁, u₂)) = (reverse(map(key, domain(W).spaces, u₂)), reverse(map(key, codomain(W).spaces, u₁))) - return map(u -> FusionTreeBlock(groups[u]), sort!(collect(keys(groups)); by = rank)) + us = collect(keys(groups)) + return map(i -> FusionTreeBlock(groups[us[i]]), sortperm(map(rank, us))) # each rank computed once end function diagonalblockstructure(W::HomSpace) From 3929e81604fc229db53eccfce15d4e33589b2343 Mon Sep 17 00:00:00 2001 From: leburgel Date: Wed, 28 Oct 2026 10:59:01 +0100 Subject: [PATCH 10/12] Keep the blocks in the order in which their uncoupled sectors first occur Co-Authored-By: Claude Opus 5.5 --- src/spaces/homspace.jl | 16 +++++++--------- test/symmetries/doubletree.jl | 4 +--- 2 files changed, 8 insertions(+), 12 deletions(-) diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index 30c364fca..dfd775052 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -211,17 +211,15 @@ function fusionblocks(W::HomSpace) I = sectortype(W) N₁, N₂ = numout(W), numin(W) F = fusiontreetype(I, N₁, N₂) - # group the (cached) fusion trees of `W` by their uncoupled sectors - groups = Dict{Tuple{NTuple{N₁, I}, NTuple{N₂, I}}, Vector{F}}() + # group the (cached) fusion trees of `W` by their uncoupled sectors, in the order in which + # these first occur + groups = Vector{F}[] + group_of = Dict{Tuple{NTuple{N₁, I}, NTuple{N₂, I}}, Vector{F}}() for (f₁, f₂) in fusiontrees(W) - push!(get!(Vector{F}, groups, (f₁.uncoupled, f₂.uncoupled)), (f₁, f₂)) + group = get!(() -> push!(groups, F[])[end], group_of, (f₁.uncoupled, f₂.uncoupled)) + push!(group, (f₁, f₂)) end - # in the order of `sectors(domain(W))` (outer) and `sectors(codomain(W))` (inner), each of - # which runs fastest over its first space (the sectors of a dual space are sorted by their dual) - key(V, s) = isdual(V) ? dual(s) : s - rank((u₁, u₂)) = (reverse(map(key, domain(W).spaces, u₂)), reverse(map(key, codomain(W).spaces, u₁))) - us = collect(keys(groups)) - return map(i -> FusionTreeBlock(groups[us[i]]), sortperm(map(rank, us))) # each rank computed once + return map(FusionTreeBlock, groups) end function diagonalblockstructure(W::HomSpace) diff --git a/test/symmetries/doubletree.jl b/test/symmetries/doubletree.jl index 908e0497a..180e19bdc 100644 --- a/test/symmetries/doubletree.jl +++ b/test/symmetries/doubletree.jl @@ -245,9 +245,7 @@ using TensorKitSectors @test fusiontrees(block) == filter(f -> (f[1].uncoupled, f[2].uncoupled) == block.uncoupled, trees) end @test sum(length, blocks) == length(trees) - # the blocks are ordered as `sectors(domain(W))` (outer) and `sectors(codomain(W))` (inner) - uncoupled = [(u₁, u₂) for u₂ in sectors(domain(W)) for u₁ in sectors(codomain(W))] - @test map(block -> block.uncoupled, blocks) == filter(u -> any(f -> (f[1].uncoupled, f[2].uncoupled) == u, trees), uncoupled) + @test allunique(map(block -> block.uncoupled, blocks)) end @testset "Double fusion tree: planar trace" begin From 6993e67d7840cece2628c217c993ee2f1ef91132 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 28 Oct 2026 15:19:34 -0400 Subject: [PATCH 11/12] Fill the vectors of the fusion tree blocks directly Co-Authored-By: Claude Opus 5.5 --- src/spaces/homspace.jl | 18 +++++++++++------- 1 file changed, 11 insertions(+), 7 deletions(-) diff --git a/src/spaces/homspace.jl b/src/spaces/homspace.jl index dfd775052..38bf7eb60 100644 --- a/src/spaces/homspace.jl +++ b/src/spaces/homspace.jl @@ -211,15 +211,19 @@ function fusionblocks(W::HomSpace) I = sectortype(W) N₁, N₂ = numout(W), numin(W) F = fusiontreetype(I, N₁, N₂) - # group the (cached) fusion trees of `W` by their uncoupled sectors, in the order in which - # these first occur - groups = Vector{F}[] - group_of = Dict{Tuple{NTuple{N₁, I}, NTuple{N₂, I}}, Vector{F}}() + # group the (cached) fusion trees of `W` by their uncoupled sectors, + # in the order in which these first occur + blocks = FusionTreeBlock{I, N₁, N₂, F}[] + trees_of = Dict{Tuple{NTuple{N₁, I}, NTuple{N₂, I}}, Vector{F}}() for (f₁, f₂) in fusiontrees(W) - group = get!(() -> push!(groups, F[])[end], group_of, (f₁.uncoupled, f₂.uncoupled)) - push!(group, (f₁, f₂)) + trees = get!(trees_of, (f₁.uncoupled, f₂.uncoupled)) do + newtrees = F[] + push!(blocks, FusionTreeBlock(newtrees)) + return newtrees + end + push!(trees, (f₁, f₂)) end - return map(FusionTreeBlock, groups) + return blocks end function diagonalblockstructure(W::HomSpace) From 3f74b72a6f8b4510637b6912ff5f4cc632fd84d8 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 28 Oct 2026 15:19:34 -0400 Subject: [PATCH 12/12] Reuse the fusion tree iterator of the longer side Co-Authored-By: Claude Opus 5.5 --- src/fusiontrees/fusiontrees.jl | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl index 6dcfd2bdf..40a2961b6 100644 --- a/src/fusiontrees/fusiontrees.jl +++ b/src/fusiontrees/fusiontrees.jl @@ -157,10 +157,16 @@ function FusionTreeBlock{I}( us = N₁ <= N₂ ? uncoupled[2] : uncoupled[1] UnitStyle(I) isa SimpleUnit || N₁ == 0 || N₂ == 0 || (leftunit(c) == leftunit(first(us)) && rightunit(c) == rightunit(last(us))) || continue - longer = N₁ <= N₂ ? fusiontrees(uncoupled[2], c, isdual[2]) : fusiontrees(uncoupled[1], c, isdual[1]) - isempty(longer) && continue - for f₂ in fusiontrees(uncoupled[2], c, isdual[2]), - f₁ in fusiontrees(uncoupled[1], c, isdual[1]) + if N₁ <= N₂ + ftrees₂ = fusiontrees(uncoupled[2], c, isdual[2]) + isempty(ftrees₂) && continue + ftrees₁ = fusiontrees(uncoupled[1], c, isdual[1]) + else + ftrees₁ = fusiontrees(uncoupled[1], c, isdual[1]) + isempty(ftrees₁) && continue + ftrees₂ = fusiontrees(uncoupled[2], c, isdual[2]) + end + for f₂ in ftrees₂, f₁ in ftrees₁ push!(trees, (f₁, f₂)) end end