File Coverage

blib/lib/Math/PlanePath/SierpinskiCurveStair.pm
Criterion Covered Total %
statement 169 196 86.2
branch 50 66 75.7
condition 18 22 81.8
subroutine 23 24 95.8
pod 8 8 100.0
total 268 316 84.8


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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              
20              
21             # math-image --path=SierpinskiCurveStair --lines --scale=10
22             #
23             # math-image --path=SierpinskiCurveStair,diagonal_length=1 --all --output=numbers_dash --offset=-10,-7 --size=78x30
24              
25              
26              
27             package Math::PlanePath::SierpinskiCurveStair;
28 1     1   8950 use 5.004;
  1         10  
29 1     1   6 use strict;
  1         2  
  1         38  
30 1     1   7 use List::Util 'min','max';
  1         1  
  1         138  
31              
32 1     1   7 use vars '$VERSION', '@ISA';
  1         2  
  1         63  
33             $VERSION = 128;
34 1     1   690 use Math::PlanePath;
  1         3  
  1         29  
35 1     1   395 use Math::PlanePath::Base::NSEW;
  1         3  
  1         41  
36             @ISA = ('Math::PlanePath::Base::NSEW',
37             'Math::PlanePath');
38              
39             use Math::PlanePath::Base::Generic
40 1         47 'is_infinite',
41 1     1   6 'round_nearest';
  1         2  
42             use Math::PlanePath::Base::Digits
43 1         70 'round_up_pow',
44 1     1   444 'round_down_pow';
  1         3  
45             *_divrem_mutate = \&Math::PlanePath::_divrem_mutate;
46              
47             # uncomment this to run the ### lines
48             #use Smart::Comments;
49              
50              
51 1     1   7 use constant n_start => 0;
  1         2  
  1         149  
52             sub x_negative {
53 6     6 1 81 my ($self) = @_;
54 6         15 return ($self->{'arms'} >= 3);
55             }
56             sub y_negative {
57 6     6 1 314 my ($self) = @_;
58 6         19 return ($self->{'arms'} >= 5);
59             }
60              
61 1         67 use constant parameter_info_array =>
62             [
63             { name => 'diagonal_length',
64             display => 'Diagonal Length',
65             type => 'integer',
66             minimum => 1,
67             default => 1,
68             width => 1,
69             description => 'Length of the diagonal in the base pattern.',
70             },
71             { name => 'arms',
72             share_key => 'arms_8',
73             display => 'Arms',
74             type => 'integer',
75             minimum => 1,
76             maximum => 8,
77             default => 1,
78             width => 1,
79             },
80 1     1   8 ];
  1         2  
81              
82 1     1   534 use Math::PlanePath::SierpinskiCurve;
  1         2  
  1         65  
83             *x_negative_at_n = \&Math::PlanePath::SierpinskiCurve::x_negative_at_n;
84             *y_negative_at_n = \&Math::PlanePath::SierpinskiCurve::y_negative_at_n;
85             *x_minimum = \&Math::PlanePath::SierpinskiCurve::x_minimum;
86             *sumxy_minimum = \&Math::PlanePath::SierpinskiCurve::sumxy_minimum;
87 1     1   7 use constant sumabsxy_minimum => 1;
  1         2  
  1         79  
88             *diffxy_minimum = \&Math::PlanePath::SierpinskiCurve::diffxy_minimum;
89 1     1   12 use constant absdiffxy_minimum => 1; # X=Y never occurs
  1         2  
  1         46  
90 1     1   5 use constant rsquared_minimum => 1; # minimum X=1,Y=0
  1         2  
  1         39  
91 1     1   7 use constant turn_any_straight => 0; # never straight
  1         2  
  1         1457  
92              
93              
94             #------------------------------------------------------------------------------
95              
96             sub new {
97 41     41 1 4201 my $self = shift->SUPER::new(@_);
98 41   100     248 $self->{'arms'} = max(1, min(8, $self->{'arms'} || 1));
99 41   100     126 $self->{'diagonal_length'} ||= 1;
100 41         93 return $self;
101             }
102              
103             # 20--21
104             # | |
105             # 18--19 22--23
106             # | |
107             # 16--17 24--25
108             # | |
109             # 15--14 27--26
110             # | |
111             # 4---5 13--12 29--28 36--37
112             # | | | | | |
113             # 2---3 6---7 10--11 30--31 34--35 38--39 42--43
114             # | | | | | | |
115             # 0---1 8---9 32--33 40--41
116              
117             # len=5
118             # N=0 to 9 is 10
119             # next N=0 to 41 is 42=4*10+2
120             # next is 4*42+2=166
121             # points(level) = 4*points(level-1)+2
122             #
123             # or side 5 points
124             # points(level) = 4*points(level-1)+1
125             # = 4*(4*points(level-2)+1)+1
126             # = 16*points(level-2) + 4 + 1
127             # = 64*points(level-3) + 16 + 4 + 1
128             # = 5 * 4^level + 1+...+4^(level-1)
129             # = 5 * 4^level + (4^level - 1) / 3
130             # = (15 * 4^level + 4^level - 1) / 3
131             # = (16 * 4^level - 1) / 3
132             # = (4^(level+2) - 1) / 3
133             # level=0 (16*1-1)/3=5
134             # level=1 (16*4-1)/3=21
135             # level=2 (16*16-1)/3=85
136             #
137             # n = (16 * 4^level - 1) / 3
138             # 3n+1 = 16 * 4^level
139             # 4^level = (3n+1)/16
140             # level = log4 ( (3n+1)/16)
141             # = log4(3n+1) - 2
142             # N=21 log4(64)-2=3-2=1
143             #
144             # nlen=4^(level+2)
145             # n = (nlen-1)/3
146             # next_n = (nlen/4-1)/3
147             # = (nlen-4)/3 /4
148             # = ((nlen-1)/3 -1) /4
149             #
150             # len=2,6,14
151             # len(k)=2*len(k-1) + 2
152             # = 2^k + 2*(2^(k-1)-1)
153             # = 2^k + 2^k - 2
154             # = 2*(2^k - 1)
155             # k=1 len=2*(2-1) = 2
156             # k=2 len=2*(4-1) = 6
157             # k=3 len=2*(8-1) = 14
158              
159             # len(k)-2=2*len(k-1)
160             # (len(k)-2)/2=len(k-1)
161             # len(k-1) = (len(k)-2)/2
162             # = len(k)/2-1
163             #
164             # ---------
165             # with P=2*L+1 points per side
166             # points(level) = 64*points(level-3) + 16 + 4 + 1
167             # = P*4^level + 1+...+4^(level-1)
168             # = P*4^level + (4^level - 1) / 3
169             # = (3P*4^level + 4^level - 1) / 3
170             # = ((3P+1)*4^level - 1) / 3
171             # = ((3*(2L+1)+1)*4^level - 1) / 3
172             # = ((6L+3+1)*4^level - 1) / 3
173             # = ((6L+4)*4^level - 1) / 3
174             # n = ((6L+4)*4^level - 1) / 3
175             # 3n+1 = (6L+4)*4^level
176             #
177             # len(k) = 2*len(k-1) + 2
178             # = 2*len(k-2) + 2 + 4
179             # = 2*len(k-3) + 2 + 4 + 8
180             # = 2^(k-1)*L + 2^k - 2
181             # = (L+2)*2^(k-1) - 2
182             # L=2 k=3 len=(2+2)*2^2-2=14
183             #
184             # ----------
185             # Nlevel = ((6L+4)*4^level - 1) / 3 - 1
186             # = ((6L+4)*4^level - 4) / 3
187             # Xlevel = (L+2)*2^level - 2 + 1
188             # = (L+2)*2^level - 1
189             #
190             # fill = Nlevel / (Xlevel*(Xlevel-1)/2)
191             # = (((6L+4)*4^level - 1) / 3 - 1) / (((L+2)*2^level - 1)*((L+2)*2^level - 2))
192             # -> (((6L+4)*4^level) / 3) / ((L+2)*2^level)^2
193             # = ((6L+4)*4^level) / ((L+2)^2*4^level) *2/3
194             # = ((6L+4)) / ((L+2)^2) * 2/3
195             # = 2*(3L+2) / ((L+2)^2) * 2/3
196             # = 4/3 * (3L+2)/(L+2)^2
197             # = (12L+8) / (3*L^2+12L+12)
198             # L=1 (12+8)/(3+12+12) = 20/27
199              
200              
201             sub n_to_xy {
202 19     19 1 959 my ($self, $n) = @_;
203             ### SierpinskiCurveStair n_to_xy(): $n
204              
205 19 50       54 if ($n < 0) {
206 0         0 return;
207             }
208 19 50       57 if (is_infinite($n)) {
209 0         0 return ($n,$n);
210             }
211              
212 19         34 my $frac;
213             {
214 19         38 my $int = int($n);
  19         31  
215 19         27 $frac = $n - $int; # inherit possible BigFloat
216 19 50       39 if ($frac) {
217 0         0 my ($x1,$y1) = $self->n_to_xy($int);
218 0         0 my ($x2,$y2) = $self->n_to_xy($int+$self->{'arms'});
219              
220 0         0 my $dx = $x2-$x1;
221 0         0 my $dy = $y2-$y1;
222 0         0 return ($frac*$dx + $x1, $frac*$dy + $y1);
223             }
224 19         30 $n = $int; # BigFloat int() gives BigInt, use that
225             }
226             ### $frac
227 19         29 my $zero = ($n * 0); # inherit bignum 0
228              
229 19         45 my $arm = _divrem_mutate ($n, $self->{'arms'});
230              
231 19         31 my $diagonal_length = $self->{'diagonal_length'};
232 19         31 my $diagonal_div = 6*$diagonal_length + 4;
233              
234 19         53 my ($nlen,$level) = round_down_pow ((3*$n+1)/$diagonal_div, 4);
235             ### $nlen
236             ### $level
237 19 50       42 if (is_infinite($level)) {
238 0         0 return $level;
239             }
240              
241 19         31 my $x = $zero;
242 19         28 my $y = $zero;
243 19         24 my $dx = 1;
244 19         27 my $dy = 0;
245              
246             # (L+2)*2^(level-1) - 2
247 19         34 my $len = ($diagonal_length+2)*2**$level - 2;
248 19         34 $nlen = ($diagonal_div*$nlen-1)/3;
249              
250 19         43 while ($level-- >= 0) {
251             ### at: "n=$n xy=$x,$y nlen=$nlen len=$len"
252              
253 52 100       127 if ($n < 2*$nlen+1) {
254 19 100       32 if ($n < $nlen) {
255             ### part 0 ...
256             } else {
257             ### part 1 ...
258 18         31 $x += ($len+1)*$dx - $len*$dy;
259 18         36 $y += ($len+1)*$dy + $len*$dx;
260 18         37 ($dx,$dy) = ($dy,-$dx); # rotate -90
261 18         29 $n -= $nlen;
262             }
263             } else {
264 33         49 $n -= 2*$nlen+1;
265 33 50       51 if ($n < $nlen) {
266             ### part 2 ...
267 0         0 $x += (2*$len+2)*$dx - $dy;
268 0         0 $y += (2*$len+2)*$dy + $dx;
269 0         0 ($dx,$dy) = (-$dy,$dx); # rotate +90
270             } else {
271             ### part 3 ...
272 33         54 $x += ($len+2)*$dx - ($len+2)*$dy;
273 33         46 $y += ($len+2)*$dy + ($len+2)*$dx;
274 33         49 $n -= $nlen;
275             }
276             }
277              
278 52         79 $nlen = ($nlen-1)/4;
279 52         105 $len = $len/2-1;
280             }
281              
282 19         34 my $lowdigit_x = int(($n+1)/2);
283 19 100       38 if ($n == 2*$diagonal_length+1) { $lowdigit_x -= 2; }
  5         8  
284 19         28 my $lowdigit_y = int($n/2);
285              
286             ### final: "n=$n xy=$x,$y dxdy=$dx,$dy"
287             ### $lowdigit_x
288             ### $lowdigit_y
289              
290 19         32 $x += $lowdigit_x*$dx - $lowdigit_y*$dy + 1; # +1 start at x=1,y=0
291 19         25 $y += $lowdigit_x*$dy + $lowdigit_y*$dx;
292              
293 19 50       39 if ($arm & 1) {
294 0         0 ($x,$y) = ($y,$x); # mirror 45
295             }
296 19 50       37 if ($arm & 2) {
297 0         0 ($x,$y) = (-1-$y,$x); # rotate +90
298             }
299 19 50       37 if ($arm & 4) {
300 0         0 $x = -1-$x; # rotate 180
301 0         0 $y = -1-$y;
302             }
303              
304 19         50 return ($x,$y);
305             }
306              
307             sub xy_to_n {
308 19776     19776 1 198110 my ($self, $x, $y) = @_;
309             ### SierpinskiCurveStair xy_to_n(): "$x, $y"
310              
311 19776         38298 $x = round_nearest($x);
312 19776         37190 $y = round_nearest($y);
313              
314 19776         29056 my $arm = 0;
315 19776 100       35355 if ($y < 0) {
316 2720         3592 $arm = 4;
317 2720         3773 $x = -1-$x; # rotate -180
318 2720         3590 $y = -1-$y;
319             }
320 19776 100       32220 if ($x < 0) {
321 4608         6113 $arm += 2;
322 4608         7873 ($x,$y) = ($y, -1-$x); # rotate -90
323             }
324 19776 100       33935 if ($y > $x) { # second octant
325 9344         12315 $arm++;
326 9344         17622 ($x,$y) = ($y,$x); # mirror 45
327             }
328              
329 19776         29834 my $arms = $self->{'arms'};
330 19776 100       33547 if ($arm >= $arms) {
331 3680         6822 return undef;
332             }
333              
334 16096         22097 $x -= 1;
335 16096 100 100     46439 if ($x < 0 || $x < $y) {
336 1008         2061 return undef;
337             }
338             ### x adjust to zero: "$x,$y"
339             ### assert: $x >= 0
340             ### assert: $y >= 0
341              
342             # len=2*(2^level - 1)
343             # len/2+1 = 2^level
344             # 2^level = len/2+1
345             # 2^(level+1) = len+2
346              
347             # len=(L+2)*2^(level-1) - 2
348             # (len+2)/(L+2) = 2^(level-1)
349              
350 15088         22792 my $diagonal_length = $self->{'diagonal_length'};
351 15088         36259 my ($len,$level) = round_down_pow (($x+1)/($diagonal_length+2), 2);
352             ### $level
353             ### $len
354 15088 50       31220 if (is_infinite($level)) {
355 0         0 return $level;
356             }
357              
358 15088         25997 my $n = 0;
359 15088         28033 my $nlen = ((6*$diagonal_length+4)*$len*$len-1)/3;
360 15088         24492 $len *= ($self->{'diagonal_length'}+2);
361             ### $len
362             ### $nlen
363              
364 15088         19556 my $n_last_1;
365 15088         26321 foreach (0 .. $level) {
366             ### at: "loop=$_ x=$x,y=$y n=$n nlen=$nlen len=$len diag cmp ".(2*$len-2)
367             ### assert: $x >= 0
368             ### assert: $y >= 0
369              
370 37510 100       67682 if ($x+$y <= 2*$len-2) {
371             ### part 0 or 1...
372 20665 100       34596 if ($x < $len-1) {
373             ### part 0 ...
374 6256         8659 $n_last_1 = 0;
375             } else {
376             ### part 1 ...
377 14409         25408 ($x,$y) = ($len-2-$y, $x-($len-1)); # shift then rotate +90
378 14409         20608 $n += $nlen;
379 14409         19932 $n_last_1 = 1;
380             }
381             } else {
382 16845         26814 $n += 2*$nlen + 1; # +1 for middle point
383             ### part 2 or 3 ...
384 16845 100       27619 if ($y < $len) {
385             ### part 2...
386 8671         15518 ($x,$y) = ($y-1, 2*$len-2-$x); # shift y-1 then rotate -90
387 8671         11735 $n_last_1 = 0;
388             } else {
389             #### digit 3...
390 8174         10803 $x -= $len;
391 8174         10723 $y -= $len;
392 8174         10709 $n += $nlen;
393             }
394 16845 100       29107 if ($x < 0) {
395 81         198 return undef;
396             }
397             }
398 37429         51888 $len /= 2;
399 37429         58753 $nlen = ($nlen-1)/4;
400             }
401              
402             ### end at: "x=$x,y=$y n=$n last2=$n_last_1"
403             ### assert: $x >= 0
404             ### assert: $y >= 0
405              
406 15007 100 100     47601 if ($x == $y || $x == $y+1) {
    100 100        
      100        
407 8932         12884 $n += $x+$y;
408             } elsif ($n_last_1 && $x == $diagonal_length-1 && $y == $diagonal_length) {
409             # in between diagonals
410 458         780 $n += 2*$diagonal_length+1;
411             } else {
412 5617         13000 return undef;
413             }
414              
415 9390         19620 return $n*$arms + $arm;
416             }
417              
418             # not exact
419             sub rect_to_n_range {
420 32     32 1 192 my ($self, $x1,$y1, $x2,$y2) = @_;
421             ### SierpinskiCurveStair rect_to_n_range(): "$x1,$y1 $x2,$y2"
422              
423 32         71 $x1 = round_nearest ($x1);
424 32         87 $x2 = round_nearest ($x2);
425 32         72 $y1 = round_nearest ($y1);
426 32         70 $y2 = round_nearest ($y2);
427 32 50       87 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
428 32 50       85 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
429              
430             # x2
431             # y2 +-------+ *
432             # | | *
433             # y1 +-------+ *
434             # *
435             # *
436             # *
437             # ------------------
438             #
439             #
440             # *
441             # x1 * x2 *
442             # +-----*-+y2*
443             # | *| *
444             # | * *
445             # | |* *
446             # | | **
447             # +-------+y1*
448             # ----------------
449             #
450 32         61 my $arms = $self->{'arms'};
451 32 100 33     248 if (($arms <= 4
    50 33        
452             ? ($y2 < 0 # y2 negative, nothing ...
453             || ($arms == 1 && $x2 <= $y1)
454             || ($arms == 2 && $x2 < 0)
455             || ($arms == 3 && $x2 < -$y2))
456              
457             # arms >= 5
458             : ($y2 < 0
459             && (($arms == 5 && $x1 >= $y2)
460             || ($arms == 6 && $x1 >= 0)
461             || ($arms == 7 && $x1 > 3-$y2))))) {
462             ### rect outside octants, for arms: $arms
463             ### $x1
464             ### $y2
465 0         0 return (1,0);
466             }
467              
468 32         63 my $max = $x2; # arms 1,8 using X, starting at X=1
469 32 100       79 if ($arms >= 2) {
470             # arms 2,3 upper using Y, starting at Y=1
471 28         80 _apply_max ($max, $y2);
472              
473 28 100       76 if ($arms >= 4) {
474             # arms 4,5 right using X, starting at X=-2
475 20         59 _apply_max ($max, -1-$x1);
476              
477 20 100       62 if ($arms >= 6) {
478             # arms 6,7 down using Y, starting at Y=-2
479 12         55 _apply_max ($max, -1-$y1);
480             }
481             }
482             }
483             ### $max
484              
485              
486             # points(level) = (4^(level+2) - 1) / 3
487             # Nlast(level) = (4^(level+2) - 1) / 3 - 1
488             # = (4^(level+2) - 4) / 3
489             # then + arms-1 for last of arms
490             # Nhi = Nlast(level) * arms + arms-1
491             # = (Nlast(level + 1)) * arms - 1
492             # = ((4^(level+2) - 4) / 3 + 1) * arms - 1
493             # = ((4^(level+2) - 1) / 3) * arms - 1
494             #
495             # len(level) = = (L+2)*2^(level-1) - 2
496             # points(level) = ((3*P+1)*4^level - 1) / 3
497             #
498 32         110 my ($pow,$level) = round_down_pow ($max/($self->{'diagonal_length'}+2),
499             2);
500             return (0,
501 32         132 ((6*$self->{'diagonal_length'}+4)*4*$pow*$pow - 1) / 3
502             * $arms - 1);
503             }
504              
505             # set $_[0] to the max of $_[0] and $_[1]
506             sub _apply_max {
507             ### _apply_max(): "$_[0] cf $_[1]"
508 60 100   60   145 unless ($_[0] > $_[1]) {
509 28         56 $_[0] = $_[1];
510             }
511             }
512              
513              
514             #------------------------------------------------------------------------------
515              
516             # Nlevel = ((3L+2)*4^level - 5) / 3
517             # LevelPoints = Nlevel+1
518             # Nlevel(arms) = (Nlevel+1)*arms - 1
519             #
520             # Eg. L=1 level=1 (5*4-5)/3 = 5
521             # arms=8 ((5*4-5)/3+1)*8 - 1 = 47
522             #
523              
524             sub level_to_n_range {
525 12     12 1 719 my ($self, $level) = @_;
526             return (0,
527             (4**$level * (3*$self->{'diagonal_length'}+2) - 2) / 3
528 12         55 * $self->{'arms'} - 1);
529             }
530             sub n_to_level {
531 0     0 1   my ($self, $n) = @_;
532 0 0         if ($n < 0) { return undef; }
  0            
533 0 0         if (is_infinite($n)) { return $n; }
  0            
534 0           $n = round_nearest($n);
535 0           _divrem_mutate ($n, $self->{'arms'});
536 0           my $diagonal_div = 3*$self->{'diagonal_length'} + 2;
537 0           my ($pow,$exp) = round_up_pow ((3*$n+3) / (3*$self->{'diagonal_length'}+2),
538             4);
539 0           return $exp;
540             }
541              
542             #------------------------------------------------------------------------------
543             1;
544             __END__