Math::PlanePath::KochCurve(3) 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 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 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, ...  (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)".

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