1Math::PlanePath::CompleUxsMeirnuCso(n3t)ributed Perl DocMuamtehn:t:aPtliaonnePath::ComplexMinus(3)
2
6 Math::PlanePath::ComplexMinus -- i-1 and other complex number bases i-r
7
9 use Math::PlanePath::ComplexMinus;
10 my $path = Math::PlanePath::ComplexMinus->new (realpart => 1);
11 my ($x, $y) = $path->n_to_xy (123);
12
14 This path traverses points by a complex number base i-r for given
15 integer r. The default is base i-1 as per
16 Solomon I. Khmelnik "Specialized Digital Computer for Operations
18 with Complex Numbers" (in Russian), Questions of Radio Electronics,
19 volume 12, number 2, 1964.
20 <http://lib.izdatelstwo.com/Papers2/s4.djvu>
21 Walter Penney, "A 'Binary' System for Complex Numbers", Journal of
23 the ACM, volume 12, number 2, April 1965, pages 247-248.
24 <http://www.nsa.gov/Portals/70/documents/news-features/declassified-documents/tech-journals/a-binary-system.pdf>
25 When continued to a power-of-2 extent, this is sometimes called the
27 "twindragon" since the shape (and fractal limit) corresponds to two
28 DragonCurve back-to-back. But the numbering of points in the
29 twindragon is not the same as here.
30 26 27 10 11 3
32 24 25 8 9 2
33 18 19 30 31 2 3 14 15 1
34 16 17 28 29 0 1 12 13 <- Y=0
35 22 23 6 7 58 59 42 43 -1
36 20 21 4 5 56 57 40 41 -2
37 50 51 62 63 34 35 46 47 -3
38 48 49 60 61 32 33 44 45 -4
39 54 55 38 39 -5
40 52 53 36 37 -6
41 ^
43 -5 -4 -3 -2 -1 X=0 1 2 3 4 5 6 7
44 A complex integer can be represented as a set of powers,
46 X+Yi = a[n]*b^n + ... + a[2]*b^2 + a[1]*b + a[0]
48 base b=i-1
49 digits a[n] to a[0] each = 0 or 1
50 N = a[n]*2^n + ... + a[2]*2^2 + a[1]*2 + a[0]
52 N is the a[i] digits as bits and X,Y is the resulting complex number.
54 It can be shown that this is a one-to-one mapping so every integer X,Y
55 of the plane is visited once each.
56 The shape of points N=0 to 2^level-1 repeats as N=2^level to
58 2^(level+1)-1. For example N=0 to N=7 is repeated as N=8 to N=15, but
59 starting at position X=2,Y=2 instead of the origin. That position 2,2
60 is because b^3 = 2+2i. There's no rotations or mirroring etc in this
61 replication, just position offsets.
62 N=0 to N=7 N=8 to N=15 repeat shape
64 2 3 10 11
66 0 1 8 9
67 6 7 14 15
68 4 5 12 13
69 For b=i-1 each N=2^level point starts at X+Yi=(i-1)^level. The
71 powering of that b means the start position rotates around by +135
72 degrees each time and outward by a radius factor sqrt(2) each time. So
73 for example b^3 = 2+2i is followed by b^4 = -4, which is 135 degrees
74 around and radius |b^3|=sqrt(8) becomes |b^4|=sqrt(16).
75 Real Part
77 The "realpart => $r" option gives a complex base b=i-r for a given
78 integer r>=1. For example "realpart => 2" is
79 20 21 22 23 24 4
81 15 16 17 18 19 3
82 10 11 12 13 14 2
83 5 6 7 8 9 1
84 45 46 47 48 49 0 1 2 3 4 <- Y=0
85 40 41 42 43 44 -1
86 35 36 37 38 39 -2
87 30 31 32 33 34 -3
88 70 71 72 73 74 25 26 27 28 29 -4
89 65 66 67 68 69 -5
90 60 61 62 63 64 -6
91 55 56 57 58 59 -7
92 50 51 52 53 54 -8
93 ^
94 -8 -7 -6 -5 -4 -3 -2 -1 X=0 1 2 3 4 5 6 7 8 9 10
95 N is broken into digits of base=norm=r*r+1, ie. digits 0 to r*r
97 inclusive. This makes horizontal runs of r*r+1 many points, such as
98 N=5 to N=9 etc above. In the default r=1 these runs are 2 long whereas
99 for r=2 they're 2*2+1=5 long, or r=3 would be 3*3+1=10, etc.
100 The offset back for each run like N=5 shown is the r in i-r, then the
102 next level is (i-r)^2 = (-2r*i + r^2-1) so N=25 begins at Y=-2*2=-4,
103 X=2*2-1=3.
104 The successive replications tile the plane for any r, though the N
106 values needed to rotate all the way around become big if norm=r*r+1 is
107 big.
108 Fractal
110 The i-1 twindragon is usually conceived as taking fractional N like
111 0.abcde in binary and giving fractional complex X+iY. The twindragon
112 is then all the points of the complex plane reached by such fractional
113 N. This set of points can be shown to be connected and to fill a
114 certain radius around the origin.
115 The code here might be pressed into use for that to some finite number
117 of bits by multiplying up to make an integer N
118 Nint = Nfrac * 256^k
120 Xfrac = Xint / 16^k
121 Yfrac = Yint / 16^k
122 256 is a good power because b^8=16 is a positive real and so there's no
124 rotations to apply to the resulting X,Y, only a power-of-16 division
125 (b^8)^k=16^k each. Using b^4=-4 for a multiplier 16^k and divisor
126 (-4)^k would be almost as easy too, requiring just sign changes if k
127 odd.
128
130 See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
131 classes.
132 "$path = Math::PlanePath::ComplexMinus->new ()"
134 "$path = Math::PlanePath::ComplexMinus->new (realpart => $r)"
135 Create and return a new path object.
136 "($x,$y) = $path->n_to_xy ($n)"
138 Return the X,Y coordinates of point number $n on the path. Points
139 begin at 0 and if "$n < 0" then the return is an empty list.
140 $n should be an integer, it's unspecified yet what will be done for
142 a fraction.
143 Level Methods
145 "($n_lo, $n_hi) = $path->level_to_n_range($level)"
146 Return "(0, 2**$level - 1)", or with "realpart" option return "(0,
147 $norm**$level - 1)" where norm=realpart^2+1.
148
150 Various formulas and pictures etc for the i-1 case can be found in the
151 author's long mathematical write-up (section "Complex Base i-1")
152 <http://user42.tuxfamily.org/dragon/index.html>
154 X,Y to N
156 A given X,Y representing X+Yi can be turned into digits of N by
157 successive complex divisions by i-r. Each digit of N is a real
158 remainder 0 to r*r inclusive from that division.
159 The base formula above is
161 X+Yi = a[n]*b^n + ... + a[2]*b^2 + a[1]*b + a[0]
163 and want the a[0]=digit to be a real 0 to r*r inclusive. Subtracting
165 a[0] and dividing by b will give
166 (X+Yi - digit) / (i-r)
168 = - (X-digit + Y*i) * (i+r) / norm
169 = (Y - (X-digit)*r)/norm
170 + i * - ((X-digit) + Y*r)/norm
171 which is
173 Xnew = Y - (X-digit)*r)/norm
175 Ynew = -((X-digit) + Y*r)/norm
176 The a[0] digit must make both Xnew and Ynew parts integers. The
178 easiest one to calculate from is the imaginary part, from which require
179 - ((X-digit) + Y*r) == 0 mod norm
181 so
183 digit = X + Y*r mod norm
185 This digit value makes the real part a multiple of norm too, as can be
187 seen from
188 Xnew = Y - (X-digit)*r
190 = Y - X*r - (X+Y*r)*r
191 = Y - X*r - X*r + Y*r*r
192 = Y*(r*r+1)
193 = Y*norm
194 Notice Ynew is the quotient from (X+Y*r)/norm rounded downwards
196 (towards negative infinity). Ie. in the division "X+Y*r mod norm"
197 which calculates the digit, the quotient is Ynew and the remainder is
198 the digit.
199 X Axis N for Realpart 1
201 For base i-1, Penney shows the N on the X axis are
202 X axis N in hexadecimal uses only digits 0, 1, C, D
204 = 0, 1, 12, 13, 16, 17, 28, 29, 192, 193, 204, 205, 208, ...
205 Those on the positive X axis have an odd number of digits and on the X
207 negative axis an even number of digits.
208 To be on the X axis the imaginary parts of the base powers b^k must
210 cancel out to leave just a real part. The powers repeat in an 8-long
211 cycle
212 k b^k for b=i-1
214 0 +1
215 1 i -1
216 2 -2i +0 \ pair cancel
217 3 2i +2 /
218 4 -4
219 5 -4i +4
220 6 8i +0 \ pair cancel
221 7 -8i -8 /
222 The k=0 and k=4 bits are always reals and can always be included. Bits
224 k=2 and k=3 have imaginary parts -2i and 2i which cancel out, so they
225 can be included together. Similarly k=6 and k=7 with 8i and -8i. The
226 two blocks k=0to3 and k=4to7 differ only in a negation so the bits can
227 be reckoned in groups of 4, which is hexadecimal. Bit 1 is digit value
228 1 and bits 2,3 together are digit value 0xC, so adding one or both of
229 those gives combinations are 0,1,0xC,0xD.
230 The high hex digit determines the sign, positive or negative, of the
232 total real part. Bits k=0 or k=2,3 are positive. Bits k=4 or k=6,7
233 are negative, so
234 N for X>0 N for X<0
236 0x01.. 0x1_.. even number of hex 0,1,C,D following
238 0x0C.. 0xC_.. "_" digit any of 0,1,C,D
239 0x0D.. 0xD_..
240 which is equivalent to X>0 is an odd number of hex digits or X<0 is an
242 even number. For example N=28=0x1C is at X=-2 since that N is X<0 form
243 "0x1_".
244 The order of the values on the positive X axis is obtained by taking
246 the digits in reverse order on alternate positions
247 0,1,C,D high digit
249 D,C,1,0
250 0,1,C,D
251 ...
252 D,C,1,0
253 0,1,C,D low digit
254 For example in the following notice the first and third digit
256 increases, but the middle digit decreases,
257 X=4to7 N=0x1D0,0x1D1,0x1DC,0x1DD
259 X=8to11 N=0x1C0,0x1C1,0x1CC,0x1CD
260 X=12to15 N=0x110,0x111,0x11C,0x11D
261 X=16to19 N=0x100,0x101,0x10C,0x10D
262 X=20to23 N=0xCD0,0xCD1,0xCDC,0xCDD
263 For the negative X axis it's the same if reading by increasing X, ie.
265 upwards toward +infinity, or the opposite way around if reading
266 decreasing X, ie. more negative downwards toward -infinity.
267 Y Axis N for Realpart 1
269 For base i-1 Penny also characterises the N values on the Y axis,
270 Y axis N in base-64 uses only
272 at even digits 0, 3, 4, 7, 48, 51, 52, 55
273 at odd digit 0, 1, 12, 13, 16, 17, 28, 29
274 = 0,3,4,7,48,51,52,55,64,67,68,71,112,115,116,119, ...
276 Base-64 means taking N in 6-bit blocks. Digit positions are counted
278 starting from the least significant digit as position 0 which is even.
279 So the low digit can be only 0,3,4,etc, then the second digit only
280 0,1,12,etc, and so on.
281 This arises from (i-1)^6 = 8i which gives a repeating pattern of 6-bit
283 blocks. The different patterns at odd and even positions are since i^2
284 = -1.
285 Boundary Length
287 The length of the boundary of unit squares for the first norm^k many
288 points, ie. N=0 to N=norm^k-1 inclusive, is calculated in
289 William J. Gilbert, "The Fractal Dimension of Sets Derived From
291 Complex Bases", Canadian Mathematical Bulletin, volume 29, number
292 4, 1986.
293 <http://www.math.uwaterloo.ca/~wgilbert/Research/GilbertFracDim.pdf>
294 The boundary formula is a 3rd-order recurrence. For the twindragon
296 case it is
297 for realpart=1
299 boundary[k] = boundary[k-1] + 2*boundary[k-3]
300 = 4, 6, 10, 18, 30, 50, 86, 146, 246, 418, ... (2*A003476)
301 4 + 2*x + 4*x^2
303 generating function ---------------
304 1 - x - 2*x^3
305 The first three boundaries are as follows. Then the recurrence gives
307 the next boundary[3] = 10+2*4 = 18.
308 k area boundary[k]
310 --- ---- -----------
311 +---+
312 0 2^k = 1 4 | 0 |
313 +---+
314 +---+---+
316 1 2^k = 2 6 | 0 1 |
317 +---+---+
318 +---+---+
320 | 2 3 |
321 2 2^k = 4 10 +---+ +---+
322 | 0 1 |
323 +---+---+
324 Gilbert calculates for any i-r by taking the boundary in three parts
326 A,B,C and showing how in the next replication level those boundary
327 parts transform into multiple copies of the preceding level parts. The
328 replication is easier to visualize for a bigger "r" than for i-1
329 because in bigger r it's clearer how the A, B and C parts differ. The
330 length replications are
331 A -> A * (2*r-1) + C * 2*r
333 B -> A * (r^2-2*r+2) + C * (r-1)^2
334 C -> B
335 starting from
337 A = 2*r
338 B = 2
339 C = 2 - 2*r
340 total boundary = A+B+C
342 For i-1 realpart=1 these A,B,C are already in the form of a recurrence
344 A->A+2*C, B->A, C->B, per the formula above. For other real parts a
345 little matrix rearrangement turns the A,B,C parts into recurrence
346 boundary[k] = boundary[k-1] * (2*r - 1)
348 + boundary[k-2] * (norm - 2*r)
349 + boundary[k-3] * norm
350 starting from
352 boundary[0] = 4 # single square cell
353 boundary[1] = 2*norm + 2 # oblong of norm many cells
354 boundary[2] = 2*(norm-1)*(r+2) + 4
355 For example
357 for realpart=2
359 boundary[k] = 3*boundary[k-1] + 1*boundary[k-2] + 5*boundary[k-3]
360 = 4, 12, 36, 140, 516, 1868, 6820, 24908, ...
361 4 - 4*x^2
363 generating function ---------------------
364 1 - 3*x - x^2 - 5*x^3
365 If calculating for large k values then the matrix form can be powered
367 up rather than repeated additions. (As usual for all such linear
368 recurrences.)
369
371 Entries in Sloane's Online Encyclopedia of Integer Sequences related to
372 this path include
373 <http://oeis.org/A318438> (etc)
375 realpart=1 (base i-1, the default)
377 A318438 X coordinate
378 A318439 Y coordinate
379 A318479 norm X^2 + Y^2
380 A066321 N on X>=0 axis, or N/2 of North-West diagonal
381 A271472 same in binary
382 A066323 number of 1 bits
383 A256441 N on negative X axis, X<=0
384 A073791 X axis X sorted by N
385 A320283 Y axis Y sorted by N
386 A137426 dX/2 at N=2^(k+2)-1
388 dY at N=2^k-1 (step to next replication level)
389 A066322 diffs (N at X=16k+4) - (N at X=16k+3)
390 A340566 permutation N by diagonals +/-
391 A003476 boundary length / 2
393 recurrence a(n) = a(n-1) + 2*a(n-3)
394 A203175 boundary length, starting from 4
395 (believe its conjectured recurrence is true)
396 A052537 boundary length part A, B or C, per Gilbert's paper
397 A193239 reverse-add steps to N binary palindrome
399 A193240 reverse-add trajectory of binary 110
400 A193241 reverse-add trajectory of binary 10110
401 A193306 reverse-subtract steps to 0 (plain-rev)
402 A193307 reverse-subtract steps to 0 (rev-plain)
403 realpart=1 (base i-2)
405 A011658 turn 0=straight, 1=not straight, repeating 0,0,0,1,1
406
408 Math::PlanePath, Math::PlanePath::DragonCurve,
409 Math::PlanePath::ComplexPlus
410
412 <http://user42.tuxfamily.org/math-planepath/index.html>
413
415 Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020,
416 2021 Kevin Ryde
417 Math-PlanePath is free software; you can redistribute it and/or modify
419 it under the terms of the GNU General Public License as published by
420 the Free Software Foundation; either version 3, or (at your option) any
421 later version.
422 Math-PlanePath is distributed in the hope that it will be useful, but
424 WITHOUT ANY WARRANTY; without even the implied warranty of
425 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
426 General Public License for more details.
427 You should have received a copy of the GNU General Public License along
429 with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
430perl v5.36.0 2022-07-22 Math::PlanePath::ComplexMinus(3)