Provided by: libmath-planepath-perl_122-1_all 
NAME
Math::PlanePath::ComplexMinus -- twindragon and other complex number base i-r
SYNOPSIS
use Math::PlanePath::ComplexMinus;
my $path = Math::PlanePath::ComplexMinus->new (realpart=>1);
my ($x, $y) = $path->n_to_xy (123);
DESCRIPTION
This path traverses points by a complex number base i-r for given integer r. The default
is base i-1 as per
Walter Penny, "A 'Binary' System for Complex Numbers", Journal of the ACM, volume 12,
number 2, April 1965, pages 247-248.
When continued to a power-of-2 extent this is called the "twindragon" shape.
26 27 10 11 3
24 25 8 9 2
18 19 30 31 2 3 14 15 1
16 17 28 29 0 1 12 13 <- Y=0
22 23 6 7 58 59 42 43 -1
20 21 4 5 56 57 40 41 -2
50 51 62 63 34 35 46 47 -3
48 49 60 61 32 33 44 45 -4
54 55 38 39 -5
52 53 36 37 -6
^
-5 -4 -3 -2 -1 X=0 1 2 3 4 5 6 7
A complex integer can be represented as a set of powers,
X+Yi = a[n]*b^n + ... + a[2]*b^2 + a[1]*b + a[0]
base b=i-1
digits a[n] to a[0] each = 0 or 1
N = a[n]*2^n + ... + a[2]*2^2 + a[1]*2 + a[0]
N is the a[i] digits as bits and X,Y is the resulting complex number. It can be shown
that this is a one-to-one mapping so every integer X,Y of the plane is visited once each.
The shape of points N=0 to N=2^level-1 repeats as N=2^level to N=2^(level+1)-1. For
example N=0 to N=7 is repeated as N=8 to N=15, but starting at position X=2,Y=2 instead of
the origin. That position 2,2 is because b^3 = 2+2i. There's no rotations or mirroring
etc in this replication, just position offsets.
N=0 to N=7 N=8 to N=15 repeat shape
2 3 10 11
0 1 8 9
6 7 14 15
4 5 12 13
For b=i-1 each N=2^level point starts at X+Yi=(i-1)^level. The powering of that b means
the start position rotates around by +135 degrees each time and outward by a radius factor
sqrt(2) each time. So for example b^3 = 2+2i is followed by b^4 = -4, which is 135
degrees around and radius |b^3|=sqrt(8) becomes |b^4|=sqrt(16).
Real Part
The "realpart => $r" option gives a complex base b=i-r for a given integer r>=1. For
example "realpart => 2" is
20 21 22 23 24 4
15 16 17 18 19 3
10 11 12 13 14 2
5 6 7 8 9 1
45 46 47 48 49 0 1 2 3 4 <- Y=0
40 41 42 43 44 -1
35 36 37 38 39 -2
30 31 32 33 34 -3
70 71 72 73 74 25 26 27 28 29 -4
65 66 67 68 69 -5
60 61 62 63 64 -6
55 56 57 58 59 -7
50 51 52 53 54 -8
^
-8 -7 -6 -5 -4 -3 -2 -1 X=0 1 2 3 4 5 6 7 8 9 10
N is broken into digits of base=norm=r*r+1, ie. digits 0 to r*r inclusive. This makes
horizontal runs of r*r+1 many points, such as N=5 to N=9 etc above. In the default r=1
these runs are 2 long whereas for r=2 they're 2*2+1=5 long, or r=3 would be 3*3+1=10, etc.
The offset back for each run like N=5 shown is the r in i-r, then the next level is
(i-r)^2 = (-2r*i + r^2-1) so N=25 begins at Y=-2*2=-4, X=2*2-1=3.
The successive replications tile the plane for any r, though the N values needed to rotate
around and do so become large if norm=r*r+1 is large.
Fractal
The i-1 twindragon is usually conceived as taking fractional N like 0.abcde in binary and
giving fractional complex X+iY. The twindragon is then all the points of the complex
plane reached by such fractional N. This set of points can be shown to be connected and
to fill a certain radius around the origin.
The code here might be pressed into use for that to some finite number of bits by
multiplying up to make an integer N
Nint = Nfrac * 256^k
Xfrac = Xint / 16^k
Yfrac = Yint / 16^k
256 is a good power because b^8=16 is a positive real and so there's no rotations to apply
to the resulting X,Y, only a power-of-16 division (b^8)^k=16^k each. Using b^4=-4 for a
multiplier 16^k and divisor (-4)^k would be almost as easy too, requiring just sign
changes if k odd.
FUNCTIONS
See "FUNCTIONS" in Math::PlanePath for behaviour common to all path classes.
"$path = Math::PlanePath::ComplexMinus->new ()"
"$path = Math::PlanePath::ComplexMinus->new (realpart => $r)"
Create and return a new path object.
"($x,$y) = $path->n_to_xy ($n)"
Return the X,Y coordinates of point number $n on the path. Points begin at 0 and if
"$n < 0" then the return is an empty list.
$n should be an integer, it's unspecified yet what will be done for a fraction.
Level Methods
"($n_lo, $n_hi) = $path->level_to_n_range($level)"
Return "(0, 2**$level - 1)", or with "realpart" option return "(0, $norm**$level - 1)"
where norm=realpart^2+1.
FORMULAS
X,Y to N
A given X,Y representing X+Yi can be turned into digits of N by successive complex
divisions by i-r. Each digit of N is a real remainder 0 to r*r inclusive from that
division.
The base formula above is
X+Yi = a[n]*b^n + ... + a[2]*b^2 + a[1]*b + a[0]
and we want the a[0]=digit to be a real 0 to r*r inclusive. Subtracting a[0] and dividing
by b will give
(X+Yi - digit) / (i-r)
= - (X-digit + Y*i) * (i+r) / norm
= (Y - (X-digit)*r)/norm
+ i * - ((X-digit) + Y*r)/norm
which is
Xnew = Y - (X-digit)*r)/norm
Ynew = -((X-digit) + Y*r)/norm
The a[0] digit must make both Xnew and Ynew parts integers. The easiest one to calculate
from is the imaginary part, from which require
- ((X-digit) + Y*r) == 0 mod norm
so
digit = X + Y*r mod norm
This digit value makes the real part a multiple of norm too, as can be seen from
Xnew = Y - (X-digit)*r
= Y - X*r - (X+Y*r)*r
= Y - X*r - X*r + Y*r*r
= Y*(r*r+1)
= Y*norm
Notice Ynew is the quotient from (X+Y*r)/norm rounded downwards (towards negative
infinity). Ie. in the division "X+Y*r mod norm" which calculates the digit, the quotient
is Ynew and the remainder is the digit.
X Axis N for Realpart 1
For base i-1, Penney shows the N on the X axis are
X axis N in hexadecimal uses only digits 0, 1, C, D
= 0, 1, 12, 13, 16, 17, 28, 29, 192, 193, 204, 205, 208, ...
Those on the positive X axis have an odd number of digits and on the X negative axis an
even number of digits.
To be on the X axis the imaginary parts of the base powers b^k must cancel out to leave
just a real part. The powers repeat in an 8-long cycle
k b^k for b=i-1
0 +1
1 i -1
2 -2i +0 \ pair cancel
3 2i +2 /
4 -4
5 -4i +4
6 8i +0 \ pair cancel
7 -8i -8 /
The k=0 and k=4 bits are always reals and can always be included. Bits k=2 and k=3 have
imaginary parts -2i and 2i which cancel out, so they can be included together. Similarly
k=6 and k=7 with 8i and -8i. The two blocks k=0to3 and k=4to7 differ only in a negation
so the bits can be reckoned in groups of 4, which is hexadecimal. Bit 1 is digit value 1
and bits 2,3 together are digit value 0xC, so adding one or both of those gives
combinations are 0,1,0xC,0xD.
The high hex digit determines the sign, positive or negative, of the total real part.
Bits k=0 or k=2,3 are positive. Bits k=4 or k=6,7 are negative, so
N for X>0 N for X<0
0x01.. 0x1_.. even number of hex 0,1,C,D following
0x0C.. 0xC_.. "_" digit any of 0,1,C,D
0x0D.. 0xD_..
which is equivalent to X>0 is an odd number of hex digits or X<0 is an even number. For
example N=28=0x1C is at X=-2 since that N is X<0 form "0x1_".
The order of the values on the positive X axis is obtained by taking the digits in reverse
order on alternate positions
0,1,C,D high digit
D,C,1,0
0,1,C,D
...
D,C,1,0
0,1,C,D low digit
For example in the following notice the first and third digit increases, but the middle
digit decreases,
X=4to7 N=0x1D0,0x1D1,0x1DC,0x1DD
X=8to11 N=0x1C0,0x1C1,0x1CC,0x1CD
X=12to15 N=0x110,0x111,0x11C,0x11D
X=16to19 N=0x100,0x101,0x10C,0x10D
X=20to23 N=0xCD0,0xCD1,0xCDC,0xCDD
For the negative X axis it's the same if reading by increasing X, ie. upwards toward
+infinity, or the opposite way around if reading decreasing X, ie. more negative downwards
toward -infinity.
Y Axis N for Realpart 1
For base i-1 Penny also characterises the N values on the Y axis,
Y axis N in base-64 uses only
at even digits 0, 3, 4, 7, 48, 51, 52, 55
at odd digit 0, 1, 12, 13, 16, 17, 28, 29
= 0,3,4,7,48,51,52,55,64,67,68,71,112,115,116,119, ...
Base-64 means taking N in 6-bit blocks. Digit positions are counted starting from the
least significant digit as position 0 which is even. So the low digit can be only
0,3,4,etc, then the second digit only 0,1,12,etc, and so on.
This arises from (i-1)^6 = 8i which gives a repeating pattern of 6-bit blocks. The
different patterns at odd and even positions are since i^2 = -1.
Boundary Length
The length of the boundary of unit squares for the first norm^k many points, ie. N=0 to
N=norm^k-1 inclusive, is calculated in
William J. Gilbert, "The Fractal Dimension of Sets Derived From Complex Bases",
Canadian Mathematical Bulletin, volume 29, number 4, 1986.
<http://www.math.uwaterloo.ca/~wgilbert/Research/GilbertFracDim.pdf>
The boundary formula is a 3rd-order recurrence. For the twindragon case it is
for realpart=1
boundary[k] = boundary[k-1] + 2*boundary[k-3]
= 4, 6, 10, 18, 30, 50, 86, 146, 246, 418, 710, ...
4 + 2*x + 4*x^2
generating function ---------------
1 - x - 2*x^3
The first three boundaries are as follows. Then the recurrence gives the next boundary[3]
= 10+2*4 = 18.
k area boundary[k]
--- ---- -----------
+---+
0 2^k = 1 4 | 0 |
+---+
+---+---+
1 2^k = 2 6 | 0 1 |
+---+---+
+---+---+
| 2 3 |
2 2^k = 4 10 +---+ +---+
| 0 1 |
+---+---+
Gilbert calculates for any i-r by taking the boundary in three parts A,B,C and showing how
in the next replication level those boundary parts transform into multiple copies of the
preceding level parts. The replication is easier to visualize for a bigger "r" than for
the twindragon because in bigger r it's clearer how the A, B and C parts differ. The
length replications are
A -> A * (2*r-1) + C * 2*r
B -> A * (r^2-2*r+2) + C * (r-1)^2
C -> B
starting from
A = 2*r
B = 2
C = 2 - 2*r
total boundary = A+B+C
For the twindragon realpart=1 these A,B,C are already in the form of a recurrence
A->A+2*C, B->A, C->B, per the formula above. For other real parts a little matrix
rearrangement turns the A,B,C parts into recurrence
boundary[k] = boundary[k-1] * (2*r - 1)
+ boundary[k-2] * (norm - 2*r)
+ boundary[k-3] * norm
starting from
boundary[0] = 4 # single square cell
boundary[1] = 2*norm + 2 # oblong of norm many cells
boundary[2] = 2*(norm-1)*(r+2) + 4
For example
for realpart=2
boundary[k] = 3*boundary[k-1] + 1*boundary[k-2] + 5*boundary[k-3]
= 4, 12, 36, 140, 516, 1868, 6820, 24908, 90884, ...
4 - 4*x^2
generating function ---------------------
1 - 3*x - x^2 - 5*x^3
If calculating for large k values then the matrix form can be powered up rather than
repeated additions. (As usual for all such linear recurrences.)
OEIS
Entries in Sloane's Online Encyclopedia of Integer Sequences related to this path include
<http://oeis.org/A066321> (etc)
realpart=1 (twindragon, the default)
A066321 N on X axis, being the base i-1 positive reals
A066323 N on X axis, in binary
A066322 diffs (N at X=16k+4) - (N at X=16k+3)
A003476 boundary length / 2
recurrence a(n) = a(n-1) + 2*a(n-3)
A203175 boundary length, starting from 4
(believe its conjectured recurrence is true)
A052537 boundary length part A, B or C, per Gilbert's paper
SEE ALSO
Math::PlanePath, Math::PlanePath::DragonCurve, Math::PlanePath::ComplexPlus
HOME PAGE
<http://user42.tuxfamily.org/math-planepath/index.html>
LICENSE
Copyright 2011, 2012, 2013, 2014, 2015 Kevin Ryde
Math-PlanePath is free software; you can redistribute it and/or modify it under the terms
of the GNU General Public License as published by the Free Software Foundation; either
version 3, or (at your option) any later version.
Math-PlanePath is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with Math-
PlanePath. If not, see <http://www.gnu.org/licenses/>.