How to call Hilbert curve encode C routines

Question

I was trying to run a Hilbert curve code written in C, which I found here http://www.tddft.org/svn/octopus/trunk/src/grid/hilbert.c

The code runs, but the results I get form the output are not correct. I made a simple driver routine, that takes 3 values as arguments from the command line and passes them to a Hilbert curve encode, decode routines. More precisely, I can't decode back the original coordinates (x,y,z).

One of my problems was to understand what the variable nbits is doing. I assume it is the size of the encoded Hilbert value. To check this I tried to modify the original definition of one of the functions from

void InttoTranspose(const int dim, const long long int h, int * x)

to

void InttoTranspose(const int dim, const long long int h, int * x, int* size)

Where I assign to *size the bit count variable ifbit. Anyway, this all didn't work. Therefore I would like to ask for your help.

The modified code is here:

#include<stdio.h>

//#include <config.h>

#include <assert.h>

/* This code is based on the implementation of the Hilbert curves
   presented in:

   J. Skilling, Programming the Hilbert curve, AIP Conf. Proc. 707, 381 (2004); http://dx.doi.org/10.1063/1.1751381

*/

*/ The int* size was included later in an attempt to find the proper size  /*
/* void InttoTranspose(const int dim, const long long int h, int * x)*/

void InttoTranspose(const int dim, const long long int h, int * x, int* size){
  /* the code uses some funny way of storing the bits */

  int idir, ibit, ifbit;

  for(idir = 0; idir < dim; idir++) x[idir] = 0;

  ifbit = 0;
  for(ibit = 0; ibit < 21; ibit++){
    for(idir = dim - 1; idir >= 0; idir--){
      x[idir] += (((h>>ifbit)&1)<<ibit);
      ifbit++;
    }
  }
 *size=ifbit; // I think that this should be nbits
}


void TransposetoAxes(int* X, int b, int n ){ /* position, #bits, dimension */

  int N = 2 << (b-1), P, Q, t;
  int i;

  /* Gray decode by H ^ (H/2) */
  t = X[n-1] >> 1;
  for(i = n - 1; i > 0; i--) X[i] ^= X[i-1];
  X[0] ^= t;

  /* Undo excess work */
  for( Q = 2; Q != N; Q <<= 1 ) {
    P = Q - 1;
    for( i = n-1; i >= 0 ; i-- ){
      if( X[i] & Q ) {
    X[0] ^= P; /* invert */
      } else{ 
    t = (X[0]^X[i]) & P; X[0] ^= t; X[i] ^= t; /* exchange */
      }
    }
  }

}

//void FC_FUNC_(hilbert_index_to_point, HILBERT_INDEX_TO_POINT)(const int * dim, const int * nbits, const long long int * index, int * point){

 // InttoTranspose(*dim, *index, point);
 // TransposetoAxes(point, *nbits, *dim);

//}


int main(int argc,char* argv[]){

  long long int hilbert;
  int i=0,x[2],m;

  while(argc--){
    if(m=atoi(*argv++)) x[i++]=m, printf("--> %5d %5d\n",m,argc);
  }


  printf("x= %5d y= %5d  z= %5d \n",x[0],x[1],x[2]);

  InttoTranspose(3, hilbert, x,&m);
  printf("hilbert encoded --> %llu size -->%d \n",hilbert,m);
  TransposetoAxes(x,m, 3 );
  printf("x= %5d y= %5d  z= %5d \n",x[0],x[1],x[2]);

  return 0;
}