1Math::PlanePath::SquareUSspeirraClo(n3t)ributed Perl DocMuamtehn:t:aPtliaonnePath::SquareSpiral(3)
2
3
4
6 Math::PlanePath::SquareSpiral -- integer points drawn around a square
7 (or rectangle)
8
10 use Math::PlanePath::SquareSpiral;
11 my $path = Math::PlanePath::SquareSpiral->new;
12 my ($x, $y) = $path->n_to_xy (123);
13
15 This path makes a square spiral,
16
17 37--36--35--34--33--32--31 3
18 | |
19 38 17--16--15--14--13 30 2
20 | | | |
21 39 18 5---4---3 12 29 1
22 | | | | | |
23 40 19 6 1---2 11 28 ... <- Y=0
24 | | | | | |
25 41 20 7---8---9--10 27 52 -1
26 | | | |
27 42 21--22--23--24--25--26 51 -2
28 | |
29 43--44--45--46--47--48--49--50 -3
30
31 ^
32 -3 -2 -1 X=0 1 2 3 4
33
34 See examples/square-numbers.pl for a simple program printing these
35 numbers.
36
37 Ulam Spiral
38 This path is well known from Stanislaw Ulam finding interesting
39 straight lines when plotting the prime numbers on it.
40
41 Stein, Ulam and Wells, "A Visual Display of Some Properties of the
42 Distribution of Primes", American Mathematical Monthly, volume 71,
43 number 5, May 1964, pages 516-520.
44 <http://www.jstor.org/stable/2312588>
45
46 The cover of Scientific American March 1964 featured this spiral,
47
48 <http://www.nature.com/scientificamerican/journal/v210/n3/covers/index.html>
49
50 <http://oeis.org/A143861/a143861.jpg>
51
52 See examples/ulam-spiral-xpm.pl for a standalone program, or see math-
53 image using this "SquareSpiral" to draw this pattern and more.
54
55 Stein, Ulam and Wells also considered primes on the
56 Math::PlanePath::Corner path, and on a half-plane like two corners.
57
58 Straight Lines
59 The perfect squares 1,4,9,16,25 fall on two diagonals with the even
60 perfect squares going to the upper left and the odd squares to the
61 lower right. The pronic numbers 2,6,12,20,30,42 etc k^2+k half way
62 between the squares fall on similar diagonals to the upper right and
63 lower left. The decagonal numbers 10,27,52,85 etc 4*k^2-3*k go
64 horizontally to the right at Y=-1.
65
66 In general straight lines and diagonals are 4*k^2 + b*k + c. b=0 is
67 the even perfect squares up to the left, then incrementing b is an
68 eighth turn anti-clockwise, or clockwise if negative. So b=1 is
69 horizontal West, b=2 diagonally down South-West, b=3 down South, etc.
70
71 Honaker's prime-generating polynomial 4*k^2 + 4*k + 59 goes down to the
72 right, after the first 30 or so values loop around a bit.
73
74 Wider
75 An optional "wider" parameter makes the path wider, becoming a
76 rectangle spiral instead of a square. For example
77
78 wider => 3
79
80 29--28--27--26--25--24--23--22 2
81 | |
82 30 11--10-- 9-- 8-- 7-- 6 21 1
83 | | | |
84 31 12 1-- 2-- 3-- 4-- 5 20 <- Y=0
85 | | |
86 32 13--14--15--16--17--18--19 -1
87 |
88 33--34--35--36-... -2
89
90 ^
91 -4 -3 -2 -1 X=0 1 2 3
92
93 The centre horizontal 1 to 2 is extended by "wider" many further
94 places, then the path loops around that shape. The starting point 1 is
95 shifted to the left by ceil(wider/2) places to keep the spiral centred
96 on the origin X=0,Y=0.
97
98 Widening doesn't change the nature of the straight lines which arise,
99 it just rotates them around. For example in this wider=3 example the
100 perfect squares are still on diagonals, but the even squares go towards
101 the bottom left (instead of top left when wider=0) and the odd squares
102 to the top right (instead of the bottom right).
103
104 Each loop is still 8 longer than the previous, as the widening is
105 basically a constant amount in each loop.
106
107 N Start
108 The default is to number points starting N=1 as shown above. An
109 optional "n_start" can give a different start with the same shape. For
110 example to start at 0,
111
112 n_start => 0
113
114 16-15-14-13-12 ...
115 | | |
116 17 4--3--2 11 28
117 | | | | |
118 18 5 0--1 10 27
119 | | | |
120 19 6--7--8--9 26
121 | |
122 20-21-22-23-24-25
123
124 The only effect is to push the N values around by a constant amount.
125 It might help match coordinates with something else zero-based.
126
127 Corners
128 Other spirals can be formed by cutting the corners of the square so as
129 to go around faster. See the following modules,
130
131 Corners Cut Class
132 ----------- -----
133 1 HeptSpiralSkewed
134 2 HexSpiralSkewed
135 3 PentSpiralSkewed
136 4 DiamondSpiral
137
138 The "PyramidSpiral" is a re-shaped "SquareSpiral" looping at the same
139 rate. It shifts corners but doesn't cut them.
140
142 See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
143 classes.
144
145 "$path = Math::PlanePath::SquareSpiral->new ()"
146 "$path = Math::PlanePath::SquareSpiral->new (wider => $integer, n_start
147 => $n)"
148 Create and return a new square spiral object. An optional "wider"
149 parameter widens the spiral path, it defaults to 0 which is no
150 widening.
151
152 "($x,$y) = $path->n_to_xy ($n)"
153 Return the X,Y coordinates of point number $n on the path.
154
155 For "$n < 1" the return is an empty list, as the path starts at 1.
156
157 "$n = $path->xy_to_n ($x,$y)"
158 Return the point number for coordinates "$x,$y". $x and $y are
159 each rounded to the nearest integer, which has the effect of
160 treating each N in the path as centred in a square of side 1, so
161 the entire plane is covered.
162
164 N to X,Y
165 There's a few ways to break an N into a side and offset into the side.
166 One convenient way is to treat a loop as starting at the bottom right
167 turn, so N=2,10,26,50,etc, If the first loop at N=2 is reckoned loop
168 number d=1 then the loop starts at
169
170 Nbase = 4*d^2 - 4*d + 2
171 = 2,10,26,50,... for d=1,2,3,4,...
172 (A069894 but it going from d=0)
173
174 For example d=3 is Nbase=4*3^2-4*3+2=26 at X=3,Y=-2. The biggest d
175 with Nbase <= N can be found by inverting with the usual quadratic
176 formula
177
178 d = floor (1/2 + sqrt(N/4 - 1/4))
179
180 For Perl it's good to keep the sqrt argument an integer (when a UV
181 integer is bigger than an NV float, and for BigRat accuracy), so
182 rearranging to
183
184 d = floor ((1+sqrt(N-1)) / 2)
185
186 So Nbase from this d leaves a remainder which is an offset into the
187 loop
188
189 Nrem = N - Nbase
190 = N - (4*d^2 - 4*d + 2)
191
192 The loop starts at X=d,Y=d-1 and has sides length 2d, 2d+1, 2d+1 and
193 2d+2,
194
195 2d
196 +------------+ <- Y=d
197 | |
198 2d | | 2d-1
199 | . |
200 | |
201 | + X=d,Y=-d+1
202 |
203 +---------------+ <- Y=-d
204 2d+1
205
206 ^
207 X=-d
208
209 The X,Y for an Nrem is then
210
211 side Nrem range X,Y result
212 ---- ---------- ----------
213 right Nrem <= 2d-1 X = d
214 Y = -d+1+Nrem
215 top 2d-1 <= Nrem <= 4d-1 X = d-(Nrem-(2d-1)) = 3d-1-Nrem
216 Y = d
217 left 4d-1 <= Nrem <= 6d-1 X = -d
218 Y = d-(Nrem-(4d-1)) = 5d-1-Nrem
219 bottom 6d-1 <= Nrem X = -d+(Nrem-(6d-1)) = -7d+1+Nrem
220 Y = -d
221
222 The corners Nrem=2d-1, Nrem=4d-1 and Nrem=6d-1 get the same result from
223 the two sides that meet so it doesn't matter if the high comparison is
224 "<" or "<=".
225
226 The bottom edge runs through to Nrem < 8d, but there's no need to check
227 that since d=floor(sqrt()) above ensures Nrem is within the loop.
228
229 A small simplification can be had by subtracting an extra 4d-1 from
230 Nrem to make negatives for the right and top sides and positives for
231 the left and bottom.
232
233 Nsig = N - Nbase - (4d-1)
234 = N - (4*d^2 - 4*d + 2) - (4d-1)
235 = N - (4*d^2 + 1)
236
237 side Nsig range X,Y result
238 ---- ---------- ----------
239 right Nsig <= -2d X = d
240 Y = d+(Nsig+2d) = 3d+Nsig
241 top -2d <= Nsig <= 0 X = -d-Nsig
242 Y = d
243 left 0 <= Nsig <= 2d X = -d
244 Y = d-Nsig
245 bottom 2d <= Nsig X = -d+1+(Nsig-(2d+1)) = Nsig-3d
246 Y = -d
247
248 This calculation can be found as an exercise in Graham, Knuth and
249 Patashnik "Concrete Mathematics", chapter 3 "Integer Functions",
250 exercise 40, page 99. They start the spiral from 0, and vertically so
251 their x is -Y here. Their formula for x(n) tests a floor(2*sqrt(N)) to
252 decide whether on a horizontal and so whether to apply the equivalent
253 of Nrem to the result.
254
255 N to X,Y with Wider
256 With the "wider" parameter stretching the spiral loops the formulas
257 above become
258
259 Nbase = 4*d^2 + (-4+2w)*d + 2-w
260
261 d = floor ((2-w + sqrt(4N + w^2 - 4)) / 4)
262
263 Notice for Nbase the w is a term 2*w*d, being an extra 2*w for each
264 loop.
265
266 The left offset ceil(w/2) described above ("Wider") for the N=1
267 starting position is written here as wl, and the other half wr arises
268 too,
269
270 wl = ceil(w/2)
271 wr = floor(w/2) = w - wl
272
273 The horizontal lengths increase by w, and positions shift by wl or wr,
274 but the verticals are unchanged.
275
276 2d+w
277 +------------+ <- Y=d
278 | |
279 2d | | 2d-1
280 | . |
281 | |
282 | + X=d+wr,Y=-d+1
283 |
284 +---------------+ <- Y=-d
285 2d+1+w
286
287 ^
288 X=-d-wl
289
290 The Nsig formulas then have w, wl or wr variously inserted. In all
291 cases if w=wl=wr=0 then they simplify to the plain versions.
292
293 Nsig = N - Nbase - (4d-1+w)
294 = N - ((4d + 2w)*d + 1)
295
296 side Nsig range X,Y result
297 ---- ---------- ----------
298 right Nsig <= -(2d+w) X = d+wr
299 Y = d+(Nsig+2d+w) = 3d+w+Nsig
300 top -(2d+w) <= Nsig <= 0 X = -d-wl-Nsig
301 Y = d
302 left 0 <= Nsig <= 2d X = -d-wl
303 Y = d-Nsig
304 bottom 2d <= Nsig X = -d+1-wl+(Nsig-(2d+1)) = Nsig-wl-3d
305 Y = -d
306
307 Rectangle to N Range
308 Within each row the minimum N is on the X=Y diagonal and N values
309 increases monotonically as X moves away to the left or right.
310 Similarly in each column there's a minimum N on the X=-Y opposite
311 diagonal, or X=-Y+1 diagonal when X negative, and N increases
312 monotonically as Y moves away from there up or down. When wider>0 the
313 location of the minimum changes, but N is still monotonic moving away
314 from the minimum.
315
316 On that basis the maximum N in a rectangle is at one of the four
317 corners,
318
319 |
320 x1,y2 M---|----M x2,y2 corner candidates
321 | | | for maximum N
322 -------O---------
323 | | |
324 | | |
325 x1,y1 M---|----M x1,y1
326 |
327
329 This path is in Sloane's Online Encyclopedia of Integer Sequences in
330 various forms. Summary at
331
332 <http://oeis.org/A068225/a068225.html>
333
334 And various sequences,
335
336 <http://oeis.org/A174344> (etc),
337 <https://oeis.org/wiki/Ulam's_spiral>
338
339 wider=0 (the default)
340 A174344 X coordinate
341 A274923 Y coordinate
342 A214526 abs(X)+abs(Y) "Manhattan" distance
343
344 A079813 abs(dY), being k 0s followed by k 1s
345 A063826 direction 1=right,2=up,3=left,4=down
346
347 A027709 boundary length of N unit squares
348 A078633 grid sticks to make N unit squares
349
350 A033638 N turn positions (extra initial 1, 1)
351 A172979 N turn positions which are primes too
352
353 A054552 N values on X axis (East)
354 A054556 N values on Y axis (North)
355 A054567 N values on negative X axis (West)
356 A033951 N values on negative Y axis (South)
357 A054554 N values on X=Y diagonal (NE)
358 A054569 N values on negative X=Y diagonal (SW)
359 A053755 N values on X=-Y opp diagonal X<=0 (NW)
360 A016754 N values on X=-Y opp diagonal X>=0 (SE)
361 A200975 N values on all four diagonals
362 A317186 N on Y axis positive and negative
363 A267682 N on Y axis positive and negative (origin twice)
364
365 A137928 N values on X=-Y+1 opposite diagonal
366 A002061 N values on X=Y diagonal pos and neg
367 A016814 (4k+1)^2, every second N on south-east diagonal
368
369 A143856 N values on ENE slope dX=2,dY=1
370 A143861 N values on NNE slope dX=1,dY=2
371 A215470 N prime and >=4 primes among its 8 neighbours
372
373 A214664 X coordinate of prime N (Ulam's spiral)
374 A214665 Y coordinate of prime N (Ulam's spiral)
375 A214666 -X \ reckoning spiral starting West
376 A214667 -Y /
377
378 A053999 prime[N] on X=-Y opp diagonal X>=0 (SE)
379 A054551 prime[N] on the X axis (E)
380 A054553 prime[N] on the X=Y diagonal (NE)
381 A054555 prime[N] on the Y axis (N)
382 A054564 prime[N] on X=-Y opp diagonal X<=0 (NW)
383 A054566 prime[N] on negative X axis (W)
384
385 A090925 permutation N at rotate +90
386 A090928 permutation N at rotate +180
387 A090929 permutation N at rotate +270
388 A090930 permutation N at clockwise spiralling
389 A020703 permutation N at rotate +90 and go clockwise
390 A090861 permutation N at rotate +180 and go clockwise
391 A090915 permutation N at rotate +270 and go clockwise
392 A185413 permutation N at 1-X,Y
393 being rotate +180, offset X+1, clockwise
394
395 A068225 permutation N to the N to its right, X+1,Y
396 A121496 run lengths of consecutive N in that permutation
397 A068226 permutation N to the N to its left, X-1,Y
398 A020703 permutation N at transpose Y,X
399 (clockwise <-> anti-clockwise)
400
401 A033952 digits on negative Y axis
402 A033953 digits on negative Y axis, starting 0
403 A033988 digits on negative X axis, starting 0
404 A033989 digits on Y axis, starting 0
405 A033990 digits on X axis, starting 0
406
407 A062410 total sum previous row or column
408
409 wider=1
410 A069894 N on South-West diagonal
411
412 The following have "offset 0" in the OEIS and therefore are based on
413 starting from N=0.
414
415 n_start=0
416 A180714 X+Y coordinate sum
417 A053615 abs(X-Y), runs n to 0 to n, distance to nearest pronic
418
419 A001107 N on X axis
420 A033991 N on Y axis
421 A033954 N on negative Y axis, second 10-gonals
422 A002939 N on X=Y diagonal North-East
423 A016742 N on North-West diagonal, 4*k^2
424 A002943 N on South-West diagonal
425 A156859 N on Y axis positive and negative
426
428 Math::PlanePath, Math::PlanePath::PyramidSpiral
429
430 Math::PlanePath::DiamondSpiral, Math::PlanePath::PentSpiralSkewed,
431 Math::PlanePath::HexSpiralSkewed, Math::PlanePath::HeptSpiralSkewed
432
433 Math::PlanePath::CretanLabyrinth
434
435 Math::NumSeq::SpiroFibonacci
436
437 X11 cursor font "box spiral" cursor which is this style (but going
438 clockwise).
439
441 <http://user42.tuxfamily.org/math-planepath/index.html>
442
444 Copyright 2010, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019
445 Kevin Ryde
446
447 This file is part of Math-PlanePath.
448
449 Math-PlanePath is free software; you can redistribute it and/or modify
450 it under the terms of the GNU General Public License as published by
451 the Free Software Foundation; either version 3, or (at your option) any
452 later version.
453
454 Math-PlanePath is distributed in the hope that it will be useful, but
455 WITHOUT ANY WARRANTY; without even the implied warranty of
456 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
457 General Public License for more details.
458
459 You should have received a copy of the GNU General Public License along
460 with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
461
462
463
464perl v5.30.1 2020-01-30 Math::PlanePath::SquareSpiral(3)