.. _program_listing_file_include_hilbert.h:

Program Listing for File hilbert.h
==================================

|exhale_lsh| :ref:`Return to documentation for file <file_include_hilbert.h>` (``include/hilbert.h``)

.. |exhale_lsh| unicode:: U+021B0 .. UPWARDS ARROW WITH TIP LEFTWARDS

.. code-block:: cpp

   #ifndef hilbert_h
   #define hilbert_h
   #include <cstdlib>
   #include <stdint.h>
   #include "exafmm_t.h"
   #define EXAFMM_HILBERT 0 
   
   namespace exafmm_t {
     inline uint64_t levelOffset(int level) {
       return (((uint64_t)1 << 3 * level) - 1) / 7;
     }
   
     int getLevel(uint64_t i) {
       int level = -1;
       uint64_t offset = 0;
       while (i >= offset) {
         level++;
         offset += (uint64_t)1 << 3 * level;
       }
       return level;
     }
   
     uint64_t getParent(uint64_t i) {
       int level = getLevel(i);
       return (i - levelOffset(level)) / 8 + levelOffset(level-1);
     }
   
     uint64_t getChild(uint64_t i) {
       int level = getLevel(i);
       return (i - levelOffset(level)) * 8 + levelOffset(level+1);
     }
     
     int getOctant(uint64_t key, bool offset=true) {
       int level = getLevel(key);
       if (offset) key -= levelOffset(level);
       return key & 7;
     }
   
     uint64_t getKey(ivec3 iX, int level, bool offset=true) {
   #if EXAFMM_HILBERT
       int M = 1 << (level - 1);
       for (int Q=M; Q>1; Q>>=1) {
         int R = Q - 1;
         for (int d=0; d<3; d++) {
           if (iX[d] & Q) iX[0] ^= R;
           else {
             int t = (iX[0] ^ iX[d]) & R;
             iX[0] ^= t;
             iX[d] ^= t;
           }
         }
       }
       for (int d=1; d<3; d++) iX[d] ^= iX[d-1];
       int t = 0;
       for (int Q=M; Q>1; Q>>=1)
         if (iX[2] & Q) t ^= Q - 1;
       for (int d=0; d<3; d++) iX[d] ^= t;
   #endif
       uint64_t i = 0;
       for (int l=0; l<level; l++) {
         i |= (iX[2] & (uint64_t)1 << l) << 2*l;
         i |= (iX[1] & (uint64_t)1 << l) << (2*l + 1);
         i |= (iX[0] & (uint64_t)1 << l) << (2*l + 2);
       }
       if (offset) i += levelOffset(level);
       return i;
     }
   
     ivec3 get3DIndex(uint64_t i) {
       int level = getLevel(i);
       i -= levelOffset(level);
       ivec3 iX = 0;
       for (int l=0; l<level; l++) {
         iX[2] |= (i & (uint64_t)1 << 3*l) >> 2*l;
         iX[1] |= (i & (uint64_t)1 << (3*l + 1)) >> (2*l + 1);
         iX[0] |= (i & (uint64_t)1 << (3*l + 2)) >> (2*l + 2);
       }
   #if EXAFMM_HILBERT
       int N = 2 << (level - 1);
       int t = iX[2] >> 1;
       for (int d=2; d>0; d--) iX[d] ^= iX[d-1];
       iX[0] ^= t;
       for (int Q=2; Q!=N; Q<<=1) {
         int R = Q - 1;
         for (int d=2; d>=0; d--) {
           if (iX[d] & Q) iX[0] ^= R;
           else {
             t = (iX[0] ^ iX[d]) & R;
             iX[0] ^= t;
             iX[d] ^= t;
           }
         }
       }
   #endif
       return iX;
     }
   
     ivec3 get3DIndex(uint64_t i, int level) {
       ivec3 iX = 0;
       for (int l=0; l<level; l++) {
         iX[2] |= (i & (uint64_t)1 << 3*l) >> 2*l;
         iX[1] |= (i & (uint64_t)1 << (3*l + 1)) >> (2*l + 1);
         iX[0] |= (i & (uint64_t)1 << (3*l + 2)) >> (2*l + 2);
       }
   #if EXAFMM_HILBERT
       int N = 2 << (level - 1);
       int t = iX[2] >> 1;
       for (int d=2; d>0; d--) iX[d] ^= iX[d-1];
       iX[0] ^= t;
       for (int Q=2; Q!=N; Q<<=1) {
         int R = Q - 1;
         for (int d=2; d>=0; d--) {
           if (iX[d] & Q) iX[0] ^= R;
           else {
             t = (iX[0] ^ iX[d]) & R;
             iX[0] ^= t;
             iX[d] ^= t;
           }
         }
       }
   #endif
       return iX;
     }
   
     ivec3 get3DIndex(vec3 X, int level, vec3 x0, real_t r0) {
       vec3 Xmin = x0 - r0;
       real_t dx = 2 * r0 / (1 << level);
       ivec3 iX;
       for (int d=0; d<3; d++) {
         iX[d] = floor((X[d] - Xmin[d]) / dx);
       }
       return iX;
     }
   
     vec3 getCoordinates(ivec3 iX, int level, vec3 x0, real_t r0) {
       vec3 Xmin = x0 - r0;
       real_t dx = 2 * r0 / (1 << level);
       vec3 X;
       for (int d=0; d<3; d++) {
         X[d] = (iX[d] + 0.5) * dx + Xmin[d];
       }
       return X;
     }
   }
   #endif