File Coverage

blib/lib/Math/PlanePath/PythagoreanTree.pm
Criterion Covered Total %
statement 239 391 61.1
branch 93 190 48.9
condition 49 82 59.7
subroutine 33 61 54.1
pod 25 25 100.0
total 439 749 58.6


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             # math-image --path=PythagoreanTree --all --scale=3
20              
21             # http://sunilchebolu.wordpress.com/pythagorean-triples-and-the-integer-points-on-a-hyperboloid/
22              
23             # http://www.math.uconn.edu/~kconrad/blurbs/ugradnumthy/pythagtriple.pdf
24             #
25             # http://www.math.ou.edu/~dmccullough/teaching/pythagoras1.pdf
26             # http://www.math.ou.edu/~dmccullough/teaching/pythagoras2.pdf
27             #
28             # http://www.microscitech.com/pythag_eigenvectors_invariants.pdf
29             #
30              
31              
32             package Math::PlanePath::PythagoreanTree;
33 2     2   3552 use 5.004;
  2         8  
34 2     2   12 use strict;
  2         4  
  2         45  
35 2     2   11 use Carp 'croak';
  2         4  
  2         95  
36              
37 2     2   11 use vars '$VERSION', '@ISA';
  2         3  
  2         168  
38             $VERSION = 127;
39 2     2   1370 use Math::PlanePath;
  2         6  
  2         134  
40             *_divrem = \&Math::PlanePath::_divrem;
41             *_sqrtint = \&Math::PlanePath::_sqrtint;
42             @ISA = ('Math::PlanePath');
43              
44             #use List::Util 'min','max';
45             *min = \&Math::PlanePath::_min;
46             *max = \&Math::PlanePath::_max;
47              
48             use Math::PlanePath::Base::Generic
49 2         131 'is_infinite',
50 2     2   12 'round_nearest';
  2         4  
51             use Math::PlanePath::Base::Digits
52 2         98 'round_down_pow',
53             'digit_split_lowtohigh',
54 2     2   637 'digit_join_lowtohigh';
  2         5  
55 2     2   1099 use Math::PlanePath::GrayCode;
  2         5  
  2         64  
56              
57             # uncomment this to run the ### lines
58             # use Smart::Comments;
59              
60 2     2   12 use constant class_x_negative => 0;
  2         4  
  2         95  
61 2     2   11 use constant class_y_negative => 0;
  2         4  
  2         77  
62 2     2   11 use constant tree_num_children_list => (3); # complete ternary tree
  2         4  
  2         73  
63 2     2   10 use constant tree_n_to_subheight => undef; # complete tree, all infinity
  2         4  
  2         150  
64              
65 2         817 use constant parameter_info_array =>
66             [ { name => 'tree_type',
67             share_key => 'tree_type_uadfb',
68             display => 'Tree Type',
69             type => 'enum',
70             default => 'UAD',
71             choices => ['UAD','UArD','FB','UMT'],
72             },
73             { name => 'coordinates',
74             share_key => 'coordinates_abcpqsm',
75             display => 'Coordinates',
76             type => 'enum',
77             default => 'AB',
78             choices => ['AB','AC','BC','PQ', 'SM','SC','MC',
79             # 'BA'
80             # 'UV', # q from x=y diagonal down at 45-deg
81             # 'RS','ST', # experimental
82             ],
83             },
84             { name => 'digit_order',
85             display => 'Digit Order',
86             type => 'enum',
87             default => 'HtoL',
88             choices => ['HtoL','LtoH'],
89             },
90 2     2   11 ];
  2         4  
91              
92             {
93             my %UAD_coordinates_always_right = (PQ => 1,
94             AB => 1,
95             AC => 1);
96             sub turn_any_left {
97 0     0 1 0 my ($self) = @_;
98             return ! ($self->{'tree_type'} eq 'UAD'
99 0   0     0 && $UAD_coordinates_always_right{$self->{'coordinates'}});
100             }
101             }
102             {
103             my %UAD_coordinates_always_left = (BC => 1);
104             sub turn_any_right {
105 0     0 1 0 my ($self) = @_;
106             return ! ($self->{'tree_type'} eq 'UAD'
107 0   0     0 && $UAD_coordinates_always_left{$self->{'coordinates'}});
108             }
109             }
110             {
111             my %UMT_coordinates_any_straight = (BC => 1, # UMT at N=5
112             PQ => 1); # UMT at N=5
113             sub turn_any_straight {
114 0     0 1 0 my ($self) = @_;
115             return ($self->{'tree_type'} eq 'UMT'
116 0   0     0 && $UMT_coordinates_any_straight{$self->{'coordinates'}});
117             }
118             }
119              
120              
121             #------------------------------------------------------------------------------
122             {
123             my %coordinate_minimum = (A => 3,
124             B => 4,
125             C => 5,
126             P => 2,
127             Q => 1,
128             S => 3,
129             M => 4,
130             );
131             sub x_minimum {
132 0     0 1 0 my ($self) = @_;
133 0         0 return $coordinate_minimum{substr($self->{'coordinates'},0,1)};
134             }
135             sub y_minimum {
136 0     0 1 0 my ($self) = @_;
137 0         0 return $coordinate_minimum{substr($self->{'coordinates'},1)};
138             }
139             }
140             {
141             my %diffxy_minimum = (PQ => 1, # octant X>=Y+1 so X-Y>=1
142             );
143             sub diffxy_minimum {
144 0     0 1 0 my ($self) = @_;
145 0         0 return $diffxy_minimum{$self->{'coordinates'}};
146             }
147             }
148             {
149             my %diffxy_maximum = (AC => -2, # C>=A+2 so X-Y<=-2
150             BC => -1, # C>=B+1 so X-Y<=-1
151             SM => -1, # S
152             SC => -2, # S
153             MC => -1, # M
154             );
155             sub diffxy_maximum {
156 0     0 1 0 my ($self) = @_;
157 0         0 return $diffxy_maximum{$self->{'coordinates'}};
158             }
159             }
160             {
161             my %absdiffxy_minimum = (PQ => 1,
162             AB => 1, # X=Y never occurs
163             BA => 1, # X=Y never occurs
164             AC => 2, # C>=A+2 so abs(X-Y)>=2
165             BC => 1,
166             SM => 1, # X=Y never occurs
167             SC => 2, # X<=Y-2
168             MC => 1, # X=Y never occurs
169             );
170             sub absdiffxy_minimum {
171 0     0 1 0 my ($self) = @_;
172 0         0 return $absdiffxy_minimum{$self->{'coordinates'}};
173             }
174             }
175 2     2   20 use constant gcdxy_maximum => 1; # no common factor
  2         5  
  2         4038  
176              
177             {
178             my %absdx_minimum = ('AB,UAD' => 2,
179             'AB,FB' => 2,
180             'AB,UMT' => 2,
181              
182             'AC,UAD' => 2,
183             'AC,FB' => 2,
184             'AC,UMT' => 2,
185              
186             'BC,UAD' => 4, # at N=37
187             'BC,FB' => 4, # at N=2 X=12,Y=13
188             'BC,UMT' => 4, # at N=2 X=12,Y=13
189              
190             'PQ,UAD' => 0,
191             'PQ,FB' => 0,
192             'PQ,UMT' => 0,
193              
194             'SM,UAD' => 1,
195             'SM,FB' => 1,
196             'SM,UMT' => 2,
197              
198             'SC,UAD' => 1,
199             'SC,FB' => 1,
200             'SC,UMT' => 1,
201              
202             'MC,UAD' => 3,
203             'MC,FB' => 3,
204             'MC,UMT' => 1,
205             );
206             sub absdx_minimum {
207 0     0 1 0 my ($self) = @_;
208 0   0     0 return $absdx_minimum{"$self->{'coordinates'},$self->{'tree_type'}"} || 0;
209             }
210             }
211             {
212             my %absdy_minimum = ('AB,UAD' => 4,
213             'AB,FB' => 4,
214             'AB,UMT' => 4,
215              
216             'AC,UAD' => 4,
217             'AC,FB' => 4,
218             'BC,UAD' => 4,
219             'BC,FB' => 4,
220             'PQ,UAD' => 0,
221             'PQ,FB' => 1,
222              
223             'SM,UAD' => 3,
224             'SM,FB' => 3,
225             'SM,UMT' => 1,
226              
227             'SC,UAD' => 4,
228             'SC,FB' => 4,
229             'MC,UAD' => 4,
230             'MC,FB' => 4,
231             );
232             sub absdy_minimum {
233 0     0 1 0 my ($self) = @_;
234 0   0     0 return $absdy_minimum{"$self->{'coordinates'},$self->{'tree_type'}"} || 0;
235             }
236             }
237              
238             {
239             my %dir_minimum_dxdy = (# AB apparent minimum dX=16,dY=8
240             'AB,UAD' => [16,8],
241             'AC,UAD' => [1,1], # it seems
242             # 'BC,UAD' => [1,0], # infimum
243             # 'SM,UAD' => [1,0], # infimum
244             # 'SC,UAD' => [1,0], # N=255 dX=7,dY=0
245             # 'MC,UAD' => [1,0], # infimum
246              
247             # 'SM,FB' => [1,0], # infimum
248             # 'SC,FB' => [1,0], # infimum
249             # 'SM,FB' => [1,0], # infimum
250              
251             'AB,UMT' => [6,12], # it seems
252              
253             # N=ternary 1111111122 dx=118,dy=40
254             # in general dx=3*4k-2 dy=4k
255             'AC,UMT' => [3,1], # infimum
256             #
257             # 'BC,UMT' => [1,0], # N=31 dX=72,dY=0
258             'PQ,UMT' => [1,1], # N=1
259             'SM,UMT' => [1,0], # infiumum dX=big,dY=3
260             'SC,UMT' => [3,1], # like AC
261             # 'MC,UMT' => [1,0], # at N=31
262             );
263             sub dir_minimum_dxdy {
264 0     0 1 0 my ($self) = @_;
265 0 0       0 return @{$dir_minimum_dxdy{"$self->{'coordinates'},$self->{'tree_type'}"}
  0         0  
266             || [1,0] };
267             }
268             }
269             {
270             # AB apparent maximum dX=-6,dY=-12 at N=3
271             # AC apparent maximum dX=-6,dY=-12 at N=3 same
272             # PQ apparent maximum dX=-1,dY=-1
273             my %dir_maximum_dxdy = ('AB,UAD' => [-6,-12],
274             'AC,UAD' => [-6,-12],
275             # 'BC,UAD' => [0,0],
276             'PQ,UAD' => [-1,-1],
277             # 'SM,UAD' => [0,0], # supremum
278             # 'SC,UAD' => [0,0], # supremum
279             # 'MC,UAD' => [0,0], # supremum
280              
281             # 'AB,FB' => [0,0],
282             # 'AC,FB' => [0,0],
283             'BC,FB' => [1,-1],
284             # 'PQ,FB' => [0,0],
285             # 'SM,FB' => [0,0], # supremum
286             # 'SC,FB' => [0,0], # supremum
287             # 'MC,FB' => [0,0], # supremum
288              
289             # N=ternary 1111111122 dx=118,dy=-40
290             # in general dx=3*4k-2 dy=-4k
291             'AB,UMT' => [3,-1], # supremum
292             #
293             'AC,UMT' => [-10,-20], # at N=9 apparent maximum
294             # 'BC,UMT' => [0,0], # apparent approach
295             'PQ,UMT' => [1,-1], # N=2
296             # 'SM,UMT' => [0,0], # supremum dX=big,dY=-1
297             'SC,UMT' => [-3,-5], # apparent approach
298             # 'MC,UMT' => [0,0], # supremum dX=big,dY=-small
299             );
300             sub dir_maximum_dxdy {
301 0     0 1 0 my ($self) = @_;
302 0 0       0 return @{$dir_maximum_dxdy{"$self->{'coordinates'},$self->{'tree_type'}"}
  0         0  
303             || [0,0]};
304             }
305             }
306              
307             #------------------------------------------------------------------------------
308              
309             sub _noop {
310 1340     1340   3320 return @_;
311             }
312             my %xy_to_pq = (AB => \&_ab_to_pq,
313             AC => \&_ac_to_pq,
314             BC => \&_bc_to_pqa, # ignoring extra $a return
315             PQ => \&_noop,
316             SM => \&_sm_to_pq,
317             SC => \&_sc_to_pq,
318             MC => \&_mc_to_pq,
319             UV => \&_uv_to_pq,
320             RS => \&_rs_to_pq,
321             ST => \&_st_to_pq,
322             );
323             my %pq_to_xy = (AB => \&_pq_to_ab,
324             AC => \&_pq_to_ac,
325             BC => \&_pq_to_bc,
326             PQ => \&_noop,
327             SM => \&_pq_to_sm,
328             SC => \&_pq_to_sc,
329             MC => \&_pq_to_mc,
330             UV => \&_pq_to_uv,
331             RS => \&_pq_to_rs,
332             ST => \&_pq_to_st,
333             );
334              
335             my %tree_types = (UAD => 1,
336             UArD => 1,
337             FB => 1,
338             UMT => 1);
339             my %digit_orders = (HtoL => 1,
340             LtoH => 1);
341             sub new {
342 21     21 1 2939 my $self = shift->SUPER::new (@_);
343             {
344 21   50     103 my $digit_order = ($self->{'digit_order'} ||= 'HtoL');
345 21 50       58 $digit_orders{$digit_order}
346             || croak "Unrecognised digit_order option: ",$digit_order;
347             }
348             {
349 21   100     43 my $tree_type = ($self->{'tree_type'} ||= 'UAD');
  21         61  
350 21 50       51 $tree_types{$tree_type}
351             || croak "Unrecognised tree_type option: ",$tree_type;
352             }
353             {
354 21   100     31 my $coordinates = ($self->{'coordinates'} ||= 'AB');
  21         38  
  21         60  
355 21   33     65 $self->{'xy_to_pq'} = $xy_to_pq{$coordinates}
356             || croak "Unrecognised coordinates option: ",$coordinates;
357 21         46 $self->{'pq_to_xy'} = $pq_to_xy{$coordinates};
358             }
359 21         60 return $self;
360             }
361              
362             sub n_to_xy {
363 56     56 1 5483 my ($self, $n) = @_;
364             ### PythagoreanTree n_to_xy(): $n
365              
366 56 50       140 if ($n < 1) { return; }
  0         0  
367 56 50       141 if (is_infinite($n)) { return ($n,$n); }
  0         0  
368              
369             {
370 56         99 my $int = int($n);
  56         82  
371 56 50       115 if ($n != $int) {
372 0         0 my $frac = $n - $int; # inherit possible BigFloat/BigRat
373 0         0 my ($x1,$y1) = $self->n_to_xy($int);
374 0         0 my ($x2,$y2) = $self->n_to_xy($int+1);
375 0         0 my $dx = $x2-$x1;
376 0         0 my $dy = $y2-$y1;
377 0         0 return ($frac*$dx + $x1, $frac*$dy + $y1);
378             }
379             }
380              
381 56         116 return &{$self->{'pq_to_xy'}}(_n_to_pq($self,$n));
  56         120  
382             }
383              
384             # maybe similar n_to_rsquared() as C^2=(P^2+Q^2)^2
385             sub n_to_radius {
386 0     0 1 0 my ($self, $n) = @_;
387              
388 0 0 0     0 if (($self->{'coordinates'} eq 'AB'
      0        
389             || $self->{'coordinates'} eq 'BA'
390             || $self->{'coordinates'} eq 'SM')
391             && $n == int($n)) {
392 0 0       0 if ($n < 1) { return undef; }
  0         0  
393 0 0       0 if (is_infinite($n)) { return $n; }
  0         0  
394 0         0 my ($p,$q) = _n_to_pq($self,$n);
395 0         0 return $p*$p + $q*$q; # C=P^2+Q^2
396             }
397              
398 0         0 return $self->SUPER::n_to_radius($n);
399             }
400              
401             sub _n_to_pq {
402 56     56   96 my ($self, $n) = @_;
403              
404 56         109 my $ndigits = _n_to_digits_lowtohigh($n);
405             ### $ndigits
406              
407 56 50       143 if ($self->{'tree_type'} eq 'UArD') {
408 0         0 Math::PlanePath::GrayCode::_digits_to_gray_reflected($ndigits,3);
409             ### gray: $ndigits
410             }
411 56 50       107 if ($self->{'digit_order'} eq 'HtoL') {
412 56         115 @$ndigits = reverse @$ndigits;
413             ### reverse: $ndigits
414             }
415              
416 56         87 my $zero = $n * 0;
417              
418 56         81 my $p = 2 + $zero;
419 56         83 my $q = 1 + $zero;
420              
421 56 100       115 if ($self->{'tree_type'} eq 'FB') {
    50          
422             ### FB ...
423              
424 26         82 foreach my $digit (@$ndigits) { # high to low, possibly $digit=undef
425             ### $p
426             ### $q
427             ### $digit
428              
429 42 100       73 if ($digit) {
430 28 100       65 if ($digit == 1) {
431 14         20 $q = $p-$q; # (2p, p-q) M2
432 14         27 $p *= 2;
433             } else {
434             # ($p,$q) = (2*$p, $p+$q);
435 14         46 $q += $p; # (p+q, 2q) M3
436 14         28 $p *= 2;
437             }
438             } else { # $digit == 0
439             # ($p,$q) = ($p+$q, 2*$q);
440 14         23 $p += $q; # (p+q, 2q) M1
441 14         23 $q *= 2;
442             }
443             }
444             } elsif ($self->{'tree_type'} eq 'UMT') {
445             ### UMT ...
446              
447 0         0 foreach my $digit (@$ndigits) { # high to low, possibly $digit=undef
448             ### $p
449             ### $q
450             ### $digit
451              
452 0 0       0 if ($digit) {
453 0 0       0 if ($digit == 1) {
454 0         0 $q = $p-$q; # (2p, p-q) M2
455 0         0 $p *= 2;
456             } else { # $digit == 2
457 0         0 $p += 3*$q; # T
458 0         0 $q *= 2;
459             }
460             } else { # $digit == 0
461             # ($p,$q) = ($p+$q, 2*$q);
462 0         0 ($p,$q) = (2*$p-$q, $p); # "U" = (2p-q, p)
463             }
464             }
465             } else {
466             ### UAD or UArD ...
467             ### assert: $self->{'tree_type'} eq 'UAD' || $self->{'tree_type'} eq 'UArD'
468              
469             # # Could optimize high zeros as repeated U
470             # # high zeros as repeated U: $depth-scalar(@$ndigits)
471             # # U^0 = p, q
472             # # U^1 = 2p-q, p eg. P=2,Q=1 is 2*2-1,2 = 3,2
473             # # U^2 = 3p-2q, 2p-q eg. P=2,Q=1 is 3*2-2*1,2*2-1 = 4,3
474             # # U^3 = 4p-3q, 3p-2q
475             # # U^k = (k+1)p-kq, kp-(k-1)q for k>=2
476             # # = p + k*(p-q), k*(p-q)+q
477             # # and with initial p=2,q=1
478             # # U^k = 2+k, 1+k
479             # #
480             # $q = $depth - $#ndigits + $zero; # count high zeros + 1
481             # $p = $q + 1 + $zero;
482              
483 30         61 foreach my $digit (@$ndigits) { # high to low, possibly $digit=undef
484             ### $p
485             ### $q
486             ### $digit
487              
488 52 100       93 if ($digit) {
489 34 100       55 if ($digit == 1) {
490 18         68 ($p,$q) = (2*$p+$q, $p); # "A" = (2p+q, p)
491             } else {
492 16         30 $p += 2*$q; # "D" = (p+2q, q)
493             }
494             } else { # $digit==0
495 18         38 ($p,$q) = (2*$p-$q, $p); # "U" = (2p-q, p)
496             }
497             }
498              
499             }
500              
501             ### final pq: "$p, $q"
502              
503 56         124 return ($p, $q);
504             }
505              
506             # _n_to_digits_lowtohigh() returns an arrayref $ndigits which is a list of
507             # ternary digits 0,1,2 from low to high which are the position of $n within
508             # its row of the tree.
509             # The length of the array is the depth.
510             #
511             # depth N N%3 2*N-1 (N-2)/3*2+1
512             # 0 1 1 1 1/3
513             # 1 2 2 3 1
514             # 2 5 2 9 3
515             # 3 14 2 27 9
516             # 4 41 2 81 27 28 + (28/2-1) = 41
517             #
518             # (N-2)/3*2+1 rounded down to pow=3^k gives depth=k+1 and base=pow+(pow+1)/2
519             # is the start of the row base=1,2,5,14,41 etc.
520             #
521             # An easier calculation is 2*N-1 rounded down to pow=3^d gives depth=d and
522             # base=2*pow-1, but 2*N-1 and 2*pow-1 might overflow an integer. Though
523             # just yet round_down_pow() goes into floats and so doesn't preserve 64-bit
524             # integer. So the technique here helps 53-bit float integers, but not right
525             # up to 64-bits.
526             #
527             sub _n_to_digits_lowtohigh {
528 76     76   1586 my ($n) = @_;
529             ### _n_to_digits_lowtohigh(): $n
530              
531 76         130 my @ndigits;
532 76 100       162 if ($n >= 2) {
533 69         192 my ($pow) = _divrem($n-2, 3);
534 69         206 ($pow, my $depth) = round_down_pow (2*$pow+1, 3);
535             ### $depth
536             ### base: $pow + ($pow+1)/2
537             ### offset: $n - $pow - ($pow+1)/2
538 69         203 @ndigits = digit_split_lowtohigh ($n - $pow - ($pow+1)/2, 3);
539 69         223 push @ndigits, (0) x ($depth - $#ndigits); # pad to $depth with 0s
540             }
541             ### @ndigits
542 76         164 return \@ndigits;
543              
544              
545             # {
546             # my ($pow, $depth) = round_down_pow (2*$n-1, 3);
547             #
548             # ### h: 2*$n-1
549             # ### $depth
550             # ### $pow
551             # ### base: ($pow + 1)/2
552             # ### rem n: $n - ($pow + 1)/2
553             #
554             # my @ndigits = digit_split_lowtohigh ($n - ($pow+1)/2, 3);
555             # $#ndigits = $depth-1; # pad to $depth with undefs
556             # ### @ndigits
557             #
558             # return \@ndigits;
559             # }
560             }
561              
562             #------------------------------------------------------------------------------
563             # xy_to_n()
564              
565             # Nrow(depth+1) - Nrow(depth)
566             # = (3*pow+1)/2 - (pow+1)/2
567             # = (3*pow + 1 - pow - 1)/2
568             # = (2*pow)/2
569             # = pow
570             #
571             sub xy_to_n {
572 5213     5213 1 30478 my ($self, $x, $y) = @_;
573 5213         9339 $x = round_nearest ($x);
574 5213         9662 $y = round_nearest ($y);
575             ### PythagoreanTree xy_to_n(): "$x, $y"
576              
577 5213 100       7630 my ($p,$q) = &{$self->{'xy_to_pq'}}($x,$y)
  5213         8914  
578             or return undef; # not a primitive A,B,C
579              
580 1369 100 100     3907 unless ($p >= 2 && $q >= 1) { # must be P > Q >= 1
581 328         577 return undef;
582             }
583 1041 50       2058 if (is_infinite($p)) {
584 0         0 return $p; # infinity
585             }
586 1041 50       2190 if (is_infinite($q)) {
587 0         0 return $q; # infinity
588             }
589 1041 100       2616 if ($p%2 == $q%2) { # must be opposite parity, not same parity
590 480         917 return undef;
591             }
592              
593 561         815 my @ndigits; # low to high
594 561 100       1160 if ($self->{'tree_type'} eq 'FB') {
    50          
595 276         362 for (;;) {
596 885 100 100     2230 unless ($p > $q && $q >= 1) {
597 114         242 return undef;
598             }
599 771 100 100     1627 last if $q <= 1 && $p <= 2;
600              
601 609 100       967 if ($q % 2) {
602             ### q odd, p even, digit 1 or 2 ...
603 323         423 $p /= 2;
604 323 100       519 if ($q > $p) {
605             ### digit 2, M3 ...
606 119         179 push @ndigits, 2;
607 119         167 $q -= $p; # opp parity of p, and < new p
608             } else {
609             ### digit 1, M2 ...
610 204         328 push @ndigits, 1;
611 204         300 $q = $p - $q; # opp parity of p, and < p
612             }
613             } else {
614             ### q even, p odd, digit 0, M1 ...
615 286         424 push @ndigits, 0;
616 286         399 $q /= 2;
617 286         416 $p -= $q; # opp parity of q
618             }
619             ### descend: "$q / $p"
620             }
621              
622             } elsif ($self->{'tree_type'} eq 'UMT') {
623 0         0 for (;;) {
624             ### at: "p=$p q=$q"
625 0         0 my $qmod2 = $q % 2;
626 0 0 0     0 unless ($p > $q && $q >= 1) {
627 0         0 return undef;
628             }
629 0 0 0     0 last if $q <= 1 && $p <= 2;
630              
631 0 0       0 if ($p < 2*$q) {
    0          
632 0         0 ($p,$q) = ($q, 2*$q-$p); # U
633 0         0 push @ndigits, 0;
634             } elsif ($qmod2) {
635 0         0 $p /= 2; # M2
636 0         0 $q = $p - $q;
637 0         0 push @ndigits, 1;
638             } else {
639 0         0 $q /= 2; # T
640 0         0 $p -= 3*$q;
641 0         0 push @ndigits, 2;
642             }
643             }
644              
645             } else {
646             ### UAD or UArD ...
647             ### assert: $self->{'tree_type'} eq 'UAD' || $self->{'tree_type'} eq 'UArD'
648 285         388 for (;;) {
649             ### $p
650             ### $q
651 1065 100 66     3620 if ($q <= 0 || $p <= 0 || $p <= $q) {
      100        
652 119         262 return undef;
653             }
654 946 100 100     1952 last if $q <= 1 && $p <= 2;
655              
656 780 100       1238 if ($p > 2*$q) {
657 317 100       499 if ($p > 3*$q) {
658             ### digit 2 ...
659 230         366 push @ndigits, 2;
660 230         338 $p -= 2*$q;
661             } else {
662             ### digit 1
663 87         125 push @ndigits, 1;
664 87         159 ($p,$q) = ($q, $p - 2*$q);
665             }
666              
667             } else {
668             ### digit 0 ...
669 463         761 push @ndigits, 0;
670 463         801 ($p,$q) = ($q, 2*$q-$p);
671             }
672             ### descend: "$q / $p"
673             }
674             }
675             ### @ndigits
676              
677 328 50       661 if ($self->{'digit_order'} eq 'LtoH') {
678 0         0 @ndigits = reverse @ndigits;
679             ### unreverse: @ndigits
680             }
681 328 50       568 if ($self->{'tree_type'} eq 'UArD') {
682 0         0 Math::PlanePath::GrayCode::_digits_from_gray_reflected(\@ndigits,3);
683             ### ungray: @ndigits
684             }
685              
686 328         546 my $zero = $x*0*$y;
687             ### offset: digit_join_lowtohigh(\@ndigits,3,$zero)
688             ### depth: scalar(@ndigits)
689             ### Nrow: $self->tree_depth_to_n($zero + scalar(@ndigits))
690              
691 328         829 return ($self->tree_depth_to_n($zero + scalar(@ndigits))
692             + digit_join_lowtohigh(\@ndigits,3,$zero)); # offset into row
693             }
694              
695             # numprims(H) = how many with hypot < H
696             # limit H->inf numprims(H) / H -> 1/2pi
697             #
698             # not exact
699             sub rect_to_n_range {
700 64     64 1 5958 my ($self, $x1,$y1, $x2,$y2) = @_;
701             ### PythagoreanTree rect_to_n_range(): "$x1,$y1 $x2,$y2"
702              
703 64         167 $x1 = round_nearest ($x1);
704 64         135 $y1 = round_nearest ($y1);
705 64         124 $x2 = round_nearest ($x2);
706 64         120 $y2 = round_nearest ($y2);
707              
708 64         101 my $zero = ($x1 * 0 * $y1 * $x2 * $y2); # inherit bignum
709              
710 64 50       147 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
711 64 50       126 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
712             ### x2: "$x2"
713             ### y2: "$y2"
714              
715 64 50       148 if ($self->{'coordinates'} eq 'BA') {
716 0         0 ($x2,$y2) = ($y2,$x2);
717             }
718 64 50       137 if ($self->{'coordinates'} eq 'SM') {
719 0 0       0 if ($x2 > $y2) { # both max
720 0         0 $y2 = $x2;
721             } else {
722 0         0 $x2 = $y2;
723             }
724             }
725              
726 64 100       119 if ($self->{'coordinates'} eq 'PQ') {
727 28 50 33     108 if ($x2 < 2 || $y2 < 1) {
728 0         0 return (1,0);
729             }
730             # P > Q so reduce y2 to at most x2-1
731 28 50       51 if ($y2 >= $x2) {
732 0         0 $y2 = $x2-1; # $y2 = min ($y2, $x2-1);
733             }
734              
735 28 50       54 if ($y2 < $y1) {
736             ### PQ y range all above X=Y diagonal ...
737 0         0 return (1,0);
738             }
739             } else {
740             # AB,AC,BC, SM,SC,MC
741 36 50 33     131 if ($x2 < 3 || $y2 < 0) {
742 0         0 return (1,0);
743             }
744             }
745              
746 64         93 my $depth;
747 64 100       120 if ($self->{'tree_type'} eq 'FB') {
748             ### FB ...
749 30 100       59 if ($self->{'coordinates'} eq 'PQ') {
750 14         23 $x2 *= 3;
751             }
752 30         79 my ($pow, $exp) = round_down_pow ($x2, 2);
753 30         48 $depth = 2*$exp;
754             } else {
755             ### UAD or UArD, and UMT ...
756 34 100       63 if ($self->{'coordinates'} eq 'PQ') {
757             ### PQ ...
758             # P=k+1,Q=k diagonal N=100..000 first of row is depth=P-2
759             # anything else in that X=P column is smaller depth
760 14         29 $depth = $x2 - 2;
761             } else {
762 20         49 my $xdepth = int (($x2+1) / 2);
763 20         31 my $ydepth = int (($y2+31) / 4);
764 20         51 $depth = min($xdepth,$ydepth);
765             }
766             }
767             ### depth: "$depth"
768 64         165 return (1, $self->tree_depth_to_n_end($zero+$depth));
769             }
770              
771             #------------------------------------------------------------------------------
772 2     2   23 use constant tree_num_roots => 1;
  2         6  
  2         2943  
773              
774             sub tree_n_children {
775 7     7 1 343 my ($self, $n) = @_;
776 7 50       17 unless ($n >= 1) {
777 0         0 return;
778             }
779 7         10 $n *= 3;
780 7         25 return ($n-1, $n, $n+1);
781             }
782             sub tree_n_num_children {
783 0     0 1 0 my ($self, $n) = @_;
784 0 0       0 return ($n >= 1 ? 3 : undef);
785             }
786             sub tree_n_parent {
787 13     13 1 651 my ($self, $n) = @_;
788 13 100       30 unless ($n >= 2) {
789 1         4 return undef;
790             }
791 12         30 return int(($n+1)/3);
792             }
793             sub tree_n_to_depth {
794 0     0 1 0 my ($self, $n) = @_;
795             ### PythagoreanTree tree_n_to_depth(): $n
796 0 0       0 unless ($n >= 1) {
797 0         0 return undef;
798             }
799 0         0 my ($pow, $depth) = round_down_pow (2*$n-1, 3);
800 0         0 return $depth;
801             }
802              
803             sub tree_depth_to_n {
804 328     328 1 570 my ($self, $depth) = @_;
805 328 50       1185 return ($depth >= 0
806             ? (3**$depth + 1)/2
807             : undef);
808             }
809             # (3^(d+1)+1)/2-1 = (3^(d+1)-1)/2
810             sub tree_depth_to_n_end {
811 64     64 1 110 my ($self, $depth) = @_;
812 64 50       235 return ($depth >= 0
813             ? (3**($depth+1) - 1)/2
814             : undef);
815             }
816             sub tree_depth_to_n_range {
817 0     0 1 0 my ($self, $depth) = @_;
818 0 0       0 if ($depth >= 0) {
819 0         0 my $n_lo = (3**$depth + 1) / 2; # same as tree_depth_to_n()
820 0         0 return ($n_lo, 3*$n_lo-2);
821             } else {
822 0         0 return;
823             }
824             }
825             sub tree_depth_to_width {
826 0     0 1 0 my ($self, $depth) = @_;
827 0 0       0 return ($depth >= 0
828             ? 3**$depth
829             : undef);
830             }
831              
832             #------------------------------------------------------------------------------
833              
834             # Maybe, or abc_to_pq() perhaps with two of three values.
835             #
836             # @EXPORT_OK = ('ab_to_pq','pq_to_ab');
837             #
838             # =item C<($p,$q) = Math::PlanePath::PythagoreanTree::ab_to_pq($a,$b)>
839             #
840             # Return the P,Q coordinates for C<$a,$b>. As described above this is
841             #
842             # P = sqrt((C+A)/2) where C=sqrt(A^2+B^2)
843             # Q = sqrt((C-A)/2)
844             #
845             # The returned P,Q are integers PE=0,QE=0, but the further
846             # conditions for the path (namely PEQE=1 and no common factor) are
847             # not enforced.
848             #
849             # If P,Q are not integers or if BE0 then return an empty list. This
850             # ensures A,B is a Pythagorean triple, ie. that C=sqrt(A^2+B^2) is an
851             # integer, but it might not be a primitive triple and might not have A odd B
852             # even.
853             #
854             # =item C<($a,$b) = Math::PlanePath::PythagoreanTree::pq_to_ab($p,$q)>
855             #
856             # Return the A,B coordinates for C<$p,$q>. This is simply
857             #
858             # $a = $p*$p - $q*$q
859             # $b = 2*$p*$q
860             #
861             # This is intended for use with C<$p,$q> satisfying PEQE=1 and no
862             # common factor, but that's not enforced.
863              
864              
865             # a=p^2-q^2, b=2pq, c=p^2+q^2
866             # Done as a=(p-q)*(p+q) for one multiply instead of two squares, and to work
867             # close to a=UINT_MAX.
868             #
869             sub _pq_to_ab {
870 27     27   45 my ($p, $q) = @_;
871 27         73 return (($p-$q)*($p+$q), 2*$p*$q);
872             }
873              
874             # C=(p-q)^2+B for one squaring instead of two.
875             # Also possible is C=(p+q)^2-B, but prefer "+B" so as not to round-off in
876             # floating point if (p+q)^2 overflows an integer.
877             sub _pq_to_bc {
878 1     1   4 my ($p, $q) = @_;
879 1         3 my $b = 2*$p*$q;
880 1         4 $p -= $q;
881 1         3 return ($b, $p*$p+$b);
882             }
883              
884             # a=p^2-q^2, b=2pq, c=p^2+q^2
885             # Could a=(p-q)*(p+q) to avoid overflow if p^2 exceeds an integer as per
886             # _pq_to_ab(), but c overflows in that case anyway.
887             sub _pq_to_ac {
888 2     2   5 my ($p, $q) = @_;
889 2         6 $p *= $p;
890 2         3 $q *= $q;
891 2         9 return ($p-$q, $p+$q);
892             }
893              
894             # a=p^2-q^2, b=2pq, c=p^2+q^2
895             # a
896             # p^2-q^2 < 2pq
897             # p^2 + 2pq - q^2 < 0
898             # (p+q)^2 - 2*q^2 < 0
899             # (p+q + sqrt(2)*q)*(p+q - sqrt(2)*q) < 0
900             # (p+q - sqrt(2)*q) < 0
901             # p + (1-sqrt(2))*q < 0
902             # p < (sqrt(2)-1)*q
903             #
904             sub _pq_to_sc {
905 0     0   0 my ($p, $q) = @_;
906 0         0 my $b = 2*$p*$q;
907 0         0 my $p_plus_q = $p + $q;
908 0         0 $p -= $q;
909 0         0 return (min($p_plus_q*$p, $b), # A = P^2-Q^2 = (P+Q)*(P-Q)
910             $p*$p+$b); # C = P^2+Q^2 = (P-Q)^2 + 2*P*Q
911             }
912             sub _pq_to_mc {
913 0     0   0 my ($p, $q) = @_;
914 0         0 my $b = 2*$p*$q;
915 0         0 my $p_plus_q = $p + $q;
916 0         0 $p -= $q;
917 0         0 return (max($p_plus_q*$p, $b), # A = P^2-Q^2 = (P+Q)*(P-Q)
918             $p*$p+$b); # C = P^2+Q^2 = (P-Q)^2 + 2*P*Q
919             }
920             sub _pq_to_sm {
921 0     0   0 my ($p, $q) = @_;
922 0         0 my ($a, $b) = _pq_to_ab($p,$q);
923 0 0       0 return ($a < $b ? ($a, $b) : ($b, $a));
924             }
925              
926             # u = p+q, v=p-q
927             # at given p, vertical q
928             # u=p,v=p on diagonal then p+q,p-q is diagonal down
929             # so mirror p axis to x=y diagonal and measure down diagonal from there
930             sub _pq_to_uv {
931 0     0   0 my ($p, $q) = @_;
932 0         0 return ($p+$q, $p-$q);
933             }
934              
935             # r = b+c = 2pq+p^2+q^2 = (p+q)^2
936             # s = c-a = p^2+q^2 - (p^2-q^2) = 2*q^2
937             sub _pq_to_rs {
938 0     0   0 my ($p, $q) = @_;
939 0         0 return (($p+$q)**2, 2*$q*$q);
940             }
941              
942             # s = c-a = p^2+q^2 - (p^2-q^2) = 2*q^2
943             # t = a+b-c = p^2-q^2 + 2pq - (p^2+q^2) = 2pq-2q^2 = 2(p-q)q
944             sub _pq_to_st {
945 0     0   0 my ($p, $q) = @_;
946 0         0 my $q2 = 2*$q;
947 0         0 return ($q2*$q, ($p-$q)*$q2);
948             }
949              
950             #------------------------------------------------------------------------------
951              
952             # a = p^2 - q^2
953             # b = 2pq
954             # c = p^2 + q^2
955             #
956             # q = b/2p
957             # a = p^2 - (b/2p)^2
958             # = p^2 - b^2/4p^2
959             # 4ap^2 = 4p^4 - b^2
960             # 4(p^2)^2 - 4a*p^2 - b^2 = 0
961             # p^2 = [ 4a +/- sqrt(16a^2 + 16*b^2) ] / 2*4
962             # = [ a +/- sqrt(a^2 + b^2) ] / 2
963             # = (a +/- c) / 2 where c=sqrt(a^2+b^2)
964             # p = sqrt((a+c)/2) since c>a
965             #
966             # a = (a+c)/2 - q^2
967             # q^2 = (a+c)/2 - a
968             # = (c-a)/2
969             # q = sqrt((c-a)/2)
970             #
971             # if c^2 = a^2+b^2 is a perfect square then a,b,c is a pythagorean triple
972             # p^2 = (a+c)/2
973             # = (a + sqrt(a^2+b^2))/2
974             # 2p^2 = a + sqrt(a^2+b^2)
975             #
976             # p>q so a>0
977             # a+c even is a odd, c odd or a even, c even
978             # if a odd then c=a^2+b^2 is opp of b parity, must have b even to make c+a even
979             # if a even then c=a^2+b^2 is same as b parity, must have b even to c+a even
980             #
981             # a=6,b=8 is c=sqrt(6^2+8^2)=10
982             # a=0,b=4 is c=sqrt(0+4^4)=4 p^2=(a+c)/2 = 2 not a square
983             # a+c even, then (a+c)/2 == 0,1 mod 4 so a+c==0,2 mod 4
984             #
985             sub _ab_to_pq {
986 5009     5009   48403 my ($a, $b) = @_;
987             ### _ab_to_pq(): "A=$a, B=$b"
988              
989 5009 100 100     15274 unless ($b >= 4 && ($a%2) && !($b%2)) { # A odd, B even
      100        
990 3931         8251 return;
991             }
992              
993             # This used to be $c=hypot($a,$b) and check $c==int($c), but libm hypot()
994             # on Darwin 8.11.0 is somehow a couple of bits off being an integer, for
995             # example hypot(57,176)==185 but a couple of bits out so $c!=int($c).
996             # Would have thought hypot() ought to be exact on integer inputs and a
997             # perfect square sum :-(. Check for a perfect square by multiplying back
998             # instead.
999             #
1000             # The condition is "$csquared != $c*$c" with operands that way around
1001             # since the other way is bad for Math::BigInt::Lite 0.14.
1002             #
1003 1078         2392 my $c;
1004             {
1005 1078         1380 my $csquared = $a*$a + $b*$b;
  1078         1651  
1006 1078         2689 $c = _sqrtint($csquared);
1007             ### $csquared
1008             ### $c
1009             # since A odd and B even should have C odd, but floating point rounding
1010             # might prevent that
1011 1078 100       2222 unless ($csquared == $c*$c) {
1012             ### A^2+B^2 not a perfect square ...
1013 1010         2195 return;
1014             }
1015             }
1016 68         217 return _ac_to_pq($a,$c);
1017             }
1018              
1019             sub _bc_to_pqa {
1020 1290     1290   2063 my ($b, $c) = @_;
1021             ### _bc_to_pqa(): "B=$b C=$c"
1022              
1023 1290 100 100     3325 unless ($c > $b && $b >= 4 && !($b%2) && ($c%2)) { # B even, C odd
      100        
      100        
1024 1216         3089 return;
1025             }
1026              
1027 74         100 my $a;
1028             {
1029 74         101 my $asquared = $c*$c - $b*$b;
  74         120  
1030 74 50       122 unless ($asquared > 0) {
1031 0         0 return;
1032             }
1033 74         145 $a = _sqrtint($asquared);
1034             ### $asquared
1035             ### $a
1036 74 100       152 unless ($asquared == $a*$a) {
1037 64         185 return;
1038             }
1039             }
1040              
1041             # If $c is near DBL_MAX can have $a overflow to infinity, leaving A>C.
1042             # _ac_to_pq() will detect that.
1043 10 100       23 my ($p,$q) = _ac_to_pq($a,$c) or return;
1044 8         23 return ($p,$q,$a);
1045             }
1046              
1047             sub _ac_to_pq {
1048 1369     1369   2130 my ($a, $c) = @_;
1049             ### _ac_to_pq(): "A=$a C=$c"
1050              
1051 1369 100 100     3841 unless ($c > $a && $a >= 3 && ($a%2) && ($c%2)) { # A odd, C odd
      100        
      100        
1052 1224         3122 return;
1053             }
1054 145         676 $a = ($a-1)/2;
1055 145         576 $c = ($c-1)/2;
1056             ### halved to: "a=$a c=$c"
1057              
1058 145         499 my $p;
1059             {
1060             # If a,b,c is a triple but not primitive then can have psquared not an
1061             # integer. Eg. a=9,b=12 has c=15 giving psquared=(9+15)/2=12 is not a
1062             # perfect square. So notice that here.
1063             #
1064 145         191 my $psquared = $c+$a+1;
  145         256  
1065 145         487 $p = _sqrtint($psquared);
1066             ### $psquared
1067             ### $p
1068 145 100       376 unless ($psquared == $p*$p) {
1069             ### P^2=A+C not a perfect square ...
1070 72         191 return;
1071             }
1072             }
1073              
1074 73         219 my $q;
1075             {
1076             # If a,b,c is a triple but not primitive then can have qsquared not an
1077             # integer. Eg. a=15,b=36 has c=39 giving qsquared=(39-15)/2=12 is not a
1078             # perfect square. So notice that here.
1079             #
1080 73         128 my $qsquared = $c-$a;
  73         120  
1081 73         222 $q = _sqrtint($qsquared);
1082             ### $qsquared
1083             ### $q
1084 73 100       245 unless ($qsquared == $q*$q) {
1085 6         19 return;
1086             }
1087             }
1088              
1089             # Might have a common factor between P,Q here. Eg.
1090             # A=27 = 3*3*3, B=36 = 4*3*3
1091             # A=45 = 3*3*5, B=108 = 4*3*3*3
1092             # A=63, B=216
1093             # A=75 =3*5*5 B=100 = 4*5*5
1094             # A=81, B=360
1095             #
1096 67         256 return ($p, $q);
1097             }
1098              
1099             sub _sm_to_pq {
1100 0     0   0 my ($s, $m) = @_;
1101 0 0       0 unless ($s < $m) {
1102 0         0 return;
1103             }
1104 0 0       0 return _ab_to_pq($s % 2
1105             ? ($s,$m) # s odd is A
1106             : ($m,$s)); # s even is B
1107             }
1108              
1109              
1110             # s^2+m^2=c^2
1111             # if s odd then a=s
1112             # ac_to_pq
1113             # b = 2pq check isn't smaller than s
1114             #
1115             # p^2=(c+a)/2
1116             # q^2=(c-a)/2
1117              
1118             sub _sc_to_pq {
1119 2     2   203 my ($s, $c) = @_;
1120 2         4 my ($p,$q);
1121 2 100       6 if ($s % 2) {
1122 1 50       4 ($p,$q) = _ac_to_pq($s,$c) # s odd is A
1123             or return;
1124 1 50       4 if ($s > 2*$p*$q) { return; } # if s>B then s is not the smaller one
  0         0  
1125             } else {
1126 1 50       4 ($p,$q,$a) = _bc_to_pqa($s,$c) # s even is B
1127             or return;
1128 1 50       3 if ($s > $a) { return; } # if s>A then s is not the smaller one
  1         3  
1129             }
1130 1         3 return ($p,$q);
1131             }
1132              
1133             sub _mc_to_pq {
1134 0     0     my ($m, $c) = @_;
1135             ### _mc_to_pq() ...
1136 0           my ($p,$q);
1137 0 0         if ($m % 2) {
1138             ### m odd is A ...
1139 0 0         ($p,$q) = _ac_to_pq($m,$c)
1140             or return;
1141 0 0         if ($m < 2*$p*$q) { return; } # if m
  0            
1142             } else {
1143             ### m even is B ...
1144 0 0         ($p,$q,$a) = _bc_to_pqa($m,$c)
1145             or return;
1146             ### $a
1147 0 0         if ($m < $a) { return; } # if m
  0            
1148             }
1149 0           return ($p,$q);
1150             }
1151              
1152             # u = p+q, v=p-q
1153             # u+v=2p p = (u+v)/2
1154             # u-v=2q q = (u-v)/2
1155             sub _uv_to_pq {
1156 0     0     my ($u, $v) = @_;
1157 0           return (($u+$v)/2, ($u-$v)/2);
1158             }
1159              
1160             # r = (p+q)^2
1161             # s = 2*q^2 so q = sqrt(r/2)
1162             sub _rs_to_pq {
1163 0     0     my ($r, $s) = @_;
1164              
1165 0 0         return if $s % 2;
1166 0           $s /= 2;
1167 0 0         return unless $s >= 1;
1168 0           my $q = _sqrtint($s);
1169 0 0         return unless $q*$q == $s;
1170              
1171 0 0         return unless $r >= 1;
1172 0           my $p_plus_q = _sqrtint($r);
1173 0 0         return unless $p_plus_q*$p_plus_q == $r;
1174              
1175 0           return ($p_plus_q - $q, $q);
1176             }
1177              
1178             # s = 2*q^2
1179             # t = a+b-c = p^2-q^2 + 2pq - (p^2+q^2) = 2pq-2q^2 = 2(p-q)q
1180             #
1181             # p=2,q=1 s=2 t=2.1.1=2
1182             #
1183             sub _st_to_pq {
1184 0     0     my ($s, $t) = @_;
1185              
1186             ### _st_to_pq(): "$s, $t"
1187 0 0         return if $s % 2;
1188 0           $s /= 2;
1189 0 0         return unless $s >= 1;
1190 0           my $q = _sqrtint($s);
1191             ### $q
1192 0 0         return unless $q*$q == $s;
1193              
1194 0 0         return if $t % 2;
1195 0           $t /= 2;
1196             ### rem: $t % $q
1197 0 0         return if $t % $q;
1198 0           $t /= $q; # p-q
1199              
1200             ### pq: ($t+$q).", $q"
1201              
1202 0           return ($t+$q, $q);
1203             }
1204              
1205             1;
1206             __END__