# Interleaving bits (Z-order curve,morton code).Transpose binary matrix

**URL:** <https://discourse.julialang.org/t/interleaving-bits-z-order-curve-morton-code-transpose-binary-matrix/18458>\
**Category:** Performance\
**Created:** [December 8, 2018, 10:14am UTC](https://discourse.julialang.org/t/interleaving-bits-z-order-curve-morton-code-transpose-binary-matrix/18458 "2018-12-08T10:14:53Z")\
**Posts on this page:** 1\
**Showing post:** 6

<div class="post-metadata">

**Author:** ![bennedich](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bennedich/32/4894_2.png) [@bennedich](https://discourse.julialang.org/u/bennedich)\
**Post date:** [December 13, 2018, 5:10pm UTC](https://discourse.julialang.org/t/interleaving-bits-z-order-curve-morton-code-transpose-binary-matrix/18458/6 "2018-12-13T17:10:17Z")

</div>

Great solution (and website)!

Some other ideas. First, if you could pack your bits as 4x16 (with 3 bits padding) instead of 4x13, you could probably speed things up a bit.

Second, consider a look-up-table. You can use any number of bits for the table. I’ve chosen 13 below, quite arbitrarily, since it divides 52, but other sizes might perform better depending on how it fits into cache:

```julia
const PERM13 = zeros(Int64, 4*2^13)

for p = 0:3, n = 0:2^13-1
    PERM13[n+1 + p*2^13] = uninterleavenibble52(n << (p*13))
end

function permute_lookup(x)
    x & 0xfff0000000000000 | # remove this if you don't need to preserve the top 12 bits
        PERM13[1 + (x & 0b1111111111111)] |
        PERM13[8193 + ((x >> 13) & 0b1111111111111)] |
        PERM13[16385 + ((x >> 26) & 0b1111111111111)] |
        PERM13[24577 + ((x >> 39) & 0b1111111111111)]
end

```

In benchmarks this performs great, presumably because it fits nicely in cache, but in practice that might not be the case, so you’d need to test it out.

Third, consider using the `pext` instruction available on Haswell and newer Intel processors (after 2013) to greatly simplify/speed up the permutation:

```julia
@inline pext(x::UInt64, y::UInt64) = ccall("llvm.x86.bmi.pext.64", llvmcall, UInt64, (UInt64, UInt64), x, y)

function permute_pext(x)
    x & 0xfff0000000000000 | # remove this if you don't need to preserve the top 12 bits
        pext(x, 0x1111111111111) |
        pext(x, 0x2222222222222) << 13 |
        pext(x, 0x4444444444444) << 26 |
        pext(x, 0x8888888888888) << 39
end

```

To figure out if the `pext` instruction is available, you could do something like this:

```julia
CPUInfo = zeros(Int32, 4)
ccall(:jl_cpuidex, Cvoid, (Ptr{Cint}, Cint, Cint), CPUInfo, 7, 0)
has_pext = CPUInfo[2] & 0x100 != 0

```

Sample timings below. Note that I’ve put `@noinline` in front of all permute functions to prevent inlining.

```julia
julia> using BenchmarkTools

julia> A = rand(UInt64, 100000);

julia> @btime (for n = $A uninterleavenibble52(n) end)
  702.064 μs (0 allocations: 0 bytes)

julia> @btime (for n = $A permute_lookup(n) end)
  371.155 μs (0 allocations: 0 bytes)

julia> @btime (for n = $A permute_pext(n) end)
  200.630 μs (0 allocations: 0 bytes)

```

This corresponds to ~20.4 cycles for the bit permute implementation, ~10.8 cycles for the look-up-table  
and ~5.8 cycles for the pext implementation.

---

_[View the full topic](https://discourse.julialang.org/t/interleaving-bits-z-order-curve-morton-code-transpose-binary-matrix/18458)._