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

NAME

       Math::PlanePath::KochCurve -- horizontal Koch curve

SYNOPSIS

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

DESCRIPTION

       This is an integer version of the self-similar Koch curve,

           Helge von Koch, "Une Methode Geometrique Elementaire pour l'Etude de Certaines
           Questions de la Theorie des Courbes Planes", Acta Arithmetica, volume 30, 1906, pages
           145-174.  <http://archive.org/details/actamathematica11lefgoog>

       It goes along the X axis and makes triangular excursions upwards.

                                      8                                   3
                                    /  \
                             6---- 7     9----10                18-...    2
                              \              /                    \
                    2           5          11          14          17     1
                  /  \        /              \        /  \        /
            0----1     3---- 4                12----13    15----16    <- Y=0

            ^
           X=0   2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19

       The replicating shape is the initial N=0 to N=4,

                   *
                  / \
             *---*   *---*

       which is rotated and repeated 3 times in the same pattern to give sections N=4 to N=8, N=8
       to N=12, and N=12 to N=16.  Then that N=0 to N=16 is itself replicated three times at the
       angles of the base pattern, and so on infinitely.

       The X,Y coordinates are arranged on a square grid using every second point, per
       "Triangular Lattice" in Math::PlanePath.  The result is flattened triangular segments with
       diagonals at a 45 degree angle.

   Level Ranges
       Each replication adds 3 copies of the existing points and is thus 4 times bigger, so if
       N=0 to N=4 is reckoned as level 1 then a given replication level goes from

           Nstart = 0
           Nlevel = 4^level   (inclusive)

       Each replication is 3 times the width.  The initial N=0 to N=4 figure is 6 wide and in
       general a level runs from

           Xstart = 0
           Xlevel = 2*3^level   at N=Nlevel

       The highest Y is 3 times greater at each level similarly.  The peak is at the midpoint of
       each level,

           Npeak = (4^level)/2
           Ypeak = 3^level
           Xpeak = 3^level

       It can be seen that the N=6 point backtracks horizontally to the same X as the start of
       its section N=4 to N=8.  This happens in the further replications too and is the maximum
       extent of the backtracking.

       The Nlevel is multiplied by 4 to get the end of the next higher level.  The same 4*N can
       be applied to all points N=0 to N=Nlevel to get the same shape but a factor of 3 bigger
       X,Y coordinates.  The in-between points 4*N+1, 4*N+2 and 4*N+3 are then new finer
       structure in the higher level.

   Fractal
       Koch conceived the curve as having a fixed length and infinitely fine structure, making it
       continuous everywhere but differentiable nowhere.  The code here can be pressed into use
       for that sort of construction for a given level of granularity by scaling

           X/3^level
           Y/3^level

       which makes it a fixed 2 wide by 1 high.  Or for unit-side equilateral triangles then
       apply further factors 1/2 and sqrt(3)/2, as noted in "Triangular Lattice" in
       Math::PlanePath.

           (X/2) / 3^level
           (Y*sqrt(3)/2) / 3^level

   Area
       The area under the curve to a given level can be calculated from its self-similar nature.
       The curve at level+1 is 3 times wider and higher and adds a triangle of unit area onto
       each line segment.  So reckoning the line segment N=0 to N=1 as level=0 (which is
       area[0]=0),

           area[level] = 9*area[level-1] + 4^(level-1)
                       = 4^(level-1) + 9*4^(level-2) + ... + 9^(level-1)*4^0

                         9^level - 4^level
                       = -----------------
                                 5

                       = 0, 1, 13, 133, 1261, 11605, 105469, ...

       The sides are 6 different angles.  The triangles added there are always the same shape
       either pointing up or down.  Base width=2 and height=1 gives area=1.

              *            *-----*   ^
             / \            \   /    | height=1
            /   \            \ /     |
           *-----*            *      v     triangle area = 2*1/2 = 1

           <-----> width=2

       If the Y coordinates are stretched to make equilateral triangles then the number of
       triangles is not changed and so the area increases by a factor of the area of the
       equilateral triangle, sqrt(3)/4.

FUNCTIONS

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

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

       "($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)"
           The returned range is exact, meaning $n_lo and $n_hi are the smallest and biggest in
           the rectangle.

       "$n = $path->n_start()"
           Return 0, the first N in the path.

FORMULAS

   N to Turn
       The curve always turns either +60 degrees or -120 degrees, it never goes straight ahead.
       In the base 4 representation of N the lowest non-zero digit gives the turn.  The first
       turn is at N=1 so there's always a non-zero digit in N.

          low digit
           base 4         turn
          ---------   ------------
             1         +60 degrees (left)
             2        -120 degrees (right)
             3         +60 degrees (left)

       For example N=8 is 20 base 4, so lowest nonzero "2" means turn -120 degrees for the next
       segment.

       If the least significant digit is non-zero then it determines the turn, making the base
       N=0 to N=4 shape.  If the least significant is zero then the next level up is in control,
       eg. N=0,4,8,12,16, making a turn according to the base shape again at that higher level.
       The first and last segments of the base shape are "straight" so there's no extra
       adjustment to apply in those higher digits.

       This base 4 digit rule is equivalent to counting low 0-bits.  A low base-4 digit 1 or 3 is
       an even number of low 0-bits and a low digit 2 is an odd number of low 0-bits.

           count low 0-bits         turn
           ----------------     ------------
                even             +60 degrees (left)
                odd             -120 degrees (right)

       For example N=8 in binary "1000" has 3 low 0-bits and 3 is odd so turn -120 degrees
       (right).

       See "Turn" in Math::PlanePath::GrayCode for a similar turn sequence arising from binary
       Gray code.

   N to Next Turn
       The turn at N+1, ie the next turn, can be found from the base-4 digits by considering how
       the digits of N change when 1 is added, and the low-digit turn calculation is applied on
       those changed digits.

       Adding 1 means low digit 0, 1 or 2 will become non-zero.  Any low 3s wrap around to become
       low 0s.  So the turn at N+1 can be found from the digits of N by seeking the lowest non-3

          lowest non-3       turn
           digit of N       at N+1
          ------------   ------------
               0          +60 degrees (left)
               1         -120 degrees (right)
               2          +60 degrees (left)

   N to Direction
       The total turn at a given N can be found by counting digits 1 and 2 in base 4.

           direction = ((count of 1-digits in base 4)
                        - (count of 2-digits in base 4)) * 60 degrees

       For example N=11 is "23" in base 4, so 60*(0-1) = -60 degrees.

       In this formula the count of 1s and 2s can go past 360 degrees, representing a spiralling
       around which occurs at progressively higher replication levels.  The direction can be
       taken mod 360 degrees, or the count mod 6, for a direction 0 to 5 if desired.

   N to abs(dX),abs(dY)
       The direction expressed as abs(dX) and abs(dY) can be calculated simply from N modulo 3.
       abs(dX) is a repeating pattern 2,1,1 and abs(dY) repeating 0,1,1.

           N mod 3     abs(dX),abs(dY)
           -------     ---------------
              0             2,0            horizontal, East or West
              1             1,1            slope North-East or South-West
              2             1,1            slope North-West or South-East

       This works because the direction calculation above corresponds to N mod 3.  Each N digit
       in base 4 becomes

           N digit
           base 4    direction add
           -------   -------------
              0            0
              1            1
              2           -1
              3            0

       Notice that direction == Ndigit mod 3.  Then because 4==1 mod 3 the power-of-4 for each
       digit reduces down to 1,

           N = 4^k * digit_k + ... 4^0 * digit_0
           N mod 3 = 1 * digit_k + ... 1 * digit_0
                   = digit_k + ... digit_0
           same as
           direction = digit_k + ... + digit_0    taken mod 3

   Rectangle to N Range -- Level
       An easy over-estimate of the N values in a rectangle can be had from the Xlevel formula
       above.  If Xlevel>rectangleX then Nlevel is past the rectangle extent.

           X = 2*3^level

       so

           floorlevel = floor log_base_3(X/2)
           Nhi = 4^(floorlevel+1) - 1

       For example a rectangle extending to X=13 has floorlevel = floor(log3(13/2))=1 and so
       Nhi=4^(1+1)-1=15.

       The rounding-down of the log3 ensures a point such as X=18 which is the first in the next
       Nlevel will give that next level.  So floorlevel=log3(18/2)=2 (exactly) and
       Nhi=4^(2+1)-1=63.

       The worst case for this over-estimate is when rectangleX==Xlevel, ie. the rectangle is
       just into the next level.  In that case Nhi is nearly a factor 4 bigger than it needs to
       be.

   Rectangle to N Range -- Exact
       The exact Nlo and Nhi in a rectangle can be found by searching along the curve.  For Nlo
       search forward from the origin N=0.  For Nhi search backward from the Nlevel over-estimate
       described above.

       At a given digit position in the prospective N the sub-part of the curve comprising the
       lower digits has a certain triangular extent.  If it's outside the target rectangle then
       step to the next digit value, and to the next of the digit above when past digit=3 (or
       below digit=0 when searching backwards).

       There's six possible orientations for the curve sub-part.  In the following pictures "o"
       is the start and the surrounding lines show the triangular extent.  There's just four
       curve parts shown in each, but these triangles bound a sub-curve of any level.

          rot=0   -+-               +-----------------+
                --   --              - .-+-*   *-+-o -
              --   *   --             --    \ /    --
            --    / \    --             --   *   --
           - o-+-*   *-+-. -              --   --
          +-----------------+       rot=3   -+-

          rot=1
          +---------+               rot=4    /+
          |      . /                        / |
          |     / /                        / o|
          |*-+-* /                        / / |
          | \   /                        / *  |
          |  * /                        /   \ |
          | / /                        / *-+-*|
          |o /                        / /     |
          | /                        / .      |
          +/                        +---------+

          +\  rot=2                 +---------+
          | \                        \ o      |
          |. \                        \ \     |
          | \ \                        \ *-+-*|
          |  * \                        \   / |
          | /   \                        \ *  |
          |*-+-* \                        \ \ |
          |     \ \                        \ .|
          |      o \                rot=5   \ |
          +---------+                        \+

       The "." is the start of the next sub-curve.  It belongs to the next digit value and so can
       be excluded.  For rot=0 and rot=3 this means simply shortening the X range permitted.  For
       rot=1 and rot=4 similarly the Y range.  For rot=2 and rot=5 it would require a separate
       test.

       Tight sub-part extent checking reduces the sub-parts which are examined, but it works
       perfectly well with a looser check, such as a square box for the sub-curve extents.  Doing
       that might be easier if the target region is not a rectangle but instead some trickier
       shape.

OEIS

       The Koch curve is in Sloane's Online Encyclopedia of Integer Sequences in various forms,

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

           A177702   abs(dX) from N=1 onwards, being 1,1,2 repeating
           A011655   abs(dY), being 0,1,1 repeating
           A035263   turn 1=left,0=right, by morphism
           A096268   turn 0=left,1=right, by morphism
           A056832   turn 1=left,2=right, by replicate and flip last
           A029883   turn +/-1=left,0=right, Thue-Morse first differences
           A089045   turn +/-1=left,0=right, by +/- something

           A003159   N positions of left turns, ending even number 0 bits
           A036554   N positions of right turns, ending odd number 0 bits

           A020988   number of left turns N=0 to N < 4^k, being 2*(4^k-1)/3
           A002450   number of right turns N=0 to N < 4^k, being (4^k-1)/3
           A016153   area under the curve, (9^n-4^n)/5

       For reference, A217586 is not quite the same as A096268 right turn.  A217586 differs by a
       0<->1 flip at N=2^k due to different initial a(1)=1.

SEE ALSO

       Math::PlanePath, Math::PlanePath::PeanoCurve, Math::PlanePath::HilbertCurve,
       Math::PlanePath::KochPeaks, Math::PlanePath::KochSnowflakes, Math::PlanePath::CCurve

       Math::Fractal::Curve

HOME PAGE

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

LICENSE

       Copyright 2011, 2012, 2013 Kevin Ryde

       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/>.