File Coverage

blib/lib/Math/PlanePath/UlamWarburtonQuarter.pm
Criterion Covered Total %
statement 176 237 74.2
branch 67 112 59.8
condition 16 26 61.5
subroutine 21 28 75.0
pod 15 15 100.0
total 295 418 70.5


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             package Math::PlanePath::UlamWarburtonQuarter;
20 2     2   1410 use 5.004;
  2         8  
21 2     2   10 use strict;
  2         5  
  2         46  
22 2     2   9 use Carp 'croak';
  2         4  
  2         89  
23 2     2   12 use List::Util 'sum';
  2         3  
  2         198  
24              
25 2     2   14 use vars '$VERSION', '@ISA';
  2         2  
  2         121  
26             $VERSION = 127;
27 2     2   767 use Math::PlanePath;
  2         5  
  2         118  
28             @ISA = ('Math::PlanePath');
29             *_divrem_mutate = \&Math::PlanePath::_divrem_mutate;
30              
31             use Math::PlanePath::Base::Generic
32 2         105 'is_infinite',
33 2     2   14 'round_nearest';
  2         4  
34             use Math::PlanePath::Base::Digits
35 2         192 'round_down_pow',
36             'bit_split_lowtohigh',
37             'digit_split_lowtohigh',
38 2     2   548 'digit_join_lowtohigh';
  2         3  
39              
40             # uncomment this to run the ### lines
41             # use Smart::Comments;
42              
43              
44 2         132 use constant parameter_info_array =>
45             [
46             { name => 'parts',
47             share_key => 'parts_ulamwarburton_quarter',
48             display => 'Parts',
49             type => 'enum',
50             default => '1',
51             choices => ['1','octant','octant_up' ],
52             choices_display => ['1','Octant','Octant Up' ],
53             description => 'Which parts of the plane to fill.',
54             },
55             Math::PlanePath::Base::Generic::parameter_info_nstart1(),
56 2     2   15 ];
  2         5  
57              
58 2     2   14 use constant class_x_negative => 0;
  2         4  
  2         93  
59 2     2   12 use constant class_y_negative => 0;
  2         3  
  2         2319  
60              
61             sub diffxy_minimum {
62 0     0 1 0 my ($self) = @_;
63 0 0       0 return ($self->{'parts'} eq 'octant' ? 0 : undef);
64             }
65             sub diffxy_maximum {
66 0     0 1 0 my ($self) = @_;
67 0 0       0 return ($self->{'parts'} eq 'octant_up' ? 0 : undef);
68             }
69              
70             # Minimum dir=0 at N=13 dX=2,dY=0.
71             # Maximum dir seems dX=13,dY=-9 at N=149 going top-left part to new bottom
72             # right diagonal.
73             my %dir_maximum_dxdy = (1 => [13,-9],
74             octant => [1,-1], # South-East
75             octant_up => [0,-1], # South
76             );
77             sub dir_maximum_dxdy {
78 0     0 1 0 my ($self) = @_;
79 0         0 return @{$dir_maximum_dxdy{$self->{'parts'}}};
  0         0  
80             }
81              
82             sub tree_num_children_list {
83 0     0 1 0 my ($self) = @_;
84 0 0       0 return ($self->{'parts'} =~ /octant/
85             ? (0, 1, 2, 3)
86             : (0, 1, 3));
87             }
88              
89             #------------------------------------------------------------------------------
90             sub new {
91 12     12 1 1450 my $self = shift->SUPER::new(@_);
92 12 100       43 if (! defined $self->{'n_start'}) {
93 8         37 $self->{'n_start'} = $self->default_n_start;
94             }
95 12   50     55 my $parts = ($self->{'parts'} ||= '1');
96 12 50       32 if (! exists $dir_maximum_dxdy{$parts}) {
97 0         0 croak "Unrecognised parts option: ", $parts;
98             }
99 12         28 return $self;
100             }
101              
102             # 7 7 7 7
103             # 6 6
104             # 7 5 5 7
105             # 4
106             # 3 3 5 7
107             # 2 6
108             # 1 3 7 7
109             #
110             # 1+1+3=5
111             # 5+1+3*5=21
112             # 1+3 = 4
113             # 1+3+3+9 = 16
114             #
115             # 0
116             # 1 0 +1
117             # 2 1 +1 <- 1
118             # 3 2 +3
119             # 4 5 +1 <- 1 + 4 = 5
120             # 5 6 +3
121             # 6 9 +3
122             # 7 12 +9
123             # 8 21 <- 1 + 4 + 16 = 21
124              
125             # 1+3 = 4 power 2
126             # 1+3+3+9 = 16 power 3
127             # 1+3+3+9+3+9+9+27 = 64 power 4
128             #
129             # (1+4+16+...+4^(l-1)) = (4^l-1)/3
130             # l=1 total=(4-1)/3 = 1
131             # l=2 total=(16-1)/3 = 5
132             # l=3 total=(64-1)/3=63/3 = 21
133             #
134             # n = 1 + (4^l-1)/3
135             # n-1 = (4^l-1)/3
136             # 3n-3 = (4^l-1)
137             # 3n-2 = 4^l
138             #
139             # 3^0+3^1+3^1+3^2 = 1+3+3+9=16
140             # x+3x+3x+9x = 16x = 256
141             #
142             # 22
143             # 20 19 18 17
144             # 12 11
145             # 21 9 8 16
146             # 6
147             # 5 4 7 15
148             # 2 10
149             # 1 3 13 14
150             #
151              
152             sub n_to_xy {
153 117     117 1 6819 my ($self, $n) = @_;
154             ### UlamWarburtonQuarter n_to_xy(): $n
155              
156 117 50       287 if ($n < $self->{'n_start'}) { return; }
  0         0  
157 117 50       317 if (is_infinite($n)) { return ($n,$n); }
  0         0  
158              
159             {
160 117         208 my $int = int($n);
  117         169  
161             ### $int
162             ### $n
163 117 50       213 if ($n != $int) {
164 0         0 my ($x1,$y1) = $self->n_to_xy($int);
165 0         0 my ($x2,$y2) = $self->n_to_xy($int+1);
166 0         0 my $frac = $n - $int; # inherit possible BigFloat
167 0         0 my $dx = $x2-$x1;
168 0         0 my $dy = $y2-$y1;
169 0         0 return ($frac*$dx + $x1, $frac*$dy + $y1);
170             }
171 117         188 $n = $int; # BigFloat int() gives BigInt, use that
172             }
173              
174 117         190 $n = $n - $self->{'n_start'} + 1; # N=1 basis
175 117 100       226 if ($n == 1) { return (0,0); }
  5         11  
176              
177 112 50       276 my ($depthsum, $nrem, $rowwidth) = _n1_to_depthsum_rem_width($self,$n)
178             or return ($n,$n); # N==nan or N==+inf
179              
180             ### assert: $nrem >= 0
181             ### assert: $nrem < $width
182 112 50       271 if ($self->{'parts'} eq 'octant_up') {
183 0         0 $nrem += ($rowwidth-1)/2;
184             ### assert: $nrem < $width
185             }
186              
187 112         261 my @ndigits = digit_split_lowtohigh($nrem,3);
188 112         234 my $dhigh = shift(@$depthsum) - 1; # highest term
189 112         163 my $x = 0;
190 112         169 my $y = 0;
191 112         212 foreach my $depthsum (reverse @$depthsum) { # depth terms low to high
192 257         359 my $ndigit = shift @ndigits; # N digits low to high
193             ### $depthsum
194             ### $ndigit
195              
196 257         404 $x += $depthsum;
197 257         336 $y += $depthsum;
198             ### depthsum to xy: "$x,$y"
199              
200 257 100       395 if ($ndigit) {
201 176 100       320 if ($ndigit == 2) {
202 88         184 ($x,$y) = (-$y,$x); # rotate +90
203             }
204             } else {
205             # digit==0 (or undef when run out of @ndigits)
206 81         156 ($x,$y) = ($y,-$x); # rotate -90
207             }
208             ### rotate to: "$x,$y"
209             }
210              
211             ### final: "$x,$y"
212 112         303 return ($dhigh + $x, $dhigh + $y);
213             }
214              
215             sub xy_to_n {
216 560     560 1 8124 my ($self, $x, $y) = @_;
217             ### UlamWarburtonQuarter xy_to_n(): "$x, $y"
218              
219 560         1042 $x = round_nearest ($x);
220 560         1076 $y = round_nearest ($y);
221 560         906 my $parts = $self->{'parts'};
222 560 50 100     1774 if ($y < 0
    100 33        
      66        
223             || $x < ($parts eq 'octant' ? $y : 0)
224             || ($parts eq 'octant_up' && $x > $y)) {
225 329         593 return undef;
226             }
227 231 100 100     496 if ($x == 0 && $y == 0) {
228 6         20 return $self->{'n_start'};
229             }
230 225         296 $x += 1; # pushed away by 1 ...
231 225         297 $y += 1;
232              
233 225         529 my ($len, $exp) = round_down_pow ($x + $y, 2);
234 225 50       472 if (is_infinite($exp)) { return $exp; }
  0         0  
235              
236 225         440 my $depth
237             = my $n
238             = ($x * 0 * $y); # inherit bignum 0
239 225         323 my $rowwidth = $depth + 1;
240              
241 225         410 while ($exp-- >= 0) {
242             ### at: "$x,$y n=$n len=$len"
243              
244             # first quadrant square
245             ### assert: $x >= 0
246             ### assert: $y >= 0
247             # ### assert: $x < 2*$len
248             # ### assert: $y < 2*$len
249              
250 857 100 100     1945 if ($x >= $len || $y >= $len) {
251             # one of three quarters away from origin
252             # +---+---+
253             # | 2 | 1 |
254             # +---+---+
255             # | | 0 |
256             # +---+---+
257              
258 545         747 $x -= $len;
259 545         689 $y -= $len;
260             ### shift to: "$x,$y"
261              
262 545 100       821 if ($x) {
263 366 100       646 unless ($y) {
264 46         105 return undef; # x==0, y!=0, nothing
265             }
266             } else {
267 179 100       321 if ($y) {
268 45         105 return undef; # x!=0, y-=0, nothing
269             }
270             }
271              
272 454         588 $depth += $len;
273 454 100 66     948 if ($x || $y) {
274 320         421 $rowwidth *= 3;
275 320         404 $n *= 3;
276 320 100       566 if ($y < 0) {
    100          
277             ### bottom right, digit 0 ...
278 101         201 ($x,$y) = (-$y,$x); # rotate +90
279             } elsif ($x >= 0) {
280             ### top right, digit 1 ...
281 104         143 $n += 1;
282             } else {
283             ### top left, digit 2 ...
284 115         196 ($x,$y) = ($y,-$x); # rotate -90
285 115         162 $n += 2;
286             }
287             }
288             }
289              
290 766         1398 $len /= 2;
291             }
292              
293             ### $n
294             ### $depth
295              
296 134 50       268 if ($self->{'parts'} eq 'octant_up') {
297 0         0 $n -= ($rowwidth-1)/2;
298             }
299              
300 134         279 return $n + $self->tree_depth_to_n($depth-1);
301             }
302              
303             # not exact
304             sub rect_to_n_range {
305 50     50 1 4148 my ($self, $x1,$y1, $x2,$y2) = @_;
306             ### UlamWarburtonQuarter rect_to_n_range(): "$x1,$y1 $x2,$y2"
307              
308 50         129 $x1 = round_nearest ($x1);
309 50         100 $y1 = round_nearest ($y1);
310 50         104 $x2 = round_nearest ($x2);
311 50         99 $y2 = round_nearest ($y2);
312              
313 50 50       119 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
314 50 50       87 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
315              
316 50 50 33     191 if ($x2 < 0 || $y2 < 0) {
317 0         0 return (1, 0); # all outside first quadrant
318             }
319              
320 50 50       86 if ($x1 < 0) { $x1 *= 0; }
  0         0  
321 50 50       97 if ($y1 < 0) { $y1 *= 0; }
  0         0  
322              
323             # level numbers
324 50 100       102 my $dlo = ($x1 > $y1 ? $x1 : $y1)+1;
325 50 100       106 my $dhi = ($x2 > $y2 ? $x2 : $y2);
326             ### $dlo
327             ### $dhi
328              
329             # round down to level=2^k numbers
330 50 50       99 if ($dlo) {
331 50         134 ($dlo) = round_down_pow ($dlo,2);
332             }
333 50         109 ($dhi) = round_down_pow ($dhi,2);
334              
335             ### rounded to pow2: "$dlo ".(2*$dhi)
336              
337 50         147 return ($self->tree_depth_to_n($dlo-1),
338             $self->tree_depth_to_n(2*$dhi-1));
339             }
340              
341             #------------------------------------------------------------------------------
342 2     2   16 use constant tree_num_roots => 1;
  2         4  
  2         2198  
343              
344             # ENHANCE-ME: step by the bits, not by X,Y
345             sub tree_n_children {
346 9     9 1 712 my ($self, $n) = @_;
347 9 50       26 if ($n < $self->{'n_start'}) {
348 0         0 return;
349             }
350 9         24 my ($x,$y) = $self->n_to_xy($n);
351 9         23 my @ret;
352 9         16 my $dx = 1;
353 9         11 my $dy = 1;
354 9         20 foreach (1 .. 4) {
355 36 100       72 if (defined (my $n_child = $self->xy_to_n($x+$dx,$y+$dy))) {
356 19 100       38 if ($n_child > $n) {
357 11         20 push @ret, $n_child;
358             }
359             }
360 36         80 ($dx,$dy) = (-$dy,$dx); # rotate +90
361             }
362 9         40 return sort {$a<=>$b} @ret;
  4         17  
363             }
364             sub tree_n_parent {
365 12     12 1 912 my ($self, $n) = @_;
366 12 100       38 if ($n <= $self->{'n_start'}) {
367 1         3 return undef;
368             }
369 11         26 my ($x,$y) = $self->n_to_xy($n);
370 11         19 my $dx = 1;
371 11         17 my $dy = 1;
372 11         25 foreach (1 .. 4) {
373 33 100       77 if (defined (my $n_parent = $self->xy_to_n($x+$dx,$y+$dy))) {
374 24 100       45 if ($n_parent < $n) {
375 11         29 return $n_parent;
376             }
377             }
378 22         53 ($dx,$dy) = (-$dy,$dx); # rotate +90
379             }
380 0         0 return undef;
381             }
382              
383             # level = depth+1 = 2^a + 2^b + 2^c + 2^d ... a>b>c>d...
384             # Ndepth = 1 + (-1
385             # + 4^a
386             # + 3 * 4^b
387             # + 3^2 * 4^c
388             # + 3^3 * 4^d + ...) / 3
389             sub tree_depth_to_n {
390 895     895 1 2047 my ($self, $depth) = @_;
391             ### tree_depth_to_n(): $depth
392 895 50       2014 if (is_infinite($depth)) {
393 0         0 return $depth;
394             }
395 895 50       2000 unless ($depth >= 0) {
396 0         0 return undef;
397             }
398 895         1398 my $n = $depth*0; # inherit bignum 0
399 895         1380 my $pow3 = 1 + $n; # inherit bignum 1
400 895         2084 foreach my $bit (reverse bit_split_lowtohigh($depth+1)) { # high to low
401 3831         5262 $n *= 4;
402 3831 100       6674 if ($bit) {
403 1990         2656 $n += $pow3;
404 1990         3019 $pow3 *= 3;
405             }
406             }
407 895 50       2307 if ($self->{'parts'} =~ /octant/) {
408 0         0 $n = ($n + (3*$depth-1))/6;
409             } else {
410 895         1491 $n = ($n-1)/3;
411             }
412 895         2330 return $n + $self->{'n_start'};
413             }
414              
415             sub tree_n_to_depth {
416 0     0 1 0 my ($self, $n) = @_;
417              
418 0         0 $n = int($n - $self->{'n_start'} + 1); # N=1 basis
419 0 0       0 if ($n < 1) {
420 0         0 return undef;
421             }
422 0 0       0 (my $depthsum, $n) = _n1_to_depthsum_rem_width($self,$n)
423             or return $n; # N==nan or N==+infinity
424 0         0 return sum(-1, @$depthsum);
425             }
426              
427             # Return ($aref, $remaining_n).
428             # sum(@$aref) = depth starting depth=1
429             #
430             # depth+1 = 2^k
431             # Ndepth(depth) = (4^k+2)/3
432             # 3N-2 = 4^k
433             # NdepthOct(depth) = ((4^k+2)/3 + 2^k)/2
434             # 6N-2 = 4^k + 3*2^k
435             #
436             sub _n1_to_depthsum_rem_width {
437 112     112   194 my ($self, $n) = @_;
438             ### _n1_to_depthsum_rem_width(): $n
439              
440 112         206 my $octant = ($self->{'parts'} =~ /octant/);
441 112 50       325 my ($power, $exp) = round_down_pow (($octant ? 6 : 3)*$n - 2, 4);
442 112 50       243 if (is_infinite($exp)) {
443 0         0 return;
444             }
445              
446             ### $power
447             ### $exp
448             ### pow base: ($power - 1)/3 + 1
449              
450             {
451 112         189 my $sub = ($power + 2)/3; # (power-1)/3 + 1
  112         202  
452 112 50       187 if ($octant) {
453 0         0 $sub = ($sub + 2**$exp) / 2;
454             ### prospective sub: $sub
455             ### assert: $sub == ($power + 3 * 2 ** $exp + 2)/6
456              
457 0 0       0 if ($sub > $n) {
458 0         0 $exp -= 1;
459 0         0 $power /= 4;
460 0         0 $sub = ($power + 3*2**$exp + 2)/6;
461             }
462             }
463             ### assert: $sub <= $n
464 112         189 $n -= $sub;
465             }
466             ### n less pow base: $n
467              
468 112         288 my @depthsum = (2**$exp);
469              
470             # find the cumulative levelpoints total <= $n, being the start of the
471             # level containing $n
472             #
473 112         171 my $factor = 1;
474 112         216 while (--$exp >= 0) {
475 387         527 $power /= 4;
476 387         534 my $sub = $power * $factor;
477 387 50       664 if ($octant) {
478 0         0 $sub = ($sub + 2**$exp)/2;
479             }
480             ### $sub
481 387         535 my $rem = $n - $sub;
482              
483             ### $n
484             ### $power
485             ### $factor
486             ### consider subtract: $sub
487             ### $rem
488              
489 387 100       751 if ($rem >= 0) {
490 257         358 $n = $rem;
491 257         399 push @depthsum, 2**$exp;
492 257         463 $factor *= 3;
493             }
494             }
495              
496             ### _n1_to_depthsum_rem_width() result ...
497             ### @depthsum
498             ### remaining n: $n
499             ### assert: $n >= 0
500             ### assert: $n < $factor
501              
502 112         377 return (\@depthsum, $n, $factor);
503             }
504              
505              
506             # at 0,2 turn and new height limit
507             # at 1 keep existing depth limit
508             # N=30 rem=1 = 0,1 depth=11=8+2+1=1011 width=9
509             #
510             sub tree_n_to_subheight {
511 0     0 1 0 my ($self, $n) = @_;
512             ### tree_n_to_subheight(): $n
513              
514 0         0 $n = int($n - $self->{'n_start'} + 1); # N=1 basis
515 0 0       0 if ($n < 1) {
516 0         0 return undef;
517             }
518 0 0       0 my ($depthsum, $nrem, $rowwidth) = _n1_to_depthsum_rem_width($self,$n)
519             or return $n; # N==nan or N==+infinity
520             ### $depthsum
521             ### $nrem
522              
523 0 0       0 if ($self->{'parts'} eq 'octant_up') {
524 0         0 $nrem += ($rowwidth-1)/2;
525             }
526              
527 0         0 my $sub = pop @$depthsum;
528 0   0     0 while (@$depthsum && _divrem_mutate($nrem,3) == 1) {
529 0         0 $sub += pop @$depthsum;
530             }
531 0 0       0 if (@$depthsum) {
532 0         0 return $depthsum->[-1] - 1 - $sub;
533             } else {
534 0         0 return undef; # $nrem all 1-digits
535             }
536             }
537              
538             #------------------------------------------------------------------------------
539             # levels
540              
541             sub level_to_n_range {
542 27     27 1 1685 my ($self, $level) = @_;
543 27         151 return ($self->{'n_start'},
544             $self->tree_depth_to_n_end(2**($level+1) - 2));
545             }
546             sub n_to_level {
547 0     0 1   my ($self, $n) = @_;
548 0           my $depth = $self->tree_n_to_depth($n);
549 0 0         if (! defined $depth) { return undef; }
  0            
550 0           my ($pow, $exp) = round_down_pow ($depth+1, 2);
551 0           return $exp;
552             }
553              
554             #------------------------------------------------------------------------------
555             1;
556             __END__