File Coverage

blib/lib/Math/PlanePath/KochCurve.pm
Criterion Covered Total %
statement 198 241 82.1
branch 52 80 65.0
condition 34 52 65.3
subroutine 30 37 81.0
pod 6 6 100.0
total 320 416 76.9


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018 Kevin Ryde
2              
3             # This file is part of Math-PlanePath.
4             #
5             # Math-PlanePath is free software; you can redistribute it and/or modify
6             # it under the terms of the GNU General Public License as published by the
7             # Free Software Foundation; either version 3, or (at your option) any later
8             # version.
9             #
10             # Math-PlanePath is distributed in the hope that it will be useful, but
11             # WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
12             # or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13             # for more details.
14             #
15             # You should have received a copy of the GNU General Public License along
16             # with Math-PlanePath. If not, see .
17              
18              
19             # math-image --path=KochCurve --lines --scale=10
20             # math-image --path=KochCurve --all --scale=10
21              
22             # continuous but nowhere differentiable
23             #
24             # Sur une courbe continue sans tangente, obtenue par une construction
25             # géométrique élémentaire
26             #
27             # http://www.nku.edu/~curtin/grenouille.html
28             # http://www.nku.edu/~curtin/koch_171.jpg
29             #
30             # Cesàro, "Remarques sur la courbe de von Koch." Atti della
31             # R. Accad. della Scienze fisiche e matem. Napoli 12, No. 15, 1-12,
32             # 1905. Reprinted as §228 in Opere scelte, a cura dell'Unione matematica
33             # italiana e col contributo del Consiglio nazionale delle ricerche, Vol. 2:
34             # Geometria, analisi, fisica matematica. Rome: Edizioni Cremonese,
35             # pp. 464-479, 1964.
36             #
37             # Thue-Morse count 1s mod 2 is net direction
38             # Toeplitz first diffs is turn sequence +1 or -1
39             #
40             # J. Ma and J.A. Holdener. When Thue-Morse Meets Koch. In Fractals:
41             # Complex Geometry, Patterns, and Scaling in Nature and Society, volume 13,
42             # pages 191-206, 2005.
43             # http://personal.kenyon.edu/holdenerj/StudentResearch/WhenThueMorsemeetsKochJan222005.pdf
44             #
45             # F.M. Dekking. On the distribution of digits in arithmetic sequences. In
46             # Seminaire de Theorie des Nombres de Bordeaux, volume 12, pages 3201-3212,
47             # 1983.
48             #
49              
50              
51              
52             package Math::PlanePath::KochCurve;
53 4     4   7716 use 5.004;
  4         18  
54 4     4   22 use strict;
  4         5  
  4         151  
55 4     4   22 use List::Util 'sum','first';
  4         12  
  4         381  
56              
57 4     4   26 use vars '$VERSION', '@ISA';
  4         6  
  4         235  
58             $VERSION = 127;
59 4     4   604 use Math::PlanePath;
  4         10  
  4         210  
60             @ISA = ('Math::PlanePath');
61             *_divrem_mutate = \&Math::PlanePath::_divrem_mutate;
62              
63             use Math::PlanePath::Base::Generic
64 4         204 'is_infinite',
65 4     4   23 'round_nearest';
  4         8  
66             use Math::PlanePath::Base::Digits
67 4         232 'round_down_pow',
68             'round_up_pow',
69             'digit_split_lowtohigh',
70 4     4   410 'digit_join_lowtohigh';
  4         8  
71              
72             # uncomment this to run the ### lines
73             # use Smart::Comments;
74              
75              
76 4     4   24 use constant n_start => 0;
  4         7  
  4         193  
77 4     4   20 use constant class_x_negative => 0;
  4         7  
  4         198  
78 4     4   22 use constant class_y_negative => 0;
  4         8  
  4         161  
79 4     4   29 use constant diffxy_minimum => 0; # X>=Y octant so X-Y>=0
  4         8  
  4         187  
80 4     4   20 use constant dx_minimum => -2;
  4         7  
  4         164  
81 4     4   20 use constant dx_maximum => 2;
  4         8  
  4         176  
82 4     4   29 use constant dy_minimum => -1;
  4         8  
  4         198  
83 4     4   21 use constant dy_maximum => 1;
  4         8  
  4         288  
84             *_UNDOCUMENTED__dxdy_list = \&Math::PlanePath::_UNDOCUMENTED__dxdy_list_six;
85 4     4   25 use constant absdx_minimum => 1; # never vertical
  4         12  
  4         173  
86 4     4   21 use constant dsumxy_minimum => -2; # diagonals
  4         7  
  4         178  
87 4     4   21 use constant dsumxy_maximum => 2;
  4         8  
  4         175  
88 4     4   27 use constant ddiffxy_minimum => -2;
  4         8  
  4         163  
89 4     4   21 use constant ddiffxy_maximum => 2;
  4         15  
  4         201  
90 4     4   24 use constant dir_maximum_dxdy => (1,-1); # South-East
  4         8  
  4         201  
91 4     4   22 use constant turn_any_straight => 0; # never straight
  4         6  
  4         7084  
92              
93              
94             #------------------------------------------------------------------------------
95              
96             sub n_to_xy {
97 1456     1456 1 42674 my ($self, $n) = @_;
98             ### KochCurve n_to_xy(): $n
99              
100             # secret negatives to -.5
101 1456 50       2246 if (2*$n < -1) { return; }
  0         0  
102 1456 50       2249 if (is_infinite($n)) { return ($n,$n); }
  0         0  
103              
104 1456         2066 my $x;
105             my $y;
106             {
107 1456         1532 my $int = int($n);
  1456         1593  
108 1456         1712 $x = 2 * ($n - $int); # usually positive, but n=-0.5 gives x=-0.5
109 1456         1525 $y = $x * 0; # inherit possible bigrat 0
110 1456         1648 $n = $int; # BigFloat int() gives BigInt, use that
111             }
112              
113 1456         1634 my $len = $y+1; # inherit bignum 1
114 1456         2247 foreach my $digit (digit_split_lowtohigh($n,4)) {
115             ### at: "$x,$y digit=$digit"
116              
117 6138 100       9077 if ($digit == 0) {
    100          
    100          
118              
119             } elsif ($digit == 1) {
120 1691         2828 ($x,$y) = (($x-3*$y)/2 + 2*$len, # rotate +60
121             ($x+$y)/2);
122              
123             } elsif ($digit == 2) {
124 1636         2899 ($x,$y) = (($x+3*$y)/2 + 3*$len, # rotate -60
125             ($y-$x)/2 + $len);
126              
127             } else {
128             ### assert: $digit==3
129 1598         1747 $x += 4*$len;
130             }
131 6138         7048 $len *= 3;
132             }
133              
134             ### final: "$x,$y"
135 1456         2814 return ($x,$y);
136             }
137              
138             sub xy_to_n {
139 5645     5645 1 31707 my ($self, $x, $y) = @_;
140             ### KochPeaks xy_to_n(): "$x, $y"
141              
142 5645         8048 $x = round_nearest ($x);
143 5645         8320 $y = round_nearest ($y);
144 5645 100 100     18286 if ($y < 0 || $x < 0 || (($x ^ $y) & 1)) {
      66        
145             ### neg y or parity different ...
146 517         701 return undef;
147             }
148 5128   100     11432 my ($len,$level) = round_down_pow(($x/2)||1, 3);
149             ### $level
150             ### $len
151 5128 50       8058 if (is_infinite($level)) {
152 0         0 return $level;
153             }
154              
155 5128         6706 my $n = 0;
156 5128         6841 foreach (0 .. $level) {
157 16522         16584 $n *= 4;
158             ### at: "level=$level len=$len x=$x,y=$y n=$n"
159 16522 100       20931 if ($x < 3*$len) {
160 8604 100       10731 if ($x < 2*$len) {
161             ### digit 0 ...
162             } else {
163             ### digit 1 ...
164 2437         2653 $x -= 2*$len;
165 2437         3681 ($x,$y) = (($x+3*$y)/2, # rotate -60
166             ($y-$x)/2);
167 2437         2764 $n += 1;
168             }
169             } else {
170 7918         8411 $x -= 4*$len;
171             ### digit 2 or 3 to: "x=$x"
172 7918 100       9433 if ($x < $y) { # before diagonal
173             ### digit 2...
174 3460         3577 $x += $len;
175 3460         3493 $y -= $len;
176 3460         5286 ($x,$y) = (($x-3*$y)/2, # rotate +60
177             ($x+$y)/2);
178 3460         3762 $n += 2;
179             } else {
180             #### digit 3...
181 4458         4667 $n += 3;
182             }
183             }
184 16522         19472 $len /= 3;
185             }
186             ### end at: "x=$x,y=$y n=$n"
187 5128 100 100     8242 if ($x != 0 || $y != 0) {
188 4947         7921 return undef;
189             }
190 181         318 return $n;
191             }
192              
193             # level extends to x= 2*3^level
194             # level = log3(x/2)
195             #
196             # exact
197             sub rect_to_n_range {
198 29     29 1 2587 my ($self, $x1,$y1, $x2,$y2) = @_;
199             ### KochCurve rect_to_n_range(): "$x1,$y1 $x2,$y2"
200              
201 29         62 $x1 = round_nearest ($x1);
202 29         44 $x2 = round_nearest ($x2);
203 29         70 $y1 = round_nearest ($y1);
204 29         49 $y2 = round_nearest ($y2);
205 29 50       55 if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  0         0  
206 29 50       39 if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }
  0         0  
207              
208 29 50 33     119 if ($x2 < 0 || $y2 < 0
      33        
209             || 3*$y1 > $x2 ) { # above line Y=X/3
210 0         0 return (1,0);
211             }
212              
213             # \
214             # \
215             # * \
216             # / \ \
217             # o-+-* *-+-e \
218             # 0 3 6
219             #
220             # 3*Y+X/2 - (Y!=0)
221             #
222             # /
223             # *-+-*
224             # \
225             # * *
226             # / \ /
227             # o-+-* *-+-*
228             # 0 3 6 X/2
229             #
230 29         95 my ($len, $level) = round_down_pow ($x2/2, 3);
231 29         59 return _rect_to_n_range_rot ($len, $level, 0, $x1,$y1, $x2,$y2);
232              
233              
234              
235             # (undef, my $level) = round_down_pow ($x2/2, 3);
236             # ### $level
237             # return (0, 4**($level+1)-1);
238             }
239              
240              
241             my @dir6_to_dx = (2, 1,-1,-2, -1, 1);
242             my @dir6_to_dy = (0, 1, 1, 0, -1,-1);
243             my @max_digit_to_rot = (1, -2, 1, 0);
244             my @min_digit_to_rot = (0, 1, -2, 1);
245             my @max_digit_to_offset = (-1, -1, -1, 2);
246              
247             sub _rect_to_n_range_rot {
248 29     29   49 my ($initial_len, $level_max, $initial_rot, $x1,$y1, $x2,$y2) = @_;
249             ### KochCurve _rect_to_n_range_rot(): "$x1,$y1 $x2,$y2 len=$initial_len level=$level_max rot=$initial_rot"
250              
251 29         38 my ($rot, $len, $x, $y);
252             my $overlap = sub {
253             ### overlap: "$x,$y len=$len rot=$rot"
254              
255 346 100   346   479 if ($len == 1) {
256 144   100     802 return ($x >= $x1 && $x <= $x2
257             && $y >= $y1 && $y <= $y2);
258             }
259 346         216 my $len = $len / 3;
260              
261 346 100       263 if ($rot < 3) {
262 153 100       218 if ($rot == 0) {
    100          
263             # *
264             # / \
265             # o-+-* *-+-.
266 75   100     331 return ($y <= $y2 # bottom before end
267             && $y+$len >= $y1
268             && $x <= $x2
269             && $x+6*$len > $x1); # right before end, exclusive
270             } elsif ($rot == 1) {
271             # .
272             # /
273             # *-+-*
274             # \
275             # * +-----
276             # / |x1,y2
277             # o
278 55   100     189 return ($x <= $x2 # left before end
279             && $y+3*$len > $y1 # top after start, exclusive
280             && $y-$x <= $y2-$x1); # diag before corner
281             } else {
282             # . |x1,y1
283             # \ +-----
284             # *
285             # /
286             # *-+-*
287             # \
288             # o
289 23   100     75 return ($y <= $y2 # bottom before end
290             && $x-3*$len <=$x2 # left before end
291             && $y+$x >= $y1+$x1); # diag after corner
292             }
293             } else {
294 193 100       81 if ($rot == 3) {
    100          
295             # .-+-* *-+-o
296             # \ /
297             # *
298 8   66     29 return ($y >= $y1 # top after start
299             && $y-$len <= $y2 # bottom before end
300             && $x >= $x1 # right after start
301             && $x-6*$len < $x2); # left before end, exclusive
302             } elsif ($rot == 4) {
303             # x2,y1| o
304             # -----+ /
305             # *
306             # \
307             # *-+-*
308             # /
309             # .
310 3   66     15 return ($x >= $x1 # right after start
311             && $y-3*$len < $y2 # bottom before end, exclusive
312             && $y-$x >= $y1-$x2); # diag after corner
313             } else {
314             # o
315             # \
316             # *-+-*
317             # /
318             # *
319             # -----+ \
320             # x2,y2| .
321 182   100     150 return ($y >= $y1 # top after start
322             && $x+3*$len >= $x1 # right after start
323             && $y+$x <= $y2+$x2); # diag before corner
324             }
325             }
326 29         119 };
327              
328 29         41 my $zero = 0*$x1*$x2*$y1*$y2;
329 29         57 my @lens = ($initial_len);
330 29         38 my $n_hi;
331 29         32 $rot = $initial_rot;
332 29         30 $len = $initial_len;
333 29         30 $x = $zero;
334 29         29 $y = $zero;
335 29         35 my @digits = (4);
336              
337 29         32 for (;;) {
338 202         237 my $digit = --$digits[-1];
339             ### max at: "digits=".join(',',@digits)." xy=$x,$y len=$len"
340              
341 202 100       294 if ($digit < 0) {
342 2         4 pop @digits;
343 2 50       4 if (! @digits) {
344             ### nothing found to level_max ...
345 0         0 return (1, 0);
346             }
347             ### end of digits, backtrack ...
348 2         4 $len = $lens[$#digits];
349 2         4 next;
350             }
351              
352 200         221 my $offset = $max_digit_to_offset[$digit];
353 200         241 $rot = ($rot - $max_digit_to_rot[$digit]) % 6;
354 200         241 $x += $dir6_to_dx[$rot] * $offset * $len;
355 200         213 $y += $dir6_to_dy[$rot] * $offset * $len;
356              
357             ### $offset
358             ### $rot
359              
360 200 100       245 if (&$overlap()) {
361 74 100       115 if ($#digits >= $level_max) {
362             ### yes overlap, found n_hi ...
363             ### digits: join(',',@digits)
364             ### n_hi: _digit_join_hightolow (\@digits, 4, $zero)
365 29         48 $n_hi = _digit_join_hightolow (\@digits, 4, $zero);
366 29         52 last;
367             }
368             ### yes overlap, descend ...
369 45         65 push @digits, 4;
370 45   66     128 $len = ($lens[$#digits] ||= $len/3);
371             } else {
372             ### no overlap, next digit ...
373             }
374             }
375              
376 29         32 $rot = $initial_rot;
377 29         32 $x = $zero;
378 29         30 $y = $zero;
379 29         31 $len = $initial_len;
380 29         57 @digits = (-1);
381              
382 29         36 for (;;) {
383 147         161 my $digit = ++$digits[-1];
384             ### min at: "digits=".join(',',@digits)." xy=$x,$y len=$len"
385              
386 147 100       198 if ($digit > 3) {
387 1         2 pop @digits;
388 1 50       3 if (! @digits) {
389             ### oops, n_lo not found to level_max ...
390 0         0 return (1, 0);
391             }
392             ### end of digits, backtrack ...
393 1         1 $len = $lens[$#digits];
394 1         2 next;
395             }
396              
397             ### $digit
398             ### rot increment: $min_digit_to_rot[$digit]
399 146         165 $rot = ($rot + $min_digit_to_rot[$digit]) % 6;
400              
401 146 100       161 if (&$overlap()) {
402 73 100       119 if ($#digits >= $level_max) {
403             ### yes overlap, found n_lo ...
404             ### digits: join(',',@digits)
405             ### n_lo: _digit_join_hightolow (\@digits, 4, $zero)
406 29         44 return (_digit_join_hightolow (\@digits, 4, $zero),
407             $n_hi);
408             }
409             ### yes overlap, descend ...
410 44         120 push @digits, -1;
411 44   33     80 $len = ($lens[$#digits] ||= $len/3);
412              
413             } else {
414             ### no overlap, next digit ...
415 73         86 $x += $dir6_to_dx[$rot] * $len;
416 73         420 $y += $dir6_to_dy[$rot] * $len;
417             }
418             }
419             }
420              
421             # $aref->[0] high digit
422             sub _digit_join_hightolow {
423 58     58   84 my ($aref, $radix, $zero) = @_;
424 58         93 my @lowtohigh = reverse @$aref;
425 58         128 return digit_join_lowtohigh(\@lowtohigh, $radix, $zero);
426             }
427              
428              
429             my @digit_to_dir = (0, 1, -1, 0);
430             my @digit_to_nextturn = (1, # digit=1 (with +1 for "next" N)
431             -2, # digit=2
432             1); # digit=3
433             sub n_to_dxdy {
434 20     20 1 90 my ($self, $n) = @_;
435             ### n_to_dxdy(): $n
436              
437 20 50       35 if ($n < 0) {
438 0         0 return; # first direction at N=0
439             }
440 20 50       31 if (is_infinite($n)) {
441 0         0 return ($n,$n);
442             }
443              
444 20         33 my $int = int($n);
445 20         25 $n -= $int;
446 20         35 my @ndigits = digit_split_lowtohigh($int,4);
447              
448 20         36 my $dir6 = sum(0, map {$digit_to_dir[$_]} @ndigits) % 6;
  111         166  
449 20         31 my $dx = $dir6_to_dx[$dir6];
450 20         26 my $dy = $dir6_to_dy[$dir6];
451              
452 20 100       30 if ($n) {
453             # fraction part
454              
455             # lowest non-3 digit, or zero if all 3s (0 above high digit)
456 17     23   70 $dir6 += $digit_to_nextturn[ first {$_!=3} @ndigits, 0 ];
  23         33  
457 17         38 $dir6 %= 6;
458 17         30 $dx += $n*($dir6_to_dx[$dir6] - $dx);
459 17         24 $dy += $n*($dir6_to_dy[$dir6] - $dy);
460             }
461 20         48 return ($dx, $dy);
462             }
463              
464             sub _UNTESTED__n_to_dir6 {
465 0     0     my ($self, $n) = @_;
466 0 0         if ($n < 0) {
467 0           return undef; # first direction at N=0
468             }
469 0 0         if (is_infinite($n)) {
470 0           return ($n,$n);
471             }
472 0   0       return (sum (map {$digit_to_dir[$_]} digit_split_lowtohigh($n,4))
473             || 0) # if empty
474             % 6;
475             }
476              
477             my @n_to_turn6 = (undef,
478             1, # +60 degrees
479             -2, # -120 degrees
480             1); # +60 degrees
481             sub _UNTESTED__n_to_turn6 {
482 0     0     my ($self, $n) = @_;
483 0 0         if (is_infinite($n)) {
484 0           return undef;
485             }
486 0           while ($n) {
487 0           my $digit = _divrem_mutate($n,4);
488 0 0         if ($digit) {
489             # lowest non-zero digit
490 0           return $n_to_turn6[$digit];
491             }
492             }
493 0           return 0;
494             }
495             sub _UNTESTED__n_to_turn_LSR {
496 0     0     my ($self, $n) = @_;
497 0   0       my $turn6 = $self->_UNTESTED__n_to_turn6($n) || return undef;
498 0 0         return ($turn6 > 0 ? 1 : -1);
499             }
500             sub _UNTESTED__n_to_turn_left {
501 0     0     my ($self, $n) = @_;
502 0   0       my $turn6 = $self->_UNTESTED__n_to_turn6($n) || return undef;
503 0 0         return ($turn6 > 0 ? 1 : 0);
504             }
505             sub _UNTESTED__n_to_turn_right {
506 0     0     my ($self, $n) = @_;
507 0   0       my $turn6 = $self->_UNTESTED__n_to_turn6($n) || return undef;
508 0 0         return ($turn6 < 0 ? 1 : 0);
509             }
510              
511             #------------------------------------------------------------------------------
512             # levels
513              
514             sub level_to_n_range {
515 0     0 1   my ($self, $level) = @_;
516 0           return (0, 4**$level);
517             }
518             sub n_to_level {
519 0     0 1   my ($self, $n) = @_;
520 0 0         if ($n < 0) { return undef; }
  0            
521 0 0         if (is_infinite($n)) { return $n; }
  0            
522 0           $n = round_nearest($n);
523 0           my ($pow, $exp) = round_up_pow ($n, 4);
524 0           return $exp;
525             }
526              
527              
528             #------------------------------------------------------------------------------
529             1;
530             __END__