📄 hilbert.c
字号:
coord_t* X, // I Transpose [n]
int b, // I # bits
int n) // I dimension
{
coord_t j, p, M;
int i, q;
M = 1 << (b - 1);
q = 0; p = M;
for( i = 0; i < n; i++ )
{
Line[i] = 0;
for( j = M; j; j >>= 1 )
{
if( X[q] & p )
Line[i] |= j;
if( ++q == n )
{
q = 0; p >>= 1;
}
}
}
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Functions: TransposetoAxes
// AxestoTranspose
//
// Purpose: Transform between Hilbert transpose and geometrical axes
//
// Example: b=5 bits for each of n=3 coordinates
// Hilbert transpose
// X[0] = A D b e 3 X[1]|
// X[1] = B E c 1 4 <-------> | /X[2]
// X[2] = C a d 2 5 axes | /
// high low |/______
// X[0]
// Axes are stored conventially as b-bit integers.
//
// History: John Skilling 20 Apr 2001, 3 Sep 2003, 14 Oct 2003
//-----------------------------------------------------------------------------
//
static void TransposetoAxes(
coord_t* X, // I O position [n]
int b, // I # bits
int n) // I dimension
{
coord_t M, P, Q, t;
int i;
// Gray decode by H ^ (H/2)
t = X[n-1] >> 1;
for( i = n-1; i; i-- )
X[i] ^= X[i-1];
X[0] ^= t;
// Undo excess work
M = 2 << (b - 1);
for( Q = 2; Q != M; Q <<= 1 )
{
P = Q - 1;
for( i = n-1; i; i-- )
if( X[i] & Q ) X[0] ^= P; // invert
else{ t = (X[0] ^ X[i]) & P; X[0] ^= t; X[i] ^= t; } // exchange
if( X[0] & Q ) X[0] ^= P; // invert
}
}
static void AxestoTranspose(
coord_t* X, // I O position [n]
int b, // I # bits
int n) // I dimension
{
coord_t P, Q, t;
int i;
// Inverse undo
for( Q = 1 << (b - 1); Q > 1; Q >>= 1 )
{
P = Q - 1;
if( X[0] & Q ) X[0] ^= P; // invert
for( i = 1; i < n; i++ )
if( X[i] & Q ) X[0] ^= P; // invert
else{ t = (X[0] ^ X[i]) & P; X[0] ^= t; X[i] ^= t; } // exchange
}
// Gray encode (inverse of decode)
for( i = 1; i < n; i++ )
X[i] ^= X[i-1];
t = X[n-1];
for( i = 1; i < b; i <<= 1 )
X[n-1] ^= X[n-1] >> i;
t ^= X[n-1];
for( i = n-2; i >= 0; i-- )
X[i] ^= t;
}
//=============================================================================
// Linked-list library
//=============================================================================
// atom___________________________________
// Base| | | Ptr | Ptr | Ptr | Ptr |
// Free| Ptr | ... | Ptr | Ptr | Ptr | Ptr |
// |_____|_____|_____|_____|_____|_____|
// / / / / / / / / /
// 0 / / / / / / / / /
// : / / / / / / / / /
// : / vacant nodes
// ______/
// | Node |
// |depth3|
// |______|
// /: \.
// /: \.
// /: \.
// /: \.
// /: \.
// /: \.
// /: \.
// /: ______
// /: | Node |
// /: |depth2|
// /: |______|
// /: /: \.
// /: /: \.
// ______ /: ______
// | Node | /: | Node |
// |depth1| /: |depth1|
// |______| /: |______|
// /: \ /: /: \.
// /: \ /: /: \.
// /: \ /: /: \.
// Base______ ______ ______ ______ ______
// | Node | | Node | | Node | | Node | | Node |
// 0----|depth0|----|depth0|----|depth0|----|depth0|----|depth0|----0
// |______| |______| |______| |______| |______|
// : : : : :
// : : : : :
// : : : : :
// atom_____ _____ _____ _____ _____
// |extra| |Atom2| |Atom0| |Atom3| |Atom1|
// | 0 | |Label| |Label| |Label| |Label|
// |atom | | etc.| | etc.| | etc.| | etc.|
// |_____| |_____| |_____| |_____| |_____|
// <---------- random-order storage --------->
//
// 0 <= < < < < < < Label < < < < <
//
// Inserted atoms have serial numbers 0, 1, 2, ..., N-1 with labels x >= 0,
// and occupy consecutive storage. To maintain this efficiently under random
// insertion and deletion demands that they be stored in dynamically changing
// order. Hence coupling to their neighbours in x must be by a linked list,
// and to avoid ambiguity all labels must be different.
//
// The binary tree allows insertion of a new atom with known label x,
// and deletion of a random atom, at a cost that is asymptotically
// logarithmic O(log N) because of the search.
// After each insertion, the newly inserted atom is at the end of storage.
// and after each deletion, the newly deleted atom is just beyond the end.
// After each insertion or deletion, the tree is re-balanced to keep its
// maximum depth close to log2(N).
//
// There are 2N+1 nodes of the tree, which only occupy consecutive storage
// when N is at a historic maximum: when N falls away through deletions,
// holes appear in memory. These are tracked and re-used by a consecutive
// vector of controlling pointers (which is conveniently stored in the
// unused atoms to save memory).
//
// Memory = sizeof(Atom) + 2 * sizeof(Node), per atom.
//
// CPU = asymptotically O(log N) per insertion or deletion.
// Tests with a range of practical N give roughly
// =======================================
// CPU = 10 N (log10(N) + 5) multiply-adds per (insertion and deletion) cycle.
// =======================================
// (Compare 10 N log10(N) multiply-adds for a single FFT.)
//
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: SetLink
//
// Purpose: Allocate and initialise memory,
// using guess for maximum number of atoms needed.
//
// See also ResetLink and FreeLink
//
// History: John Skilling 9 Mar 1996, 6 Feb 1998, 12 Apr 2001, 10 Sep 2001,
// 22 Jan 2003
//=============================================================================
// STRUCTURE OF ATOMS
// Node = apex of linked list
// \.
// \.
// \.
// List \______
// 0 | Label| (unused)
// | flux | (unused)
// | base | Link up to origin base of tree
// | free | Point to apex of tree
// |______|
// 1 | Label| x ./|\.
// | flux | extra properties of atom |
// | base | Link up to base of tree |
// | free | (unused) |
// |______| | Surviving atoms
// | | |
// ........ | in scrambled
// |______| |
// Natoms | | | storage order
// | .... | |
// | | |
// | .... | |
// |______| \|/
// Natoms+1 | Label| (unused)
// | flux | (unused)
// | base | Point to vacant node address
// | free | Point to vacant node address
// |______|
// | |
// ........
// |______|
// NMAX | |
// | .... |
// | |
// | .... |
// |______| Reallocate memory if Natoms attempts to exceed NMAX
//
//=============================================================================
// STRUCTURE OF NODES
//
// Link node ______
// | depth| >= 1
// |number| >= 2
// |parent| Link up to parent (NULL if apex)
// | lo | Link down to "next" child
// | hi | Link down to "previous" child
// | atom | Point down "previous" branches to atom
// |______|
//
// Base node ______
// | depth| 0
// |number| 1
// |parent| Link up to parent (NULL if apex of empty tree)
// | lo | Link across to "next" neighbour (NULL if last)
// | hi | Link across to "previous" neighbour (NULL if first)
// | atom | Link down to atom
// |______|
//
//-----------------------------------------------------------------------------
//
int SetLink( // O 0, or -ve error code
int Ndim, // I # words per label
Node* psLink) // I O Linked list
{
int NMAX = 100; // Initial number of allowed insertions of an atom
pNode Nodes = NULL; // Begin tree storage
pAtom List = NULL; // List of atoms
int i; // Local counter
int CALLvalue = 0;
psLink->atom = NULL;
psLink->parent = NULL;
// Allocate storage for tree
CALLOC( Nodes, NMAX+NMAX+1, Node )
CALLOC( List, NMAX+1, Atom )
CALLOC( List->Label, Ndim * (NMAX+1), coord_t )
// Initialise artificial header atom
List->Base = Nodes;
List->Free = Nodes;
// Initialise pointers to available nodes
for( i = 1; i <= NMAX; ++i )
{
List[i].Base = &Nodes[i+i-1];
List[i].Free = &Nodes[i+i];
List[i].Label = List[i-1].Label + Ndim;
}
// Initialise top of tree
Nodes->depth = 0;
Nodes->number = 1;
Nodes->parent = NULL;
Nodes->hi = NULL;
Nodes->lo = NULL;
Nodes->atom = List;
psLink->depth = NMAX;
psLink->number = Ndim;
psLink->atom = List;
psLink->parent = Nodes;
return 0;
Exit:
FREE( List->Label )
FREE( List )
FREE( Nodes )
return CALLvalue;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: FreeLink
//
// Purpose: Free memory allocated by SetLink.
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -