File Coverage

blib/lib/Math/PlanePath/ToothpickUpist.pm
Criterion Covered Total %
statement 90 249 36.1
branch 13 90 14.4
condition 4 37 10.8
subroutine 23 33 69.7
pod 12 12 100.0
total 142 421 33.7


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              
19             #
20             # A151567 four copies of leftist toothpicks
21             # becomes 2*left(n)+2*left(n+1)-4n-1 undoubling diagonals
22             #
23             # A151565 ,1,1,2,2,2,2, 4, 4, 2, 2,4,4,4,4,8,8,2,2,4,4,4,4,8,8,4,4,8,8,8,8,16,
24             # A151566 ,0,1,2,4,6,8,10,14,18,20,22,26,30,34,38,46,54,56,58,62,66,70,74,82,90
25              
26             # A175099,A160018 leftist closed rectangles
27              
28             package Math::PlanePath::ToothpickUpist;
29 1     1   1348 use 5.004;
  1         3  
30 1     1   5 use strict;
  1         1  
  1         25  
31              
32 1     1   6 use vars '$VERSION', '@ISA';
  1         1  
  1         58  
33             $VERSION = 18;
34 1     1   1136 use Math::PlanePath;
  1         6293  
  1         148  
35             @ISA = ('Math::PlanePath');
36              
37              
38             # return $remainder, modify $n
39             # the scalar $_[0] is modified, but if it's a BigInt then a new BigInt is made
40             # and stored there, the bigint value is not changed
41             sub _divrem_mutate {
42 7     7   10 my $d = $_[1];
43 7         55 my $rem;
44 7 50 33     22 if (ref $_[0] && $_[0]->isa('Math::BigInt')) {
45 0         0 ($_[0], $rem) = $_[0]->copy->bdiv($d); # quot,rem in array context
46 0 0 0     0 if (! ref $d || $d < 1_000_000) {
47 0         0 return $rem->numify; # plain remainder if fits
48             }
49             } else {
50 7         9 $rem = $_[0] % $d;
51 7         12 $_[0] = int(($_[0]-$rem)/$d); # exact division stays in UV
52             }
53 7         13 return $rem;
54             }
55              
56              
57             use Math::PlanePath::Base::Generic
58 1         61 '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         72 'round_up_pow',
62             'round_down_pow',
63             'bit_split_lowtohigh',
64 1     1   749 'digit_join_lowtohigh';
  1         1644  
65              
66             # uncomment this to run the ### lines
67             # use Smart::Comments;
68              
69              
70 1     1   5 use constant default_n_start => 0;
  1         2  
  1         49  
71 1     1   5 use constant class_x_negative => 1;
  1         2  
  1         42  
72 1     1   5 use constant class_y_negative => 0;
  1         2  
  1         39  
73 1     1   5 use constant x_negative_at_n => 2;
  1         3  
  1         44  
74 1     1   5 use constant sumxy_minimum => 0; # triangular X>=-Y
  1         2  
  1         38  
75 1     1   5 use constant diffxy_maximum => 0; # triangular X<=Y so X-Y<=0
  1         1  
  1         40  
76 1     1   4 use constant dy_minimum => 0; # across rows dY=0
  1         2  
  1         43  
77 1     1   5 use constant dy_maximum => 1; # then up dY=1 at end
  1         1  
  1         57  
78 1     1   5 use constant tree_num_children_list => (0,1,2);
  1         2  
  1         62  
79 1     1   4 use constant dir_maximum_dxdy => (-1,0); # West
  1         2  
  1         792  
80              
81              
82             #------------------------------------------------------------------------------
83             sub new {
84 4     4 1 645 my $self = shift->SUPER::new(@_);
85 4 50       32 if (! defined $self->{'n_start'}) {
86 4         21 $self->{'n_start'} = $self->default_n_start;
87             }
88 4         11 return $self;
89             }
90              
91             sub n_to_xy {
92 0     0 1 0 my ($self, $n) = @_;
93             ### ToothpickUpist n_to_xy(): $n
94              
95             # written as $n-n_start() rather than "-=" so as to provoke an
96             # uninitialized value warning if $n==undef
97 0         0 $n = $n - $self->{'n_start'}; # N=0 basis
98              
99 0 0       0 if ($n < 0) {
100 0         0 return;
101             }
102 0 0 0     0 if ($n == 0 || is_infinite($n)) {
103 0         0 return ($n,$n);
104             }
105              
106             # this frac behaviour unspecified yet
107             {
108 0         0 my $int = int($n);
  0         0  
109             ### $int
110             ### $n
111 0 0       0 if ($n != $int) {
112 0         0 my $frac = $n - $int; # inherit possible BigFloat
113 0         0 $int += $self->{'n_start'};
114 0         0 my ($x1,$y1) = $self->n_to_xy($int);
115 0         0 my ($x2,$y2) = $self->n_to_xy($int+1);
116 0         0 my $dx = $x2-$x1;
117 0         0 my $dy = $y2-$y1;
118 0         0 return ($frac*$dx + $x1, $frac*$dy + $y1);
119             }
120 0         0 $n = $int; # BigFloat int() gives BigInt, use that
121             }
122             ### $n
123              
124 0         0 my ($depthbits, $lowbit, $ndepth) = _n0_to_depthbits($n);
125             ### $depthbits
126             ### $ndepth
127             ### n remainder: $n-$ndepth
128              
129 0         0 my @nbits = bit_split_lowtohigh($n-$ndepth); # offset into row
130              
131             ### @nbits
132             ### $lowbit
133              
134             # Where there's a 0-bit in the depth remains a 0-bit.
135             # Where there's a 1-bit in the depth takes a bit from Noffset.
136             # Small Noffset has less bits than the depth 1s, hence "|| 0".
137             #
138 0 0 0     0 my @xbits = map {$_ && (shift @nbits || 0)} @$depthbits;
  0         0  
139             ### @xbits
140              
141 0         0 my $zero = $n * 0;
142 0         0 my $x = digit_join_lowtohigh (\@xbits, 2, $zero);
143 0         0 my $y = digit_join_lowtohigh ($depthbits, 2, $zero);
144              
145             ### Y without lowbit: $y
146              
147 0         0 return (2*$x-$y, # triangular style
148             $y + $lowbit);
149             }
150              
151             sub xy_to_n {
152 0     0 1 0 my ($self, $x, $y) = @_;
153             ### ToothpickUpist xy_to_n(): "$x, $y"
154              
155 0         0 $y = round_nearest ($y);
156 0         0 $x = round_nearest($x);
157              
158             # odd points X!=Ymod2 are the second copy of the triangle, go to Y-1 for them
159 0         0 $x += $y;
160 0         0 my $lowbit = _divrem_mutate ($x, 2);
161 0         0 $y -= $lowbit;
162             ### odd adjusted xy: "$x,$y"
163              
164 0         0 return _right_xy_to_n ($self, $x,$y, $lowbit);
165             }
166              
167             # with X,Y in the align="right" style,
168             #
169             # |
170             sub _right_xy_to_n {
171 7     7   9 my ($self, $x, $y, $lowbit) = @_;
172             ### _right_xy_to_n(): "x=$x y=$y lowbit=$lowbit"
173              
174 7 50 33     50 unless ($x >= 0 && $x <= $y && $y >= 0) {
      33        
175             ### outside horizontal row range ...
176 0         0 return undef;
177             }
178 7 50       16 if (is_infinite($y)) {
179 0         0 return $y;
180             }
181              
182 7         40 my $zero = ($y * 0);
183 7         9 my $n = $zero; # inherit bignum 0
184 7         11 my $npower = $zero+2; # inherit bignum 2
185              
186 7         15 my @xbits = bit_split_lowtohigh($x);
187 7         40 my @depthbits = bit_split_lowtohigh($y);
188              
189 7         73 my @nbits; # N offset into row
190 7         14 foreach my $i (0 .. $#depthbits) { # x,y bits low to high
191 28 100       46 if ($depthbits[$i]) {
192 7         7 $n = 2*$n + $npower;
193 7   50     32 push @nbits, $xbits[$i] || 0; # low to high
194             } else {
195 21 50       42 if ($xbits[$i]) {
196 0         0 return undef;
197             }
198             }
199 28         39 $npower *= 3;
200             }
201              
202 7 50       15 if ($lowbit) {
203 0         0 push @nbits, 1;
204             }
205              
206             ### n at left end of y row: $n
207             ### n offset for x: @nbits
208             ### total: $n + digit_join_lowtohigh(\@nbits,2,$zero) + $self->{'n_start'}
209              
210 7         19 return $n + digit_join_lowtohigh(\@nbits,2,$zero) + $self->{'n_start'};
211             }
212              
213             # not exact
214             sub rect_to_n_range {
215 0     0 1 0 my ($self, $x1,$y1, $x2,$y2) = @_;
216             ### ToothpickUpist rect_to_n_range(): "$x1,$y1, $x2,$y2"
217              
218 0         0 $y1 = round_nearest ($y1);
219 0         0 $y2 = round_nearest ($y2);
220 0 0       0 if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1) }
  0         0  
221              
222 0         0 $x1 = round_nearest ($x1);
223 0         0 $x2 = round_nearest ($x2);
224 0 0       0 if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1) }
  0         0  
225              
226 0 0       0 if ($y2 < 0) {
227             ### all negative ...
228 0         0 return (1, 0);
229             }
230 0         0 $y1 -= 1;
231 0 0       0 if ($y1 < 0) {
232 0         0 $y1 = 0;
233             }
234              
235             ### range using: "y1=$y1 y2=$y2"
236              
237 0         0 return (_right_xy_to_n($self, 0,$y1, 0),
238             _right_xy_to_n($self, $y2,$y2, 1));
239             }
240              
241              
242             #------------------------------------------------------------------------------
243 1     1   6 use constant tree_num_roots => 1;
  1         1  
  1         1308  
244              
245             sub tree_n_num_children {
246 0     0 1 0 my ($self, $n) = @_;
247              
248 0         0 $n = $n - $self->{'n_start'}; # N=0 basis
249 0 0 0     0 if (is_infinite($n) || $n < 0) {
250 0         0 return undef;
251             }
252              
253 0         0 my ($depthbits, $lowbit, $ndepth) = _n0_to_depthbits($n);
254 0 0       0 if (! $lowbit) {
255 0         0 return 1;
256             }
257 0 0       0 unless (shift @$depthbits) { # low bit above $lowbit doubling
258             # Depth even (or zero), two children under every point.
259 0         0 return 2;
260             }
261              
262             # Depth odd, single child under some or all points.
263             # When depth==1mod4 it's all points, when depth has more than one
264             # trailing 1-bit then it's only some points.
265             #
266 0         0 $n -= $ndepth; # Noffset into row
267 0         0 my $repbit = _divrem_mutate($n,2);
268 0         0 while (shift @$depthbits) { # low to high
269 0 0       0 if (_divrem_mutate($n,2) != $repbit) {
270 0         0 return 0;
271             }
272             }
273 0         0 return 1;
274             }
275              
276             sub tree_n_children {
277 0     0 1 0 my ($self, $n) = @_;
278             ### tree_n_children(): $n
279              
280 0         0 $n = $n - $self->{'n_start'}; # N=0 basis
281 0 0 0     0 if (is_infinite($n) || $n < 0) {
282 0         0 return;
283             }
284              
285 0         0 my ($depthbits, $lowbit, $ndepth, $nwidth) = _n0_to_depthbits($n);
286 0 0       0 if (! $lowbit) {
287             ### doubled to children at nwidth below ...
288 0         0 return ($n + $nwidth);
289             }
290              
291 0         0 $n -= $ndepth; # Noffset into row
292              
293 0 0       0 if (shift @$depthbits) {
294             # Depth odd, single child under some or all points.
295             # When depth==1mod4 it's all points, when depth has more than one
296             # trailing 1-bit then it's only some points.
297 0         0 while (shift @$depthbits) { # depth==3mod4 or more low 1s
298 0         0 my $repbit = _divrem_mutate($n,2);
299 0 0       0 if (($n % 2) != $repbit) {
300 0         0 return;
301             }
302             }
303 0         0 return $n + $ndepth+$nwidth + $self->{'n_start'};
304              
305             } else {
306             # Depth even (or zero), two children under every point.
307 0         0 $n = 2*$n + $ndepth+$nwidth + $self->{'n_start'};
308 0         0 return ($n,$n+1);
309             }
310             }
311              
312             sub tree_n_parent {
313 0     0 1 0 my ($self, $n) = @_;
314              
315 0 0       0 my ($x,$y) = $self->n_to_xy($n)
316             or return undef;
317              
318 0 0       0 if (($x%2) != ($y%2)) {
319             ### odd, directly down ...
320 0         0 return $self->xy_to_n($x,$y-1);
321             }
322              
323             ### even, to one side or the other ...
324 0         0 my $n_parent = $self->xy_to_n($x-1, $y);
325 0 0       0 if (defined $n_parent) {
326 0         0 return $n_parent;
327             }
328 0         0 return $self->xy_to_n($x+1,$y);
329             }
330              
331             sub tree_n_to_depth {
332 1     1 1 41 my ($self, $n) = @_;
333             ### ToothpickUpist n_to_depth(): $n
334 1         2 $n = $n - $self->{'n_start'};
335 1 50       4 unless ($n >= 0) {
336 0         0 return undef; # negatives, -infinity, NaN
337             }
338 1 50       3 if (is_infinite($n)) {
339 1         8 return $n; # +infinity
340             }
341 0         0 my ($depthbits, $lowbit) = _n0_to_depthbits($n);
342 0         0 unshift @$depthbits, $lowbit;
343 0         0 return digit_join_lowtohigh ($depthbits, 2, $n*0);
344             }
345             sub tree_depth_to_n {
346 8     8 1 59 my ($self, $depth) = @_;
347             ### tree_depth_to_n(): $depth
348 8 50       16 if ($depth >= 0) {
349             # $depth==+infinity becomes nan from divrem, prefer to return N=+infinity
350             # for +inf depth
351 8 100       21 if (is_infinite($depth)) {
352 1         8 return $depth;
353             }
354 7         45 my $lowbit = _divrem_mutate($depth,2);
355 7         15 return _right_xy_to_n($self,0,$depth, $lowbit);
356             } else {
357 0         0 return undef;
358             }
359             }
360              
361             sub tree_n_to_subheight {
362 0     0 1 0 my ($self, $n) = @_;
363             ### ToothpickUpist tree_n_to_subheight(): $n
364              
365 0         0 $n = $n - $self->{'n_start'};
366 0 0 0     0 if (is_infinite($n) || $n < 0) {
367 0         0 return undef;
368             }
369 0         0 my ($depthbits, $lowbit, $ndepth) = _n0_to_depthbits($n);
370 0         0 $n -= $ndepth; # remaining offset into row
371 0         0 my @nbits = bit_split_lowtohigh($n);
372              
373             ### $lowbit
374             ### $depthbits
375              
376 0   0     0 my $target = $nbits[0] || 0;
377 0         0 foreach my $i (0 .. $#$depthbits) {
378 0 0       0 unless ($depthbits->[$i] ^= 1) { # flip 0<->1, at original==1 take nbit
379 0 0 0     0 if ((shift @nbits || 0) != $target) {
380 0         0 unshift @$depthbits, 1-$lowbit;
381 0         0 $#$depthbits = $i;
382             ### $depthbits
383 0         0 return digit_join_lowtohigh($depthbits, 2, $n*0);
384             }
385             }
386             }
387 0         0 return undef; # first or last of row, infinite
388             }
389              
390             sub _EXPERIMENTAL__tree_n_to_leafdist {
391 0     0   0 my ($self, $n) = @_;
392             ### _EXPERIMENTAL__tree_n_to_leafdist(): $n
393              
394 0         0 $n = $n - $self->{'n_start'}; # N=0 basis
395 0 0 0     0 if (is_infinite($n) || $n < 0) {
396 0         0 return undef;
397             }
398              
399             # depth bits leafdist
400             # 0 0,0 7
401             # 1 0,1 6
402             # 2 1,0 5
403             # 3 1,1 4
404             # 4 1,0,0 3
405             # 5 1,0,1 2
406             # 6 1,1,0 1 or 9
407             # 7 1,1,1 0 or 8
408             # ignore $lowbit until last, bits above same as SierpinskiTriangle
409             #
410 0         0 my ($depthbits, $lowbit, $ndepth) = _n0_to_depthbits($n);
411 0         0 $lowbit = 1-$lowbit;
412              
413 0   0     0 my $ret = 6 - 2*((shift @$depthbits)||0);
414 0 0       0 if (shift @$depthbits) { $ret -= 4; }
  0         0  
415             ### $ret
416 0 0       0 if ($ret) {
417 0         0 return $ret + $lowbit;
418             }
419              
420 0         0 $n -= $ndepth;
421             ### Noffset into row: $n
422              
423             # Low bits of Nrem unchanging while trailing 1-bits in @depthbits,
424             # to distinguish between leaf or non-leaf. Same as tree_n_children().
425             #
426 0         0 my $repbit = _divrem_mutate($n,2); # low bit of $n
427             ### $repbit
428 0         0 do {
429             ### next bit: $n%2
430 0 0       0 if (_divrem_mutate($n,2) != $repbit) { # bits of $n offset low to high
431 0         0 return $lowbit; # is a leaf
432             }
433             } while (shift @$depthbits);
434 0         0 return 8+$lowbit; # is a non-leaf
435             }
436              
437             # Ndepth = 2 * ( 3^a first N at this depth
438             # + 2 * 3^b
439             # + 2^2 * 3^c
440             # + 2^3 * 3^d
441             # + ... )
442              
443             sub _n0_to_depthbits {
444 0     0   0 my ($n) = @_;
445             ### _n0_to_depthbits(): $n
446              
447 0 0       0 if ($n == 0) {
448 0         0 return ([], 0, 0, 1);
449             }
450              
451 0         0 my ($nwidth, $bitpos) = round_down_pow ($n/2, 3);
452             ### nwidth power-of-3: $nwidth
453             ### $bitpos
454              
455 0         0 $nwidth *= 2; # two of each row
456              
457 0         0 my @depthbits;
458 0         0 my $ndepth = 0;
459 0         0 for (;;) {
460             ### at: "n=$n nwidth=$nwidth bitpos=$bitpos depthbits=".join(',',map{$_||0}@depthbits)
461              
462 0 0       0 if ($n >= $ndepth + $nwidth) {
463 0         0 $depthbits[$bitpos] = 1;
464 0         0 $ndepth += $nwidth;
465 0         0 $nwidth *= 2;
466             } else {
467 0         0 $depthbits[$bitpos] = 0;
468             }
469 0 0       0 last unless --$bitpos >= 0;
470 0         0 $nwidth /= 3;
471             }
472              
473             # Nwidth = 2**count1bits(depth)
474             ### assert: $nwidth == 2*(1 << scalar(grep{$_}@depthbits))
475              
476             # first or second of the two of each row
477 0         0 $nwidth /= 2;
478 0 0       0 my $lowbit = ($n >= $ndepth + $nwidth ? 1 : 0);
479 0 0       0 if ($lowbit) {
480 0         0 $ndepth += $nwidth;
481             }
482             ### final depthbits: join(',',@depthbits)
483              
484 0         0 return (\@depthbits, $lowbit, $ndepth, $nwidth);
485             }
486              
487             #------------------------------------------------------------------------------
488             # levels
489              
490             sub level_to_n_range {
491 11     11 1 617 my ($self, $level) = @_;
492 11         30 return (0, 2* 3**$level - 1);
493             }
494             sub n_to_level {
495 0     0 1   my ($self, $n) = @_;
496 0 0         if ($n < 0) { return undef; }
  0            
497 0 0         if (is_infinite($n)) { return $n; }
  0            
498 0           $n = round_nearest($n);
499 0           _divrem_mutate ($n, 2);
500 0           my ($pow, $exp) = round_up_pow ($n+1, 3);
501 0           return $exp;
502             }
503              
504             #------------------------------------------------------------------------------
505             1;
506             __END__