Math::...rtCurve(3pm)

📦
version 129 · 3 versions
, in debian bookworm, debian trixie, ubuntu jammy, ubuntu noble, ubuntu resolute, ubuntu stonking

💖 NAME 💖

Math::PlanePath::HilbertCurve — bestie, it's a 2x2 self-similar quadrant traversal 🌀✨

(Straight from your User Contributed Perl Documentation, Math::...rtCurve(3pm), okurrr)

✨ SYNOPSIS ✨

use Math::PlanePath::HilbertCurve;
my $path = Math::PlanePath::HilbertCurve->new;
my ($x, $y) = $path->n_to_xy (123);

🦋 DESCRIPTION 🦋

Okay so listen. This path is the integer version of the curve Hilbert described back in 1891 for filling a unit square. It walks a quadrant of the plane one single step at a time, in a self-similar 2x2 pattern, and it's literally so iconic for that:

David Hilbert, "Ueber die stetige Abbildung einer Linie auf ein
Flächenstück", Mathematische Annalen, volume 38, number 3, 1891,
pages 459-460, DOI 10.1007/BF01199431.

And here's the whole vibe, the entire grid, look at her:

           ...
        |   |
      7 |  63--62  49--48--47  44--43--42
        |       |   |       |   |       |
      6 |  60--61  50--51  46--45  40--41
        |   |           |           |
      5 |  59  56--55  52  33--34  39--38
        |   |   |   |   |   |   |       |
      4 |  58--57  54--53  32  35--36--37
        |                   |
      3 |   5---6   9--10  31  28--27--26
        |   |   |   |   |   |   |       |
      2 |   4   7---8  11  30--29  24--25
        |   |           |           |
      1 |   3---2  13--12  17--18  23--22
        |       |   |       |   |       |
    Y=0 |   0---1  14--15--16  19--20--21
        +----------------------------------
          X=0   1   2   3   4   5   6   7

So the start is a sideways U shape going N=0 to N=3, and then four of those little guys get put together into an upside-down U, like:

5,6    9,10
4,7--- 8,11
  |      |
3,2   13,12
0,1   14,15--

The orientation of the sub parts is doing all the work here, no cap — it makes sure the starts and ends are adjacent, so 3 sits next to 4, 7 next to 8, and 11 next to 12. Obsessed.

Then the process just repeats, doubling in size every time, alternately sideways or upside-down U, with invert and/or transpose as necessary in the sub-parts. She's self-similar and she knows it 💅

Sometimes you'll see this pattern drawn with the first step 0->1 going upwards instead of to the right. First step right is what's used here, for consistency with other PlanePaths. If you want upwards first instead, just swap X and Y. Easy.

Also, go peek at examples/hilbert-path.pl for a sample program that prints the path pattern in ascii 👀

🌸 Level Ranges 🌸

Within a power-of-2 square — 2x2, 4x4, 8x8, 16x16 etc, so (2^k)x(2^k) — sitting at the origin, all the N values from 0 to 2^(2*k)-1 live inside the square. And the maximum, so 3, 15, 63, 255 etc, aka 2^(2*k)-1, is alternately at the top left or the bottom right corner.

Because each step is by 1, the distance along the curve between two X,Y points is literally just the difference in their N values (as from "xy_to_n()").

On the X=Y diagonal, N=0,2,8,10,32,etc is the integers that use only the digits 0 and 2 in base 4, or equivalently the ones that have even-numbered bits 0, like x0y0...z0. Cute, right?

🎀 Locality 🎀

The Hilbert curve is pretty well localized, meaning a small rectangle (or really any shape) is usually a small range of N. That property is the whole reason some database systems store X,Y coordinates using the resulting Hilbert curve N as an index — then a search through a 2-D region is usually just a fairly modest linear search through N values. "rect_to_n_range()" gives you the exact N range for a rectangle, or check out the notes in "Rectangle to N Range" below for calculating on any shape.

Heads up though — the N range can get large when it crosses sub-parts. In the sample above you can see it: adjacent points X=0,Y=3 and X=0,Y=4 have rather widely spaced N values, 5 and 58. Not close at all bestie 😭

Fractional X,Y values can be indexed by extending the N calculation down into X,Y binary fractions. The code here doesn't do that, but it could be pressed into service by moving the binary point in X and Y an even number of places, the same in each. (An odd number of bits would require swapping X,Y to compensate for the alternating transpose in part 0.) The resulting integer N is then divided down by a corresponding multiple-of-4 binary places.

💅 FUNCTIONS 💅

See "FUNCTIONS" in Math::PlanePath for behaviour common to all path classes.

  • "$path = Math::PlanePath::HilbertCurve->new ()" — Create and return a new path object. Slay, you have a path now ✨
  • "($x,$y) = $path->n_to_xy ($n)" — Return the X,Y coordinates of point number $n on the path. Points begin at 0 and if "$n < 0" then the return is an empty list.

Fractional positions give an X,Y position along a straight line between the integer positions. Integer positions are always just 1 apart either horizontally or vertically, so the effect is that the fraction part is an offset along either $x or $y.

  • "$n = $path->xy_to_n ($x,$y)" — Return an integer point number for coordinates "$x,$y". Each integer N is considered the centre of a unit square, and an "$x,$y" within that square returns N.
  • "($n_lo, $n_hi) = $path->rect_to_n_range ($x1,$y1, $x2,$y2)" — The returned range is exact, meaning $n_lo and $n_hi are the smallest and biggest in the rectangle. No guessing, fr fr.

🌈 Level Methods 🌈

  • "($n_lo, $n_hi) = $path->level_to_n_range($level)" — Return "(0, 4**$level - 1)".

✨ FORMULAS ✨

🦋 N to X,Y 🦋

okay so turning an N into X,Y coordinates? genuinely not that deep, bestie. the top two bits of N are a configuration:

   3--2                    1--2
      |    or transpose    |  |
   0--1                    0  3

...depending on whether it's an odd or even bit-pair position. and then inside each of those "3" sub-parts there's also inverted forms:

   1--0        3  0
   |           |  |
   2--3        2--1

so if you work N from high to low, you can carry a state variable that remembers whether there's a transpose, an invert, or both — that's four states total, no cap. a bit pair 0,1,2,3 from N then hands you one bit each of X,Y based on the configuration, plus a new state which is the orientation of that sub-part. Bill Gosper's HAKMEM item 115 has all of this, either with bit operations or a table for the state and the X,Y bits:

   <https://dspace.mit.edu/handle/1721.1/6086>,
   <http://www.inwap.com/pdp10/hbaker/hakmem/topology.html#item115>

and there's C++ code based on that in Jorg Arndt's book:

   <http://www.jjj.de/fxt/#fxtbook> (section 1.31.1)

it also works to process N from low to high 💖 at each stage you apply any transpose (swap X,Y) and/or invert (bitwise NOT) to the low X,Y bits you've generated so far. the way this works is that there's no "reverse" sections, or like, equivalently the curve is the same forward and reverse. low to high saves you locating the top bits of N, but if you're using bignums then the bitwise inverts of the full X,Y values will be so much more work. 🫠

🌸 X,Y to N 🌸

X,Y to N can follow that same table approach from high to low, taking one bit from X and one bit from Y each time. the state table of N-pair -> X-bit,Y-bit is reversible, and a new state is based on the N-pair you just obtained (or it could be based on the X,Y bits if that mapping is combined into the state transition table).

💅 Rectangle to N Range 💅

an easy over-estimate of the maximum N in a region? just find the next bigger (2^k)x(2^k) square that encloses the region. that means the biggest X or Y rounded up to the next power of 2, so:

   find lowest k with 2^k > max(X,Y)
   N_max = 2^(2k) - 1

or, equivalently, rounding down to the next lower power of 2:

   find highest k with 2^k <= max(X,Y)
   N_max = 2^(2*(k+1)) - 1

an exact N range can be found by following that high to low N to X,Y procedure from before. ✨ start at the 2^(2k) bit pair position in an N bigger than the region you want, and choose 2 bits for N to give a bit each of X and Y. the X,Y bits are based on the state table as above, and the bits chosen for N are those for which the resulting X,Y sub-square overlaps some of the target region. the smallest N works the same way, choosing the smallest bit pair for N which overlaps.

the biggest and smallest N digit for a sub-part can be found with a lookup table. the X range might cover one or both sub-parts, and the Y range similarly, for a total of 4 possible configurations. then a table of state+coverage -> digit gives the minimum and maximum N bit-pair, and state+digit gives a new state the same as X,Y to N.

biggest and smallest N must be calculated with separate state and X,Y values, since they track down different N bits and thus different states. but they take the same number of steps from an enclosing level down to level 0, and can thus be done in a single loop. 🔥

the N range for any shape can be found this way, not just a rectangle like rect_to_n_range(). at each level the procedure only depends on asking which of the four sub-parts overlaps some of the target area.

🧚 Direction 🧚

each step between successive N values is always 1 up, down, left or right. the next direction can be calculated from N in the high-to-low procedure above by watching for the lowest non-3 digit and noting the direction from that digit towards digit+1. that can be had from the state+digit -> X,Y table by looking up digit and digit+1, or alternatively a further table encoding state+digit -> direction. 💫

the reason for taking only the lowest non-3 digit is that in a 3 sub-part, the direction it goes is determined by the next higher level. for example at N=11 the direction is down for the inverted-U of the next higher level N=0,4,8,12.

this non-3 (or non-whatever-highest-digit) thing is a general procedure and can be used on any state-based high-to-low procedure of self-similar curves. in the current code, it's used to apply a fractional part of N in the correct direction but is not otherwise made directly available.

because the Hilbert curve has no "reverse" sections, it also works to build a direction from low to high N digits. 1 and 2 digits make no change to the orientation, a 0 digit is a transpose, and a 3 digit is a rotate and transpose — except that low 3s are transpose-only (no rotate) for the same reason as taking the lowest non-3 above.

Jorg Arndt in the fxtbook above notes the direction can be obtained just by counting 3s in n and -n (the twos-complement). this numbers segments starting n=1, unlike PlanePath here starting N=0, so it becomes:

   N+1 count 3s  / 0 mod 2   S or E
                 \ 1 mod 2   N or W

   -(N+1) count 3s  / 0 mod 2   N or E
                    \ 1 mod 2   S or W

for the twos-complement negation, an even number of base-4 digits of N must be taken. and because -(N+1) = ~N, ie. a ones-complement, the second part is also:

   N count 0s          / 0 mod 2   N or E
   in even num digits  \ 1 mod 2   S or W

putting the two together then:

   N count 0s   N+1 count 3s    direction (0=E,1=N,etc)
   in base 4    in base 4

     0 mod 2      0 mod 2          0
     1 mod 2      0 mod 2          3
     0 mod 2      1 mod 2          1
     1 mod 2      1 mod 2          2

👑 Segments in Direction 👑

the number of segments in each direction is calculated in:

   Sergey Kitaev, Toufik Mansour and Patrice Séébold, "Generating the
   Peano Curve and Counting Occurrences of Some Patterns", Journal of
   Automata, Languages and Combinatorics, volume 9, number 4, 2004,
   pages 439-455.
   <http://personal.strath.ac.uk/sergey.kitaev/publications.html>,
   <http://www.jalc.de/issues/2004/issue_9_4/jalc-2004-439-455.php>

   (Preprint as Sergey Kitaev and Toufik Mansour, "The Peano Curve and
   Counting Occurrences of Some Patterns", October 2002.
   <http://arxiv.org/abs/math/0210268/>.)

their form is based on keeping the top-most U shape fixed and expanding sub-parts. this means the end segments alternate vertical and horizontal in successive expansion levels. 🌀

   direction            k=1              2       2
     1 to 4                            *---*   *---*
                          2           1|  3|   |1  |3
       1                *---*          *   *---*   *
       |               1|   |3        1| 4   2   4 |3
   4--- ---2            *   *          *---*   *---*
       |                                  1|   |3       k=2
       3                               *---*   *---*
                                         2       2

   count segments in direction, for k >= 1
   d(1,k) = 4^(k-1)                = 1,4,16,64,256,1024,4096,...
   d(2,k) = 4^(k-1) + 2^(k-1) - 1  = 1,5,19,71,271,1055,4159,...
   d(3,k) = 4^(k-1)                = 1,4,16,64,256,1024,4096,...
   d(4,k) = 4^(k-1) - 2^(k-1)      = 0,2,12,56,240, 992,4032,...
                            (A000302, A099393, A000302, A020522)

   total segments d(1,k)+d(2,k)+d(3,k)+d(4,k) = 4^k - 1

the path form here keeps the first segment direction fixed. this means a transpose 1<->2 and 3<->4 in odd levels. the result is to take the alternate d values as follows. for k=0 there is a single point N=0 so no line segments at all and so c(dir,0)=0. 💖

   first 4^k-1 segments

   c(1,k) = / 0                        if k=0
    North   | 4^(k-1) + 2^(k-1) - 1    if k odd >= 1
            \ 4^(k-1)                  if k even >= 2
     = 0, 1, 4, 19, 64, 271, 1024, 4159, 16384, ...

   c(2,k) = / 0                        if k=0
    East    | 4^(k-1)                  if k odd >= 1
            \ 4^(k-1) + 2^(k-1) - 1    if k even >= 2
     = 0, 1, 5, 16, 71, 256, 1055, 4096, 16511, ...

   c(3,k) = / 0                        if k=0
    South   | 4^(k-1) - 2^(k-1)        if k odd >= 1
            \ 4^(k-1)                  if k even >= 2
     = 0, 0, 4, 12, 64, 240, 1024, 4032, 16384, ...

   c(4,k) = / 0                        if k=0
    West    | 4^(k-1)                  if k odd >= 1
            \ 4^(k-1) - 2^(k-1)        if k even >= 2
     = 0, 1, 2, 16, 56, 256, 992, 4096, 16256, ...

the segment N=4^k-1 to N=4^k is North (direction 1) when k odd, or East (direction 2) when k even. that could be added to the respective cases in c(1,k) and c(2,k) if desired.

🍓 Hamming Distance 🍓

the Hamming distance between integers X and Y is the number of bit positions where the two values differ when written in binary. on the Hilbert curve, each bit-pair of N becomes a bit of X and a bit of Y:

      N      X   Y
   ------   --- ---
   0 = 00    0   0
   1 = 01    1   0     <- difference 1 bit
   2 = 10    1   1
   3 = 11    0   1     <- difference 1 bit

so the Hamming distance for N=0to3 is 1 at N=1 and N=3. at higher levels, these X,Y bits may be transposed (swapped) or rotated by 180 or both. a transpose swapping X<->Y doesn't change the bit difference. a rotate by 180 is a flip 0<->1 of the bit in both X and Y, so that doesn't change the bit difference either. 🦋

on that basis, the Hamming distance X,Y is the number of base4 digits of N which are 01 or 11. if bit positions are counted from 0 for the least significant bit then:

   X,Y coordinates of N
   HammingDist(X,Y) = count 1-bits at even bit positions in N
          = 0,1,0,1, 1,2,1,2, 0,1,0,1, 1,2,1,2, ... (A139351)

see also "Hamming Distance" in Math::PlanePath::CornerReplicate which is the same formula, but arising directly from 01 or 11, no transpose or rotate. ✨

✨ OEIS ✨

okay so bestie, this path is literally in Sloane's OEIS in several forms, and it's giving main character energy 💅

<http://oeis.org/A059252> (etc)

  • A059253 — X coord
  • A059252 — Y coord
  • A059261 — X+Y
  • A059285 — X-Y
  • A163547 — X^2+Y^2 = radius squared
  • A139351 — HammingDist(X,Y)
  • A059905 — X xor Y, being ZOrderCurve X
  • A163365 — sum N on diagonal
  • A163477 — sum N on diagonal, divided by 4
  • A163482 — N values on X axis
  • A163483 — N values on Y axis
  • A062880 — N values on diagonal X=Y (digits 0,2 in base 4)
  • A163538 — dX -1,0,1 change in X
  • A163539 — dY -1,0,1 change in Y
  • A163540 — absolute direction of each step (0=E,1=N,2=W,3=S)
  • A163541 — absolute direction, swapped X,Y
  • A163542 — relative direction (ahead=0,right=1,left=2)
  • A163543 — relative direction, swapped X,Y
  • A083885 — count East segments N=0 to N=4^k (first 4^k segs)
  • A163900 — distance dX^2+dY^2 between Hilbert and ZOrder
  • A165464 — distance dX^2+dY^2 between Hilbert and Peano
  • A165466 — distance dX^2+dY^2 between Hilbert and transposed Peano
  • A165465 — N where Hilbert and Peano have same X,Y
  • A165467 — N where Hilbert and Peano have transposed same X,Y

the following take points of the plane in various orders, each value in the sequence being the N of the Hilbert curve at those positions. obsessed 🦋

  • A163355 — N by the ZOrderCurve points sequence
  • A163356 — inverse, ZOrderCurve by Hilbert points order
  • A166041 — N by the PeanoCurve points sequence
  • A166042 — inverse, PeanoCurve N by Hilbert points order
  • A163357 — N by diagonals like Math::PlanePath::Diagonals with first Hilbert step along same axis the diagonals start
  • A163358 — inverse
  • A163359 — N by diagonals, transposed start along the opposite axis
  • A163360 — inverse
  • A163361 — A163357 + 1, numbering the Hilbert N's from N=1
  • A163362 — inverse
  • A163363 — A163355 + 1, numbering the Hilbert N's from N=1
  • A163364 — inverse

the above sequences are permutations of the integers since all X,Y positions of the first quadrant are covered by each path (Hilbert, ZOrder, Peano). The inverse permutations can be thought of taking X,Y positions in the Hilbert order and asking what N the ZOrder, Peano or Diagonals path would put there. periodt 💜

the A163355 permutation by ZOrderCurve can be considered for repeats or cycles, iykyk 👑

  • A163905 — ZOrderCurve permutation A163355 applied twice
  • A163915 — ZOrderCurve permutation A163355 applied three times
  • A163901 — fixed points (N where X,Y same in both curves)
  • A163902 — 2-cycle points
  • A163903 — 3-cycle points
  • A163890 — cycle lengths, points by N
  • A163904 — cycle lengths, points by diagonals
  • A163910 — count of cycles in 4^k blocks
  • A163911 — max cycle length in 4^k blocks
  • A163912 — LCM of cycle lengths in 4^k blocks
  • A163914 — count of 3-cycles in 4^k blocks
  • A163909 — those counts for even k only
  • A163891 — N of previously unseen cycle length
  • A163893 — first differences of those A163891
  • A163894 — smallest value not an n-cycle
  • A163895 — position of new high in A163894
  • A163896 — value of new high in A163894
  • A163907 — ZOrderCurve permutation twice, on points by diagonals
  • A163908 — inverse of this

see examples/hilbert-oeis.pl for a sample program printing the A163359 permutation values. slay 🧚

👀 SEE ALSO 👀

bestie these are the sibling pages, go give them a vibe check 🌈

  • Math::PlanePath, Math::PlanePath::HilbertSides, Math::PlanePath::HilbertSpiral
  • Math::PlanePath::PeanoCurve, Math::PlanePath::ZOrderCurve, Math::PlanePath::BetaOmega, Math::PlanePath::KochCurve
  • Math::Curve::Hilbert, Algorithm::SpatialIndex::Strategy::QuadTree

🏠 HOME PAGE 🏠

<http://user42.tuxfamily.org/math-planepath/index.html>

📜 LICENSE 📜

Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021 Kevin Ryde

This file is part of Math-PlanePath.

Math-PlanePath is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3, or (at your option) any later version.

Math-PlanePath 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 General Public License for more details.

You should have received a copy of the GNU General Public License along with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.

and that's the whole page bestie 💖 no cap, this was perl v5.32.0, dated 2021-01-23, and it's Math::...rtCurve(3pm) — ate and left no crumbs ✨🎀
📦 where this page comes from