From 83693efa63ee6d0ef893d7cb10487cb73327de2a Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 7 Oct 2026 10:12:21 -0400 Subject: [PATCH 1/7] Reuse stride orders when planning mapreduce loops --- ext/StridedGPUArraysExt.jl | 2 +- src/mapreduce.jl | 74 ++++++++++++++++++++++---------------- 2 files changed, 44 insertions(+), 32 deletions(-) diff --git a/ext/StridedGPUArraysExt.jl b/ext/StridedGPUArraysExt.jl index 7525b44..aadfef0 100644 --- a/ext/StridedGPUArraysExt.jl +++ b/ext/StridedGPUArraysExt.jl @@ -128,7 +128,7 @@ end function Strided._mapreduce_block!( f, op, initop, dims::Dims{N}, - strides, offsets, costs, + strides, offsets, costs, strideorders, arrays::Tuple{GPUStridedView{TO, N}, Vararg{GPUStridedView{<:Any, N}}} ) where {TO, N} diff --git a/src/mapreduce.jl b/src/mapreduce.jl index 48b3e66..77959aa 100644 --- a/src/mapreduce.jl +++ b/src/mapreduce.jl @@ -118,26 +118,6 @@ function _mapreducedim!( return arrays[1] end -# Cache-friendly reordering of the loop dimensions: put the most important -# (smallest-stride, most contiguous) dimension first and size-1 dimensions last. -# The importance of each dimension is modelled from the `indexorder` of every -# array's strides, with the output array (`strides[1]`) weighted by a factor 2. -# `strides` is the per-array tuple-of-tuples; returns the reordered `(dims, strides)`. -function _sortdims(dims, strides) - M = length(strides) - N = length(dims) - # ceil(Int, log2(M+2)) # to account for the fact that there are M arrays, where the first one is counted with a factor 2 - g = 8 * sizeof(Int) - leading_zeros(M + 1) - importance = 2 .* (1 .<< (g .* (N .- indexorder(strides[1])))) - for k in 2:M - importance = importance .+ (1 .<< (g .* (N .- indexorder(strides[k])))) - end - - importance = importance .* (dims .> 1) # put dims 1 at the back - p = TupleTools.sortperm(importance; rev = true) - return TupleTools.getindices(dims, p), broadcast(TupleTools.getindices, strides, (p,)) -end - # Fuse dimensions if possible: merge dimension `i` into `i-1` when the two are # contiguous in *every* array (`s[i] == dims[i-1] * s[i-1]`). # Merged dimensions fold their extent into `i-1` and become size 1. @@ -158,10 +138,46 @@ function _fusedims(dims, allstrides) return dims end +# Compute stride orders once: fusing only changes dimensions, so the surviving +# dimensions retain their importance order and singleton dimensions move last. +@inline function _mapreduce_plan(dims::Dims, allstrides::Tuple) + orders = map(indexorder, allstrides) + p1 = _sortperm_desc(_importance(dims, orders)) + strides1 = map(s -> TupleTools.getindices(s, p1), allstrides) + dims1 = _fusedims(TupleTools.getindices(dims, p1), strides1) + p2 = _sortperm_desc(map(d -> Int(d > 1), dims1)) + p = TupleTools.getindices(p1, p2) + return TupleTools.getindices(dims1, p2), + map(s -> TupleTools.getindices(s, p2), strides1), + map(o -> TupleTools.getindices(o, p), orders) +end + function order_and_fuse_dims(dims, allstrides) - dims, allstrides = _sortdims(dims, allstrides) - dims = _fusedims(dims, allstrides) - return _sortdims(dims, allstrides) + dims, allstrides, _ = _mapreduce_plan(dims, allstrides) + return dims, allstrides +end + +# Weight the output array twice when choosing a cache-friendly loop order. +function _importance(dims::NTuple{N, Int}, orders::NTuple{M, NTuple{N, Int}}) where {N, M} + g = 8 * sizeof(Int) - leading_zeros(M + 1) + importance = 2 .* (1 .<< (g .* (N .- orders[1]))) + for k in 2:M + importance = importance .+ (1 .<< (g .* (N .- orders[k]))) + end + return importance .* (dims .> 1) +end + +# Stable descending insertion sort for short tuples. +function _sortperm_desc(v::NTuple{N, Int}) where {N} + p = ntuple(identity, Val(N)) + @inbounds for i in 2:N + j = i + while j > 1 && v[p[j - 1]] < v[p[j]] + p = TupleTools.setindex(TupleTools.setindex(p, p[j], j - 1), p[j - 1], j) + j -= 1 + end + end + return p end # Per-dimension cost used by the blocking and thread-splitting heuristics, @@ -170,14 +186,14 @@ _computecosts(strides) = map(a -> ifelse(iszero(a), 1, a << 1), map(min, strides # Pipeline entry point: order → fuse → order → block → kernel. # Fusing needs ordered entries, and final order pass brings remaining dim 1 to end -function _mapreduce_order!( +@inline function _mapreduce_order!( f::F, op::OP, initop::I, dims::Dims, arrays::Tuple{Vararg{StridedView}} ) where {F, OP, I} isempty(dims) && return _mapreduce_scalar!(f, op, initop, arrays) - dims, allstrides = order_and_fuse_dims(dims, map(strides, arrays)) + dims, allstrides, strideorders = _mapreduce_plan(dims, map(strides, arrays)) offsets = map(offset, arrays) costs = _computecosts(allstrides) - return _mapreduce_block!(f, op, initop, dims, allstrides, offsets, costs, arrays) + return _mapreduce_block!(f, op, initop, dims, allstrides, offsets, costs, strideorders, arrays) end # 0-dimensional fast path: bypass @generated kernel @@ -196,15 +212,11 @@ end const MINTHREADLENGTH = 1 << 15 # minimal length before any kind of threading is applied function _mapreduce_block!( - f::F, op::OP, initop::I, dims, strides, offsets, costs, arrays + f::F, op::OP, initop::I, dims, strides, offsets, costs, strideorders, arrays ) where {F, OP, I} bytestrides = map((s, stride) -> s .* stride, sizeof.(eltype.(arrays)), strides) - strideorders = map(indexorder, strides) blocks = _computeblocks(dims, costs, bytestrides, strideorders) - # t = @elapsed _computeblocks(dims, costs, bytestrides, strideorders) - # println("_computeblocks time: $t") - if get_num_threads() == 1 || prod(dims) <= MINTHREADLENGTH _mapreduce_kernel!(f, op, initop, dims, blocks, arrays, strides, offsets) elseif op !== nothing && _length(dims, strides[1]) == 1 # complete reduction From a565fb373199b796fb95f6f9cb6b601b1b0d6310 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 7 Oct 2026 10:32:28 -0400 Subject: [PATCH 2/7] Keep execution boundaries unchanged when optimizing dimension order --- ext/StridedGPUArraysExt.jl | 2 +- src/mapreduce.jl | 28 +++++++++++++--------------- 2 files changed, 14 insertions(+), 16 deletions(-) diff --git a/ext/StridedGPUArraysExt.jl b/ext/StridedGPUArraysExt.jl index aadfef0..7525b44 100644 --- a/ext/StridedGPUArraysExt.jl +++ b/ext/StridedGPUArraysExt.jl @@ -128,7 +128,7 @@ end function Strided._mapreduce_block!( f, op, initop, dims::Dims{N}, - strides, offsets, costs, strideorders, + strides, offsets, costs, arrays::Tuple{GPUStridedView{TO, N}, Vararg{GPUStridedView{<:Any, N}}} ) where {TO, N} diff --git a/src/mapreduce.jl b/src/mapreduce.jl index 77959aa..33fe3a7 100644 --- a/src/mapreduce.jl +++ b/src/mapreduce.jl @@ -138,23 +138,17 @@ function _fusedims(dims, allstrides) return dims end -# Compute stride orders once: fusing only changes dimensions, so the surviving -# dimensions retain their importance order and singleton dimensions move last. -@inline function _mapreduce_plan(dims::Dims, allstrides::Tuple) +# Compute stride orders once for dimension ordering: fusing only changes sizes, +# so surviving dimensions retain their importance order and singleton dimensions +# move last. +function order_and_fuse_dims(dims::Dims, allstrides::Tuple) orders = map(indexorder, allstrides) p1 = _sortperm_desc(_importance(dims, orders)) strides1 = map(s -> TupleTools.getindices(s, p1), allstrides) dims1 = _fusedims(TupleTools.getindices(dims, p1), strides1) p2 = _sortperm_desc(map(d -> Int(d > 1), dims1)) - p = TupleTools.getindices(p1, p2) return TupleTools.getindices(dims1, p2), - map(s -> TupleTools.getindices(s, p2), strides1), - map(o -> TupleTools.getindices(o, p), orders) -end - -function order_and_fuse_dims(dims, allstrides) - dims, allstrides, _ = _mapreduce_plan(dims, allstrides) - return dims, allstrides + map(s -> TupleTools.getindices(s, p2), strides1) end # Weight the output array twice when choosing a cache-friendly loop order. @@ -186,14 +180,14 @@ _computecosts(strides) = map(a -> ifelse(iszero(a), 1, a << 1), map(min, strides # Pipeline entry point: order → fuse → order → block → kernel. # Fusing needs ordered entries, and final order pass brings remaining dim 1 to end -@inline function _mapreduce_order!( +function _mapreduce_order!( f::F, op::OP, initop::I, dims::Dims, arrays::Tuple{Vararg{StridedView}} ) where {F, OP, I} isempty(dims) && return _mapreduce_scalar!(f, op, initop, arrays) - dims, allstrides, strideorders = _mapreduce_plan(dims, map(strides, arrays)) + dims, allstrides = order_and_fuse_dims(dims, map(strides, arrays)) offsets = map(offset, arrays) costs = _computecosts(allstrides) - return _mapreduce_block!(f, op, initop, dims, allstrides, offsets, costs, strideorders, arrays) + return _mapreduce_block!(f, op, initop, dims, allstrides, offsets, costs, arrays) end # 0-dimensional fast path: bypass @generated kernel @@ -212,11 +206,15 @@ end const MINTHREADLENGTH = 1 << 15 # minimal length before any kind of threading is applied function _mapreduce_block!( - f::F, op::OP, initop::I, dims, strides, offsets, costs, strideorders, arrays + f::F, op::OP, initop::I, dims, strides, offsets, costs, arrays ) where {F, OP, I} bytestrides = map((s, stride) -> s .* stride, sizeof.(eltype.(arrays)), strides) + strideorders = map(indexorder, strides) blocks = _computeblocks(dims, costs, bytestrides, strideorders) + # t = @elapsed _computeblocks(dims, costs, bytestrides, strideorders) + # println("_computeblocks time: $t") + if get_num_threads() == 1 || prod(dims) <= MINTHREADLENGTH _mapreduce_kernel!(f, op, initop, dims, blocks, arrays, strides, offsets) elseif op !== nothing && _length(dims, strides[1]) == 1 # complete reduction From c1d4d8c426b74a32e9bb8cb3e5a33d8815421b25 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 7 Oct 2026 11:17:52 -0400 Subject: [PATCH 3/7] Clarify loop planning and use branchless comparison counting --- src/mapreduce.jl | 76 +++++++++++++++++++++++++++++++----------------- 1 file changed, 50 insertions(+), 26 deletions(-) diff --git a/src/mapreduce.jl b/src/mapreduce.jl index 33fe3a7..547e39d 100644 --- a/src/mapreduce.jl +++ b/src/mapreduce.jl @@ -138,40 +138,65 @@ function _fusedims(dims, allstrides) return dims end -# Compute stride orders once for dimension ordering: fusing only changes sizes, -# so surviving dimensions retain their importance order and singleton dimensions -# move last. +# Fusing only changes sizes, so surviving dimensions retain their importance +# order. The second permutation only moves singleton dimensions to the back. function order_and_fuse_dims(dims::Dims, allstrides::Tuple) - orders = map(indexorder, allstrides) - p1 = _sortperm_desc(_importance(dims, orders)) - strides1 = map(s -> TupleTools.getindices(s, p1), allstrides) - dims1 = _fusedims(TupleTools.getindices(dims, p1), strides1) - p2 = _sortperm_desc(map(d -> Int(d > 1), dims1)) - return TupleTools.getindices(dims1, p2), - map(s -> TupleTools.getindices(s, p2), strides1) + stride_orders = map(indexorder, allstrides) + loop_order = _sortperm_desc(_importance(dims, stride_orders)) + ordered_dims = TupleTools.getindices(dims, loop_order) + ordered_strides = map(s -> TupleTools.getindices(s, loop_order), allstrides) + + fused_dims = _fusedims(ordered_dims, ordered_strides) + active_order = _sortperm_desc(map(d -> Int(d > 1), fused_dims)) + return TupleTools.getindices(fused_dims, active_order), + map(s -> TupleTools.getindices(s, active_order), ordered_strides) end -# Weight the output array twice when choosing a cache-friendly loop order. -function _importance(dims::NTuple{N, Int}, orders::NTuple{M, NTuple{N, Int}}) where {N, M} - g = 8 * sizeof(Int) - leading_zeros(M + 1) - importance = 2 .* (1 .<< (g .* (N .- orders[1]))) +# Each stride rank gets enough bits to hold all array votes without carries. +# The output array gets two votes; each input gets one. +function _importance(dims::NTuple{N, Int}, stride_orders::NTuple{M, NTuple{N, Int}}) where {N, M} + bits_per_rank = 8 * sizeof(Int) - leading_zeros(M + 1) + importance = 2 .* (1 .<< (bits_per_rank .* (N .- stride_orders[1]))) for k in 2:M - importance = importance .+ (1 .<< (g .* (N .- orders[k]))) + importance = importance .+ (1 .<< (bits_per_rank .* (N .- stride_orders[k]))) end return importance .* (dims .> 1) end -# Stable descending insertion sort for short tuples. -function _sortperm_desc(v::NTuple{N, Int}) where {N} - p = ntuple(identity, Val(N)) - @inbounds for i in 2:N - j = i - while j > 1 && v[p[j - 1]] < v[p[j]] - p = TupleTools.setindex(TupleTools.setindex(p, p[j], j - 1), p[j - 1], j) - j -= 1 +# Tiny tuples need at most one swap; longer tuples use comparison counting. +_sortperm_desc(values::NTuple{N, Int}) where {N} = + N <= 2 ? _sortperm_desc_insertion(values) : _sortperm_desc_ranked(values) + +function _sortperm_desc_insertion(values::NTuple{N, Int}) where {N} + permutation = ntuple(identity, Val(N)) + @inbounds for source in 2:N + position = source + while position > 1 && values[permutation[position - 1]] < values[permutation[position]] + permutation = TupleTools.setindex( + TupleTools.setindex(permutation, permutation[position], position - 1), + permutation[position - 1], position + ) + position -= 1 + end + end + return permutation +end + +# Stable descending permutation: count larger values and earlier equal values. +function _sortperm_desc_ranked(values::NTuple{N, Int}) where {N} + permutation = ntuple(identity, Val(N)) + issorted(values; rev = true) && return permutation + + @inbounds for source in 1:N + value = values[source] + position = 1 + @simd for other in 1:N + other_value = values[other] + position += (other_value > value) | ((other_value == value) & (other < source)) end + permutation = TupleTools.setindex(permutation, source, position) end - return p + return permutation end # Per-dimension cost used by the blocking and thread-splitting heuristics, @@ -556,8 +581,7 @@ end end function indexorder(strides::NTuple{N, Int}) where {N} - # returns order such that strides[i] is the order[i]th smallest element of strides, not - # counting zero strides zero strides have order 1 + # Rank absolute nonzero strides. Ties share a rank; zero strides get rank 1. return ntuple(Val(N)) do i si = abs(strides[i]) iszero(si) && return 1 From 46d8bf23bd99a7ffd2793228c9c907dd2e6f6b35 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 7 Oct 2026 13:16:24 -0400 Subject: [PATCH 4/7] Dispatch small stride permutations directly --- src/mapreduce.jl | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/src/mapreduce.jl b/src/mapreduce.jl index 547e39d..d265cbd 100644 --- a/src/mapreduce.jl +++ b/src/mapreduce.jl @@ -163,27 +163,22 @@ function _importance(dims::NTuple{N, Int}, stride_orders::NTuple{M, NTuple{N, In return importance .* (dims .> 1) end -# Tiny tuples need at most one swap; longer tuples use comparison counting. -_sortperm_desc(values::NTuple{N, Int}) where {N} = - N <= 2 ? _sortperm_desc_insertion(values) : _sortperm_desc_ranked(values) - -function _sortperm_desc_insertion(values::NTuple{N, Int}) where {N} - permutation = ntuple(identity, Val(N)) - @inbounds for source in 2:N - position = source - while position > 1 && values[permutation[position - 1]] < values[permutation[position]] - permutation = TupleTools.setindex( - TupleTools.setindex(permutation, permutation[position], position - 1), - permutation[position - 1], position - ) - position -= 1 - end +# Small tuples have direct stable permutations; ties keep their original order. +_sortperm_desc(::Tuple{}) = () +_sortperm_desc(::Tuple{Int}) = (1,) +_sortperm_desc(values::NTuple{2, Int}) = ifelse(values[1] < values[2], (2, 1), (1, 2)) + +function _sortperm_desc(values::NTuple{3, Int}) + a, b, c = values + if a >= b + return b >= c ? (1, 2, 3) : a >= c ? (1, 3, 2) : (3, 1, 2) + else + return a >= c ? (2, 1, 3) : b >= c ? (2, 3, 1) : (3, 2, 1) end - return permutation end # Stable descending permutation: count larger values and earlier equal values. -function _sortperm_desc_ranked(values::NTuple{N, Int}) where {N} +function _sortperm_desc(values::NTuple{N, Int}) where {N} permutation = ntuple(identity, Val(N)) issorted(values; rev = true) && return permutation From 1b937eba511db5bca66f7a53913521330c9871e7 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Wed, 7 Oct 2026 13:19:55 -0400 Subject: [PATCH 5/7] Spell out descending permutation helper name --- src/mapreduce.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/mapreduce.jl b/src/mapreduce.jl index d265cbd..5c7637c 100644 --- a/src/mapreduce.jl +++ b/src/mapreduce.jl @@ -142,12 +142,12 @@ end # order. The second permutation only moves singleton dimensions to the back. function order_and_fuse_dims(dims::Dims, allstrides::Tuple) stride_orders = map(indexorder, allstrides) - loop_order = _sortperm_desc(_importance(dims, stride_orders)) + loop_order = _sortperm_descending(_importance(dims, stride_orders)) ordered_dims = TupleTools.getindices(dims, loop_order) ordered_strides = map(s -> TupleTools.getindices(s, loop_order), allstrides) fused_dims = _fusedims(ordered_dims, ordered_strides) - active_order = _sortperm_desc(map(d -> Int(d > 1), fused_dims)) + active_order = _sortperm_descending(map(d -> Int(d > 1), fused_dims)) return TupleTools.getindices(fused_dims, active_order), map(s -> TupleTools.getindices(s, active_order), ordered_strides) end @@ -163,12 +163,13 @@ function _importance(dims::NTuple{N, Int}, stride_orders::NTuple{M, NTuple{N, In return importance .* (dims .> 1) end -# Small tuples have direct stable permutations; ties keep their original order. -_sortperm_desc(::Tuple{}) = () -_sortperm_desc(::Tuple{Int}) = (1,) -_sortperm_desc(values::NTuple{2, Int}) = ifelse(values[1] < values[2], (2, 1), (1, 2)) +# Return indices that order values from largest to smallest, preserving ties. +# For example, (2, 5, 2) gives (2, 1, 3). Small tuples use direct permutations. +_sortperm_descending(::Tuple{}) = () +_sortperm_descending(::Tuple{Int}) = (1,) +_sortperm_descending(values::NTuple{2, Int}) = ifelse(values[1] < values[2], (2, 1), (1, 2)) -function _sortperm_desc(values::NTuple{3, Int}) +function _sortperm_descending(values::NTuple{3, Int}) a, b, c = values if a >= b return b >= c ? (1, 2, 3) : a >= c ? (1, 3, 2) : (3, 1, 2) @@ -178,7 +179,7 @@ function _sortperm_desc(values::NTuple{3, Int}) end # Stable descending permutation: count larger values and earlier equal values. -function _sortperm_desc(values::NTuple{N, Int}) where {N} +function _sortperm_descending(values::NTuple{N, Int}) where {N} permutation = ntuple(identity, Val(N)) issorted(values; rev = true) && return permutation From e3a99b6d2e0500dff1c6178688b2beb04c0f8a75 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 07:59:02 -0400 Subject: [PATCH 6/7] Use TupleTools for descending dimension permutations --- src/mapreduce.jl | 36 ++---------------------------------- 1 file changed, 2 insertions(+), 34 deletions(-) diff --git a/src/mapreduce.jl b/src/mapreduce.jl index 5c7637c..d361419 100644 --- a/src/mapreduce.jl +++ b/src/mapreduce.jl @@ -142,12 +142,12 @@ end # order. The second permutation only moves singleton dimensions to the back. function order_and_fuse_dims(dims::Dims, allstrides::Tuple) stride_orders = map(indexorder, allstrides) - loop_order = _sortperm_descending(_importance(dims, stride_orders)) + loop_order = TupleTools.sortperm(_importance(dims, stride_orders); rev = true) ordered_dims = TupleTools.getindices(dims, loop_order) ordered_strides = map(s -> TupleTools.getindices(s, loop_order), allstrides) fused_dims = _fusedims(ordered_dims, ordered_strides) - active_order = _sortperm_descending(map(d -> Int(d > 1), fused_dims)) + active_order = TupleTools.sortperm(map(d -> Int(d > 1), fused_dims); rev = true) return TupleTools.getindices(fused_dims, active_order), map(s -> TupleTools.getindices(s, active_order), ordered_strides) end @@ -163,38 +163,6 @@ function _importance(dims::NTuple{N, Int}, stride_orders::NTuple{M, NTuple{N, In return importance .* (dims .> 1) end -# Return indices that order values from largest to smallest, preserving ties. -# For example, (2, 5, 2) gives (2, 1, 3). Small tuples use direct permutations. -_sortperm_descending(::Tuple{}) = () -_sortperm_descending(::Tuple{Int}) = (1,) -_sortperm_descending(values::NTuple{2, Int}) = ifelse(values[1] < values[2], (2, 1), (1, 2)) - -function _sortperm_descending(values::NTuple{3, Int}) - a, b, c = values - if a >= b - return b >= c ? (1, 2, 3) : a >= c ? (1, 3, 2) : (3, 1, 2) - else - return a >= c ? (2, 1, 3) : b >= c ? (2, 3, 1) : (3, 2, 1) - end -end - -# Stable descending permutation: count larger values and earlier equal values. -function _sortperm_descending(values::NTuple{N, Int}) where {N} - permutation = ntuple(identity, Val(N)) - issorted(values; rev = true) && return permutation - - @inbounds for source in 1:N - value = values[source] - position = 1 - @simd for other in 1:N - other_value = values[other] - position += (other_value > value) | ((other_value == value) & (other < source)) - end - permutation = TupleTools.setindex(permutation, source, position) - end - return permutation -end - # Per-dimension cost used by the blocking and thread-splitting heuristics, # derived from the minimum (doubled) stride across arrays. _computecosts(strides) = map(a -> ifelse(iszero(a), 1, a << 1), map(min, strides...)) From 2e553dfc454108c69e6145e0245c790a15bd2790 Mon Sep 17 00:00:00 2001 From: lkdvos Date: Thu, 8 Oct 2026 11:00:15 -0400 Subject: [PATCH 7/7] add explanatory comment --- src/mapreduce.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mapreduce.jl b/src/mapreduce.jl index d361419..77adf57 100644 --- a/src/mapreduce.jl +++ b/src/mapreduce.jl @@ -155,7 +155,7 @@ end # Each stride rank gets enough bits to hold all array votes without carries. # The output array gets two votes; each input gets one. function _importance(dims::NTuple{N, Int}, stride_orders::NTuple{M, NTuple{N, Int}}) where {N, M} - bits_per_rank = 8 * sizeof(Int) - leading_zeros(M + 1) + bits_per_rank = 8 * sizeof(Int) - leading_zeros(M + 1) # ceil(Int, log2(M + 2)) importance = 2 .* (1 .<< (bits_per_rank .* (N .- stride_orders[1]))) for k in 2:M importance = importance .+ (1 .<< (bits_per_rank .* (N .- stride_orders[k])))