Typedefs
typedef unsigned int	byte

Functions
static void	calctables (int n)

void	hilbert_i2c (int n, int m, long int r, a)

void	hilbert_c2i (int n, int m, a, long int *r)

Variables
static char	rcsid [] = "$Header: /l/spencer/src/urt/lib/RCS/hilbert.c,v 3.0.1.1 1992/04/30 14:07:11 spencer Exp $"

static int	nbits = 0

static byte	p_to_s [512]

static byte	s_to_p [512]

static byte	p_to_J [512]

static byte	bit [9]

static byte	circshift [512][9]

static byte	parity [512]

static byte	bitof [512][9]

Typedef Documentation

typedef unsigned int byte

Definition at line 47 of file hilbert.c.

Function Documentation

static void calctables ( int n )

static

Definition at line 62 of file hilbert.c.

 {
     register int i, b;
     int two_n = 1 << n;
 
     if ( nbits == n )
     return;
 
     nbits = n;
     /* bit array is easy. */
     for ( b = 0; b < n; b++ )
     bit[b] = 1 << (n - b - 1);
 
     /* Next, do bitof. */
     for ( i = 0; i < two_n; i++ )
     for ( b = 0; b < n; b++ )
         bitof[i][b] = (i & bit[b]) ? 1 : 0;
 
     /* circshift is independent of the others. */
     for ( i = 0; i < two_n; i++ )
     for ( b = 0; b < n; b++ )
         circshift[i][b] = (i >> (b)) |
         ((i << (n - b)) & (two_n - 1));
 
     /* So is parity. */
     parity[0] = 0;
     for ( i = 1, b = 1; i < two_n; i++ )
     {
     /* Parity of i is opposite of the number you get when you
      * knock the high order bit off of i.
      */
     if ( i == b * 2 )
         b *= 2;
     parity[i] = !parity[i - b];
     }
 
     /* Now do p_to_s, s_to_p, and p_to_J. */
     for ( i = 0; i < two_n; i++ )
     {
     int s;
 
     s = i & bit[0];
     for ( b = 1; b < n; b++ )
         if ( bitof[i][b] ^ bitof[i][b-1] )
         s |= bit[b];
     p_to_s[i] = s;
     s_to_p[s] = i;
 
     p_to_J[i] = n - 1;
     for ( b = 0; b < n; b++ )
         if ( bitof[i][b] != bitof[i][n-1] )
         p_to_J[i] = b;
     }
 }

void hilbert_c2i	(	int	n,
		int	m,
		a	,
		long int *	r
	)

Definition at line 245 of file hilbert.c.

 {
     byte rho[9], J, sigma, tau,
          sigmaT, tauT, tauT1 = 0, omega, omega1 = 0, alpha[9];
     register int i, b;
     int Jsum;
     long int rl;
 
     calctables(n);
 
     /* Unpack the coordinates into alpha[i]. */
     /* First, zero out the alphas. */
     switch ( m ) {
     case 9: alpha[8] = 0;
     case 8: alpha[7] = 0;
     case 7: alpha[6] = 0;
     case 6: alpha[5] = 0;
     case 5: alpha[4] = 0;
     case 4: alpha[3] = 0;
     case 3: alpha[2] = 0;
     case 2: alpha[1] = 0;
     case 1: alpha[0] = 0;
     }
 
     /* The loop that unpacks bits of a[b] into alpha[i] has been unrolled.
      * The high-order bit of a[b] has to go into alpha[0], so we pre-shift
      * a[b] so that its high-order bit is always in the 0x100 position.
      */
     for ( b = 0; b < n; b++ )
     {
     register int bt = bit[b], t = a[b] << (9 - m);
 
     switch (m)
     {
     case 9: if ( t & 0x01 ) alpha[8] |= bt;
     case 8: if ( t & 0x02 ) alpha[7] |= bt;
     case 7: if ( t & 0x04 ) alpha[6] |= bt;
     case 6: if ( t & 0x08 ) alpha[5] |= bt;
     case 5: if ( t & 0x10 ) alpha[4] |= bt;
     case 4: if ( t & 0x20 ) alpha[3] |= bt;
     case 3: if ( t & 0x40 ) alpha[2] |= bt;
     case 2: if ( t & 0x80 ) alpha[1] |= bt;
     case 1: if ( t & 0x100 ) alpha[0] |= bt;
     }
     }
 
     Jsum = 0;
     for ( i = 0; i < m; i++ )
     {
     /* Compute omega[i] = omega[i-1] xor tauT[i-1]. */
     if ( i == 0 )
         omega = 0;
     else
         omega = omega1 ^ tauT1;
 
     sigmaT = alpha[i] ^ omega;
     /* sigma[i] is the left circular shift of sigmaT[i]. */
     if ( Jsum != 0 )
         sigma = circshift[sigmaT][n - Jsum];
     else
         sigma = sigmaT;
 
     rho[i] = s_to_p[sigma];
 
     /* Now we can get the principle position. */
     J = p_to_J[rho[i]];
 
     /* And compute tau[i] and tauT[i]. */
     /* tau[i] complements low bit of sigma[i], and bit at J[i] if
      * necessary to make even parity.
      */
     tau = sigma ^ 1;
     if ( parity[tau] )
         tau ^= bit[J];
 
     /* tauT[i] is right circular shift of tau[i]. */
     if ( Jsum != 0 )
         tauT = circshift[tau][Jsum];
     else
         tauT = tau;
     Jsum += J;
     if ( Jsum >= n )
         Jsum -= n;
 
     /* Carry forth the "i-1" values. */
     tauT1 = tauT;
     omega1 = omega;
     }
 
     /* Pack rho values into r. */
     rl = 0;
     for ( i = 0; i < m; i++ )
     rl = (rl << n) | rho[i];
     *r = rl;
 }

void hilbert_i2c	(	int	n,
		int	m,
		long int	r,
		a
	)

Definition at line 136 of file hilbert.c.

 {
     byte rho[9], rh, J, sigma, tau,
          sigmaT, tauT, tauT1 = 0, omega, omega1 = 0, alpha[9];
     register int i, b;
     int Jsum;
 
     /* Initialize bit twiddle tables. */
     calctables(n);
 
     /* Distribute bits from r into rho. */
     for ( i = m - 1; i >= 0; i-- )
     {
     rho[i] = r & ((1 << n) - 1);
     r >>= n;
     }
 
     /* Loop over bytes. */
     Jsum = 0;
     for ( i = 0; i < m; i++ )
     {
     rh = rho[i];
     /* J[i] is principle position of rho[i]. */
     J = p_to_J[rh];
 
     /* sigma[i] is derived from rho[i] by exclusive-oring adjacent bits. */
     sigma = p_to_s[rh];
 
     /* tau[i] complements low bit of sigma[i], and bit at J[i] if
      * necessary to make even parity.
      */
     tau = sigma ^ 1;
     if ( parity[tau] )
         tau ^= bit[J];
     
     /* sigmaT[i] is circular shift of sigma[i] by sum of J[0..i-1] */
     /* tauT[i] is same circular shift of tau[i]. */
     if ( Jsum > 0 )
     {
         sigmaT = circshift[sigma][Jsum];
         tauT = circshift[tau][Jsum];
     }
     else
     {
         sigmaT = sigma;
         tauT = tau;
     }
 
     Jsum += J;
     if ( Jsum >= n )
         Jsum -= n;
 
     /* omega[i] is xor of omega[i-1] and tauT[i-1]. */
     if ( i == 0 )
         omega = 0;
     else
         omega = omega1 ^ tauT1;
     omega1 = omega;
     tauT1 = tauT;
 
     /* alpha[i] is xor of omega[i] and sigmaT[i] */
     alpha[i] = omega ^ sigmaT;
     }
 
     /* Build coordinates by taking bits from alphas. */
     for ( b = 0; b < n; b++ )
     {
     register int ab, bt;
     ab = 0;
     bt = bit[b];
     /* Unroll the loop that stuffs bits into ab.
      * The result is shifted left by 9-m bits.
      */
     switch( m )
     {
     case 9: if ( alpha[8] & bt) ab |= 0x01;
     case 8: if ( alpha[7] & bt) ab |= 0x02;
     case 7: if ( alpha[6] & bt) ab |= 0x04;
     case 6: if ( alpha[5] & bt) ab |= 0x08;
     case 5: if ( alpha[4] & bt) ab |= 0x10;
     case 4: if ( alpha[3] & bt) ab |= 0x20;
     case 3: if ( alpha[2] & bt) ab |= 0x40;
     case 2: if ( alpha[1] & bt) ab |= 0x80;
     case 1: if ( alpha[0] & bt) ab |= 0x100;
     }
     a[b] = ab >> (9 - m);
     }
 }