/*
 * This software is copyrighted as noted below.  It may be freely copied,
 * modified, and redistributed, provided that the copyright notice is 
 * preserved on all copies.
 * 
 * There is no warranty or other guarantee of fitness for this software,
 * it is provided solely "as is".  Bug reports or fixes may be sent
 * to the author, who may or may not act on them as he desires.
 *
 * You may not include this software in a program or other software product
 * without supplying the source, or without informing the end-user that the 
 * source is available for no extra charge.
 *
 * If you modify this software, you should include a notice giving the
 * name of the person performing the modification, the date of modification,
 * and the reason for such modification.
 */
/* 
 * hilbert.c - Computes Hilbert curve coordinates from position and v.v.
 * 
 * Author:	Spencer W. Thomas
 * 		EECS Dept.
 * 		University of Michigan
 * Date:	Thu Feb  7 1991
 * Copyright (c) 1991, University of Michigan
 */
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 $";

/*
 * Lots of tables to simplify calculations.  Notation: p#i means bit i
 * in byte p (high order bit first).
 * p_to_s:	Output s is a byte from input p such that
 * 		s#i = p#i xor p#(i-1)
 * s_to_p:	The inverse of the above.
 * p_to_J:	Output J is "principle position" of input p.  The
 * 		principle position is the last bit s.t.
 * 		p#J != p#(n-1) (or n-1 if all bits are equal).
 * bit:		bit[i] == (1 << (n - i))
 * circshift:	circshift[b][i] is a right circular shift of b by i
 * 		bits in n bits.
 * parity:	Parity[i] is 1 or 0 depending on the parity of i (1 is odd).
 * bitof:	bitof[b][i] is b#i.
 * nbits:	The value of n for which the above tables have been
 * 		calculated.
 */

typedef unsigned int byte;

static int nbits = 0;

static byte
    p_to_s[512],
    s_to_p[512],
    p_to_J[512],
    bit[9],
    circshift[512][9],
    parity[512],
    bitof[512][9];

/* Calculate the above tables when nbits changes. */
static void
calctables(n)
int n;
{
    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;
    }
}

/*****************************************************************
 * TAG( hilbert_i2c )
 * 
 * Convert an index into a Hilbert curve to a set of coordinates.
 * Inputs:
 * 	n:	Number of coordinate axes.
 * 	m:	Number of bits per axis.
 * 	r:	The index, contains n*m bits (so n*m must be <= 32).
 * Outputs:
 * 	a:	The list of n coordinates, each with m bits.
 * Assumptions:
 * 	n*m < (sizeof r) * (bits_per_byte), n <= 8, m <= 9.
 * Algorithm:
 * 	From A. R. Butz, "Alternative Algorithm for Hilbert's
 * 		Space-Filling Curve", IEEE Trans. Comp., April, 1971,
 * 		pp 424-426.
 */
void
hilbert_i2c( n, m, r, a)
int n;
int m;
long int r;
int a[];
{
    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);
    }
}

/*****************************************************************
 * TAG( hilbert_c2i )
 * 
 * Convert coordinates of a point on a Hilbert curve to its index.
 * Inputs:
 * 	n:	Number of coordinates.
 * 	m:	Number of bits/coordinate.
 * 	a:	Array of n m-bit coordinates.
 * Outputs:
 * 	r:	Output index value.  n*m bits.
 * Assumptions:
 * 	n*m <= 32, n <= 8, m <= 9.
 * Algorithm:
 * 	Invert the above.
 */
void
hilbert_c2i( n, m, a, r)
int n;
int m;
int a[];
long int *r;
{
    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;
}


#ifdef test
#include <stdio.h>

main()
{
    int a[9], n, m, i;
    long int r, r1;

    for (;;)
    {
	printf( "Enter n, m: " );
	scanf( "%d", &n );
	if ( n == 0 )
	    break;
	scanf( "%d", &m );
	while ( (i = getchar()) != '\n' && i != EOF )
	    ;
	if ( i == EOF )
	    break;
	for ( r = 0; r < 1 << (n*m); r++ )
	{
	    hilbert_i2c( n, m, r, a );
	    if ( n*m <= 6 )
	    {
		printf( "a = " );
		for ( i = 0; i < n; i++ )
		    printf( "0x%0*x ", (m+3)/4, a[i] );
	    }
	    hilbert_c2i( n, m, a, &r1 );
	    if ( n*m <= 6 )
		printf( "r = 0x%0*x\n", (n*m+3)/4, r1 );
	    else if ( r != r1 )
		printf( "r = 0x%0*x; r1 = 0x%0*x\n", (n*m+3)/4, r,
			(n*m+3)/4, r1 );
	}
    }
}

p1t( tbl, n )
byte tbl[];
int n;
{
    int i;

    for ( i = 0; i < n; i++ )
	printf( "%02x ", tbl[i] );
    putchar( '\n' );
}

p2t( tbl, n, m )
byte tbl[][9];
int n;
{
    int i;

    for ( i = 0; i < n; i++ )
    {
	printf( "%3d: ", i );
	p1t( tbl[i], m );
    }
}


#endif