File Coverage

blib/lib/Math/PlanePath/SquareSpiral.pm
Criterion Covered Total %
statement 88 114 77.1
branch 27 36 75.0
condition 2 2 100.0
subroutine 13 19 68.4
pod 7 7 100.0
total 137 178 76.9


line stmt bran cond sub pod time code
1             # Copyright 2010, 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             # http://d4maths.lowtech.org/mirage/ulam.htm
20             # http://d4maths.lowtech.org/mirage/img/ulam.gif
21             # sample gif of primes made by APL or something
22             #
23             # http://www.sciencenews.org/view/generic/id/2696/title/Prime_Spirals
24             # Ulam's spiral of primes
25             #
26             # http://web.archive.org/web/20160424015805id_/http://yoyo.cc.monash.edu.au/%7Ebunyip/primes/primeSpiral.htm
27             # http://yoyo.cc.monash.edu.au/%7Ebunyip/primes/triangleUlam.htm
28             # Pulchritudinous Primes of Ulam spiral.
29              
30             # http://mathworld.wolfram.com/PrimeSpiral.html
31             #
32             # Mark C. Chu-Carroll "The Surprises Never Eend: The Ulam Spiral of Primes"
33             # http://scienceblogs.com/goodmath/2010/06/the_surprises_never_eend_the_u.php
34             #
35             # http://yoyo.cc.monash.edu.au/%7Ebunyip/primes/index.html
36             # including image highlighting the lines
37              
38             # S. M. Ellerstein, The Square Spiral, Journal of Recreational
39             # Mathematics 29 (#3, 1998) 188; 30 (#4, 1999-2000), 246-250.
40             #
41             # M. Stein and S. M. Ulam. "An Observation on the
42             # Distribution of Primes." Amer. Math. Monthly 74, 43-44,
43             # 1967.
44             #
45             # M.L. Stein; S. M. Ulam; and M. B. Wells. "A Visual
46             # Display of Some Properties of the Distribution of Primes."
47             # Amer. Math. Monthly 71, 516-520, 1964.
48              
49             # cf sides alternately prime and Fibonacci
50             # A160790 corner N
51             # A160791 side lengths, alternately integer and triangular adding that integer
52             # A160792 corner N
53             # A160793 side lengths, alternately integer and sum primes
54             # A160794 corner N
55             # A160795 side lengths, alternately primes and fibonaccis
56              
57              
58             package Math::PlanePath::SquareSpiral;
59 8     8   2077 use 5.004;
  8         24  
60 8     8   39 use strict;
  8         13  
  8         288  
61             #use List::Util 'max';
62             *max = \&Math::PlanePath::_max;
63              
64 8     8   59 use vars '$VERSION', '@ISA';
  8         15  
  8         493  
65             $VERSION = 128;
66 8     8   1080 use Math::PlanePath;
  8         14  
  8         260  
67             *_sqrtint = \&Math::PlanePath::_sqrtint;
68 8     8   3527 use Math::PlanePath::Base::NSEW;
  8         19  
  8         288  
69             @ISA = ('Math::PlanePath::Base::NSEW',
70             'Math::PlanePath');
71              
72             use Math::PlanePath::Base::Generic
73 8     8   43 'round_nearest';
  8         15  
  8         440  
74              
75             # uncomment this to run the ### lines
76             #use Smart::Comments '###';
77              
78              
79             # Note: this shared by other paths
80 8         391 use constant parameter_info_array =>
81             [
82             { name => 'wider',
83             display => 'Wider',
84             type => 'integer',
85             minimum => 0,
86             default => 0,
87             width => 3,
88             description => 'Wider path.',
89             },
90             Math::PlanePath::Base::Generic::parameter_info_nstart1(),
91 8     8   38 ];
  8         14  
92              
93 8     8   40 use constant xy_is_visited => 1;
  8         15  
  8         1284  
94              
95             # 2w+4 -- 2w+3 ----- w+2
96             # | |
97             # 2w+5 0------- w+1
98             # |
99             # 2w+6 ---
100             # ^
101             # X=0
102             #
103             sub x_negative_at_n {
104 0     0 1 0 my ($self) = @_;
105 0 0       0 return $self->n_start + ($self->{'wider'} ? 0 : 4);
106             }
107             sub y_negative_at_n {
108 0     0 1 0 my ($self) = @_;
109 0         0 return $self->n_start + 2*$self->{'wider'} + 6;
110             }
111             sub _UNDOCUMENTED__dxdy_list_at_n {
112 0     0   0 my ($self) = @_;
113 0         0 return $self->n_start + 2*$self->{'wider'} + 4;
114             }
115              
116 8     8   59 use constant turn_any_right => 0; # only left or straight
  8         14  
  8         7127  
117              
118             sub _UNDOCUMENTED__turn_any_left_at_n {
119 0     0   0 my ($self) = @_;
120 0         0 return $self->n_start + $self->{'wider'} + 1;
121             }
122              
123              
124             #------------------------------------------------------------------------------
125              
126             sub new {
127 18     18 1 2177 my $self = shift->SUPER::new (@_);
128              
129             # parameters
130 18   100     81 $self->{'wider'} ||= 0; # default
131 18 100       39 if (! defined $self->{'n_start'}) {
132 16         49 $self->{'n_start'} = $self->default_n_start;
133             }
134              
135 18         38 return $self;
136             }
137              
138             # wider==0
139             # base from bottom-right corner
140             # d = [ 1, 2, 3, 4 ]
141             # N = [ 2, 10, 26, 50 ]
142             # N = (4 d^2 - 4 d + 2)
143             # d = 1/2 + sqrt(1/4 * $n + -4/16)
144             #
145             # wider==1
146             # base from bottom-right corner
147             # d = [ 1, 2, 3, 4 ]
148             # N = [ 3, 13, 31, 57 ]
149             # N = (4 d^2 - 2 d + 1)
150             # d = 1/4 + sqrt(1/4 * $n + -3/16)
151             #
152             # wider==2
153             # base from bottom-right corner
154             # d = [ 1, 2, 3, 4 ]
155             # N = [ 4, 16, 36, 64 ]
156             # N = (4 d^2)
157             # d = 0 + sqrt(1/4 * $n + 0)
158             #
159             # wider==3
160             # base from bottom-right corner
161             # d = [ 1, 2, 3 ]
162             # N = [ 5, 19, 41 ]
163             # N = (4 d^2 + 2 d - 1)
164             # d = -1/4 + sqrt(1/4 * $n + 5/16)
165             #
166             # N = 4*d^2 + (-4+2*w)*d + (2-w)
167             # = 4*$d*$d + (-4+2*$w)*$d + (2-$w)
168             # d = 1/2-w/4 + sqrt(1/4*$n + b^2-4ac)
169             # (b^2-4ac)/(2a)^2 = [ (2w-4)^2 - 4*4*(2-w) ] / 64
170             # = [ 4w^2 - 16w + 16 - 32 + 16w ] / 64
171             # = [ 4w^2 - 16 ] / 64
172             # = [ w^2 - 4 ] / 16
173             # d = 1/2-w/4 + sqrt(1/4*$n + (w^2 - 4) / 16)
174             # = 1/4 * (2-w + sqrt(4*$n + w^2 - 4))
175             # = 0.25 * (2-$w + sqrt(4*$n + $w*$w - 4))
176             #
177             # then offset the base by +4*$d+$w-1 for top left corner for +/- remainder
178             # rem = $n - (4*$d*$d + (-4+2*$w)*$d + (2-$w) + 4*$d + $w - 1)
179             # = $n - (4*$d*$d + (-4+2*$w)*$d + 2 - $w + 4*$d + $w - 1)
180             # = $n - (4*$d*$d + (-4+2*$w)*$d + 1 - $w + 4*$d + $w)
181             # = $n - (4*$d*$d + (-4+2*$w)*$d + 1 + 4*$d)
182             # = $n - (4*$d*$d + (2*$w)*$d + 1)
183             # = $n - ((4*$d + 2*$w)*$d + 1)
184             #
185              
186             sub n_to_xy {
187 372     372 1 24439 my ($self, $n) = @_;
188             #### SquareSpiral n_to_xy: $n
189              
190 372         737 $n = $n - $self->{'n_start'}; # starting $n==0, warn if $n==undef
191 372 50       752 if ($n < 0) {
192             #### before n_start ...
193 0         0 return;
194             }
195              
196 372         556 my $w = $self->{'wider'};
197 372         626 my $w_right = int($w/2);
198 372         489 my $w_left = $w - $w_right;
199 372 100       660 if ($n <= $w+1) {
200             #### centre horizontal
201             # n=0 at w_left
202             # x = $n - int(($w+1)/2)
203             # = $n - int(($w+1)/2)
204 43         109 return ($n - $w_left, # n=0 at w_left
205             0);
206             }
207              
208 329         800 my $d = int ((2-$w + _sqrtint(4*$n + $w*$w)) / 4);
209             #### d frac: ((2-$w + sqrt(int(4*$n) + $w*$w)) / 4)
210             #### $d
211              
212             #### base: 4*$d*$d + (-4+2*$w)*$d + (2-$w)
213 329         570 $n -= ((4*$d + 2*$w)*$d);
214             #### remainder: $n
215              
216 329 100       614 if ($n >= 0) {
217 144 100       239 if ($n <= 2*$d) {
218             ### left vertical
219 75         197 return (-$d - $w_left,
220             -$n + $d);
221             } else {
222             ### bottom horizontal
223 69         187 return ($n - $w_left - 3*$d,
224             -$d);
225             }
226             } else {
227 185 100       665 if ($n >= -2*$d-$w) {
228             ### top horizontal
229 109         313 return (-$n - $d - $w_left,
230             $d);
231             } else {
232             ### right vertical
233 76         186 return ($d + $w_right,
234             $n + 3*$d + $w);
235             }
236             }
237             }
238              
239             sub xy_to_n {
240 27     27 1 1557 my ($self, $x, $y) = @_;
241              
242 27         41 my $w = $self->{'wider'};
243 27         57 my $w_right = int($w/2);
244 27         41 my $w_left = $w - $w_right;
245 27         69 $x = round_nearest ($x);
246 27         58 $y = round_nearest ($y);
247             ### xy_to_n: "x=$x, y=$y"
248             ### $w_left
249             ### $w_right
250              
251 27         39 my $d;
252 27 100       60 if (($d = $x - $w_right) > abs($y)) {
253             ### right vertical
254             ### $d
255             #
256             # base bottom right per above
257             ### BR: 4*$d*$d + (-4+2*$w)*$d + (2-$w)
258             # then +$d-1 for the y=0 point
259             # N_Y0 = 4*$d*$d + (-4+2*$w)*$d + (2-$w) + $d-1
260             # = 4*$d*$d + (-3+2*$w)*$d + (2-$w) + -1
261             # = 4*$d*$d + (-3+2*$w)*$d + 1-$w
262             ### N_Y0: (4*$d + -3 + 2*$w)*$d + 1-$w
263             #
264 6         18 return (4*$d + -3 + 2*$w)*$d - $w + $y + $self->{'n_start'};
265             }
266              
267 21 100       40 if (($d = -$x - $w_left) > abs($y)) {
268             ### left vertical
269             ### $d
270             #
271             # top left per above
272             ### TL: 4*$d*$d + (2*$w)*$d + 1
273             # then +$d for the y=0 point
274             # N_Y0 = 4*$d*$d + (2*$w)*$d + 1 + $d
275             # = 4*$d*$d + (1 + 2*$w)*$d + 1
276             ### N_Y0: (4*$d + 1 + 2*$w)*$d + 1
277             #
278 4         14 return (4*$d + 1 + 2*$w)*$d - $y + $self->{'n_start'};
279             }
280              
281 17         23 $d = abs($y);
282 17 100       30 if ($y > 0) {
283             ### top horizontal
284             ### $d
285             #
286             # top left per above
287             ### TL: 4*$d*$d + (2*$w)*$d + 1
288             # then -($d+$w_left) for the x=0 point
289             # N_X0 = 4*$d*$d + (2*$w)*$d + 1 + -($d+$w_left)
290             # = 4*$d*$d + (-1 + 2*$w)*$d + 1 - $w_left
291             ### N_Y0: (4*$d - 1 + 2*$w)*$d + 1 - $w_left
292             #
293 8         21 return (4*$d - 1 + 2*$w)*$d - $w_left - $x + $self->{'n_start'};
294             }
295              
296             ### bottom horizontal, and centre y=0
297             ### $d
298             #
299             # top left per above
300             ### TL: 4*$d*$d + (2*$w)*$d + 1
301             # then +2*$d to bottom left, +$d+$w_left for the x=0 point
302             # N_X0 = 4*$d*$d + (2*$w)*$d + 1 + 2*$d + $d+$w_left)
303             # = 4*$d*$d + (3 + 2*$w)*$d + 1 + $w_left
304             ### N_Y0: (4*$d + 3 + 2*$w)*$d + 1 + $w_left
305             #
306 9         38 return (4*$d + 3 + 2*$w)*$d + $w_left + $x + $self->{'n_start'};
307             }
308              
309             # hi is exact but lo is not
310             # not exact
311             sub rect_to_n_range {
312 0     0 1 0 my ($self, $x1,$y1, $x2,$y2) = @_;
313              
314 0         0 $x1 = round_nearest ($x1);
315 0         0 $y1 = round_nearest ($y1);
316 0         0 $x2 = round_nearest ($x2);
317 0         0 $y2 = round_nearest ($y2);
318              
319             # ENHANCE-ME: find actual minimum if rect doesn't cover 0,0
320 0         0 return ($self->{'n_start'},
321             max ($self->xy_to_n($x1,$y1),
322             $self->xy_to_n($x2,$y1),
323             $self->xy_to_n($x1,$y2),
324             $self->xy_to_n($x2,$y2)));
325              
326             # my $w = $self->{'wider'};
327             # my $w_right = int($w/2);
328             # my $w_left = $w - $w_right;
329             #
330             # my $d = 1 + max (abs($y1),
331             # abs($y2),
332             # $x1 - $w_right, -$x1 - $w_left,
333             # $x2 - $w_right, -$x2 - $w_left,
334             # 1);
335             # ### $d
336             # ### is: $d*$d
337             #
338             # # ENHANCE-ME: find actual minimum if rect doesn't cover 0,0
339             # return (1,
340             # (4*$d - 4 + 2*$w)*$d + 2); # bottom-right
341             }
342              
343              
344             # [ 1, 2, 3, 4, 5 ],
345             # [ 1, 3, 7, 13, 21 ]
346             # N = (d^2 - d + 1)
347             # = ($d**2 - $d + 1)
348             # = (($d - 1)*$d + 1)
349             # d = 1/2 + sqrt(1 * $n + -3/4)
350             # = (1 + sqrt(4*$n - 3)) / 2
351             #
352             # wider=3
353             # [ 2, 3, 4, 5 ],
354             # [ 6, 13, 22, 33 ]
355             # N = (d^2 + 2 d - 2)
356             # = ($d**2 + 2*$d - 2)
357             # = (($d + 2)*$d - 2)
358             # d = -1 + sqrt(1 * $n + 3)
359             #
360             # wider=5
361             # [ 2, 3, 4, 5 ],
362             # [ 8, 17, 28, 41 ]
363             # N = (d^2 + 4 d - 4)
364             # = ($d**2 + 4*$d - 4)
365             # = (($d + 4)*$d - 4)
366             # d = -2 + sqrt(1 * $n + 8)
367             #
368             # wider=7
369             # [ 2, 3, 4, 5 ],
370             # [ 10, 21, 34, 49 ]
371             # N = (d^2 + 6 d - 6)
372             # = ($d**2 + 6*$d - 6)
373             # = (($d + 6)*$d - 6)
374             # d = -3 + sqrt(1 * $n + 15)
375             #
376             #
377             # N = (d^2 + (w-1)*d + 1-w)
378             # d = (1-w)/2 + sqrt($n + (w^2 + 2w - 3)/4)
379             # = (1-w + sqrt(4*$n + (w-3)(w+1))) / 2
380             #
381             # extra subtract d+w-1
382             # Nbase = (d^2 + (w-1)*d + 1-w) + d+w-1
383             # = d^2 + w*d
384              
385             sub n_to_dxdy {
386 118     118 1 558 my ($self, $n) = @_;
387             ### n_to_dxdy(): $n
388              
389 118         180 $n = $n - $self->{'n_start'}; # starting $n==0, warn if $n==undef
390 118 100       215 if ($n < 0) {
391             #### before n_start ...
392 2         17 return;
393             }
394              
395 116         169 my $w = $self->{'wider'};
396 116         277 my $d = int((1-$w + _sqrtint(4*$n + ($w+2)*$w+1)) / 2);
397              
398 116         176 my $int = int($n);
399 116         162 $n -= $int; # fraction 0 <= $n < 1
400 116         177 $int -= ($d+$w)*$d-1;
401              
402             ### $d
403             ### $w
404             ### $n
405             ### $int
406              
407 116         160 my ($dx, $dy);
408 116 100       206 if ($int <= 0) {
409 78 100       122 if ($int < 0) {
410             ### horizontal ...
411 73         96 $dx = 1;
412 73         95 $dy = 0;
413             } else {
414             ### corner horiz to vert ...
415 5         7 $dx = 1-$n;
416 5         10 $dy = $n;
417             }
418             } else {
419 38 100       70 if ($int < $d) {
420             ### vertical ...
421 34         49 $dx = 0;
422 34         51 $dy = 1;
423             } else {
424             ### corner vert to horiz ...
425 4         6 $dx = -$n;
426 4         10 $dy = 1-$n;
427             }
428             }
429              
430 116 100       224 unless ($d % 2) {
431             ### rotate +180 for even d ...
432 50         74 $dx = -$dx;
433 50         66 $dy = -$dy;
434             }
435              
436             ### result: "$dx, $dy"
437 116         256 return ($dx,$dy);
438             }
439              
440              
441              
442             # old bit:
443             #
444             # wider==0
445             # base from two-way diagonal top-right and bottom-left
446             # s even for top-right diagonal doing top leftwards then left downwards
447             # s odd for bottom-left diagonal doing bottom rightwards then right pupwards
448             # s = [ 0, 1, 2, 3, 4, 5, 6 ]
449             # N = [ 1, 1, 3, 7, 13, 21, 31 ]
450             # +0 +2 +4 +6 +8 +10
451             # 2 2 2 2 2
452             #
453             # n = (($d - 1)*$d + 1)
454             # s = 1/2 + sqrt(1 * $n + -3/4)
455             # = .5 + sqrt ($n - .75)
456             #
457             #
458              
459             #------------------------------------------------------------------------------
460              
461             sub _NOTDOCUMENTED_n_to_figure_boundary {
462 0     0     my ($self, $n) = @_;
463             ### _NOTDOCUMENTED_n_to_figure_boundary(): $n
464              
465             # adjust to N=1 at origin X=0,Y=0
466 0           $n = $n - $self->{'n_start'} + 1;
467              
468 0 0         if ($n < 1) {
469 0           return undef;
470             }
471              
472 0           my $wider = $self->{'wider'};
473 0 0         if ($n <= $wider) {
474             # single block row
475             # +---+-----+----+
476             # | 1 | ... | $n | boundary = 2*N + 2
477             # +---+-----+----+
478 0           return 2*$n + 2;
479             }
480              
481 0           my $d = int((_sqrtint(4*$n + $wider*$wider - 2) - $wider) / 2);
482             ### $d
483             ### $wider
484             ### cmp: $d*($d+1+$wider) + $wider + 1
485              
486 0 0         if ($n > $d*($d+1+$wider)) {
487 0           $wider++;
488             ### increment for +2 after turn ...
489             }
490 0           return 4*$d + 2*$wider + 2;
491             }
492              
493             #------------------------------------------------------------------------------
494             1;
495             __END__