Provided by: libmath-planepath-perl_122-1_all bug

NAME

       Math::PlanePath::Flowsnake -- self-similar path through hexagons

SYNOPSIS

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

DESCRIPTION

       This path is an integer version of the flowsnake curve by William Gosper.

       A single arm of the curve fills 1/3 of the plane, spiralling around anti-clockwise ever
       fatter and with jagged self-similar edges.

                                39----40----41                        8
                                  \           \
                 32----33----34    38----37    42                     7
                   \           \        /     /
                    31----30    35----36    43    47----48            6
                         /                    \     \     \
                 28----29    17----16----15    44    46    49--..     5
                /              \           \     \  /
              27    23----22    18----19    14    45                  4
                \     \     \        /     /
                 26    24    21----20    13    11----10               3
                   \  /                    \  /     /
                    25     4---- 5---- 6    12     9                  2
                            \           \         /
                              3---- 2     7---- 8                     1
                                  /
                           0---- 1                                  Y=0

            X=-4 -3 -2 -1  0  1  2  3  4  5  6  7  8  9 10 11

       The points are spread out on every second X coordinate to make little triangles with
       integer coordinates, per "Triangular Lattice" in Math::PlanePath.

       The base pattern is the seven points 0 to 6,

           4---- 5---- 6
            \           \
              3---- 2
                  /
           0---- 1

       This repeats at 7-fold increasing scale, with sub-sections rotated according to the edge
       direction and the 1, 2 and 6 sections in reverse.  The next level can be seen at the
       multiple-of-7 points N=0,7,14,21,28,35,42,49.

                                         42
                             -----------    ---
                          35                   ---
              -----------                         ---
           28                                        49 ---
             ---
                ----                  14
                    ---   -----------  |
                       21              |
                                      |
                                     |
                                     |
                               ---- 7
                          -----
                     0 ---

       Notice this is the same shape as N=0 to N=6, but rotated by atan(1/sqrt(7)) = 20.68
       degrees anti-clockwise.  Each level rotates further and for example after about 18 levels
       it goes all the way around and back to the first quadrant.

       The rotation doesn't fill the plane though, only 1/3 of it.  The shape fattens as it curls
       around, but leaves a spiral gap beneath (see "Arms" below).

   Tiling
       The base pattern corresponds to a tiling by hexagons as follows, with the "***" lines
       being the base figure.

                 .     .
                / \   / \
               /   \ /   \
              .     .     .
              |     |     |
              |     |     |
              4*****5*****6
             /*\   / \   /*\
            / * \ /   \ / * \
           .   * .     .   * .
           |   * |     |    *|
           |    *|     |    *|
           .     3*****2     7...
            \   / \   /*\   /
             \ /   \ / * \ /
              .     . *   .
              |     | *   |
              |     |*    |
              0*****1     .
               \   / \   /
                \ /   \ /
                 .     .

       In the next level the parts corresponding to 1, 2 and 6 are reversed because they have
       their hexagon tile to the right of the line segment, rather than to the left.

   Arms
       The optional "arms" parameter can give up to three copies of the flowsnake, each advancing
       successively.  For example "arms=>3" is as follows.

           arms => 3                        51----48----45                 5
                                              \           \
                             ...   69----66    54----57    42              4
                               \     \     \        /     /
              28----25----22    78    72    63----60    39    33----30     3
                \           \     \  /                    \  /     /
                 31----34    19    75    12----15----18    36    27        2
                      /     /              \           \        /
              40----37    16     4---- 1     9---- 6    21----24           1
             /              \     \              /
           43    55----58    13     7     0---- 3    74----77---...    <- Y=0
             \     \     \     \  /                    \
              46    52    61    10     2     8----11    71----68          -1
                \  /     /              \  /     /           /
                 49    64    70----73     5    14    62----65             -2
                         \  /     /           /     /
                          67    76    20----17    59    53----50          -3
                               /     /              \  /     /
                             ...   23    35----38    56    47             -4
                                     \     \     \        /
                                      26    32    41----44                -5
                                        \  /
                                         29                               -6

                                          ^
              -9 -8 -7 -6 -5 -4 -3 -2 -1 X=0 1  2  3  4  5  6  7  8  9

       The N=3*k points are the plain curve, N=3*k+1 is a copy rotated by 120 degrees (1/3
       around), and N=3*k+2 is a copy rotated by 240 degrees (2/3 around).  The initial N=1 of
       the second arm and N=2 of the third correspond to N=3 of the first arm, rotated around.

       Essentially the flowsnake fills an ever expanding hexagon with one corner at the origin,
       and wiggly sides.

                   *---*
                  /     \
             *---*   A   *         Plain curve in A.
            /     \     /          Room for two more arms in B and C,
           *   B   O---*           rotated 120 and 240 degrees.
            \     /     \
             *---*   C   *
                  \     /
                   *---*

       The sides of these "hexagons" are not straight lines but instead wild wiggly spiralling S
       shapes, and the endpoints rotate around by the angle described above at each level.
       Opposing sides are symmetric, so they mesh perfectly and with three arms fill the plane.

   Fractal
       The flowsnake can also be thought of as successively subdividing line segments with
       suitably scaled copies of the 0 to 7 figure (or its reversal).

       The code here could be used for that by taking points N=0 to N=7^level.  The Y coordinates
       should be multiplied by sqrt(3) to make proper equilateral triangles, then a rotation and
       scaling to make the endpoint come out at some desired point, such as X=1,Y=0.  With such a
       scaling the path is confined to a finite fractal boundary.

FUNCTIONS

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

       "$path = Math::PlanePath::Flowsnake->new ()"
       "$path = Math::PlanePath::Flowsnake->new (arms => $a)"
           Create and return a new flowsnake path object.

           The optional "arms" parameter gives between 1 and 3 copies of the curve successively
           advancing.

       "($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.

       "($n_lo, $n_hi) = $path->rect_to_n_range ($x1,$y1, $x2,$y2)"
           In the current code the returned range is exact, meaning $n_lo and $n_hi are the
           smallest and biggest in the rectangle, but don't rely on that yet since finding the
           exact range is a touch on the slow side.  (The advantage of which though is that it
           helps avoid very big ranges from a simple over-estimate.)

   Level Methods
       "($n_lo, $n_hi) = $path->level_to_n_range($level)"
           Return "(0, 7**$level)", or for multiple arms return "(0, $arms * 7**$level)".

           There are 7^level + 1 points in a level, numbered starting from 0.  On the second and
           third arms the origin is omitted (so as not to repeat that point) and so just 7^level
           for them, giving 7^level+1 + (arms-1)*7^level = arms*7^level + 1 many points starting
           from 0.

FORMULAS

   N to X,Y
       The position of a given N can be calculated from the base-7 digits of N from high to low.
       At a given digit position the state maintained is

           direction 0 to 5, multiple of 60-degrees
           plain or reverse pattern

       It's convenient to work in the "i,j" coordinates per "Triangular Calculations" in
       Math::PlanePath.  This represents a point in the triangular grid as i+j*w where
       w=1/2+I*sqrt(3)/2 the a complex sixth root of unity at +60 degrees.

           foreach base-7 digit high to low
             (i,j) = (2i-j, i+3j)   # multiply by 2+w
             (i,j) += position of digit in plain or reverse,
                      and rotated by "direction"

       The multiply by 2+w scales up i,j by that vector, so for instance i=1,j=0 becomes i=2,j=1.
       This spreads the points as per the multiple-of-7 diagram shown above, so what was at N
       scales up to 7*N.

       The digit is then added as either the plain or reversed base figure,

             plain             reverse

           4-->5-->6
           ^       ^
            \       \
             3-->2   *             6<---*
                /                  ^
               v                  /
           0-->1             0   5<--4
                              \       \
                               v       v
                               1<--2<--3

       The arrow shown in each part is whether the state becomes plain or reverse.  For example
       in plain state at digit=1 the arrow points backwards so if digit=1 then the state changes
       to reverse for the next digit.  The direction likewise follows the direction of each
       segment in the pattern.

       Notice the endpoint "*" is at at 2+w in both patterns.  When considering the rotation it's
       convenient to reckon the direction by this endpoint.

       The combination of direction and plain/reverse is a total of 14 different states, and for
       each there's 7 digit values (0 to 6) for a total 84 i,j.

   X,Y to N
       The current approach uses the "FlowsnakeCentres" code.  The tiling in "Flowsnake" and
       "FlowsnakeCentres" is the same so the X,Y of a given N are no more than 1 away in the grid
       of the two forms.

       The way the two lowest shapes are arranged in fact means that if the Flowsnake N is at X,Y
       then the same N in "FlowsnakeCentres" is at one of three locations

           X, Y         same
           X-2, Y       left      (+180 degrees)
           X-1, Y-1     left down (+240 degrees)

       This is true even when the rotated "arms" multiple paths (the same number of arms in both
       paths).

       Is there an easy way to know which of the three offsets is right?  The current approach is
       to put each through "FlowsnakeCentres" to make an N, and put that N back through Flowsnake
       "n_to_xy()" to see if it's the target $n.

   Rectangle to N Range
       The current code calculates an exact "rect_to_n_range()" by searching for the highest and
       lowest Ns which are in the rectangle.

       The curve at a given level is bounded by the Gosper island shape but the wiggly sides make
       it difficult to calculate, so a bounding radius sqrt(7)^level, plus a bit, is used.  The
       portion of the curve comprising a given set of high digits of N can be excluded if the N
       point is more than that radius away from the rectangle.

       When a part of the curve is excluded it prunes a whole branch of the digits tree.  When
       the lowest digit is reached then a check for that point being actually within the
       rectangle is made.  The radius calculation is a bit rough, and it doesn't take into
       account the direction of the curve, so it's a rather large over-estimate, but it works.

       The same sort of search can be applied for highest and lowest N in a non-rectangular
       shapes, calculating a radial distance away from the shape.  The distance calculation
       doesn't have to be exact either, it can go from something bounding the shape until the
       lowest digit is reached and an individual X,Y is considered as a candidate for high or low
       N.

   N to Turn
       The turn made by the curve at a point N>=1 can be calculated from the lowest non-0 digit
       and the plain/reverse state per the lowest non-3 above there.

          N digits in base 7
          strip low 0-digits, zcount many of them
          zdigit = take next digit (lowest non-0)
          strip 3-digits
          tdigit = take next digit (0 if no digits left)
          plain if tdigit=0,4,5, reverse if tdigit=1,2,6

                                                zcount
                ---+--------+--------+--------+--------+
           high    | tdigit |   3s   | zdigit |   0s   |   low
                ---+--------+--------+--------+--------+

                    if zcount=0       if zcount>=1
                    ie. no low 0s     ie. some low 0s
          zdigit   plain reverse     plain   reverse
          ------   ----- -------     -----   -------
             1        1      1          1        1
             2        2      0          1       -1       turn left
             3       -1      2         -1        1       multiple of
             4       -2      1         -1        1       60 degrees
             5        0     -2          1       -1
             6        1     -1         -1       -1

       For example N=9079 is base-7 "35320" so a single low 0 for zcount=1 and strip it to
       "3532".  Take zdigit=2 leaving "353".  Skip low 3s leaving "35".  Take tdigit=5 which is
       "plain".  So table "plain" with zcount>=1 is the third column and there zdigit=2 is
       turn=+1.

       The turns in the zcount=0 "no low 0s" columns follow the turns of the base pattern shown
       above.  For example zdigit=1 is as per N=1 turning 120 degrees left, so +2.  For the
       reverse pattern the turns are negated and the zdigit value reversed, so the "reverse"
       column read 6 to 1 is the same as the plain column negated and read 1 to 6.

       Low 0s are stripped because the turn at a point such as N=7 ("10" in base-7) is determined
       by the pattern above it, the self-similar multiple-of-7s shape.  But when there's low 0s
       in the way there's an adjustment to apply because the last segment of the base pattern is
       not in the same direction as the first, but instead at -60 degrees.  Likewise the first
       segment of the reverse pattern.  At some zdigit positions those two cancel out, such as at
       zdigit=1 where a plain and reverse meet, but others it's not so and hence separate table
       columns for with or without low 0s.

       The plain or reverse pattern is determined by the lowest non-3 digit.  This works because
       the digit=0, digit=4, and digit=5 segments of the base pattern have their sub-parts
       "plain" in all cases, both the plain and reverse forms.  Conversely digit=1, digit=2 and
       digit=6 segments are "reverse" in all cases, both plain and reverse forms.  But the
       digit=3 part is plain in plain and reverse in reverse, so it inherits the orientation of
       the digit above and it's therefore necessary to skip across any 3s.

       When taking digits, N is treated as having infinite 0-digits at the high end.  This only
       affects the tdigit plain/reverse step.  If N has a single non-zero such as "5000" then
       it's taken as zdigit=5 and above that for the plain/reverse a tdigit=0 is then assumed.
       The first turn is at N=1 so there's always at least one non-0 for the zdigit.

OEIS

       Entries in Sloane's Online Encyclopedia of Integer Sequences related to this path include

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

           A229214   direction 1,2,3,-1,-2,-3 (clockwise)

SEE ALSO

       Math::PlanePath, Math::PlanePath::FlowsnakeCentres, Math::PlanePath::GosperIslands

       Math::PlanePath::KochCurve, Math::PlanePath::HilbertCurve, Math::PlanePath::PeanoCurve,
       Math::PlanePath::ZOrderCurve

       Fukuda, Shimizu and Nakamura, "New Gosper Space Filling Curves", on flowsnake variations
       in bigger hexagons (with wiggly sides too).

           <http://kilin.clas.kitasato-u.ac.jp/museum/gosperex/343-024.pdf> or
           <http://web.archive.org/web/20070630031400/http://kilin.u-shizuoka-ken.ac.jp/museum/gosperex/343-024.pdf>

HOME PAGE

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

LICENSE

       Copyright 2010, 2011, 2012, 2013, 2014, 2015 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/>.