How to sort points along a Hilbert curve without using Hilbert indices?

Question

I'm trying to implement the algorithm described in the paper Fast Hilbert Sort Algorithm Without Using Hilbert Indices (https://www.researchgate.net/profile/Takeshi_Shinohara/publication/313074453_Fast_Hilbert_Sort_Algorithm_Without_Using_Hilbert_Indices/links/5b8468bd299bf1d5a72b9a0c/Fast-Hilbert-Sort-Algorithm-Without-Using-Hilbert-Indices.pdf?origin=publication_detail), but I can't get the right results.

Below is my python code (For bitset and it's member functions flip and test in C++ , please refer to https://en.cppreference.com/w/cpp/utility/bitset):

N=9 # 9 points
n=2 # 2 dimension 
m=3 # order of Hilbert curve
b=m-1
def BitTest(x,od,maxlen=3):
    bit=format(x,'b').zfill(maxlen)
    return int(bit[maxlen-1-od])

def BitFlip(b,pos,):
    b ^= 1 << pos
    return b
def partition(A,st,en,od,ax,di):
    i = st
    j = en
    while True:
        while i < j and BitTest(A[i][ax],od)==di:
            i = i + 1
        while i < j and BitTest(A[j][ax],od)!=di:
            j = j - 1
        if i >= j:
            return i
        A[i], A[j] = A[j], A[i]

def HSort(A,st,en,od,c,e,d,di,cnt):
    if en<=st: 
        return
    p =partition(A,st,en,od,(d+c)%n,BitTest(e,(d+c)%n))
    if c==n-1:
        if b==0:
            return
        d2= (d+n+n-(di if(di==2) else cnt+2))%n
        e=BitFlip(e,d2)
        e=BitFlip(e,(d+c)%n)
        HSort(A,st,p-1,b-1,0,e,d2,False,0)
        
        e=BitFlip(e,(d+c)%n)
        e=BitFlip(e,d2)
        d2= (d+n+n-(di if(di==cnt+2) else 2))%n
        HSort(A,p+1,en,b-1,0,e,d2,False,0)
    else:
        HSort(A,st,p-1,b,c+1,e,d,False,(di if(di==1) else cnt+1))
        e=BitFlip(e,(d+c)%n)
        e=BitFlip(e,(d+c+1)%n)
        HSort(A,p+1,en,b,c+1,e,d,True,(di if(di==cnt+1) else 1))
        e=BitFlip(e,(d+c+1)%n)
        e=BitFlip(e,(d+c)%n)
        
array = [[2,2],[2,4],[3,4],[2,5],[3,5],[1,6],[3,6],[5,6],[3,7]]
HSort(array,st=0,en=N-1,od=m-1,c=0,e=0,d=0,di=False,cnt=0)
print(array)

Answer 1

That document has a typo, the constant "b" should be replaced with "od". Here is a working code in c++:

#include <iostream>
#include <vector>
#include <array>

constexpr std::int32_t m = 3;
constexpr std::int32_t n = 2;

bool test_bit(std::int32_t value, std::int32_t pos)
{
    const auto result = value & (1 << pos);

    return result;
}

void flip_bit(std::int32_t &value, std::int32_t pos)
{
    value ^= 1 << pos;
}

std::int32_t partition(std::vector<std::array<std::int32_t, 2>> &A, std::size_t st, std::size_t en, std::int32_t od, std::int32_t ax, bool di)
{
    std::int32_t i = st - 1;
    std::int32_t j = en + 1;

    while(true)
    {
        do
            i = i + 1;
        while(i < j && test_bit(A[i][ax], od) == di);

        do
            j = j - 1;
        while(i < j && test_bit(A[j][ax], od) != di);

        if(j <= i)
            return i; //partition is complete

        std::swap(A[i], A[j]);
    }
}

void hilbert_sort(std::vector<std::array<std::int32_t, 2>> &A, std::size_t st, std::size_t en, std::int32_t od, std::int32_t c, std::int32_t &e, std::int32_t 
d, bool di, std::int32_t cnt)
{
    std::int32_t p;
    std::int32_t d2;

    if(en <= st)
        return;

    p = partition(A, st, en, od, (d + c) % n, test_bit(e, (d + c) % n));

    if(c == n - 1)
    {
        if(od == 0)
            return;

        d2 = (d + n + n - (di ? 2 : cnt + 2)) % n;

        flip_bit(e, d2);
        flip_bit(e, (d + c) % n);

        hilbert_sort(A, st, p - 1, od - 1, 0, e, d2, false, 0);

        flip_bit(e, (d + c) % n);
        flip_bit(e, d2);

        d2 = (d + n + n - (di ? cnt + 2 : 2)) % n;

        hilbert_sort(A, p, en, od - 1, 0, e, d2, false, 0);
    }
    else
    {
        hilbert_sort(A, st, p - 1, od, c + 1, e, d, false, di ? 1 : cnt + 1); 

        flip_bit(e, (d + c) % n);
        flip_bit(e, (d + c + 1) % n);

        hilbert_sort(A, p, en, od, c + 1, e, d, true, di ? cnt + 1 : 1); 

        flip_bit(e, (d + c + 1) % n);
        flip_bit(e, (d + c) % n);
    }
}

int main()
{
    std::vector<std::array<std::int32_t, 2>> points = {{2,2},{2,4},{3,4},{2,5},{3,5},{1,6},{3,6},{5,6},{3,7}};

    std::int32_t e = 0;

    hilbert_sort(points, 0, points.size() - 1, m - 1, 0, e, 0, false , 0);

    for(const auto &point : points)
        std::clog << "(" << point[0] << ", " << point[1] << ")\n";

    return 0;
}

You also seems to have a typo "p+1" it should be just "p". Here is a working python code:

N=9 # 9 points
n=2 # 2 dimension 
m=3 # order of Hilbert curve
def BitTest(x,od):
    result = x & (1 << od)
    return int(bool(result))

def BitFlip(b,pos):
    b ^= 1 << pos
    return b
def partition(A,st,en,od,ax,di):
    i = st
    j = en
    while True:
        while i < j and BitTest(A[i][ax],od) == di:
            i = i + 1
        while i < j and BitTest(A[j][ax],od) != di:
            j = j - 1
            
        if j <= i:
            return i
        
        A[i], A[j] = A[j], A[i]

def HSort(A,st,en,od,c,e,d,di,cnt):
    if en<=st: 
        return
    p = partition(A,st,en,od,(d+c)%n,BitTest(e,(d+c)%n))

    if c==n-1:
        if od==0:
            return
        
        d2= (d+n+n-(2 if di else cnt + 2)) % n
        e=BitFlip(e,d2)
        e=BitFlip(e,(d+c)%n)
        HSort(A,st,p-1,od-1,0,e,d2,False,0)
        
        e=BitFlip(e,(d+c)%n)
        e=BitFlip(e,d2)
        d2= (d+n+n-(cnt + 2 if di else 2))%n
        HSort(A,p,en,od-1,0,e,d2,False,0)
    else:
        HSort(A,st,p-1,od,c+1,e,d,False,(1 if di else cnt+1))
        e=BitFlip(e,(d+c)%n)
        e=BitFlip(e,(d+c+1)%n)
        HSort(A,p,en,od,c+1,e,d,True,(cnt+1 if di else 1))
        e=BitFlip(e,(d+c+1)%n)
        e=BitFlip(e,(d+c)%n)
        
array = [[2,2],[2,4],[3,4],[2,5],[3,5],[1,6],[3,6],[5,6],[3,7]]
HSort(array,st=0,en=N-1,od=m-1,c=0,e=0,d=0,di=False,cnt=0)
print(array)