希尔伯特变换代码

/* See LICENSE below for information on rights to use modify and distribute
   this code. */

/* 
 * hilbert.c - Computes Hilbert space-filling curve coordinates without
 * recursion from integer index and vice versa and other Hilbert-related
 * calculations.  Also known as Pi-order or Peano scan.
 * 
 * Author:      Doug Moore
 *              Dept. of Computational and Applied Math
 *              Rice University
 *              http://www.caam.rice.edu/~dougm
 * Date:        Sun Feb 20 2000
 * Copyright （c） 1998-2000 Rice University
 *
 * Acknowledgement:
 * This implementation is based on the work of A. R. Butz （“Alternative
 * Algorithm for Hilbert‘s Space-Filling Curve“ IEEE Trans. Comp. April
 * 1971 pp 424-426） and its interpretation by Spencer W. Thomas University
 * of Michigan （http://www-personal.umich.edu/~spencer/Home.html） in his widely
 * available C software.  While the implementation here differs considerably
 * from his the first two interfaces and the style of some comments are very
 * much derived from his work. */


#include “hilbert.h“

/* implementation of the hilbert functions */

#define adjust_rotation（rotationnDimsbits）                            \
do {                                                                    \
      /* rotation = （rotation + 1 + ffs（bits）） % nDims; */              \
      bits &= -bits & nd1Ones;                                          \
      while （bits）                                                      \
        bits >>= 1 ++rotation;                                         \
      if （ ++rotation >= nDims ）                                        \
        rotation -= nDims;                                              \
} while （0）

#define ones（Tk） （（（（T）2） << （k-1）） - 1）

#define rdbit（wk） （（（w） >> （k）） & 1）
     
#define rotateRight（arg nRots nDims）                                  \
（（（（arg） >> （nRots）） | （（arg） << （（nDims）-（nRots）））） & ones（bitmask_tnDims））

#define rotateLeft（arg nRots nDims）                                   \
（（（（arg） << （nRots）） | （（arg） >> （（nDims）-（nRots）））） & ones（bitmask_tnDims））

#define DLOGB_BIT_TRANSPOSE
static bitmask_t
bitTranspose（unsigned nDims unsigned nBits bitmask_t inCoords）
#if defined（DLOGB_BIT_TRANSPOSE）
{
  unsigned const nDims1 = nDims-1;
  unsigned inB = nBits;
  unsigned utB;
  bitmask_t inFieldEnds = 1;
  bitmask_t inMask = ones（bitmask_tinB）;
  bitmask_t coords = 0;

  while （（utB = inB / 2））
    {
      unsigned const shiftAmt = nDims1 * utB;
      bitmask_t const utFieldEnds =
	inFieldEnds | （inFieldEnds << （shiftAmt+utB））;
      bitmask_t const utMask =
	（utFieldEnds << utB） - utFieldEnds;
      bitmask_t utCoords = 0;
      unsigned d;
      if （inB & 1）
	{
	  bitmask_t const inFieldStarts = inFieldEnds << （inB-1）;
	  unsigned oddShift = 2*shiftAmt;
	  for （d = 0; d < nDims; ++d）
	    {
	      bitmask_t in = inCoords & inMask;
	      inCoords >>= inB;
	      coords |= （in & inFieldStarts） <<	oddShift++;