📄 hilbert.c
字号:
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Filename: hilbert.c
//
// Purpose: Hilbert and Linked-list utility procedures for BayeSys3.
//
// History: TreeSys.c 17 Apr 1996 - 31 Dec 2002
// Peano.c 10 Apr 2001 - 11 Jan 2003
// merged 1 Feb 2003
// Arith debug 28 Aug 2003
// Hilbert.c 14 Oct 2003
// 2 Dec 2003
//-----------------------------------------------------------------------------
/*
Copyright (c) 1996-2003 Maximum Entropy Data Consultants Ltd,
114c Milton Road, Cambridge CB4 1XE, England
This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.
This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Lesser General Public License for more details.
You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
#include "license.txt"
*/
#include <stdlib.h>
#include "hilbert.h"
typedef Atom* pAtom; // Atom address
typedef Node* pNode; // Node address
// Internal prototypes
static void LinetoTranspose(coord_t*, coord_t*, int, int);
static void TransposetoLine(coord_t*, coord_t*, int, int);
static void TransposetoAxes(coord_t*, int, int);
static void AxestoTranspose(coord_t*, int, int);
static int ResetLink (pNode);
static void Balance (pNode);
// Internal macros
#undef CALLOC // allocates vector p[0:n-1] of type t
#define CALLOC(p,n,t) {p=NULL;\
if((n)>0&&!(p=(t*)calloc((size_t)(n),sizeof(t))))\
{CALLvalue=E_MALLOC;goto Exit;}/*printf("%p %d\n",p,(size_t)(n)*sizeof(t));*/}
#undef FREE // frees CALLOC or REALLOC or NULL vector p[0:*], sets p=NULL
#define FREE(p) {if(p){/*printf("%p -1\n",p);*/(void)free((void*)p);} p=NULL;}
#undef CALL // catches negative error codes
#define CALL(x) {if( (CALLvalue = (x)) < 0 ) goto Exit;}
//=============================================================================
// Composite-integer arithmetic library
//=============================================================================
//
// A composite-integer is a multi-word unsigned integer "Label" stored
// "big endian" in N conventional unsigned integers with [0] high.
// ___________________________________________________
// | | | | |
// | Label[0] | Label[1] | .... | Label[N-1] |
// |____________|____________|____________|____________|
// high low
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: CmpLabel
//
// Purpose: Compare labels by generating
// +1 if u > v,
// sign(u - v) = 0 if u = v,
// -1 if u < v.
//
// History: John Skilling 12 Apr 2001
//-----------------------------------------------------------------------------
//
int CmpLabel( // O comparison
coord_t* u, // I composite integer ([0] high)
coord_t* v, // I composite integer ([0] high)
int Ndim) // I dimension
{
int j;
for( j = 0; j < Ndim; j++ )
if( u[j] < v[j] )
return -1;
else if( u[j] > v[j] )
return 1;
return 0;
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: AddLabel
//
// Purpose: Set w = u + v
//
// History: JS 28 Jan 2002, 31 Dec 2002
// Julian Center 28 Aug 2003 debug
//-----------------------------------------------------------------------------
//
void AddLabel(
coord_t* w, // O w = u + v [Ndim]
coord_t* u, // I can be overwritten [Ndim]
coord_t* v, // I must not be overwritten [Ndim]
int Ndim) // I dimension
{
int carry = 0;
int i;
for( i = Ndim-1; i >= 0; i-- )
{
w[i] = u[i] + v[i];
carry = carry ? (++w[i] <= v[i]) : (w[i] < v[i]);
}
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Function: SubLabel
//
// Purpose: Set w = u - v
//
// History: JS 28 Jan 2002, 31 Dec 2002
// Julian Center 28 Aug 2003 debug
//-----------------------------------------------------------------------------
//
void SubLabel(
coord_t* w, // O w = u - v [Ndim]
coord_t* u, // I can be overwritten [Ndim]
coord_t* v, // I must not be overwritten [Ndim]
int Ndim) // I dimension
{
int carry = 0;
int i;
for( i = Ndim-1; i >= 0; i-- )
{
w[i] = u[i] - v[i];
carry = carry ? (--w[i] >= u[i]) : (w[i] > u[i]);
}
}
//=============================================================================
// Hilbert-curve (a space-filling Peano curve) library
//=============================================================================
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Functions: LinetoAxes
// AxestoLine
//
// Purpose: Serial Hilbert length <----> multidimensional Axes position.
//
// Space = n-dimensional hypercube of side R = 2^b
// Number of cells = N = R^n = 2^(n*b)
//
// Line = serial number of cell along Hilbert curve through hypercube
// = extended integer of n*b bits ranging from 0 to N-1,
// stored as vector of n unsigned b-bit integers with [0] high.
//
// Axes = Geometrical position of cell
// = n b-bit integers representing coordinates.
//
// Example: side R = 16, dimension n = 2, number of cells = N = 256.
// Line = 9, stored in base-16 words as
// Line[0] = 0 (high), Line[1] = 9 (low),
// corresponds to position (2,3) as in diagram, stored as
// Axes[0] = 2, Axes[1] = 3.
//
// |
// 15 | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@
// | | | | | | | | | | | | | | | | |
// | @ @---@ @ @ @---@ @ @ @---@ @ @ @---@ @
// | | | | | | | | |
// | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@
// | | | | | | | | |
// | @---@ @---@---@---@ @---@ @---@ @---@---@---@ @---@
// | | | | |
// | @ @---@---@ @---@---@ @ @ @---@---@ @---@---@ @
// | | | | | | | | | | | | |
// Axes[1]| @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@
// | | | | |
// | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@
// | | | | | | | | | | | | |
// | @ @---@---@ @---@---@ @---@ @---@---@ @---@---@ @
// | | |
// | @---@ @---@---@ @---@---@ @---@---@ @---@---@ @---@
// | | | | | | | | | | |
// | @---@ @---@ @---@ @---@ @---@ @---@ @---@ @---@
// | | | | | | |
// | @ @---@ @ @---@ @---@ @---@ @---@ @ @---@ @
// | | | | | | | | | | | | | | |
// | @---@ @---@ @ @---@---@ @---@---@ @ @---@ @---@
// | | |
// 3 | 5---6 9---@ @ @---@---@ @---@---@ @ @---@ @---@
// | | | | | | | | | | | | | | |
// 2 | 4 7---8 @ @---@ @---@ @---@ @---@ @ @---@ @
// | | | | | | |
// 1 | 3---2 @---@ @---@ @---@ @---@ @---@ @---@ @---@
// | | | | | | | | | | |
// 0 | 0---1 @---@---@ @---@---@ @---@---@ @---@---@ @--255
// |
// -------------------------------------------------------------------
// 0 1 2 3 ---> Axes[0] 15
//
// Notes: (1) Unit change in Line yields single unit change in Axes position:
// the Hilbert curve is maximally local.
// (2) CPU proportional to total number of bits, = b * n.
//
// History: John Skilling 20 Apr 2001, 11 Jan 2003, 3 Sep 2003
//-----------------------------------------------------------------------------
//
void LinetoAxes(
coord_t* Axes, // O multidimensional geometrical axes [n]
coord_t* Line, // I linear serial number, stored as [n]
int b, // I # bits used in each word
int n) // I dimension
{
if( n <= 1 ) // trivial case
*Axes = *Line;
else
{
LinetoTranspose(Axes, Line, b, n);
TransposetoAxes(Axes, b, n);
}
}
void AxestoLine(
coord_t* Line, // O linear serial number, stored as [n]
coord_t* Axes, // I multidimensional geometrical axes [n]
int b, // I # bits used in each word
int n) // I dimension
{
coord_t store[1024]; // avoid overwriting Axes
int i; // counter
if( n <= 1 ) // trivial case
*Line = *Axes;
else if( n <= 1024 ) // surely the usual case
{
for( i = 0; i < n; ++i )
store[i] = Axes[i];
AxestoTranspose( store, b, n);
TransposetoLine(Line, store, b, n);
}
else // must do in place at greater cost
{
AxestoTranspose( Axes, b, n);
TransposetoLine(Line, Axes, b, n);
TransposetoAxes( Axes, b, n);
}
}
//+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
// Functions: LinetoTranspose
// TransposetoLine
//
// Purpose: Recover Hilbert integer by bit-transposition
//
// Example: b=5 bits for each of n=3 coordinates
// 15-bit Hilbert integer = A B C D E a b c d e 1 2 3 4 5
// X[0]..... X[1]..... X[2].....
// transposed to
// X[0](high) = A D b e 3
// X[1] = B E c 1 4
// X[2](low) = C a d 2 5
// high low
//
// History: John Skilling 20 Apr 2001, 3 Sep 2003, 14 Oct 2003
//-----------------------------------------------------------------------------
//
static void LinetoTranspose(
coord_t* X, // O Transpose [n]
coord_t* Line, // I Hilbert integer [n]
int b, // I # bits
int n) // I dimension
{
coord_t j, p, M;
int i, q;
M = 1 << (b - 1);
for( i = 0; i < n; i++ )
X[i] = 0;
q = 0; p = M;
for( i = 0; i < n; i++ )
{
for( j = M; j; j >>= 1 )
{
if( Line[i] & j )
X[q] |= p;
if( ++q == n )
{
q = 0; p >>= 1;
}
}
}
}
static void TransposetoLine(
coord_t* Line, // O Hilbert integer [n]
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -