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 selfsimilar Koch curve,
 Helge von Koch, ``Une Méthode Géométrique Élémentaire pour l'Étude de Certaines Questions de la Théorie des Courbes Planes'', Acta Arithmetica, volume 30, 1906, pages 145174. <http://archive.org/details/actamathematica11lefgoog>
It goes along the X axis and makes triangular excursions upwards.
8 3 / \ 6 7 910 18... 2 \ / \ 2 5 11 14 17 1 / \ / \ / \ / 01 3 4 1213 1516 < 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 inbetween 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 unitside 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 selfsimilar 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[level1] + 4^(level1) = 4^(level1) + 9*4^(level2) + ... + 9^(level1)*4^0 9^level  4^level =  5 = 0, 1, 13, 133, 1261, 11605, 105469, ... (A016153)
The sides are 6 different angles. The triangles added on the sides 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.
Level Methods
 "($n_lo, $n_hi) = $path>level_to_n_range($level)"
 Return "(0, 4**$level)".
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 nonzero digit gives the turn. The first turn is at N=1 so there's always a nonzero 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 nonzero 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 0bits. A low base4 digit 1 or 3 is an even number of low 0bits and a low digit 2 is an odd number of low 0bits.
count low 0bits turn   even +60 degrees (left) odd 120 degrees (right)
For example N=8 in binary ``1000'' has 3 low 0bits 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 base4 digits by considering how the digits of N change when 1 is added, and the lowdigit turn calculation is applied on those changed digits.Adding 1 means low digit 0, 1 or 2 will become nonzero. 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 non3
lowest non3 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 1digits in base 4)  (count of 2digits in base 4)) * 60 degrees
For example N=11 is ``23'' in base 4, so 60*(01) = 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 NorthEast or SouthWest 2 1,1 slope NorthWest or SouthEast
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 powerof4 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 overestimate 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 roundingdown 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 overestimate 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 overestimate described above.At a given digit position in the prospective N the subpart 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 subpart. 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 subcurve 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 subcurve. 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 subpart extent checking reduces the subparts which are examined, but it works perfectly well with a looser check, such as a square box for the subcurve 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, ThueMorse 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^k1)/3 A002450 number of right turns N=0 to N < 4^k, being (4^k1)/3 A016153 area under the curve, (9^n4^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.
LICENSE
Copyright 2011, 2012, 2013, 2014, 2015, 2016 Kevin RydeMathPlanePath 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.
MathPlanePath 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 MathPlanePath. If not, see <http://www.gnu.org/licenses/>.