NAME
Math::PlanePath::CCurve -- Levy C curve
SYNOPSIS
use Math::PlanePath::CCurve;
my $path = Math::PlanePath::CCurve->new;
my ($x, $y) = $path->n_to_xy (123);DESCRIPTION
This is an integer version of the Levy "C" curve.
11-----10-----9,7-----6------5 3
| | |
13-----12 8 4------3 2
| |
19---14,18----17 2 1
| | | |
21-----20 15-----16 0------1 <- Y=0
|
22 -1
|
25,23---24 -2
|
26 35-----34-----33 -3
| | |
27,37--28,36 32 -4
| | |
38 29-----30-----31 -5
|
39,41---40 -6
|
42 ... -7
| |
43-----44 49-----48 64-----63 -8
| | | |
45---46,50----47 62 -9
| |
51-----52 56 60-----61 -10
| | |
53-----54----55,57---58-----59 -11
^
-7 -6 -5 -4 -3 -2 -1 X=0 1The initial segment N=0 to N=1 is repeated with a turn +90 degrees left to give N=1 to N=2. Then N=0to2 is repeated likewise turned +90 degrees and placed at N=2 to make N=2to4. And so on doubling each time.
The 90 degree rotation is the same at each repetition, so the segment at N=2^k is the initial N=0to1 turned +90 degrees. This means at N=1,2,4,8,16,etc the direction is always upwards.
If 2^k is the highest 1-bit in N then the X,Y position can be written in complex numbers as
XY(N) = XY(2^k) + i*XY(r) N = 2^k + r with r<2^k
= (1+i)^k + i*XY(r)The effect is a change of base from binary to base 1+i but with a power of i on each term. Suppose the 1-bits in N are at positions k, k1, k2, etc, then
XY(N) = b^k N= 2^k + 2^(k1) + 2^(k2) + ... in binary
+ b^k1 * i base b=1+i
+ b^k2 * i^2
+ b^k3 * i^3
+ ...Notice the power of i is not the bit position k, but rather the count of how many 1-bits are above the position. This calculation is straightforward but the resulting structure of boundary and shapes enclosed has many different parts.
Level Ranges 4^k
The X,Y extents of the path through to Nlevel=2^k can be expressed as a width and height measured relative to the endpoints.
*------------------* <-+
| | |
*--* *--* | height h[k]
| | |
* N=4^k N=0 * <-+
| | | | | below l[k]
*--*--* *--*--* <-+
^-----^ ^-----^ Extents to N=4^k
width 2^k width
w[k] w[k]
<------------------------>
total width -> 2N=4^k is on either the X or Y axis and for the extents here it's taken rotated as necessary to be horizontal. k=2 N=4^2=16 shown above is already horizontal. The next level k=3 N=64=4^3 would be rotated -90 degrees to be horizontal.
The width w[k] is measured from the N=0 and N=4^k endpoints. It doesn't include the 2^k length between those endpoints. The two ends are symmetric so the extent is the same for each.
h[k] = 2^k - 1 0,1,3,7,15,31,etc
w[k] = / 0 for k=0
\ 2^(k-1) - 1 for k>=1 0,0,1,3,7,15,etc
l[k] = / 0 for k<=1
\ 2^(k-2) - 1 for k>=2 0,0,0,1,3,7,etcThe initial N=0 to N=0 to N=64 shown above is k=3. h[3]=7 is the X=-7 horizontal. l[3]=1 is the X=1 horizontal. w[3]=3 is the vertical Y=3, and also Y=-11 which is 3 below the endpoint N=64 at Y=8.
Expressed as a fraction of the 2^k distance between the endpoints the extents approach total 2 wide by 1.25 high,
*------------------* <-+
| | | 1
*--* *--* | total
| | | height
* N=4^k N=0 * <-+ 1+1/4
| | | | | 1/4
*--*--* *--*--* <-+
^-----^ ^-----^
1/2 1 1/2 total width 2The extent formulas can be found by considering the self-similar blocks. The initial k=0 is a single line segment and all its extents are 0.
h[0] = 0
N=1 ----- N=0
l[0] = 0
w[0] = 0Thereafter the replication overlap as
+-------+---+-------+
| | | |
+------+ | | +------+
| | D | | C | | B | | <-+
| +-------+---+-------+ | | 2^(k-1)
| | | | | previous
| | | | | level ends
| E | | A | <-+
+------+ +------+
^---------------^
2^k this level ends
w[k] = max (h[k-1], w[k-1]) # right of A,B
h[k] = 2^(k-1) + max (h[k-1], w[k-1]) # above B,C,D
l[k] = max w[k-1], l[k-1]-2^(k-1) # below A,ESince h[k]=2^(k-1)+w[k] have h[k] > w[k] for k>=1 and with the initial h[0]=w[k]=0 have h[k]>=w[k] always. So the max of those two is h.
h[k] = 2^(k-1) + h[k-1] giving h[k] = 2^k-1 for k>=1
w[k] = h[k-1] giving w[k] = 2^(k-1)-1 for k>=1The max for l[k] is always w[k-1] as l[k] is never big enough that the parts B-C and C-D can extend down past their 2^(k-1) vertical position. (l[0]=w[0]=0 and thereafter by induction l[k]<=w[k].)
l[k] = w[k-1] giving l[k] = 2^(k-2)-1 for k>=2Repeated Points
The curve crosses itself and repeats some X,Y positions up to 4 times. The first doubled, tripled and quadrupled points are
visits first X,Y N
--------- --------- ----------------------
2 -2, 3 7, 9
3 18, -7 189, 279, 281
4 -32, 55 1727, 1813, 2283, 2369Each line segment between integer points is traversed at most 2 times, once forward and once backward. There's 4 such lines reaching each integer point and so the points are visited at most 4 times.
As per "Direction" below the direction of the curve is given by the count of 1-bits in N. Since no line is repeated each of the N values at a given X,Y have a different count 1-bits mod 4. For example N=7 is 3 1-bits and N=9 is 2 1-bits. The full counts need not be consecutive, as for example N=1727 is 9 1-bits and N=2369 is 4 1-bits.
The maximum 2 segment traversals can be seen from the way the curve replicates. Suppose the entire plane had all line segments traversed forward and backward.
v | v |
-- <-------- <-
[0,1] [1,1] [X,Y] = integer points
-> --------> -- each edge traversed
| ^ | ^ forward and backward
| | | |
| | | |
v | v |
-- <-------- <--
[0,0] [1,0]
-> --------> --
| ^ | ^Then when each line segment expands on the right the result is the same pattern of traversals when viewed rotated by 45-degrees and scaled by factor sqrt(2).
\ v / v \ v / v
[0,1] [1,1]
/ / ^ \ ^ / ^ \
/ / \ \ / / \ \
\ \ / /
\ v / v
[1/2,1/2]
^ / ^ \
/ / \ \
\ \ / / \ \ / /
\ v / v \ v / v
[0,0] 1,0
^ / ^ \ ^ / ^ \The curve is a subset of this pattern. It begins as a single line segment which has this pattern and thereafter the pattern preserves itself. Hence at most 2 segment traversals in the curve.
Tiling
The segment traversal argument above can also be made by taking the line segments as triangles which are a quarter of a unit square with peak pointing to the right of the traversal direction.
to *
^\
| \
| \ triangle peak
| /
| /
|/
from *These triangles in the two directions tile the plane. On expansion each splits into 2 halves in new positions. Those parts don't overlap and the plane is still tiled. See for example Larry Riddle's pages
http://www.agnesscott.edu/lriddle/levy.html http://www.agnesscott.edu/lriddle/tiling.html
For the integer version of the curve this kind of tiling can be used to combine copies of the curve so that each every point is visited precisely 4 times. The h[k], w[k] and l[k] extents above are less than the 2^k endpoint length, so a square of side 2^k can be fully tiled with copies of the curve at each corner,
| ^ | ^
| | | | 24 copies of the curve
| | | | to visit all points of the
v | v | inside square precisely
<------- <-------- <-------- 4 times each
* *
--------> --------> --------> points N=0 to N=4^k-1
| ^ | ^ rotated and shifted
| | | | suitably
| | | |
v | v |
<-------- <-------- <--------
* *
-------- --------> -------->
| ^ | ^
| | | |
| | | |
v | v |The four innermost copies of the curve cover most of the inside square, but the other copies surrounding them loop into the square and fill in the remainder to make 4 visits at every point.
It's interesting to note that a set of 8 curves at the origin only covers the axes with 4-fold visits,
_ _ _
| ^ 8 arms at the origin
| | cover only X,Y axes
v | with 4-visits
<-------- <--------
0,0 away from the axes
-------- --------> some points < 4 visits
| ^
| |
v |The "_ _ _" line shown which is part of the 24-pattern above but omitted here. This line is at Y=2^k. The extents described above mean that it extends down to Y=2^k - h[k] = 2^k-(2^k-1)=1, so it visits some points in row Y=1 and higher. Omitting the curve means there are Y>=1 not visited 4 times. Similarly Y<=-1 and X<-1 and X>=+1.
This means that if the path had some sort of "arms" of multiple curves extending from the origin then it would visit all points on the axes X=0 Y=0 a full 4 times, but there would be infinitely many points off the axes without full 4 visits.
FUNCTIONS
See "FUNCTIONS" in Math::PlanePath for the behaviour common to all path classes.
$path = Math::PlanePath::CCurve->new ()-
Create and return a new path object.
($x,$y) = $path->n_to_xy ($n)-
Return the X,Y coordinates of point number
$non the path. Points begin at 0 and if$n < 0then the return is an empty list.Fractional positions give an X,Y position along a straight line between the integer positions.
$n = $path->xy_to_n ($x,$y)-
Return the point number for coordinates
$x,$y. If there's nothing at$x,$ythen returnundef. If$x,$yis visited more than once then return the smallest$nwhich visits it. @n_list = $path->xy_to_n_list ($x,$y)-
Return a list of N point numbers at coordinates
$x,$y. If there's nothing at$x,$ythen return an empty list.A given
$x,$yis visited at most 4 times so the returned list is at most 4 values. $n = $path->n_start()-
Return 0, the first N in the path.
FORMULAS
Direction
The direction or net turn of the curve is the count of 1 bits in N,
direction = count_1_bits(N) * 90degreesFor example N=11 is binary 1011 has three 1 bits, so direction 3*90=270 degrees, ie. to the south.
This bit count is because at each power-of-2 position the curve is a copy of the lower bits but turned +90 degrees, so +90 for each 1-bit.
For powers-of-2 N=2,4,8,16, etc, there's only a single 1-bit so the direction is always +90 degrees there, ie. always upwards.
Turn
At each point N the curve can turn in any direction: left, right, straight, or 180 degrees back. The turn is given by the number of low 0-bits of N,
turn right = (count_low_0_bits(N) - 1) * 90degreesFor example N=8 is binary 0b100 which is 2 low 0-bits for turn=(2-1)*90=90 degrees to the right.
When N is odd there's no low zero bits and the turn is always (0-1)*90=-90 to the right, so every second turn is 90 degrees to the left.
Next Turn
The turn at the point following N, ie. at N+1, can be calculated by counting the low 1-bits of N,
next turn right = (count_low_1_bits(N) - 1) * 90degreesFor example N=11 is binary 0b1011 which is 2 low one bits for nextturn=(2-1)*90=90 degrees to the right at the following point, ie. at N=12.
This works simply because low 1-bits like ..0111 increment to low 0-bits ..1000 to become N+1. The low 1-bits at N are thus the low 0-bits at N+1.
N to dX,dY
n_to_dxdy() is implemented using the direction described above. For integer N the count mod 4 gives the direction for dX,dY.
dir = count_1_bits(N) mod 4
dx = dir_to_dx[dir] # table 0 to 3
dy = dir_to_dy[dir]For fractional N the direction at int(N)+1 can be obtained from combining the direction at int(N) and the turn at int(N)+1, that being the low 1-bits of N per "Next Turn" above. Those two directions can then be combined as described in "N to dX,dY -- Fractional" in Math::PlanePath.
# apply turn to make direction at Nint+1
turn = count_low_1_bits(N) - 1 # N integer part
dir = (dir - turn) mod 4 # direction at N+1
# adjust dx,dy by fractional amount in this direction
dx += Nfrac * (dir_to_dx[dir] - dx)
dy += Nfrac * (dir_to_dy[dir] - dy)A small optimization can be made by working the "-1" of the turn formula into a +90 degree rotation of the dir_to_dx[] and dir_to_dy[] parts by swap and sign change,
turn_plus_1 = count_low_1_bits(N) # on N integer part
dir = (dir - turn_plus_1) mod 4 # direction-1 at N+1
# adjustment including extra +90 degrees on dir
dx -= $n*(dir_to_dy[dir] + dx)
dy += $n*(dir_to_dx[dir] - dy)X,Y to N
The N values at a given X,Y can be found by taking terms low to high from the complex number formula (as given above),
X+iY = b^k N = 2^k + 2^(k1) + 2^(k2) + ... in binary
+ b^k1 * i base b=1+i
+ b^k2 * i^2
+ ...If the lowest term is b^0 then X+iY has X+Y odd. If the lowest term is not b^0 but instead some power b^n then X+iY has X+Y even. This is because a multiple of b=1+i,
X+iY = (x+iy)*(1+i)
= (x-y) + (x+y)i
so X=x-y Y=x+y
sum X+Y = 2x is even if X+iY a multiple of 1+iSo the lowest bit of N is found by
bit = (X+Y) mod 2If bit=1 then a power i^p is to be subtracted from X+iY. p is how many 1-bits are above that point, and this is not yet known. It represents a direction to move X,Y to put it on an even position. It's also the direction of the step N-2^l to N, where 2^l is the lowest 1-bit of N.
The reduction should be attempted with p as each of the four possible directions N,S,E,W. Some or all will lead to an N. For quadrupled points (such as X=-32, Y=55 described above) all four will lead to an N.
for p 0 to 3
dX,dY = i^p # directions [1,0] [0,1] [-1,0] [0,-1]
loop until X,Y = [0,0] or [1,0] or [-1,0] or [0,1] or [0,-1]
{
bit = X+Y mod 2 # bits of N from low to high
if bit == 1 {
X -= dX # move to "even" X+Y == 0 mod 2
Y -= dY
(dX,dY) = (dY,-dX) # rotate -90
}
(X,Y) = (X+Y)/2, (Y-X)/2 # divide (X+iY)/(1+i)
}
if not (dX=1 and dY=0)
wrong final direction, try next p
if X=dX and Y=dY
further high 1-bit for N
found an N
if X=0 and Y=0
found an NThe loop ends at one of the five points
0,1
|
-1,0 -- 0,0 -- 1,0
|
0,-1It's not possible to wait for X=0,Y=0 to be reached because some dX,dY directions will step infinitely among the four non-zeros. Only the case X=dX,Y=dY is sure to reach 0,0.
The successive p decrements which are dX,dY rotate -90 must end at p == 0 mod 4 for highest term in the X+iY formula having i^0=1. This means must end dX=1,dY=0 East.
The number of 1-bits in N is == p mod 4. So the order the N values are obtained follows the order the p directions are attempted. In general the N values will not be smallest to biggest N so a little sort is necessary if that's desired.
It can be seen that sum X+Y is used for the bit calculation and then again in the divide by 1+i. It's convenient to write the whole loop in terms of sum S=X+Y and difference D=Y-X.
for dS = +1 or -1 # four directions
for dD = +1 or -1 #
S = X+Y
D = Y-X
loop until -1 <= S <= 1 and -1 <= D <= 1 {
bit = S mod 2 # bits of N from low to high
if bit == 1 {
S -= dS # move to "even" S+D == 0 mod 2
D -= dD
(dS,dD) = (dD,-dS) # rotate -90
}
(S,D) = (S+D)/2, (D-S)/2 # divide (S+iD)/(1+i)
}
if not (dS=1 and dD=-1)
wrong final direction, try next dS,dD direction
if S=dS and D=dD
further high 1-bit for N
found an N
if S=0 and D=0
found an NThe effect of S=X+Y, D=Y-D is to rotate by -45 degrees and use every second point of the plane.
D= 2 X=0,Y=2 . rotate -45
D= 1 X=0,Y=1 . X=1,Y=2 .
D= 0 X=0,Y=0 . X=1,Y=1 . X=2,Y=2
D=-1 X=1,Y=0 . X=2,Y=1 .
D=-2 X=2,Y=0 .
S=0 S=1 S=2 S=3 S=4The final five points described above are then in a 3x3 block at the origin. The four in-between points S=0,D=1 etc don't occur so ranges tests -1<=S<=1 and -1<=D<=1 can be used.
S=-1,D=1 . S=1,D=1
. S=0,D=0 .
S=-1,D=-1 . S=1,D=-1OEIS
Entries in Sloane's Online Encyclopedia of Integer Sequences related to this path include
http://oeis.org/A179868 (etc)
A010059 abs(dX), count1bits(N) mod 2
A010060 abs(dY), count1bits(N)+1 mod 2, being Thue-Morse
A000120 direction, being total turn, count 1-bits
A179868 direction 0to3, count 1-bits mod 4
A035263 turn 0=straight or 180, 1=left or right,
being (count low 0-bits + 1) mod 2
A096268 next turn 1=straight or 180, 0=left or right,
being count low 1-bits mod 2
A007814 turn-1 to the right,
being count low 0-bits
A003159 N positions of left or right turn, ends even num 0 bits
A036554 N positions of straight or 180 turn, ends odd num 0 bits
A146559 X at N=2^k, being Re((i+1)^k)
A009545 Y at N=2^k, being Im((i+1)^k)
A191689 fractal dimension of the boundarySEE ALSO
Math::PlanePath, Math::PlanePath::DragonCurve, Math::PlanePath::AlternatePaper, Math::PlanePath::KochCurve
ccurve(6x) back-end of xscreensaver(1) displaying the C curve (and various other dragon curve and Koch curves).
HOME PAGE
http://user42.tuxfamily.org/math-planepath/index.html
LICENSE
Copyright 2011, 2012, 2013, 2014 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/>.