File Coverage

blib/lib/Math/PlanePath/HypotOctant.pm
Criterion Covered Total %
statement 98 145 67.5
branch 16 40 40.0
condition 2 14 14.2
subroutine 14 18 77.7
pod 7 7 100.0
total 137 224 61.1


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             # Circle drop splash rings from
20             # math-image --path=HypotOctant --values=DigitProductSteps,values_type=count
21             # math-image --path=Hypot --values=DigitProduct
22             # math-image --path=Hypot --values=DigitCount
23             # math-image --path=Hypot --values=Modulo,modulus=1000
24             # http://stefan.guninski.com/oeisposter/
25             #
26             # pi*r^2 - pi*(r-1)^2 = pi*(2r-1)
27             # octant is 1/8 of that pi*(2x-1)/8
28             # pi*(2x-1)/8=100k
29             # 2x-1 = 100k*8/pi
30             # x = 100*4/pi*k
31             #
32             # A000328 Number of points of norm <= n^2 in square lattice.
33             # 1, 5, 13, 29, 49, 81, 113, 149, 197, 253, 317, 377, 441, 529, 613, 709, 797
34             # a(n) = 1 + 4 * sum(j=0, n^2 / 4, n^2 / (4*j+1) - n^2 / (4*j+3) )
35             #
36             # A057655 num points norm <= n in square lattice.
37             #
38             # A036702 num points |z=a+bi| <= n with 0<=a, 0<=b<=a, so octant
39             # A036703 num points n-1 < z <= n, first diffs?
40              
41              
42              
43             package Math::PlanePath::HypotOctant;
44 1     1   7814 use 5.004;
  1         10  
45 1     1   5 use strict;
  1         2  
  1         26  
46 1     1   5 use Carp 'croak';
  1         2  
  1         39  
47              
48 1     1   5 use vars '$VERSION', '@ISA';
  1         1  
  1         53  
49             $VERSION = 127;
50 1     1   588 use Math::PlanePath;
  1         2  
  1         35  
51             @ISA = ('Math::PlanePath');
52              
53             use Math::PlanePath::Base::Generic
54 1         53 'is_infinite',
55 1     1   6 'round_nearest';
  1         1  
56              
57             # uncomment this to run the ### lines
58             #use Smart::Comments;
59              
60              
61 1         44 use constant parameter_info_array =>
62             [ { name => 'points',
63             share_key => 'points_aeo',
64             display => 'Points',
65             type => 'enum',
66             default => 'all',
67             choices => ['all','even','odd'],
68             choices_display => ['All','Even','Odd'],
69             description => 'Which X,Y points visit, either all of them or just X+Y even or X+Y odd.',
70             },
71 1     1   5 ];
  1         1  
72              
73 1     1   5 use constant class_x_negative => 0;
  1         1  
  1         32  
74 1     1   5 use constant class_y_negative => 0;
  1         2  
  1         167  
75              
76             sub x_minimum {
77 6     6 1 12 my ($self) = @_;
78 6 100       25 return ($self->{'points'} eq 'odd'
79             ? 1 # odd, line X=Y not included
80             : 0); # octant Y<=X so X-Y>=0
81             }
82             # points=odd X=1,Y=0
83             # otherwise X=0,Y=0
84             *sumabsxy_minimum = \&x_minimum;
85             *diffxy_minimum = \&x_minimum; # X>=Y so X-Y>=0
86             *absdiffxy_minimum = \&x_minimum;
87             *rsquared_minimum = \&x_minimum;
88              
89             sub absdy_minimum {
90 0     0 1 0 my ($self) = @_;
91 0 0       0 return ($self->{'points'} eq 'all'
92             ? 0
93             : 1); # never same Y
94             }
95              
96             sub dir_minimum_dxdy {
97 0     0 1 0 my ($self) = @_;
98 0 0       0 return ($self->{'points'} eq 'all'
99             ? (1,0) # all i=1 to X=1,Y=0
100             : (1,1)); # odd,even always at least NE
101             }
102             # max direction SE diagonal as anything else is at most tangent to the
103             # eighth of a circle
104 1     1   5 use constant dir_maximum_dxdy => (1,-1); # South-East
  1         2  
  1         906  
105              
106              
107             #------------------------------------------------------------------------------
108              
109             # my @n_to_x = (undef, 0);
110             # my @n_to_y = (undef, 0);
111             # my @hypot_to_n = (1);
112             # my @y_next_x = (1, 1);
113             # my @y_next_hypot = (1, 2);
114              
115             sub new {
116 6     6 1 1134 my $self = shift->SUPER::new(@_);
117              
118 6   100     28 my $points = ($self->{'points'} ||= 'all');
119 6 100       25 if ($points eq 'all') {
    100          
    50          
120 3         6 $self->{'n_to_x'} = [undef];
121 3         6 $self->{'n_to_y'} = [undef];
122 3         5 $self->{'hypot_to_n'} = [];
123 3         6 $self->{'y_next_x'} = [0];
124 3         5 $self->{'y_next_hypot'} = [0];
125 3         5 $self->{'x_inc'} = 1;
126 3         6 $self->{'x_inc_factor'} = 2;
127 3         4 $self->{'x_inc_squared'} = 1;
128 3         6 $self->{'opposite_parity'} = -1;
129              
130             } elsif ($points eq 'even') {
131 1         5 $self->{'n_to_x'} = [undef, 0];
132 1         3 $self->{'n_to_y'} = [undef, 0];
133 1         3 $self->{'hypot_to_n'} = [1];
134 1         3 $self->{'y_next_x'} = [2, 1];
135 1         4 $self->{'y_next_hypot'} = [4, 2];
136 1         2 $self->{'x_inc'} = 2;
137 1         4 $self->{'x_inc_factor'} = 4;
138 1         3 $self->{'x_inc_squared'} = 4;
139 1         3 $self->{'opposite_parity'} = 1;
140              
141             } elsif ($points eq 'odd') {
142 2         4 $self->{'n_to_x'} = [undef];
143 2         5 $self->{'n_to_y'} = [undef];
144 2         5 $self->{'hypot_to_n'} = [undef];
145 2         5 $self->{'y_next_x'} = [1];
146 2         5 $self->{'y_next_hypot'} = [1];
147 2         5 $self->{'x_inc'} = 2;
148 2         7 $self->{'x_inc_factor'} = 4;
149 2         3 $self->{'x_inc_squared'} = 4;
150 2         4 $self->{'opposite_parity'} = 0;
151              
152             } else {
153 0         0 croak "Unrecognised points option: ", $points;
154             }
155 6         12 return $self;
156             }
157              
158              
159             # at h=x^2+y^2
160             # step to (x+k)^2+y^2
161             # is add 2*x*k+k*k
162              
163             sub _extend {
164 2173     2173   2660 my ($self) = @_;
165             ### _extend() n: scalar(@{$self->{'n_to_x'}})
166              
167 2173         2442 my $n_to_x = $self->{'n_to_x'};
168 2173         2265 my $n_to_y = $self->{'n_to_y'};
169 2173         2302 my $hypot_to_n = $self->{'hypot_to_n'};
170 2173         2258 my $y_next_x = $self->{'y_next_x'};
171 2173         2360 my $y_next_hypot = $self->{'y_next_hypot'};
172              
173 2173         2642 my @y = (0);
174 2173         2367 my $hypot = $y_next_hypot->[0];
175 2173         3165 for (my $i = 1; $i < @$y_next_x; $i++) {
176 63482 100       115625 if ($hypot == $y_next_hypot->[$i]) {
    100          
177 1157         1901 push @y, $i;
178             } elsif ($hypot > $y_next_hypot->[$i]) {
179 4158         4899 @y = ($i);
180 4158         6239 $hypot = $y_next_hypot->[$i];
181             }
182             }
183              
184 2173 100       3044 if ($y[-1] == $#$y_next_x) {
185 134         149 my $y = scalar(@$y_next_x);
186 134         190 my $x = $y + ($self->{'points'} eq 'odd');
187 134         178 $y_next_x->[$y] = $x;
188 134         186 $y_next_hypot->[$y] = $x*$x+$y*$y;
189             ### assert: $y_next_hypot->[$y] == $y**2 + $y_next_x->[$y]**2
190             }
191              
192             ### store: join(' ',map{"$n_to_x->[$_],$n_to_y->[$_]"} 0 .. $#$n_to_x)
193             ### at n: scalar(@$n_to_x)
194             ### hypot_to_n: "h=$hypot n=".scalar(@$n_to_x)
195              
196 2173         3123 $hypot_to_n->[$hypot] = scalar(@$n_to_x);
197 2173         2803 push @$n_to_y, @y;
198             push @$n_to_x,
199             map {
200 2173         2676 my $x = $y_next_x->[$_];
  2999         3291  
201 2999         3378 $y_next_x->[$_] += $self->{'x_inc'};
202             $y_next_hypot->[$_]
203 2999         3662 += $self->{'x_inc_factor'} * $x + $self->{'x_inc_squared'};
204             ### assert: $y_next_hypot->[$_] == $_**2 + $y_next_x->[$_]**2
205 2999         6213 $x
206             } @y;
207              
208             # ### hypot_to_n now: join(' ',map {defined($hypot_to_n->[$_]) && "h=$_,n=$hypot_to_n->[$_]"} 0 .. $#$hypot_to_n)
209             }
210              
211             sub n_to_xy {
212 3000     3000 1 17590 my ($self, $n) = @_;
213             ### Hypot n_to_xy(): $n
214              
215 3000 50       4173 if ($n < 1) { return; }
  0         0  
216 3000 50       4341 if (is_infinite($n)) { return ($n,$n); }
  0         0  
217              
218             {
219 3000         3811 my $int = int($n);
  3000         3298  
220 3000 50       4260 if ($n != $int) {
221 0         0 my $frac = $n - $int; # inherit possible BigFloat/BigRat
222 0         0 my ($x1,$y1) = $self->n_to_xy($int);
223 0         0 my ($x2,$y2) = $self->n_to_xy($int+1);
224 0         0 my $dx = $x2-$x1;
225 0         0 my $dy = $y2-$y1;
226 0         0 return ($frac*$dx + $x1, $frac*$dy + $y1);
227             }
228             }
229              
230 3000         3613 my $n_to_x = $self->{'n_to_x'};
231 3000         3176 my $n_to_y = $self->{'n_to_y'};
232              
233 3000         4201 while ($n > $#$n_to_x) {
234 2173         2749 _extend($self);
235             }
236              
237 3000         5461 return ($n_to_x->[$n], $n_to_y->[$n]);
238             }
239              
240             sub xy_to_n {
241 0     0 1   my ($self, $x, $y) = @_;
242             ### Hypot xy_to_n(): "$x, $y"
243             ### hypot_to_n last: $#{$self->{'hypot_to_n'}}
244              
245 0           $x = round_nearest ($x);
246 0           $y = round_nearest ($y);
247              
248 0 0         if ((($x%2) ^ ($y%2)) == $self->{'opposite_parity'}) {
249 0           return undef;
250             }
251              
252 0           my $hypot = $x*$x + $y*$y;
253 0 0         if (is_infinite($hypot)) {
254 0           return $hypot;
255             }
256              
257 0 0 0       if ($x < 0 || $y < 0 || $y > $x) {
      0        
258             ### outside first octant ...
259 0           return undef;
260             }
261              
262 0           my $hypot_to_n = $self->{'hypot_to_n'};
263 0           while ($hypot > $#$hypot_to_n) {
264 0           _extend($self);
265             }
266              
267 0           my $n_to_x = $self->{'n_to_x'};
268 0           my $n_to_y = $self->{'n_to_y'};
269              
270 0           my $n = $hypot_to_n->[$hypot];
271 0           for (;;) {
272 0 0 0       if ($x == $n_to_x->[$n] && $y == $n_to_y->[$n]) {
273 0           return $n;
274             }
275 0           $n += 1;
276              
277 0 0         if ($n_to_x->[$n]**2 + $n_to_y->[$n]**2 != $hypot) {
278             ### oops, hypot_to_n no good ...
279 0           return undef;
280             }
281             }
282             }
283              
284             # not exact
285             sub rect_to_n_range {
286 0     0 1   my ($self, $x1,$y1, $x2,$y2) = @_;
287              
288 0           $x1 = round_nearest ($x1);
289 0           $y1 = round_nearest ($y1);
290 0           $x2 = round_nearest ($x2);
291 0           $y2 = round_nearest ($y2);
292 0 0         if ($x1 > $x2) { ($x1,$x2) = ($x2,$x1); }
  0            
293 0 0         if ($y1 > $y2) { ($y1,$y2) = ($y2,$y1); }
  0            
294              
295 0 0 0       if ($x2 < 0 || $y2 < 0) {
296 0           return (1, 0);
297             }
298              
299             # circle area pi*r^2, with r^2 = $x2**2 + $y2**2
300 0           return (1, 1 + int (3.2/8 * (($x2+1)**2 + ($y2+1)**2)));
301             }
302              
303             1;
304             __END__