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 Penny, "A 'Binary' System for Complex Numbers", Journal of
23 the ACM, volume 12, number 2, April 1965, pages 247-248.
24 When continued to a power-of-2 extent this is sometimes called the
26 "twindragon" since the shape (and fractal limit) limit correspond to
27 two DragonCurve back-to-back. But the numbering of points in the
28 twindragon is not the same as here.
29 26 27 10 11 3
31 24 25 8 9 2
32 18 19 30 31 2 3 14 15 1
33 16 17 28 29 0 1 12 13 <- Y=0
34 22 23 6 7 58 59 42 43 -1
35 20 21 4 5 56 57 40 41 -2
36 50 51 62 63 34 35 46 47 -3
37 48 49 60 61 32 33 44 45 -4
38 54 55 38 39 -5
39 52 53 36 37 -6
40 ^
42 -5 -4 -3 -2 -1 X=0 1 2 3 4 5 6 7
43 A complex integer can be represented as a set of powers,
45 X+Yi = a[n]*b^n + ... + a[2]*b^2 + a[1]*b + a[0]
47 base b=i-1
48 digits a[n] to a[0] each = 0 or 1
49 N = a[n]*2^n + ... + a[2]*2^2 + a[1]*2 + a[0]
51 N is the a[i] digits as bits and X,Y is the resulting complex number.
53 It can be shown that this is a one-to-one mapping so every integer X,Y
54 of the plane is visited once each.
55 The shape of points N=0 to 2^level-1 repeats as N=2^level to
57 2^(level+1)-1. For example N=0 to N=7 is repeated as N=8 to N=15, but
58 starting at position X=2,Y=2 instead of the origin. That position 2,2
59 is because b^3 = 2+2i. There's no rotations or mirroring etc in this
60 replication, just position offsets.
61 N=0 to N=7 N=8 to N=15 repeat shape
63 2 3 10 11
65 0 1 8 9
66 6 7 14 15
67 4 5 12 13
68 For b=i-1 each N=2^level point starts at X+Yi=(i-1)^level. The
70 powering of that b means the start position rotates around by +135
71 degrees each time and outward by a radius factor sqrt(2) each time. So
72 for example b^3 = 2+2i is followed by b^4 = -4, which is 135 degrees
73 around and radius |b^3|=sqrt(8) becomes |b^4|=sqrt(16).
74 Real Part
76 The "realpart => $r" option gives a complex base b=i-r for a given
77 integer r>=1. For example "realpart => 2" is
78 20 21 22 23 24 4
80 15 16 17 18 19 3
81 10 11 12 13 14 2
82 5 6 7 8 9 1
83 45 46 47 48 49 0 1 2 3 4 <- Y=0
84 40 41 42 43 44 -1
85 35 36 37 38 39 -2
86 30 31 32 33 34 -3
87 70 71 72 73 74 25 26 27 28 29 -4
88 65 66 67 68 69 -5
89 60 61 62 63 64 -6
90 55 56 57 58 59 -7
91 50 51 52 53 54 -8
92 ^
93 -8 -7 -6 -5 -4 -3 -2 -1 X=0 1 2 3 4 5 6 7 8 9 10
94 N is broken into digits of base=norm=r*r+1, ie. digits 0 to r*r
96 inclusive. This makes horizontal runs of r*r+1 many points, such as
97 N=5 to N=9 etc above. In the default r=1 these runs are 2 long whereas
98 for r=2 they're 2*2+1=5 long, or r=3 would be 3*3+1=10, etc.
99 The offset back for each run like N=5 shown is the r in i-r, then the
101 next level is (i-r)^2 = (-2r*i + r^2-1) so N=25 begins at Y=-2*2=-4,
102 X=2*2-1=3.
103 The successive replications tile the plane for any r, though the N
105 values needed to rotate around and do so become large if norm=r*r+1 is
106 large.
107 Fractal
109 The i-1 twindragon is usually conceived as taking fractional N like
110 0.abcde in binary and giving fractional complex X+iY. The twindragon
111 is then all the points of the complex plane reached by such fractional
112 N. This set of points can be shown to be connected and to fill a
113 certain radius around the origin.
114 The code here might be pressed into use for that to some finite number
116 of bits by multiplying up to make an integer N
117 Nint = Nfrac * 256^k
119 Xfrac = Xint / 16^k
120 Yfrac = Yint / 16^k
121 256 is a good power because b^8=16 is a positive real and so there's no
123 rotations to apply to the resulting X,Y, only a power-of-16 division
124 (b^8)^k=16^k each. Using b^4=-4 for a multiplier 16^k and divisor
125 (-4)^k would be almost as easy too, requiring just sign changes if k
126 odd.
127
129 See "FUNCTIONS" in Math::PlanePath for behaviour common to all path
130 classes.
131 "$path = Math::PlanePath::ComplexMinus->new ()"
133 "$path = Math::PlanePath::ComplexMinus->new (realpart => $r)"
134 Create and return a new path object.
135 "($x,$y) = $path->n_to_xy ($n)"
137 Return the X,Y coordinates of point number $n on the path. Points
138 begin at 0 and if "$n < 0" then the return is an empty list.
139 $n should be an integer, it's unspecified yet what will be done for
141 a fraction.
142 Level Methods
144 "($n_lo, $n_hi) = $path->level_to_n_range($level)"
145 Return "(0, 2**$level - 1)", or with "realpart" option return "(0,
146 $norm**$level - 1)" where norm=realpart^2+1.
147
149 Various formulas and pictures etc for the i-1 case can be found in the
150 author's long mathematical write-up (section "Complex Base i-1")
151 <http://user42.tuxfamily.org/dragon/index.html>
153 X,Y to N
155 A given X,Y representing X+Yi can be turned into digits of N by
156 successive complex divisions by i-r. Each digit of N is a real
157 remainder 0 to r*r inclusive from that division.
158 The base formula above is
160 X+Yi = a[n]*b^n + ... + a[2]*b^2 + a[1]*b + a[0]
162 and want the a[0]=digit to be a real 0 to r*r inclusive. Subtracting
164 a[0] and dividing by b will give
165 (X+Yi - digit) / (i-r)
167 = - (X-digit + Y*i) * (i+r) / norm
168 = (Y - (X-digit)*r)/norm
169 + i * - ((X-digit) + Y*r)/norm
170 which is
172 Xnew = Y - (X-digit)*r)/norm
174 Ynew = -((X-digit) + Y*r)/norm
175 The a[0] digit must make both Xnew and Ynew parts integers. The
177 easiest one to calculate from is the imaginary part, from which require
178 - ((X-digit) + Y*r) == 0 mod norm
180 so
182 digit = X + Y*r mod norm
184 This digit value makes the real part a multiple of norm too, as can be
186 seen from
187 Xnew = Y - (X-digit)*r
189 = Y - X*r - (X+Y*r)*r
190 = Y - X*r - X*r + Y*r*r
191 = Y*(r*r+1)
192 = Y*norm
193 Notice Ynew is the quotient from (X+Y*r)/norm rounded downwards
195 (towards negative infinity). Ie. in the division "X+Y*r mod norm"
196 which calculates the digit, the quotient is Ynew and the remainder is
197 the digit.
198 X Axis N for Realpart 1
200 For base i-1, Penney shows the N on the X axis are
201 X axis N in hexadecimal uses only digits 0, 1, C, D
203 = 0, 1, 12, 13, 16, 17, 28, 29, 192, 193, 204, 205, 208, ...
204 Those on the positive X axis have an odd number of digits and on the X
206 negative axis an even number of digits.
207 To be on the X axis the imaginary parts of the base powers b^k must
209 cancel out to leave just a real part. The powers repeat in an 8-long
210 cycle
211 k b^k for b=i-1
213 0 +1
214 1 i -1
215 2 -2i +0 \ pair cancel
216 3 2i +2 /
217 4 -4
218 5 -4i +4
219 6 8i +0 \ pair cancel
220 7 -8i -8 /
221 The k=0 and k=4 bits are always reals and can always be included. Bits
223 k=2 and k=3 have imaginary parts -2i and 2i which cancel out, so they
224 can be included together. Similarly k=6 and k=7 with 8i and -8i. The
225 two blocks k=0to3 and k=4to7 differ only in a negation so the bits can
226 be reckoned in groups of 4, which is hexadecimal. Bit 1 is digit value
227 1 and bits 2,3 together are digit value 0xC, so adding one or both of
228 those gives combinations are 0,1,0xC,0xD.
229 The high hex digit determines the sign, positive or negative, of the
231 total real part. Bits k=0 or k=2,3 are positive. Bits k=4 or k=6,7
232 are negative, so
233 N for X>0 N for X<0
235 0x01.. 0x1_.. even number of hex 0,1,C,D following
237 0x0C.. 0xC_.. "_" digit any of 0,1,C,D
238 0x0D.. 0xD_..
239 which is equivalent to X>0 is an odd number of hex digits or X<0 is an
241 even number. For example N=28=0x1C is at X=-2 since that N is X<0 form
242 "0x1_".
243 The order of the values on the positive X axis is obtained by taking
245 the digits in reverse order on alternate positions
246 0,1,C,D high digit
248 D,C,1,0
249 0,1,C,D
250 ...
251 D,C,1,0
252 0,1,C,D low digit
253 For example in the following notice the first and third digit
255 increases, but the middle digit decreases,
256 X=4to7 N=0x1D0,0x1D1,0x1DC,0x1DD
258 X=8to11 N=0x1C0,0x1C1,0x1CC,0x1CD
259 X=12to15 N=0x110,0x111,0x11C,0x11D
260 X=16to19 N=0x100,0x101,0x10C,0x10D
261 X=20to23 N=0xCD0,0xCD1,0xCDC,0xCDD
262 For the negative X axis it's the same if reading by increasing X, ie.
264 upwards toward +infinity, or the opposite way around if reading
265 decreasing X, ie. more negative downwards toward -infinity.
266 Y Axis N for Realpart 1
268 For base i-1 Penny also characterises the N values on the Y axis,
269 Y axis N in base-64 uses only
271 at even digits 0, 3, 4, 7, 48, 51, 52, 55
272 at odd digit 0, 1, 12, 13, 16, 17, 28, 29
273 = 0,3,4,7,48,51,52,55,64,67,68,71,112,115,116,119, ...
275 Base-64 means taking N in 6-bit blocks. Digit positions are counted
277 starting from the least significant digit as position 0 which is even.
278 So the low digit can be only 0,3,4,etc, then the second digit only
279 0,1,12,etc, and so on.
280 This arises from (i-1)^6 = 8i which gives a repeating pattern of 6-bit
282 blocks. The different patterns at odd and even positions are since i^2
283 = -1.
284 Boundary Length
286 The length of the boundary of unit squares for the first norm^k many
287 points, ie. N=0 to N=norm^k-1 inclusive, is calculated in
288 William J. Gilbert, "The Fractal Dimension of Sets Derived From
290 Complex Bases", Canadian Mathematical Bulletin, volume 29, number
291 4, 1986.
292 <http://www.math.uwaterloo.ca/~wgilbert/Research/GilbertFracDim.pdf>
293 The boundary formula is a 3rd-order recurrence. For the twindragon
295 case it is
296 for realpart=1
298 boundary[k] = boundary[k-1] + 2*boundary[k-3]
299 = 4, 6, 10, 18, 30, 50, 86, 146, 246, 418, ... (2*A003476)
300 4 + 2*x + 4*x^2
302 generating function ---------------
303 1 - x - 2*x^3
304 The first three boundaries are as follows. Then the recurrence gives
306 the next boundary[3] = 10+2*4 = 18.
307 k area boundary[k]
309 --- ---- -----------
310 +---+
311 0 2^k = 1 4 | 0 |
312 +---+
313 +---+---+
315 1 2^k = 2 6 | 0 1 |
316 +---+---+
317 +---+---+
319 | 2 3 |
320 2 2^k = 4 10 +---+ +---+
321 | 0 1 |
322 +---+---+
323 Gilbert calculates for any i-r by taking the boundary in three parts
325 A,B,C and showing how in the next replication level those boundary
326 parts transform into multiple copies of the preceding level parts. The
327 replication is easier to visualize for a bigger "r" than for i-1
328 because in bigger r it's clearer how the A, B and C parts differ. The
329 length replications are
330 A -> A * (2*r-1) + C * 2*r
332 B -> A * (r^2-2*r+2) + C * (r-1)^2
333 C -> B
334 starting from
336 A = 2*r
337 B = 2
338 C = 2 - 2*r
339 total boundary = A+B+C
341 For i-1 realpart=1 these A,B,C are already in the form of a recurrence
343 A->A+2*C, B->A, C->B, per the formula above. For other real parts a
344 little matrix rearrangement turns the A,B,C parts into recurrence
345 boundary[k] = boundary[k-1] * (2*r - 1)
347 + boundary[k-2] * (norm - 2*r)
348 + boundary[k-3] * norm
349 starting from
351 boundary[0] = 4 # single square cell
352 boundary[1] = 2*norm + 2 # oblong of norm many cells
353 boundary[2] = 2*(norm-1)*(r+2) + 4
354 For example
356 for realpart=2
358 boundary[k] = 3*boundary[k-1] + 1*boundary[k-2] + 5*boundary[k-3]
359 = 4, 12, 36, 140, 516, 1868, 6820, 24908, ...
360 4 - 4*x^2
362 generating function ---------------------
363 1 - 3*x - x^2 - 5*x^3
364 If calculating for large k values then the matrix form can be powered
366 up rather than repeated additions. (As usual for all such linear
367 recurrences.)
368
370 Entries in Sloane's Online Encyclopedia of Integer Sequences related to
371 this path include
372 <http://oeis.org/A066321> (etc)
374 realpart=1 (base i-1, the default)
376 A066321 N on X axis, X>=0
377 A271472 and in binary
378 A066322 diffs (N at X=16k+4) - (N at X=16k+3)
379 A066323 N on X axis, number of 1 bits
380 A137426 dX/2 at N=2^(k+2)-1
381 dY at N=2^k-1 (step to next replication level)
382 A256441 N on negative X axis, X<=0
383 A003476 boundary length / 2
385 recurrence a(n) = a(n-1) + 2*a(n-3)
386 A203175 boundary length, starting from 4
387 (believe its conjectured recurrence is true)
388 A052537 boundary length part A, B or C, per Gilbert's paper
389 A193239 reverse-add steps to N binary palindrome
391 A193240 reverse-add trajectory of binary 110
392 A193241 reverse-add trajectory of binary 10110
393 A193306 reverse-subtract steps to 0 (plain-rev)
394 A193307 reverse-subtract steps to 0 (rev-plain)
395
397 Math::PlanePath, Math::PlanePath::DragonCurve,
398 Math::PlanePath::ComplexPlus
399
401 <http://user42.tuxfamily.org/math-planepath/index.html>
402
404 Copyright 2011, 2012, 2013, 2014, 2015, 2016, 2017 Kevin Ryde
405 Math-PlanePath is free software; you can redistribute it and/or modify
407 it under the terms of the GNU General Public License as published by
408 the Free Software Foundation; either version 3, or (at your option) any
409 later version.
410 Math-PlanePath is distributed in the hope that it will be useful, but
412 WITHOUT ANY WARRANTY; without even the implied warranty of
413 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
414 General Public License for more details.
415 You should have received a copy of the GNU General Public License along
417 with Math-PlanePath. If not, see <http://www.gnu.org/licenses/>.
418perl v5.28.0 2017-12-03 Math::PlanePath::ComplexMinus(3)