1Math::PlanePath::TheodoUrsuesrSpCiornatlr(i3b)uted PerlMDaotchu:m:ePnltaanteiPoanth::TheodorusSpiral(3)
2
6 Math::PlanePath::TheodorusSpiral -- right-angle unit step spiral
7
9 use Math::PlanePath::TheodorusSpiral;
10 my $path = Math::PlanePath::TheodorusSpiral->new;
11 my ($x, $y) = $path->n_to_xy (123);
12
14 This path puts points on the spiral of Theodorus, also called the
15 square root spiral.
16 61 6
18 60
19 27 26 25 24 5
20 28 23 59
21 29 22 58 4
22 30 21 57 3
24 31 20
25 4 56 2
26 32 5 3 19
27 6 2 55 1
28 33 18
29 7 0 1 54 <- Y=0
30 34 17
31 8 53 -1
32 35 16
33 9 52 -2
34 36 15
35 10 14 51 -3
36 37 11 12 13 50
37 -4
38 38 49
39 39 48 -5
40 40 47
41 41 46 -6
42 42 43 44 45
43 ^
46 -6 -5 -4 -3 -2 -1 X=0 1 2 3 4 5 6 7
47 Each step is a unit distance at right angles to the previous radial
49 spoke. So for example,
50 3 <- Y=1+1/sqrt(2)
52 \
53 \
54 ..2 <- Y=1
55 .. |
56 . |
57 0-----1 <- Y=0
58 ^
60 X=0 X=1
61 1 to 2 is a unit step at right angles to the 0 to 1 radial. Then 2 to
63 3 steps at a right angle to radial 0 to 2 which is 45 degrees, etc.
64 The radial distance 0 to 2 is sqrt(2), 0 to 3 is sqrt(3), and in
66 general
67 R = sqrt(N)
69 because each step is a right triangle with radius(N+1)^2 =
71 radius(N)^2 + 1. The resulting shape is very close to an Archimedean
72 spiral with successive loops increasing in radius by pi = 3.14159 or
73 thereabouts each time.
74 X,Y positions returned are fractional and each integer N position is
76 exactly 1 away from the previous. Fractional N values give positions
77 on the straight line between the integer points. (An analytic
78 continuation for a rounded curve between points is possible, but not
79 currently implemented.)
80 Each loop is just under 2*pi^2 = 19.7392 many N points longer than the
82 previous. This means quadratic values 9.8696*k^2 for integer k are an
83 almost straight line. Quadratics close to 9.87 (or a square multiple
84 of that) nearly line up. For example the 22-polygonal numbers have
85 10*k^2 and at low values are nearly straight because 10 is close to
86 9.87, but then spiral away.
87
89 See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
90 classes.
91 The code is currently implemented by adding unit steps in X,Y
93 coordinates, so it's not particularly fast. The last X,Y is saved in
94 the object anticipating successively higher N (not necessarily
95 consecutive), and previous positions 1000 apart are saved for re-use or
96 to go back.
97 "$path = Math::PlanePath::TheodorusSpiral->new ()"
99 Create and return a new Theodorus spiral object.
100 "($x,$y) = $path->n_to_xy ($n)"
102 Return the X,Y coordinates of point number $n on the path.
103 $n can be any value "$n >= 0" and fractions give positions on the
105 spiral in between the integer points.
106 For "$n < 0" the return is an empty list, it being currently
108 considered there are no negative points in the spiral. (The
109 analytic continuation by Davis would be a possibility, though the
110 resulting "inner spiral" makes positive and negative points overlap
111 a bit. A spiral starting at X=-1 would fit in between the positive
112 points.)
113 "$rsquared = $path->n_to_rsquared ($n)"
115 Return the radial distance R^2 of point $n, or "undef" if $n is
116 negative. For integer $n this is simply $n itself.
117 "$n = $path->xy_to_n ($x,$y)"
119 Return an integer point number for coordinates "$x,$y". Each
120 integer N is considered the centre of a circle of diameter 1 and an
121 "$x,$y" within that circle returns N.
122 The unit steps of the spiral means those unit circles don't
124 overlap, but the loops are roughly 3.14 apart so there's gaps in
125 between. If "$x,$y" is not within one of the unit circles then the
126 return is "undef".
127 "$str = $path->figure ()"
129 Return string "circle".
130
132 N to RSquared
133 For integer N the spiral has radius R=sqrt(N) and the square is simply
134 RSquared=R^2=N. For fractional N, the point is on a straight line at
135 right angles to the integer position, so
136 R = hypot(sqrt(Ninteger), Nfrac)
138 RSquared = (sqrt(Ninteger))^2 + Nfrac^2
139 = Ninteger + Nfrac^2
140 X,Y to N
142 For a given X,Y the radius R=hypot(X,Y) determines the N position as
143 N=R^2. An N point up to 0.5 away radially might cover X,Y, so the
144 range of N to consider is
145 Nlo = (R-.5)^2
147 Nhi = (R+.5)^2
148 A simple search is made through those N's seeking which, if any, covers
150 X,Y. The number of N's searched is Nhi-Nlo = 2*R+1 which is about 1/3
151 of a loop around the spiral (2*R/2*pi*R ~= 1/3). Actually 0.51 is used
152 so as to guard against floating point round-off, which is then about
153 4*.51 = 2.04*R many points.
154 The angle of the X,Y position determines which part of the spiral is
156 intersected, but using that doesn't seem particularly easy. The angle
157 for a given N is an arctan sum and don't currently have a good closed-
158 form or converging series to invert, or apply some Newton's method, or
159 whatever.
160 Rectangle to N Range
162 For rect_to_n_range() the corner furthest from the origin determines
163 the high N. For that corner
164 Rhi = hypot(xhi,yhi)
166 Nhi = (Rhi+.5)^2
167 The extra .5 is since a unit circle figure centred as much as .5
169 further out might intersect the xhi,yhi. The square root hypot() can
170 be avoided by the following over-estimate, and ceil can keep it in
171 integers for integer Nhi.
172 Nhi = Rhi^2 + Rhi + 1/4
174 <= Xhi^2+Yhi^2 + Xhi+Yhi + 1 # since Rhi<=Xhi+Yhi
175 = Xhi*(Xhi+1) + Yhi*(Yhi+1) + 1
176 <= ceilXhi*(ceilXhi+1) + ceilYhi*(ceilYhi+1) + 1
177 With either formula the worst case is when Nhi doesn't intersect the
179 xhi,yhi corner but is just before it, anti-clockwise. Nhi is then a
180 full revolution bigger than it need be, depending where the other
181 corners fall.
182 Similarly for the corner or axis crossing nearest the origin (when the
184 origin itself isn't covered by the rectangle),
185 Rlo = hypot(Xlo,Ylo)
187 Nlo = (Rlo-.5)^2, or 0 if origin covered by rectangle
188 And again in integers without a square root if desired,
190 Nlo = Rlo^2 - Rlo + 1/4
192 >= Xlo^2+Ylo^2 - (Xlo+Ylo) # since Xlo+Ylo>=Rlo
193 = Xlo*(Xlo-1) + Ylo*(Ylo-1)
194 >= floorXlo*(floorXlo-1) + floorYlo(floorYlo-1)
195 The worst case is when this Nlo doesn't intersect the xlo,ylo corner
197 but is just after it anti-clockwise, so Nlo is a full revolution
198 smaller than it need be.
199
201 Entries in Sloane's Online Encyclopedia of Integer Sequences related to
202 this path include
203 <http://oeis.org/A072895> (etc)
205 A072895 N just below X axis
207 A137515 N-1 just below X axis
208 counting num points for n revolutions
209 A172164 loop length increases
210 A164102 2*pi^2
211
213 Math::PlanePath, Math::PlanePath::ArchimedeanChords,
214 Math::PlanePath::SacksSpiral, Math::PlanePath::MultipleRings
215
217 <http://user42.tuxfamily.org/math-planepath/index.html>
218
220 Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019,
221 2020 Kevin Ryde
222 This file is part of Math-PlanePath.
224 Math-PlanePath is free software; you can redistribute it and/or modify
226 it under the terms of the GNU General Public License as published by
227 the Free Software Foundation; either version 3, or (at your option) any
228 later version.
229 Math-PlanePath is distributed in the hope that it will be useful, but
231 WITHOUT ANY WARRANTY; without even the implied warranty of
232 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
233 General Public License for more details.
234 You should have received a copy of the GNU General Public License along
236 with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
237perl v5.36.0 2023-01-20Math::PlanePath::TheodorusSpiral(3)