Hashing for Tolerant Index-Of
Roger Hui
Abstract
We consider the problem of x⍳y where x and y are real or complex, with tolerant comparison. We show how to use hashing to solve the problem.
Hashing algorithms for real and complex arguments have been implemented in Dyalog v13.0. APL models and test cases are provided. Benchmarks demonstrate the improvement over the previous implementation.
Acknowledgment
I am grateful to Jay Foad for several comments on the paper, in particular the elucidation of the region of tolerant equality.
0. Introduction
⎕io←0 is used throughout this text.
Two numbers x and y are tolerantly equal
relative to ⎕ct < 1 if
(|x-y) ≤ ⎕ct×(|x)⌈(|y)
(the ≤ in the definition is not tolerant)
[0].
The examples presented here are constructed
with ⎕ct←1e¯14 unless specified otherwise.
x⍳y is the index-of function, the smallest index in x of each element of y , modelled as follows:
iota←{(0⌷⍴⍺)-+/∨\⍵∘.=⍺}
3 1 4 1 5 9 iota 2 3⍴⍳6
6 1 6
0 2 4
3 1 4 1 5 9 iota 1+1e¯15 1e¯13
1 6
x⍳y is tolerant because the comparison used (the = in iota) is tolerant.
x⍳y is commonly not implemented as in iota because it runs in time of order (⍴x)×(⍴y) . SHARP APL and Dyalog APL v12.1 use a method which sorts x , then for each element of y do a binary search to find a tolerant match, followed by a linear search on all the tolerant matches [1] [2]. This has time of order (⍟⍴x)×(⍴x)+(⍴y) on the average but (⍴x)×(⍴y) in the worst case.
3 rand 0.5
430.609375 ¯191.5625 ¯30.82421875
x←8e4 rand 0.5
5 timer 'x⍳x' ⍝ timing in Dyalog v12.1
0.0908
monster←{1+(0.0001×⎕ct)×?⍵⍴100000}
x←monster 8e4
5 timer 'x⍳x' ⍝ timing in Dyalog v12.1
54.0406
| n | Typical | Monster | ||
| 1e4 | 0.0186 | 0.6718 | ||
| 2e4 | 0.0250 | 2.5282 | ||
| 4e4 | 0.0406 | 10.6000 | ||
| 8e4 | 0.0908 | 54.0406 |
(rand , timer , and other functions used in the text are defined in the Appendix.)
Hashing algorithms also run in time of order (⍴x)×(⍴y)
in the worst case
but order (⍴x)+(⍴y) on the average.
We show how hashing can be used for tolerant x⍳y
on real and complex arguments.
1. Real Arguments
The definition of tolerant equality implies that the closed interval of numbers tolerantly equal to a real number t has endpoints t×1-⎕ct and t÷1-⎕ct .
As the diagram (exaggeratingly) illustrates, the interval is not symmetric around t . The diagram is reflected about the origin if t is negative.
Koss[3] made an observation in 1985 that enabled hashing for tolerant index-of: The sign, exponent, and leading digits of one or both of t×1-⎕ct and t÷1-⎕ct are the same as those of any number in the interval, and the number of such digits is a function of ⎕ct . The following examples illustrate the point in base 10:
intv←{⍺←⎕ct ⋄ ⍵×(1-⍺),1,÷1-⍺}
⍪ 0.0001 intv ○1
3.141278494
3.141592654
3.141906844
The digits of an endpoint (but not both endpoints) can be completely different from those of t if t is near a power of the base.
⍪ 0.0001 intv 1e20 9.99900000E19 1.00000000E20 1.00010001E20 ⍪ 0.0001 intv ¯1e¯20+1e¯27 ¯9.99899900E¯21 ¯9.99999900E¯21 ¯1.00009991E¯20
An IEEE-754 64-bit floating point number has a sign bit, followed by an exponent, followed by a mantissa. (IEEE-754 is used by most current machines.) For such numbers the observation simplifies to saying that the leading bits of one or both of t×1-⎕ct and t÷1-⎕ct are the same as those of any number in the interval.
hexf ⍪ intv ○1 ⍝ hexf is from the dfns workspace 400921FB54442CD1 400921FB54442D18 400921FB54442D5F
Again, if t is near a power of the base, the bits of t×1-⎕ct or t÷1-⎕ct can be substantially different from those of t .
hexf ⍪ intv 4 400FFFFFFFFFFFA6 4010000000000000 401000000000002D hexf ⍪ intv 4 - 1.2e¯14 400FFFFFFFFFFF8B 400FFFFFFFFFFFE5 401000000000001F
A hashing algorithm for x⍳y for real x and y :
| 0. | Compute a mask based on ⎕ct (function hashmask). | |
| 1. | Initialize a hash table h of size prime q larger than 2×⍴x . | |
| 2. | For each element i⌷x , mask out the trailing bits, then compute a hash index, then insert into the hash table at that index. Collisions of elements that are not exactly equal are resolved by inserting into the next hash table entry, wrapping back to the front as necessary. Exactly equal elements are inserted just once. | |
| 3. | For each element i⌷y , for each of (i⌷y)×1-⎕ct and (i⌷y)÷1-⎕ct , mask out the trailing bits, then compute a hash index, then look in the hash table starting at that index, using tolerant equality. |
The hash index for a 64-bit masked floating point
number t is
Other than the parts involving ⎕ct and the masking of the trailing bits, this is Algorithm L of Knuth [4] with a load factor of 0.5. The algorithm is modelled in APL as diota , written as a “tradfn” to facilitate translation into C.
x←1+(0.25×⎕ct)×?200⍴851 y←1+(0.25×⎕ct)×?300⍴951 (x iota y) ≡ x diota y 1
2. Complex Arguments
The region of numbers tolerantly equal to a complex number t is a near-circle with radius r←⎕ct×|t .
More precisely, it is the union of two circular regions: The smaller is centered on t with radius r . The larger derives as follows. A point u on the boundary with magnitude ≥ |t satisfies (|u-t)=⎕ct×|u . That is, the ratio between |u-t and |u is ⎕ct , a constant; therefore, that boundary is a circle of Apollonius with foci t and the origin. The larger circle circumscribes the two points on the smaller circle with magnitude |t and the point with magnitude (|t)÷1-⎕ct collinear with t and the origin. (A more precise diagram than the one above can be found in Appendix A0.)
The union of two circles (or even a circle) is inconvenient to deal with computationally. For present purposes it is good enough to work with a square that contains the region, “good enough” in the sense that even though numbers outside the region will be hashed to the same index, such numbers will be rejected when they are tested for tolerant equality. We choose as the delta (half the length of a side of the square) ⎕ct × √2 × the larger of the magnitudes of the real and imaginary parts of t . This choice avoids the problem and expense of computing the magnitude of t , requisite for a tighter fitting square.
delta←{⍺←⎕ct ⋄ ⍺×(2*0.5)×⌈/|9 11○⍵}
delta 3j4
5.656854249E¯14
delta 3e6j4
4.242640687E¯8
delta 3j4e6
5.656854249E¯8
Like the real case, the leading bits of the four corners are representative of those of points in the square. However, unlike the real case, the number of leading bits will depend on ⎕ct and t and (in general) are different for the real and imaginary parts.
hexf 3 4 4008000000000000 4010000000000000 hexf 3 4+○1e¯14 4008000000000047 4010000000000023 hexf 3 1e¯20 4008000000000000 3BC79CA10C924223 hexf 3, ○¯1e¯21 4008000000000000 BBADABE7AA9FDCC0
The first example shows that if the limbs (the real and imaginary parts) are similar in length, then the numbers in the region can have many common leading bits in both limbs. The second example shows that if one limb is much longer than the other, then two tolerantly equal numbers can have completely different bits in their shorter limbs.
For a complex number t , the number of leading bits for the longer limb is computed as before, a function of ⎕ct alone (function hashmask). The number of leading bits for the shorter limb s with a delta of d is the length of the common prefix for s and s+d or for s and s-d , whichever is longer (function hashmask2).
A hashing algorithm for x⍳y for complex x and y :
| 0. | Compute a mask based on ⎕ct (function hashmask). | |||||
| 1. | Initialize a hash table h of size prime q larger than 2×⍴x . | |||||
| 2. | For each element i⌷x , compute d←delta i⌷x and m2←s hashmask2 d where s is the length of the shorter limb, then mask out the trailing bits separately for the real and imaginary parts, then compute a hash index, then insert into the hash table at that index. Collisions of elements that are not exactly equal are resolved by inserting into the next hash table entry, wrapping back to the front as necessary. Exactly equal elements are inserted just once. | |||||
| 3. | For each element i⌷y ,
compute d←delta i⌷y
and m2←s hashmask2 d where s
is the length of the shorter limb,
then for each of the “four corners”
|
The hash index for a 128-bit masked complex
number t is
The algorithm is modelled in APL as ziota , written as a “tradfn” to facilitate translation into C.
x←*0j1×1+(0.25×⎕ct)×?199⍴223 y←*0j1×1+(0.25×⎕ct)×?307⍴239 (x iota y)≡x ziota y 1
3. Implementation
The algorithms modelled by diota and ziota
have been implemented in C in Dyalog v. 13.0.
The same algorithms underlie the implementation
of x∊y x~y x∩y x∪y ∪x .
The implementation has additional optimizations
for ⎕ct←0 and
for x⍳x and ∪x
where the initial hashing of x (step 2)
is merged into the main loop (step 3).
The hash table is retained for reuse in f←x∘⍳
and in x∘⍳¨ (also ∊∘y , etc.).
4. Tests
assert p is a monad which does nothing if ~0∊p and signals assertion failure if 0∊p . For example:
x←1000 rand 0.5 assert (x⍳x)≡x iota x assert 3 = 3×1+2×⎕ct assertion failure assert 3=3×1+2×⎕ct ∧
The test suites use assert extensively. A test function returns a 1 if it runs successfully, or suspends with an assertion failure at the first point of failure. Four main functions are included (see Appendix A2):
test_iotaD ⍝ x⍳y for real x and y 1 test_iotaD_retained ⍝ f←x∘⍳ and x∘⍳¨ for real x 1 test_iotaZ ⍝ x⍳y for complex x and y 1 test_iotaZ_retained ⍝ f←x∘⍳ and x∘⍳¨ for complex x 1
5. Timings
The following table compares the timings in Dyalog v13.0 vs. v12.1.
{x y←⍵ ⍵ rand¨0.5 ⋄ 5 timer'x⍳y'}¨1e6×1 2 4 8
{x ←⍵ rand 0.5 ⋄ 5 timer'x⍳x'}¨1e6×1 2 4 8
| Real x⍳y | Real x⍳x | ||||||||||||
| n | v13.0 | v12.1 | Ratio | v13.0 | v12.1 | Ratio | |||||||
| 1e6 | 0.5812 | 2.3438 | 4.03 | 0.3282 | 2.1282 | 6.48 | |||||||
| 2e6 | 1.2438 | 5.2124 | 4.19 | 0.6906 | 5.0874 | 7.37 | |||||||
| 4e6 | 2.2938 | 11.1250 | 4.85 | 1.3438 | 11.1844 | 8.32 | |||||||
| 8e6 | 4.9750 | 25.9906 | 5.22 | 2.9562 | 24.8844 | 8.42 | |||||||
The following table presents timings done only in Dyalog v13.0 as they are impossible or impractical in v12.1.
{x y←⍵ ⍵ rand¨0j1 ⋄ 5 timer'x⍳y'}¨1e6×1 2 4 8
{x ←⍵ rand 0j1 ⋄ 5 timer'x⍳x'}¨1e6×1 2 4 8
{x y←monster¨⍵ ⍵ ⋄ 5 timer'x⍳y'}¨1e6×1 2 4 8
{x ←monster ⍵ ⋄ 5 timer'x⍳x'}¨1e6×1 2 4 8
| Complex | Monster | |||||||
| n | x⍳y | x⍳x | x⍳y | x⍳x | ||||
| 1e6 | 1.6500 | 0.9532 | 1.8032 | 1.6750 | ||||
| 2e6 | 3.2750 | 2.1436 | 3.5032 | 3.2314 | ||||
| 4e6 | 5.1468 | 3.6906 | 6.8906 | 6.8562 | ||||
| 8e6 | 9.7032 | 7.2250 | 13.3812 | 13.3436 | ||||
References
| 0. | ISO/IEC, Programming Language Extended APL, ISO/IEC 13751:2001(E), 2001-02-01; Section 5.2.5, Tolerantly-Equal, page 19. | |
| 1. | Bernecky, Robert, e-mail, 2010-05-11 05:04. | |
| 2. | Streeter, Geoff, e-mail, 2010-05-13 00:25. | |
| 3. | Koss, Steven L., Tolerant Hash Functions, Minnowbrook, 1985-10. | |
| 4. | Knuth, Donald E., The Art of Computer Programming, Volume 3, Searching and Sorting, Addison-Wesley, 1973; Section 6.4. |
Appendices
A0. Region of Tolerant Equality
The region of numbers tolerantly equal to a complex number t is a near-circle with radius r←⎕ct×|t .
More precisely, it is the union of two circular regions: The smaller (in black) is centered on t with radius r . The larger derives as follows. A point u on the boundary with magnitude ≥ |t satisfies (|u-t)=⎕ct×|u . That is, the ratio between |u-t and |u is ⎕ct , a constant; therefore, that boundary is a circle of Apollonius with foci t and the origin. The larger circle (in red) circumscribes the two points p and q on the smaller circle with magnitude |t (in green) and the point s with magnitude (|t)÷1-⎕ct collinear with t and the origin.
 |
r←⎕ct×|t (|u-t) = ⎕ct×|u (|p) = |t (|q) = |t (|s) = (|t)÷1-⎕ct |
If t is real, the region simplifies to the closed interval with endpoints t×1-⎕ct and t÷1-⎕ct .