1Math::PlanePath::DragonUCsuerrveC(o3n)tributed Perl DocuMmaetnht:a:tPiloannePath::DragonCurve(3)
2
6 Math::PlanePath::DragonCurve -- dragon curve
7
9 use Math::PlanePath::DragonCurve;
10 my $path = Math::PlanePath::DragonCurve->new;
11 my ($x, $y) = $path->n_to_xy (123);
12
14 This is the dragon or paper folding curve by Heighway, Harter, et al.
15 9----8 5---4 2
17 | | | |
18 10--11,7---6 3---2 1
19 | |
20 17---16 13---12 0---1 <- Y=0
21 | | |
22 18-19,15-14,22-23 -1
23 | | |
24 20--21,25-24 -2
25 |
26 26---27 -3
27 |
28 --32 29---28 -4
29 | |
30 31---30 -5
31 ^ ^ ^ ^ ^ ^ ^
33 -5 -4 -3 -2 -1 X=0 1 ...
34 The curve visits "inside" X,Y points twice. The first of these is
36 X=-2,Y=1 which is N=7 and also N=11. The segments N=6,7,8 and
37 N=10,11,12 have touched, but the path doesn't cross itself. The
38 doubled vertices are all like this, touching but not crossing and no
39 edges repeating.
40 Arms
42 The curve fills a quarter of the plane and four copies mesh together
43 perfectly when rotated by 90, 180 and 270 degrees. The "arms"
44 parameter can choose 1 to 4 curve arms successively advancing.
45 For example arms=4 begins as follows, with N=0,4,8,12,etc being the
47 first arm, N=1,5,9,13 the second, N=2,6,10,14 the third and N=3,7,11,15
48 the fourth.
49 arms => 4
51 20 ------ 16
53 |
54 9 ------5/12 ----- 8 23
55 | | | |
56 17 --- 13/6 --- 0/1/2/3 --- 4/15 --- 19
57 | | | |
58 21 10 ----- 14/7 ----- 11
59 |
60 18 ------ 22
61 With four arms every X,Y point is visited twice (except the origin 0,0
63 where all four begin) and every edge between the points is traversed
64 once.
65 Level Angle
67 The first step N=1 is to the right along the X axis and the path then
68 slowly spirals anti-clockwise and progressively fatter. The end of
69 each replication is N=2^level which is at level*45 degrees around,
70 N X,Y angle radial dist
72 ---- ----- ----- -----------
73 1 1,0 0 1
74 2 1,1 45 sqrt(2)
75 4 0,2 90 sqrt(4)=2
76 8 -2,2 135 sqrt(8)
77 16 -4,0 180 sqrt(16)=4
78 32 -4,-4 225 sqrt(32)
79 ...
80 Here's points N=0 to N=2^9=512. "0" is the origin and "+" is N=512.
82 Notice it's spiralled around full-circle to angle 45 degrees up again,
83 like the initial N=2.
84 * * * *
86 * * * * * *
87 * * * * * * * * *
88 * * * * * * * * *
89 * * * * * * * *
90 * * * * * * * + * *
91 * * * * * * * *
92 * * * * * * *
93 * * * * * * * *
94 * * * * * *
95 * * * *
96 * * * * * * *
97 * * * * * * * * * *
98 * * * * * * * * * * *
99 * * * * * * * * * *
100 * * * * * * * * * * * * * * *
101 * * * * * * * * * * * * * * * *
102 * * * * * * * * * * * * * *
103 * * * * * * * * * * * *
104 * * * * * * * * * * * * * * *
105 * * * * * * * * * * * * * * * *
106 * * * * 0 * * * * * * * * * * * * * * *
107 * * * * * * * * * * * * * * *
108 * * * * * * * * * * * * * *
109 * * * * * * * * * * * * * *
110 * * * * * * * * * * * * * * * * *
111 * * * * * * * * * * * * * * * * *
112 * * * * * * * * * * * * * * * * *
113 * * * * * * * * * *
114 * * * * * * * * * * * * * *
115 * * * * * * * * * *
116 * * * * * * * *
117 At a power of two Nlevel=2^level for N=2 or higher, the curve always
119 goes upward from Nlevel-1 to Nlevel, and then goes to the left for
120 Nlevel+1. For example at N=16 the curve goes up N=15 to N=16, then
121 left for N=16 to N=17. Likewise at N=32, etc. The spiral curls around
122 ever further but the self-similar twist back means the Nlevel endpoint
123 is always at this same up/left orientation. See "Total Turn" below for
124 the net direction in general.
125 Level Ranges
127 The X,Y extents of the path through to Nlevel=2^level can be expressed
128 as a "length" in the direction of the Xlevel,Ylevel endpoint and a
129 "width" across.
130 level even, so endpoint is a straight line
132 k = level/2
133 +--+ <- Lmax
135 | |
136 | E <- Lend = 2^k at Nlevel=2^level
137 |
138 +-----+
139 |
140 O | <- Lstart=0
141 | |
142 +--+ <- Lmin
143 ^ ^
145 Wmin Wmax
146 Lmax = (7*2^k - 4)/6 if k even
148 (7*2^k - 2)/6 if k odd
149 Lmin = - (2^k - 1)/3 if k even
151 - (2^k - 2)/3 if k odd
152 Wmax = (2*2^k - 2) / 3 if k even
154 (2*2^k - 1) / 3 if k odd
155 Wmin = Lmin
157 For example level=2 is to Nlevel=2^2=4 and k=level/2=1 is odd so it
159 measures as follows,
160 4 <- Lmax = (7*2^1 - 2)/6 = 2
162 |
163 3--2
164 |
165 0--1 <- Lmin = -(2^1 - 2)/3 = 0
166 ^ ^Wmax = (2*2^1 - 1)/3 = 1
168 |
169 Wmin = Lmin = 0
170 Or level=4 is to Nlevel=2^4=16 and k=4/2=2 is even. It measures as
172 follows. The lengthways "L" measures are in the direction of the N=16
173 endpoint and the "W" measures are across.
174 9----8 5---4 <- Wmax = (2*2^2 - 2)/3 = 2
176 | | | |
177 10--11,7---6 3---2
178 | |
179 16 13---12 0---1
180 | |
181 15---14 <- Wmin = -(2^2 - 1)/3 = -1
182 ^ ^Lmin = Wmin = -1
184 |
185 Lmax = (7*2^2 - 4)/6 = 4
186 The formulas are all integer values, but the fractions 7/6, 1/3 and 2/3
188 show the limits as the level increases. If scaled so that length
189 Lend=2^k is reckoned as 1 unit then Lmax extends 1/6 past the end, Lmin
190 and Wmin extend 1/3, and Wmax extends across 2/3.
191 +--------+ --
193 | - | 1/6 total length
194 || | | = 1/6+1+1/3 = 3/2
195 || E | --
196 || |
197 || |
198 | \ | 1
199 | \ |
200 | --\ |
201 | \ |
202 | ||
203 | O || --
204 | | ||
205 | | || 1/3
206 | ---- |
207 +--------+ --
208 1/3| 2/3
209 total width = 1/3+2/3 = 1
211 Paper Folding
213 The path is called a paper folding curve because it can be generated by
214 thinking of a long strip of paper folded in half repeatedly and then
215 unfolded so each crease is a 90 degree angle. The effect is that the
216 curve repeats in successive doublings turned by 90 degrees and
217 reversed.
218 The first segment unfolds, pivoting at the "1",
220 2
222 -> |
223 unfold / |
224 ===> | |
225 |
226 0-------1 0-------1
227 Then the same again with that L shape, pivoting at the "2", then after
229 that pivoting at the "4", and so on.
230 4
232 |
233 |
234 |
235 3--------2
236 2 |
237 | unfold ^ |
238 | ===> \_ |
239 | |
240 0------1 0--------1
241 It can be shown that this unfolding doesn't overlap itself but the
243 corners may touch, such as at the X=-2,Y=1 etc noted above.
244
246 See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
247 classes.
248 "$path = Math::PlanePath::DragonCurve->new ()"
250 "$path = Math::PlanePath::DragonCurve->new (arms => $int)"
251 Create and return a new path object.
252 The optional "arms" parameter can make 1 to 4 copies of the curve,
254 each arm successively advancing.
255 "($x,$y) = $path->n_to_xy ($n)"
257 Return the X,Y coordinates of point number $n on the path. Points
258 begin at 0 and if "$n < 0" then the return is an empty list.
259 Fractional positions give an X,Y position along a straight line
261 between the integer positions.
262 "$n = $path->xy_to_n ($x,$y)"
264 Return the point number for coordinates "$x,$y". If there's
265 nothing at "$x,$y" then return "undef".
266 The curve visits an "$x,$y" twice for various points (all the
268 "inside" points). The smaller of the two N values is returned.
269 "@n_list = $path->xy_to_n_list ($x,$y)"
271 Return a list of N point numbers for coordinates "$x,$y".
272 The origin 0,0 has "arms_count()" many N since it's the starting
274 point for each arm. Other points have up to two Ns for a given
275 "$x,$y". If arms=4 then every "$x,$y" except the origin has
276 exactly two Ns.
277 "$n = $path->n_start()"
279 Return 0, the first N in the path.
280 Level Methods
282 "($n_lo, $n_hi) = $path->level_to_n_range($level)"
283 Return "(0, 2**$level)", or for multiple arms return "(0, $arms *
284 2**$level + ($arms-1))".
285 There are 2^level segments in a curve level, so 2^level+1 points
287 numbered from 0. For multiple arms there are arms*(2^level+1)
288 points, numbered from 0 so n_hi = arms*(2^level+1)-1.
289
291 Various formulas for coordinates, lengths and area can be found in the
292 author's long mathematical write-up
293 <http://user42.tuxfamily.org/dragon/index.html>
295 X,Y to N
297 The current code uses the "DragonMidpoint" "xy_to_n()" by rotating -45
298 degrees and offsetting to the midpoints of the four edges around the
299 target X,Y. The "DragonMidpoint" algorithm then gives four candidate N
300 values and those which convert back to the desired X,Y in the
301 "DragonCurve" "n_to_xy()" are the results for "xy_to_n_list()".
302 Xmid,Ymid = X+Y, Y-X # rotate -45 degrees
304 for dx = 0 or -1
305 for dy = 0 or 1
306 N candidate = DragonMidpoint xy_to_n(Xmid+dx,Ymid+dy)
307 Since there's at most two "DragonCurve" Ns at a given X,Y the loop can
309 stop when two Ns are found.
310 Only the "leaving" edges will convert back to the target N, so only two
312 of the four edges actually need to be considered. Is there a way to
313 identify them? For arm 1 and 3 the leaving edges are up,down on odd
314 points (meaning sum X+Y odd) and right,left for even points (meaning
315 sum X+Y even). But for arm 2 and 4 it's the other way around. Without
316 an easy way to determine the arm this doesn't seem to help.
317 X,Y is Visited
319 Whether a given X,Y is visited by the curve can be determined from one
320 or two segments (rather then up to four for X,Y to N).
321 | S midpoint Xmid = X+Y
323 | Ymid = Y-X
324 *---T--X,Y--S---*
325 | T midpoint Xmid-1
326 | Ymid+1
327 Segment S is to the East of X,Y. The arm it falls on can be determined
329 as per "X,Y to N" in Math::PlanePath::DragonMidpoint. Numbering arm(S)
330 = 0,1,2,3 then
331 X,Y Visited
333 -----------
334 if arms_count()==4 yes # whole plane
335 if arm(S) < arms_count() yes
336 if arm(S)==2 and arms_count()==1 no
337 if arm(T) < arms_count() yes
338 This works because when two arms touch they approach and leave by a
340 right angle, without crossing. So two opposite segments S and T
341 identify the two possible arms coming to the X,Y point.
342 |
344 |
345 \
346 ---- ----
347 \
348 |
349 |
350 An arm only touches its immediate neighbour, ie. arm-1 or arm+1 mod 4.
352 This means if arm(S)==2 then arm(T) can only be 1,2,3, not 0. So if
353 "arms_count()" is 1 then arm(T) cannot be on the curve and no need to
354 run its segment check.
355 The only exception to the right-angle touching rule is at the origin
357 X,Y = 0,0. In that case Xmid,Ymid = 0,0 is on the first arm and X,Y is
358 correctly determined to be on the curve. If S was not to the East but
359 some other direction away from X,Y then this wouldn't be so.
360 Turn
362 At each point the curve always turns either left or right, it never
363 goes straight ahead. The bit above the lowest 1-bit in N gives the
364 turn direction.
365 N = 0b...z10000 (possibly no trailing 0s)
367 z bit Turn
369 ----- ----
370 0 left
371 1 right
372 For example N=12 is binary 0b1100, the lowest 1 bit is 0b_1__ and the
374 bit above that is a 1, which means turn to the right. Or N=18 is
375 binary 0b10010, the lowest 1 is 0b___1_ and the bit above that is 0, so
376 turn left there.
377 This z bit can be picked out with some bit twiddling
379 $mask = $n & -$n; # lowest 1 bit, 000100..00
381 $z = $n & ($mask << 1); # the bit above it
382 $turn = ($z == 0 ? 'left' : 'right');
383 This sequence is in Knuth volume 2 "Seminumerical Algorithms" answer to
385 section 4.5.3 question 41 and is called the "dragon sequence". It's
386 expressed there recursively as
387 d(0) = 1 # unused, since first turn at N=1
389 d(2N) = d(N) # shift down looking for low 1-bit
390 d(4N+1) = 0 # bit above lowest 1-bit is 0
391 d(4N+3) = 1 # bit above lowest 1-bit is 1
392 Next Turn
394 The bits also give the turn after next by looking at the bit above the
395 lowest 0-bit. This works because 011..11 + 1 = 100..00 so the bit
396 above the lowest 0 becomes the bit above the lowest 1.
397 N = 0b...w01111 (possibly no trailing 1s)
399 w bit Next Turn
401 ---- ---------
402 0 left
403 1 right
404 For example at N=12=0b1100 the lowest 0 is the least significant bit
406 0b___0, and above that is a 0 too, so at N=13 the turn is to the left.
407 Or for N=18=0b10010 the lowest 0 is again the least significant bit,
408 but above it is a 1, so at N=19 the turn is to the right.
409 This too can be found with some bit twiddling, as for example
411 $mask = $n ^ ($n+1); # low one and below 000111..11
413 $w = $n & ($mask + 1); # the bit above there
414 $turn = ($w == 0 ? 'left' : 'right');
415 Total Turn
417 The total turn is the count of 0<->1 transitions in the runs of bits of
418 N, which is the same as how many bit pairs of N (including overlaps)
419 are different so "01" or "10".
420 This can be seen from the segment replacements resulting from bits of
422 N,
423 N bits from high to low, start in "plain" state
425 plain state
427 0 bit -> no change
428 1 bit -> count left, and go to reversed state
429 reversed state
431 0 bit -> count left, and go to plain state
432 1 bit -> no change
433 The 0 or 1 counts are from the different side a segment expands on in
435 plain or reversed state. Segment A to B expands to an "L" shape bend
436 which is on the right in plain state, but on the left in reversed
437 state.
438 plain state reverse state
440 A = = = = B +
442 \ ^ 0bit / \
443 \ / turn / \ 1bit
444 0bit \ / 1bit left/ \
445 \ / turn / v
446 + left A = = = = B
447 In both cases a rotate of +45 degrees keeps the very first segment of
449 the whole curve in a fixed direction (along the X axis), which means
450 the south-east slope shown is no-change. This is the 0 of plain or the
451 1 of reversed. And the north-east slope which is the other new edge is
452 a turn towards the left.
453 It can be seen the "state" above is simply the previous bit, so the
455 effect for the bits of N is to count a left turn at each transition
456 from 0->1 or 1->0. Initial "plain" state means the infinite zero bits
457 at the high end of N are included. For example N=9 is 0b1001 so three
458 left turns for curve direction south to N=10 (as can be seen in the
459 diagram above).
460 1 00 1 N=9
462 ^ ^ ^
463 +-+--+---three transitions,
464 so three left turns for direction south
465 The transitions can also be viewed as a count of how many runs of
467 contiguous 0s or 1s, up to the highest 1-bit.
468 1 00 1 three blocks of 0s and 1s
470 This can be calculated by some bit twiddling with a shift and xor to
472 turn transitions into 1-bits which can then be counted, as per Jorg
473 Arndt (fxtbook section 1.31.3.1 "The Dragon Curve").
474 total turn = count_1_bits ($n ^ ($n >> 1))
476 The reversing structure of the curve shows up in the total turn at each
478 point. The total turns for a block of 2^N is followed by its own
479 reversal plus 1. For example,
480 ------->
482 N=0 to N=7 0, 1, 2, 1, 2, 3, 2, 1
483 N=15 to N=8 1, 2, 3, 2, 3, 4, 3, 2 each is +1
485 <-------
486 N to dX,dY
488 "n_to_dxdy()" is the "total turn" per above, or for fractional N then
489 an offset according to the "next turn" above. If using the bit
490 twiddling operators described then the two can be calculated
491 separately.
492 The current "n_to_dxdy()" code tries to support floating point or other
494 number types without bitwise XOR etc by processing bits high to low
495 with a state table which combines the calculations for total turn and
496 next turn. The state encodes
497 total turn 0 to 3
499 next turn 0 or 1
500 previous bit 0 or 1 (the bit above the current bit)
501 The "next turn" remembers the bit above lowest 0 seen so far (or 0
503 initially). The "total turn" counts 0->1 or 1->0 transitions. The
504 "previous bit" shows when there's a transition, or what bit is above
505 when a 0 is seen. It also works not to have this previous bit in the
506 state but instead pick out two bits each time.
507 At the end of bit processing any "previous bit" in state is no longer
509 needed and can be masked out to lookup the final four dx, dy, next dx,
510 next dy.
511
513 The Dragon curve is in Sloane's Online Encyclopedia of Integer
514 Sequences in many forms (and see "DragonMidpoint" for its forms too),
515 <http://oeis.org/A014577> (etc)
517 A038189 turn, 0=left,1=right, bit above lowest 1, extra 0
519 A089013 same, but initial extra 1
520 A082410 turn, 1=left,0=right, reversing complement, extra 0
521 A099545 turn, 1=left,3=right, as [odd part n] mod 4
522 so turn by 90 degrees * 1 or 3
523 A034947 turn, 1=left,-1=right, Jacobi (-1/n)
524 A112347 turn, 1=left,-1=right, Kronecker (-1/n), extra 0
525 A121238 turn, 1=left,-1=right, -1^(n + some partitions) extra 1
526 A014577 next turn, 1=left,0=right
527 A014707 next turn, 0=left,1=right
528 A014709 next turn, 1=left,2=right
529 A014710 next turn, 2=left,1=right
530 A143347 paperfolding constant, 0=left,1=right in binary
531 A090678 1=same turn as previous, 0=different
532 These numerous turn sequences differ only in having left or right
534 represented as 0, 1, -1, etc, and possibly "extra" initial 0 or 1 at
535 n=0 arising from the definitions and the first turn being at n=N=1.
536 The "next turn" forms begin at n=0 for turn at N=1 and so are the turn
537 at N=n+1.
538 A005811 direction total turn
540 A088748 direction total turn + 1
541 A037834 direction total turn - 1
542 A136004 direction total turn + 4
543 A246960 direction 0,1,2,3 = total turn mod 4
544 A173318 cumulative(total turn)
545 A164910 cumulative(total turn + 1)
546 A166242 2^(total turn), by double/halving
547 A000975 N of new maximum total turn, binary 10101...
548 A268411 direction of horizontals, 0=East, 1=West
549 A043724 N of East
550 A043725 N of North
551 A043726 N of West
552 A043727 N of South
553 A088431 turn sequence run lengths
555 A007400 2*runlength
556 A091072 N positions of the left turns, being odd part form 4K+1
558 A091067 N positions of the right turns, being odd part form 4K+3
559 A255068 N positions where next turn right
560 A060833 N positions where previous turn right
561 A106836 N steps between right turns
562 A088742 N steps between left turns
563 A255070 num right turns 1 to N
564 A236840 2* num right turns 1 to N
565 A003460 turns N=1 to N=2^n-1 packed as bits 1=left,0=right
566 low to high, then written in octal
567 A126937 points numbered like SquareSpiral (start N=0 and flip Y)
568 A038503 num segments East in level k
570 A038504 num segments North in level k
571 A038505 num segments West in level k
572 A000749 num segments South in level k
573 A146559 X at N=2^k, for k>=1, being Re((i+1)^k)
575 A009545 Y at N=2^k, for k>=1, being Im((i+1)^k)
576 A227036 boundary length N=0 to N=2^k
578 also right boundary length to N=2^(k+1)
579 A203175 left boundary length N=0 to N=2^k
580 = differences of total boundary
581 = squares on left boundary
582 A003476 squares on right boundary
584 = single points N=0 to N=2^(k-1) inclusive
585 A164395 single points N=0 to N=2^k-1 inclusive, for k=4 up
586 A003230 area enclosed N=0 to N=2^k, for k=4 up
588 same as double points
589 A003478 area enclosed by left side,
590 also area increment
591 A003477 area of each blob (touching enclosed unit squares)
592 A003479 join area between N=2^k replications
594 A003229 join area increment,
595 also area left side extra over doubling
596 A077949 same
597 A289265 growth rate r = 1.695 of boundaries etc
599 A272031 fractal dimension log(r)/log(sqrt(2))
600 For reference, "dragon-like" A059125 is similar to the turn sequence
602 A014707, but differs in having the "middle" values for each replication
603 come from successive values of the sequence itself, or some such.
604 A088431 and A007400
606 The run lengths A088431 and A007400 are from a continued fraction
607 expansion of an infinite sum
608 1 1 1 1 1 1
610 1 + - + - + -- + --- + ----- + ... + ------- + ...
611 2 4 16 256 65536 2^(2^k)
612 Jeffrey Shallit and independently M. Kmosek show how continued fraction
614 terms repeated in reverse give rise to this sort of power sum,
615 Jeffrey Shallit, "Simple Continued Fractions for Some Irrational
617 Numbers", Journal of Number Theory, volume 11, 1979, pages 209-217.
618 <http://www.cs.uwaterloo.ca/~shallit/papers.html>
619 <http://www.cs.uwaterloo.ca/~shallit/Papers/scf.ps>
620 (And which appears in Knuth "Art of Computer Programming", volume
622 2, section 4.5.3 exercise 41.)
623 A126937
625 The A126937 "SquareSpiral" numbering has the dragon curve and square
626 spiralling with their Y axes in opposite directions, as shown in its
627 a126937.pdf. So the dragon curve turns up towards positive Y but the
628 square spiral is numbered down towards negative Y (or vice versa).
629 "PlanePath" code for this starting at "$i=0" would be
630 my $dragon = Math::PlanePath::DragonCurve->new;
632 my $square = Math::PlanePath::SquareSpiral->new (n_start => 0);
633 my ($x, $y) = $dragon->n_to_xy ($i);
634 my $A126937_of_i = $square->xy_to_n ($x, -$y);
635
637 Math::PlanePath, Math::PlanePath::DragonRounded,
638 Math::PlanePath::DragonMidpoint, Math::PlanePath::R5DragonCurve,
639 Math::PlanePath::TerdragonCurve
640 Math::PlanePath::ComplexMinus, Math::PlanePath::ComplexPlus,
642 Math::PlanePath::CCurve, Math::PlanePath::AlternatePaper
643 Graph::Maker::TwindragonAreaTree
645 <http://rosettacode.org/wiki/Dragon_curve>
647
649 <http://user42.tuxfamily.org/math-planepath/index.html>
650
652 Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 Kevin Ryde
653 Math-PlanePath is free software; you can redistribute it and/or modify
655 it under the terms of the GNU General Public License as published by
656 the Free Software Foundation; either version 3, or (at your option) any
657 later version.
658 Math-PlanePath is distributed in the hope that it will be useful, but
660 WITHOUT ANY WARRANTY; without even the implied warranty of
661 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
662 General Public License for more details.
663 You should have received a copy of the GNU General Public License along
665 with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
666perl v5.30.1 2020-01-30 Math::PlanePath::DragonCurve(3)