From 843af52c7ba703ce1070404b4d62f4dd99f7f31d Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Wed, 14 Oct 2026 12:59:05 +0200
Subject: [PATCH 01/13] Keep the intermediate blocks of `fsbraid` in production
 order
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

The bends and Artin braids of a fusion tree block used to construct their destination block in canonical order, together with a dictionary from trees to rows, at every step. The intermediate blocks of `fsbraid` now keep their trees in the order in which they are produced (`TreeRows`, `canonical = Val(false)` for `bendleft`, `bendright` and `artin_braid`, `repartition_noncanonical`), and only the result is brought into canonical order (`canonicalize`). Small blocks find rows by a linear search. In `artin_braid`, the channels are obtained by filtering `a ⊗ d` instead of intersecting two fusion products, and the symbols that do not depend on the channel are computed once.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/braiding_manipulations.jl | 63 +++++++++++++----------
 src/fusiontrees/duality_manipulations.jl  | 41 +++++++++------
 src/fusiontrees/fusiontrees.jl            | 49 ++++++++++++++++++
 3 files changed, 110 insertions(+), 43 deletions(-)

diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index b737f354f..e61a66f44 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -78,7 +78,16 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N}
     return f′ => coeff
 end
 
-function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where {I, N}
+# the inner lines `c′` of `a ⊗ d → c′ → c′ ⊗ b → e`, i.e. `intersect(a ⊗ d, e ⊗ conj(b))`,
+# without constructing the second fusion product for sectors with a simple unit
+function _artin_channels(a::I, b::I, d::I, e::I) where {I <: Sector}
+    UnitStyle(I) isa SimpleUnit || return intersect(a ⊗ d, e ⊗ conj(b))
+    return Iterators.filter(c′ -> !iszero(Nsymbol(c′, b, e)), a ⊗ d)
+end
+
+function artin_braid(
+        src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false, canonical::Val = Val(true)
+    ) where {I, N}
     1 <= i < N || throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs"))
     isempty(src) && return src => zeros(sectorscalartype(I), 0, 0)
 
@@ -91,12 +100,11 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
     isdual = src.isdual[1]
     isdual′ = TupleTools.setindex(isdual, isdual[i], i + 1)
     isdual′ = TupleTools.setindex(isdual′, isdual[i + 1], i)
-    dst = FusionTreeBlock{I}((uncoupled′, ()), (isdual′, ()); sizehint = length(src))
+    rows = TreeRows(I, (uncoupled′, ()), (isdual′, ()), length(src), canonical)
 
     oneT = one(sectorscalartype(I))
 
-    indexmap = treeindex_map(dst)
-    U = zeros(sectorscalartype(I), length(dst), length(src))
+    U = zeros(sectorscalartype(I), max(length(rows.trees), length(src)), length(src)) # see `bendright`
 
     if isunit(a) || isunit(b) # braiding with trivial sector: simple and always possible
         for (col, (f, f₂)) in enumerate(fusiontrees(src))
@@ -113,10 +121,10 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
                 vertices′ = TupleTools.setindex(vertices′, vertices[i - 1], i)
             end
             f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′)
-            row = indexmap[treeindex_data((f′, f₂))]
+            row = row!(rows, (f′, f₂))
             @inbounds U[row, col] = oneT
         end
-        return dst => U
+        return FusionTreeBlock(rows.trees) => U
     end
 
     BraidingStyle(I) isa NoBraiding &&
@@ -132,7 +140,7 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
             if FusionStyle(I) isa MultiplicityFreeFusion
                 R = oftype(oneT, (inv ? conj(Rsymbol(b, a, c)) : Rsymbol(a, b, c)))
                 f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices)
-                row = indexmap[treeindex_data((f′, f₂))]
+                row = row!(rows, (f′, f₂))
                 @inbounds U[row, col] = R
             else # GenericFusion
                 μ = vertices[1]
@@ -142,7 +150,7 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
                     iszero(R) && continue
                     vertices′ = TupleTools.setindex(vertices, ν, 1)
                     f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices′)
-                    row = indexmap[treeindex_data((f′, f₂))]
+                    row = row!(rows, (f′, f₂))
                     @inbounds U[row, col] = R
                 end
             end
@@ -155,29 +163,30 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
         c = inner_extended[i]
         e = inner_extended[i + 1]
         if FusionStyle(I) isa MultiplicityFreeFusion
-            for c′ in intersect(a ⊗ d, e ⊗ conj(b))
+            R = inv ? Rsymbol(d, c, e) : Rsymbol(c, d, e)
+            for c′ in _artin_channels(a, b, d, e)
                 coeff = if inv
-                    conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′)
+                    conj(R * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′)
                 else
-                    Rsymbol(c, d, e) * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′))
+                    R * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′))
                 end
                 iszero(coeff) && continue
                 inner′ = TupleTools.setindex(inner, c′, i - 1)
                 f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′)
-                row = indexmap[treeindex_data((f′, f₂))]
+                row = row!(rows, (f′, f₂))
                 @inbounds U[row, col] = coeff
             end
         else # GenericFusion
-            for c′ in intersect(a ⊗ d, e ⊗ conj(b))
-                Rmat1 = inv ? Rsymbol(d, c, e)' : Rsymbol(c, d, e)
+            Rmat1 = inv ? Rsymbol(d, c, e)' : Rsymbol(c, d, e)
+            μ = vertices[i - 1]
+            ν = vertices[i]
+            for c′ in _artin_channels(a, b, d, e)
                 Rmat2 = inv ? Rsymbol(d, a, c′)' : Rsymbol(a, d, c′)
                 Fmat = Fsymbol(d, a, b, e, c′, c)
-                μ = vertices[i - 1]
-                ν = vertices[i]
-                for σ in 1:Nsymbol(a, d, c′)
-                    for λ in 1:Nsymbol(c′, b, e)
+                for σ in axes(Rmat2, 1)
+                    for λ in axes(Fmat, 2)
                         coeff = zero(oneT)
-                        for ρ in 1:Nsymbol(d, c, e), κ in 1:Nsymbol(d, a, c′)
+                        for ρ in axes(Rmat1, 2), κ in axes(Rmat2, 2)
                             coeff += Rmat1[ν, ρ] * conj(Fmat[κ, λ, μ, ρ]) *
                                 conj(Rmat2[σ, κ])
                         end
@@ -186,7 +195,7 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
                         vertices′ = TupleTools.setindex(vertices′, λ, i)
                         inner′ = TupleTools.setindex(inner, c′, i - 1)
                         f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′)
-                        row = indexmap[treeindex_data((f′, f₂))]
+                        row = row!(rows, (f′, f₂))
                         @inbounds U[row, col] = coeff
                     end
                 end
@@ -194,7 +203,7 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
         end
     end
 
-    return dst => U
+    return FusionTreeBlock(rows.trees) => U
 end
 
 # braid fusion tree
@@ -323,13 +332,15 @@ end
         p = linearizepermutation(p1, p2, numout(src), numin(src))
         levels = (l1..., reverse(l2)...)
 
-        dst, U = repartition(src, numind(src))
+        # the intermediate blocks keep their trees in the order in which they are produced,
+        # and only the result is brought into canonical order
+        dst, U = repartition_noncanonical(src, Val(numind(src)))
 
         braid_style = BraidingStyle(I)
         for s in permutation2swaps(p)
             _check_levels(braid_style, levels, s)
             inv = levels[s] > levels[s + 1]
-            dst, U_tmp = artin_braid(dst, s; inv)
+            dst, U_tmp = artin_braid(dst, s; inv, canonical = Val(false))
             U = U_tmp * U
             l = levels[s]
             levels = TupleTools.setindex(levels, levels[s + 1], s)
@@ -337,11 +348,11 @@ end
         end
 
         if N₂ == 0
-            return dst => U
+            return canonicalize(dst, U)
         else
-            dst, U_tmp = repartition(dst::FusionTreeBlock{I}, N₁)
+            dst, U_tmp = repartition_noncanonical(dst::FusionTreeBlock{I}, Val(N₁))
             U = U_tmp * U
-            return dst => U
+            return canonicalize(dst, U)
         end
     end
 end
diff --git a/src/fusiontrees/duality_manipulations.jl b/src/fusiontrees/duality_manipulations.jl
index 9e613d50e..48cb106a2 100644
--- a/src/fusiontrees/duality_manipulations.jl
+++ b/src/fusiontrees/duality_manipulations.jl
@@ -66,7 +66,7 @@ function bendright((f₁, f₂)::FusionTreePair)
 
     return (f₁′, f₂′) => coeff
 end
-function bendright(src::FusionTreeBlock)
+function bendright(src::FusionTreeBlock, canonical::Val = Val(true))
     I = sectortype(src)
     N₁ = numout(src)
     N₂ = numin(src)
@@ -80,9 +80,10 @@ function bendright(src::FusionTreeBlock)
         (src.isdual[2]..., !(src.isdual[1][N₁])),
     )
 
-    dst = FusionTreeBlock{I}(uncoupled_dst, isdual_dst; sizehint = length(src))
-    indexmap = treeindex_map(dst)
-    U = zeros(fusionscalartype(I), length(dst), length(src))
+    rows = TreeRows(I, uncoupled_dst, isdual_dst, length(src), canonical)
+    # the canonical destination block has more trees than `src` if `src` holds only some of the
+    # trees with its uncoupled sectors (as in `planar_trace`)
+    U = zeros(fusionscalartype(I), max(length(rows.trees), length(src)), length(src))
 
     for (col, (f₁, f₂)) in enumerate(fusiontrees(src))
         (a, b, c), (f₁′, f₂′) = _bendright_treepair((f₁, f₂))
@@ -92,7 +93,7 @@ function bendright(src::FusionTreeBlock)
         end
         if FusionStyle(I) isa MultiplicityFreeFusion
             coeff = coeff₀ * Bsymbol(a, b, c)
-            row = indexmap[treeindex_data((f₁′, f₂′))]
+            row = row!(rows, (f₁′, f₂′))
             @inbounds U[row, col] = coeff
         else
             Bmat = Bsymbol(a, b, c)
@@ -106,12 +107,12 @@ function bendright(src::FusionTreeBlock)
                 iszero(coeff) && continue
                 vertices₂ = N₂ > 0 ? (f₂.vertices..., ν) : ()
                 f₂′ = FusionTree(uncoupled₂, coupled₂, isdual₂, inner₂, vertices₂)
-                row = indexmap[treeindex_data((f₁′, f₂′))]
+                row = row!(rows, (f₁′, f₂′))
                 @inbounds U[row, col] = coeff
             end
         end
     end
-    return dst => U
+    return FusionTreeBlock(rows.trees) => U
 end
 
 @doc """
@@ -145,7 +146,7 @@ end
 
 # !! note that this is more or less a copy of bendright through
 # (f1, f2) => conj(coeff) for ((f2, f1), coeff) in bendleft(src)
-function bendleft(src::FusionTreeBlock)
+function bendleft(src::FusionTreeBlock, canonical::Val = Val(true))
     I = sectortype(src)
     N₁ = numout(src)
     N₂ = numin(src)
@@ -159,9 +160,10 @@ function bendleft(src::FusionTreeBlock)
         TupleTools.front(src.isdual[2]),
     )
 
-    dst = FusionTreeBlock{I}(uncoupled_dst, isdual_dst; sizehint = length(src))
-    indexmap = treeindex_map(dst)
-    U = zeros(fusionscalartype(I), length(dst), length(src))
+    rows = TreeRows(I, uncoupled_dst, isdual_dst, length(src), canonical)
+    # the canonical destination block has more trees than `src` if `src` holds only some of the
+    # trees with its uncoupled sectors (as in `planar_trace`)
+    U = zeros(fusionscalartype(I), max(length(rows.trees), length(src)), length(src))
 
     for (col, (f₁, f₂)) in enumerate(fusiontrees(src))
         (a, b, c), (f₂′, f₁′) = _bendright_treepair((f₂, f₁))
@@ -171,7 +173,7 @@ function bendleft(src::FusionTreeBlock)
         end
         if FusionStyle(I) isa MultiplicityFreeFusion
             coeff = coeff₀ * Bsymbol(a, b, c)
-            row = indexmap[treeindex_data((f₁′, f₂′))]
+            row = row!(rows, (f₁′, f₂′))
             @inbounds U[row, col] = conj(coeff)
         else
             Bmat = Bsymbol(a, b, c)
@@ -185,12 +187,12 @@ function bendleft(src::FusionTreeBlock)
                 iszero(coeff) && continue
                 vertices₁ = N₁ > 0 ? (f₁.vertices..., ν) : ()
                 f₁′ = FusionTree(uncoupled₁, coupled₁, isdual₁, inner₁, vertices₁)
-                row = indexmap[treeindex_data((f₁′, f₂′))]
+                row = row!(rows, (f₁′, f₂′))
                 @inbounds U[row, col] = conj(coeff)
             end
         end
     end
-    return dst => U
+    return FusionTreeBlock(rows.trees) => U
 end
 
 @doc """
@@ -478,7 +480,7 @@ return dst, U
 # it, so that it is visible to the generator's world age. Otherwise a precompile workload that
 # first triggers `repartition` fails with `UndefVarError: _repartition_body` (at runtime the
 # first call happens at a later world age, which masks the ordering issue).
-function _repartition_body(N)
+function _repartition_body(N, canonical = true)
     if N == 0
         ex = quote
             T = fusionscalartype(sectortype(src))
@@ -491,13 +493,14 @@ function _repartition_body(N)
         end
     else
         f = N < 0 ? bendleft : bendright
+        args = canonical ? () : (Val(false),)
         ex_rep = Expr(:block)
         for _ in 1:(abs(N) - 1)
-            push!(ex_rep.args, :((dst, Utmp) = $f(dst)))
+            push!(ex_rep.args, :((dst, Utmp) = $f(dst, $(args...))))
             push!(ex_rep.args, :(U = Utmp * U))
         end
         ex = quote
-            dst, U = $f(src)
+            dst, U = $f(src, $(args...))
             $ex_rep
             return dst => U
         end
@@ -507,6 +510,10 @@ end
 @generated function repartition(src::Union{FusionTreePair, FusionTreeBlock}, ::Val{N}) where {N}
     return _repartition_body(numout(src) - N)
 end
+# the same for a block, with the trees of the result in the order in which they are produced
+@generated function repartition_noncanonical(src::FusionTreeBlock, ::Val{N}) where {N}
+    return _repartition_body(numout(src) - N, false)
+end
 
 """
     transpose((f₁, f₂)::FusionTreePair{I}, p::Index2Tuple{N₁, N₂}) where {I, N₁, N₂}
diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl
index 4ea3ab915..a08156022 100644
--- a/src/fusiontrees/fusiontrees.jl
+++ b/src/fusiontrees/fusiontrees.jl
@@ -117,6 +117,8 @@ const FusionTreePair{I, N₁, N₂} = Tuple{FusionTree{I, N₁}, FusionTree{I, N
 
 Collection of fusion-splitting tree pairs that share the same uncoupled sectors.
 Mostly internal structure to speed up non-`UniqueFusion` fusiontree manipulation computations.
+Braiding a block requires it to be complete, i.e. to contain all the fusion-splitting tree pairs
+with its uncoupled sectors and duality flags, as `FusionTreeBlock{I}(uncoupled, isdual)` does.
 """
 struct FusionTreeBlock{I, N₁, N₂, F <: FusionTreePair{I, N₁, N₂}}
     trees::Vector{F}
@@ -283,6 +285,53 @@ function treeindex_map(fs::FusionTreeBlock)
     return fusiontreedict(I)(treeindex_data(f) => ind for (ind, f) in enumerate(fusiontrees(fs)))
 end
 
+# The trees of the destination block of a block manipulation, and the row of each of them in the
+# transformation matrix. With `Val(true)`, these are the trees of the destination block in their
+# canonical order; with `Val(false)`, the trees are added in the order in which they are produced,
+# which avoids constructing the destination block. The latter is used for the intermediate blocks
+# of a composite manipulation, of which only the result is brought into canonical order.
+# With `Val(false)`, the source block has to be complete: the destination block then has as many
+# trees (`n`), and an incomplete source block raises an error if its destination block has more.
+struct TreeRows{F, D}
+    trees::Vector{F}
+    index::Union{D, Nothing} # `nothing`: few trees, the rows are found by a linear search
+    n::Int
+end
+function TreeRows(::Type{I}, uncoupled, isdual, n::Int, ::Val{canonical}) where {I, canonical}
+    trees = if canonical
+        fusiontrees(FusionTreeBlock{I}(uncoupled, isdual; sizehint = n))
+    else
+        sizehint!(Vector{fusiontreetype(I, length(uncoupled[1]), length(uncoupled[2]))}(), n)
+    end
+    K = Base.promote_op(treeindex_data, eltype(trees))
+    index = n > 32 ? Dict{K, Int}(treeindex_data(f) => i for (i, f) in enumerate(trees)) : nothing
+    return TreeRows{eltype(trees), Dict{K, Int}}(trees, index, n)
+end
+function row!(rows::TreeRows, f)
+    key = treeindex_data(f)
+    i = if isnothing(rows.index)
+        findfirst(f′ -> treeindex_data(f′) == key, rows.trees)
+    else
+        get(rows.index, key, nothing)
+    end
+    isnothing(i) || return i
+    length(rows.trees) < rows.n ||
+        throw(ArgumentError("braiding requires a complete `FusionTreeBlock`"))
+    push!(rows.trees, f)
+    isnothing(rows.index) || (rows.index[key] = length(rows.trees))
+    return length(rows.trees)
+end
+
+# bring a block whose trees are not in canonical order, together with the transformation matrix
+# whose rows correspond to these trees, into canonical order
+function canonicalize(src::FusionTreeBlock{I}, U::AbstractMatrix) where {I}
+    length(src) == 1 && return src => U
+    rows = TreeRows(I, src.uncoupled, src.isdual, length(src), Val(true))
+    U′ = zeros(eltype(U), length(rows.trees), size(U, 2))
+    U′[map(f -> row!(rows, f), fusiontrees(src)), :] = U
+    return FusionTreeBlock(rows.trees) => U′
+end
+
 
 # Facilitate getting correct fusion tree types
 Base.@assume_effects :foldable function fusiontreetype(::Type{I}, N::Int) where {I <: Sector}

From bd99713eb04f1c2ec64a36704b479ad8f030a3a1 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Wed, 14 Oct 2026 12:59:05 +0200
Subject: [PATCH 02/13] Recycle tree vectors and matrices between the braids in
 `fsbraid`

The Artin braids keep the size of a block, so `fsbraid` passes the tree vector and the matrices of the previous step on to `artin_braid` and multiplies in place.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/braiding_manipulations.jl | 18 +++++++++++++-----
 src/fusiontrees/fusiontrees.jl            | 10 +++++++---
 2 files changed, 20 insertions(+), 8 deletions(-)

diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index e61a66f44..6a1de9d69 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -86,7 +86,8 @@ function _artin_channels(a::I, b::I, d::I, e::I) where {I <: Sector}
 end
 
 function artin_braid(
-        src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false, canonical::Val = Val(true)
+        src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false, canonical::Val = Val(true),
+        buffers = (nothing, nothing) # a tree vector and a matrix that may be overwritten
     ) where {I, N}
     1 <= i < N || throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs"))
     isempty(src) && return src => zeros(sectorscalartype(I), 0, 0)
@@ -100,11 +101,12 @@ function artin_braid(
     isdual = src.isdual[1]
     isdual′ = TupleTools.setindex(isdual, isdual[i], i + 1)
     isdual′ = TupleTools.setindex(isdual′, isdual[i + 1], i)
-    rows = TreeRows(I, (uncoupled′, ()), (isdual′, ()), length(src), canonical)
+    rows = TreeRows(I, (uncoupled′, ()), (isdual′, ()), length(src), canonical, buffers[1])
 
     oneT = one(sectorscalartype(I))
 
-    U = zeros(sectorscalartype(I), max(length(rows.trees), length(src)), length(src)) # see `bendright`
+    U = buffers[2] isa Matrix{typeof(oneT)} ? fill!(buffers[2], zero(oneT)) :
+        zeros(sectorscalartype(I), max(length(rows.trees), length(src)), length(src)) # see `bendright`
 
     if isunit(a) || isunit(b) # braiding with trivial sector: simple and always possible
         for (col, (f, f₂)) in enumerate(fusiontrees(src))
@@ -336,12 +338,18 @@ end
         # and only the result is brought into canonical order
         dst, U = repartition_noncanonical(src, Val(numind(src)))
 
+        # the swaps keep the size of the block, so that the tree vector and the matrices of the
+        # previous step are recycled (except for the trees of `src`, which are part of the key)
+        buffers = (nothing, nothing)
+        U_spare = nothing
         braid_style = BraidingStyle(I)
         for s in permutation2swaps(p)
             _check_levels(braid_style, levels, s)
             inv = levels[s] > levels[s + 1]
-            dst, U_tmp = artin_braid(dst, s; inv, canonical = Val(false))
-            U = U_tmp * U
+            dst′, U_tmp = artin_braid(dst, s; inv, canonical = Val(false), buffers)
+            buffers = (dst.trees === src.trees ? nothing : dst.trees, U_tmp)
+            dst = dst′
+            U_spare, U = U, U_spare isa typeof(U_tmp) ? mul!(U_spare, U_tmp, U) : U_tmp * U
             l = levels[s]
             levels = TupleTools.setindex(levels, levels[s + 1], s)
             levels = TupleTools.setindex(levels, l, s + 1)
diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl
index a08156022..33faca9c4 100644
--- a/src/fusiontrees/fusiontrees.jl
+++ b/src/fusiontrees/fusiontrees.jl
@@ -289,7 +289,8 @@ end
 # transformation matrix. With `Val(true)`, these are the trees of the destination block in their
 # canonical order; with `Val(false)`, the trees are added in the order in which they are produced,
 # which avoids constructing the destination block. The latter is used for the intermediate blocks
-# of a composite manipulation, of which only the result is brought into canonical order.
+# of a composite manipulation, of which only the result is brought into canonical order; the tree
+# vector of an intermediate block that is no longer needed can then be passed on as `buffer`.
 # With `Val(false)`, the source block has to be complete: the destination block then has as many
 # trees (`n`), and an incomplete source block raises an error if its destination block has more.
 struct TreeRows{F, D}
@@ -297,11 +298,14 @@ struct TreeRows{F, D}
     index::Union{D, Nothing} # `nothing`: few trees, the rows are found by a linear search
     n::Int
 end
-function TreeRows(::Type{I}, uncoupled, isdual, n::Int, ::Val{canonical}) where {I, canonical}
+function TreeRows(
+        ::Type{I}, uncoupled, isdual, n::Int, ::Val{canonical}, buffer = nothing
+    ) where {I, canonical}
+    F = fusiontreetype(I, length(uncoupled[1]), length(uncoupled[2]))
     trees = if canonical
         fusiontrees(FusionTreeBlock{I}(uncoupled, isdual; sizehint = n))
     else
-        sizehint!(Vector{fusiontreetype(I, length(uncoupled[1]), length(uncoupled[2]))}(), n)
+        buffer isa Vector{F} ? empty!(buffer) : sizehint!(Vector{F}(), n)
     end
     K = Base.promote_op(treeindex_data, eltype(trees))
     index = n > 32 ? Dict{K, Int}(treeindex_data(f) => i for (i, f) in enumerate(trees)) : nothing

From 6a47af83a257a314d0b5bbbe7d9e53cc5e27ecb4 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Wed, 14 Oct 2026 12:59:06 +0200
Subject: [PATCH 03/13] Braid single-tree blocks as a product of scalars

A block with a single tree stays a single tree under every step of `fsbraid`, so its braid is computed with the per-tree manipulations, without blocks or matrices. Their bodies move into private `_artin_braid_single` and `_bendright_single`, which take the single coefficient of each symbol with `only`; the public methods keep their `UniqueFusion` assertion.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/braiding_manipulations.jl | 30 +++++++++++++++++++----
 src/fusiontrees/duality_manipulations.jl  | 21 +++++++++++++---
 2 files changed, 43 insertions(+), 8 deletions(-)

diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index 6a1de9d69..4428e233e 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -16,8 +16,13 @@ tree with non-zero coefficient, namely `f` with coefficient `1`. This keyword ha
 if `BraidingStyle(sectortype(f)) isa SymmetricBraiding`.
 """
 function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N}
-    1 <= i < N || throw(ArgumentError(lazy"Cannot swap outputs i=$i and i+1 out of only $N outputs"))
     @assert FusionStyle(I) === UniqueFusion()
+    return _artin_braid_single(f, i, inv)
+end
+# also for a tree that is the only one with its uncoupled sectors, whose braid is again a single
+# tree (`only` throws otherwise)
+function _artin_braid_single(f::FusionTree{I, N}, i, inv::Bool) where {I, N}
+    1 <= i < N || throw(ArgumentError(lazy"Cannot swap outputs i=$i and i+1 out of only $N outputs"))
 
     uncoupled = f.uncoupled
     a, b = uncoupled[i], uncoupled[i + 1]
@@ -53,7 +58,7 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N}
 
     if i == 1
         c = N > 2 ? inner[1] : coupled′
-        R = oftype(oneT, (inv ? conj(Rsymbol(b, a, c)) : Rsymbol(a, b, c)))
+        R = oftype(oneT, (inv ? conj(only(Rsymbol(b, a, c))) : only(Rsymbol(a, b, c))))
         f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner, vertices)
         return f′ => R
     end
@@ -64,13 +69,13 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N}
     a = inner_extended[i - 1]
     c = inner_extended[i]
     e = inner_extended[i + 1]
-    c′ = only(a ⊗ d)
+    c′ = FusionStyle(I) isa UniqueFusion ? only(a ⊗ d) : only(_artin_channels(a, b, d, e))
     coeff = oftype(
         oneT,
         if inv
-            conj(Rsymbol(d, c, e) * Fsymbol(d, a, b, e, c′, c)) * Rsymbol(d, a, c′)
+            conj(only(Rsymbol(d, c, e)) * only(Fsymbol(d, a, b, e, c′, c))) * only(Rsymbol(d, a, c′))
         else
-            Rsymbol(c, d, e) * conj(Fsymbol(d, a, b, e, c′, c) * Rsymbol(a, d, c′))
+            only(Rsymbol(c, d, e)) * conj(only(Fsymbol(d, a, b, e, c′, c)) * only(Rsymbol(a, d, c′)))
         end
     )
     inner′ = TupleTools.setindex(inner, c′, i - 1)
@@ -334,6 +339,21 @@ end
         p = linearizepermutation(p1, p2, numout(src), numin(src))
         levels = (l1..., reverse(l2)...)
 
+        if length(src) == 1 # then so is every intermediate block: the braid is a product of scalars
+            (f, f₀), coeff = _bend_single(_bendleft_single, only(fusiontrees(src)), Val(numin(src)))
+            for s in permutation2swaps(p)
+                _check_levels(BraidingStyle(I), levels, s)
+                inv = levels[s] > levels[s + 1]
+                f, c = _artin_braid_single(f, s, inv)
+                coeff *= c
+                l = levels[s]
+                levels = TupleTools.setindex(levels, levels[s + 1], s)
+                levels = TupleTools.setindex(levels, l, s + 1)
+            end
+            (f₁′, f₂′), c = _bend_single(_bendright_single, (f, f₀), Val(N₂))
+            return FusionTreeBlock([(f₁′, f₂′)]) => fill(coeff * c, 1, 1)
+        end
+
         # the intermediate blocks keep their trees in the order in which they are produced,
         # and only the result is brought into canonical order
         dst, U = repartition_noncanonical(src, Val(numind(src)))
diff --git a/src/fusiontrees/duality_manipulations.jl b/src/fusiontrees/duality_manipulations.jl
index 48cb106a2..78a7784c6 100644
--- a/src/fusiontrees/duality_manipulations.jl
+++ b/src/fusiontrees/duality_manipulations.jl
@@ -54,15 +54,18 @@ function _bendright_treepair((f₁, f₂)::FusionTreePair)
 end
 
 function bendright((f₁, f₂)::FusionTreePair)
-    I = sectortype((f₁, f₂))
+    @assert FusionStyle((f₁, f₂)) === UniqueFusion()
+    return _bendright_single((f₁, f₂))
+end
+# also for a final fusion vertex without multiplicity (`only` throws otherwise)
+function _bendright_single((f₁, f₂)::FusionTreePair)
     N₁ = numout((f₁, f₂))
-    @assert FusionStyle(I) === UniqueFusion()
     (a, b, c), (f₁′, f₂′) = _bendright_treepair((f₁, f₂))
 
     # compute the coefficient
     coeff₀ = sqrtdim(c) * invsqrtdim(a)
     f₁.isdual[N₁] && (coeff₀ *= conj(frobenius_schur_phase(dual(b))))
-    coeff = coeff₀ * Bsymbol(a, b, c)
+    coeff = coeff₀ * only(Bsymbol(a, b, c))
 
     return (f₁′, f₂′) => coeff
 end
@@ -143,6 +146,18 @@ function bendleft((f₁, f₂)::FusionTreePair)
     (f₂′, f₁′), coeff = bendright((f₂, f₁))
     return (f₁′, f₂′) => conj(coeff)
 end
+function _bendleft_single((f₁, f₂)::FusionTreePair)
+    (f₂′, f₁′), coeff = _bendright_single((f₂, f₁))
+    return (f₁′, f₂′) => conj(coeff)
+end
+
+# `n` bends of a tree pair that stays a single tree pair, and the product of their coefficients
+_bend_single(bend, f, ::Val{0}) = f => one(sectorscalartype(sectortype(f)))
+function _bend_single(bend, f, ::Val{n}) where {n}
+    f′, c′ = bend(f)
+    f″, c″ = _bend_single(bend, f′, Val(n - 1))
+    return f″ => c′ * c″
+end
 
 # !! note that this is more or less a copy of bendright through
 # (f1, f2) => conj(coeff) for ((f2, f1), coeff) in bendleft(src)

From b279d3a3ae77b4e5102b96fc5577b8d4e65e75ee Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Wed, 14 Oct 2026 14:47:21 +0200
Subject: [PATCH 04/13] Test braiding a block with a single tree; changelog

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 docs/src/Changelog.md         |  1 +
 test/symmetries/doubletree.jl | 25 +++++++++++++++++++++++++
 2 files changed, 26 insertions(+)

diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md
index aa3d1b297..9c90c2649 100644
--- a/docs/src/Changelog.md
+++ b/docs/src/Changelog.md
@@ -86,6 +86,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - Fusion tree block manipulations of (named) product sectors are computed as Kronecker products of the manipulations of their parts, splitting off the `UniqueFusion` sectors, such that they share the caches of the individual sectors ([#592](https://github.com/QuantumKitHub/TensorKit.jl/pull/592))
 - The TensorOperations ChainRules rules no longer copy or retain `C` when `β = Zero()`, as is the case for every tensor allocated by `@tensor` ([#555](https://github.com/QuantumKitHub/TensorKit.jl/pull/555))
 - The Enzyme `norm` rule reuses the cached primal and the `tr` rule avoids an unnecessary copy ([#554](https://github.com/QuantumKitHub/TensorKit.jl/pull/554))
+- Braiding fusion tree blocks no longer constructs and indexes the canonical intermediate blocks of every step: these keep their trees in the order in which they are produced, their buffers are recycled, and blocks with a single tree are braided as a product of scalars (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 0af5db19c..1c0f51652 100644
--- a/test/symmetries/doubletree.jl
+++ b/test/symmetries/doubletree.jl
@@ -235,6 +235,31 @@ using TensorKitSectors
         end
     end
 
+    FusionStyle(I) isa UniqueFusion || !(BraidingStyle(I) isa HasBraiding) ||
+        @testset "Double fusion tree: braiding a block with a single tree" begin
+        a, b = random_fusion(I, Val(2))
+        cs = filter(c -> Nsymbol(a, b, c) == 1, collect(a ⊗ b))
+        if !isempty(cs)
+            src1 = FusionTreeBlock{I}(((a, b), (rand(cs),)), ((rand(Bool), rand(Bool)), (rand(Bool),)))
+            @test length(src1) == 1
+            for n in 0:3
+                p = (randperm(3)...,)
+                p1, p2 = p[1:n], p[(n + 1):3]
+                ip = invperm(p)
+                ilevels = TupleTools.getindices((1, 2, 3), p)
+                dst, U = @constinferred TensorKit.braid(src1, (p1, p2), ((1, 2), (3,)))
+                @test length(dst) == 1 && size(U) == (1, 1)
+                dst′, U′ = @constinferred TensorKit.braid(dst, (ip[1:2], ip[3:3]), (ilevels[1:n], ilevels[(n + 1):3]))
+                @test dst′ == src1
+                @test _isone(U * U′)
+                if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
+                    A1 = fusiontensor(only(fusiontrees(src1)))
+                    @test permutedims(A1, (p1..., p2...)) ≈ only(U) * fusiontensor(only(fusiontrees(dst)))
+                end
+            end
+        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

From 698d70100f945449973f677eb50cc63c0475c8f8 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Tue, 20 Oct 2026 11:30:33 +0200
Subject: [PATCH 05/13] Assert the block type in the block branch of `fsbraid`
 for JET

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/braiding_manipulations.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index 4428e233e..1a1925852 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -335,6 +335,7 @@ end
         return (f₁′, f₂′) => coeff1 * coeff2 * coeff3
     else
         src, (p1, p2), (l1, l2) = key
+        src = src::FusionTreeBlock{I}
 
         p = linearizepermutation(p1, p2, numout(src), numin(src))
         levels = (l1..., reverse(l2)...)

From e824b1fe1c85c600384b6c1dd150988baf9e6e1f Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Fri, 23 Oct 2026 08:27:21 +0200
Subject: [PATCH 06/13] Fill a given matrix in `artin_braid!`, and choose the
 single braid channel by dispatch

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/braiding_manipulations.jl | 36 +++++++++++++++--------
 1 file changed, 23 insertions(+), 13 deletions(-)

diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index 1a1925852..94797d22a 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -69,7 +69,7 @@ function _artin_braid_single(f::FusionTree{I, N}, i, inv::Bool) where {I, N}
     a = inner_extended[i - 1]
     c = inner_extended[i]
     e = inner_extended[i + 1]
-    c′ = FusionStyle(I) isa UniqueFusion ? only(a ⊗ d) : only(_artin_channels(a, b, d, e))
+    c′ = _artin_channel(a, b, d, e)
     coeff = oftype(
         oneT,
         if inv
@@ -89,13 +89,24 @@ function _artin_channels(a::I, b::I, d::I, e::I) where {I <: Sector}
     UnitStyle(I) isa SimpleUnit || return intersect(a ⊗ d, e ⊗ conj(b))
     return Iterators.filter(c′ -> !iszero(Nsymbol(c′, b, e)), a ⊗ d)
 end
-
-function artin_braid(
-        src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false, canonical::Val = Val(true),
-        buffers = (nothing, nothing) # a tree vector and a matrix that may be overwritten
+# the only such line, for a fusion tree that braids into a single one
+_artin_channel(a, b, d, e) = _artin_channel(FusionStyle(a), a, b, d, e)
+_artin_channel(::UniqueFusion, a, b, d, e) = only(a ⊗ d)
+_artin_channel(::FusionStyle, a, b, d, e) = first(_artin_channels(a, b, d, e))
+
+function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where {I, N}
+    n = length(src) # braiding a complete block keeps its size
+    return artin_braid!(Matrix{sectorscalartype(I)}(undef, n, n), src, i; inv) # zeroed by `artin_braid!`
+end
+# `U` is overwritten with the braid matrix of `src`, a complete block; the trees of the result are
+# in canonical order or in the order in which they are produced, which recycles `trees` if given
+function artin_braid!(
+        U::AbstractMatrix, src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false,
+        canonical::Val = Val(true), trees = nothing
     ) where {I, N}
     1 <= i < N || throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs"))
-    isempty(src) && return src => zeros(sectorscalartype(I), 0, 0)
+    size(U) == (length(src), length(src)) || throw(DimensionMismatch(lazy"braid matrix of size $(size(U)) for a block of $(length(src)) trees"))
+    isempty(src) && return src => U
 
     uncoupled = src.uncoupled[1]
     a, b = uncoupled[i], uncoupled[i + 1]
@@ -106,12 +117,10 @@ function artin_braid(
     isdual = src.isdual[1]
     isdual′ = TupleTools.setindex(isdual, isdual[i], i + 1)
     isdual′ = TupleTools.setindex(isdual′, isdual[i + 1], i)
-    rows = TreeRows(I, (uncoupled′, ()), (isdual′, ()), length(src), canonical, buffers[1])
+    rows = TreeRows(I, (uncoupled′, ()), (isdual′, ()), length(src), canonical, trees)
 
     oneT = one(sectorscalartype(I))
-
-    U = buffers[2] isa Matrix{typeof(oneT)} ? fill!(buffers[2], zero(oneT)) :
-        zeros(sectorscalartype(I), max(length(rows.trees), length(src)), length(src)) # see `bendright`
+    fill!(U, zero(oneT))
 
     if isunit(a) || isunit(b) # braiding with trivial sector: simple and always possible
         for (col, (f, f₂)) in enumerate(fusiontrees(src))
@@ -361,14 +370,15 @@ end
 
         # the swaps keep the size of the block, so that the tree vector and the matrices of the
         # previous step are recycled (except for the trees of `src`, which are part of the key)
-        buffers = (nothing, nothing)
+        trees = nothing
+        U_tmp = Matrix{sectorscalartype(I)}(undef, length(dst), length(dst))
         U_spare = nothing
         braid_style = BraidingStyle(I)
         for s in permutation2swaps(p)
             _check_levels(braid_style, levels, s)
             inv = levels[s] > levels[s + 1]
-            dst′, U_tmp = artin_braid(dst, s; inv, canonical = Val(false), buffers)
-            buffers = (dst.trees === src.trees ? nothing : dst.trees, U_tmp)
+            dst′, _ = artin_braid!(U_tmp, dst, s; inv, canonical = Val(false), trees)
+            trees = dst.trees === src.trees ? nothing : dst.trees
             dst = dst′
             U_spare, U = U, U_spare isa typeof(U_tmp) ? mul!(U_spare, U_tmp, U) : U_tmp * U
             l = levels[s]

From 5a0905dcb80136f0ae903c5f62491217965e1af2 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Thu, 29 Oct 2026 06:51:13 +0100
Subject: [PATCH 07/13] Give `TreeRows` an inner constructor, so that Aqua
 finds no unbound type parameter

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/fusiontrees.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl
index 33faca9c4..b23598043 100644
--- a/src/fusiontrees/fusiontrees.jl
+++ b/src/fusiontrees/fusiontrees.jl
@@ -297,6 +297,8 @@ struct TreeRows{F, D}
     trees::Vector{F}
     index::Union{D, Nothing} # `nothing`: few trees, the rows are found by a linear search
     n::Int
+    # no default outer constructor, of which `D` would be unbound if `index === nothing`
+    TreeRows{F, D}(trees::Vector{F}, index, n::Int) where {F, D} = new{F, D}(trees, index, n)
 end
 function TreeRows(
         ::Type{I}, uncoupled, isdual, n::Int, ::Val{canonical}, buffer = nothing

From 6322e4dd34587c0d0413abeedaa6d6dd52dfa1b5 Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Thu, 29 Oct 2026 14:12:20 +0100
Subject: [PATCH 08/13] Braid single trees through `artin_braid` itself, take
 the only braid channel with `only`, and filter the braid channels for every
 unit style

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/braiding_manipulations.jl | 22 ++++++++--------------
 1 file changed, 8 insertions(+), 14 deletions(-)

diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index 94797d22a..437b17737 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -14,14 +14,11 @@ applying `artin_braid(f′, i; inv = true)` to all the outputs `f′` of
 `artin_braid(f, i; inv = false)` and collecting the results should yield a single fusion
 tree with non-zero coefficient, namely `f` with coefficient `1`. This keyword has no effect
 if `BraidingStyle(sectortype(f)) isa SymmetricBraiding`.
+
+For sectors without `UniqueFusion`, the braid of `f` has to be a single fusion tree, which is the
+case if `f` is the only fusion tree with its uncoupled and coupled sectors; this throws otherwise.
 """
 function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N}
-    @assert FusionStyle(I) === UniqueFusion()
-    return _artin_braid_single(f, i, inv)
-end
-# also for a tree that is the only one with its uncoupled sectors, whose braid is again a single
-# tree (`only` throws otherwise)
-function _artin_braid_single(f::FusionTree{I, N}, i, inv::Bool) where {I, N}
     1 <= i < N || throw(ArgumentError(lazy"Cannot swap outputs i=$i and i+1 out of only $N outputs"))
 
     uncoupled = f.uncoupled
@@ -84,15 +81,12 @@ function _artin_braid_single(f::FusionTree{I, N}, i, inv::Bool) where {I, N}
 end
 
 # the inner lines `c′` of `a ⊗ d → c′ → c′ ⊗ b → e`, i.e. `intersect(a ⊗ d, e ⊗ conj(b))`,
-# without constructing the second fusion product for sectors with a simple unit
-function _artin_channels(a::I, b::I, d::I, e::I) where {I <: Sector}
-    UnitStyle(I) isa SimpleUnit || return intersect(a ⊗ d, e ⊗ conj(b))
-    return Iterators.filter(c′ -> !iszero(Nsymbol(c′, b, e)), a ⊗ d)
-end
-# the only such line, for a fusion tree that braids into a single one
+# without constructing the second fusion product
+_artin_channels(a, b, d, e) = Iterators.filter(c′ -> !iszero(Nsymbol(c′, b, e)), a ⊗ d)
+# the only such line, for a fusion tree that braids into a single one (`only` throws otherwise)
 _artin_channel(a, b, d, e) = _artin_channel(FusionStyle(a), a, b, d, e)
 _artin_channel(::UniqueFusion, a, b, d, e) = only(a ⊗ d)
-_artin_channel(::FusionStyle, a, b, d, e) = first(_artin_channels(a, b, d, e))
+_artin_channel(::FusionStyle, a, b, d, e) = only(_artin_channels(a, b, d, e))
 
 function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where {I, N}
     n = length(src) # braiding a complete block keeps its size
@@ -354,7 +348,7 @@ end
             for s in permutation2swaps(p)
                 _check_levels(BraidingStyle(I), levels, s)
                 inv = levels[s] > levels[s + 1]
-                f, c = _artin_braid_single(f, s, inv)
+                f, c = artin_braid(f, s; inv)
                 coeff *= c
                 l = levels[s]
                 levels = TupleTools.setindex(levels, levels[s + 1], s)

From d13e2ee0fea6e03b0fdb55e37f702e7bbebf854d Mon Sep 17 00:00:00 2001
From: leburgel <lander.burgelman@gmail.com>
Date: Thu, 29 Oct 2026 14:12:20 +0100
Subject: [PATCH 09/13] Set the inner line and the first vertex of a braided
 tree outside the loops over the multiplicities

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/braiding_manipulations.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index 437b17737..28237a7bb 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -193,7 +193,9 @@ function artin_braid!(
             for c′ in _artin_channels(a, b, d, e)
                 Rmat2 = inv ? Rsymbol(d, a, c′)' : Rsymbol(a, d, c′)
                 Fmat = Fsymbol(d, a, b, e, c′, c)
+                inner′ = TupleTools.setindex(inner, c′, i - 1)
                 for σ in axes(Rmat2, 1)
+                    vertices_σ = TupleTools.setindex(vertices, σ, i - 1)
                     for λ in axes(Fmat, 2)
                         coeff = zero(oneT)
                         for ρ in axes(Rmat1, 2), κ in axes(Rmat2, 2)
@@ -201,9 +203,7 @@ function artin_braid!(
                                 conj(Rmat2[σ, κ])
                         end
                         iszero(coeff) && continue
-                        vertices′ = TupleTools.setindex(vertices, σ, i - 1)
-                        vertices′ = TupleTools.setindex(vertices′, λ, i)
-                        inner′ = TupleTools.setindex(inner, c′, i - 1)
+                        vertices′ = TupleTools.setindex(vertices_σ, λ, i)
                         f′ = FusionTree{I}(uncoupled′, coupled′, isdual′, inner′, vertices′)
                         row = row!(rows, (f′, f₂))
                         @inbounds U[row, col] = coeff

From 70246c45670b242de5181a2bfb7281db8075f473 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 29 Oct 2026 17:20:52 -0400
Subject: [PATCH 10/13] Braid a single tree pair through the tree manipulations
 themselves

`bendright` and `bendleft` of a tree pair no longer assert `UniqueFusion`, but
rely on `only` to throw if the bent vertex has a multiplicity, which replaces
their `_single` variants. The shortcut of `braid` for a symmetric braiding is
restricted to `UniqueFusion`, where the `Fsymbol`s are trivial, such that a
block with a single tree is braided as its tree pair.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/braiding_manipulations.jl | 42 ++++++++++-------------
 src/fusiontrees/duality_manipulations.jl  | 21 ++----------
 2 files changed, 20 insertions(+), 43 deletions(-)

diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index 28237a7bb..6bd973414 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -235,9 +235,8 @@ braid(f::FusionTree{I, N}, p::IndexTuple{N}, levels::IndexTuple{N}) where {I, N}
     braid(f, (p, ()), (levels, ()))
 function braid(f::FusionTree{I, N}, (p, _)::Index2Tuple{N, 0}, (levels, _)::Index2Tuple{N, 0}) where {I, N}
     TupleTools.isperm(p) || throw(ArgumentError(lazy"not a valid permutation: $p"))
-    @assert FusionStyle(I) isa UniqueFusion
     braid_style = BraidingStyle(I)
-    if braid_style isa SymmetricBraiding # this assumes Fsymbols are 1!
+    if FusionStyle(I) isa UniqueFusion && braid_style isa SymmetricBraiding # Fsymbols are 1
         coeff = one(sectorscalartype(I))
         for i in 1:N
             for j in 1:(i - 1)
@@ -329,35 +328,20 @@ end
 
 @cached function fsbraid(key::K)::_fsdicttype(K) where {I, N₁, N₂, K <: Union{FSPBraidKey{I, N₁, N₂}, FSBBraidKey{I, N₁, N₂}}}
     if K <: FSPBraidKey
-        ((f₁, f₂), (p1, p2), (l1, l2)) = key
-        p = linearizepermutation(p1, p2, length(f₁), length(f₂))
-        levels = (l1..., reverse(l2)...)
-        (f, f0), coeff1 = repartition((f₁, f₂), N₁ + N₂)
-        f′, coeff2 = braid(f, p, levels)
-        (f₁′, f₂′), coeff3 = repartition((f′, f0), N₁)
-        return (f₁′, f₂′) => coeff1 * coeff2 * coeff3
+        f, p, levels = key
+        return _braid_single(f::FusionTreePair{I}, p, levels)
     else
         src, (p1, p2), (l1, l2) = key
         src = src::FusionTreeBlock{I}
 
-        p = linearizepermutation(p1, p2, numout(src), numin(src))
-        levels = (l1..., reverse(l2)...)
-
         if length(src) == 1 # then so is every intermediate block: the braid is a product of scalars
-            (f, f₀), coeff = _bend_single(_bendleft_single, only(fusiontrees(src)), Val(numin(src)))
-            for s in permutation2swaps(p)
-                _check_levels(BraidingStyle(I), levels, s)
-                inv = levels[s] > levels[s + 1]
-                f, c = artin_braid(f, s; inv)
-                coeff *= c
-                l = levels[s]
-                levels = TupleTools.setindex(levels, levels[s + 1], s)
-                levels = TupleTools.setindex(levels, l, s + 1)
-            end
-            (f₁′, f₂′), c = _bend_single(_bendright_single, (f, f₀), Val(N₂))
-            return FusionTreeBlock([(f₁′, f₂′)]) => fill(coeff * c, 1, 1)
+            f′, coeff = _braid_single(only(fusiontrees(src)), (p1, p2), (l1, l2))
+            return FusionTreeBlock([f′]) => fill(coeff, 1, 1)
         end
 
+        p = linearizepermutation(p1, p2, numout(src), numin(src))
+        levels = (l1..., reverse(l2)...)
+
         # the intermediate blocks keep their trees in the order in which they are produced,
         # and only the result is brought into canonical order
         dst, U = repartition_noncanonical(src, Val(numind(src)))
@@ -390,6 +374,16 @@ end
     end
 end
 
+# also for a tree pair that is the only one with its uncoupled sectors (`only` throws otherwise)
+function _braid_single((f₁, f₂)::FusionTreePair, (p1, p2)::Index2Tuple{N₁, N₂}, (l1, l2)::Index2Tuple) where {N₁, N₂}
+    p = linearizepermutation(p1, p2, length(f₁), length(f₂))
+    levels = (l1..., reverse(l2)...)
+    (f, f0), coeff1 = repartition((f₁, f₂), N₁ + N₂)
+    f′, coeff2 = braid(f, p, levels)
+    (f₁′, f₂′), coeff3 = repartition((f′, f0), N₁)
+    return (f₁′, f₂′) => coeff1 * coeff2 * coeff3
+end
+
 CacheStyle(::typeof(fsbraid), k::Union{FSPBraidKey{I}, FSBBraidKey{I}}) where {I} =
     FusionStyle(I) isa UniqueFusion ? NoCache() : GlobalLRUCache()
 
diff --git a/src/fusiontrees/duality_manipulations.jl b/src/fusiontrees/duality_manipulations.jl
index 78a7784c6..bec996e76 100644
--- a/src/fusiontrees/duality_manipulations.jl
+++ b/src/fusiontrees/duality_manipulations.jl
@@ -53,12 +53,8 @@ function _bendright_treepair((f₁, f₂)::FusionTreePair)
     return (a, b, c), (f₁′, f₂′)
 end
 
+# also for a final splitting vertex without multiplicity (`only` throws otherwise)
 function bendright((f₁, f₂)::FusionTreePair)
-    @assert FusionStyle((f₁, f₂)) === UniqueFusion()
-    return _bendright_single((f₁, f₂))
-end
-# also for a final fusion vertex without multiplicity (`only` throws otherwise)
-function _bendright_single((f₁, f₂)::FusionTreePair)
     N₁ = numout((f₁, f₂))
     (a, b, c), (f₁′, f₂′) = _bendright_treepair((f₁, f₂))
 
@@ -142,22 +138,9 @@ See also [`bendright`](@ref).
 """ bendleft
 
 function bendleft((f₁, f₂)::FusionTreePair)
-    @assert FusionStyle((f₁, f₂)) === UniqueFusion()
     (f₂′, f₁′), coeff = bendright((f₂, f₁))
     return (f₁′, f₂′) => conj(coeff)
 end
-function _bendleft_single((f₁, f₂)::FusionTreePair)
-    (f₂′, f₁′), coeff = _bendright_single((f₂, f₁))
-    return (f₁′, f₂′) => conj(coeff)
-end
-
-# `n` bends of a tree pair that stays a single tree pair, and the product of their coefficients
-_bend_single(bend, f, ::Val{0}) = f => one(sectorscalartype(sectortype(f)))
-function _bend_single(bend, f, ::Val{n}) where {n}
-    f′, c′ = bend(f)
-    f″, c″ = _bend_single(bend, f′, Val(n - 1))
-    return f″ => c′ * c″
-end
 
 # !! note that this is more or less a copy of bendright through
 # (f1, f2) => conj(coeff) for ((f2, f1), coeff) in bendleft(src)
@@ -499,7 +482,7 @@ function _repartition_body(N, canonical = true)
     if N == 0
         ex = quote
             T = fusionscalartype(sectortype(src))
-            if FusionStyle(src) === UniqueFusion()
+            if src isa FusionTreePair
                 return src => one(T)
             else
                 U = copyto!(zeros(T, length(src), length(src)), LinearAlgebra.I)

From 520e55bb35bf292f994d0faf8a94040d3a07f404 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 29 Oct 2026 16:08:53 -0400
Subject: [PATCH 11/13] Keep the trees of every fusion tree block manipulation
 in production order

The destination blocks of bends, folds and braids are only used together with
their transformation matrix, so they no longer construct and index the canonical
destination block, but keep their trees in the order in which they are produced.
Since `planar_trace` now also works on complete blocks, the canonical mode of
`TreeRows`, `canonicalize`, `repartition_noncanonical` and `treeindex_map` are
removed. The tests compare the results of manipulations up to the order of
their trees.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 docs/src/Changelog.md                     |  2 +-
 src/fusiontrees/braiding_manipulations.jl | 23 +++----
 src/fusiontrees/duality_manipulations.jl  | 78 ++++++-----------------
 src/fusiontrees/fusiontrees.jl            | 51 ++++-----------
 test/setup.jl                             | 13 +++-
 test/symmetries/doubletree.jl             | 20 +++---
 6 files changed, 61 insertions(+), 126 deletions(-)

diff --git a/docs/src/Changelog.md b/docs/src/Changelog.md
index 9c90c2649..8f88c8ad6 100644
--- a/docs/src/Changelog.md
+++ b/docs/src/Changelog.md
@@ -86,7 +86,7 @@ When releasing a new version, move the "Unreleased" changes to a new version sec
 - Fusion tree block manipulations of (named) product sectors are computed as Kronecker products of the manipulations of their parts, splitting off the `UniqueFusion` sectors, such that they share the caches of the individual sectors ([#592](https://github.com/QuantumKitHub/TensorKit.jl/pull/592))
 - The TensorOperations ChainRules rules no longer copy or retain `C` when `β = Zero()`, as is the case for every tensor allocated by `@tensor` ([#555](https://github.com/QuantumKitHub/TensorKit.jl/pull/555))
 - The Enzyme `norm` rule reuses the cached primal and the `tr` rule avoids an unnecessary copy ([#554](https://github.com/QuantumKitHub/TensorKit.jl/pull/554))
-- Braiding fusion tree blocks no longer constructs and indexes the canonical intermediate blocks of every step: these keep their trees in the order in which they are produced, their buffers are recycled, and blocks with a single tree are braided as a product of scalars (internal)
+- Fusion tree block manipulations no longer construct and index the canonical destination block, but keep the trees of the result in the order in which they are produced; the buffers of the intermediate blocks of a braid are recycled, and blocks with a single tree are braided as their tree pair (internal)
 
 ## [0.17.2](https://github.com/QuantumKitHub/TensorKit.jl/compare/v0.17.1...v0.17.2) - 2026-09-20
 
diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index 6bd973414..60f968c29 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -93,10 +93,9 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
     return artin_braid!(Matrix{sectorscalartype(I)}(undef, n, n), src, i; inv) # zeroed by `artin_braid!`
 end
 # `U` is overwritten with the braid matrix of `src`, a complete block; the trees of the result are
-# in canonical order or in the order in which they are produced, which recycles `trees` if given
+# stored in `trees` if given
 function artin_braid!(
-        U::AbstractMatrix, src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false,
-        canonical::Val = Val(true), trees = nothing
+        U::AbstractMatrix, src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false, trees = nothing
     ) where {I, N}
     1 <= i < N || throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs"))
     size(U) == (length(src), length(src)) || throw(DimensionMismatch(lazy"braid matrix of size $(size(U)) for a block of $(length(src)) trees"))
@@ -111,7 +110,7 @@ function artin_braid!(
     isdual = src.isdual[1]
     isdual′ = TupleTools.setindex(isdual, isdual[i], i + 1)
     isdual′ = TupleTools.setindex(isdual′, isdual[i + 1], i)
-    rows = TreeRows(I, (uncoupled′, ()), (isdual′, ()), length(src), canonical, trees)
+    rows = TreeRows(fusiontreetype(I, N, 0), length(src), trees)
 
     oneT = one(sectorscalartype(I))
     fill!(U, zero(oneT))
@@ -342,9 +341,7 @@ end
         p = linearizepermutation(p1, p2, numout(src), numin(src))
         levels = (l1..., reverse(l2)...)
 
-        # the intermediate blocks keep their trees in the order in which they are produced,
-        # and only the result is brought into canonical order
-        dst, U = repartition_noncanonical(src, Val(numind(src)))
+        dst, U = repartition(src, Val(numind(src)))
 
         # the swaps keep the size of the block, so that the tree vector and the matrices of the
         # previous step are recycled (except for the trees of `src`, which are part of the key)
@@ -355,7 +352,7 @@ end
         for s in permutation2swaps(p)
             _check_levels(braid_style, levels, s)
             inv = levels[s] > levels[s + 1]
-            dst′, _ = artin_braid!(U_tmp, dst, s; inv, canonical = Val(false), trees)
+            dst′, _ = artin_braid!(U_tmp, dst, s; inv, trees)
             trees = dst.trees === src.trees ? nothing : dst.trees
             dst = dst′
             U_spare, U = U, U_spare isa typeof(U_tmp) ? mul!(U_spare, U_tmp, U) : U_tmp * U
@@ -364,13 +361,9 @@ end
             levels = TupleTools.setindex(levels, l, s + 1)
         end
 
-        if N₂ == 0
-            return canonicalize(dst, U)
-        else
-            dst, U_tmp = repartition_noncanonical(dst::FusionTreeBlock{I}, Val(N₁))
-            U = U_tmp * U
-            return canonicalize(dst, U)
-        end
+        N₂ == 0 && return dst => U
+        dst, U_tmp = repartition(dst::FusionTreeBlock{I}, Val(N₁))
+        return dst => U_tmp * U
     end
 end
 
diff --git a/src/fusiontrees/duality_manipulations.jl b/src/fusiontrees/duality_manipulations.jl
index bec996e76..171ef2ba7 100644
--- a/src/fusiontrees/duality_manipulations.jl
+++ b/src/fusiontrees/duality_manipulations.jl
@@ -65,24 +65,14 @@ function bendright((f₁, f₂)::FusionTreePair)
 
     return (f₁′, f₂′) => coeff
 end
-function bendright(src::FusionTreeBlock, canonical::Val = Val(true))
+function bendright(src::FusionTreeBlock)
     I = sectortype(src)
     N₁ = numout(src)
     N₂ = numin(src)
     @assert N₁ > 0
-    uncoupled_dst = (
-        TupleTools.front(src.uncoupled[1]),
-        (src.uncoupled[2]..., dual(src.uncoupled[1][N₁])),
-    )
-    isdual_dst = (
-        TupleTools.front(src.isdual[1]),
-        (src.isdual[2]..., !(src.isdual[1][N₁])),
-    )
-
-    rows = TreeRows(I, uncoupled_dst, isdual_dst, length(src), canonical)
-    # the canonical destination block has more trees than `src` if `src` holds only some of the
-    # trees with its uncoupled sectors (as in `planar_trace`)
-    U = zeros(fusionscalartype(I), max(length(rows.trees), length(src)), length(src))
+
+    rows = TreeRows(fusiontreetype(I, N₁ - 1, N₂ + 1), length(src))
+    U = zeros(fusionscalartype(I), length(src), length(src))
 
     for (col, (f₁, f₂)) in enumerate(fusiontrees(src))
         (a, b, c), (f₁′, f₂′) = _bendright_treepair((f₁, f₂))
@@ -144,24 +134,14 @@ end
 
 # !! note that this is more or less a copy of bendright through
 # (f1, f2) => conj(coeff) for ((f2, f1), coeff) in bendleft(src)
-function bendleft(src::FusionTreeBlock, canonical::Val = Val(true))
+function bendleft(src::FusionTreeBlock)
     I = sectortype(src)
     N₁ = numout(src)
     N₂ = numin(src)
     @assert N₂ > 0
-    uncoupled_dst = (
-        (src.uncoupled[1]..., dual(src.uncoupled[2][N₂])),
-        TupleTools.front(src.uncoupled[2]),
-    )
-    isdual_dst = (
-        (src.isdual[1]..., !(src.isdual[2][N₂])),
-        TupleTools.front(src.isdual[2]),
-    )
-
-    rows = TreeRows(I, uncoupled_dst, isdual_dst, length(src), canonical)
-    # the canonical destination block has more trees than `src` if `src` holds only some of the
-    # trees with its uncoupled sectors (as in `planar_trace`)
-    U = zeros(fusionscalartype(I), max(length(rows.trees), length(src)), length(src))
+
+    rows = TreeRows(fusiontreetype(I, N₁ + 1, N₂ - 1), length(src))
+    U = zeros(fusionscalartype(I), length(src), length(src))
 
     for (col, (f₁, f₂)) in enumerate(fusiontrees(src))
         (a, b, c), (f₂′, f₁′) = _bendright_treepair((f₂, f₁))
@@ -236,17 +216,11 @@ function foldright((f₁, f₂)::FusionTreePair)
 end
 
 function foldright(src::FusionTreeBlock)
-    uncoupled_dst = (
-        Base.tail(src.uncoupled[1]),
-        (dual(first(src.uncoupled[1])), src.uncoupled[2]...),
-    )
-    isdual_dst = (Base.tail(src.isdual[1]), (!first(src.isdual[1]), src.isdual[2]...))
     I = sectortype(src)
     N₁ = numout(src)
     @assert N₁ > 0
 
-    dst = FusionTreeBlock{I}(uncoupled_dst, isdual_dst; sizehint = length(src))
-    indexmap = treeindex_map(dst)
+    rows = TreeRows(fusiontreetype(I, N₁ - 1, numin(src) + 1), length(src))
 
     f₁, f₂ = first(fusiontrees(src))
     a::I = f₁.uncoupled[1]
@@ -260,7 +234,7 @@ function foldright(src::FusionTreeBlock)
     c::I = f₁.coupled
     cache₃ = Dict((b, c) => Asymbol(a, b, c))
 
-    U = zeros(eltype(coeff₁), length(dst), length(src))
+    U = zeros(eltype(coeff₁), length(src), length(src))
     for (col, (f₁, f₂)) in enumerate(fusiontrees(src))
         f₁′s, coeffs₁ = get!(cache₁, f₁) do
             multi_Fmove(f₁)
@@ -280,12 +254,12 @@ function foldright(src::FusionTreeBlock)
                 if isduala
                     coeff *= κₐ
                 end
-                row = indexmap[treeindex_data((f₁′, f₂′))]
+                row = row!(rows, (f₁′, f₂′))
                 @inbounds U[row, col] += coeff
             end
         end
     end
-    return dst => U
+    return FusionTreeBlock(rows.trees) => U
 end
 
 
@@ -319,21 +293,12 @@ function foldleft((f₁, f₂)::FusionTreePair)
 end
 
 function foldleft(src::FusionTreeBlock)
-    uncoupled_dst = (
-        (dual(first(src.uncoupled[2])), src.uncoupled[1]...),
-        Base.tail(src.uncoupled[2]),
-    )
-    isdual_dst = (
-        (!first(src.isdual[2]), src.isdual[1]...),
-        Base.tail(src.isdual[2]),
-    )
     I = sectortype(src)
     N₁ = numin(src)
     N₂ = numout(src)
     @assert N₁ > 0
 
-    dst = FusionTreeBlock{I}(uncoupled_dst, isdual_dst; sizehint = length(src))
-    indexmap = treeindex_map(dst)
+    rows = TreeRows(fusiontreetype(I, N₂ + 1, N₁ - 1), length(src))
 
     f₁, f₂ = first(fusiontrees(src))
     a::I = f₂.uncoupled[1]
@@ -347,7 +312,7 @@ function foldleft(src::FusionTreeBlock)
     c::I = f₂.coupled
     cache₃ = Dict((b, c) => Asymbol(a, b, c))
 
-    U = zeros(eltype(coeff₂), length(dst), length(src))
+    U = zeros(eltype(coeff₂), length(src), length(src))
     for (col, (f₁, f₂)) in enumerate(fusiontrees(src))
         f₂′s, coeffs₂ = get!(cache₂, f₂) do
             multi_Fmove(f₂)
@@ -367,12 +332,12 @@ function foldleft(src::FusionTreeBlock)
                 if isduala
                     coeff *= conj(κₐ)
                 end
-                row = indexmap[treeindex_data((f₁′, f₂′))]
+                row = row!(rows, (f₁′, f₂′))
                 @inbounds U[row, col] += coeff
             end
         end
     end
-    return dst => U
+    return FusionTreeBlock(rows.trees) => U
 end
 
 # clockwise cyclic permutation while preserving (N₁, N₂): foldright & bendleft
@@ -478,7 +443,7 @@ return dst, U
 # it, so that it is visible to the generator's world age. Otherwise a precompile workload that
 # first triggers `repartition` fails with `UndefVarError: _repartition_body` (at runtime the
 # first call happens at a later world age, which masks the ordering issue).
-function _repartition_body(N, canonical = true)
+function _repartition_body(N)
     if N == 0
         ex = quote
             T = fusionscalartype(sectortype(src))
@@ -491,14 +456,13 @@ function _repartition_body(N, canonical = true)
         end
     else
         f = N < 0 ? bendleft : bendright
-        args = canonical ? () : (Val(false),)
         ex_rep = Expr(:block)
         for _ in 1:(abs(N) - 1)
-            push!(ex_rep.args, :((dst, Utmp) = $f(dst, $(args...))))
+            push!(ex_rep.args, :((dst, Utmp) = $f(dst)))
             push!(ex_rep.args, :(U = Utmp * U))
         end
         ex = quote
-            dst, U = $f(src, $(args...))
+            dst, U = $f(src)
             $ex_rep
             return dst => U
         end
@@ -508,10 +472,6 @@ end
 @generated function repartition(src::Union{FusionTreePair, FusionTreeBlock}, ::Val{N}) where {N}
     return _repartition_body(numout(src) - N)
 end
-# the same for a block, with the trees of the result in the order in which they are produced
-@generated function repartition_noncanonical(src::FusionTreeBlock, ::Val{N}) where {N}
-    return _repartition_body(numout(src) - N, false)
-end
 
 """
     transpose((f₁, f₂)::FusionTreePair{I}, p::Index2Tuple{N₁, N₂}) where {I, N₁, N₂}
diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl
index b23598043..ace229b7c 100644
--- a/src/fusiontrees/fusiontrees.jl
+++ b/src/fusiontrees/fusiontrees.jl
@@ -117,8 +117,10 @@ const FusionTreePair{I, N₁, N₂} = Tuple{FusionTree{I, N₁}, FusionTree{I, N
 
 Collection of fusion-splitting tree pairs that share the same uncoupled sectors.
 Mostly internal structure to speed up non-`UniqueFusion` fusiontree manipulation computations.
-Braiding a block requires it to be complete, i.e. to contain all the fusion-splitting tree pairs
-with its uncoupled sectors and duality flags, as `FusionTreeBlock{I}(uncoupled, isdual)` does.
+Manipulating a block requires it to be complete, i.e. to contain all the fusion-splitting tree
+pairs with its uncoupled sectors and duality flags, as `FusionTreeBlock{I}(uncoupled, isdual)`
+does. The resulting block holds its trees in the order in which they are produced, rather than in
+the canonical order of `FusionTreeBlock{I}(uncoupled, isdual)`.
 """
 struct FusionTreeBlock{I, N₁, N₂, F <: FusionTreePair{I, N₁, N₂}}
     trees::Vector{F}
@@ -280,19 +282,11 @@ function treeindex_data((f₁, f₂))
         return false
     end
 end
-function treeindex_map(fs::FusionTreeBlock)
-    I = sectortype(fs)
-    return fusiontreedict(I)(treeindex_data(f) => ind for (ind, f) in enumerate(fusiontrees(fs)))
-end
 
-# The trees of the destination block of a block manipulation, and the row of each of them in the
-# transformation matrix. With `Val(true)`, these are the trees of the destination block in their
-# canonical order; with `Val(false)`, the trees are added in the order in which they are produced,
-# which avoids constructing the destination block. The latter is used for the intermediate blocks
-# of a composite manipulation, of which only the result is brought into canonical order; the tree
-# vector of an intermediate block that is no longer needed can then be passed on as `buffer`.
-# With `Val(false)`, the source block has to be complete: the destination block then has as many
-# trees (`n`), and an incomplete source block raises an error if its destination block has more.
+# The trees of the destination block of a block manipulation, in the order in which they are
+# produced, and the row of each of them in the transformation matrix. The source block has to be
+# complete, such that the destination block has as many trees (`n`); the tree vector of a block
+# that is no longer needed can be passed on as `buffer`.
 struct TreeRows{F, D}
     trees::Vector{F}
     index::Union{D, Nothing} # `nothing`: few trees, the rows are found by a linear search
@@ -300,18 +294,11 @@ struct TreeRows{F, D}
     # no default outer constructor, of which `D` would be unbound if `index === nothing`
     TreeRows{F, D}(trees::Vector{F}, index, n::Int) where {F, D} = new{F, D}(trees, index, n)
 end
-function TreeRows(
-        ::Type{I}, uncoupled, isdual, n::Int, ::Val{canonical}, buffer = nothing
-    ) where {I, canonical}
-    F = fusiontreetype(I, length(uncoupled[1]), length(uncoupled[2]))
-    trees = if canonical
-        fusiontrees(FusionTreeBlock{I}(uncoupled, isdual; sizehint = n))
-    else
-        buffer isa Vector{F} ? empty!(buffer) : sizehint!(Vector{F}(), n)
-    end
-    K = Base.promote_op(treeindex_data, eltype(trees))
-    index = n > 32 ? Dict{K, Int}(treeindex_data(f) => i for (i, f) in enumerate(trees)) : nothing
-    return TreeRows{eltype(trees), Dict{K, Int}}(trees, index, n)
+function TreeRows(::Type{F}, n::Int, buffer = nothing) where {F}
+    trees = buffer isa Vector{F} ? empty!(buffer) : sizehint!(Vector{F}(), n)
+    K = Base.promote_op(treeindex_data, F)
+    index = n > 32 ? Dict{K, Int}() : nothing
+    return TreeRows{F, Dict{K, Int}}(trees, index, n)
 end
 function row!(rows::TreeRows, f)
     key = treeindex_data(f)
@@ -322,22 +309,12 @@ function row!(rows::TreeRows, f)
     end
     isnothing(i) || return i
     length(rows.trees) < rows.n ||
-        throw(ArgumentError("braiding requires a complete `FusionTreeBlock`"))
+        throw(ArgumentError("manipulating a block requires it to be complete"))
     push!(rows.trees, f)
     isnothing(rows.index) || (rows.index[key] = length(rows.trees))
     return length(rows.trees)
 end
 
-# bring a block whose trees are not in canonical order, together with the transformation matrix
-# whose rows correspond to these trees, into canonical order
-function canonicalize(src::FusionTreeBlock{I}, U::AbstractMatrix) where {I}
-    length(src) == 1 && return src => U
-    rows = TreeRows(I, src.uncoupled, src.isdual, length(src), Val(true))
-    U′ = zeros(eltype(U), length(rows.trees), size(U, 2))
-    U′[map(f -> row!(rows, f), fusiontrees(src)), :] = U
-    return FusionTreeBlock(rows.trees) => U′
-end
-
 
 # Facilitate getting correct fusion tree types
 Base.@assume_effects :foldable function fusiontreetype(::Type{I}, N::Int) where {I <: Sector}
diff --git a/test/setup.jl b/test/setup.jl
index dca3a6c4e..d44366ae9 100644
--- a/test/setup.jl
+++ b/test/setup.jl
@@ -9,12 +9,12 @@ export sectorlist, fast_sectorlist
 export default_spacelist, factorization_spacelist, ad_spacelist
 export VIBM, VIBMRepA4
 export test_ad_rrule
-export _isunitary, _isone
+export _isunitary, _isone, _inorder
 
 using Random
 using Test: @test
 using TensorKit
-using TensorKit: ℙ, PlanarTrivial
+using TensorKit: ℙ, PlanarTrivial, FusionTreeBlock
 using TensorKitSectors
 using TensorOperations: IndexTuple, Index2Tuple
 using Base.Iterators: take, product
@@ -149,6 +149,15 @@ _isunitary(x::Number; kwargs...) = isapprox(x * x', one(x); kwargs...)
 _isunitary(x; kwargs...) = isunitary(x; kwargs...)
 _isone(x; kwargs...) = isapprox(x, one(x); kwargs...)
 
+# the rows of the transformation matrix `U` onto the trees of `dst`, in the order of the trees of
+# `ref`, as block manipulations keep their trees in the order in which they are produced
+_inorder(U::Number, dst, ref) = (@test dst == ref; U)
+function _inorder(U::AbstractMatrix, dst::FusionTreeBlock, ref::FusionTreeBlock)
+    rows = indexin(fusiontrees(ref), fusiontrees(dst))
+    @test length(dst) == length(ref) && !any(isnothing, rows)
+    return U[something.(rows, 1), :]
+end
+
 
 uniquefusionsectorlist = (
     Z2Irrep, Z3Irrep, Z4Irrep, Z3Irrep ⊠ Z4Irrep, U1Irrep,
diff --git a/test/symmetries/doubletree.jl b/test/symmetries/doubletree.jl
index 1c0f51652..1d51b6730 100644
--- a/test/symmetries/doubletree.jl
+++ b/test/symmetries/doubletree.jl
@@ -57,15 +57,13 @@ using TensorKitSectors
         # single bend
         dst, U = @constinferred TK.bendright(src)
         dst2, U2 = @constinferred TK.bendleft(dst)
-        @test src == dst2
-        @test _isone(U2 * U)
+        @test _isone(_inorder(U2 * U, dst2, src))
         # double bend
         dst1, U1 = @constinferred TK.bendleft(src)
         dst2, U2 = @constinferred TK.bendleft(dst1)
         dst3, U3 = @constinferred TK.bendright(dst2)
         dst4, U4 = @constinferred TK.bendright(dst3)
-        @test src == dst4
-        @test _isone(U4 * U3 * U2 * U1)
+        @test _isone(_inorder(U4 * U3 * U2 * U1, dst4, src))
 
         if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
             all_inds = (ntuple(identity, numout(src))..., reverse(ntuple(i -> i + numout(src), numin(src)))...)
@@ -88,15 +86,13 @@ using TensorKitSectors
         # single bend
         dst, U = @constinferred TK.foldleft(src)
         dst2, U2 = @constinferred TK.foldright(dst)
-        @test src == dst2
-        @test _isone(U2 * U)
+        @test _isone(_inorder(U2 * U, dst2, src))
         # double bend
         dst1, U1 = @constinferred TK.foldright(src)
         dst2, U2 = @constinferred TK.foldright(dst1)
         dst3, U3 = @constinferred TK.foldleft(dst2)
         dst4, U4 = @constinferred TK.foldleft(dst3)
-        @test src == dst4
-        @test _isone(U4 * U3 * U2 * U1)
+        @test _isone(_inorder(U4 * U3 * U2 * U1, dst4, src))
 
         if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
             all_inds = TupleTools.circshift((ntuple(identity, numout(src))..., reverse(ntuple(i -> i + numout(src), numin(src)))...), -2)
@@ -121,7 +117,7 @@ using TensorKitSectors
             # @test _isunitary(U)
 
             dst′, U′ = repartition(dst, N)
-            @test _isone(U * U′)
+            @test _isone(U * _inorder(U′, dst′, src))
 
             if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
                 all_inds = (ntuple(identity, numout(src))..., reverse(ntuple(i -> i + numout(src), numin(src)))...)
@@ -152,11 +148,11 @@ using TensorKitSectors
 
             dst, U = @constinferred transpose(src, (p1, p2))
             dst′, U′ = @constinferred transpose(dst, (ip1, ip2))
-            @test _isone(U * U′)
+            @test _isone(U * _inorder(U′, dst′, src))
 
             if BraidingStyle(I) isa Bosonic
                 dst″, U″ = permute(src, (p1, p2))
-                @test U″ ≈ U
+                @test _inorder(U″, dst″, dst) ≈ U
             end
 
             if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)
@@ -217,7 +213,7 @@ using TensorKitSectors
             else
                 dst′, U′ = @constinferred TensorKit.braid(dst, (ip1, ip2), (il1, il2))
             end
-            @test _isone(U * U′)
+            @test _isone(U * _inorder(U′, dst′, src))
 
             # check fusiontensor compatibility
             if (BraidingStyle(I) isa Bosonic) && hasfusiontensor(I)

From 238ee049d84b67b6ad17e42a40356fee062a90a5 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Wed, 4 Nov 2026 18:23:42 -0500
Subject: [PATCH 12/13] Drop the recycling of the tree vectors in `fsbraid`,
 which gives no measurable speedup

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/braiding_manipulations.jl | 17 +++++------------
 src/fusiontrees/fusiontrees.jl            |  7 +++----
 2 files changed, 8 insertions(+), 16 deletions(-)

diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index 60f968c29..a621c43f0 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -92,11 +92,8 @@ function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where
     n = length(src) # braiding a complete block keeps its size
     return artin_braid!(Matrix{sectorscalartype(I)}(undef, n, n), src, i; inv) # zeroed by `artin_braid!`
 end
-# `U` is overwritten with the braid matrix of `src`, a complete block; the trees of the result are
-# stored in `trees` if given
-function artin_braid!(
-        U::AbstractMatrix, src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false, trees = nothing
-    ) where {I, N}
+# `U` is overwritten with the braid matrix of `src`, a complete block
+function artin_braid!(U::AbstractMatrix, src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where {I, N}
     1 <= i < N || throw(ArgumentError("Cannot swap outputs i=$i and i+1 out of only $N outputs"))
     size(U) == (length(src), length(src)) || throw(DimensionMismatch(lazy"braid matrix of size $(size(U)) for a block of $(length(src)) trees"))
     isempty(src) && return src => U
@@ -110,7 +107,7 @@ function artin_braid!(
     isdual = src.isdual[1]
     isdual′ = TupleTools.setindex(isdual, isdual[i], i + 1)
     isdual′ = TupleTools.setindex(isdual′, isdual[i + 1], i)
-    rows = TreeRows(fusiontreetype(I, N, 0), length(src), trees)
+    rows = TreeRows(fusiontreetype(I, N, 0), length(src))
 
     oneT = one(sectorscalartype(I))
     fill!(U, zero(oneT))
@@ -343,18 +340,14 @@ end
 
         dst, U = repartition(src, Val(numind(src)))
 
-        # the swaps keep the size of the block, so that the tree vector and the matrices of the
-        # previous step are recycled (except for the trees of `src`, which are part of the key)
-        trees = nothing
+        # the swaps keep the size of the block, so that the matrices of the previous step are recycled
         U_tmp = Matrix{sectorscalartype(I)}(undef, length(dst), length(dst))
         U_spare = nothing
         braid_style = BraidingStyle(I)
         for s in permutation2swaps(p)
             _check_levels(braid_style, levels, s)
             inv = levels[s] > levels[s + 1]
-            dst′, _ = artin_braid!(U_tmp, dst, s; inv, trees)
-            trees = dst.trees === src.trees ? nothing : dst.trees
-            dst = dst′
+            dst, _ = artin_braid!(U_tmp, dst, s; inv)
             U_spare, U = U, U_spare isa typeof(U_tmp) ? mul!(U_spare, U_tmp, U) : U_tmp * U
             l = levels[s]
             levels = TupleTools.setindex(levels, levels[s + 1], s)
diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl
index ace229b7c..32695e210 100644
--- a/src/fusiontrees/fusiontrees.jl
+++ b/src/fusiontrees/fusiontrees.jl
@@ -285,8 +285,7 @@ end
 
 # The trees of the destination block of a block manipulation, in the order in which they are
 # produced, and the row of each of them in the transformation matrix. The source block has to be
-# complete, such that the destination block has as many trees (`n`); the tree vector of a block
-# that is no longer needed can be passed on as `buffer`.
+# complete, such that the destination block has as many trees (`n`).
 struct TreeRows{F, D}
     trees::Vector{F}
     index::Union{D, Nothing} # `nothing`: few trees, the rows are found by a linear search
@@ -294,8 +293,8 @@ struct TreeRows{F, D}
     # no default outer constructor, of which `D` would be unbound if `index === nothing`
     TreeRows{F, D}(trees::Vector{F}, index, n::Int) where {F, D} = new{F, D}(trees, index, n)
 end
-function TreeRows(::Type{F}, n::Int, buffer = nothing) where {F}
-    trees = buffer isa Vector{F} ? empty!(buffer) : sizehint!(Vector{F}(), n)
+function TreeRows(::Type{F}, n::Int) where {F}
+    trees = sizehint!(Vector{F}(), n)
     K = Base.promote_op(treeindex_data, F)
     index = n > 32 ? Dict{K, Int}() : nothing
     return TreeRows{F, Dict{K, Int}}(trees, index, n)

From 5a6dcfe4c115cad61becc360347c009d749180b9 Mon Sep 17 00:00:00 2001
From: lkdvos <ldevos98@gmail.com>
Date: Thu, 5 Nov 2026 10:21:08 -0500
Subject: [PATCH 13/13] Inline the single braid channel of a tree, and presize
 the row index of large blocks

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
---
 src/fusiontrees/braiding_manipulations.jl | 6 +-----
 src/fusiontrees/fusiontrees.jl            | 2 +-
 2 files changed, 2 insertions(+), 6 deletions(-)

diff --git a/src/fusiontrees/braiding_manipulations.jl b/src/fusiontrees/braiding_manipulations.jl
index a621c43f0..b7fe74ce6 100644
--- a/src/fusiontrees/braiding_manipulations.jl
+++ b/src/fusiontrees/braiding_manipulations.jl
@@ -66,7 +66,7 @@ function artin_braid(f::FusionTree{I, N}, i; inv::Bool = false) where {I, N}
     a = inner_extended[i - 1]
     c = inner_extended[i]
     e = inner_extended[i + 1]
-    c′ = _artin_channel(a, b, d, e)
+    c′ = only(_artin_channels(a, b, d, e))
     coeff = oftype(
         oneT,
         if inv
@@ -83,10 +83,6 @@ end
 # the inner lines `c′` of `a ⊗ d → c′ → c′ ⊗ b → e`, i.e. `intersect(a ⊗ d, e ⊗ conj(b))`,
 # without constructing the second fusion product
 _artin_channels(a, b, d, e) = Iterators.filter(c′ -> !iszero(Nsymbol(c′, b, e)), a ⊗ d)
-# the only such line, for a fusion tree that braids into a single one (`only` throws otherwise)
-_artin_channel(a, b, d, e) = _artin_channel(FusionStyle(a), a, b, d, e)
-_artin_channel(::UniqueFusion, a, b, d, e) = only(a ⊗ d)
-_artin_channel(::FusionStyle, a, b, d, e) = only(_artin_channels(a, b, d, e))
 
 function artin_braid(src::FusionTreeBlock{I, N, 0}, i; inv::Bool = false) where {I, N}
     n = length(src) # braiding a complete block keeps its size
diff --git a/src/fusiontrees/fusiontrees.jl b/src/fusiontrees/fusiontrees.jl
index 32695e210..dcc63ae42 100644
--- a/src/fusiontrees/fusiontrees.jl
+++ b/src/fusiontrees/fusiontrees.jl
@@ -296,7 +296,7 @@ end
 function TreeRows(::Type{F}, n::Int) where {F}
     trees = sizehint!(Vector{F}(), n)
     K = Base.promote_op(treeindex_data, F)
-    index = n > 32 ? Dict{K, Int}() : nothing
+    index = n > 32 ? sizehint!(Dict{K, Int}(), n) : nothing
     return TreeRows{F, Dict{K, Int}}(trees, index, n)
 end
 function row!(rows::TreeRows, f)