File Coverage

blib/lib/Math/PlanePath/Hypot.pm
Criterion Covered Total %
statement 121 186 65.0
branch 19 50 38.0
condition 5 17 29.4
subroutine 11 18 61.1
pod 9 9 100.0
total 165 280 58.9


line stmt bran cond sub pod time code
1             # Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020 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              
20             # A000328 Number of points of norm <= n^2 in square lattice.
21             # 1, 5, 13, 29, 49, 81, 113, 149, 197, 253, 317, 377, 441, 529, 613, 709, 797
22             # a(n) = 1 + 4 * sum(j=0, n^2 / 4, n^2 / (4*j+1) - n^2 / (4*j+3) )
23             # A014200 num points norm <= n^2, excluding 0, divided by 4
24             #
25             # A046109 num points norm == n^2
26             #
27             # A057655 num points x^2+y^2 <= n
28             # A014198 = A057655 - 1
29             #
30             # A004018 num points x^2+y^2 == n
31             #
32             # A057962 hypot count x-1/2,y-1/2 <= n
33             # is last point of each hypot in points=odd
34             #
35             # A057961 hypot count as radius increases
36             #
37              
38             # points="square_horiz"
39             # points="square_vert"
40             # points="square_centre"
41             # A199015 square_centred partial sums
42             #
43              
44              
45             package Math::PlanePath::Hypot;
46 1     1   1211 use 5.004;
  1         4  
47 1     1   8 use strict;
  1         7  
  1         25  
48 1     1   5 use Carp 'croak';
  1         2  
  1         68  
49              
50 1     1   7 use vars '$VERSION', '@ISA';
  1         2  
  1         69  
51             $VERSION = 129;
52 1     1   796 use Math::PlanePath;
  1         3  
  1         41  
53             @ISA = ('Math::PlanePath');
54              
55             use Math::PlanePath::Base::Generic
56 1         69 'is_infinite',
57 1     1   7 'round_nearest';
  1         2  
58              
59             # uncomment this to run the ### lines
60             # use Smart::Comments;
61              
62              
63 1         250 use constant parameter_info_array =>
64             [ { name => 'points',
65             share_key => 'points_aeo',
66             display => 'Points',
67             type => 'enum',
68             default => 'all',
69             choices => ['all','even','odd'],
70             choices_display => ['All','Even','Odd'],
71             description => 'Which X,Y points visit, either all of them or just X+Y=even or odd.',
72             },
73             Math::PlanePath::Base::Generic::parameter_info_nstart1(),
74 1     1   5 ];
  1         2  
75              
76             {
77             my %x_negative_at_n = (all => 3,
78             even => 2,
79             odd => 2);
80             sub x_negative_at_n {
81 0     0 1 0 my ($self) = @_;
82 0         0 return $self->n_start + $x_negative_at_n{$self->{'points'}};
83             }
84             }
85             {
86             my %y_negative_at_n = (all => 4,
87             even => 3,
88             odd => 3);
89             sub y_negative_at_n {
90 0     0 1 0 my ($self) = @_;
91 0         0 return $self->n_start + $y_negative_at_n{$self->{'points'}};
92             }
93             }
94             sub rsquared_minimum {
95 0     0 1 0 my ($self) = @_;
96 0 0       0 return ($self->{'points'} eq 'odd'
97             ? 1 # odd at X=1,Y=0
98             : 0); # even,all at X=0,Y=0
99             }
100             # points=even includes X=Y so abs(X-Y)>=0
101             # points=odd doesn't include X=Y so abs(X-Y)>=1
102             *absdiffxy_minimum = \&rsquared_minimum;
103             *sumabsxy_minimum = \&rsquared_minimum;
104              
105 1     1   8 use constant turn_any_right => 0; # always left or straight
  1         11  
  1         1491  
106             sub turn_any_straight {
107 0     0 1 0 my ($self) = @_;
108 0         0 return ($self->{'points'} ne 'all'); # points=all is left always
109             }
110              
111              
112             #------------------------------------------------------------------------------
113              
114             sub new {
115 10     10 1 567 my $self = shift->SUPER::new(@_);
116              
117 10 100       43 if (! defined $self->{'n_start'}) {
118 1         10 $self->{'n_start'} = $self->default_n_start;
119             }
120              
121 10   100     40 my $points = ($self->{'points'} ||= 'all');
122 10 100       43 if ($points eq 'all') {
    100          
    50          
    0          
123 4         14 $self->{'n_to_x'} = [0];
124 4         9 $self->{'n_to_y'} = [0];
125 4         9 $self->{'hypot_to_n'} = [0];
126 4         8 $self->{'y_next_x'} = [1, 1];
127 4         9 $self->{'y_next_hypot'} = [1, 2];
128 4         11 $self->{'x_inc'} = 1;
129 4         7 $self->{'x_inc_factor'} = 2;
130 4         9 $self->{'x_inc_squared'} = 1;
131 4         7 $self->{'y_factor'} = 2;
132 4         8 $self->{'opposite_parity'} = -1;
133              
134             } elsif ($points eq 'even') {
135 3         11 $self->{'n_to_x'} = [0];
136 3         7 $self->{'n_to_y'} = [0];
137 3         8 $self->{'hypot_to_n'} = [0];
138 3         7 $self->{'y_next_x'} = [2, 1];
139 3         5 $self->{'y_next_hypot'} = [4, 2];
140 3         10 $self->{'x_inc'} = 2;
141 3         6 $self->{'x_inc_factor'} = 4;
142 3         7 $self->{'x_inc_squared'} = 4;
143 3         4 $self->{'y_factor'} = 2;
144 3         7 $self->{'opposite_parity'} = 1;
145              
146             } elsif ($points eq 'odd') {
147 3         10 $self->{'n_to_x'} = [];
148 3         8 $self->{'n_to_y'} = [];
149 3         8 $self->{'hypot_to_n'} = [];
150 3         10 $self->{'y_next_x'} = [1];
151 3         7 $self->{'y_next_hypot'} = [1];
152 3         13 $self->{'x_inc'} = 2;
153 3         6 $self->{'x_inc_factor'} = 4;
154 3         5 $self->{'x_inc_squared'} = 4;
155 3         6 $self->{'y_factor'} = 2;
156 3         7 $self->{'opposite_parity'} = 0;
157              
158             } elsif ($points eq 'square_centred') {
159 0         0 $self->{'n_to_x'} = [];
160 0         0 $self->{'n_to_y'} = [];
161 0         0 $self->{'hypot_to_n'} = [];
162 0         0 $self->{'y_next_x'} = [undef,1];
163 0         0 $self->{'y_next_hypot'} = [undef,2];
164 0         0 $self->{'x_inc'} = 2;
165 0         0 $self->{'x_inc_factor'} = 4; # ((x+2)^2 - x^2) = 4*x+4
166 0         0 $self->{'x_inc_squared'} = 4;
167 0         0 $self->{'y_start'} = 1;
168 0         0 $self->{'y_inc'} = 2;
169 0         0 $self->{'opposite_parity'} = -1;
170              
171             } else {
172 0         0 croak "Unrecognised points option: ", $points;
173             }
174 10         27 return $self;
175             }
176              
177             sub _extend {
178 1017     1017   1575 my ($self) = @_;
179             ### _extend() n: scalar(@{$self->{'n_to_x'}})
180             ### y_next_x: $self->{'y_next_x'}
181              
182 1017         1356 my $n_to_x = $self->{'n_to_x'};
183 1017         1400 my $n_to_y = $self->{'n_to_y'};
184 1017         1375 my $hypot_to_n = $self->{'hypot_to_n'};
185 1017         1376 my $y_next_x = $self->{'y_next_x'};
186 1017         1356 my $y_next_hypot = $self->{'y_next_hypot'};
187 1017   50     2361 my $y_start = $self->{'y_start'} || 0;
188 1017   50     2245 my $y_inc = $self->{'y_inc'} || 1;
189              
190             # set @y to the Y with the smallest $y_next_hypot[$y], and if there's some
191             # Y's with equal smallest hypot then all those Y's
192 1017         1512 my @y = ($y_start);
193 1017   50     1825 my $hypot = $y_next_hypot->[$y_start] || 99;
194 1017         1916 for (my $y = $y_start+$y_inc; $y < @$y_next_x; $y += $y_inc) {
195 10749 100       24531 if ($hypot == $y_next_hypot->[$y]) {
    100          
196 324         613 push @y, $y;
197             } elsif ($hypot > $y_next_hypot->[$y]) {
198 1512         2240 @y = ($y);
199 1512         2743 $hypot = $y_next_hypot->[$y];
200             }
201             }
202              
203             ### chosen y list: @y
204              
205             # if the endmost of the @$y_next_x, @$y_next_hypot arrays are used then
206             # extend them by one
207 1017 100       1892 if ($y[-1] == $#$y_next_x) {
208             ### grow y_next_x ...
209 141         209 my $y = $#$y_next_x + $y_inc;
210 141         224 my $x = $y + ($self->{'points'} eq 'odd');
211 141         245 $y_next_x->[$y] = $x;
212 141         246 $y_next_hypot->[$y] = $x*$x+$y*$y;
213             ### $y_next_x
214             ### $y_next_hypot
215             ### assert: $y_next_hypot->[$y] == $y**2 + $x*$x
216             }
217              
218             # @x is the $y_next_x[$y] for each of the @y smallests, and step those
219             # selected elements next X and hypot for that new X,Y
220             my @x = map {
221 1017         1644 my $y = $_;
  1236         1575  
222 1236         1626 my $x = $y_next_x->[$y];
223 1236         1645 $y_next_x->[$y] += $self->{'x_inc'};
224             $y_next_hypot->[$y]
225 1236         1853 += $self->{'x_inc_factor'} * $x + $self->{'x_inc_squared'};
226             ### assert: $y_next_hypot->[$y] == ($x+$self->{'x_inc'})**2 + $y**2
227 1236         2451 $x
228             } @y;
229             ### $hypot
230             ### base octant: join(' ',map{"$x[$_],$y[$_]"} 0 .. $#x)
231              
232             # transpose X,Y to Y,X
233             {
234 1017         1438 my @base_x = @x;
235 1017         1412 my @base_y = @y;
236 1017 100       1822 unless ($y[0]) { # no transpose of x,0
237 126         172 shift @base_x;
238 126         169 shift @base_y;
239             }
240 1017 100       1762 if ($x[-1] == $y[-1]) { # no transpose of x,x
241 87         118 pop @base_x;
242 87         109 pop @base_y;
243             }
244 1017         1527 push @x, reverse @base_y;
245 1017         1613 push @y, reverse @base_x;
246             }
247             ### with transpose q1: join(' ',map{"$x[$_],$y[$_]"} 0 .. $#x)
248              
249             # rotate +90 quadrant 1 into quadrant 2
250             {
251 1017         1347 my @base_y = @y;
  1017         1327  
  1017         1442  
252 1017         1469 push @y, @x;
253 1017         1438 push @x, map {-$_} @base_y;
  2259         3716  
254             }
255             ### with rotate q2: join(' ',map{"$x[$_],$y[$_]"} 0 .. $#x)
256              
257             # rotate +180 quadrants 1+2 into quadrants 2+3
258 1017         1564 push @x, map {-$_} @x;
  4518         6183  
259 1017         1540 push @y, map {-$_} @y;
  4518         6167  
260              
261             ### store: join(' ',map{"$x[$_],$y[$_]"} 0 .. $#x)
262             ### at n: scalar(@$n_to_x)
263             ### hypot_to_n: "h=$hypot n=".scalar(@$n_to_x)
264 1017         1844 $hypot_to_n->[$hypot] = scalar(@$n_to_x);
265 1017         2374 push @$n_to_x, @x;
266 1017         3766 push @$n_to_y, @y;
267              
268             # ### hypot_to_n now: join(' ',map {defined($hypot_to_n->[$_]) && "h=$_,n=$hypot_to_n->[$_]"} 0 .. $#$hypot_to_n)
269              
270              
271             # my $x = $y_next_x->[0];
272             #
273             # $x = $y_next_x->[$y];
274             # $n_to_x->[$next_n] = $x;
275             # $n_to_y->[$next_n] = $y;
276             # $xy_to_n{"$x,$y"} = $next_n++;
277             #
278             # $y_next_x->[$y]++;
279             # $y_next_hypot->[$y] = $y*$y + $y_next_x->[$y]**2;
280             }
281              
282             sub n_to_xy {
283 9009     9009 1 106381 my ($self, $n) = @_;
284             ### Hypot n_to_xy(): $n
285              
286 9009         12569 $n = $n - $self->{'n_start'}; # starting $n==0, warn if $n==undef
287 9009 50       15353 if ($n < 0) { return; }
  0         0  
288 9009 50       15207 if (is_infinite($n)) { return ($n,$n); }
  0         0  
289              
290 9009         14432 my $int = int($n);
291 9009         11185 $n -= $int; # fraction part
292              
293 9009         12520 my $n_to_x = $self->{'n_to_x'};
294 9009         11649 my $n_to_y = $self->{'n_to_y'};
295              
296 9009         16011 while ($int >= $#$n_to_x) {
297 1017         1585 _extend($self);
298             }
299              
300 9009         12289 my $x = $n_to_x->[$int];
301 9009         11344 my $y = $n_to_y->[$int];
302 9009         21480 return ($x + $n * ($n_to_x->[$int+1] - $x),
303             $y + $n * ($n_to_y->[$int+1] - $y));
304             }
305              
306             sub xy_is_visited {
307 0     0 1   my ($self, $x, $y) = @_;
308              
309 0 0         if ($self->{'opposite_parity'} >= 0) {
310 0           $x = round_nearest ($x);
311 0           $y = round_nearest ($y);
312 0 0         if ((($x%2) ^ ($y%2)) == $self->{'opposite_parity'}) {
313 0           return 0;
314             }
315             }
316 0 0         if ($self->{'points'} eq 'square_centred') {
317 0 0 0       unless (($y%2) && ($x%2)) {
318 0           return 0;
319             }
320             }
321 0           return 1;
322             }
323              
324             sub xy_to_n {
325 0     0 1   my ($self, $x, $y) = @_;
326             ### Hypot xy_to_n(): "$x, $y"
327             ### hypot_to_n last: $#{$self->{'hypot_to_n'}}
328              
329 0           $x = round_nearest ($x);
330 0           $y = round_nearest ($y);
331              
332 0 0         if ((($x%2) ^ ($y%2)) == $self->{'opposite_parity'}) {
333 0           return undef;
334             }
335 0 0         if ($self->{'points'} eq 'square_centred') {
336 0 0 0       unless (($y%2) && ($x%2)) {
337 0           return undef;
338             }
339             }
340              
341 0           my $hypot = $x*$x + $y*$y;
342 0 0         if (is_infinite($hypot)) {
343             ### infinity
344 0           return undef;
345             }
346              
347 0           my $n_to_x = $self->{'n_to_x'};
348 0           my $n_to_y = $self->{'n_to_y'};
349              
350 0           my $hypot_to_n = $self->{'hypot_to_n'};
351 0           while ($hypot > $#$hypot_to_n) {
352 0           _extend($self);
353             }
354              
355 0           my $n = $hypot_to_n->[$hypot];
356 0           for (;;) {
357 0 0 0       if ($x == $n_to_x->[$n] && $y == $n_to_y->[$n]) {
358 0           return $n + $self->{'n_start'};
359             }
360 0           $n += 1;
361              
362 0 0         if ($n_to_x->[$n]**2 + $n_to_y->[$n]**2 != $hypot) {
363             ### oops, hypot_to_n no good ...
364 0           return undef;
365             }
366             }
367              
368             # if ($x < 0 || $y < 0) {
369             # return undef;
370             # }
371             # my $h = $x*$x + $y*$y;
372             #
373             # while ($y_next_x[$y] <= $x) {
374             # _extend($self);
375             # }
376             # return $xy_to_n{"$x,$y"};
377             }
378              
379             # not exact
380             sub rect_to_n_range {
381 0     0 1   my ($self, $x1,$y1, $x2,$y2) = @_;
382              
383 0           $x1 = abs (round_nearest ($x1));
384 0           $y1 = abs (round_nearest ($y1));
385 0           $x2 = abs (round_nearest ($x2));
386 0           $y2 = abs (round_nearest ($y2));
387              
388 0 0         if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  0            
389 0 0         if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }
  0            
390              
391             # circle area pi*r^2, with r^2 = $x2**2 + $y2**2
392             return ($self->{'n_start'},
393 0           $self->{'n_start'} + int (3.2 * (($x2+1)**2 + ($y2+1)**2)));
394             }
395              
396             1;
397             __END__