File Coverage

blib/lib/Math/PlanePath/ToothpickReplicate.pm
Criterion Covered Total %
statement 109 226 48.2
branch 31 110 28.1
condition 1 15 6.6
subroutine 15 20 75.0
pod 6 6 100.0
total 162 377 42.9


line stmt bran cond sub pod time code
1             # Copyright 2012, 2013, 2014, 2015 Kevin Ryde
2              
3             # This file is part of Math-PlanePath-Toothpick.
4             #
5             # Math-PlanePath-Toothpick is free software; you can redistribute it and/or
6             # modify it under the terms of the GNU General Public License as published
7             # by the Free Software Foundation; either version 3, or (at your option) any
8             # later version.
9             #
10             # Math-PlanePath-Toothpick is distributed in the hope that it will be
11             # useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
12             # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
13             # Public License for more details.
14             #
15             # You should have received a copy of the GNU General Public License along
16             # with Math-PlanePath-Toothpick. If not, see .
17              
18             # block_order => 'AB123'
19             # block_order => 'A1B32' is depth first and finite parts first,
20             # in parts=1 where single infinite spine
21             #
22             # maybe tree methods same structure as ToothpickTree
23             #
24              
25             # cf A175262 odd binary length and middle digit 1
26             # A175263 odd binary length and middle digit 0
27             #
28              
29             package Math::PlanePath::ToothpickReplicate;
30 1     1   1939 use 5.004;
  1         3  
31 1     1   6 use strict;
  1         2  
  1         43  
32             #use List::Util 'max';
33             *max = \&Math::PlanePath::_max;
34              
35 1     1   4 use vars '$VERSION', '@ISA';
  1         2  
  1         61  
36             $VERSION = 18;
37 1     1   1111 use Math::PlanePath;
  1         6453  
  1         135  
38             @ISA = ('Math::PlanePath');
39              
40              
41             # return ($quotient, $remainder)
42             sub _divrem {
43 14     14   20 my ($n, $d) = @_;
44 14 50 33     32 if (ref $n && $n->isa('Math::BigInt')) {
45 0         0 my ($quot,$rem) = $n->copy->bdiv($d);
46 0 0 0     0 if (! ref $d || $d < 1_000_000) {
47 0         0 $rem = $rem->numify; # plain remainder if fits
48             }
49 0         0 return ($quot, $rem);
50             }
51 14         19 my $rem = $n % $d;
52 14         30 return (int(($n-$rem)/$d), # exact division stays in UV
53             $rem);
54             }
55              
56              
57             use Math::PlanePath::Base::Generic
58 1         51 'is_infinite',
59 1     1   6 'round_nearest';
  1         2  
60             use Math::PlanePath::Base::Digits 119 # v.119 for round_up_pow()
61 1         61 'round_up_pow',
62 1     1   758 'round_down_pow';
  1         1609  
63              
64             # uncomment this to run the ### lines
65             # use Smart::Comments;
66              
67 1     1   1205 use Math::PlanePath::ToothpickTree;
  1         3  
  1         158  
68             *new = \&Math::PlanePath::ToothpickTree::new;
69             *x_negative = \&Math::PlanePath::ToothpickTree::x_negative;
70             *y_negative = \&Math::PlanePath::ToothpickTree::y_negative;
71             *rect_to_n_range = \&Math::PlanePath::ToothpickTree::rect_to_n_range;
72             *x_minimum = \&Math::PlanePath::ToothpickTree::x_minimum;
73             *y_minimum = \&Math::PlanePath::ToothpickTree::y_minimum;
74             *sumxy_minimum = \&Math::PlanePath::ToothpickTree::sumxy_minimum;
75             *sumabsxy_minimum = \&Math::PlanePath::ToothpickTree::sumabsxy_minimum;
76             *rsquared_minimum = \&Math::PlanePath::ToothpickTree::rsquared_minimum;
77              
78 1         71 use constant parameter_info_array =>
79             [ { name => 'parts',
80             share_key => 'parts_toothpickreplicate',
81             display => 'Parts',
82             type => 'enum',
83             default => '4',
84             choices => ['4','3','2','1'],
85             choices_display => ['4','3','2','1'],
86             description => 'Which parts of the pattern to generate.',
87             },
88 1     1   7 ];
  1         2  
89              
90 1     1   5 use constant n_start => 0;
  1         2  
  1         42  
91 1     1   4 use constant class_x_negative => 1;
  1         2  
  1         45  
92 1     1   5 use constant class_y_negative => 1;
  1         2  
  1         150  
93              
94             {
95             my @x_negative_at_n = (undef,
96             undef, # 1
97             3, # 2
98             6, # 3
99             5, # 4
100             );
101             sub x_negative_at_n {
102 0     0 1 0 my ($self) = @_;
103 0         0 return $x_negative_at_n[$self->{'parts'}];
104             }
105             }
106             {
107             my @y_negative_at_n = (undef,
108             undef, # 1
109             undef, # 2
110             2, # 3
111             2, # 4
112             );
113             sub y_negative_at_n {
114 0     0 1 0 my ($self) = @_;
115 0         0 return $y_negative_at_n[$self->{'parts'}];
116             }
117             }
118              
119             # parts=1 same as parts=4
120             # parts=2 same as parts=4
121             # parts=3 same as parts=4
122             # parts=4 33,-12
123             # 133,-30
124             # 333,-112
125             # 1333,-230
126             # 3332,-1112 -> 3,-1
127 1     1   5 use constant dir_maximum_dxdy => (3,-1);
  1         1  
  1         1728  
128              
129             #------------------------------------------------------------------------------
130             # Fraction covered
131             # Xlevel = 2^(level+1) - 1
132             # Ylevel = 2^(level+1)
133             # Nend = (2*4^(level+1) + 1)/3 - 1
134             #
135             # Nend / (Xlevel*Ylevel)
136             # -> ((2*4^(level+1) + 1)/3 - 1) / 4^(level+1)
137             # -> (2*4^(level+1) + 1)/3 / 4^(level+1)
138             # -> 2*4^(level+1)/3 / 4^(level+1)
139             # -> 2/3
140              
141             # Leading diagonal 1,3, 7,11,
142             # 23,25,29,43, +22,22,22,32
143             # 87,89,93,97, +86,86,86,86
144             # 109,111,115,171, +86,128
145             # 343
146             # part2start = (4^level + 5)/3 = 3,7,23,87,343
147             # sums of part2start(level), but +2 in second half of each
148             # (3)/3=1
149             # (3+ 1+5)/3=3
150             # (3+ 1+5 + 4+5)/3=9
151              
152              
153             # v v
154             # | -> | part 3
155             # +---h h---+
156             #
157             # +---v h
158             # | -> | part 1 rot then part 3
159             # h +---v
160             #
161             # v v
162             # | -> | part 3 then part 3 again
163             # h---+ +---h
164             #
165              
166             # v +---v
167             # | -> | part 1
168             # +---h h
169             #
170             # v v---+
171             # | -> | part 3 then part 1 rot is +90
172             # h---+ h
173              
174             # N = (2*4^level + 1)/3 + 1 is first of "level"
175             # 3N-3 = 2*4^level + 1
176             # 2*4^level = 3N-4
177             # 4^(level+1) = 6N-8
178             #
179             # part = (2*4^level - 2)/3 many points in "level"
180             # above = (2*4^(level+1) - 2)/3
181             # = (4*2*4^level - 2)/3
182             # = 4*(2*4^level - 2/4)/3
183             # = 4*(2*4^level - 2)/3 + 4*(+ 2 - 2/4)/3
184             # = 4*(2*4^level - 2)/3 + 2
185             # = 4*part + 2
186             # part = (above-2)/4
187              
188             my @quadrant_to_hdx = (1,-1, -1,1);
189             my @quadrant_to_vdy = (1, 1, -1,-1);
190              
191             sub n_to_xy {
192 54     54 1 4131 my ($self, $n) = @_;
193             ### ToothpickReplicate n_to_xy(): $n
194              
195 54 50       112 if ($n < 0) { return; }
  0         0  
196 54 50       120 if (is_infinite($n)) { return ($n,$n); }
  0         0  
197              
198             {
199 54         316 my $int = int($n);
  54         69  
200             ### $int
201             ### $n
202 54 50       99 if ($n != $int) {
203 0         0 my ($x1,$y1) = $self->n_to_xy($int);
204 0         0 my ($x2,$y2) = $self->n_to_xy($int+1);
205 0         0 my $frac = $n - $int; # inherit possible BigFloat
206 0         0 my $dx = $x2-$x1;
207 0         0 my $dy = $y2-$y1;
208 0         0 return ($frac*$dx + $x1, $frac*$dy + $y1);
209             }
210 54         71 $n = $int; # BigFloat int() gives BigInt, use that
211             }
212              
213 54         71 my $parts = $self->{'parts'};
214 54         60 my $x = 0;
215 54         57 my $y = 0;
216 54         57 my $hdx = 1;
217 54         60 my $hdy = 0;
218 54         54 my $vdx = 0;
219 54         58 my $vdy = 1;
220              
221 54 50       164 if ($parts eq '2') {
    100          
    100          
222 0 0       0 if ($n == 0) {
223 0         0 return (0,1);
224             }
225              
226             # first of a replication level
227             # Nlevel = 2*(2*4^level - 2)/3 + 1
228             # = (4*4^level - 4)/3 + 1
229             # = (4*4^level - 4 + 3)/3
230             # = (4*4^level - 1)/3 = 5,21
231             # 3N = 4*4^level - 1
232             # 4^(level+1) = 3N+1
233              
234 0         0 my ($len,$level) = round_down_pow(3*$n+1, 4);
235 0         0 my $three_parts = $len/2;
236              
237             ### $len
238             ### $level
239             ### $three_parts
240             ### start this level: ($len-1)/3
241             ### n reduced: $n-($len-1)/3
242              
243 0         0 (my $quadrant, $n) = _divrem ($n-($len-1)/3, $three_parts);
244             ### $quadrant
245             ### n remainder: $n
246             ### assert: $quadrant >= 0
247             ### assert: $quadrant <= 1
248              
249 0         0 $n += ($len/2-2)/3;
250 0 0       0 if ($quadrant) { $hdx = -1; }
  0         0  
251             ### n in quarter: $n
252              
253             } elsif ($parts == 3) {
254 8 100       16 if ($n <= 1) {
255 2         5 return (0,$n);
256             }
257             # Nend = 3*(2*4^level - 2)/3 + 2
258             # = (2*4^level - 2) + 2
259             # = 2*4^level = 2,8,32
260             # N-1 = 2*4^level
261             # 4^(level+1) = 2N-2
262              
263 6         17 my ($len,$level) = round_down_pow(2*$n, 4);
264 6         62 my $three_parts = $len/2;
265              
266             ### $len
267             ### $level
268             ### $three_parts
269             ### start this level: ($len/2+1)
270             ### n reduced: $n-($len/2+1)
271              
272 6         15 (my $quadrant, $n) = _divrem ($n-$len/2, $three_parts);
273             ### $quadrant
274             ### n remainder: $n
275             ### assert: $quadrant >= 0
276             ### assert: $quadrant <= 2
277              
278 6         12 $n += ($len/2-2)/3;
279             ### n in quarter: $n
280              
281 6 100       21 if ($quadrant == 0) {
    100          
282 2         2 $hdx = 0; # rotate -90
283 2         3 $hdy = -1;
284 2         2 $vdx = 1;
285 2         33 $vdy = 0;
286 2         3 $x = -1; # offset
287             } elsif ($quadrant == 2) {
288 2         4 $hdx = -1; # mirror
289             }
290              
291             } elsif ($parts == 4) {
292 11 100       21 if ($n <= 2) {
293 3 100       7 if ($n == 0) { return (0,0); }
  1         3  
294 2 100       5 if ($n == 1) { return (0,1); }
  1         3  
295 1         3 return (0,-1); # N==2
296             }
297             # first of a replication level
298             # Nlevel = 4*(2*4^level - 2)/3 + 3
299             # = (8*4^level - 8)/3 + 3
300             # = (8*4^level - 8 + 9)/3
301             # = (8*4^level+1)/3 11,43,171
302             # 3N = 8*4^level+1
303             # 8*4^level = 3N-1
304             # 4^(level+2) = 6N-2
305             #
306             # first of this level, using level+2
307             # Nlevel = (4^(level+2)/2+1)/3
308             # = (4^(level+2)+2)/6
309             #
310             # three count = 3*(2*4^level - 2)/3 + 2
311             # = 2*4^level
312             # 43-11 = 32
313             # 172-44 = 128
314              
315             # getting level+2 and len = 4^(level+2)
316 8         25 my ($len,$level) = round_down_pow(6*$n-2, 4);
317 8         77 my $three_parts = $len/8;
318              
319             ### all breakdown ...
320             ### $level
321             ### $len
322             ### $three_parts
323             ### Nlevel base: ($len+2)/6
324              
325 8         21 (my $quadrant, $n) = _divrem ($n-($len+2)/6, $three_parts);
326             ### $quadrant
327             ### n remainder: $n
328             ### assert: $quadrant >= 0
329             ### assert: $quadrant <= 3
330              
331             # quarter middle
332             # Nquarter = (2*4^level - 2)/3 = 2,10,42
333 8         13 $n += ($len/8-2)/3;
334 8         11 $hdx = $quadrant_to_hdx[$quadrant];
335 8         13 $vdy = $quadrant_to_vdy[$quadrant];
336             ### n in quarter: $n
337             }
338              
339             # quarter first of a replication level
340             # Nlevel = 4*(2*4^level - 2)/3 + 2
341             # = (8*4^level - 8)/3 + 2
342             # = (8*4^level - 8 + 6)/3
343             # = (8*4^level - 2)/3 2,10,42
344             # 3N = 8*4^level-2
345             # 8*4^level = 3N+2
346             # 4^(level+2) = 6N+4
347             #
348             # using level+1
349             # Nlevel = (8*4^level - 2)/3
350             # = (2*4^(level+1) - 2)/3
351              
352              
353             # getting level+2 and 16*len
354 49         142 my ($len,$level) = round_down_pow(6*$n+4, 4);
355 49         458 my $part_n = (2*$len-2)/3;
356             ### $level
357             ### $part_n
358              
359 49         54 $len = 2**$level;
360 49         106 for ( ;
361             $level-- >= 0;
362             $len /= 2, $part_n = ($part_n-2)/4) {
363              
364             ### at: "x=$x,y=$y level=$level hxy=$hdx,$hdy vxy=$vdx,$vdy n=$n"
365             ### $len
366             ### $part_n
367             ### assert: $len == 2 ** ($level+1)
368             ### assert: $part_n == (2 * 4 ** ($level+1) - 2)/3
369              
370 145 100       263 if ($n < $part_n) {
371             ### part 0, no change ...
372 55         136 next;
373             }
374              
375 90         100 $n -= $part_n;
376 90         117 $x += $len * ($hdx + $vdx); # diagonal
377 90         104 $y += $len * ($hdy + $vdy);
378              
379 90 100       155 if ($n == 0) {
380             ### toothpick A ...
381 25         32 last;
382             }
383 65 100       114 if ($n == 1) {
384             ### toothpick B ...
385 24         25 $x += $vdx;
386 24         27 $y += $vdy;
387 24         26 last;
388             }
389 41         49 $n -= 2;
390              
391 41 100       73 if ($n < $part_n) {
392             ### part 1, rotate ...
393 16         18 $x -= $hdx; # offset
394 16         17 $y -= $hdy;
395 16         29 ($hdx,$hdy, $vdx,$vdy) # rotate 90 in direction v toward h
396             = (-$vdx,-$vdy, $hdx,$hdy);
397 16         41 next;
398             }
399 25         26 $n -= $part_n;
400              
401 25 100       47 if ($n < $part_n) {
402             ### part 2 ...
403 9         22 next;
404             }
405 16         16 $n -= $part_n;
406              
407             ### part 3, mirror ...
408 16         18 $hdx = -$hdx;
409 16         54 $hdy = -$hdy;
410             }
411              
412             ### assert: $n == 0 || $n == 1
413              
414             ### final: "x=$x y=$y"
415 49         111 return ($x,$y);
416             }
417              
418             sub xy_to_n {
419 0     0 1 0 my ($self, $x, $y) = @_;
420             ### ToothpickReplicate xy_to_n(): "$x, $y"
421              
422 0         0 $x = round_nearest ($x);
423 0         0 $y = round_nearest ($y);
424              
425 0         0 my $parts = $self->{'parts'};
426 0   0     0 my $rotated = ($parts == 3 && $x >= 0 && $y < 0);
427 0 0       0 if ($rotated) {
428 0         0 ($x,$y) = (-$y,$x+1); # rotate +90 and shift up
429             ### rotated: "x=$x y=$y"
430             }
431              
432 0         0 my ($len,$level) = round_down_pow (max(abs($x), abs($y)-1),
433             2);
434 0 0       0 if (is_infinite($level)) {
435 0         0 return $level;
436             }
437             ### $level
438             ### $len
439              
440 0         0 my $zero = $x * 0 * $y;
441 0         0 my $n = $zero;
442              
443 0 0       0 if ($parts == 2) {
    0          
    0          
444 0 0       0 if ($x == 0) {
445 0 0       0 if ($y == 1) { return 0; }
  0         0  
446             }
447 0         0 $n += (2*$len*$len+1)/3; # +1,+3,+11,+43
448 0 0       0 if ($x < 0) {
449 0         0 $x = -$x;
450 0         0 $n += 2*$len*$len; # second quad, +2,+8,+32
451             }
452              
453             } elsif ($parts == 3) {
454             ### 3/4 ...
455 0 0       0 if ($x == 0) {
456 0 0       0 if ($y == 0) { return 0; }
  0         0  
457 0 0       0 if ($y == 1) { return 1; }
  0         0  
458             }
459 0         0 $n += (10*$len*$len+2)/3; # +4,+14,+54,+214,+854,+3414
460 0 0       0 if ($rotated) {
    0          
461 0         0 $n -= 2*$len*$len; # fourth quad, -2, -8, -32
462             } elsif ($x < 0) {
463 0         0 $x = -$x;
464 0 0       0 if ($y > 0) {
465 0         0 $n += 2*$len*$len; # second quad, +2, +8, +32
466             } else {
467 0         0 return undef; # third quad, empty
468             }
469             }
470             } elsif ($parts == 4) {
471 0 0       0 if ($x == 0) {
472 0 0       0 if ($y == 0) { return 0; }
  0         0  
473 0 0       0 if ($y == 1) { return 1; }
  0         0  
474 0 0       0 if ($y == -1) { return 2; }
  0         0  
475             }
476 0         0 $n += (2*$len*$len+1);
477 0 0       0 if ($x < 0) {
478 0         0 $x = -$x;
479 0 0       0 if ($y > 0) {
480 0         0 $n += 2*$len*$len; # second quad, +2, +8, +32
481             } else {
482 0         0 $n += 4*$len*$len; # third quad, +4,+16
483 0         0 $y = -$y;
484             }
485             } else {
486 0 0       0 if ($y < 0) {
487 0         0 $n += 6*$len*$len; # fourth quad
488 0         0 $y = -$y;
489             }
490             }
491             }
492              
493             # 2^(level+1)-1
494             # v
495             # +-----------+---------+
496             # | | | <- 2^(level+1)
497             # | 3 2 |
498             # | mirror same |
499             # | --B-- | <- 2^level + 1
500             # | | |
501             # +-- A --+ <- 2^level
502             # | |
503             # 1 |
504             # rot |
505             # 0 +90 |
506             # | |
507             # +-----------+
508             # ^
509             # 2^level
510              
511 0         0 my $part_n = (2*$len*$len - 2) / 3;
512             ### $part_n
513              
514 0         0 while ($level-- > 0) {
515             ### at: "x=$x,y=$y len=$len part_n=$part_n n=$n"
516             ### assert: $len == 2 ** ($level+1)
517             ### assert: $part_n == (2 * 4 ** ($level+1) - 2)/3
518              
519 0 0       0 if ($x == $len) {
520 0 0       0 if ($y == $len) {
521             ### toothpick A ...
522 0         0 return $n + $part_n;
523             }
524 0 0       0 if ($y == $len+1) {
525             ### toothpick B ...
526 0         0 return $n + $part_n + 1;
527             }
528             }
529              
530 0 0       0 if ($y <= $len) {
531 0 0       0 if ($x < $len) {
532             ### part 0 ...
533             } else {
534             ### part 1, rotate ...
535 0         0 $n += $part_n + 2;
536 0         0 ($x,$y) = ($len-$y,$x-$len+1); # shift, rotate +90
537             }
538             } else {
539 0         0 $y -= $len;
540 0 0       0 if ($x > $len) {
541             ### part 2 ...
542 0         0 $n += 2*$part_n + 2;
543 0         0 $x -= $len;
544             } else {
545             ### part 3 ...
546 0         0 $n += 3*$part_n + 2;
547 0         0 $x = $len-$x; # mirror
548             }
549             }
550              
551 0         0 $len /= 2;
552 0         0 $part_n = ($part_n-2)/4;
553             }
554              
555             ### end loop: "x=$x y=$y n=$n"
556              
557 0 0       0 if ($x == 1) {
558 0 0       0 if ($y == 1) {
    0          
559 0         0 return $n;
560             } elsif ($y == 2) {
561 0         0 return $n + 1;
562             }
563             }
564              
565 0         0 return undef;
566             }
567              
568             #------------------------------------------------------------------------------
569             # levels
570              
571             # parts=1
572             # LevelPoints[k] = 4*LevelPoints[k] + 2 starting LevelPoints[0] = 2
573             # LevelPoints[k] = 2 + 2*4 + 2*4^2 + ... + 2*4^(k-1) + 4^k*LevelPoints[0]
574             # LevelPoints[k] = 2 + 2*4 + 2*4^2 + ... + 2*4^(k-1) + 2*4^k
575             # LevelPoints[k] = 2*(4^(k+1) - 1)/3
576              
577             {
578             my %level_to_n_range = (4 => -2,
579             3 => -3,
580             2 => -4,
581             1 => -5,
582             );
583             sub level_to_n_range {
584 9     9 1 402 my ($self, $level) = @_;
585             return (0,
586             (4**($level+1) * (2*$self->{'parts'})
587 9         37 + $level_to_n_range{$self->{'parts'}}) / 3);
588             }
589             }
590             {
591             # $level_to_n_range{} and _divrem_mutate() rounded up
592             my %n_to_level = (4 => 2 + 2*4-1,
593             3 => 3 + 2*3-1,
594             2 => 4 + 2*2-1,
595             1 => 5 + 2-1,
596             );
597             sub n_to_level {
598 0     0 1   my ($self, $n) = @_;
599 0 0         if ($n < 0) { return undef; }
  0            
600 0 0         if (is_infinite($n)) { return $n; }
  0            
601 0           $n = round_nearest($n);
602 0           $n *= 3;
603 0           $n += $n_to_level{$self->{'parts'}};
604 0           _divrem_mutate ($n, 2*$self->{'parts'});
605 0           my ($pow, $exp) = round_down_pow ($n-1, 4);
606 0           return $exp;
607             }
608             }
609              
610             # return $remainder, modify $n
611             # the scalar $_[0] is modified, but if it's a BigInt then a new BigInt is made
612             # and stored there, the bigint value is not changed
613             sub _divrem_mutate {
614 0     0     my $d = $_[1];
615 0           my $rem;
616 0 0 0       if (ref $_[0] && $_[0]->isa('Math::BigInt')) {
617 0           ($_[0], $rem) = $_[0]->copy->bdiv($d); # quot,rem in array context
618 0 0 0       if (! ref $d || $d < 1_000_000) {
619 0           return $rem->numify; # plain remainder if fits
620             }
621             } else {
622 0           $rem = $_[0] % $d;
623 0           $_[0] = int(($_[0]-$rem)/$d); # exact division stays in UV
624             }
625 0           return $rem;
626             }
627              
628             #------------------------------------------------------------------------------
629             1;
630             __END__