File Coverage

blib/lib/Math/PlanePath/SacksSpiral.pm
Criterion Covered Total %
statement 94 103 91.2
branch 19 30 63.3
condition 0 3 0.0
subroutine 28 28 100.0
pod 4 4 100.0
total 145 168 86.3


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             # could loop by more or less, eg. 4*n^2 each time like a square spiral
20             # (Kevin Vicklund at the_surprises_never_eend_the_u.php)
21              
22              
23             package Math::PlanePath::SacksSpiral;
24 15     15   1635 use 5.004;
  15         61  
25 15     15   120 use strict;
  15         32  
  15         403  
26 15     15   2861 use Math::Libm 'hypot';
  15         28604  
  15         942  
27 15     15   4144 use POSIX 'floor';
  15         45962  
  15         147  
28             #use List::Util 'max';
29             *max = \&Math::PlanePath::_max;
30              
31 15     15   12795 use Math::PlanePath;
  15         80  
  15         444  
32 15     15   8164 use Math::PlanePath::MultipleRings;
  15         39  
  15         535  
33              
34 15     15   109 use vars '$VERSION', '@ISA';
  15         32  
  15         1168  
35             $VERSION = 129;
36             @ISA = ('Math::PlanePath');
37              
38              
39             # uncomment this to run the ### lines
40             #use Smart::Comments;
41              
42              
43 15     15   100 use constant n_start => 0;
  15         32  
  15         796  
44 15     15   88 use constant figure => 'circle';
  15         39  
  15         707  
45 15     15   88 use constant x_negative_at_n => 2;
  15         28  
  15         751  
46 15     15   94 use constant y_negative_at_n => 3;
  15         30  
  15         687  
47              
48 15     15   86 use constant 1.02; # for leading underscore
  15         235  
  15         581  
49 15     15   87 use constant _TWO_PI => 4*atan2(1,0);
  15         49  
  15         1383  
50              
51             # at N=k^2 polygon of 2k+1 sides R=k
52             # dX -> sin(2pi/(2k+1))*k
53             # -> 2pi/(2k+1) * k
54             # -> pi
55 15     15   106 use constant dx_minimum => - 2*atan2(1,0); # -pi
  15         30  
  15         922  
56 15     15   99 use constant dx_maximum => 2*atan2(1,0); # +pi
  15         29  
  15         1009  
57 15     15   95 use constant dy_minimum => - 2*atan2(1,0);
  15         35  
  15         854  
58 15     15   95 use constant dy_maximum => 2*atan2(1,0);
  15         29  
  15         742  
59              
60 15     15   95 use constant turn_any_right => 0; # left always
  15         46  
  15         873  
61 15     15   103 use constant turn_any_straight => 0; # left always
  15         29  
  15         1433  
62              
63              
64             #------------------------------------------------------------------------------
65             # sub _as_float {
66             # my ($x) = @_;
67             # if (ref $x) {
68             # if ($x->isa('Math::BigInt')) {
69             # return Math::BigFloat->new($x);
70             # }
71             # if ($x->isa('Math::BigRat')) {
72             # return $x->as_float;
73             # }
74             # }
75             # return $x;
76             # }
77              
78             # Note: this is "use Math::BigFloat" not "require Math::BigFloat" because
79             # BigFloat 1.997 does some setups in its import() needed to tie-in to the
80             # BigInt back-end, or something.
81             use constant::defer _bigfloat => sub {
82 1 50   1   84 eval "use Math::BigFloat; 1" or die $@;
  1         8  
  1         2  
  1         5  
83 1         5 return "Math::BigFloat";
84 15     15   8244 };
  15         12418  
  15         158  
85              
86             sub n_to_xy {
87 55     55 1 4004 my ($self, $n) = @_;
88 55 50       130 if ($n < 0) {
89 0         0 return;
90             }
91 55         87 my $two_pi = _TWO_PI();
92              
93 55 50       102 if (ref $n) {
94 0 0       0 if ($n->isa('Math::BigInt')) {
95 0         0 $n = _bigfloat()->new($n);
96             }
97 0 0       0 if ($n->isa('Math::BigRat')) {
98 0         0 $n = $n->as_float;
99             }
100 0 0       0 if ($n->isa('Math::BigFloat')) {
101 0   0     0 $two_pi = 2 * Math::BigFloat->bpi ($n->accuracy
102             || $n->precision
103             || $n->div_scale);
104             }
105             }
106              
107 55         122 my $r = sqrt($n);
108 55         119 my $theta = $two_pi * ($r - int($r)); # 0 <= $theta < 2*pi
109 55         201 return ($r * cos($theta),
110             $r * sin($theta));
111              
112             }
113              
114             sub n_to_rsquared {
115 3     3 1 9 my ($self, $n) = @_;
116 3 50       12 if ($n < 0) { return undef; }
  0         0  
117 3         9 return $n; # exactly RSquared=$n
118             }
119              
120             sub xy_to_n {
121 5     5 1 311 my ($self, $x, $y) = @_;
122             ### SacksSpiral xy_to_n(): "$x, $y"
123              
124 5         15 my $theta_frac = Math::PlanePath::MultipleRings::_xy_to_angle_frac($x,$y);
125             ### assert: 0 <= $theta_frac && $theta_frac < 1
126              
127             # the nearest arc, integer
128 5         27 my $s = floor (hypot($x,$y) - $theta_frac + 0.5);
129              
130             # the nearest N on the arc
131 5         16 my $n = floor ($s*$s + $theta_frac * (2*$s + 1) + 0.5);
132              
133             # check within 0.5 radius
134 5         11 my ($nx, $ny) = $self->n_to_xy($n);
135              
136             ### $theta_frac
137             ### raw hypot: hypot($x,$y)
138             ### $s
139             ### $n
140             ### hypot: hypot($nx-$x, $ny-$y)
141 5 50       19 if (hypot($nx-$x,$ny-$y) <= 0.5) {
142 5         16 return $n;
143             } else {
144 0         0 return undef;
145             }
146             }
147              
148             # r^2 = x^2 + y^2
149             # (r+1)^2 = r^2 + 2r + 1
150             # r < x+y
151             # (r+1)^2 < x^2+y^2 + x + y + 1
152             # < (x+.5)^2 + (y+.5)^2 + 1
153             # (x+1)^2 + (y+1)^2 = x^2+y^2 + 2x+2y+2
154             #
155             # (x+1)^2 + (y+1)^2 - (r+1)^2
156             # = x^2+y^2 + 2x+2y+2 - (r^2 + 2r + 1)
157             # = x^2+y^2 + 2x+2y+2 - x^2-y^2 - 2*sqrt(x^2+y^2) - 1
158             # = 2x+2y+1 - 2*sqrt(x^2+y^2)
159             # >= 2x+2y+1 - 2*(x+y)
160             # = 1
161             #
162             # (x+e)^2 + (y+e)^2 - (r+e)^2
163             # = x^2+y^2 + 2xe+2ye + 2e^2 - (r^2 + 2re + e^2)
164             # = x^2+y^2 + 2xe+2ye + 2e^2 - x^2-y^2 - 2*e*sqrt(x^2+y^2) - e^2
165             # = 2xe+2ye + e^2 - 2*e*sqrt(x^2+y^2)
166             # >= 2xe+2ye + e^2 - 2*e*(x+y)
167             # = e^2
168             #
169             # x+1,y+1 increases the radius by at least 1 thus pushing it to the outside
170             # of a ring. Actually it's more, as much as sqrt(2)=1.4142 on the leading
171             # diagonal X=Y. But the over-estimate is close enough for now.
172             #
173              
174             # r = hypot(xmin,ymin)
175             # Nlo = (r-1/2)^2
176             # = r^2 - r + 1/4
177             # >= x^2+y^2 - (x+y) because x+y >= r
178             # = x(x-1) + y(y-1)
179             # >= floorx(floorx-1) + floory(floory-1)
180             # in integers if round down to x=0 then x*(x-1)=0 too, so not negative
181             #
182             # r = hypot(xmax,ymax)
183             # Nhi = (r+1/2)^2
184             # = r^2 + r + 1/4
185             # <= x^2+y^2 + (x+y) + 1
186             # = x(x+1) + y(y+1) + 1
187             # <= ceilx(ceilx+1) + ceily(ceily+1) + 1
188              
189             # Note: this code shared by TheodorusSpiral. If start using the polar angle
190             # for more accuracy here then unshare it first.
191             #
192             # not exact
193             sub rect_to_n_range {
194 52     52 1 1106 my ($self, $x1,$y1, $x2,$y2) = @_;
195 52         107 ($x1,$y1, $x2,$y2) = _rect_to_radius_corners ($x1,$y1, $x2,$y2);
196              
197             ### $x_min
198             ### $y_min
199             ### N min: $x_min*($x_min-1) + $y_min*($y_min-1)
200              
201             ### $x_max
202             ### $y_max
203             ### N max: $x_max*($x_max+1) + $y_max*($y_max+1) + 1
204              
205 52         162 return ($x1*($x1-1) + $y1*($y1-1),
206             $x2*($x2+1) + $y2*($y2+1) + 1);
207             }
208              
209             #------------------------------------------------------------------------------
210             # generic
211              
212             # $x1,$y1, $x2,$y2 is a rectangle.
213             # Return ($xmin,$ymin, $xmax,$ymax).
214             #
215             # The two points are respectively minimum and maximum radius from the
216             # origin, rounded down or up to integers.
217             #
218             # If the rectangle is entirely one quadrant then the points are two opposing
219             # corners. But if an axis is crossed then the minimum is on that axis and
220             # if the origin is covered then the minimum is 0,0.
221             #
222             # Currently the return is abs() absolute values of the places. Could change
223             # that if there was any significance to the quadrant containing the min/max
224             # corners.
225             #
226             sub _rect_to_radius_corners {
227 919     919   1913 my ($x1,$y1, $x2,$y2) = @_;
228              
229 919 100       2368 ($x1,$x2) = ($x2,$x1) if $x1 > $x2;
230 919 100       2097 ($y1,$y2) = ($y2,$y1) if $y1 > $y2;
231              
232 919 100       4247 return (int($x2 < 0 ? -$x2
    100          
    100          
    100          
233             : $x1 > 0 ? $x1
234             : 0),
235             int($y2 < 0 ? -$y2
236             : $y1 > 0 ? $y1
237             : 0),
238              
239             max(_ceil(abs($x1)), _ceil(abs($x2))),
240             max(_ceil(abs($y1)), _ceil(abs($y2))));
241             }
242              
243             sub _ceil {
244 3676     3676   6031 my ($x) = @_;
245 3676         5221 my $int = int($x);
246 3676 100       9638 return ($x > $int ? $int+1 : $int);
247             }
248              
249             # FIXME: prefer to stay in integers if possible
250             # return ($rlo,$rhi) which is the radial distance range found in the rectangle
251             sub _rect_to_radius_range {
252 867     867   6248 my ($x1,$y1, $x2,$y2) = @_;
253              
254 867         2231 ($x1,$y1, $x2,$y2) = _rect_to_radius_corners ($x1,$y1, $x2,$y2);
255 867         4584 return (hypot($x1,$y1),
256             hypot($x2,$y2));
257             }
258              
259             1;
260             __END__