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 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, ...
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 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.
SEE ALSO
Math::PlanePath, Math::PlanePath::PeanoCurve, Math::PlanePath::HilbertCurve, Math::PlanePath::KochPeaks, Math::PlanePath::KochSnowflakes, Math::PlanePath::CCurve
HOME PAGE
http://user42.tuxfamily.org/mathplanepath/index.html
LICENSE
Copyright 2011, 2012, 2013, 2014 Kevin Ryde
MathPlanePath 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/>.