File Coverage

blib/lib/Math/PlanePath/HTree.pm
Criterion Covered Total %
statement 34 213 15.9
branch 0 80 0.0
condition 0 33 0.0
subroutine 12 24 50.0
pod 12 12 100.0
total 58 362 16.0


line stmt bran cond sub pod time code
1             # Copyright 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              
19             package Math::PlanePath::HTree;
20 1     1   1268 use 5.004;
  1         3  
21 1     1   5 use strict;
  1         2  
  1         45  
22             #use List::Util 'max';
23             *max = \&Math::PlanePath::_max;
24              
25             use Math::PlanePath::Base::Generic
26 1         71 'is_infinite',
27 1     1   768 'round_nearest';
  1         828  
28             use Math::PlanePath::Base::Digits 119 # v.119 for round_up_pow()
29 1         70 'round_up_pow',
30             'round_down_pow',
31             'bit_split_lowtohigh',
32 1     1   766 'digit_join_lowtohigh';
  1         1624  
33              
34 1     1   846 use Math::PlanePath::LCornerTree;
  1         5  
  1         55  
35             *_divrem = \&Math::PlanePath::LCornerTree::_divrem;
36              
37 1     1   7 use vars '$VERSION', '@ISA';
  1         2  
  1         53  
38             $VERSION = 18;
39 1     1   5 use Math::PlanePath;
  1         2  
  1         150  
40             @ISA = ('Math::PlanePath');
41              
42             # uncomment this to run the ### lines
43             # use Smart::Comments;
44              
45              
46             # return $remainder, modify $n
47             # the scalar $_[0] is modified, but if it's a BigInt then a new BigInt is made
48             # and stored there, the bigint value is not changed
49             sub _divrem_mutate {
50 0     0   0 my $d = $_[1];
51 0         0 my $rem;
52 0 0 0     0 if (ref $_[0] && $_[0]->isa('Math::BigInt')) {
53 0         0 ($_[0], $rem) = $_[0]->copy->bdiv($d); # quot,rem in array context
54 0 0 0     0 if (! ref $d || $d < 1_000_000) {
55 0         0 return $rem->numify; # plain remainder if fits
56             }
57             } else {
58 0         0 $rem = $_[0] % $d;
59 0         0 $_[0] = int(($_[0]-$rem)/$d); # exact division stays in UV
60             }
61 0         0 return $rem;
62             }
63              
64              
65 1     1   5 use constant n_start => 1;
  1         1  
  1         53  
66 1     1   5 use constant class_x_negative => 0;
  1         1  
  1         52  
67 1     1   5 use constant class_y_negative => 0;
  1         2  
  1         42  
68 1     1   5 use constant tree_num_children_list => (0,1,2);
  1         1  
  1         1660  
69              
70             sub n_to_xy {
71 0     0 1 0 my ($self, $n) = @_;
72             ### HTree n_to_xy(): $n
73              
74 0 0       0 if ($n < 1) { return; }
  0         0  
75 0 0       0 if (is_infinite($n)) { return ($n,$n); }
  0         0  
76             {
77 0         0 my $int = int($n);
  0         0  
78             ### $int
79             ### $n
80 0 0       0 if ($n != $int) {
81 0         0 my ($x1,$y1) = $self->n_to_xy($int);
82 0         0 my ($x2,$y2) = $self->n_to_xy($int+1);
83 0         0 my $frac = $n - $int; # inherit possible BigFloat
84 0         0 my $dx = $x2-$x1;
85 0         0 my $dy = $y2-$y1;
86 0         0 return ($frac*$dx + $x1, $frac*$dy + $y1);
87             }
88 0         0 $n = $int; # BigFloat int() gives BigInt, use that
89             }
90              
91 0         0 my $zero = $n * 0;
92 0         0 my @nbits = bit_split_lowtohigh($n);
93 0         0 my $len = (2 + $zero) ** int($#nbits/2); # half power
94             ### $len
95              
96             # eg. N=9 "up" block, nbits even, $#nbits odd, dX=1 dY=0 East
97             # High 1-bit of sub-tree position does rotate +90 to North for
98             # initial step North N=8 to N=9.
99             #
100             # eg. N=17 "right" block, nbits odd, $#nbits even, dX=0 dY=-1 South
101             # High 1-bit of sub-tree position does rotate +90 to East for
102             # initial step East N=16 to N=17.
103             #
104 0         0 my $dx = ($#nbits % 2);
105 0         0 my $dy = $dx - 1;
106             ### initial direction: "$dx, $dy"
107              
108 0         0 my $x = my $y = $len;
109 0 0       0 if ($dx) {
110 0         0 $y *= 2;
111             }
112 0         0 $x -= 1;
113 0         0 $y -= 1;
114             ### initial xy: "$x, $y"
115              
116             ### assert: $nbits[-1] == 1
117 0         0 pop @nbits; # strip high 1-bit which is N=2^k spine bit
118             ### strip high spine 1-bit to: @nbits
119             # N=10001xxx
120             # ^^^^^^^ leaving these
121              
122             # Strip high 0-bits of sub-tree.
123             # Could have @nbits all 0-bits if N=2^k spine point.
124 0   0     0 while (@nbits && ! $nbits[-1]) {
125 0         0 pop @nbits;
126             }
127             ### strip high zeros to: @nbits
128             # N=10001xxx
129             # ^^^^ leaving these
130             # or if an N=1000000 spine point then @nbits now empty and no move $x,$y.
131              
132 0         0 foreach my $bit (reverse @nbits) { # high to low
133             ### at: "$x,$y bit=$bit len=$len dir=$dx,$dy"
134              
135 0         0 ($dx,$dy) = ($dy,-$dx); # rotate -90 for $bit==0
136 0 0       0 if ($bit) {
137 0         0 $dx = -$dx; # rotate 180 to give rotate +90 for $bit==1
138 0         0 $dy = -$dy;
139             }
140             ### turn to: "dir=$dx,$dy"
141              
142 0 0       0 if ($dx) { $len /= 2; } # halve when going horizontal
  0         0  
143 0         0 $x += $dx * $len;
144 0         0 $y += $dy * $len;
145             }
146              
147             ### return: "$x,$y"
148 0         0 return ($x,$y);
149             }
150              
151             # | [37] | 34=10,0010
152             # 13 | 4,13--5,13--6,13 35=10,0011
153             # | | | |
154             # 12 | | 36=10,0100
155             # | [33] | 37=10,0101
156             # 11 [35]1,7---------3,7----------5,7[34]
157             # | | |
158             # 10 0,10 | | 4,6 |
159             # | | | | | | |
160             # 9 0,9-----1,9---2,9 | 4,5---5,5---6,5
161             # | | | | [36] |
162             # 8 0,8 | 4,4
163             # |
164             # 7 [32]3,7-------------------------------
165             # |
166             # 6 0,6 2,6 | [30] [29] 17=1,0001
167             # | [9] | | | [19] | 18=1,0010
168             # 5 0,5------1,5---2,5 | [23]---5,5---[22] 20=1,0100
169             # | | | | | | | 24=1,1000
170             # 4 0,4 | 2,4 | [31] | [28]
171             # [15] | | |
172             # 3 [8]1,3---------3,3-----------5,3[17]
173             # | [16] |
174             # 2 0,2[3] | 2,2[7] [24] | [27]
175             # | | | | | |
176             # 1 0,1[2]---1,1---2,1[5] [20]4,1---5,1---6,1[21] 21=1,0101
177             # | 4 | | [18] |
178             # 0 0,0[1] 2,0[6] [25] [26]
179             #
180             # 0 1 2 3 4 5 6
181              
182             sub xy_to_n {
183 0     0 1 0 my ($self, $x, $y) = @_;
184             ### HTree xy_to_n(): "$x,$y"
185              
186 0         0 $x = round_nearest ($x);
187 0         0 $y = round_nearest ($y);
188 0 0       0 if (is_infinite($x)) {
189 0         0 return $x;
190             }
191 0 0       0 if (is_infinite($y)) {
192 0         0 return $y;
193             }
194              
195 0         0 my ($len,$exp);
196 0         0 my $n;
197 0         0 my ($xlen,$xexp) = round_down_pow($x+1, 2);
198 0         0 my ($ylen,$yexp) = round_down_pow($y+1, 2);
199 0 0       0 if ($yexp > $xexp) {
200             ### Y bigger ...
201 0         0 $len = $ylen/2;
202 0         0 $exp = $yexp;
203 0         0 $n = $len*$len*2;
204              
205 0         0 $y -= $ylen;
206             ### to: "$x,$y len=$len"
207 0 0 0     0 if ($x == $len-1 && $y == -1) {
208             ### spine ...
209 0         0 return $n;
210             }
211              
212             } else {
213             ### X bigger, fake initial X bit ...
214 0         0 $n = $xlen;
215 0         0 $len = $xlen;
216 0         0 $exp = $xexp;
217 0         0 $n = $len*$len;
218              
219 0 0 0     0 if ($x == $len-1 && $y == $len-1) {
220             ### spine ...
221 0         0 return $n;
222             }
223             }
224              
225 0 0 0     0 if ($x < 0 || $y < 0) {
226 0         0 return undef;
227             }
228              
229             ### $n
230 0         0 my @nbits; # high to low
231              
232              
233 0         0 while ($exp-- >= 0) {
234             ### X at: "$x,$y len=$len"
235             ### assert: $len >= 1
236             ### assert: $x >= 0
237             ### assert: $y >= 0
238             ### assert: $x <= 2*$len
239             ### assert: $y <= 2*$len
240              
241 0 0 0     0 if ($x == $len-1 && $y == $len-1) {
242             ### midpoint X ...
243              
244             # ### nbits HtoL: join('',@nbits)
245             # ### shift off high nbits ...
246             # shift @nbits;
247              
248 0         0 last;
249             }
250 0 0       0 if ($x >= $len) {
251             ### move left, digit 0 ...
252 0         0 $x -= $len;
253 0         0 push @nbits, 0;
254             } else {
255             ### rotate 180, digit 0: "$x,$y"
256 0         0 push @nbits, 1;
257 0         0 $x = $len-2 - $x;
258 0         0 $y = 2*$len-2 - $y;
259 0 0 0     0 if ($x < 0 || $y < 0) {
260             ### outside: "$x,$y"
261 0         0 return undef;
262             }
263             }
264              
265             ### Y at: "$x,$y len=$len"
266             ### assert: $x >= 0
267             ### assert: $y >= 0
268             ### assert: $x <= $len
269             ### assert: $y <= 2*$len
270              
271 0 0 0     0 if ($y == $len-1 && $x == $len/2-1) {
272             ### midpoint Y ...
273              
274 0         0 last;
275             }
276 0 0       0 if ($y >= $len) {
277             ### move down only, digit 1 ...
278 0         0 $y -= $len;
279 0         0 push @nbits, 1;
280             } else {
281             ### rotate 180, digit 0 ...
282 0         0 push @nbits, 0;
283 0         0 $x = $len-2 - $x;
284 0         0 $y = $len-2 - $y;
285 0 0 0     0 if ($x < 0 || $y < 0) {
286             ### outside: "$x,$y"
287 0         0 return undef;
288             }
289             }
290              
291 0         0 $len /= 2;
292             }
293              
294 0 0       0 if ($yexp > $xexp) {
295             } else {
296             ### nbits HtoL: join('',@nbits)
297             ### shift off high nbits ...
298 0         0 shift @nbits;
299             }
300              
301             ### nbits HtoL: join('',@nbits)
302              
303 0 0       0 if ($yexp > $xexp) {
304             } else {
305             }
306              
307              
308 0         0 @nbits = reverse @nbits;
309 0         0 push @nbits, 1;
310              
311             ### nbits HtoL: join('',reverse @nbits)
312 0         0 return $n + digit_join_lowtohigh(\@nbits,2);
313             }
314              
315              
316             # 7
317             # |
318             # 6 | | | | |
319             # 5 *---*---* | *---*---*
320             # 4 | | | | | | |
321             # | | |
322             # 3 *------*------*
323             # | |
324             # 2 | | | | | |
325             # 1 *---*---* *---*---*
326             # 0 | | | |
327             # 0 1 2 3 4 5 6
328             #
329             # not exact
330             sub rect_to_n_range {
331 0     0 1 0 my ($self, $x1,$y1, $x2,$y2) = @_;
332             ### HTree rect_to_n_range(): "$x1,$y1 $x2,$y2"
333              
334 0         0 $x1 = round_nearest($x1);
335 0         0 $x2 = round_nearest($x2);
336 0         0 $y1 = round_nearest($y1);
337 0         0 $y2 = round_nearest($y2);
338              
339 0         0 $x2 = max($x1,$x2);
340 0         0 $y2 = max($y1,$y2);
341 0 0 0     0 if ($x2 < 0 || $y2 < 0) {
342             ### all outside first quadrant ...
343 0         0 return (1,0);
344             }
345              
346 0         0 my ($pow) = round_down_pow(max(2*$x2, $y2) + 1,
347             2);
348 0         0 return (1, $pow*$pow);
349              
350             # my ($xpow) = round_down_pow($x2+1, 2);
351             # my ($ypow) = round_down_pow($y2+1, 2);
352             # if ($xpow > $ypow) {
353             # return (1, 2*$xpow*$xpow);
354             # } else {
355             # return (1, $ypow*$ypow);
356             # }
357             }
358              
359             # 5*2^n, 6*2^n successively
360             # being start of last two rows of each sub-tree
361             # depth=13 160 10100000
362             # depth=14 192 11000000
363             #
364             sub tree_depth_to_n {
365 0     0 1 0 my ($self, $depth) = @_;
366             ### HTree depth_to_n(): $depth
367 0         0 $depth = int($depth);
368 0 0       0 if ($depth < 3) {
369 0 0       0 if ($depth < 0) { return undef; }
  0         0  
370 0         0 return $depth + 1;
371             }
372 0         0 ($depth, my $rem) = _divrem ($depth-3, 2);
373 0         0 return (5+$rem) * 2**$depth;
374             }
375              
376             # spine 2^n
377             #
378             sub tree_depth_to_n_end {
379 0     0 1 0 my ($self, $depth) = @_;
380             ### HTree depth_to_n(): $depth
381 0 0       0 if ($depth < 0) {
382 0         0 return undef;
383             }
384 0         0 return 2**int($depth);
385             }
386              
387             # 0 N=1
388             # |
389             # 1 N=2--
390             # | \
391             # 2 N=3 N=4-------
392             # | \
393             # 3 N=5 N=8 ------------
394             # / \ | \
395             # 4 N=6 N=7 N=9 N=16----------
396             # / \ | \
397             # 5 N=10 N=11 N=17 N=32------
398             # / \ / \ / \ | \
399             # 6 N=12 N=13 N=14 N=15 N=18 N=19 N=33 N=64---
400             # / \ / \ | \
401             # N=20 [4of] N=34 N=35 N=65 N=128
402             # 0 1 = 1
403             # 1 1 = 1
404             # 2 2 = 1 + 1
405             # 3 2 = 1 + 1
406             # 4 4 = 2 + 1 + 1
407             # 5 4 = 2 + 1 + 1
408             # 6 8 = 4 + 2 + 1 + 1
409             # 7 8
410             # 8 16
411             # 9 16
412             #
413             # 1 + sum i=0 to floor(depth/2)-1 of 2^i
414             # = 2^floor(depth/2)
415              
416             sub tree_depth_to_width {
417 0     0 1 0 my ($self, $depth) = @_;
418             ### HTree tree_n_to_subheight(): $depth
419              
420 0 0       0 if ($depth < 0) {
421 0         0 return undef;
422             }
423 0         0 $depth = int($depth/2);
424 0         0 return 2**$depth;
425             }
426              
427              
428             sub tree_n_to_depth {
429 0     0 1 0 my ($self, $n) = @_;
430             ### HTree n_to_depth(): $n
431              
432 0 0       0 if ($n < 1) { return undef; }
  0         0  
433 0         0 $n = int($n);
434 0 0       0 if (is_infinite($n)) { return $n; }
  0         0  
435 0         0 my ($pow,$depth) = round_down_pow($n,2);
436 0 0       0 if ($n -= $pow) {
437 0         0 ($pow, my $exp) = round_down_pow($n,2);
438 0         0 $depth += $exp+1;
439             }
440 0         0 return $depth;
441             }
442              
443              
444             # (n-pow)*2 + pow
445             # = 2*n-2*pow+pow
446             # = 2*n-pow
447             # (n-pow)*2 < pow
448             # 2*n-2*pow < pow
449             # 2*n-pow < 2*pow
450             # 1011 -> 011 -> 110 -> 1110
451             sub tree_n_children {
452 0     0 1 0 my ($self, $n) = @_;
453             ### HTree tree_n_children(): $n
454              
455 0 0       0 if ($n < 1) {
456 0         0 return;
457             }
458              
459 0         0 my ($pow) = round_down_pow($n,2);
460 0 0       0 if ($pow == 1) {
461 0         0 return $n+1;
462             }
463 0 0       0 if ($n == $pow) {
464 0         0 return ($n+1, 2*$n);
465             }
466              
467 0         0 $n *= 2;
468 0         0 $n -= $pow;
469 0 0       0 if ($n < 2*$pow) {
470 0         0 return ($n, $n+1);
471             } else {
472 0         0 return;
473             }
474             }
475              
476             # 1
477             # 2 3
478             # 4 5, 6,7
479             # 101 11x
480             # 8 9, 10,11, 12,13,14,15
481             # 1001 1010,1011 1100,1101
482             #
483             # (n-pow)/2 + pow
484             # = (n-pow+2*pow)/2
485             # = (n+pow)/2
486             #
487             sub tree_n_parent {
488 0     0 1 0 my ($self, $n) = @_;
489             ### HTree tree_n_parent(): $n
490              
491 0 0       0 if ($n < 2) {
492 0         0 return undef;
493             }
494 0         0 my ($pow) = round_down_pow($n,2);
495 0 0       0 if ($n == $pow) {
496 0         0 return $n/2;
497             }
498 0         0 return ($n + $pow - ($n%2)) / 2;
499             }
500              
501             # length of the run of 0s immediately below the high 1-bit
502             # 10001111
503             # ^^^------ 3 0-bits
504             #
505             sub tree_n_to_subheight {
506 0     0 1 0 my ($self, $n) = @_;
507             ### HTree tree_n_to_subheight(): $n
508              
509 0 0       0 if ($n < 2) {
510 0         0 return undef;
511             }
512 0         0 my ($pow,$exp) = round_down_pow($n,2);
513 0 0       0 if ($n == $pow) {
514 0         0 return undef; # spine, infinite
515             }
516 0         0 $n -= $pow;
517 0         0 ($pow, my $exp2) = round_down_pow($n,2);
518 0         0 return $exp - $exp2 - 1;
519             }
520              
521              
522             #------------------------------------------------------------------------------
523             # levels
524              
525             sub level_to_n_range {
526 4     4 1 665 my ($self, $level) = @_;
527 4         13 return (1, 2**(2*$level+1) - 1);
528             }
529             sub n_to_level {
530 0     0 1   my ($self, $n) = @_;
531 0 0         if ($n < 1) { return undef; }
  0            
532 0 0         if (is_infinite($n)) { return $n; }
  0            
533 0           $n = round_nearest($n);
534 0           _divrem_mutate ($n, 2);
535 0           my ($pow, $exp) = round_up_pow ($n+1, 4);
536 0           return $exp;
537             }
538              
539             #------------------------------------------------------------------------------
540             1;
541             __END__