From b30f2055eb4661e87c97024fad03beb62f0f8d58 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 08:26:39 +0200
Subject: [PATCH 01/21] specialise dim to storage type

---
 src/TensorKit.jl          |  2 +-
 src/spaces/gradedspace.jl | 17 ++++++++++++++---
 2 files changed, 15 insertions(+), 4 deletions(-)

diff --git a/src/TensorKit.jl b/src/TensorKit.jl
index 87a3a2380..929dd804e 100644
--- a/src/TensorKit.jl
+++ b/src/TensorKit.jl
@@ -44,7 +44,7 @@ export infimum, supremum, isisomorphic, ismonomorphic, isepimorphic
 export sectortype, sectors, hassector
 export unit, rightunit, leftunit, allunits, isunit, otimes, deligneproduct, timereversed
 export Nsymbol, Fsymbol, Rsymbol, Bsymbol, frobenius_schur_phase, frobenius_schur_indicator, twist, fusiontensor
-export sectorscalartype, fusionscalartype, braidingscalartype
+export sectorscalartype, fusionscalartype, braidingscalartype, dimscalartype
 
 # Export methods for fusion trees
 export fusiontrees, braid, permute, transpose
diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index 53f48dafd..d469648cc 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -89,9 +89,20 @@ GradedSpace(g::AbstractDict; dual::Bool = false) = GradedSpace(g...; dual = dual
 field(::Type{<:GradedSpace}) = ℂ
 InnerProductStyle(::Type{<:GradedSpace}) = EuclideanInnerProduct()
 
-function dim(V::GradedSpace)
-    init = 0 * dim(first(allunits(sectortype(V))))
-    return sum(c -> dim(c) * dim(V, c), sectors(V); init = init)
+function dim(V::GradedSpace{I, <:AbstractDict}) where {I <: Sector}
+    init = zero(dimscalartype(I))
+    return sum(((c, d),) -> dim(c) * d, V.dims; init)
+end
+function dim(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N}
+    init = zero(dimscalartype(I))
+    D = init
+    vals = values(I)
+    @inbounds for n in 1:N
+        d = V.dims[n]
+        iszero(d) && continue
+        D += dim(vals[n]) * d # dim(c) = dim(dual(c))
+    end
+    return D
 end
 function dim(V::GradedSpace{I, <:AbstractDict}, c::I) where {I <: Sector}
     return get(V.dims, isdual(V) ? dual(c) : c, 0)

From 41f43fdf997580b4128ee22cf9291570967f826e Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 12:13:17 +0200
Subject: [PATCH 02/21] oplus and ominus

---
 src/spaces/gradedspace.jl | 72 +++++++++++++++++++++++++++++++--------
 1 file changed, 57 insertions(+), 15 deletions(-)

diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index d469648cc..6064eeaef 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -137,25 +137,67 @@ function unitspace(S::Type{<:GradedSpace{I}}) where {I <: Sector}
 end
 zerospace(S::Type{<:GradedSpace}) = S()
 
-# TODO: the following methods can probably be implemented more efficiently for
-# `FiniteGradedSpace`, but we don't expect them to be used often in hot loops, so
-# these generic definitions (which are still quite efficient) are good for now.
-function ⊕(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector}
+function ⊕(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
+    dual1 = isdual(V₁)
+    dual1 == isdual(V₂) || throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist"))
+    k1, k2 = V₁.dims.keys, V₂.dims.keys # already sorted
+    v1, v2 = V₁.dims.values, V₂.dims.values
+    n1, n2 = length(k1), length(k2)
+    ks, vs = Vector{I}(), Vector{Int}()
+    sizehint!(ks, n1 + n2)
+    sizehint!(vs, n1 + n2)
+    i, j = 1, 1
+    @inbounds while i <= n1 && j <= n2
+        if k1[i] == k2[j]
+            push!(ks, k1[i]); push!(vs, v1[i] + v2[j]); i += 1; j += 1
+        elseif k1[i] < k2[j]
+            push!(ks, k1[i]); push!(vs, v1[i]); i += 1
+        else
+            push!(ks, k2[j]); push!(vs, v2[j]); j += 1
+        end
+    end
+    @inbounds while i <= n1
+        push!(ks, k1[i]); push!(vs, v1[i]); i += 1
+    end
+    @inbounds while j <= n2
+        push!(ks, k2[j]); push!(vs, v2[j]); j += 1
+    end
+    return typeof(V₁)(SectorDict{I, Int}(ks, vs), dual1)
+end
+function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector}
     dual1 = isdual(V₁)
     dual1 == isdual(V₂) ||
         throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist"))
-    dims = SectorDict{I, Int}()
-    for c in union(sectors(V₁), sectors(V₂))
-        cout = ifelse(dual1, dual(c), c)
-        dims[cout] = dim(V₁, c) + dim(V₂, c)
+    newdims = map(+, V₁.dims, V₂.dims)
+    return typeof(V₁)(newdims, dual1)
+end
+function ⊖(V::GradedSpace{I, <: Tuple}, W::GradedSpace{I, <: Tuple}) where {I <: Sector}
+    dualV = isdual(V)
+    V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)"))
+    newdims = map(-, V.dims, W.dims)
+    return typeof(V)(newdims, dualV)
+end
+function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) where {I <: Sector}
+    dualV = isdual(V)
+    V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)"))
+    kv, kw = V.dims.keys, W.dims.keys # already sorted
+    vv, vw = V.dims.values, W.dims.values
+    ks, vs = Vector{I}(), Vector{Int}()
+    nv, nw = length(kv), length(kw)
+    sizehint!(ks, nv)
+    sizehint!(vs, nv)
+    j = 1
+    @inbounds for i in eachindex(kv) # keys(W) ⊆ keys(V)
+        d = vv[i]
+        if j <= nw && kw[j] == kv[i]
+            d -= vw[j]
+            j += 1
+        end
+        if !iszero(d)
+            push!(ks, kv[i]); push!(vs, d)
+        end
     end
-    return typeof(V₁)(dims; dual = dual1)
-end
-function ⊖(V::GradedSpace{I}, W::GradedSpace{I}) where {I <: Sector}
-    dual = isdual(V)
-    V ≿ W && dual == isdual(W) ||
-        throw(SpaceMismatch("$(W) is not a subspace of $(V)"))
-    return typeof(V)(c => dim(V, c) - dim(W, c) for c in sectors(V); dual)
+    return typeof(V)(SectorDict{I, Int}(ks, vs), dualV)
 end
 
 function fuse(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector}

From 619804db43336769f4031ba76413944befa24d87 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 12:21:28 +0200
Subject: [PATCH 03/21] infimum and supremum

---
 src/spaces/gradedspace.jl | 76 ++++++++++++++++++++++++++++++---------
 1 file changed, 60 insertions(+), 16 deletions(-)

diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index 6064eeaef..46c94eaf7 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -171,7 +171,7 @@ function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where
     newdims = map(+, V₁.dims, V₂.dims)
     return typeof(V₁)(newdims, dual1)
 end
-function ⊖(V::GradedSpace{I, <: Tuple}, W::GradedSpace{I, <: Tuple}) where {I <: Sector}
+function ⊖(V::GradedSpace{I, <:Tuple}, W::GradedSpace{I, <:Tuple}) where {I <: Sector}
     dualV = isdual(V)
     V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)"))
     newdims = map(-, V.dims, W.dims)
@@ -210,23 +210,67 @@ function fuse(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector}
     return typeof(V₁)(dims)
 end
 
-function infimum(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector}
+function infimum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector}
     Visdual = isdual(V₁)
-    Visdual == isdual(V₂) ||
-        throw(SpaceMismatch("Infimum of space and dual space does not exist"))
-    return typeof(V₁)(
-        (Visdual ? dual(c) : c) => min(dim(V₁, c), dim(V₂, c))
-            for c in intersect(sectors(V₁), sectors(V₂)); dual = Visdual
-    )
-end
-function supremum(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector}
+    Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist"))
+    newdims = map(min, V₁.dims, V₂.dims)
+    return typeof(V₁)(newdims, Visdual)
+end
+function infimum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
+    Visdual = isdual(V₁)
+    Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist"))
+    k1, k2 = V₁.dims.keys, V₂.dims.keys
+    v1, v2 = V₁.dims.values, V₂.dims.values
+    n1, n2 = length(k1), length(k2)
+    ks, vs = Vector{I}(), Vector{Int}()
+    i, j = 1, 1
+    @inbounds while i <= n1 && j <= n2
+        if k1[i] == k2[j]
+            m = min(v1[i], v2[j])
+            if !iszero(m)
+                push!(ks, k1[i]); push!(vs, m)
+            end
+            i += 1; j += 1
+        elseif k1[i] < k2[j]
+            i += 1
+        else
+            j += 1
+        end
+    end
+    return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual)
+end
+function supremum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector}
     Visdual = isdual(V₁)
-    Visdual == isdual(V₂) ||
-        throw(SpaceMismatch("Supremum of space and dual space does not exist"))
-    return typeof(V₁)(
-        (Visdual ? dual(c) : c) => max(dim(V₁, c), dim(V₂, c))
-            for c in union(sectors(V₁), sectors(V₂)); dual = Visdual
-    )
+    Visdual == isdual(V₂) || throw(SpaceMismatch("Supremum of space and dual space does not exist"))
+    newdims = map(max, V₁.dims, V₂.dims)
+    return typeof(V₁)(newdims, Visdual)
+end
+function supremum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
+    Visdual = isdual(V₁)
+    Visdual == isdual(V₂) || throw(SpaceMismatch("Supremum of space and dual space does not exist"))
+    k1, k2 = V₁.dims.keys, V₂.dims.keys
+    v1, v2 = V₁.dims.values, V₂.dims.values
+    n1, n2 = length(k1), length(k2)
+    ks, vs = Vector{I}(), Vector{Int}()
+    sizehint!(ks, n1 + n2)
+    sizehint!(vs, n1 + n2)
+    i, j = 1, 1
+    @inbounds while i <= n1 && j <= n2
+        if k1[i] == k2[j]
+            push!(ks, k1[i]); push!(vs, max(v1[i], v2[j])); i += 1; j += 1
+        elseif k1[i] < k2[j]
+            push!(ks, k1[i]); push!(vs, v1[i]); i += 1
+        else
+            push!(ks, k2[j]); push!(vs, v2[j]); j += 1
+        end
+    end
+    @inbounds while i <= n1
+        push!(ks, k1[i]); push!(vs, v1[i]); i += 1
+    end
+    @inbounds while j <= n2
+        push!(ks, k2[j]); push!(vs, v2[j]); j += 1
+    end
+    return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual)
 end
 
 hassector(V::GradedSpace{I}, s::I) where {I <: Sector} = dim(V, s) != 0

From 0967bd071455860e189687dfb347f267736d9e1e Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 12:52:02 +0200
Subject: [PATCH 04/21] fuse

---
 src/spaces/gradedspace.jl | 47 ++++++++++++++++++++++++++++++++++-----
 1 file changed, 41 insertions(+), 6 deletions(-)

diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index 46c94eaf7..d32358405 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -200,14 +200,49 @@ function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) w
     return typeof(V)(SectorDict{I, Int}(ks, vs), dualV)
 end
 
-function fuse(V₁::GradedSpace{I}, V₂::GradedSpace{I}) where {I <: Sector}
-    dims = SectorDict{I, Int}()
-    for a in sectors(V₁), b in sectors(V₂)
-        for c in a ⊗ b
-            dims[c] = get(dims, c, 0) + Nsymbol(a, b, c) * dim(V₁, a) * dim(V₂, b)
+function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
+    dual1, dual2 = isdual(V₁), isdual(V₂)
+    acc = Dict{I, Int}() # SectorDict `get` within the double for loop accumulates O(N^2) ` findindex` calls -> sort afterwards
+    k1, k2 = V₁.dims.keys, V₂.dims.keys
+    v1, v2 = V₁.dims.values, V₂.dims.values
+    @inbounds for n1 in eachindex(k1)
+        a0 = k1[n1]; d1 = v1[n1]
+        a = dual1 ? dual(a0) : a0
+        for n2 in eachindex(k2)
+            b0 = k2[n2]; d2 = v2[n2]
+            b = dual2 ? dual(b0) : b0
+            dab = d1 * d2
+            for c in a ⊗ b
+                acc[c] = get(acc, c, 0) + Nsymbol(a, b, c) * dab
+            end
+        end
+    end
+    ks = sort!(collect(keys(acc)))
+    vs = [acc[k] for k in ks]
+    return typeof(V₁)(SectorDict{I, Int}(ks, vs), false)
+end
+function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N}
+    vals = values(I)
+    dual1, dual2 = isdual(V₁), isdual(V₂)
+    newdims = zeros(Int, N) #TODO: is there a way to avoid dense storage even for sparse results?
+    @inbounds for n1 in 1:N
+        d1 = V₁.dims[n1]
+        iszero(d1) && continue
+        a0 = vals[n1] # avoid call to sectors(V₁)
+        a = dual1 ? dual(a0) : a0
+        for n2 in 1:N
+            d2 = V₂.dims[n2]
+            iszero(d2) && continue
+            b0 = vals[n2] # idem for V₂
+            b = dual2 ? dual(b0) : b0
+            dab = d1 * d2
+            for c in a ⊗ b
+                nc = findindex(vals, c)
+                newdims[nc] += Nsymbol(a, b, c) * dab
+            end
         end
     end
-    return typeof(V₁)(dims)
+    return typeof(V₁)(ntuple(i -> newdims[i], Val(N)), false)
 end
 
 function infimum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector}

From 60186d55213b442c3f46c1d0c3753e2c660e1050 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 13:14:49 +0200
Subject: [PATCH 05/21] truncate_space

---
 src/factorizations/truncation.jl | 25 +++++++++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl
index 2facfd2e9..67dac8e56 100644
--- a/src/factorizations/truncation.jl
+++ b/src/factorizations/truncation.jl
@@ -37,6 +37,31 @@ _blocklength(ax::Base.OneTo, ind::AbstractVector{Bool}) = count(ind)
 function truncate_space(V::ElementarySpace, inds)
     return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in pairs(inds))
 end
+function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N}
+    vals = values(I)
+    dualV = isdual(V)
+    newdims = zeros(Int, N)
+    for (c, ind) in pairs(inds)
+        n_read = findindex(vals, dualV ? dual(c) : c) # dual-adjusted index for reading V.dims
+        n_write = findindex(vals, c) # output is never dual, so c is fine as-is
+        newdims[n_write] = _blocklength(V.dims[n_read], ind) # dim(c) = dim(dual(c))
+    end
+    return typeof(V)(ntuple(i -> newdims[i], Val(N)), false)
+end
+function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector}
+    dualV = isdual(V)
+    ks, vs = Vector{I}(), Vector{Int}() # accumulate and sort once at the end
+    for (c, ind) in pairs(inds)
+        d = get(V.dims, dualV ? dual(c) : c, 0)
+        len = _blocklength(d, ind)
+        if !iszero(len)
+            push!(ks, c)
+            push!(vs, len)
+        end
+    end
+    perm = sortperm(ks)
+    return typeof(V)(SectorDict{I, Int}(ks[perm], vs[perm]), false)
+end
 
 function truncate_domain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds)
     for (c, b) in blocks(tdst)

From 8f90be2ddaf82a2d25de10c88abcc859d5d904d2 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 13:20:36 +0200
Subject: [PATCH 06/21] restore binary search for sectordicts

---
 src/auxiliary/dicts.jl | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl
index 626599c01..d495103eb 100644
--- a/src/auxiliary/dicts.jl
+++ b/src/auxiliary/dicts.jl
@@ -89,14 +89,14 @@ end
 Base.empty(::SortedVectorDict, ::Type{K}, ::Type{V}) where {K, V} = SortedVectorDict{K, V}()
 Base.empty!(d::SortedVectorDict) = (empty!(d.keys); empty!(d.values); return d)
 
-# _searchsortedfirst(v::Vector, k) = searchsortedfirst(v, k)
-function _searchsortedfirst(v::Vector, k)
-    i = 1
-    @inbounds while i <= length(v) && isless(v[i], k)
-        i += 1
-    end
-    return i
-end
+_searchsortedfirst(v::Vector, k) = searchsortedfirst(v, k)
+# function _searchsortedfirst(v::Vector, k)
+#     i = 1
+#     @inbounds while i <= length(v) && isless(v[i], k)
+#         i += 1
+#     end
+#     return i
+# end
 
 function Base.delete!(d::SortedVectorDict{K}, k) where {K}
     key = convert(K, k)

From ae10f5b565cc93c26caa33eab00d4c4ca8635e2f Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 13:54:36 +0200
Subject: [PATCH 07/21] refactor sorted merge procedure

---
 src/auxiliary/dicts.jl    | 49 ++++++++++++++++++++++
 src/spaces/gradedspace.jl | 87 +++------------------------------------
 2 files changed, 55 insertions(+), 81 deletions(-)

diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl
index d495103eb..17297177e 100644
--- a/src/auxiliary/dicts.jl
+++ b/src/auxiliary/dicts.jl
@@ -186,6 +186,55 @@ function Base.:(==)(d1::SortedVectorDict, d2::SortedVectorDict)
     return true
 end
 
+# merge over two SORTED vector pairs representing keys and values
+# - combine(v1,v2): value for a key present in both operands
+# - only1(v1) / only2(v2): value for a key present in only one operand;
+# pass `nothing` to drop such keys entirely (e.g. for an intersection)
+# zero results are dropped (either from `combine` or `only1`/`only2`), matching how GradedSpace never stores an explicit zero dimension
+# k1 and k2 originate from GradedSpace.dims.keys, which are guaranteed to be sorted
+function _sortedmerge(k1::Vector{I}, v1::Vector{Int}, k2::Vector{I}, v2::Vector{Int}, combine, only1, only2) where {I}
+    n1, n2 = length(k1), length(k2)
+    ks, vs = Vector{I}(), Vector{Int}()
+    sizehint!(ks, n1 + n2)
+    sizehint!(vs, n1 + n2)
+    i, j = 1, 1
+    @inbounds while i <= n1 && j <= n2
+        if k1[i] == k2[j]
+            d = combine(v1[i], v2[j])
+            if !iszero(d)
+                push!(ks, k1[i])
+                push!(vs, d)
+            end
+            i += 1
+            j += 1
+        elseif k1[i] < k2[j]
+            _mergeonly!(ks, vs, k1[i], v1[i], only1)
+            i += 1
+        else
+            _mergeonly!(ks, vs, k2[j], v2[j], only2)
+            j += 1
+        end
+    end
+    @inbounds while i <= n1
+        _mergeonly!(ks, vs, k1[i], v1[i], only1)
+        i += 1
+    end
+    @inbounds while j <= n2
+        _mergeonly!(ks, vs, k2[j], v2[j], only2)
+        j += 1
+    end
+    return ks, vs
+end
+@inline _mergeonly!(ks, vs, k, v, ::Nothing) = nothing
+@inline function _mergeonly!(ks, vs, k, v, f)
+    d = f(v)
+    if !iszero(d)
+        push!(ks, k)
+        push!(vs, d)
+    end
+    return nothing
+end
+
 """
     Hashed(value, hashfunction = Base.hash, isequal = Base.isequal)
 
diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index d32358405..5ca39f912 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -140,28 +140,7 @@ zerospace(S::Type{<:GradedSpace}) = S()
 function ⊕(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
     dual1 = isdual(V₁)
     dual1 == isdual(V₂) || throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist"))
-    k1, k2 = V₁.dims.keys, V₂.dims.keys # already sorted
-    v1, v2 = V₁.dims.values, V₂.dims.values
-    n1, n2 = length(k1), length(k2)
-    ks, vs = Vector{I}(), Vector{Int}()
-    sizehint!(ks, n1 + n2)
-    sizehint!(vs, n1 + n2)
-    i, j = 1, 1
-    @inbounds while i <= n1 && j <= n2
-        if k1[i] == k2[j]
-            push!(ks, k1[i]); push!(vs, v1[i] + v2[j]); i += 1; j += 1
-        elseif k1[i] < k2[j]
-            push!(ks, k1[i]); push!(vs, v1[i]); i += 1
-        else
-            push!(ks, k2[j]); push!(vs, v2[j]); j += 1
-        end
-    end
-    @inbounds while i <= n1
-        push!(ks, k1[i]); push!(vs, v1[i]); i += 1
-    end
-    @inbounds while j <= n2
-        push!(ks, k2[j]); push!(vs, v2[j]); j += 1
-    end
+    ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, +, identity, identity)
     return typeof(V₁)(SectorDict{I, Int}(ks, vs), dual1)
 end
 function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector}
@@ -180,23 +159,7 @@ end
 function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) where {I <: Sector}
     dualV = isdual(V)
     V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)"))
-    kv, kw = V.dims.keys, W.dims.keys # already sorted
-    vv, vw = V.dims.values, W.dims.values
-    ks, vs = Vector{I}(), Vector{Int}()
-    nv, nw = length(kv), length(kw)
-    sizehint!(ks, nv)
-    sizehint!(vs, nv)
-    j = 1
-    @inbounds for i in eachindex(kv) # keys(W) ⊆ keys(V)
-        d = vv[i]
-        if j <= nw && kw[j] == kv[i]
-            d -= vw[j]
-            j += 1
-        end
-        if !iszero(d)
-            push!(ks, kv[i]); push!(vs, d)
-        end
-    end
+    ks, vs = _sortedmerge(V.dims.keys, V.dims.values, W.dims.keys, W.dims.values, -, identity, nothing)
     return typeof(V)(SectorDict{I, Int}(ks, vs), dualV)
 end
 
@@ -217,14 +180,14 @@ function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorD
             end
         end
     end
-    ks = sort!(collect(keys(acc)))
+    ks = sort!(collect(keys(acc))) #TODO: sortperm?
     vs = [acc[k] for k in ks]
     return typeof(V₁)(SectorDict{I, Int}(ks, vs), false)
 end
 function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N}
     vals = values(I)
     dual1, dual2 = isdual(V₁), isdual(V₂)
-    newdims = zeros(Int, N) #TODO: is there a way to avoid dense storage even for sparse results?
+    newdims = zeros(Int, N)
     @inbounds for n1 in 1:N
         d1 = V₁.dims[n1]
         iszero(d1) && continue
@@ -254,24 +217,7 @@ end
 function infimum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
     Visdual = isdual(V₁)
     Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist"))
-    k1, k2 = V₁.dims.keys, V₂.dims.keys
-    v1, v2 = V₁.dims.values, V₂.dims.values
-    n1, n2 = length(k1), length(k2)
-    ks, vs = Vector{I}(), Vector{Int}()
-    i, j = 1, 1
-    @inbounds while i <= n1 && j <= n2
-        if k1[i] == k2[j]
-            m = min(v1[i], v2[j])
-            if !iszero(m)
-                push!(ks, k1[i]); push!(vs, m)
-            end
-            i += 1; j += 1
-        elseif k1[i] < k2[j]
-            i += 1
-        else
-            j += 1
-        end
-    end
+    ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, min, nothing, nothing)
     return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual)
 end
 function supremum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector}
@@ -283,28 +229,7 @@ end
 function supremum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
     Visdual = isdual(V₁)
     Visdual == isdual(V₂) || throw(SpaceMismatch("Supremum of space and dual space does not exist"))
-    k1, k2 = V₁.dims.keys, V₂.dims.keys
-    v1, v2 = V₁.dims.values, V₂.dims.values
-    n1, n2 = length(k1), length(k2)
-    ks, vs = Vector{I}(), Vector{Int}()
-    sizehint!(ks, n1 + n2)
-    sizehint!(vs, n1 + n2)
-    i, j = 1, 1
-    @inbounds while i <= n1 && j <= n2
-        if k1[i] == k2[j]
-            push!(ks, k1[i]); push!(vs, max(v1[i], v2[j])); i += 1; j += 1
-        elseif k1[i] < k2[j]
-            push!(ks, k1[i]); push!(vs, v1[i]); i += 1
-        else
-            push!(ks, k2[j]); push!(vs, v2[j]); j += 1
-        end
-    end
-    @inbounds while i <= n1
-        push!(ks, k1[i]); push!(vs, v1[i]); i += 1
-    end
-    @inbounds while j <= n2
-        push!(ks, k2[j]); push!(vs, v2[j]); j += 1
-    end
+    ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, max, identity, identity)
     return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual)
 end
 

From 47a5f3437054d7533c454276ca95f604899e6b0a Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 15:21:40 +0200
Subject: [PATCH 08/21] speed up fuse slightly with sortperm

---
 src/spaces/gradedspace.jl | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index 5ca39f912..0a79293c2 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -180,9 +180,10 @@ function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorD
             end
         end
     end
-    ks = sort!(collect(keys(acc))) #TODO: sortperm?
-    vs = [acc[k] for k in ks]
-    return typeof(V₁)(SectorDict{I, Int}(ks, vs), false)
+    ks0 = collect(keys(acc))
+    vs0 = collect(values(acc))
+    perm = sortperm(ks0)
+    return typeof(V₁)(SectorDict{I, Int}(ks0[perm], vs0[perm]), false)
 end
 function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N}
     vals = values(I)

From 411dc4472b20d214bb851e53b1468a5938a7c29d Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 11 Aug 2026 15:21:55 +0200
Subject: [PATCH 09/21] import thing

---
 src/factorizations/factorizations.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl
index fbe87a63a..80d96ff9a 100644
--- a/src/factorizations/factorizations.jl
+++ b/src/factorizations/factorizations.jl
@@ -6,7 +6,7 @@ module Factorizations
 export copy_oftype, factorisation_scalartype, one!, truncspace
 
 using ..TensorKit
-using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector,
+using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, findindex,
     blocktype, foreachblock, one!,
     similar_diagonal, similarstoragetype
 

From 0adec3ba4357537db2857dac9c25125673302a0d Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Wed, 12 Aug 2026 11:47:24 +0200
Subject: [PATCH 10/21] splat with type annotation above Val

---
 src/factorizations/truncation.jl |  2 +-
 src/spaces/gradedspace.jl        | 14 +++++++-------
 2 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl
index 67dac8e56..f3c64ed7c 100644
--- a/src/factorizations/truncation.jl
+++ b/src/factorizations/truncation.jl
@@ -46,7 +46,7 @@ function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sec
         n_write = findindex(vals, c) # output is never dual, so c is fine as-is
         newdims[n_write] = _blocklength(V.dims[n_read], ind) # dim(c) = dim(dual(c))
     end
-    return typeof(V)(ntuple(i -> newdims[i], Val(N)), false)
+    return typeof(V)((newdims...,)::NTuple{N, Int}, false)
 end
 function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector}
     dualV = isdual(V)
diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index 0a79293c2..3e363d19c 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -30,17 +30,17 @@ end
 sectortype(::Type{<:GradedSpace{I}}) where {I <: Sector} = I
 
 function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N}
-    d = ntuple(n -> 0, N)
-    isset = ntuple(n -> false, N)
+    d = zeros(Int, N)
+    isset = falses(N)
     for (c, dc) in dims
         k = convert(I, c)
         i = findindex(values(I), k)
-        k = dc < 0 && throw(ArgumentError(lazy"Sector $k has negative dimension $dc"))
+        dc < 0 && throw(ArgumentError(lazy"Sector $k has negative dimension $dc"))
         isset[i] && throw(ArgumentError(lazy"Sector $c appears multiple times"))
-        isset = TupleTools.setindex(isset, true, i)
-        d = TupleTools.setindex(d, dc, i)
+        isset[i] = true
+        d[i] = dc
     end
-    return GradedSpace{I, NTuple{N, Int}}(d, dual)
+    return GradedSpace{I, NTuple{N, Int}}((d...,)::NTuple{N, Int}, dual)
 end
 function GradedSpace{I, NTuple{N, Int}}(dims::Pair; dual::Bool = false) where {I, N}
     return GradedSpace{I, NTuple{N, Int}}((dims,); dual = dual)
@@ -206,7 +206,7 @@ function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{
             end
         end
     end
-    return typeof(V₁)(ntuple(i -> newdims[i], Val(N)), false)
+    return typeof(V₁)((newdims...,)::NTuple{N, Int}, false)
 end
 
 function infimum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector}

From 13677f3129b19757c7ce4fd72a0129b577c64e88 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 21 Aug 2026 15:33:03 +0200
Subject: [PATCH 11/21] actually don't splat, but construct directly where
 previously a vector was made

---
 src/factorizations/truncation.jl | 2 +-
 src/spaces/gradedspace.jl        | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl
index f3c64ed7c..3ec0457a9 100644
--- a/src/factorizations/truncation.jl
+++ b/src/factorizations/truncation.jl
@@ -46,7 +46,7 @@ function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sec
         n_write = findindex(vals, c) # output is never dual, so c is fine as-is
         newdims[n_write] = _blocklength(V.dims[n_read], ind) # dim(c) = dim(dual(c))
     end
-    return typeof(V)((newdims...,)::NTuple{N, Int}, false)
+    return typeof(V)(NTuple{N, Int}(newdims), false)
 end
 function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector}
     dualV = isdual(V)
diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index 3e363d19c..ede620e62 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -40,7 +40,7 @@ function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N}
         isset[i] = true
         d[i] = dc
     end
-    return GradedSpace{I, NTuple{N, Int}}((d...,)::NTuple{N, Int}, dual)
+    return GradedSpace{I, NTuple{N, Int}}(NTuple{N, Int}(d), dual)
 end
 function GradedSpace{I, NTuple{N, Int}}(dims::Pair; dual::Bool = false) where {I, N}
     return GradedSpace{I, NTuple{N, Int}}((dims,); dual = dual)
@@ -206,7 +206,7 @@ function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{
             end
         end
     end
-    return typeof(V₁)((newdims...,)::NTuple{N, Int}, false)
+    return typeof(V₁)(NTuple{N, Int}(newdims), false)
 end
 
 function infimum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector}

From 399e121f677afbd330df396feaed54b97918c2c4 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 21 Aug 2026 16:17:06 +0200
Subject: [PATCH 12/21] make slightly more readable maybe perhaps

---
 src/spaces/gradedspace.jl | 36 ++++++++++++++++++------------------
 1 file changed, 18 insertions(+), 18 deletions(-)

diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index ede620e62..bdea064e1 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -168,13 +168,13 @@ function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorD
     acc = Dict{I, Int}() # SectorDict `get` within the double for loop accumulates O(N^2) ` findindex` calls -> sort afterwards
     k1, k2 = V₁.dims.keys, V₂.dims.keys
     v1, v2 = V₁.dims.values, V₂.dims.values
-    @inbounds for n1 in eachindex(k1)
-        a0 = k1[n1]; d1 = v1[n1]
-        a = dual1 ? dual(a0) : a0
-        for n2 in eachindex(k2)
-            b0 = k2[n2]; d2 = v2[n2]
-            b = dual2 ? dual(b0) : b0
-            dab = d1 * d2
+    @inbounds for na in eachindex(k1)
+        a₀, da = k1[na], v1[na]
+        a = dual1 ? dual(a₀) : a₀
+        for nb in eachindex(k2)
+            b₀, db = k2[nb], v2[nb]
+            b = dual2 ? dual(b₀) : b₀
+            dab = da * db
             for c in a ⊗ b
                 acc[c] = get(acc, c, 0) + Nsymbol(a, b, c) * dab
             end
@@ -189,17 +189,17 @@ function fuse(V₁::GradedSpace{I, NTuple{N, Int}}, V₂::GradedSpace{I, NTuple{
     vals = values(I)
     dual1, dual2 = isdual(V₁), isdual(V₂)
     newdims = zeros(Int, N)
-    @inbounds for n1 in 1:N
-        d1 = V₁.dims[n1]
-        iszero(d1) && continue
-        a0 = vals[n1] # avoid call to sectors(V₁)
-        a = dual1 ? dual(a0) : a0
-        for n2 in 1:N
-            d2 = V₂.dims[n2]
-            iszero(d2) && continue
-            b0 = vals[n2] # idem for V₂
-            b = dual2 ? dual(b0) : b0
-            dab = d1 * d2
+    @inbounds for na in 1:N
+        da = V₁.dims[na]
+        iszero(da) && continue
+        a₀ = vals[na] # avoid call to sectors(V₁)
+        a = dual1 ? dual(a₀) : a₀
+        for nb in 1:N
+            db = V₂.dims[nb]
+            iszero(db) && continue
+            b₀ = vals[nb] # idem for V₂
+            b = dual2 ? dual(b₀) : b₀
+            dab = da * db
             for c in a ⊗ b
                 nc = findindex(vals, c)
                 newdims[nc] += Nsymbol(a, b, c) * dab

From d2a41f16fff0ff10108dda3966433a68b1269ef0 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 25 Aug 2026 17:46:59 +0200
Subject: [PATCH 13/21] truncate_space always has non-dual entry spaces

---
 src/factorizations/truncation.jl | 10 ++++------
 1 file changed, 4 insertions(+), 6 deletions(-)

diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl
index 3ec0457a9..87971ca88 100644
--- a/src/factorizations/truncation.jl
+++ b/src/factorizations/truncation.jl
@@ -39,20 +39,18 @@ function truncate_space(V::ElementarySpace, inds)
 end
 function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N}
     vals = values(I)
-    dualV = isdual(V)
     newdims = zeros(Int, N)
     for (c, ind) in pairs(inds)
-        n_read = findindex(vals, dualV ? dual(c) : c) # dual-adjusted index for reading V.dims
-        n_write = findindex(vals, c) # output is never dual, so c is fine as-is
-        newdims[n_write] = _blocklength(V.dims[n_read], ind) # dim(c) = dim(dual(c))
+        d = dim(V, c)
+        n_write = findindex(vals, c)
+        newdims[n_write] = _blocklength(d, ind)
     end
     return typeof(V)(NTuple{N, Int}(newdims), false)
 end
 function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector}
-    dualV = isdual(V)
     ks, vs = Vector{I}(), Vector{Int}() # accumulate and sort once at the end
     for (c, ind) in pairs(inds)
-        d = get(V.dims, dualV ? dual(c) : c, 0)
+        d = dim(V, c)
         len = _blocklength(d, ind)
         if !iszero(len)
             push!(ks, c)

From d7bdb07189022986ebff1be5a6846ff1fad20775 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 25 Aug 2026 18:04:44 +0200
Subject: [PATCH 14/21] overkill iszero check in dim

---
 src/spaces/gradedspace.jl | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index bdea064e1..cde32d709 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -97,10 +97,8 @@ function dim(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N}
     init = zero(dimscalartype(I))
     D = init
     vals = values(I)
-    @inbounds for n in 1:N
-        d = V.dims[n]
-        iszero(d) && continue
-        D += dim(vals[n]) * d # dim(c) = dim(dual(c))
+    @inbounds for (c, d) in zip(values(I), V.dims)
+        D += dim(c) * d
     end
     return D
 end

From 271e96b4bfe1004661712c30652b934e7288c357 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Wed, 26 Aug 2026 09:48:22 +0200
Subject: [PATCH 15/21] assert truncate_space spaces being non-dual

---
 src/factorizations/truncation.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl
index 87971ca88..03e636298 100644
--- a/src/factorizations/truncation.jl
+++ b/src/factorizations/truncation.jl
@@ -35,9 +35,11 @@ _blocklength(ax::Base.OneTo, ind::AbstractVector{<:Integer}) = length(ind)
 _blocklength(ax::Base.OneTo, ind::AbstractVector{Bool}) = count(ind)
 
 function truncate_space(V::ElementarySpace, inds)
+    @assert !isdual(V)
     return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in pairs(inds))
 end
 function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N}
+    @assert !isdual(V)
     vals = values(I)
     newdims = zeros(Int, N)
     for (c, ind) in pairs(inds)
@@ -48,6 +50,7 @@ function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sec
     return typeof(V)(NTuple{N, Int}(newdims), false)
 end
 function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Sector}
+    @assert !isdual(V)
     ks, vs = Vector{I}(), Vector{Int}() # accumulate and sort once at the end
     for (c, ind) in pairs(inds)
         d = dim(V, c)

From 335785f282cf6ee21a856c04e5c3f479d06c13a3 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 28 Aug 2026 17:18:18 +0200
Subject: [PATCH 16/21] introduce `sectorstoragetype`

---
 src/spaces/gradedspace.jl | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index cde32d709..b7b35c8f9 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -323,12 +323,20 @@ specify `D`.
 const Vect = SpaceTable()
 Base.getindex(::SpaceTable) = ComplexSpace
 Base.getindex(::SpaceTable, ::Type{Trivial}) = ComplexSpace
-function Base.getindex(::SpaceTable, I::Type{<:Sector})
+Base.getindex(::SpaceTable, I::Type{<:Sector}) = GradedSpace{I, sectorstoragetype(I)}
+
+"""
+    sectorstoragetype(I::Type{<:Sector}) -> Type
+
+The storage type `D` used for the `dims` field of `GradedSpace{I, D}`.
+This is `NTuple{N,Int}` with `N = length(values(I))` if `I` has a finite, known length,
+or `SectorDict{I,Int}` otherwise.
+"""
+Base.@assume_effects :foldable function sectorstoragetype(::Type{I}) where {I <: Sector}
     if Base.IteratorSize(values(I)) isa Union{HasLength, HasShape}
-        N = length(values(I))
-        return GradedSpace{I, NTuple{N, Int}}
+        return NTuple{length(values(I)), Int}
     else
-        return GradedSpace{I, SectorDict{I, Int}}
+        return SectorDict{I, Int}
     end
 end
 

From d79931ec59734c0dbefb589d83010b295b586677 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 28 Aug 2026 18:12:59 +0200
Subject: [PATCH 17/21] specialise truncation code to sectorstoragetype

---
 src/factorizations/factorizations.jl |   2 +-
 src/factorizations/pullbacks.jl      |  12 ++-
 src/factorizations/truncation.jl     | 152 ++++++++++++++++++++-------
 3 files changed, 123 insertions(+), 43 deletions(-)

diff --git a/src/factorizations/factorizations.jl b/src/factorizations/factorizations.jl
index 80d96ff9a..b16a28e9b 100644
--- a/src/factorizations/factorizations.jl
+++ b/src/factorizations/factorizations.jl
@@ -8,7 +8,7 @@ export copy_oftype, factorisation_scalartype, one!, truncspace
 using ..TensorKit
 using ..TensorKit: AdjointTensorMap, SectorDict, SectorVector, findindex,
     blocktype, foreachblock, one!,
-    similar_diagonal, similarstoragetype
+    similar_diagonal, similarstoragetype, sectorstoragetype
 
 using LinearAlgebra: LinearAlgebra, BlasFloat, Diagonal,
     svdvals, svdvals!, eigen, eigen!,
diff --git a/src/factorizations/pullbacks.jl b/src/factorizations/pullbacks.jl
index b910a184c..68ac57bb4 100644
--- a/src/factorizations/pullbacks.jl
+++ b/src/factorizations/pullbacks.jl
@@ -24,16 +24,22 @@ for pullback! in (:qr_null_pullback!, :lq_null_pullback!)
         return Δt
     end
 end
-_notrunc_ind(t) = SectorDict(c => Colon() for c in blocksectors(t))
+function _notrunc_ind(t)
+    I = sectortype(t)
+    return _builddensemap(sectorstoragetype(I), I, blocks(t), Colon) do _, _
+        Colon()
+    end
+end
 
 for pullback! in (:svd_pullback!, :eig_pullback!, :eigh_pullback!)
     @eval function MAK.$pullback!(
             Δt::AbstractTensorMap, t::AbstractTensorMap, F, ΔF, inds = _notrunc_ind(t);
             kwargs...
         )
+        Isec = sectortype(t)
         foreachblock(Δt, t) do c, (Δb, b)
-            haskey(inds, c) || return nothing
-            ind = inds[c]
+            ind = _denseget(inds, Isec, c)
+            isnothing(ind) && return nothing
             Fc = block.(F, Ref(c))
             ΔFc = block.(ΔF, Ref(c))
             MAK.$pullback!(Δb, b, Fc, ΔFc, ind; kwargs...)
diff --git a/src/factorizations/truncation.jl b/src/factorizations/truncation.jl
index 03e636298..d2e455ebd 100644
--- a/src/factorizations/truncation.jl
+++ b/src/factorizations/truncation.jl
@@ -34,15 +34,55 @@ _blocklength(ax, ind) = length(ax[ind])
 _blocklength(ax::Base.OneTo, ind::AbstractVector{<:Integer}) = length(ind)
 _blocklength(ax::Base.OneTo, ind::AbstractVector{Bool}) = count(ind)
 
+# TODO: it quacks like a duck, just define a subtype of AbstractDict?
+# represent the sector-index mapping as Vector{Union{Nothing, V}} where V is the type of the index
+# mapping is indexed through findindex
+# the type V is needed because the concrete type of ind depends on the strategy (except for intersect/union)
+_densenew(::Type{I}, ::Type{V}) where {I <: Sector, V} =
+    Vector{Union{Nothing, V}}(nothing, length(values(I)))
+
+function _denseset!(v::Vector, ::Type{I}, c::I, val) where {I <: Sector}
+    v[findindex(values(I), c)] = val
+    return v
+end
+_denseget(v::Vector, ::Type{I}, c::I) where {I <: Sector} = v[findindex(values(I), c)]
+function _densepairs(v::Vector, ::Type{I}) where {I <: Sector}
+    vals = values(I)
+    return (vals[i] => x for (i, x) in enumerate(v) if !isnothing(x))
+end
+_densekeys(v::Vector, ::Type{I}) where {I <: Sector} = (c for (c, _) in _densepairs(v, I))
+
+# fallbacks to catch SectorVector/SectorDict, even for NTuple sectorstoragetype
+_denseget(v, ::Type{I}, c::I) where {I <: Sector} = get(v, c, nothing)
+_densekeys(v, ::Type{I}) where {I <: Sector} = keys(v)
+_densepairs(v, ::Type{I}) where {I <: Sector} = pairs(v)
+
+# builds either a dense Vector or SectorDict based on sectorstoragetype
+# mapping each (c, v) pair's sector c to f(c, v)
+# so every `findtruncated` method shares one output-construction path
+# `pairsiter` are c => v pairs, can be c => nothing for NoTruncation/TruncationIntersection/TruncationUnion
+function _builddensemap(f, ::Type{D}, ::Type{I}, pairsiter, ::Type{V}) where {D <: Tuple, I <: Sector, V}
+    d = _densenew(I, V)
+    for (c, v) in pairsiter
+        _denseset!(d, I, c, f(c, v))
+    end
+    return d
+end
+function _builddensemap(f, ::Type{D}, ::Type{I}, pairsiter, ::Type{V}) where {D <: SectorDict, I <: Sector, V}
+    return SectorDict(c => f(c, v) for (c, v) in pairsiter) # V unused
+end
+
 function truncate_space(V::ElementarySpace, inds)
     @assert !isdual(V)
-    return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in pairs(inds))
+    I = sectortype(V)
+    @assert I == Trivial
+    return spacetype(V)(c => _blocklength(dim(V, c), ind) for (c, ind) in _densepairs(inds, I))
 end
 function truncate_space(V::GradedSpace{I, NTuple{N, Int}}, inds) where {I <: Sector, N}
     @assert !isdual(V)
     vals = values(I)
     newdims = zeros(Int, N)
-    for (c, ind) in pairs(inds)
+    for (c, ind) in _densepairs(inds, I)
         d = dim(V, c)
         n_write = findindex(vals, c)
         newdims[n_write] = _blocklength(d, ind)
@@ -65,27 +105,30 @@ function truncate_space(V::GradedSpace{I, <:SectorDict}, inds) where {I <: Secto
 end
 
 function truncate_domain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds)
+    Isec = sectortype(tdst)
     for (c, b) in blocks(tdst)
-        I = get(inds, c, nothing)
-        @assert !isnothing(I)
+        I = _denseget(inds, Isec, c)
+        @assert !isnothing(I) # kept for safety, but should be guaranteed by _densepairs
         b′ = block(tsrc, c)
         b .= view(b′, :, I)
     end
     return tdst
 end
 function truncate_codomain!(tdst::AbstractTensorMap, tsrc::AbstractTensorMap, inds)
+    Isec = sectortype(tdst)
     for (c, b) in blocks(tdst)
-        I = get(inds, c, nothing)
-        @assert !isnothing(I)
+        I = _denseget(inds, Isec, c)
+        @assert !isnothing(I) # kept for safety, but should be guaranteed by _densepairs
         b′ = block(tsrc, c)
         b .= view(b′, I, :)
     end
     return tdst
 end
 function truncate_diagonal!(Ddst::DiagonalTensorMap, Dsrc::DiagonalTensorMap, inds)
+    Isec = sectortype(Ddst)
     for (c, b) in blocks(Ddst)
-        I = get(inds, c, nothing)
-        @assert !isnothing(I)
+        I = _denseget(inds, Isec, c)
+        @assert !isnothing(I) # kept for safety, but should be guaranteed by _densepairs
         diagview(b) .= view(diagview(block(Dsrc, c)), I)
     end
     return Ddst
@@ -140,7 +183,10 @@ end
 function MAK.truncate(
         ::typeof(left_null!), (U, S)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation
     )
-    ind = SectorDict(c => (size(b, 2) + 1):size(b, 1) for (c, b) in blocks(S))
+    I = sectortype(S)
+    ind = _builddensemap(sectorstoragetype(I), I, blocks(S), UnitRange{Int}) do _, b
+        (size(b, 2) + 1):size(b, 1)
+    end
     V_truncated = truncate_space(space(S, 1), ind)
     Ũ = similar(U, codomain(U) ← V_truncated)
     truncate_domain!(Ũ, U, ind)
@@ -149,7 +195,10 @@ end
 function MAK.truncate(
         ::typeof(right_null!), (S, Vᴴ)::NTuple{2, AbstractTensorMap}, strategy::NoTruncation
     )
-    ind = SectorDict(c => (size(b, 1) + 1):size(b, 2) for (c, b) in blocks(S))
+    I = sectortype(S)
+    ind = _builddensemap(sectorstoragetype(I), I, blocks(S), UnitRange{Int}) do _, b
+        (size(b, 1) + 1):size(b, 2)
+    end
     V_truncated = truncate_space(dual(space(S, 2)), ind)
     Ṽᴴ = similar(Vᴴ, V_truncated ← domain(Vᴴ))
     truncate_codomain!(Ṽᴴ, Vᴴ, ind)
@@ -186,7 +235,10 @@ function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationStrateg
 end
 
 function MAK.findtruncated(values::SectorVector, ::NoTruncation)
-    return SectorDict(c => Colon() for c in keys(values))
+    I = sectortype(values)
+    return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in keys(values)), Colon) do _, _
+        Colon()
+    end
 end
 
 # Need to select the first k values here after sorting across blocks, weighted by quantum dimension
@@ -232,18 +284,29 @@ MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByOrder) =
     MAK.findtruncated(values, strategy)
 
 function MAK.findtruncated(values::SectorVector, strategy::TruncationByFilter)
-    return SectorDict(c => findall(strategy.filter, d) for (c, d) in pairs(values))
+    I = sectortype(values)
+    return _builddensemap(sectorstoragetype(I), I, pairs(values), Vector{Int}) do _, v
+        findall(strategy.filter, v)
+    end
 end
 
 function MAK.findtruncated(values::SectorVector, strategy::TruncationByValue)
+    I = sectortype(values)
     atol = rtol_to_atol(values, strategy.p, strategy.atol, strategy.rtol)
     strategy′ = trunctol(; atol, strategy.by, strategy.keep_below)
-    return SectorDict(c => MAK.findtruncated(d, strategy′) for (c, d) in pairs(values))
+    V = Base.promote_op(MAK.findtruncated, valtype(values), typeof(strategy′))
+    return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do _, v
+        MAK.findtruncated(v, strategy′)
+    end
 end
 function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByValue)
+    I = sectortype(values)
     atol = rtol_to_atol(values, strategy.p, strategy.atol, strategy.rtol)
     strategy′ = trunctol(; atol, strategy.by, strategy.keep_below)
-    return SectorDict(c => MAK.findtruncated_svd(d, strategy′) for (c, d) in pairs(values))
+    V = Base.promote_op(MAK.findtruncated_svd, valtype(values), typeof(strategy′))
+    return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do _, v
+        MAK.findtruncated_svd(v, strategy′)
+    end
 end
 
 # Need to select the first k values here after sorting by error across blocks,
@@ -285,14 +348,24 @@ MAK.findtruncated_svd(values::SectorVector, strategy::TruncationByError) =
     MAK.findtruncated(values, strategy)
 
 function MAK.findtruncated(values::SectorVector, strategy::TruncationSpace)
-    sectortype(values) == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values"))
+    I = sectortype(values)
+    I == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values"))
     blockstrategy(c) = truncrank(dim(strategy.space, c); strategy.by, strategy.rev)
-    return SectorDict(c => MAK.findtruncated(d, blockstrategy(c)) for (c, d) in pairs(values))
+    Vstrategy = Base.promote_op(blockstrategy, I)
+    V = Base.promote_op(MAK.findtruncated, valtype(values), Vstrategy)
+    return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do c, v
+        MAK.findtruncated(v, blockstrategy(c))
+    end
 end
 function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationSpace)
-    sectortype(values) == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values"))
+    I = sectortype(values)
+    I == sectortype(strategy) || throw(SectorMismatch("sectortype of truncation strategy does not match values"))
     blockstrategy(c) = truncrank(dim(strategy.space, c); strategy.by, strategy.rev)
-    return SectorDict(c => MAK.findtruncated_svd(d, blockstrategy(c)) for (c, d) in pairs(values))
+    Vstrategy = Base.promote_op(blockstrategy, I)
+    V = Base.promote_op(MAK.findtruncated_svd, valtype(values), Vstrategy)
+    return _builddensemap(sectorstoragetype(I), I, pairs(values), V) do c, v
+        MAK.findtruncated_svd(v, blockstrategy(c))
+    end
 end
 
 # The implementations below assume that the `SectorDict` always contains an entry for every block sector
@@ -300,40 +373,40 @@ end
 # This is always the case in the implementations above.
 
 function MAK.findtruncated(values::SectorVector, strategy::TruncationIntersection)
+    I = sectortype(values)
     inds = map(Base.Fix1(MAK.findtruncated, values), strategy.components)
-    @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now"
-    sectors = keys(first(inds))
-    vals = map(keys(first(inds))) do c
-        mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_intersect, inds)
+    @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now"
+    sectors = collect(_densekeys(first(inds), I))
+    return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _
+        mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_intersect, inds)
     end
-    return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals)
 end
 function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationIntersection)
+    I = sectortype(values)
     inds = map(Base.Fix1(MAK.findtruncated_svd, values), strategy.components)
-    @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now"
-    sectors = keys(first(inds))
-    vals = map(keys(first(inds))) do c
-        mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_intersect, inds)
+    @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now"
+    sectors = collect(_densekeys(first(inds), I))
+    return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _
+        mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_intersect, inds)
     end
-    return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals)
 end
 function MAK.findtruncated(values::SectorVector, strategy::TruncationUnion)
+    I = sectortype(values)
     inds = map(Base.Fix1(MAK.findtruncated, values), strategy.components)
-    @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now"
-    sectors = keys(first(inds))
-    vals = map(keys(first(inds))) do c
-        mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_union, inds)
+    @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now"
+    sectors = collect(_densekeys(first(inds), I))
+    return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _
+        mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_union, inds)
     end
-    return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals)
 end
 function MAK.findtruncated_svd(values::SectorVector, strategy::TruncationUnion)
+    I = sectortype(values)
     inds = map(Base.Fix1(MAK.findtruncated_svd, values), strategy.components)
-    @assert TensorKit._allequal(keys, inds) "missing blocks are not supported right now"
-    sectors = keys(first(inds))
-    vals = map(keys(first(inds))) do c
-        mapreduce(Base.Fix2(getindex, c), MatrixAlgebraKit._ind_union, inds)
+    @assert TensorKit._allequal(v -> collect(_densekeys(v, I)), inds) "missing blocks are not supported right now"
+    sectors = collect(_densekeys(first(inds), I))
+    return _builddensemap(sectorstoragetype(I), I, (c => nothing for c in sectors), Any) do c, _
+        mapreduce(v -> _denseget(v, I, c), MatrixAlgebraKit._ind_union, inds)
     end
-    return SectorDict{eltype(sectors), eltype(vals)}(sectors, vals)
 end
 
 # Truncation error
@@ -341,7 +414,8 @@ end
 MAK.truncation_error(values::SectorVector, ind) = MAK.truncation_error!(copy(values), ind)
 
 function MAK.truncation_error!(values::SectorVector, ind)
-    for (c, ind_c) in pairs(ind)
+    Isec = sectortype(values)
+    for (c, ind_c) in _densepairs(ind, Isec)
         v = values[c]
         v[ind_c] .= zero(eltype(v))
     end

From d42aa28bea2e3a3b15c339112e21011d92d3c768 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 28 Aug 2026 18:31:07 +0200
Subject: [PATCH 18/21] put `_sortedmerge` in `Base.mergewith` and use where
 possible

---
 src/auxiliary/dicts.jl    | 67 ++++++++++++++++++++++-----------------
 src/spaces/gradedspace.jl | 12 +++----
 2 files changed, 42 insertions(+), 37 deletions(-)

diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl
index 17297177e..c84a9e1e5 100644
--- a/src/auxiliary/dicts.jl
+++ b/src/auxiliary/dicts.jl
@@ -186,17 +186,19 @@ function Base.:(==)(d1::SortedVectorDict, d2::SortedVectorDict)
     return true
 end
 
-# merge over two SORTED vector pairs representing keys and values
-# - combine(v1,v2): value for a key present in both operands
-# - only1(v1) / only2(v2): value for a key present in only one operand;
-# pass `nothing` to drop such keys entirely (e.g. for an intersection)
-# zero results are dropped (either from `combine` or `only1`/`only2`), matching how GradedSpace never stores an explicit zero dimension
-# k1 and k2 originate from GradedSpace.dims.keys, which are guaranteed to be sorted
-function _sortedmerge(k1::Vector{I}, v1::Vector{Int}, k2::Vector{I}, v2::Vector{Int}, combine, only1, only2) where {I}
+# merge over two SectorDicts
+# the intersect case for infimum is kind of tricky, so there's an extra bool
+# to indicate keeping keys that are only present in one of the two dicts
+# zero results are dropped, matching how GradedSpace never stores an explicit zero dimension
+function _sortedmerge(
+        combine, ::Val{keepunique}, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}
+    ) where {keepunique, K, V}
+    k1, v1 = d1.keys, d1.values
+    k2, v2 = d2.keys, d2.values
     n1, n2 = length(k1), length(k2)
-    ks, vs = Vector{I}(), Vector{Int}()
-    sizehint!(ks, n1 + n2)
-    sizehint!(vs, n1 + n2)
+    ks, vs = Vector{K}(), Vector{V}()
+    sizehint!(ks, keepunique ? n1 + n2 : min(n1, n2))
+    sizehint!(vs, keepunique ? n1 + n2 : min(n1, n2))
     i, j = 1, 1
     @inbounds while i <= n1 && j <= n2
         if k1[i] == k2[j]
@@ -208,33 +210,40 @@ function _sortedmerge(k1::Vector{I}, v1::Vector{Int}, k2::Vector{I}, v2::Vector{
             i += 1
             j += 1
         elseif k1[i] < k2[j]
-            _mergeonly!(ks, vs, k1[i], v1[i], only1)
+            if keepunique
+                push!(ks, k1[i])
+                push!(vs, v1[i])
+            end
             i += 1
         else
-            _mergeonly!(ks, vs, k2[j], v2[j], only2)
+            if keepunique
+                push!(ks, k2[j])
+                push!(vs, v2[j])
+            end
             j += 1
         end
     end
-    @inbounds while i <= n1
-        _mergeonly!(ks, vs, k1[i], v1[i], only1)
-        i += 1
-    end
-    @inbounds while j <= n2
-        _mergeonly!(ks, vs, k2[j], v2[j], only2)
-        j += 1
-    end
-    return ks, vs
-end
-@inline _mergeonly!(ks, vs, k, v, ::Nothing) = nothing
-@inline function _mergeonly!(ks, vs, k, v, f)
-    d = f(v)
-    if !iszero(d)
-        push!(ks, k)
-        push!(vs, d)
+    if keepunique
+        @inbounds while i <= n1
+            push!(ks, k1[i])
+            push!(vs, v1[i])
+            i += 1
+        end
+        @inbounds while j <= n2
+            push!(ks, k2[j])
+            push!(vs, v2[j])
+            j += 1
+        end
     end
-    return nothing
+    return SortedVectorDict{K, V}(ks, vs)
 end
 
+Base.mergewith(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V} =
+    _sortedmerge(combine, Val(true), d1, d2)
+
+_sortedintersect(combine, d1::SortedVectorDict{K, V}, d2::SortedVectorDict{K, V}) where {K, V} =
+    _sortedmerge(combine, Val(false), d1, d2)
+
 """
     Hashed(value, hashfunction = Base.hash, isequal = Base.isequal)
 
diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index b7b35c8f9..25bce0ec5 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -138,8 +138,7 @@ zerospace(S::Type{<:GradedSpace}) = S()
 function ⊕(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
     dual1 = isdual(V₁)
     dual1 == isdual(V₂) || throw(SpaceMismatch("Direct sum of a vector space and a dual space does not exist"))
-    ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, +, identity, identity)
-    return typeof(V₁)(SectorDict{I, Int}(ks, vs), dual1)
+    return typeof(V₁)(mergewith(+, V₁.dims, V₂.dims), dual1)
 end
 function ⊕(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector}
     dual1 = isdual(V₁)
@@ -157,8 +156,7 @@ end
 function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) where {I <: Sector}
     dualV = isdual(V)
     V ≿ W && dualV == isdual(W) || throw(SpaceMismatch("$(W) is not a subspace of $(V)"))
-    ks, vs = _sortedmerge(V.dims.keys, V.dims.values, W.dims.keys, W.dims.values, -, identity, nothing)
-    return typeof(V)(SectorDict{I, Int}(ks, vs), dualV)
+    return typeof(V)(mergewith(-, V.dims, W.dims), dualV)
 end
 
 function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
@@ -216,8 +214,7 @@ end
 function infimum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
     Visdual = isdual(V₁)
     Visdual == isdual(V₂) || throw(SpaceMismatch("Infimum of space and dual space does not exist"))
-    ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, min, nothing, nothing)
-    return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual)
+    return typeof(V₁)(_sortedintersect(min, V₁.dims, V₂.dims), Visdual)
 end
 function supremum(V₁::GradedSpace{I, <:Tuple}, V₂::GradedSpace{I, <:Tuple}) where {I <: Sector}
     Visdual = isdual(V₁)
@@ -228,8 +225,7 @@ end
 function supremum(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
     Visdual = isdual(V₁)
     Visdual == isdual(V₂) || throw(SpaceMismatch("Supremum of space and dual space does not exist"))
-    ks, vs = _sortedmerge(V₁.dims.keys, V₁.dims.values, V₂.dims.keys, V₂.dims.values, max, identity, identity)
-    return typeof(V₁)(SectorDict{I, Int}(ks, vs), Visdual)
+    return typeof(V₁)(mergewith(max, V₁.dims, V₂.dims), Visdual)
 end
 
 hassector(V::GradedSpace{I}, s::I) where {I <: Sector} = dim(V, s) != 0

From f8b0288cfd4d7bb11e51348f9f7a998530b62631 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Tue, 1 Sep 2026 18:04:01 +0200
Subject: [PATCH 19/21] introduce and use `blockdims`

---
 src/spaces/gradedspace.jl | 51 +++++++++++++++++++++++----------------
 1 file changed, 30 insertions(+), 21 deletions(-)

diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index 25bce0ec5..e4bb53b21 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -31,7 +31,7 @@ sectortype(::Type{<:GradedSpace{I}}) where {I <: Sector} = I
 
 function GradedSpace{I, NTuple{N, Int}}(dims; dual::Bool = false) where {I, N}
     d = zeros(Int, N)
-    isset = falses(N)
+    isset = falses(N) # see if this is still needed if we're restricting to small N
     for (c, dc) in dims
         k = convert(I, c)
         i = findindex(values(I), k)
@@ -111,9 +111,9 @@ end
 Base.axes(V::GradedSpace) = Base.OneTo(dim(V))
 function Base.axes(V::GradedSpace{I}, c::I) where {I <: Sector}
     offset = 0
-    for c′ in sectors(V)
+    for (c′, d′) in blockdims(V)
         c′ == c && break
-        offset += dim(c′) * dim(V, c′)
+        offset += dim(c′) * d′
     end
     return (offset + 1):(offset + dim(c) * dim(V, c))
 end
@@ -124,9 +124,9 @@ isdual(V::GradedSpace) = V.dual
 isconj(V::GradedSpace) = isdual(V)
 function flip(V::GradedSpace{I}) where {I <: Sector}
     return if isdual(V)
-        typeof(V)(c => dim(V, c) for c in sectors(V))
+        typeof(V)(blockdims(V))
     else
-        typeof(V)(dual(c) => dim(V, c) for c in sectors(V))'
+        typeof(V)(dual(c) => d for (c, d) in blockdims(V))'
     end
 end
 
@@ -160,20 +160,11 @@ function ⊖(V::GradedSpace{I, <:SectorDict}, W::GradedSpace{I, <:SectorDict}) w
 end
 
 function fuse(V₁::GradedSpace{I, <:SectorDict}, V₂::GradedSpace{I, <:SectorDict}) where {I <: Sector}
-    dual1, dual2 = isdual(V₁), isdual(V₂)
-    acc = Dict{I, Int}() # SectorDict `get` within the double for loop accumulates O(N^2) ` findindex` calls -> sort afterwards
-    k1, k2 = V₁.dims.keys, V₂.dims.keys
-    v1, v2 = V₁.dims.values, V₂.dims.values
-    @inbounds for na in eachindex(k1)
-        a₀, da = k1[na], v1[na]
-        a = dual1 ? dual(a₀) : a₀
-        for nb in eachindex(k2)
-            b₀, db = k2[nb], v2[nb]
-            b = dual2 ? dual(b₀) : b₀
-            dab = da * db
-            for c in a ⊗ b
-                acc[c] = get(acc, c, 0) + Nsymbol(a, b, c) * dab
-            end
+    acc = Dict{I, Int}() # SectorDict `get` within the double for loop accumulates O(N^2) `findindex` calls -> sort afterwards
+    for (a, da) in blockdims(V₁), (b, db) in blockdims(V₂)
+        dab = da * db
+        for c in a ⊗ b
+            acc[c] = get(acc, c, 0) + Nsymbol(a, b, c) * dab
         end
     end
     ks0 = collect(keys(acc))
@@ -238,6 +229,24 @@ function sectors(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N}
     end
 end
 
+"""
+    blockdims(V::GradedSpace)
+
+Return an iterator over the non-zero blocks of the graded space `V`.
+These blocks contain the `Sector`s and their corresponding degeneracy,
+i.e. the number of times the sector appears in the direct sum decomposition of `V`.
+"""
+function blockdims(V::GradedSpace{I, <:AbstractDict}) where {I <: Sector}
+    return ((isdual(V) ? dual(c) : c) => d for (c, d) in V.dims)
+end
+function blockdims(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N}
+    vals = values(I)
+    return (
+        (isdual(V) ? dual(vals[n]) : vals[n]) => V.dims[n]
+            for n in 1:N if !iszero(V.dims[n])
+    )
+end
+
 Base.hash(V::GradedSpace, h::UInt) = hash(V.dual, hash(V.dims, h))
 function Base.:(==)(V₁::GradedSpace, V₂::GradedSpace)
     return sectortype(V₁) == sectortype(V₂) && (V₁.dims == V₂.dims) && V₁.dual == V₂.dual
@@ -270,7 +279,7 @@ function Base.show(io::IO, V::GradedSpace)
         cls = ")"
     end
 
-    v = [c => dim(V, c) for c in sectors(V)]
+    v = collect(blockdims(V))
 
     # logic stolen from Base.show_vector
     limited = get(io, :limit, false)::Bool
@@ -301,7 +310,7 @@ function Base.show(io::IO, ::MIME"text/plain", V::GradedSpace)
     # print detailed sector information - hijack Base.Vector printing
     print(io, ":\n")
     isdual(V) && (V = dual(V))
-    print_data = [c => dim(V, c) for c in sectors(V)]
+    print_data = collect(blockdims(V))
     ioc = IOContext(io, :typeinfo => eltype(print_data))
     Base.print_matrix(ioc, print_data)
 

From 2fbf6feb2b8d7bd294ef80d1700b2591686f5be5 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Wed, 2 Sep 2026 17:06:56 +0200
Subject: [PATCH 20/21] introduce the ntuple cutoff

---
 src/spaces/gradedspace.jl | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index e4bb53b21..c2f6758f9 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -330,19 +330,22 @@ Base.getindex(::SpaceTable) = ComplexSpace
 Base.getindex(::SpaceTable, ::Type{Trivial}) = ComplexSpace
 Base.getindex(::SpaceTable, I::Type{<:Sector}) = GradedSpace{I, sectorstoragetype(I)}
 
+# based on Julia tuple unrolling range
+const _ntuple_storage_threshold = 32
+
 """
     sectorstoragetype(I::Type{<:Sector}) -> Type
 
 The storage type `D` used for the `dims` field of `GradedSpace{I, D}`.
-This is `NTuple{N,Int}` with `N = length(values(I))` if `I` has a finite, known length,
-or `SectorDict{I,Int}` otherwise.
+This is `NTuple{N,Int}` with `N = length(values(I))` if `I` has a finite, known length
+of at most `$_ntuple_storage_threshold`, or `SectorDict{I,Int}` otherwise.
 """
 Base.@assume_effects :foldable function sectorstoragetype(::Type{I}) where {I <: Sector}
     if Base.IteratorSize(values(I)) isa Union{HasLength, HasShape}
-        return NTuple{length(values(I)), Int}
-    else
-        return SectorDict{I, Int}
+        N = length(values(I))
+        N <= _ntuple_storage_threshold && return NTuple{N, Int}
     end
+    return SectorDict{I, Int}
 end
 
 Base.getindex(::ComplexNumbers, I::Type{<:Sector}) = Vect[I]

From e5877dce1be4fd3023b95c2d236e0c4ca42d3eb1 Mon Sep 17 00:00:00 2001
From: Boris De Vos <boris.devos@ugent.be>
Date: Fri, 4 Sep 2026 10:31:22 +0200
Subject: [PATCH 21/21] apply code suggestions

---
 src/auxiliary/dicts.jl    | 21 ++++++---------------
 src/spaces/gradedspace.jl | 13 ++-----------
 2 files changed, 8 insertions(+), 26 deletions(-)

diff --git a/src/auxiliary/dicts.jl b/src/auxiliary/dicts.jl
index c84a9e1e5..dbce97087 100644
--- a/src/auxiliary/dicts.jl
+++ b/src/auxiliary/dicts.jl
@@ -89,21 +89,12 @@ end
 Base.empty(::SortedVectorDict, ::Type{K}, ::Type{V}) where {K, V} = SortedVectorDict{K, V}()
 Base.empty!(d::SortedVectorDict) = (empty!(d.keys); empty!(d.values); return d)
 
-_searchsortedfirst(v::Vector, k) = searchsortedfirst(v, k)
-# function _searchsortedfirst(v::Vector, k)
-#     i = 1
-#     @inbounds while i <= length(v) && isless(v[i], k)
-#         i += 1
-#     end
-#     return i
-# end
-
 function Base.delete!(d::SortedVectorDict{K}, k) where {K}
     key = convert(K, k)
     if !isequal(k, key)
         return d
     end
-    i = _searchsortedfirst(d.keys, key)
+    i = searchsortedfirst(d.keys, key)
     if i <= length(d) && isequal(d.keys[i], key)
         deleteat!(d.keys, i)
         deleteat!(d.values, i)
@@ -118,7 +109,7 @@ function Base.haskey(d::SortedVectorDict{K}, k) where {K}
     if !isequal(k, key)
         return false
     end
-    i = _searchsortedfirst(d.keys, key)
+    i = searchsortedfirst(d.keys, key)
     return (i <= length(d) && isequal(d.keys[i], key))
 end
 function Base.getindex(d::SortedVectorDict{K}, k) where {K}
@@ -126,7 +117,7 @@ function Base.getindex(d::SortedVectorDict{K}, k) where {K}
     if !isequal(k, key)
         throw(KeyError(k))
     end
-    i = _searchsortedfirst(d.keys, key)
+    i = searchsortedfirst(d.keys, key)
     @inbounds if (i <= length(d) && isequal(d.keys[i], key))
         return d.values[i]
     else
@@ -138,7 +129,7 @@ function Base.setindex!(d::SortedVectorDict{K}, v, k) where {K}
     if !isequal(k, key)
         throw(ArgumentError("$k is not a valid key for type $K"))
     end
-    i = _searchsortedfirst(d.keys, key)
+    i = searchsortedfirst(d.keys, key)
     if i <= length(d) && isequal(d.keys[i], key)
         d.values[i] = v
     else
@@ -153,7 +144,7 @@ function Base.get(d::SortedVectorDict{K}, k, default) where {K}
     if !isequal(k, key)
         return default
     end
-    i = _searchsortedfirst(d.keys, key)
+    i = searchsortedfirst(d.keys, key)
     @inbounds begin
         return (i <= length(d) && isequal(d.keys[i], key)) ? d.values[i] : default
     end
@@ -163,7 +154,7 @@ function Base.get(f::Union{Function, Type}, d::SortedVectorDict{K}, k) where {K}
     if !isequal(k, key)
         return f()
     end
-    i = _searchsortedfirst(d.keys, key)
+    i = searchsortedfirst(d.keys, key)
     @inbounds begin
         return (i <= length(d) && isequal(d.keys[i], key)) ? d.values[i] : f()
     end
diff --git a/src/spaces/gradedspace.jl b/src/spaces/gradedspace.jl
index c2f6758f9..facb63811 100644
--- a/src/spaces/gradedspace.jl
+++ b/src/spaces/gradedspace.jl
@@ -89,18 +89,9 @@ GradedSpace(g::AbstractDict; dual::Bool = false) = GradedSpace(g...; dual = dual
 field(::Type{<:GradedSpace}) = ℂ
 InnerProductStyle(::Type{<:GradedSpace}) = EuclideanInnerProduct()
 
-function dim(V::GradedSpace{I, <:AbstractDict}) where {I <: Sector}
+function dim(V::GradedSpace{I}) where {I <: Sector}
     init = zero(dimscalartype(I))
-    return sum(((c, d),) -> dim(c) * d, V.dims; init)
-end
-function dim(V::GradedSpace{I, NTuple{N, Int}}) where {I <: Sector, N}
-    init = zero(dimscalartype(I))
-    D = init
-    vals = values(I)
-    @inbounds for (c, d) in zip(values(I), V.dims)
-        D += dim(c) * d
-    end
-    return D
+    return sum(((c, d),) -> dim(c) * d, blockdims(V); init)
 end
 function dim(V::GradedSpace{I, <:AbstractDict}, c::I) where {I <: Sector}
     return get(V.dims, isdual(V) ? dual(c) : c, 0)