1Math::PlanePath::GosperURseeprliCcoantter(i3b)uted PerlMDaotchu:m:ePnltaanteiPoanth::GosperReplicate(3)
2
6 Math::PlanePath::GosperReplicate -- self-similar hexagon replications
7
9 use Math::PlanePath::GosperReplicate;
10 my $path = Math::PlanePath::GosperReplicate->new;
11 my ($x, $y) = $path->n_to_xy (123);
12
14 This is a self-similar hexagonal tiling of the plane. At each level
15 the shape is the Gosper island.
16 17----16 4
18 / \
19 24----23 18 14----15 3
20 / \ \
21 25 21----22 19----20 10---- 9 2
22 \ / \
23 26----27 3---- 2 11 7---- 8 1
24 / \ \
25 31----30 4 0---- 1 12----13 <- Y=0
26 / \ \
27 32 28----29 5---- 6 45----44 -1
28 \ / \
29 33----34 38----37 46 42----43 -2
30 / \ \
31 39 35----36 47----48 -3
32 \
33 40----41 -4
34 ^
36 -7 -6 -5 -4 -3 -2 -1 X=0 1 2 3 4 5 6 7
37 Points are spread out on every second X coordinate to make a triangular
39 lattice in integer coordinates (see "Triangular Lattice" in
40 Math::PlanePath).
41 The base pattern is the inner N=0 to N=6, then six copies of that shape
43 are arranged around as the blocks N=7,14,21,28,35,42. Then six copies
44 of the resulting N=0 to N=48 shape are replicated around, etc.
45 Each point can be taken as a little hexagon, so that all points tile
47 the plane with hexagons. The innermost N=0 to N=6 are for instance,
48 * *
50 / \ / \
51 / \ / \
52 * * *
53 | 3 | 2 |
54 * * *
55 / \ / \ / \
56 / \ / \ / \
57 * * * *
58 | 4 | 0 | 1 |
59 * * * *
60 \ / \ / \ /
61 \ / \ / \ /
62 * * *
63 | 5 | 6 |
64 * * *
65 \ / \ /
66 \ / \ /
67 * *
68 The further replications are the same arrangement, but the sides become
70 ever wigglier and the centres rotate around. The rotation can be seen
71 N=7 at X=5,Y=1 which is up from the X axis.
72 The "FlowsnakeCentres" path is this same replicating shape, but
74 starting from a side instead of the middle and traversing in such as
75 way as to make each N adjacent. The "Flowsnake" curve itself is this
76 replication too, but segments across hexagons.
77 Complex Base
79 The path corresponds to expressing complex integers X+i*Y in a base
80 b = 5/2 + i*sqrt(3)/2
82 with coordinates scaled to put equilateral triangles on a square grid.
84 So for integer X,Y on the triangular grid (X,Y either both odd or both
85 even),
86 X/2 + i*Y*sqrt(3)/2 = a[n]*b^n + ... + a[2]*b^2 + a[1]*b + a[0]
88 where each digit a[i] is either 0 or a sixth root of unity encoded into
90 base-7 digits of N,
91 w6 = e^(i*pi/3) sixth root of unity, b = 2 + w6
93 = 1/2 + i*sqrt(3)/2
94 N digit a[i] complex number
96 ------- -------------------
97 0 0
98 1 w6^0 = 1
99 2 w6^1 = 1/2 + i*sqrt(3)/2
100 3 w6^2 = -1/2 + i*sqrt(3)/2
101 4 w6^3 = -1
102 5 w6^4 = -1/2 - i*sqrt(3)/2
103 6 w6^5 = 1/2 - i*sqrt(3)/2
104 7 digits suffice because
106 norm(b) = (5/2)^2 + (sqrt(3)/2)^2 = 7
108 Rotate Numbering
110 Parameter "numbering_type => 'rotate'" applies a rotation in each sub-
111 part according to its location around the preceding level.
112 The effect can be illustrated by writing N in base-7. Part 10-16 is
114 the same as the middle 0-6. Part 20-26 has a rotation by +60 degrees.
115 Part 30-36 has rotation by +120 degrees, and so on.
116 22----21
118 / / numbering_type => 'rotate'
119 31 36 23 20 26 N shown in base-7
120 / \ \ \ /
121 32 30 35 24----25 13----12
122 \ / / \
123 33----34 3---- 2 14 10----11
124 / \ \
125 46----45 4 0---- 1 15----16
126 \ \
127 41----40 44 5---- 6 64----63
128 \ / / \
129 42----43 55----54 65 60 62
130 / \ \ \ /
131 56 50 53 66 61
132 / /
133 51----52
134 Notice this means in each part the 11, 21, 31, etc, points are directed
136 away from the middle in the same way, relative to the sub-part
137 locations.
138 Working through the expansions gives the following rule for when an N
140 is on the boundary of level k,
141 write N in k many base-7 digits (empty string if k=0)
143 if any 0 digit then non-boundary
144 ignore high digit and all 1 digits
145 if any 4 or 5 digit then non-boundary
146 if any 32, 33, 66 pair then non-boundary
147 A 0 digit is the middle of a block, or 4 or 5 digit the inner side of a
149 block, for k>=1, hence non-boundary. After that the 6,1,2,3 parts
150 variously expand with rotations so that a 66 is enclosed on the
151 clockwise side and 32 and 33 on the anti-clockwise side.
152
154 See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
155 classes.
156 "$path = Math::PlanePath::GosperReplicate->new ()"
158 "$path = Math::PlanePath::GosperReplicate->new (numbering_type =>
159 $str)"
160 Create and return a new path object. The "numbering_type"
161 parameter can be
162 "fixed" (default)
164 "rotate"
165 "($x,$y) = $path->n_to_xy ($n)"
167 Return the X,Y coordinates of point number $n on the path. Points
168 begin at 0 and if "$n < 0" then the return is an empty list.
169 Level Methods
171 "($n_lo, $n_hi) = $path->level_to_n_range($level)"
172 Return "(0, 7**$level - 1)".
173
175 Axis Rotations
176 In the fixed numbering, digit positions 1,2,3,4,5,6 go around +60deg
177 each, so the N for rotation of X,Y by +60 degrees is each digit +1.
178 N = 0, 1, 2, 3, 4, 5, 6, 10, 11, 12
180 rot+60(N) = 0, 2, 3, 4, 5, 6, 1, 14, 16, 17, ... decimal
182 = 0, 2, 3, 4, 5, 6, 1, 20, 22, 23, ... base7
183 rot+120(N) = 0, 3, 4, 5, 6, 1, 2, 21, 24, 25, ... decimal
185 = 0, 3, 4, 5, 6, 1, 2, 30, 33, 34, ... base7
186 etc
188 In the rotate numbering, just adding +1 (etc) at the high digit alone
190 is rotation.
191 X,Y Extents
193 The maximum X in a given level N=0 to 7^k-1 can be calculated from the
194 replications. A given high digit 1 to 6 has sub-parts located at
195 b^k*w6^(d-1). Those sub-parts are all the same, so the one with
196 maximum real(b^k*w6^(d-1)) contains the maximum X.
197 N_xmax_digit(j) = d=1to6 where real(w6^(d-1) * b^j) is maximum
199 = 1,1,6,6,6,5,5,5,4,4,4,3,3,3,3,2,2, ...
200 k-1
202 N_xmax(k) = digits N_xmax_digit(j) low digit j=0
203 j=0
204 = 0, 1, 8, 302, 2360, 16766, 100801, ... decimal
205 = 0, 1, 11, 611, 6611, 66611, 566611, ... base7
206 k-1
208 z_xmax(k) = sum w6^d[j] * b^j
209 j=0 each d[j] with real(w6^d[j] * b^j) maximum
210 = 0, 1, 7/2+1/2*sqrt3*i, 10-sqrt3*i, 57/2-3/2*sqrt3*i,...
211 xmax(k) = 2*real(z_xmax(k))
213 = 0, 2, 7, 20, 57, 151, 387, 1070, 2833, 7106, ...
214 For computer calculation these maximums can be calculated from the
216 powers. The parts resulting can also be written in terms of the angle
217 arg(b) = atan(sqrt(3)/5) = 19.106... degrees
219 For successive k, if adding this pushes the b^k angle past +30deg then
221 the preceding digit goes past -30deg and becomes the new maximum X.
222 Write the angle as a fraction of 60deg (pi/3),
223 F = atan(sqrt(3)/5) / (pi/3) = 0.318443 ...
225 This is irrational since b^k is never on the X or Y axes. That can be
227 seen since 2/sqrt3*imag(b^k) mod 7 goes in a repeating pattern
228 1,5,4,6,2,3. Similarly 2*real(b^k) mod 7 so not on the Y axis, and
229 also anything on the Y axis would have 3*k fall on the X axis.
230 Digits low to high are successive steps back cyclically 6,5,4,3,2,1 so
232 that (with mod giving 0 to 5),
233 N_xmax_digit(j) = (-floor(F*j+1/2) mod 6) + 1
235 The +1/2 is since initial direction b^0=1 is angle 0 which is half way
237 between -30 and +30 deg.
238 Similarly for the location, using conj(w6) for rotation back
240 z_xmax_exp(j) = floor(F*j+1/2)
242 = 0,0,1,1,1,2,2,2,3,3,3,4,4,4,4,5,5,5, ...
243 z_xmax(k) = sum(j=0,k-1, conj(w6)^z_xmax_exp(j) * b^j)
244 By symmetry the maximum extent is the same in 60deg, 120deg, etc
246 directions, suitably rotated. The N in those cases has the digits
247 1,2,3,4,5,6 cycled around for the rotation. In PlanePath triangular
248 X,Y coordinates direction 60deg means when sum X+3*Y is a maximum, etc.
249 If the +1/2 in the floor is omitted then the effect is to find the
251 maximum point in direction +30deg. In the PlanePath coordinates this
252 means maximum sum S = X+Y.
253 N_smax_digit(j) = (-floor(F*j) mod 6) + 1
255 = 1,1,1,1,6,6,6,5,5,5,4,4,4,3,3, ...
256 k-1
258 N_smax(k) = digits N_smax_digit(j) low digit j=0
259 j=0
260 = 0, 1, 8, 57, 400, 14806, 115648, ... decimal
261 = 0, 1, 11, 111, 1111, 61111, 661111, ... base7
262 and also N_smax() + 1
263 z_smax_exp(j) = floor(F*j)
265 = 0,0,0,0,1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6, ...
266 z_smax(k) = sum(j=0,k-1, conj(w6)^z_smax_exp(j) * b^j)
267 = 0, 1, 7/2+1/2*sqrt3*i, 9+3*sqrt3*i, 19+12*sqrt3*i, ...
268 and also z_smax() + w6^2
269 smax(k) = 2*real(z_smax(k)) + imag(z_smax(k))*2/sqrt3
271 = 0, 2, 8, 24, 62, 172, 470, 1190, 3202, 8740, ...
272 coordinate sum X+Y max
273 In the base figure, points 1 and 2 have the same X+Y=2 and this remains
275 so in subsequent levels, so that for k>=1 N_smax(k) and N_smax(k)+1 are
276 equal maximums.
277
279 Math::PlanePath, Math::PlanePath::GosperIslands,
280 Math::PlanePath::Flowsnake, Math::PlanePath::FlowsnakeCentres,
281 Math::PlanePath::QuintetReplicate, Math::PlanePath::ComplexPlus
282
284 <http://user42.tuxfamily.org/math-planepath/index.html>
285
287 Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020
288 Kevin Ryde
289 This file is part of Math-PlanePath.
291 Math-PlanePath is free software; you can redistribute it and/or modify
293 it under the terms of the GNU General Public License as published by
294 the Free Software Foundation; either version 3, or (at your option) any
295 later version.
296 Math-PlanePath is distributed in the hope that it will be useful, but
298 WITHOUT ANY WARRANTY; without even the implied warranty of
299 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
300 General Public License for more details.
301 You should have received a copy of the GNU General Public License along
303 with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
304perl v5.38.0 2023-07-20Math::PlanePath::GosperReplicate(3)