1Math::PlanePath::CCurveU(s3e)r Contributed Perl DocumentaMtaitohn::PlanePath::CCurve(3)
2
3
4

NAME

6       Math::PlanePath::CCurve -- Levy C curve
7

SYNOPSIS

9        use Math::PlanePath::CCurve;
10        my $path = Math::PlanePath::CCurve->new;
11        my ($x, $y) = $path->n_to_xy (123);
12

DESCRIPTION

14       This is an integer version of the "C" curve by Lévy.
15
16           "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
20           <http://gallica.bnf.fr/ark:/12148/bpt6k57344323/f53.image>
21           <http://gallica.bnf.fr/ark:/12148/bpt6k57344820>
22
23       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
27                                 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
57                                                              ^
58            -7     -6     -5     -4     -3     -2     -1     X=0     1
59
60       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
65                                         4----3
66                                              |      N=0to2
67                              2               2      repeated
68                              |               |      as N=2to4
69           0----1        0----1          0----1      with turn +90
70
71       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
75       The X,Y position can be written in complex numbers as a recurrence
76
77           with N = 2^k + r      high bit 2^k, rest r<2^k
78
79           C(N) = C(2^k)  + i*C(r)
80                = (1+i)^k + i*C(r)
81
82       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
86           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
92       Notice the i power is not the bit position k, but rather how many
93       1-bits are above the position.
94
95   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
99              *------------------*       <-+
100              |                  |         |
101           *--*                  *--*      | height h[k]
102           |                        |      |
103           *   N=4^k         N=0    *    <-+
104           |     |            |     |      | below l[k]
105           *--*--*            *--*--*    <-+
106
107           ^-----^            ^-----^
108            width     2^k      width
109             w[k]               w[k]           Extents to N=4^k
110
111           <------------------------>
112           total width = 2^k + 2*w[k]
113
114       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
119       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
123           h[k] = 2^k - 1                     0,1,3,7,15,31,etc
124
125           w[k] = /  0            for k=0
126                  \  2^(k-1) - 1  for k>=1    0,0,1,3,7,15,etc
127
128           l[k] = /  0            for k<=1
129                  \  2^(k-2) - 1  for k>=2    0,0,0,1,3,7,etc
130
131       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
135       Expressed as a fraction of the 2^k distance between the endpoints the
136       extents approach total 2 wide by 1.25 high,
137
138              *------------------*       <-+
139              |                  |         |  1
140           *--*                  *--*      |         total
141           |                        |      |         height
142           *   N=4^k         N=0    *    <-+         -> 1+1/4
143           |     |            |     |      |  1/4
144           *--*--*            *--*--*    <-+
145
146           ^-----^            ^-----^
147             1/2        1       1/2     total width -> 2
148
149       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
153                                 h[0] = 0
154                 N=1 ----- N=0
155                                 l[0] = 0
156                           w[0] = 0
157
158       Thereafter the replication overlap as
159
160              +-------+---+-------+
161              |       |   |       |
162           +------+   |   |   +------+
163           |  | D |   | C |   | B |  |        <-+
164           |  +-------+---+-------+  |          | 2^(k-1)
165           |      |           |      |          | previous
166           |      |           |      |          | level ends
167           |    E |           | A    |        <-+
168           +------+           +------+
169
170                ^---------------^
171               2^k this level ends
172
173           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
177       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
180           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
183       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
187           l[k] = w[k-1]   giving l[k] = 2^(k-2)-1 for k>=2
188
189   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
193           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
199       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
204       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
210       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
214             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
227       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
231            \ 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
245       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
249   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
254              to  *
255                  ^\
256                  | \
257                  |  \ triangle peak
258                  |  /
259                  | /
260                  |/       quarter of a unit square
261             from *
262
263       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
267           Larry Riddle
268           <http://ecademy.agnesscott.edu/~lriddle/ifs/levy/levy.htm>
269           <http://ecademy.agnesscott.edu/~lriddle/ifs/levy/tiling.htm>
270
271       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
277                    | ^         | ^
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
296       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
300       It's interesting to note that a set of 8 curves at the origin only
301       covers the axes with 4-fold visits,
302
303                    | ^              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
313       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
318       See examples/c-curve-wx.pl for a wxWidgets program drawing various
319       forms and tilings of the curve.
320

FUNCTIONS

322       See "FUNCTIONS" in Math::PlanePath for the behaviour common to all path
323       classes.
324
325       "$path = Math::PlanePath::CCurve->new ()"
326           Create and return a new path object.
327
328       "($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
332           Fractional positions give an X,Y position along a straight line
333           between the integer positions.
334
335       "$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
340       "@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
344           A given "$x,$y" is visited at most 4 times so the returned list is
345           at most 4 values.
346
347       "$n = $path->n_start()"
348           Return 0, the first N in the path.
349
350   Level Methods
351       "($n_lo, $n_hi) = $path->level_to_n_range($level)"
352           Return "(0, 2**$level)".
353

FORMULAS

355       Some formulas and results can also be found in the author's
356       mathematical write-up
357
358           <http://user42.tuxfamily.org/c-curve/index.html>
359
360   Direction
361       The direction or net turn of the curve is the count of 1 bits in N,
362
363           direction = count_1_bits(N) * 90degrees
364
365       For example N=11 is binary 1011 has three 1 bits, so direction 3*90=270
366       degrees, ie. to the south.
367
368       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
371       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
374   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
379           turn right = (count_low_0_bits(N) - 1) * 90degrees
380
381       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
384       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
388   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
392           next turn right = (count_low_1_bits(N) - 1) * 90degrees
393
394       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
398       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
402   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
406           dir = count_1_bits(N) mod 4
407           dx = dir_to_dx[dir]    # table 0 to 3
408           dy = dir_to_dy[dir]
409
410       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
415           # 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
419           # 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
423       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
427           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
430           # adjustment including extra +90 degrees on dir
431           dx -= $n*(dir_to_dy[dir] + dx)
432           dy += $n*(dir_to_dx[dir] - dy)
433
434   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
438           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
443       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
447           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
452       So the lowest bit of N is found by
453
454           bit = (X+Y) mod 2
455
456       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
461       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
466           for p 0 to 3
467             dX,dY = i^p   # directions [1,0]  [0,1]  [-1,0]  [0,-1]
468
469             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
480             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
488       The "loop until" ends at one of the five points
489
490                   0,1
491                    |
492           -1,0 -- 0,0 -- 1,0
493                    |
494                   0,-1
495
496       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
500       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
505       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
510       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
514           for dS = +1 or -1      # four directions
515             for dD = +1 or -1    #
516               S = X+Y
517               D = Y-X
518
519               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
529               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
537       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
540           D= 2                      X=0,Y=2       .              rotate -45
541
542           D= 1            X=0,Y=1      .       X=1,Y=2       .
543
544           D= 0  X=0,Y=0      .      X=1,Y=1       .       X=2,Y=2
545
546           D=-1            X=1,Y=0      .       X=2,Y=1       .
547
548           D=-2                      X=2,Y=0       .
549
550                  S=0        S=1       S=2        S=3        S=4
551
552       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
556            S=-1,D=1      .      S=1,D=1
557
558               .       S=0,D=0      .
559
560            S=-1,D=-1     .      S=1,D=-1
561
562   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
566                  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
573           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
577           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
582       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
587       The total is 2^k since there are a total 2^k points from N=0 to 2^k-1
588       inclusive.
589
590           M0[k] + M1[k] + M2[k] + M3[k] = 2^k
591
592       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
595           d(0) + d(2) + d(4) + d(6) = 0
596           d(1) + d(3) + d(5) + d(7) = 0
597
598       The counts can be calculated in two ways.  Firstly they satisfy mutual
599       recurrences.  Each adds the preceding rotated M.
600
601           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
606       Geometrically this can be seen from the way each level extends by a
607       copy of the previous level rotated +90,
608
609           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
617       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
621       Some substitutions give 3rd order recurrences
622
623           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
629       The characteristic polynomial  of these recurrences is
630
631           x^3 - 4x^2 + 6x - 4
632           = (x-2) * (x - (1-i)) * (x - (1+i))
633
634       So explicit formulas can be written in powers of the roots 2, 1-i and
635       1+i,
636
637           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
642       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
650       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
656           M0[k] = /k\ + /k\ + ... + / k\      m = floor(k/4)
657                   \0/   \4/         \4m/
658
659           M1[k] = /k\ + /k\ + ... + / k  \    m = floor((k-1)/4)
660                   \1/   \5/         \4m+1/
661
662           M2[k] = /k\ + /k\ + ... + / k  \    m = floor((k-2)/4)
663                   \2/   \6/         \4m+2/
664
665           M3[k] = /k\ + /k\ + ... + / k  \    m = floor((k-3)/4)
666                   \3/   \7/         \4m+3/
667
668       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
673       The total M0+M1+M2+M3=2^k is the total binomials across a row of
674       Pascal's triangle.
675
676           /k\ + /k\ + ... + /k\ = 2^k
677           \0/   \1/         \k/
678
679       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
684   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
688           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
693           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
697           R[k] = 2*R[k-1] + R[k-2] - 4*R[k-3] + 2*R[k-4]
698
699       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
703               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
711       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
715           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
719       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
724                                             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
732       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
738       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
744                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
752       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
757            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
767       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
773                       / 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
780                =  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
785                = 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
789   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
794           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
802             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
815       The integer part is the straight sides of the hull and the sqrt(2) part
816       is the diagonal sides of the hull.
817
818       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
825                         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
836                       S[h] + Z[h-2]-Z[h-1]
837
838           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
842       When k is odd the shape is similar but Z[h] in the middle.
843
844                               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
857           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
862   Convex Hull Area
863       The area of the convex hull for points N=0 to N=2^k inclusive is
864
865                   / 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
872           = 0, 1/2, 2, 13/2, 17, 46, 102, 230, 482, 1018, 2082, 4274, ...
873
874       HA[1] and HA[3] are fractions but all others are integers.
875
876       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
880   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
885                     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
897       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
902       "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
908       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
911                   |\       /|            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
927       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
933       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
936           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
939                                      x + x^2          (common factor 1+x
940           generating function  -------------------     in numerator and
941                                (1-2*x^2) * (1-x^8)     denominator)
942
943           d[2h-1] = d[2h] = floor( 8/15 * 2^h )
944
945       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
948       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
952       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

OEIS

958       Entries in Sloane's Online Encyclopedia of Integer Sequences related to
959       this path include
960
961           <http://oeis.org/A179868> (etc)
962
963           A010059   abs(dX), count1bits(N)+1 mod 2
964           A010060   abs(dY), count1bits(N) mod 2, being Thue-Morse
965
966           A000120   direction, being total turn, count 1-bits
967           A179868   direction 0to3, count 1-bits mod 4
968
969           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
976           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
979           A146559   X at N=2^k, being Re((i+1)^k)
980           A009545   Y at N=2^k, being Im((i+1)^k)
981
982           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
986           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
991           A191689   fractal dimension of the interior boundary
992
993       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
998           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

SEE ALSO

1005       Math::PlanePath, Math::PlanePath::DragonCurve,
1006       Math::PlanePath::AlternatePaper, Math::PlanePath::KochCurve
1007
1008       ccurve(6x) back-end for xscreensaver(1) which displays the C curve (and
1009       various other dragon curve and Koch curves).
1010

HOME PAGE

1012       <http://user42.tuxfamily.org/math-planepath/index.html>
1013

LICENSE

1015       Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 Kevin Ryde
1016
1017       This file is part of Math-PlanePath.
1018
1019       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
1024       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
1029       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/>.
1031
1032
1033
1034perl v5.32.0                      2020-07-28        Math::PlanePath::CCurve(3)
Impressum