From c38839667541fcf93b17d6a05e1a37a1be5b431a Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 15 Oct 2026 10:21:19 -0400 Subject: [PATCH 1/8] Closed-form Manhattan indexing Replace the recursive enumeration in `to_manhattan_index` and `manhattan_to_multidimensional_index` by an inclusion-exclusion count of the points below a given Manhattan distance, together with a binary search for the inverse mapping. This makes both mappings O(N 2^N) instead of growing polynomially with the distance, which speeds up `findindex` and `getindex` for `values(::ProductSector)`, and with it `isless` for sectors such as SU(N) irreps whose `findindex` relies on these functions. Co-Authored-By: Claude Opus 5.5 --- src/auxiliary.jl | 105 ++++++++++++++++++++++++++++------------------- 1 file changed, 62 insertions(+), 43 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 6ed53605..f2fee3b5 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -19,76 +19,95 @@ function _kron_promote(A₁, B₁, sz₁, sz₂) end # Manhattan based distance enumeration: I is supposed to be one-based index -# TODO: is there any way to make this faster? +# Points are zero-based coordinates `x` with `0 ≤ x[i] < sz[i]`, ordered first by their +# Manhattan distance `d = sum(x)` to the origin, and then lexicographically. -# forward mapping from multidimensional to single Manhattan index -@inline function num_manhattan_points(d::Int, sz::Dims{1}) - return Int(sz[1] > d) -end -@inline function num_manhattan_points(d::Int, sz::Dims{N}) where {N} - d == 0 && return 1 +# number of points with Manhattan distance smaller than `d`, computed in closed form through +# inclusion-exclusion over the dimensions where the upper bound `sz[i]` is violated +@inline function num_manhattan_points_below(d::Int, sz::Dims{N}) where {N} + d ≤ 0 && return 0 num = 0 - for i in 1:min(sz[1], d + 1) - num += num_manhattan_points(d - i + 1, Base.tail(sz)) + for mask in 0:((1 << N) - 1) + s = 0 + sign = 1 + for i in 1:N + if isodd(mask >> (i - 1)) + s += min(sz[i], d) # clamp to avoid overflow, terms with s ≥ d vanish anyway + sign = -sign + end + end + s < d && (num += sign * _binomial(d - 1 - s + N, Val(N))) end return num end +# binomial(n, N) for n ≥ N, unrolled for small static N +@inline function _binomial(n::Int, ::Val{N}) where {N} + r = 1 + for k in 1:N + r = (r * (n - N + k)) ÷ k + end + return r +end + +# number of points with Manhattan distance exactly `d` +@inline function num_manhattan_points(d::Int, sz::Dims) + return num_manhattan_points_below(d + 1, sz) - num_manhattan_points_below(d, sz) +end + +# forward mapping from multidimensional to single Manhattan index @inline localoffset(d, I::Tuple{Int}, sz::Tuple{Int}) = 0 @inline function localoffset(d, I::NTuple{N, Int}, sz::NTuple{N, Int}) where {N} - offset = 0 - for i in 1:(I[1] - 1) - offset += num_manhattan_points(d - i + 1, Base.tail(sz)) - end - offset += localoffset(d - I[1] + 1, Base.tail(I), Base.tail(sz)) - return offset + # number of points at distance `d` with a smaller first coordinate + sz′ = Base.tail(sz) + offset = num_manhattan_points_below(d + 1, sz′) - num_manhattan_points_below(d - I[1] + 2, sz′) + return offset + localoffset(d - I[1] + 1, Base.tail(I), sz′) end to_manhattan_index(::Tuple{}, ::Tuple{}) = 1 to_manhattan_index(I::Tuple{Int}, ::Tuple{Int}) = I[1] function to_manhattan_index(I::Dims{N}, sz::Dims{N}) where {N} d = sum(I) - N # Manhattan distance to origin for one-based indices - d == 0 && return 1 - index = 1 - # count all the points with smaller Manhatten distance - for k in 0:(d - 1) - index += num_manhattan_points(k, sz) - end - return index + localoffset(d, I, sz) + return 1 + num_manhattan_points_below(d, sz) + localoffset(d, I, sz) end # inverse mapping @inline invertlocaloffset(d, offset, sz::Tuple{Int}) = (d + 1,) @inline function invertlocaloffset(d, offset, sz::Tuple{Int, Int, Vararg{Int}}) - i₁ = 1 - while i₁ < sz[1] - jump = num_manhattan_points(d - i₁ + 1, Base.tail(sz)) - if offset < jump - break + # find the largest first coordinate `x₁` such that the number of points at distance `d` + # with a smaller first coordinate does not exceed `offset` + sz′ = Base.tail(sz) + nd = num_manhattan_points_below(d + 1, sz′) + lo, hi = 0, min(sz[1] - 1, d) + while lo < hi + mid = (lo + hi + 1) >> 1 + if nd - num_manhattan_points_below(d - mid + 1, sz′) ≤ offset + lo = mid + else + hi = mid - 1 end - offset -= jump - i₁ += 1 end - return (i₁, invertlocaloffset(d - i₁ + 1, offset, Base.tail(sz))...) + offset -= nd - num_manhattan_points_below(d - lo + 1, sz′) + return (lo + 1, invertlocaloffset(d - lo, offset, sz′)...) end manhattan_to_multidimensional_index(index::Int, ::Dims{0}) = () manhattan_to_multidimensional_index(index::Int, ::Dims{1}) = (index,) function manhattan_to_multidimensional_index(index::Int, sz::Dims{N}) where {N} - # find the layer of the Manhattan distance - index == 1 && return ntuple(one, Val(N)) - offset = 1 - d = 1 - while true - currentlayer = num_manhattan_points(d, sz) - if index <= offset + currentlayer - break + # find the layer of the Manhattan distance, i.e. the smallest `d` such that there are at + # least `index` points with distance `≤ d`, by exponential followed by binary search + lo, hi = 0, 1 + while num_manhattan_points_below(hi + 1, sz) < index + lo, hi = hi + 1, 2hi + end + while lo < hi + mid = (lo + hi) >> 1 + if num_manhattan_points_below(mid + 1, sz) < index + lo = mid + 1 + else + hi = mid end - d += 1 - offset += currentlayer end - # find the position within the layer - index -= offset - return invertlocaloffset(d, index - 1, sz) + return invertlocaloffset(lo, index - 1 - num_manhattan_points_below(lo, sz), sz) end From 7880fc87d53d1a6fab47b504d8b9c9db2d74fe5f Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 15 Oct 2026 10:28:19 -0400 Subject: [PATCH 2/8] Stateful iteration over ProductSector values Iterating `values(::ProductSector)` previously called `getindex` at every step, recomputing the multidimensional index from scratch. Instead keep the current index in the iterator state and advance it to its Manhattan-order successor with `next_manhattan_index`, which is O(N) per step. Co-Authored-By: Claude Opus 5.5 --- src/auxiliary.jl | 27 +++++++++++++++++++++++++++ src/product.jl | 11 ++++++++--- 2 files changed, 35 insertions(+), 3 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index f2fee3b5..7dbe7eee 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -91,6 +91,33 @@ end return (lo + 1, invertlocaloffset(d - lo, offset, sz′)...) end +# successor of a multidimensional index in Manhattan order, assuming it exists +function next_manhattan_index(I::Dims{N}, sz::Dims{N}) where {N} + # find the last position that can be incremented while lowering the distance of the + # positions after it, and otherwise move to the next layer + r = 0 # Manhattan distance of the positions after k + k = N + while k > 0 + r > 0 && I[k] < sz[k] && break + r += I[k] - 1 + k -= 1 + end + if k > 0 + I = Base.setindex(I, I[k] + 1, k) + r -= 1 + else + r += 1 + end + # the smallest index with distance `r` for the positions after k is obtained by + # filling them up from the back + for j in N:-1:(k + 1) + x = min(r, sz[j] - 1) + I = Base.setindex(I, x + 1, j) + r -= x + end + return I +end + manhattan_to_multidimensional_index(index::Int, ::Dims{0}) = () manhattan_to_multidimensional_index(index::Int, ::Dims{1}) = (index,) function manhattan_to_multidimensional_index(index::Int, sz::Dims{N}) where {N} diff --git a/src/product.jl b/src/product.jl index b0ed5554..f89cc1e8 100644 --- a/src/product.jl +++ b/src/product.jl @@ -47,9 +47,14 @@ function findindex(P::SectorValues{I}, c::I) where {I <: ProductSector} return to_manhattan_index(findindex.(values.(_sectors(I)), Tuple(c)), _size(P)) end -function Base.iterate(P::SectorValues{I}, i = 1) where {I <: ProductSector} - Base.IteratorSize(P) != Base.IsInfinite() && i > length(P) && return nothing - return getindex(P, i), i + 1 +function Base.iterate(P::SectorValues{I}) where {I <: ProductSector} + inds = map(Returns(1), _size(P)) + return I(getindex.(values.(_sectors(I)), inds)), (1, inds) +end +function Base.iterate(P::SectorValues{I}, (i, inds)) where {I <: ProductSector} + Base.IteratorSize(P) != Base.IsInfinite() && i ≥ length(P) && return nothing + inds′ = next_manhattan_index(inds, _size(P)) + return I(getindex.(values.(_sectors(I)), inds′)), (i + 1, inds′) end function Base.convert(::Type{ProductSector{T}}, t::Tuple) where {T <: SectorTuple} From 91e934b9912dd568cc30af5926ac052de6fd6ff0 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 15 Oct 2026 10:33:04 -0400 Subject: [PATCH 3/8] Use a for loop in `next_manhattan_index` Co-Authored-By: Claude Opus 5.5 --- src/auxiliary.jl | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 7dbe7eee..ec94f919 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -96,11 +96,13 @@ function next_manhattan_index(I::Dims{N}, sz::Dims{N}) where {N} # find the last position that can be incremented while lowering the distance of the # positions after it, and otherwise move to the next layer r = 0 # Manhattan distance of the positions after k - k = N - while k > 0 - r > 0 && I[k] < sz[k] && break - r += I[k] - 1 - k -= 1 + k = 0 + for j in N:-1:1 + if r > 0 && I[j] < sz[j] + k = j + break + end + r += I[j] - 1 end if k > 0 I = Base.setindex(I, I[k] + 1, k) From 2fda13fe7d9a0b1ce8f9ba331fa6706d88243511 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 15 Oct 2026 10:33:04 -0400 Subject: [PATCH 4/8] Add `@inbounds` annotations to `next_manhattan_index` Co-Authored-By: Claude Opus 5.5 --- src/auxiliary.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index ec94f919..1c25153a 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -97,7 +97,7 @@ function next_manhattan_index(I::Dims{N}, sz::Dims{N}) where {N} # positions after it, and otherwise move to the next layer r = 0 # Manhattan distance of the positions after k k = 0 - for j in N:-1:1 + @inbounds for j in N:-1:1 if r > 0 && I[j] < sz[j] k = j break @@ -105,14 +105,14 @@ function next_manhattan_index(I::Dims{N}, sz::Dims{N}) where {N} r += I[j] - 1 end if k > 0 - I = Base.setindex(I, I[k] + 1, k) + I = @inbounds Base.setindex(I, I[k] + 1, k) r -= 1 else r += 1 end # the smallest index with distance `r` for the positions after k is obtained by # filling them up from the back - for j in N:-1:(k + 1) + @inbounds for j in N:-1:(k + 1) x = min(r, sz[j] - 1) I = Base.setindex(I, x + 1, j) r -= x From 85f7503fdf78353b58209aecef7d04bf7d121341 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 15 Oct 2026 11:00:19 -0400 Subject: [PATCH 5/8] Drop unnecessary `@inline` annotations in Manhattan indexing Benchmarks show that only `num_manhattan_points_below`, both `localoffset` methods and the recursive `invertlocaloffset` method benefit from the annotation; the recursive base case of `localoffset` is not inlined otherwise. Co-Authored-By: Claude Opus 5.5 --- src/auxiliary.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 1c25153a..9221fa1e 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -42,7 +42,7 @@ end end # binomial(n, N) for n ≥ N, unrolled for small static N -@inline function _binomial(n::Int, ::Val{N}) where {N} +function _binomial(n::Int, ::Val{N}) where {N} r = 1 for k in 1:N r = (r * (n - N + k)) ÷ k @@ -51,7 +51,7 @@ end end # number of points with Manhattan distance exactly `d` -@inline function num_manhattan_points(d::Int, sz::Dims) +function num_manhattan_points(d::Int, sz::Dims) return num_manhattan_points_below(d + 1, sz) - num_manhattan_points_below(d, sz) end @@ -72,7 +72,7 @@ function to_manhattan_index(I::Dims{N}, sz::Dims{N}) where {N} end # inverse mapping -@inline invertlocaloffset(d, offset, sz::Tuple{Int}) = (d + 1,) +invertlocaloffset(d, offset, sz::Tuple{Int}) = (d + 1,) @inline function invertlocaloffset(d, offset, sz::Tuple{Int, Int, Vararg{Int}}) # find the largest first coordinate `x₁` such that the number of points at distance `d` # with a smaller first coordinate does not exceed `offset` From 352d68489100bcfce872a3322b24b80c908f2582 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 15 Oct 2026 11:00:38 -0400 Subject: [PATCH 6/8] Explain the inclusion-exclusion count in `num_manhattan_points_below` Co-Authored-By: Claude Opus 5.5 --- src/auxiliary.jl | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 9221fa1e..dcba6f87 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -22,8 +22,16 @@ end # Points are zero-based coordinates `x` with `0 ≤ x[i] < sz[i]`, ordered first by their # Manhattan distance `d = sum(x)` to the origin, and then lexicographically. -# number of points with Manhattan distance smaller than `d`, computed in closed form through -# inclusion-exclusion over the dimensions where the upper bound `sz[i]` is violated +# Number of points `x` with `0 ≤ x[i] < sz[i]` and Manhattan distance `sum(x) < d`. +# +# Without the upper bounds, the number of points `x ≥ 0` with `sum(x) ≤ d - 1` is +# `binomial(d - 1 + N, N)` (stars and bars). The points that violate the upper bound in +# (at least) a subset `S` of the dimensions, i.e. with `x[i] ≥ sz[i]` for `i ∈ S`, are +# counted in the same way after shifting `x[i] → x[i] - sz[i]`, which gives +# `binomial(d - 1 - s + N, N)` with `s = sum(sz[i] for i in S)`. By inclusion-exclusion, +# the number of points that satisfy all bounds is then +# ∑_S (-1)^|S| binomial(d - 1 - s + N, N) +# where the subsets `S` are enumerated as the bits of `mask`, and terms with `s ≥ d` vanish. @inline function num_manhattan_points_below(d::Int, sz::Dims{N}) where {N} d ≤ 0 && return 0 num = 0 @@ -42,6 +50,7 @@ end end # binomial(n, N) for n ≥ N, unrolled for small static N +# (this is faster than a lookup table, since the divisions by constants are cheap) function _binomial(n::Int, ::Val{N}) where {N} r = 1 for k in 1:N From f0a3ad1312fc2e23c33f4c0c9c5fb1b51d886604 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 15 Oct 2026 11:00:38 -0400 Subject: [PATCH 7/8] Single-line definition of `num_manhattan_points` Co-Authored-By: Claude Opus 5.5 --- src/auxiliary.jl | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index dcba6f87..0c657f83 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -60,9 +60,7 @@ function _binomial(n::Int, ::Val{N}) where {N} end # number of points with Manhattan distance exactly `d` -function num_manhattan_points(d::Int, sz::Dims) - return num_manhattan_points_below(d + 1, sz) - num_manhattan_points_below(d, sz) -end +num_manhattan_points(d::Int, sz::Dims) = num_manhattan_points_below(d + 1, sz) - num_manhattan_points_below(d, sz) # forward mapping from multidimensional to single Manhattan index @inline localoffset(d, I::Tuple{Int}, sz::Tuple{Int}) = 0 From dc1a9f149447c281b1b4db1af186a4c6cf21c2a6 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 15 Oct 2026 11:25:13 -0400 Subject: [PATCH 8/8] Remove unused `num_manhattan_points` Co-Authored-By: Claude Opus 5.5 --- src/auxiliary.jl | 3 --- 1 file changed, 3 deletions(-) diff --git a/src/auxiliary.jl b/src/auxiliary.jl index 0c657f83..c90692cb 100644 --- a/src/auxiliary.jl +++ b/src/auxiliary.jl @@ -59,9 +59,6 @@ function _binomial(n::Int, ::Val{N}) where {N} return r end -# number of points with Manhattan distance exactly `d` -num_manhattan_points(d::Int, sz::Dims) = num_manhattan_points_below(d + 1, sz) - num_manhattan_points_below(d, sz) - # forward mapping from multidimensional to single Manhattan index @inline localoffset(d, I::Tuple{Int}, sz::Tuple{Int}) = 0 @inline function localoffset(d, I::NTuple{N, Int}, sz::NTuple{N, Int}) where {N}