3D-кодирование Мортона с использованием чередования битов, обычный набор инструкций по сравнению с BMI2

Question

Я хочу написать две функции для Morton Z-Order Encoding и Decoding на C в быстром и эффективном виде, а именно.

uint64_t morton_encode(uint32_t xindex, uint32_t yindex, uint32_t zindex);
void morton_decode(uint64_t morton_number, uint32_t *xindex, uint32_t *yindex, uint32_t *zindex);

Я ранее следовал за вопросами

как вычислить число 3-х миньонов с чередованием битов по 3 дюйма

Мое текущее решение на основе SO и открытых исходных кодов

uint64_t spread(uint64_t w)  {
    w &=                0x00000000001fffff; 
    w = (w | w << 32) & 0x001f00000000ffff;  
    w = (w | w << 16) & 0x001f0000ff0000ff;  
    w = (w | w <<  8) & 0x010f00f00f00f00f; 
    w = (w | w <<  4) & 0x10c30c30c30c30c3; 
    w = (w | w <<  2) & 0x1249249249249249;
    return w;
    }

uint64_t morton_encode(uint32_t x, uint32_t y, uint32_t z)  {
   return ((spread((uint64_t)x)) | (spread((uint64_t)y) << 1) | (spread((uint64_t)z) << 2));
   }

///////////////// For Decoding //////////////////////

uint32_t compact(uint64_t w) {
    w &=                  0x1249249249249249;
    w = (w ^ (w >> 2))  & 0x30c30c30c30c30c3;
    w = (w ^ (w >> 4))  & 0xf00f00f00f00f00f;
    w = (w ^ (w >> 8))  & 0x00ff0000ff0000ff;
    w = (w ^ (w >> 16)) & 0x00ff00000000ffff;
    w = (w ^ (w >> 32)) & 0x00000000001fffff;
    return (uint32_t)w;
    }

void morton_decode(uint64_t morton_number, uint32_t *xindex, uint32_t *yindex, uint32_t *zindex){
    *xindex = compact(code);
    *yindex = compact(code >> 1);
    *zindex = compact(code >> 2);
}

Недавно я столкнулся с таким SO вопросом (пытаясь поиграться с 2D-кодом Morton): 2-й Mort-код кодирует декодировать 64 бит

#include <immintrin.h>
#include <stdint.h>

// on GCC, compile with option -mbmi2, requires Haswell or better.

uint64_t xy_to_morton (uint32_t x, uint32_t y)
{
  return _pdep_u32(x, 0x55555555) | _pdep_u32(y,0xaaaaaaaa);
}

uint64_t morton_to_xy (uint64_t m, uint32_t *x, uint32_t *y)
{
  *x = _pext_u64(m, 0x5555555555555555);
  *y = _pext_u64(m, 0xaaaaaaaaaaaaaaaa);
}

Из того, что я понимаю, это НЕ портативное решение, но так как каждая система, на которой я (буду) запускать свой код, имеет процессор Haswell (даже на кластере HPC). Мои вопросы:

1. Как изменить этот код для 3D-системы или эти наборы инструкций BMI можно использовать для кодирования декодирования 3D-числа Мортона?
2. Является ли / будет ли более эффективно использовать эти инструкции над стандартным решением, которое я сейчас использую, учитывая случай, когда мне нужно декодировать несколько миллионов чисел за раз на каждом временном шаге, а таких временных шагов миллион?

Редактировать: В первом квартале я уже близко к решению, но все еще не могу разобраться

0x55555555 -> 0000 0000 0101 0101 0101 0101 0101 0101 0101 0101 
0xaaaaaaaa -> 0000 0000 1010 1010 1010 1010 1010 1010 1010 1010

очевидно, что маски представляют собой чередующиеся биты x и y. Так что для 3d мне нужно получить маску типа

0000 0000 01 001 001 001 001 001 001 001 001 001 001 (for x)
0000 0000 01 010 010 010 010 010 010 010 010 010 010 (for y)
0000 0000 01 100 100 100 100 100 100 100 100 100 100 (for z)
           ^

Я немного сбит с толку насчет битов до меток ^ для 64-битного кода Morton, только первые 21 бит x, y и z, которые являются 32-битными целыми числами, должны иметь значение.

Answer 1

Итак, немного повозившись, я пришел к решению, которое, как мне показалось, должно быть здесь в качестве ответа.

// on GCC, compile with option -mbmi2, requires Haswell or better.
#include <stdio.h>
#include <limits.h>
#include <immintrin.h>
#include <inttypes.h>
#include <sys/time.h>

#define maask 0x1249249249249249

/* Morton Encoding Mehtod 1 */
uint64_t Z_encode1 (uint32_t x, uint32_t y, uint32_t z)
{
  return _pdep_u32(x, maask)       | \
         _pdep_u32(y,(maask << 1)) | \
         _pdep_u32(z,(maask << 2));
}

/* Morton Decoding Method 1 */
uint64_t Z_decode1 (uint64_t m, uint32_t *x, uint32_t *y, uint32_t *z)
{
  *x = _pext_u64(m, maask);
  *y = _pext_u64(m, (maask << 1));
  *z = _pext_u64(m, (maask << 2));
}

// method 2 functions 
uint64_t spread(uint64_t w)  {
    w &=                0x00000000001fffff; 
    w = (w | w << 32) & 0x001f00000000ffff;  
    w = (w | w << 16) & 0x001f0000ff0000ff;  
    w = (w | w <<  8) & 0x010f00f00f00f00f; 
    w = (w | w <<  4) & 0x10c30c30c30c30c3; 
    w = (w | w <<  2) & 0x1249249249249249;
    return w;
    }

uint32_t compact(uint64_t w) {
    w &=                  0x1249249249249249;
    w = (w ^ (w >> 2))  & 0x30c30c30c30c30c3;
    w = (w ^ (w >> 4))  & 0xf00f00f00f00f00f;
    w = (w ^ (w >> 8))  & 0x00ff0000ff0000ff;
    w = (w ^ (w >> 16)) & 0x00ff00000000ffff;
    w = (w ^ (w >> 32)) & 0x00000000001fffff;
    return (uint32_t)w;
    }

uint64_t Z_encode2(uint32_t x, uint32_t y, uint32_t z)  {
   return ((spread((uint64_t)x)) | (spread((uint64_t)y) << 1) | (spread((uint64_t)z) << 2));
   }



void Z_decode2(uint64_t Z_code, uint32_t *xindex, uint32_t *yindex, uint32_t *zindex){
    *xindex = compact(Z_code);
    *yindex = compact(Z_code >> 1);
    *zindex = compact(Z_code >> 2);
}
int main()
{
    const int size = 1024;
    struct timeval start, stop;
    double time_encode1 = 0.0, time_encode2 = 0.0;
    double time_decode1 = 0.0, time_decode2 = 0.0;

    uint64_t Zindex = 0;
    uint32_t xindex=0,yindex=0,zindex=0;

    /* method 1 ENCODING benchmark */
    gettimeofday(&start, NULL);
    for (uint32_t i = 0; i < size; i++){
        for (uint32_t j = 0; j < size; j++) {
            for (uint32_t k = 0; k < size; k++) {
                Zindex = Z_encode1(i, j, k);
            }
        }
    }
    gettimeofday(&stop, NULL);
    time_encode1 = (double)(stop.tv_usec - start.tv_usec) / 1000000 + (double)(stop.tv_sec - start.tv_sec);

    /* method 2 ENCODING benchmark */
    gettimeofday(&start, NULL);
    for (uint32_t i = 0; i < size; i++){
        for (uint32_t j = 0; j < size; j++) {
            for (uint32_t k = 0; k < size; k++) {
                Zindex = Z_encode2(i, j, k);
            }
        }
    }
    gettimeofday(&stop, NULL);
    time_encode2 = (double)(stop.tv_usec - start.tv_usec) / 1000000 + (double)(stop.tv_sec - start.tv_sec);

    //////////////////////// DECODING ////////////////
    /* method 1 DECODING benchmark */
    gettimeofday(&start, NULL);
    for (uint64_t i = 0; i < size; i++)
        Z_decode1(i, &xindex, &yindex, &zindex);
    gettimeofday(&stop, NULL);
    time_decode1 = (double)(stop.tv_usec - start.tv_usec) / 1000000 + (double)(stop.tv_sec - start.tv_sec);

    /* method 1 DECODING benchmark */
    gettimeofday(&start, NULL);
    for (uint64_t i = 0; i < size; i++)
        Z_decode2(i, &xindex, &yindex, &zindex);
    gettimeofday(&stop, NULL);
    time_decode2 = (double)(stop.tv_usec - start.tv_usec) / 1000000 + (double)(stop.tv_sec - start.tv_sec);

    printf("Method1 -> Encoding: %f Decoding: %f\n", time_encode1, time_decode1);
    printf("Method2 -> Encoding: %f Decoding: %f\n", time_encode2, time_decode2);
    return 0;
}

Вот результаты

size = 512 ( 512x512x512 = 134217728 numbers)
======================================================
Method 1 -> Encoding: 0.600302sec Decoding: 0.000003sec
Method 2 -> Encoding: 2.778170sec Decoding: 0.000011sec

size = 1024 ( 1024x1024x1024 = 1073741824 numbers)
======================================================
Method 1 -> Encoding:  4.623594sec Decoding: 0.000006sec
Method 2 -> Encoding: 22.339238sec Decoding: 0.000022sec

size = 2048 ( 2048*2048*2048 = 8589934592 numbers)
======================================================
Method 1 -> Encoding:  36.981743sec Decoding: 0.000011sec
Method 2 -> Encoding: 178.164773sec Decoding: 0.000045sec

Вывод: кодирование обходится дороже, чем декодирование, используйте набор инструкций BMI для оптимизации производительности.

PS. - не переносим, ​​так как нуждается в процессоре Haswell или выше.