(*
LITMAX/BIGMIN computation  Pascal source code
Author: Hermann Tropf
Date: Feb 22, 2021.
Realisation hints updated: Mar 30, 2021.
This program refers to the article "Multidimensional Range Search in Dynamically Balanced Trees",
by Hermann Tropf and Helmut Herzog, Angewandte Informatik (applied informatics) Febr. 1981 pp 7177.
Download: www.visiontools.com/htropf/multidimensionalrangequery.pdf
With this post, I hope to provide a more vivid explanation of the key functions LITMAX / BIGMIN
described in the article. I try it by means of the old fashioned character graphics below.
The above article describes the use of bitwise interleaved data for efficient multidimensional range
searching in dynamic data. A key idea of the method is the usage of two symmetric functions
called LITMAX and BIGMIN to which I refer in the following.
40 Years after publication of the article, I researched for its impact and I found that the method
has found its use in some commercial databases and in a bunch of different technical application
areas  and who knows where else without being mentioned.
However, I have the impression that some people had difficulties in implementing the LITMAX/BIGMIN
functions described in the article. Assumedly some of them have resorted to other methods that are
conceptually simpler at first sight. This is a pity, because in contrast to other methods like
kdtrees or Range trees, the approach is independent of the data organisation chosen. You can use
any 1dimensional data organisation (binary search trees, skip lists, 1Drange trees, whatever).
Even simple sorted arrays can be used, skipping large portions of the data by use of LITMAX and/or
BIGMIN. Due to this independency it is easier to integrate the method into existing data bases or
technical applications.
Realising that the LITMAX/BIGMIN calculation seems to be a major hurd when going into using the
method, I decided to present the LITMAX/BIGMIN calculation directly as Pascal source code,
together with a somewhat more elaborate explanation of how it works.
At some time I worked with Hilbert Order instead of Zorder and realised that Zorder bit
interleaving is not necessarily done explicitly (US Patent 7321890 , section "Technical Remarks",
Patent dropped). Instead, data are kept as they are, they are just processed in bit interleaved
manner. So there are no problems with word length restrictions, and the code can easily be
adapted to any number of dimensions and to any key wordlength.
And, for this presentation, keeping the data as they are and processing them instead in interleaved
order, makes the explanation of how LITMAX/BIGMIN works somewhat easier.
For a low number of dimensions and short data wordlength, it may still remain advantegeous
to interleave the data, especially if ParallelBitsDeposit instructions or the like are available.
To realise any structure for the zordered Data (search trees, skip lists, sorted arrays or
whatever), you need a function comparing the Zvalue of two records for the given number of
dimensions. It is also given below (function "zvalue_compare").
Some realisation hints (including usage of floating point data) are given at the end.
Before going into the source code, at first a visualisation of LITMAX by means of an example in 2D:
LITMAX, how it works

See Fig. 6 in the article.
The data field:

Bitwise interleaving the 2D data (x,y) effects that the entire point field is alternatingly
devided up in x and y direction.
The first division is in x direchtion and splits it into a low xsection and a high xsection
These sections are divided in y direchtion; the sections are all split into a low ysection
and a high ysection
These sections are divided in x direchtion; the sections are all split into a low xsection
and a high xsection
etc etc
> x
   
   
v   
y                        
  
  
  

                       

  
  
  
                      
  
  
  
ZValues

The Zvalue of a data point is the value of its coordinates when x and y are bitwise
interleaved. As mentioned above, it is not necessary to interleave the date explicitly.
There is only needed a function to decide which of two Points has the greater Zvalue.
LITMAX calculation

F is the point found in the data
****
* * is the search rectangle
****
I=Min is the corner of the search rectangle with minimum data (upper left corner in
the illustration). Its zvalue is the lowest zvalue in the search rectangle
A=Max is the corner of the search rectangle with maximum data (lower right corner in
the illustration). Its zvalue is the highest zvalue in the search rectangle
> x
   
  I****************** 
v  *  * 
y               *         
 *  * 
 *  * 
 *  * 
*  *
        *      *         
F *  *
 *'''''''''''******A 
  
  
                      
  
  
  
At first, division in xdirection.
The splitting line (double, vertical) splits the field into a low section
and an high section
The search range is overlapping the splitting line:
Min is in the low section (MinBit = 0).
Max is in the high section (MaxBit = 1).
F is in the low section, that means Fbit =0.
Situation FBit=0 MinBit =0 MaxBit =1
All zvalues in the high section are greater than the Z Values of F, Z(F).
They are no candidates for LITMAX whose Zvalue is smaller than Z(F) by
defintion.
So LITMAX must be in the low section. The search range is shrinked the way
that it is intirely in the low section. This is done by Jumping Max in
xdirection to the highest possible value in the lower section. This is
done by loading 011111... into max, starting fom the actial bit position
(first position here). See source code below: MAX := loadxxx011111( MAX)
(Compare Fig. 11 in the article, which is for interlaced data).
Result:
> x
splitting line
v
   
  I*********** 
v  * * 
y             *           
 * * 
 * * 
 * * 
* *
        *   *            
F * *
 *''''''''''A 
  
  
                      
  
  
  
Now, the low xsection becomes the active section.
Division in ydirection.
The splitting line (double, horizontal) splits the field into a
low section and a high section
This section is now Here, At first, division in xdirection.
> x
   
  I*********** 
v  * * 
y             *           
 * * 
 * * 
 * * 
* C
splitting line > = = = = = = = = = = = = = = = = = = =           
F * *
 *''''''''''A 
  
  
                      
  
  
  
Min is in the low section (MinBit = 0).
Max is in the high section (MaxBit = 1).
F is in the upper section, that means Fbit =1.
Situation: FBit=1 MinBit =0 MaxBit =1
All zvalues in the low section (here on top) are lower than the Z(F).
They are candidates for LITMAX whose Zvalue is smaller than Z(F) by
defintion. From them, the greatest possible Zvalue is in the corner
denoted by C, where the "candidate point" is placed. This is done by
C := loadxxx011111( MAX)
F can also be in the high section.
Search goes on in the high section, the range is now restricted to the
high section.
This is done by Jumping Min in ydirection to the lowest possible value
in the upper section.
This is done by loading 10000 into Min, starting fom the actial bit
position (secodn bit position here).
See source code below: MIN := loadxxx100000 ( MIN)
The search for F in the high section can fail (this is the case if the
staircase produced by F, see article, goes strictly along the splitting
line); if it fails, then the candidate point C is litmax.
result:
> x
   
   
v   
y                        
  
  
  
C
                       
F  I***********
 *''''''''''A 
  
  
                      
  
  
  
splitting line ^
All Zvalues in the search range are greater than Z(F). Search fails, the
candidate point C is LITMAX.
The remaining situations are evident, the complete LITMAX decision table is here:
FBit=0 MinBit =0 MaxBit =0 : everything in the same section: just continue
FBit=0 MinBit =0 MaxBit =1 : see above
FBit=0 MinBit =1 MaxBit =0 : impossible, because always F(Max) > F(Min)
FBit=0 MinBit =1 MaxBit =1 : F is in the lower section, the remaining search range
is completely in the upper section. LITMAX cannot be
in the remaining search range. Therfore LITMAX=C.
(In the source code below, the LITMAX variable is
directly used for C, there no special variable for C).
FBit=1 MinBit =0 MaxBit =0 : F is in the upper section, the remaining search range
completely in the lower section. LITMAX is the greates
possible ZValue in the lower section, this is MAX.
FBit=1 MinBit =0 MaxBit =1 : see above
FBit=1 MinBit =1 MaxBit =0 : impossible, because always F(Max) > F(Min)
FBit=1 MinBit =1 MaxBit =1 : everything in the same section: just continue
The search for BIGMIN is analogue to that, with symmetries. You find the table in the article.
Here is the Pascal source code:

(I tested it, but of course there is no guarantee that it is error free)
*)
{
//Basic definitions that can be modified to adopt to one's needs:
//no. of dimensions (here: 3), dimension names, Key format (here: word):
const ndims=3; //number of dimensions
type trecord = record
case integer of
1: (x,y,z: word;);
2: (valuefield: array[1..ndims] of word;
//so you can address the data either with a.x or a.valuefield[1]
maxvalue: word); //Max of x,y,z... just to speed up zvalueComparisons
end; //record
var startbit: integer;
}
procedure adjust_maxvalue(var R: trecord);
//adjusts R.maxvalue to the max value in R
var dim, max: integer;
begin
max:= R.valuefield[1];
for dim:=2 to ndims do
if R.valuefield[dim] > max then max := R.valuefield[dim];
R.maxvalue:=max;
end; //adjust_maxvalue
function LOAD_xxx10000(loadword: word; bitindex: integer): word;
//loads a bitpattern "10000.." into loadword, starting at bitindex
var pattern_00010000: word;
var pattern_00001111: word;
var pattern_11110000: word;
begin
pattern_00010000 := 1 shl (bitindex1);
pattern_00001111 := pattern_00010000 1;
pattern_11110000 := not pattern_00001111;
result:=loadword or pattern_00010000;
result:=result and pattern_11110000;
end; //LOAD_xxx10000
function LOAD_xxx01111(loadword: word; bitindex: integer): word;
//loads a bitpattern "011111.." into loadword, starting ad bitindex
var pattern_00010000: word;
var pattern_00001111: word;
var pattern_11101111: word;
begin
pattern_00010000 := 1 shl (bitindex1);
pattern_11101111 := not pattern_00010000;
pattern_00001111 := pattern_00010000 1;
result:=loadword and pattern_11101111;
result:=result or pattern_00001111;
end; //LOAD_xxx01111
procedure getLITMAX(MIN_, MAX_, DIV_: trecord; startbit: integer; var LITMAX: trecord);
//LITMAX is the maximum possible record in search Rectanle MIN_/Max_ whose ZValue
//is smaller than Z(DIV_)
//how it works: see above
//clearly some lines can be omitted ("case not possible.... error..."); I included it to keep
//the code in accordance with the description above
var dimindex: integer; //0.. ndims1
var bitindex: integer; //1..sizeof(word)
var bitmask: word;
var isminbit, ismaxbit, isdivbit: boolean;
begin
//precondition: Zvalue of DIV_ > Zvalue of MIN_
// Zvalue of DIV_ < Zvalue of MAX_
//to check it, use function zvalue_compare given below
LITMAX:=MAX_;
//for bitindex:=sizeof(word) * 8 (*one byte has 8 bits*) downto 1 do
for bitindex:=startbit downto 1 do
begin
bitmask:= 1 shl (bitindex1); //a bit faster when beginning with 1 shl (bitindex1) and then
//shr 1 each time. Easier to read this way.
for dimindex:=1 to ndims do
//Remark: interleaving order (3D) is xyzxyzxyz.... first dimension is x, second is y, etc.
//Therefore it does not work with "ndims downto 1".
// As depicted above and in the article, this visually yields an Ncurve, which is basically
// the same. Just interchange x and y in the figures to get a Z.
begin
isdivbit:=(DIV_.valuefield[dimindex] and bitmask) = bitmask;
isminbit:=(MIN_.valuefield[dimindex] and bitmask) = bitmask;
ismaxbit:=(MAX_.valuefield[dimindex] and bitmask) = bitmask;
//decision table:
if (not isdivbit) and (not isminbit) and(not ismaxbit) // 0 0 0 > no action, continue
then
else
if (not isdivbit) and (not isminbit) and(ismaxbit) // 0 0 1 > MAX:=LOAD xxx01111 into MAX
then
begin
MAX_.valuefield[dimindex] := LOAD_xxx01111(MAX_.valuefield[dimindex], bitindex);
end
else
if (not isdivbit) and (isminbit) and(not ismaxbit) // 0 1 0 not possible
then
error('not possible')
else
if (not isdivbit) and (isminbit) and(ismaxbit) // 0 1 1 finish
then
begin
exit;
end
else
if (isdivbit) and (not isminbit) and(not ismaxbit) // 1 0 0 LITMAX:=MAX; finish
then
begin
LITMAX:=MAX_; exit;
end
else
if (isdivbit) and (not isminbit) and(ismaxbit) // 1 0 1 LITMAX:=load xxx011111 into MAX
// MIN:=load xxx10000 into MIN
then
begin
LITMAX:=MAX_; LITMAX.valuefield[dimindex] := LOAD_xxx01111(MAX_.valuefield[dimindex],
bitindex);
//Remark: Don't forget "LITMAX:=MAX_"; it has become necessary, as interleaving is
//no more done explicitely. Here, the load function handles just one key, not the
//complete (bit interleaved ) record.
MIN_.valuefield [dimindex] := LOAD_xxx10000(MIN_.valuefield[dimindex], bitindex);
end
else
if (isdivbit) and (isminbit) and(not ismaxbit) // 1 1 0 not possible
then
error('not possible')
else
if (isdivbit) and (isminbit) and(ismaxbit) // 1 1 1 no action, continue
then (*no action*);
end; //for dimindex
end; //for bitindex
end; //getLITMAX
procedure getBIGMIN(MIN_, MAX_, DIV_: trecord; startbit: integer; var BIGMIN: trecord);
//BIGMIN is the minimum possible record in search rectanle MIN_/Max_ whose ZValue is
//bigger than Z(DIV_).
//how it works: see above
//clearly some lines can be omitted ("case not possible.... error..."); I included it to
//keep the code in accordance
// with the description above
var dimindex: integer; //0.. ndims1
var bitindex: integer; //1..sizeof(word)
var bitmask: word;
var isminbit, ismaxbit, isdivbit: boolean;
begin
//precondition: Zvalue of DIV_ > Zvalue of MIN_
// Zvalue of DIV_ < Zvalue of MAX_
//to check it, use function zvalue_compare given below
BIGMIN:=MIN_;
for bitindex:=startbit downto 1 do
//for bitindex:=sizeof(word) *8 (*one byte has 8 bits*) downto 1 do
begin
bitmask:= 1 shl (bitindex1); //a bit faster when beginning with 1 shl (bitindex1) and then
//shr 1 each time. Easier to read this way.
for dimindex:=1 to ndims do
//Remark: interleaving order (3D) is xyzxyzxyz.... first dimension is x, second is y, etc.
// it does not work with "ndims downto 1"
// As depicted above and in the article, this visually yields an Ncurve, which is basically
// the same. Just interchange x and y in the figures to get a Z.
begin
isdivbit:=(DIV_.valuefield[dimindex] and bitmask) = bitmask;
isminbit:=(MIN_.valuefield[dimindex] and bitmask) = bitmask;
ismaxbit:=(MAX_.valuefield[dimindex] and bitmask) = bitmask;
//decision table:
if (not isdivbit) and (not isminbit) and(not ismaxbit) // 0 0 0 > no action, continue
then
else
if (not isdivbit) and (not isminbit) and(ismaxbit) // 0 0 1 > BIGMIN:=LOAD xxx10000 into MIN
then // MAX :=LOAD xxx01111 into MAX
begin
BIGMIN:=MIN_; BIGMIN.valuefield[dimindex] := LOAD_xxx10000(MIN_.valuefield[dimindex],
bitindex);
//Remark: Don't forget "BIGMIN:=MIN"; it has become necessary, as interleaving is
//no more done explicitely. Here, the load function handles just one key, not the
//complete (bit interleaved ) record.
MAX_.valuefield [dimindex] := LOAD_xxx01111(MAX_.valuefield[dimindex], bitindex);
end
else
if (not isdivbit) and (isminbit) and(not ismaxbit) // 0 1 0 not possible
then
error ('not possible')
else
if (not isdivbit) and (isminbit) and(ismaxbit) // 0 1 1 bigmin:=min, finish
then
begin
BIGMIN:=MIN_; exit;
end
else
if (isdivbit) and (not isminbit) and(not ismaxbit) // 1 0 0 finish
then
exit
else
if (isdivbit) and (not isminbit) and(ismaxbit) // 1 0 1 load xxx10000 into MIN
then
begin
MIN_.valuefield[dimindex] := LOAD_xxx10000(MIN_.valuefield[dimindex], bitindex);
end
else
if (isdivbit) and (isminbit) and(not ismaxbit) // 1 1 0 not possible
then
error('not possible')
else
if (isdivbit) and (isminbit) and(ismaxbit) // 1 1 1 no action, continue
then (*no action*);
end; //for dimindex
end; //for bitindex
end; //getBIGMIN
//the function zvalue_compare described in the following becomes necessary because bit
//interleaving is not done explicitly. The data are left as they are and they are processed
//instead in interleaved manner.
//(Comparing interleaved data a and b was easier: just write "if a > b" and the like).
//type tcompare = (less, equal, greater); (defined elsewhere)
function zvalue_compare(a,b: trecord): tcompare;
//works with clobal ndims = number of dimensions
var dim: integer;
var a_xor_b: array [1..ndims] of word;
var AND12, XOR12: word;
(*
Task: go through the data columnwise, a and b in parallel, until you
hit an "1" Bit, this is the decisive bit position and the relevant dimension.
There:
if aBit is "1" and bBit is "0", then Z(a)>Z(b)
if aBit is "0" and bBit is "1", then Z(b)>Z(a)
if aBit is "1" and bBit is "1", then Z(a)>Z(b).
The latter rule applies because the interleaving scheme is xyzxyz...,
that means x has proirity against y, y has proirity against z, etc.
examples:
a b
x 00000 00100
y 01000 00000
z 00000 00000
Z(a)>Z(b) because the "1" is hit at a with "0" at b (y is the relevant dimension)
a b
x 00000 01000
y 01100 00000
z 01000 00000
Z(b)>Z(a) because the "1" is hit at b with "0" at a (x is the relevant dimension)
x 01000 01000
y 01000 01000
z 01000 00000
Z(a)>Z(b) because the "1" is hit at a with "1" at b (x is the relevant dimension)
So, a straightforward but slow realisation would be as follows (ndims=no. of dimensions):
var bitpos, dim: integer;
var bitmask: word;
var a_is_1, b_is_1: boolean;
begin
result:=equal; //init for the case that all values are Zero
for bitpos:=sizeof(word) *8 downto 1 do //one word has 2 bytes, one byte has 8 bits
begin
bitmask:= 1 shl (bitpos1);
for dim:=1 to ndims do
begin
a_is_1 :=(a.valuefield[dim] and bitmask) <>0;
b_is_1 :=(b.valuefield[dim] and bitmask) <>0;
if a_is_1 then if not b_is_1 then
begin
result:=greater; //Z(a)>Z(b)
exit;
end;
if b_is_1 then if not a_is_1 then
begin
result:=less; //Z(a) XOR12
if so, both y and x bits are 1 and the relevant dimension is x.
case 2: at the decisive bit position AND12 is 0
at the decisive bit position XOR12 is 1
here, clearly, XOR12 > AND12
if so, only one of the y and x bits has "1" and the relevant dimension
is the dimension which has "1".
In order to handle n dimensions, the thing is done first for dimensions 1 and 2, the relevant
dimension found thus far is kept. Then the same thing is done with the relevant dimension thus far
and dimension 3, wherehy again the (possibly new) relevant dimension is kept. And so on.
Here is the code:
*)
var relevant_dimension, run_dimension: integer;
begin
//to accelerate, this should be done already when creating the records...
//(it is done here just to be sure)
adjust_maxvalue((*var*) a);
adjust_maxvalue((*var*) b);
//try shortcut:
if (a.maxvalue XOR b.maxvalue) > (a.maxvalue AND b.maxvalue)
then
begin
//inc (shortcutzaehler);
if a.maxvalue > b.maxvalue
//if maxvalue_A > maxvalue_B
then
begin
result:=greater;
exit;
end
else
if a.maxvalue < b.maxvalue
//if maxvalue_A < maxvalue_B //das else kann man sich sparen
then
begin
result:=less;
exit;
end
end;
//end try shortcut
result:=equal;
//prepare a_xor_b for all dimensions:
for dim:=1 to ndims do
a_xor_b[dim]:=a.valuefield[dim] xor b.valuefield[dim];
relevant_dimension:=1; //init for the first dimension loop
run_dimension:=2; //init
//dimensionloop:
while run_dimension <= ndims do
begin
AND12 := a_xor_b[relevant_dimension] AND a_xor_b[run_dimension];
XOR12 := a_xor_b[relevant_dimension] XOR a_xor_b[run_dimension];
if AND12 > XOR12 then
//dimension 1 (x) is relevant as it is the first dimension in the interlacing scheme
begin
if a.valuefield[relevant_dimension] > b.valuefield[relevant_dimension]
then result:=greater
else result:=less;
//relevant_dimension for next loop remains
end
else
if AND12 < XOR12 then
begin
if a_xor_b[relevant_dimension] > a_xor_b[run_dimension]
then
begin //dim 1 is relevant
if a.valuefield[relevant_dimension] > b.valuefield[relevant_dimension]
then result:=greater
else result:=less;
//relevant_dimension next loop remains
end
else
begin //dim 2 is relevant
if a.valuefield[run_dimension] > b.valuefield[run_dimension]
then result:=greater
else result:=less;
relevant_dimension:=run_dimension;
(*new; this is the init decisive dimension for the next dimension loop*)
end;
end;
//dimension_is_ready:
inc(run_dimension);
end; //dimensionloop
//remark: result remains "equal" only if never a[i]<>b[i], that means only if a=b
end; //zvalue_compare
(*
Realisation hints:

Data types:
The above program works with word type data (no negatives). To handle integer data
(possibly negative), just invert the sign bit. This converts the data range to
0000..0 (lowest negative) to 1111..1 (highest positive).
The meaning behind the data is not relevant, only greaterlessequal matters! Anything
sortable that can be mapped to a 0000..0 to 1111..1 range, can be range searched (maybe
buckets of whatever that could again be searched anyhow).
To achieve this mapping for floating point data is surprisingly simple: If the float
is >=0, invert the sign bit, if the float is <0, invert all bits.
The consideration behind it is as follows: The sign bit is the first significant bit.
Unfortunately, it is 0 for positive data and 1 for negative data. So just invert it.
Exponent and mantissa are both not 2's complement (as it is the case with integers)
but in a range of 0000..0 (lowest value) to 1111..1 (highest value). Increasing the
exponent shifts more to the left, which corresponds to an increasing absolute value
of the float. The same is true for the mantissa: increasing the mantissa increases the
absolute value of the float. So, what we do is sorting the float according to its
absolute value. However, a negative with great absolute value is less than a negative
with small absolute value. This can simply be handled by inverting both exponent and
nmantissa (in addition to the sign bit).
Of course, such bit manipulations are to be reverted when reporting results.
Data structures:
The LITMAX/BIGMIN method can be helpful for any data structuring scheme. In the
original article cited above, binary search trees are used (that can be dynamically
balanced when insertions or deletions take place).
For data given as an array of records, sorted by "zvalue_compare", this original
method can be realized by binary search, without realizing a search tree explicitly.
For an array of records, it is somewhat better, to split the data hierarchically
into intervals and use for each interval processed both BIGMIN (of leftmost Record)
and LITMAX (of rightmost Record), narrowing the interval of ZValues remaining in
question even further.
The same can be done with other interval based data structures, e.g. skiplists
(that also can be dynamically balanced when insertions or deletions take place).
Performance considerations:
Align the data the way that for each dimension the leading bit is at the same position.
Example:
maximum key value first dimension is 00111010
maximum key value second dimension is 00000111
Here, shift the second dimension input data 3 bits to the left, shift the second
dimension output data 3 bits to the right.
Adjust the startbit to the maximum possible key value; this makes the LILMAX/BIGMIN
calculation a bit faster.
When using a data structure where in some step the number of consecutive records in
question is known (as it is the case for sorted arrays partitioned into intervals), make
a trivial sequential search for these records if this number of records is small enough,
(with limit at, say, 100..10000, depending on the computer and data situation).
Use the LITMAX / BIGMIN stuff only for the big leaps. The limit can be determinded
beforehand by experiments with typical test data and the computer at hand. It may be
even learned from experience while processing.
Double Records:
The thing works also when double records are possible (concerning searchable keys).
This is due to the fact that, when a record F encountered is within the search range,
LITMAX/BIGMIN are not calculated (as described in the article algorithm, page 75):
LITMAX/BIGMIN(F) results are records with zvalue less/greater than the Zcodes of F,
respectively. Therefore, don't use LITMAX/BIGMIN when F is in the search range and
double records are possible. When using it in this case anyway, scan the neighbours
of F in both directions as long they have the same keys as F.
*)