File Coverage

blib/lib/Math/PlanePath/FilledRings.pm
Criterion Covered Total %
statement 150 181 82.8
branch 34 56 60.7
condition 6 11 54.5
subroutine 22 26 84.6
pod 6 6 100.0
total 218 280 77.8


line stmt bran cond sub pod time code
1             # Copyright 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             # math-image --path=FilledRings --all --output=numbers_dash --size=70x30
20             #
21             # offset 0 to 1
22             # same PixelRings fractional offset for midpoint
23             #
24             # default 0.5 |-------int-------| 0.5
25             # |------------------| 0
26             # |------------------| 0.99 or 1.0
27             #
28             # default 0 |-------int-------|--------int-------|
29             # |------------------| 0.5
30             # |------------------| 0.99 or 1.0
31             #
32             # innermost
33             # |-------0-------|-------1--------|
34             #
35             #
36             # offset -0.5 to +0.5
37             # h-0.5 < R < h+0.5
38             # h-0.5 <= R-0.5 < h+0.5
39             # h-0.5 < R+0.5 <= h+0.5
40              
41             # A036702 count Gaussian |z| <= n
42             # A036706 count Gaussian n-1/2 < |z| < n+1/2 with a>0,b>=0, so 1/4
43             # A036707 count Gaussian |z| < n+1/2 with b>=0, so 1/2 plane
44             # A036708 count Gaussian n-1/2 < |z| < n+1/2 with b>=0, so 1/4
45             #
46             # A000328 num points <= circle radius n
47             # 1, 5, 13, 29, 49, 81, 113, 149, 197, 253, 317, 377, 441,
48             # A046109 num points == circle radius n
49             # A051132 num points < circle radius n
50             # 0, 1, 9, 25, 45, 69, 109, 145, 193, 249, 305, 373, 437,
51             # A057655 num points x^2+y^2 <= n
52             # 1, 5, 9, 9, 13, 21, 21, 21, 25, 29, 37, 37, 37, 45, 45,
53             #
54              
55             package Math::PlanePath::FilledRings;
56 1     1   9534 use 5.004;
  1         14  
57 1     1   8 use strict;
  1         1  
  1         36  
58              
59 1     1   6 use vars '$VERSION', '@ISA';
  1         2  
  1         79  
60             $VERSION = 129;
61 1     1   747 use Math::PlanePath;
  1         2  
  1         54  
62             @ISA = ('Math::PlanePath');
63             *_divrem = \&Math::PlanePath::_divrem;
64              
65             use Math::PlanePath::Base::Generic
66 1         42 'is_infinite',
67 1     1   7 'round_nearest';
  1         2  
68              
69 1     1   497 use Math::PlanePath::SacksSpiral;
  1         4  
  1         83  
70             *_rect_to_radius_corners = \&Math::PlanePath::SacksSpiral::_rect_to_radius_corners;
71              
72             # uncomment this to run the ### lines
73             #use Smart::Comments;
74              
75              
76             # N(r) = 1 + 4*sum floor(r^2/(4i+1)) - floor(r^2/(4i+3))
77             #
78             # N(r+1) - N(r)
79             # = 1 + 4*sum floor((r+1)^2/(4i+1)) - floor((r+1)^2/(4i+3))
80             # - 1 + 4*sum floor(r^2/(4i+1)) - floor(r^2/(4i+3))
81             # = 4*sum floor(((r+1)^2-r^2)/(4i+1)) - floor(((r+1)^2-r^2)/(4i+3))
82             # = 4*sum floor((2r+1)/(4i+1)) - floor((2r+1)/(4i+3))
83             #
84             # _cumul[0] index=0 is r=1/2
85             # r = index+1/2
86             # 2r+1 = 2(index+1/2)+1
87             # = 2*index+1+1
88             # = 2*index+2
89             #
90             # 2r+1 >= 4i+1
91             # 2r >= 4i
92             # i <= (2*index+2)/2
93             # i <= index+1
94             #
95             # r=3.5
96             # sqrt(3*3+3*3) = 4.24 out
97             # sqrt(3*3+2*2) = 3.60 out
98             # sqrt(3*3+1*1) = 3.16 in
99             #
100             # * * *
101             # * * * * *
102             # * * * * * * *
103             # * * * o * * * 3+5+7+7+7+5+3 = 37
104             # * * * * * * *
105             # * * * * *
106             # * * *
107             #
108             # N(r) = 1 + 4*( floor(12.25/1)-floor(12.25/3)
109             # + floor(12.25/5)-floor(12.25/7)
110             # + floor(12.25/9)-floor(12.25/11) )
111             # = 37
112             #
113             # (index+1/2)^2 = index^2 + index + 1/4
114             # >= index*(index+1)
115             # (end+1 + 1/2)^2
116             # = (end+3/2)^2
117             # = end^2 + 3*end + 9/4
118             # = end*(end+3) + 2 + 1/4
119             #
120             # (r+1/2)^2 = r^2+r+1/4 floor=r*(r+1)
121             # (r-1/2)^2 = r^2-r+1/4 ceil=r*(r-1)+1
122              
123 1         57 use constant parameter_info_array =>
124             [
125             Math::PlanePath::Base::Generic::parameter_info_nstart1(),
126 1     1   7 ];
  1         2  
127              
128 1     1   6 use constant n_frac_discontinuity => 0;
  1         2  
  1         47  
129 1     1   5 use constant xy_is_visited => 1;
  1         2  
  1         45  
130              
131 1     1   6 use constant dx_minimum => -1;
  1         1  
  1         44  
132 1     1   6 use constant dx_maximum => 1;
  1         2  
  1         41  
133 1     1   6 use constant dy_minimum => -1;
  1         1  
  1         55  
134 1     1   7 use constant dy_maximum => 1;
  1         3  
  1         143  
135             sub x_negative_at_n {
136 0     0 1 0 my ($self) = @_;
137 0         0 return $self->n_start + 4;
138             }
139             sub y_negative_at_n {
140 0     0 1 0 my ($self) = @_;
141 0         0 return $self->n_start + 6;
142             }
143             *_UNDOCUMENTED__dxdy_list = \&Math::PlanePath::_UNDOCUMENTED__dxdy_list_eight;
144 1     1   7 use constant dsumxy_minimum => -2; # diagonals
  1         2  
  1         50  
145 1     1   5 use constant dsumxy_maximum => 2;
  1         2  
  1         54  
146 1     1   7 use constant ddiffxy_minimum => -2;
  1         2  
  1         57  
147 1     1   7 use constant ddiffxy_maximum => 2;
  1         2  
  1         57  
148 1     1   6 use constant dir_maximum_dxdy => (1,-1); # South-East
  1         2  
  1         1152  
149              
150             sub _UNDOCUMENTED__turn_any_right_at_n {
151 0     0   0 my ($self) = @_;
152 0         0 return $self->n_start + 40;
153             }
154              
155             #------------------------------------------------------------------------------
156              
157             sub new {
158 3     3 1 1125 my $self = shift->SUPER::new(@_);
159              
160             # parameters
161 3 50       17 if (! defined $self->{'n_start'}) {
162 3         18 $self->{'n_start'} = $self->default_n_start;
163             }
164              
165             # internals
166 3         11 $self->{'cumul'} = [ 1 ]; # N=0 basis
167 3   50     21 $self->{'offset'} ||= 0;
168              
169 3         9 return $self;
170             }
171              
172             sub _cumul_extend {
173 55     55   183766 my ($self) = @_;
174             ### _cumul_extend() ...
175              
176 55         113 my $cumul = $self->{'cumul'};
177 55         112 my $r2 = ($#$cumul + 3) * $#$cumul + 2;
178 55         92 my $c = 0;
179 55         146 for (my $d = 1; $d <= $r2; $d += 4) {
180 11737         22834 $c += int($r2/$d) - int($r2/($d+2));
181             }
182 55         161 push @$cumul, 4*$c + 1;
183             ### $cumul
184             }
185              
186             sub n_to_xy {
187 27     27 1 1770 my ($self, $n) = @_;
188             ### FilledRings n_to_xy(): $n
189              
190 27         51 $n = $n - $self->{'n_start'}; # to N=0 basis, warning if $n==undef
191 27 100       60 if ($n <= 1) {
192 5 50       13 if ($n < 0) {
193 0         0 return;
194             } else {
195 5         15 return ($n, 0); # 0<=N<=1
196             }
197             }
198 22 50       88 if (is_infinite($n)) {
199 0         0 return ($n,$n);
200             }
201              
202             {
203             # ENHANCE-ME: direction of N+1 from the cumulative lookup
204 22         40 my $int = int($n);
  22         35  
205 22 100       38 if ($n != $int) {
206 6         11 my $frac = $n - $int;
207 6         21 my ($x1,$y1) = $self->n_to_xy($int + $self->{'n_start'});
208 6         19 my ($x2,$y2) = $self->n_to_xy($int+1 + $self->{'n_start'});
209 6 100 66     29 if ($y2 == 0 && $x2 > 0) { $x2 -= 1; }
  3         4  
210 6         13 my $dx = $x2-$x1;
211 6         8 my $dy = $y2-$y1;
212 6         27 return ($frac*$dx + $x1, $frac*$dy + $y1);
213             }
214 16         29 $n = $int;
215             }
216              
217             ### search cumul for: "n=$n"
218 16         27 my $cumul = $self->{'cumul'};
219 16         22 my $r = 1;
220 16         19 for (;;) {
221 32 100       69 if ($r > $#$cumul) {
222 4         7 _cumul_extend ($self);
223             }
224 32 100       64 if ($cumul->[$r] > $n) {
225 16         21 last;
226             }
227 16         24 $r++;
228             }
229             ### $r
230              
231 16         30 $n -= $cumul->[$r-1];
232 16         34 my $len = $cumul->[$r] - $cumul->[$r-1]; # length of this ring
233              
234             ### cumul: "$cumul->[$r-1] to $cumul->[$r]"
235             ### $len
236             ### n rem: $n
237              
238 16         29 $len /= 4; # length of a quadrant of this ring
239 16         41 (my $quadrant, $n) = _divrem ($n, $len);
240              
241             ### len of quadrant: $len
242             ### $quadrant
243             ### n into quadrant: $n
244             ### assert: $quadrant >= 0
245             ### assert: $quadrant < 4
246              
247 16         23 my $rev;
248 16 100       39 if ($rev = ($n > $len/2)) {
249 3         6 $n = $len - $n;
250             }
251             ### $rev
252             ### $n
253              
254             # my $rhi = ($r+1)*$r;
255             # my $rlo = ($r-1)*$r+1;
256              
257 16         49 my $rlo = ($r-0.5+$self->{'offset'})**2;
258 16         28 my $rhi = ($r+0.5+$self->{'offset'})**2;
259 16         26 my $x = $r;
260 16         24 my $y = 0;
261 16         34 while ($n > 0) {
262             ### at: "$x,$y n=$n"
263              
264 10         13 $y++;
265             ### inc y to: $y
266              
267 10 50       25 if ($x*$x + $y*$y > $rhi) {
268 0         0 $x--;
269             ### dec x to: $x
270             ### assert: $x*$x + $y*$y <= $rhi
271             ### assert: $x*$x + $y*$y >= $rlo
272             }
273 10         15 $n--;
274 10 50       23 last if $n <= 0;
275              
276 0 0       0 if (($x-1)*($x-1) + $y*$y >= $rlo) {
277             ### another dec x to: $x
278 0         0 $x--;
279 0         0 $n--;
280 0 0       0 last if $n <= 0;
281             }
282             }
283              
284             # if ($n) {
285             # ### n frac: $n
286             # }
287              
288 16 100       36 if ($rev) {
289 3         9 ($x,$y) = ($y,$x);
290             }
291 16 100       39 if ($quadrant & 2) {
292 4         8 $x = -$x;
293 4         57 $y = -$y;
294             }
295 16 100       30 if ($quadrant & 1) {
296 8         85 ($x,$y) = (-$y, $x);
297             }
298             ### return: "$x, $y"
299 16         42 return ($x, $y);
300             }
301              
302              
303             # h=x^2+y^2
304             # h >= (r-1/2)^2
305             # sqrt(h) >= r-1/2
306             # sqrt(h)+1/2 >= r
307             # r = int (sqrt(h)+1/2)
308             # = int( (2*sqrt(h)+1)/2 }
309             # = int( (sqrt(4*h) + 1)/2 }
310              
311             sub xy_to_n {
312 7     7 1 543 my ($self, $x, $y) = @_;
313             ### FilledRings xy_to_n(): "$x, $y"
314 7         22 $x = round_nearest ($x);
315 7         16 $y = round_nearest ($y);
316              
317 7 100 66     20 if ($x == 0 && $y == 0) {
318 1         3 return $self->{'n_start'};
319             }
320              
321 6         19 my $r = int ((sqrt(4*($x*$x+$y*$y)) + 1) / 2);
322             ### $r
323 6 50       14 if (is_infinite($r)) {
324 0         0 return undef;
325             }
326              
327 6         14 my $cumul = $self->{'cumul'};
328 6         18 while ($#$cumul < $r) {
329 0         0 _cumul_extend ($self);
330             }
331 6         14 my $n = $cumul->[$r-1];
332             ### n base: $n
333              
334 6         10 my $len = $cumul->[$r] - $n;
335             ### $len
336 6         10 $len /= 4;
337             ### len/4: $len
338              
339 6 100       14 if ($y < 0) {
340             ### y neg, rotate 180
341 1         5 $y = -$y;
342 1         3 $x = -$x;
343 1         11 $n += 2*$len;
344             }
345              
346 6 100       14 if ($x < 0) {
347 2         4 $n += $len;
348 2         5 ($x,$y) = ($y,-$x);
349             ### neg x, rotate 90
350             ### n base now: $n
351             }
352              
353             ### assert: $x >= 0
354             ### assert: $y >= 0
355              
356 6         9 my $rev;
357 6 50       14 if ($rev = ($x < $y)) {
358             ### top octant, reverse: "x=$x len/4=".($len/4)." gives ".($len/4 - $x)
359 0         0 ($x,$y) = ($y,$x);
360             }
361              
362 6         8 my $offset = 0;
363 6         11 my $rhi = ($r+1)*$r;
364 6         7 my $rlo = ($r-1)*$r+1;
365             ### assert: $x*$x + $y*$y <= $rhi
366             ### assert: $x*$x + $y*$y >= $rlo
367              
368 6         10 my $tx = $r;
369 6         10 my $ty = 0;
370 6         13 while ($ty < $y) {
371             ### at: "$tx,$ty offset=$offset"
372              
373 3         6 $ty++;
374             ### inc ty to: $ty
375 3 50       9 if ($tx*$tx + $ty*$ty > $rhi) {
376 0         0 $tx--;
377             ### dec tx to: $tx
378             ### assert: $tx*$tx + $ty*$ty <= $rhi
379             ### assert: $tx*$tx + $ty*$ty >= $rlo
380             }
381 3         5 $offset++;
382 3 50 33     13 last if $x == $tx && $y == $ty;
383              
384 0 0       0 if (($tx-1)*($tx-1) + $ty*$ty >= $rlo) {
385             ### another dec tx to: "tx=$tx"
386 0         0 $tx--;
387 0         0 $offset++;
388 0 0       0 last if $y == $ty;
389             }
390             }
391              
392 6         10 $n += $self->{'n_start'};
393 6 50       15 if ($rev) {
394 0         0 return $n + $len - $offset;
395             } else {
396 6         16 return $n + $offset;
397             }
398             }
399              
400             # not exact
401             sub rect_to_n_range {
402 0     0 1   my ($self, $x1,$y1, $x2,$y2) = @_;
403             ### FilledRings rect_to_n_range(): "$x1,$y1 $x2,$y2"
404              
405 0           ($x1,$y1, $x2,$y2) = _rect_to_radius_corners ($x1,$y1, $x2,$y2);
406             ### radius range: "$x1,$y1 $x2,$y2"
407              
408 0 0         if ($x1 >= 1) { $x1 -= 1; }
  0            
409 0 0         if ($y1 >= 1) { $y1 -= 1; }
  0            
410 0           $x2 += 1;
411 0           $y2 += 1;
412              
413             return (int((21*($x1*$x1 + $y1*$y1)) / 7) + $self->{'n_start'},
414 0           int((22*($x2*$x2 + $y2*$y2)) / 7) + $self->{'n_start'} - 1);
415             }
416              
417             1;
418             __END__