1Math::PlanePath::CCurveU(s3e)r Contributed Perl DocumentaMtaitohn::PlanePath::CCurve(3)
2
6 Math::PlanePath::CCurve -- Levy C curve
7
9 use Math::PlanePath::CCurve;
10 my $path = Math::PlanePath::CCurve->new;
11 my ($x, $y) = $path->n_to_xy (123);
12
14 This is an integer version of the "C" curve by Lévy.
15 "Les Courbes Planes ou Gauches et les Surfaces Composée de Parties
17 Semblables au Tout", Journal de l'École Polytechnique, July 1938
18 pages 227-247 and October 1938 pages 249-292
19 <http://gallica.bnf.fr/ark:/12148/bpt6k57344323/f53.image>
21 <http://gallica.bnf.fr/ark:/12148/bpt6k57344820>
22 It spirals anti-clockwise, variously crossing and overlapping itself.
24 The construction is straightforward but various measurements like how
25 many distinct points are quite complicated.
26 11-----10-----9,7-----6------5 3
28 | | |
29 13-----12 8 4------3 2
30 | |
31 19---14,18----17 2 1
32 | | | |
33 21-----20 15-----16 0------1 <- Y=0
34 |
35 22 -1
36 |
37 25,23---24 -2
38 |
39 26 35-----34-----33 -3
40 | | |
41 27,37--28,36 32 -4
42 | | |
43 38 29-----30-----31 -5
44 |
45 39,41---40 -6
46 |
47 42 ... -7
48 | |
49 43-----44 49-----48 64-----63 -8
50 | | | |
51 45---46,50----47 62 -9
52 | |
53 51-----52 56 60-----61 -10
54 | | |
55 53-----54----55,57---58-----59 -11
56 ^
58 -7 -6 -5 -4 -3 -2 -1 X=0 1
59 The initial segment N=0 to N=1 is repeated with a turn +90 degrees left
61 to give N=1 to N=2. Then N=0to2 is repeated likewise turned +90
62 degrees and placed at N=2 to make N=2to4. And so on doubling each
63 time.
64 4----3
66 | N=0to2
67 2 2 repeated
68 | | as N=2to4
69 0----1 0----1 0----1 with turn +90
70 The 90 degree rotation is the same at each repetition, so the segment
72 at N=2^k is always the initial N=0to1 turned +90 degrees. This means
73 at N=1,2,4,8,16,etc the direction is always upwards.
74 The X,Y position can be written in complex numbers as a recurrence
76 with N = 2^k + r high bit 2^k, rest r<2^k
78 C(N) = C(2^k) + i*C(r)
80 = (1+i)^k + i*C(r)
81 The effect is a change from base 2 to base 1+i but with a further power
83 of i on each term. Suppose the 1-bits in N are at positions k0, k1,
84 k2, etc (high to low), then
85 C(N) = b^k0 * i^0 N= 2^k0 + 2^(k1) + 2^(k2) + ... in binary
87 + b^k1 * i^1 k0 > k1 > k2 > ...
88 + b^k2 * i^2 base b=1+i
89 + b^k3 * i^3
90 + ...
91 Notice the i power is not the bit position k, but rather how many
93 1-bits are above the position.
94 Level Ranges 4^k
96 The X,Y extents of the path through to Nlevel=2^k can be expressed as a
97 width and height measured relative to the endpoints.
98 *------------------* <-+
100 | | |
101 *--* *--* | height h[k]
102 | | |
103 * N=4^k N=0 * <-+
104 | | | | | below l[k]
105 *--*--* *--*--* <-+
106 ^-----^ ^-----^
108 width 2^k width
109 w[k] w[k] Extents to N=4^k
110 <------------------------>
112 total width = 2^k + 2*w[k]
113 N=4^k is on either the X or Y axis and for the extents here it's taken
115 rotated as necessary to be horizontal. k=2 N=4^2=16 shown above is
116 already horizontal. The next level k=3 N=64=4^3 would be rotated -90
117 degrees to be horizontal.
118 The width w[k] is measured from the N=0 and N=4^k endpoints. It
120 doesn't include the 2^k length between those endpoints. The two ends
121 are symmetric so the extent is the same at each end.
122 h[k] = 2^k - 1 0,1,3,7,15,31,etc
124 w[k] = / 0 for k=0
126 \ 2^(k-1) - 1 for k>=1 0,0,1,3,7,15,etc
127 l[k] = / 0 for k<=1
129 \ 2^(k-2) - 1 for k>=2 0,0,0,1,3,7,etc
130 The initial N=0 to N=64 shown above is k=3. h[3]=7 is the X=-7
132 horizontal. l[3]=1 is the X=1 horizontal. w[3]=3 is the vertical Y=3,
133 and also Y=-11 which is 3 below the endpoint N=64 at Y=8.
134 Expressed as a fraction of the 2^k distance between the endpoints the
136 extents approach total 2 wide by 1.25 high,
137 *------------------* <-+
139 | | | 1
140 *--* *--* | total
141 | | | height
142 * N=4^k N=0 * <-+ -> 1+1/4
143 | | | | | 1/4
144 *--*--* *--*--* <-+
145 ^-----^ ^-----^
147 1/2 1 1/2 total width -> 2
148 The extent formulas can be found by considering the self-similar
150 blocks. The initial k=0 is a single line segment and all its extents
151 are 0.
152 h[0] = 0
154 N=1 ----- N=0
155 l[0] = 0
156 w[0] = 0
157 Thereafter the replication overlap as
159 +-------+---+-------+
161 | | | |
162 +------+ | | +------+
163 | | D | | C | | B | | <-+
164 | +-------+---+-------+ | | 2^(k-1)
165 | | | | | previous
166 | | | | | level ends
167 | E | | A | <-+
168 +------+ +------+
169 ^---------------^
171 2^k this level ends
172 w[k] = max (h[k-1], w[k-1]) # right of A,B
174 h[k] = 2^(k-1) + max (h[k-1], w[k-1]) # above B,C,D
175 l[k] = max w[k-1], l[k-1]-2^(k-1) # below A,E
176 Since h[k]=2^(k-1)+w[k] have h[k] > w[k] for k>=1 and with the initial
178 h[0]=w[k]=0 have h[k]>=w[k] always. So the max of those two is h.
179 h[k] = 2^(k-1) + h[k-1] giving h[k] = 2^k-1 for k>=1
181 w[k] = h[k-1] giving w[k] = 2^(k-1)-1 for k>=1
182 The max for l[k] is always w[k-1] as l[k] is never big enough that the
184 parts B-C and C-D can extend down past their 2^(k-1) vertical position.
185 (l[0]=w[0]=0 and thereafter by induction l[k]<=w[k].)
186 l[k] = w[k-1] giving l[k] = 2^(k-2)-1 for k>=2
188 Repeated Points
190 The curve crosses itself and can repeat X,Y positions up to 4 times.
191 The first double, triple and quadruple points are at
192 visits X,Y N
194 ------ ------- ----------------------
195 2 -2, 3 7, 9
196 3 18, -7 189, 279, 281
197 4 -32, 55 1727, 1813, 2283, 2369
198 Each line segment between integer points is traversed at most 2 times,
200 once forward and once backward. There's 4 lines reaching each integer
201 point and this line traversal means the points are visited at most 4
202 times.
203 As per "Direction" below the direction of the curve is given by the
205 count of 1-bits in N. Since no line is repeated each of the N values
206 at a given X,Y have a different count-1-bits mod 4. For example N=7 is
207 3 1-bits and N=9 is 2 1-bits. The full counts need not be consecutive,
208 as for example N=1727 is 9 1-bits and N=2369 is 4 1-bits.
209 The maximum of 2 line segment traversals can be seen from the way the
211 curve replicates. Suppose the entire plane had all line segments
212 traversed forward and backward.
213 v | v |
215 -- <-------- <-
216 [0,1] [1,1] [X,Y] = integer points
217 -> --------> -- each edge traversed
218 | ^ | ^ forward and backward
219 | | | |
220 | | | |
221 v | v |
222 -- <-------- <--
223 [0,0] [1,0]
224 -> --------> --
225 | ^ | ^
226 Then when each line segment expands on the right the result is the same
228 pattern of traversals -- viewed rotated by 45-degrees and scaled by
229 factor sqrt(2).
230 \ v / v \ v / v
232 [0,1] [1,1]
233 / / ^ \ ^ / ^ \
234 / / \ \ / / \ \
235 \ \ / /
236 \ v / v
237 [1/2,1/2]
238 ^ / ^ \
239 / / \ \
240 \ \ / / \ \ / /
241 \ v / v \ v / v
242 [0,0] 1,0
243 ^ / ^ \ ^ / ^ \
244 The curve is a subset of this pattern. It begins as a single line
246 segment which has this pattern and thereafter the pattern preserves
247 itself. Hence at most 2 segment traversals in the curve.
248 Tiling
250 The segment traversal argument above can also be made by taking the
251 line segments as triangles which are a quarter of a unit square with
252 peak pointing to the right of the traversal direction.
253 to *
255 ^\
256 | \
257 | \ triangle peak
258 | /
259 | /
260 |/ quarter of a unit square
261 from *
262 These triangles in the two directions tile the plane. On expansion
264 each splits into 2 halves in new positions. Those parts don't overlap
265 and the plane is still tiled. See for example
266 Larry Riddle
268 <http://ecademy.agnesscott.edu/~lriddle/ifs/levy/levy.htm>
269 <http://ecademy.agnesscott.edu/~lriddle/ifs/levy/tiling.htm>
270 For the integer version of the curve this kind of tiling can be used to
272 combine copies of the curve so that each every point is visited
273 precisely 4 times. The h[k], w[k] and l[k] extents above are less than
274 the 2^k endpoint length, so a square of side 2^k can be fully tiled
275 with copies of the curve at each corner,
276 | ^ | ^
278 | | | | 24 copies of the curve
279 | | | | to visit all points of the
280 v | v | inside square ABCD
281 <------- <-------- <-------- precisely 4 times each
282 A B
283 --------> --------> --------> each part points
284 | ^ | ^ N=0 to N=4^k-1
285 | | | | rotated and shifted
286 | | | | suitably
287 v | v |
288 <-------- <-------- <--------
289 C D
290 -------- --------> -------->
291 | ^ | ^
292 | | | |
293 | | | |
294 v | v |
295 The four innermost copies of the curve cover most of the inside square,
297 but the other copies surrounding them loop into the square and fill in
298 the remainder to make 4 visits at every point.
299 It's interesting to note that a set of 8 curves at the origin only
301 covers the axes with 4-fold visits,
302 | ^ 8 arms at the origin
304 | | cover only X,Y axes
305 v | with 4-visits
306 <-------- <--------
307 0,0 away from the axes
308 -------- --------> some points < 4 visits
309 | ^
310 | |
311 v |
312 This means that if the path had some sort of "arms" of multiple curves
314 extending from the origin then it would visit all points on the axes
315 X=0 Y=0 a full 4 times, but off the axes there would be points without
316 full 4 visits.
317 See examples/c-curve-wx.pl for a wxWidgets program drawing various
319 forms and tilings of the curve.
320
322 See "FUNCTIONS" in Math::PlanePath for the behaviour common to all path
323 classes.
324 "$path = Math::PlanePath::CCurve->new ()"
326 Create and return a new path object.
327 "($x,$y) = $path->n_to_xy ($n)"
329 Return the X,Y coordinates of point number $n on the path. Points
330 begin at 0 and if "$n < 0" then the return is an empty list.
331 Fractional positions give an X,Y position along a straight line
333 between the integer positions.
334 "$n = $path->xy_to_n ($x,$y)"
336 Return the point number for coordinates "$x,$y". If there's
337 nothing at "$x,$y" then return "undef". If "$x,$y" is visited more
338 than once then return the smallest $n which visits it.
339 "@n_list = $path->xy_to_n_list ($x,$y)"
341 Return a list of N point numbers at coordinates "$x,$y". If
342 there's nothing at "$x,$y" then return an empty list.
343 A given "$x,$y" is visited at most 4 times so the returned list is
345 at most 4 values.
346 "$n = $path->n_start()"
348 Return 0, the first N in the path.
349 Level Methods
351 "($n_lo, $n_hi) = $path->level_to_n_range($level)"
352 Return "(0, 2**$level)".
353
355 Some formulas and results can also be found in the author's
356 mathematical write-up
357 <http://user42.tuxfamily.org/c-curve/index.html>
359 Direction
361 The direction or net turn of the curve is the count of 1 bits in N,
362 direction = count_1_bits(N) * 90degrees
364 For example N=11 is binary 1011 has three 1 bits, so direction 3*90=270
366 degrees, ie. to the south.
367 This bit count is because at each power-of-2 position the curve is a
369 copy of the lower bits but turned +90 degrees, so +90 for each 1-bit.
370 For powers-of-2 N=2,4,8,16, etc, there's only a single 1-bit so the
372 direction is always +90 degrees there, ie. always upwards.
373 Turn
375 At each point N the curve can turn in any direction: left, right,
376 straight, or 180 degrees back. The turn is given by the number of low
377 0-bits of N,
378 turn right = (count_low_0_bits(N) - 1) * 90degrees
380 For example N=8 is binary 0b100 which is 2 low 0-bits for
382 turn=(2-1)*90=90 degrees to the right.
383 When N is odd there's no low zero bits and the turn is always
385 (0-1)*90=-90 to the right, so every second turn is 90 degrees to the
386 left.
387 Next Turn
389 The turn at the point following N, ie. at N+1, can be calculated by
390 counting the low 1-bits of N,
391 next turn right = (count_low_1_bits(N) - 1) * 90degrees
393 For example N=11 is binary 0b1011 which is 2 low one bits for
395 nextturn=(2-1)*90=90 degrees to the right at the following point, ie.
396 at N=12.
397 This works simply because low 1-bits like ..0111 increment to low
399 0-bits ..1000 to become N+1. The low 1-bits at N are thus the low
400 0-bits at N+1.
401 N to dX,dY
403 "n_to_dxdy()" is implemented using the direction described above. For
404 integer N the count mod 4 gives the direction for dX,dY.
405 dir = count_1_bits(N) mod 4
407 dx = dir_to_dx[dir] # table 0 to 3
408 dy = dir_to_dy[dir]
409 For fractional N the direction at int(N)+1 can be obtained from the
411 direction at int(N) and the turn at int(N)+1, which is the low 1-bits
412 of N per "Next Turn" above. Those two directions can then be combined
413 as described in "N to dX,dY -- Fractional" in Math::PlanePath.
414 # apply turn to make direction at Nint+1
416 turn = count_low_1_bits(N) - 1 # N integer part
417 dir = (dir - turn) mod 4 # direction at N+1
418 # adjust dx,dy by fractional amount in this direction
420 dx += Nfrac * (dir_to_dx[dir] - dx)
421 dy += Nfrac * (dir_to_dy[dir] - dy)
422 A small optimization can be made by working the "-1" of the turn
424 formula into a +90 degree rotation of the "dir_to_dx[]" and
425 "dir_to_dy[]" parts by swap and sign change,
426 turn_plus_1 = count_low_1_bits(N) # on N integer part
428 dir = (dir - turn_plus_1) mod 4 # direction-1 at N+1
429 # adjustment including extra +90 degrees on dir
431 dx -= $n*(dir_to_dy[dir] + dx)
432 dy += $n*(dir_to_dx[dir] - dy)
433 X,Y to N
435 The N values at a given X,Y can be found by taking terms low to high
436 from the complex number formula (the same as given above)
437 X+iY = b^k N = 2^k + 2^(k1) + 2^(k2) + ... in binary
439 + b^k1 * i base b=1+i
440 + b^k2 * i^2
441 + ...
442 If the lowest term is b^0 then X+iY has X+Y odd. If the lowest term is
444 not b^0 but instead some power b^n then X+iY has X+Y even. This is
445 because a multiple of b=1+i,
446 X+iY = (x+iy)*(1+i)
448 = (x-y) + (x+y)i
449 so X=x-y Y=x+y
450 sum X+Y = 2x is even if X+iY a multiple of 1+i
451 So the lowest bit of N is found by
453 bit = (X+Y) mod 2
455 If bit=1 then a power i^p is to be subtracted from X+iY. p is how many
457 1-bits are above that point, and this is not yet known. It represents
458 a direction to move X,Y to put it on an even position. It's also the
459 direction of the step N-2^l to N, where 2^l is the lowest 1-bit of N.
460 The reduction should be attempted with p commencing as each of the four
462 possible directions N,S,E,W. Some or all will lead to an N. For
463 quadrupled points (such as X=-32, Y=55 described above) all four will
464 lead to an N.
465 for p 0 to 3
467 dX,dY = i^p # directions [1,0] [0,1] [-1,0] [0,-1]
468 loop until X,Y = [0,0] or [1,0] or [-1,0] or [0,1] or [0,-1]
470 {
471 bit = X+Y mod 2 # bits of N from low to high
472 if bit == 1 {
473 X -= dX # move to "even" X+Y == 0 mod 2
474 Y -= dY
475 (dX,dY) = (dY,-dX) # rotate -90 as for p-1
476 }
477 X,Y = (X+Y)/2, (Y-X)/2 # divide (X+iY)/(1+i)
478 }
479 if not (dX=1 and dY=0)
481 wrong final direction, try next p
482 if X=dX and Y=dY
483 further high 1-bit for N
484 found an N
485 if X=0 and Y=0
486 found an N
487 The "loop until" ends at one of the five points
489 0,1
491 |
492 -1,0 -- 0,0 -- 1,0
493 |
494 0,-1
495 It's not possible to wait for X=0,Y=0 to be reached because some dX,dY
497 directions will step infinitely among the four non-zeros. Only the
498 case X=dX,Y=dY is sure to reach 0,0.
499 The successive p decrements which rotate dX,dY by -90 degrees must end
501 at p == 0 mod 4 for highest term in the X+iY formula having i^0=1.
502 This means must end dX=1,dY=0 East. If this doesn't happen then there
503 is no N for that p direction.
504 The number of 1-bits in N is == p mod 4. So the order the N values are
506 obtained follows the order the p directions are attempted. In general
507 the N values will not be smallest to biggest N so a little sort is
508 necessary if that's desired.
509 It can be seen that sum X+Y is used for the bit calculation and then
511 again in the divide by 1+i. It's convenient to write the whole loop in
512 terms of sum S=X+Y and difference D=Y-X.
513 for dS = +1 or -1 # four directions
515 for dD = +1 or -1 #
516 S = X+Y
517 D = Y-X
518 loop until -1 <= S <= 1 and -1 <= D <= 1 {
520 bit = S mod 2 # bits of N from low to high
521 if bit == 1 {
522 S -= dS # move to "even" S+D == 0 mod 2
523 D -= dD
524 (dS,dD) = (dD,-dS) # rotate -90
525 }
526 (S,D) = (S+D)/2, (D-S)/2 # divide (S+iD)/(1+i)
527 }
528 if not (dS=1 and dD=-1)
530 wrong final direction, try next dS,dD direction
531 if S=dS and D=dD
532 further high 1-bit for N
533 found an N
534 if S=0 and D=0
535 found an N
536 The effect of S=X+Y, D=Y-D is to rotate by -45 degrees and use every
538 second point of the plane.
539 D= 2 X=0,Y=2 . rotate -45
541 D= 1 X=0,Y=1 . X=1,Y=2 .
543 D= 0 X=0,Y=0 . X=1,Y=1 . X=2,Y=2
545 D=-1 X=1,Y=0 . X=2,Y=1 .
547 D=-2 X=2,Y=0 .
549 S=0 S=1 S=2 S=3 S=4
551 The final five points described above are then in a 3x3 block at the
553 origin. The four in-between points S=0,D=1 etc don't occur so range
554 tests -1<=S<=1 and -1<=D<=1 can be used.
555 S=-1,D=1 . S=1,D=1
557 . S=0,D=0 .
559 S=-1,D=-1 . S=1,D=-1
561 Segments by Direction
563 In a level N=0 to N=2^k-1 inclusive, the number of segments in each
564 direction 0=East, 1=North, 2=West, 3=South are given by
565 k=0 for k >= 1
567 --- ----------
568 M0[k] = 1, 2^(k-2) + d(k+2)*2^(h-1)
569 M1[k] = 0, 2^(k-2) + d(k+0)*2^(h-1)
570 M2[k] = 0, 2^(k-2) + d(k-2)*2^(h-1)
571 M3[k] = 0, 2^(k-2) + d(k-4)*2^(h-1)
572 where h = floor(k/2)
574 and d(m) = 0 1 1 1 0 -1 -1 -1
575 for m == 0 to 7 mod 8
576 M0[k] = 1, 1, 1, 1, 2, 6, 16, 36, 72, 136, 256, ...
578 M1[k] = 0, 1, 2, 3, 4, 6, 12, 28, 64, 136, 272, ...
579 M2[k] = 0, 0, 1, 3, 6, 10, 16, 28, 56, 120, 256, ...
580 M3[k] = 0, 0, 0, 1, 4, 10, 20, 36, 64, 120, 240, ...
581 d(n) is a factor +1, -1 or 0 according to n mod 8. Each M goes as a
583 power 2^(k-2), so roughly 1/4 each, but a half power 2^(h-1) possibly
584 added or subtracted in a k mod 8 pattern. In binary this is a 2^(k-2)
585 high 1-bit with another 1-bit in the middle added or subtracted.
586 The total is 2^k since there are a total 2^k points from N=0 to 2^k-1
588 inclusive.
589 M0[k] + M1[k] + M2[k] + M3[k] = 2^k
591 It can be seen that the d(n) parts sum to 0 so the 2^(h-1) parts cancel
593 out leaving 4*2^(k-2) = 2^k.
594 d(0) + d(2) + d(4) + d(6) = 0
596 d(1) + d(3) + d(5) + d(7) = 0
597 The counts can be calculated in two ways. Firstly they satisfy mutual
599 recurrences. Each adds the preceding rotated M.
600 M0[k+1] = M0[k] + M3[k] initially M0[0] = 1 (N=0 to N=1)
602 M1[k+1] = M1[k] + M0[k] M1[0] = 0
603 M2[k+1] = M2[k] + M1[k] M2[0] = 0
604 M3[k+1] = M3[k] + M2[k] M3[0] = 0
605 Geometrically this can be seen from the way each level extends by a
607 copy of the previous level rotated +90,
608 7---6---5 Easts in N=0 to 8
610 | | = Easts in N=0 to 4
611 8 4---3 + Wests in N=0 to 4
612 | since N=4 to N=8 is
613 2 the N=0 to N=4 rotated +90
614 |
615 0---1
616 For the bits in N, level k+1 introduces a new bit either 0 or 1. In
618 M0[k+1] the a 0-bit is count M0[k] the same direction, and when a 1-bit
619 is M3[k] since one less bit mod 4. Similarly the other counts.
620 Some substitutions give 3rd order recurrences
622 for k >= 4
624 M0[k] = 4*M0[k-1] - 6*M0[k-2] + 4*M0[k-3] initial 1,1,1,1
625 M1[k] = 4*M1[k-1] - 6*M1[k-2] + 4*M1[k-3] initial 0,1,2,3
626 M2[k] = 4*M2[k-1] - 6*M2[k-2] + 4*M2[k-3] initial 0,0,1,3
627 M3[k] = 4*M3[k-1] - 6*M3[k-2] + 4*M3[k-3] initial 0,0,0,1
628 The characteristic polynomial of these recurrences is
630 x^3 - 4x^2 + 6x - 4
632 = (x-2) * (x - (1-i)) * (x - (1+i))
633 So explicit formulas can be written in powers of the roots 2, 1-i and
635 1+i,
636 M0[k] = ( 2^k + (1-i)^k + (1+i)^k )/4 for k>=1
638 M1[k] = ( 2^k + i*(1-i)^k - i*(1+i)^k )/4
639 M2[k] = ( 2^k - (1-i)^k - (1+i)^k )/4
640 M3[k] = ( 2^k - i*(1-i)^k + i*(1+i)^k )/4
641 The complex numbers 1-i and 1+i are 45 degree lines clockwise and anti-
643 clockwise respectively. The powers turn them in opposite directions so
644 the imaginary parts always cancel out. The remaining real parts can be
645 had by a half power h=floor(k/2) which is the magnitude
646 abs(1-i)=sqrt(2) projected onto the real axis. The sign selector d(n)
647 above is whether the positive or negative part of the real axis, or
648 zero when at the origin.
649 The second way to calculate is the combinatorial interpretation that
651 per "Direction" above the direction is count_1_bits(N) mod 4 so East
652 segments are all N values with count_1_bits(N) == 0 mod 4, ie. N with
653 0, 4, 8, etc many 1-bits. The number of ways to have those bit counts
654 within total k bits is k choose 0, 4, 8 etc.
655 M0[k] = /k\ + /k\ + ... + / k\ m = floor(k/4)
657 \0/ \4/ \4m/
658 M1[k] = /k\ + /k\ + ... + / k \ m = floor((k-1)/4)
660 \1/ \5/ \4m+1/
661 M2[k] = /k\ + /k\ + ... + / k \ m = floor((k-2)/4)
663 \2/ \6/ \4m+2/
664 M3[k] = /k\ + /k\ + ... + / k \ m = floor((k-3)/4)
666 \3/ \7/ \4m+3/
667 The power forms above are cases of the identity by Ramus for sums of
669 binomial coefficients in arithmetic progression like this. (See Knuth
670 volume 1 section 1.2.6 exercise 30 for a form with cosines resulting
671 from w=i+1 as 8th roots of unity.)
672 The total M0+M1+M2+M3=2^k is the total binomials across a row of
674 Pascal's triangle.
675 /k\ + /k\ + ... + /k\ = 2^k
677 \0/ \1/ \k/
678 It's interesting to note the M counts here are the same in the dragon
680 curve (Math::PlanePath::DragonCurve). The shapes of the curves are
681 different since the segments are in a different order, but the total
682 puts points N=2^k at the same X,Y position.
683 Right Boundary
685 The length of the right-side boundary of the curve, which is the
686 outside of the "C", from N=0 to N=2^k is
687 R[k] = / 7*2^h - 2k - 6 if k even
689 \ 10*2^h - 2k - 6 if k odd
690 where h = floor(k/2)
691 = 1, 2, 4, 8, 14, 24, 38, 60, 90, 136, 198, 292, 418, ...
692 R[k] = (7/2 + 5/2 * sqrt(2)) * ( sqrt(2))^k
694 + (7/2 - 5/2 * sqrt(2)) * (-sqrt(2))^k
695 - 2*k - 6
696 R[k] = 2*R[k-1] + R[k-2] - 4*R[k-3] + 2*R[k-4]
698 The length doubles until R[4]=14 which is points N=0 to N=2^4=16. At
700 k=4 the points N=7,8,9 have turned inward and closed off some of the
701 outside of the curve so the boundary less than 2x.
702 11--10--9,7--6--5 right boundary
704 | | | around "outside"
705 13--12 8 4--3 N=0 to N=2^4=16
706 | |
707 14 2 R[4]=14
708 | |
709 15--16 0--1
710 The floor(k/2) and odd/even cases are eliminated by the +/-sqrt(2)
712 powering shown. Those powers are also per the characteristic equation
713 of the recurrence,
714 x^4 - 2*X^3 - x^2 + 4*x - 2
716 = (x - 1)^2 * (x + sqrt(2)) * (x - sqrt(2))
717 roots 1, sqrt(2), -sqrt(2)
718 The right boundary comprises runs of straight lines and zig-zags. When
720 it expands the straight lines become zig-zags and the zig-zags become
721 straight lines. The straight lines all point "forward", which is anti-
722 clockwise.
723 c * a
725 / ^ / ^ / ^
726 => v \ v \ v \
727 D<----C<----B<----A D C B A
728 | ^ / ^
729 v | v \
730 straight S=3 zig-zag Z[k+1] = 2S[k]-2 = 4
731 The count Z here is both sides of each "V" shape from points "a"
733 through to "c". So Z counts the boundary length (rather than the
734 number of "V"s). Each S becomes an upward peak. The first and last
735 side of those peaks become part of the following "straight" section (at
736 A and D), hence Z[k+1]=2*S[k]-2.
737 The zigzags all point "forward" too. When they expand they close off
739 the V shape and become 2 straight lines for each V, which means 1
740 straight line for each Z side. The segment immediately before and
741 after contribute a segment to the resulting straight run too, hence
742 S[k+1]=Z[k]+2.
743 C B A *<---C<---*<---B<---*<---A<---*
745 / ^ / ^ / ^ | | | |
746 v \ v \ v \ => | | | |
747 * * * * <--* * * *<--
748 / ^
749 v \
750 zig-zag Z=4 segments straight S[k+1] = Z[k]+2 = 6
751 The initial N=0 to N=1 is a single straight segment S[0]=1 and from
753 there the runs grow. N=1 to N=3 is a straight section S[1]=2. Z[0]=0
754 represents an empty zigzag at N=1. Z[1] is the first non-empty at N=3
755 to N=5.
756 h S[h] Z[h] Z[h] = 2*S[h]-2
758 -- ---- ---- S[h+1] = Z[h]+2
759 0 1 0
760 1 2 2 S[h+1] = 2*S[h]-2+2 = 2*S[h]
761 2 4 6 so
762 3 8 14 S[h] = 2^h
763 4 16 30 Z[h] = 2*2^h-2
764 5 32 62
765 5 64 126
766 The curve N=0 to N=2^k is symmetric at each end and is made up of runs
768 S[0], Z[0], S[1], Z[1], etc, of straight and zigzag alternately at each
769 end. When k is even there's a single copy of a middle S[k/2]. When k
770 is odd there's a single middle Z[(k-1)/2] (with an S[(k-1)/2] before
771 and after). So
772 / i=h-1 \ # where h = floor(k/2)
774 R[k] = 2 * | sum S[i]+Z[i] |
775 \ i=0 /
776 + S[h]
777 + / S[h]+Z[h] if k odd
778 \ 0 if k even
779 = 2*( 1+2+4+...+2^(h-1) # S[0] to S[h-1]
781 + 2+4+8+...+2^h - 2*h) # Z[0] to Z[h-1]
782 + 2^h # S[h]
783 + if k odd (2^h + 2*2^h - 2) # possible S[h]+Z[h]
784 = 2*(2^h-1 + 2*2^h-2 - 2h) + 2^h + (k odd 3*2^h - 2)
786 = 7*2^h - 4h-6 + (if k odd then + 3*2^h - 2)
787 = 7*2^h - 2k-6 + (if k odd then + 3*2^h)
788 Convex Hull Boundary
790 A convex hull is the smallest convex polygon which contains a given set
791 of points. For the C curve the boundary length of the convex hull for
792 points N=0 to N=2^k inclusive is
793 hull boundary[k]
795 / 2 if k=0
796 | 2+sqrt(2) if k=1
797 = | 6 if k=2
798 | 6*2^(h-1) + (7*2^(h-1) - 4)*sqrt(2) if k odd >=3
799 \ 7*2^(h-1) + (3*2^(h-1) - 4)*sqrt(2) if k even >=4
800 where h = floor(k/2)
801 k hull boundary
803 --- ----------------------------
804 0 2 + 0 * sqrt(2) = 2
805 1 2 + 1 * sqrt(2) = 3.41
806 2 6 + 0 * sqrt(2) = 6
807 3 6 + 3 * sqrt(2) = 10.24
808 4 14 + 2 * sqrt(2) = 16.82
809 5 12 + 10 * sqrt(2) = 26.14
810 6 28 + 8 * sqrt(2) = 39.31
811 7 24 + 24 * sqrt(2) = 57.94
812 8 56 + 20 * sqrt(2) = 84.28
813 9 48 + 52 * sqrt(2) = 121.53
814 The integer part is the straight sides of the hull and the sqrt(2) part
816 is the diagonal sides of the hull.
817 When k is even the hull has the following shape. The sides are as per
819 the right boundary above but after Z[h-2] the curl goes inwards and so
820 parts beyond Z[h-2] are not part of the hull. Each Z stair-step
821 diagonal becomes a sqrt(2) length for the hull. Z counts both vertical
822 and horizontal of each stairstep, hence sqrt(2)*Z/2 for the hull
823 boundary.
824 S[h]
826 *--------* * Z=2
827 Z[h-1] / \ Z[h-1] | \ diagonal
828 / \ | \ sqrt(2)*Z/2
829 * * *----* = sqrt(2)
830 S[h-1] | | S[h-1]
831 | |
832 * *
833 Z[h-2] \ / Z[h-2]
834 *-- --*
835 S[h] + Z[h-2]-Z[h-1]
837 k even
839 hull boundary[k] = S[h] + 2*S[h-1] + S[h+Z[h-2]-Z[h-1]
840 + sqrt(2)*(2*Z[h-1] + 2*Z[h-2])/2
841 When k is odd the shape is similar but Z[h] in the middle.
843 S[h]
845 *----*
846 Z[h-1] / \ middle
847 * \ Z[h]
848 S[h-1] | \
849 * *
850 \ | S[h]
851 Z[h] |
852 + 2*(S[h]-S[h-1]) *
853 \ / Z[h-1]
854 *---*
855 S[h-1]
856 k odd
858 hull boundary[k] = 2*S[h] + 2*S[h-1]
859 + sqrt(2)*(Z[h]/2 + 2*Z[h-1]/2
860 + Z[h]/2 + S[h]-S[h-1]
861 Convex Hull Area
863 The area of the convex hull for points N=0 to N=2^k inclusive is
864 / 0 if k=0
866 | 1/2 if k=1
867 HA[k] = | 2 if k=2
868 | 35*2^(k-4) - 13*2^(h-1) + 2 if k odd >=3
869 \ 35*2^(k-4) - 10*2^(h-1) + 2 if k even >=4
870 where h = floor(k/2)
871 = 0, 1/2, 2, 13/2, 17, 46, 102, 230, 482, 1018, 2082, 4274, ...
873 HA[1] and HA[3] are fractions but all others are integers.
875 The area can be calculated from the shapes shown for the hull boundary
877 above. For k odd it can be noted the width and height are equal, then
878 the various corners are cut off.
879 Line Points
881 The number of points which fall on straight and diagonal lines from the
882 endpoints can be calculated by considering how the previous level
883 duplicates to make the next.
884 d d
886 c \ / c
887 b | + | b
888 \ | / \ | / curve endpoints
889 \ | / \ | / "S" start
890 \|/ \|/ "E" end
891 a------E---e---S------a
892 /|\ /|\
893 / | \ / | \
894 / | f f | \
895 h g g h
896 The curve is rotated to make the endpoints horizontal. Each "a"
898 through "h" is the number of points which fall on the respective line.
899 The curve is symmetric in left to right so the line counts are the same
900 each side in mirror image.
901 "S" start and "E" end points are not included in any of the counts.
903 "e" is the count in between S and E. The two "d" lines meet at point
904 "+" and that point is counted in d. That point is where two previous
905 level curves meet for k>=1. Points are visited up to 4 times (per
906 "Repeated Points" above) and all those multiple visits are counted.
907 The following diagram shows how curve level k+1 is made from two level
909 k curves. One is from S to M and another M to E.
910 |\ /| curve level k copies
912 | \ / | S to M and M to E
913 | c+a c+a | making curve k+1 S to E
914 | \|/ |
915 \ | --M-- | /
916 \ | /|\ | c a[k+1] = b[k]
917 c d e+g e+g d / b[k+1] = c[k]
918 \ | / \ | / c[k+1] = d[k]
919 \|/ \|/ d[k+1] = a[k]+c[k] + e[k]+g[k] + 1
920 b-------E--f---f--S-------b e[k+1] = 2*f[k]
921 /|\ /|\ f[k+1] = g[k]
922 a | g g | a g[k+1] = h[k]
923 / h \ / h \ h[k+1] = a[k]
924 / | \ / | \
925 / | | \
926 For example the line S to M is an e[k], but also the M to E contributes
928 a g[k] on that same line so e+g. Similarly c[k] and a[k] on the outer
929 sides of M. Point M itself is visited too so the grand total for
930 d[k+1] is a+c+e+g+1. The other lines are simpler, being just rotations
931 except for the middle line e[k+1] which is made of two f[k].
932 The successive g[k+1]=h[k]=a[k-1]=b[k-2]=c[k-3]=d[k-4] can be
934 substituted into the d to give a recurrence
935 d[k+1] = d[k-1] + d[k-3] + d[k-5] + 2*d[k-7] + 1
937 = 0,1,1,2,2,4,4,8,8,17,17,34,34,68,68,136,136,273,273,...
938 x + x^2 (common factor 1+x
940 generating function ------------------- in numerator and
941 (1-2*x^2) * (1-x^8) denominator)
942 d[2h-1] = d[2h] = floor( 8/15 * 2^h )
944 The recurrence begins with the single segment N=0 to N=1 and the two
946 endpoints are not included so initial all zeros a[0]=...=h[0]=0.
947 As an example, the N=0 to N=64 picture above is level k=6 and its "d"
949 line relative to those endpoints is the South-West diagonal down from
950 N=0. The points on that line are N=32,30,40,42 giving d[6]=4.
951 All the measures are relative to the endpoint direction. The points on
953 the fixed X or Y axis or diagonal can be found by taking the
954 appropriate a through h, or sum of two of them for both positive and
955 negative of a direction.
956
958 Entries in Sloane's Online Encyclopedia of Integer Sequences related to
959 this path include
960 <http://oeis.org/A179868> (etc)
962 A010059 abs(dX), count1bits(N)+1 mod 2
964 A010060 abs(dY), count1bits(N) mod 2, being Thue-Morse
965 A000120 direction, being total turn, count 1-bits
967 A179868 direction 0to3, count 1-bits mod 4
968 A035263 turn 0=straight or 180, 1=left or right,
970 being (count low 0-bits + 1) mod 2
971 A096268 next turn 1=straight or 180, 0=left or right,
972 being count low 1-bits mod 2
973 A007814 turn-1 to the right,
974 being count low 0-bits
975 A003159 N positions of left or right turn, ends even num 0 bits
977 A036554 N positions of straight or 180 turn, ends odd num 0 bits
978 A146559 X at N=2^k, being Re((i+1)^k)
980 A009545 Y at N=2^k, being Im((i+1)^k)
981 A131064 right boundary length to odd power N=2^(2k-1),
983 being 5*2^n-4n-4, skip initial 1
984 A027383 right boundary length differences
985 A038503 number of East segments in N=0 to N=2^k-1
987 A038504 number of North segments in N=0 to N=2^k-1
988 A038505 number of West segments in N=0 to N=2^k-1
989 A000749 number of South segments in N=0 to N=2^k-1
990 A191689 fractal dimension of the interior boundary
992 A191689 is the fractal dimension which roughly speaking means what
994 power r^k the boundary length grows by when each segment is taken as a
995 little triangle (or similar). There are various holes inside the
996 curling spiralling curve and they are all boundary for this purpose.
997 P. Duvall and J. Keesling, "The Dimension of the Boundary of the
999 Lévy Dragon", International Journal Math and Math Sci, volume 20,
1000 number 4, 1997, pages 627-632. (Preprint "The Hausdorff Dimension
1001 of the Boundary of the Lévy Dragon"
1002 <http://at.yorku.ca/p/a/a/h/08.htm>.)
1003
1005 Math::PlanePath, Math::PlanePath::DragonCurve,
1006 Math::PlanePath::AlternatePaper, Math::PlanePath::KochCurve
1007 ccurve(6x) back-end for xscreensaver(1) which displays the C curve (and
1009 various other dragon curve and Koch curves).
1010
1012 <http://user42.tuxfamily.org/math-planepath/index.html>
1013
1015 Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017 Kevin Ryde
1016 This file is part of Math-PlanePath.
1018 Math-PlanePath is free software; you can redistribute it and/or modify
1020 it under the terms of the GNU General Public License as published by
1021 the Free Software Foundation; either version 3, or (at your option) any
1022 later version.
1023 Math-PlanePath is distributed in the hope that it will be useful, but
1025 WITHOUT ANY WARRANTY; without even the implied warranty of
1026 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
1027 General Public License for more details.
1028 You should have received a copy of the GNU General Public License along
1030 with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
1031perl v5.28.0 2017-12-03 Math::PlanePath::CCurve(3)