Math-PlanePath
view release on metacpan or search on metacpan
lib/Math/PlanePath/ComplexMinus.pm view on Meta::CPAN
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 all the way around become big if norm=r*r+1 is big.
=head2 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.
=head1 FUNCTIONS
See L<Math::PlanePath/FUNCTIONS> for behaviour common to all path classes.
=over 4
=item C<$path = Math::PlanePath::ComplexMinus-E<gt>new ()>
=item C<$path = Math::PlanePath::ComplexMinus-E<gt>new (realpart =E<gt> $r)>
Create and return a new path object.
=item C<($x,$y) = $path-E<gt>n_to_xy ($n)>
Return the X,Y coordinates of point number C<$n> on the path. Points begin
at 0 and if C<$n E<lt> 0> then the return is an empty list.
C<$n> should be an integer, it's unspecified yet what will be done for a
fraction.
=back
=head2 Level Methods
=over
=item C<($n_lo, $n_hi) = $path-E<gt>level_to_n_range($level)>
Return C<(0, 2**$level - 1)>, or with C<realpart> option return C<(0,
$norm**$level - 1)> where norm=realpart^2+1.
=back
=head1 FORMULAS
Various formulas and pictures etc for the i-1 case can be found in the
author's long mathematical write-up (section "Complex Base i-1")
=over
L<http://user42.tuxfamily.org/dragon/index.html>
=back
=head2 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 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.
=cut
# Is this quite right ? ...
#
# =head2 Radius Range
#
# In general for base i-1 after the first few innermost levels each
# N=2^level increases the covered radius around by a factor sqrt(2), ie.
#
# N = 0 to 2^level-1
# Xmin,Ymin closest to origin
# Xmin^2+Ymin^2 approx 2^(level-7)
#
# The "level-7" is since the innermost few levels take a while to cover the
# points surrounding the origin. Notice for example X=1,Y=-1 is not reached
# until N=58. But after that it grows like N approx = pi*R^2.
=pod
=head2 X Axis N for Realpart 1
For base i-1, Penney shows the N on the X axis are
lib/Math/PlanePath/ComplexMinus.pm view on Meta::CPAN
which is equivalent to XE<gt>0 is an odd number of hex digits or XE<lt>0 is
an even number. For example N=28=0x1C is at X=-2 since that N is XE<lt>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.
=head2 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.
=head2 Boundary Length
X<Gilbert, William J.>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
=over
William J. Gilbert, "The Fractal Dimension of Sets Derived From Complex
Bases", Canadian Mathematical Bulletin, volume 29, number 4, 1986.
L<http://www.math.uwaterloo.ca/~wgilbert/Research/GilbertFracDim.pdf>
=back
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, ... (2*A003476)
4 + 2*x + 4*x^2
generating function ---------------
1 - x - 2*x^3
=cut
# GP-Test 2*4 + 10 == 18
# GP-Test 2*6 + 18 == 30
# GP-DEFINE gB1(x) = (4 + 2*x + 4*x^2) / (1 - x - 2*x^3)
# GP-Test Vec(gB1(x) - O(x^11)) == [4,6,10,18,30,50,86,146,246,418,710]
=pod
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 i-1 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 i-1 realpart=1 these A,B,C are already in the form of a recurrence
A-E<gt>A+2*C, B-E<gt>A, C-E<gt>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, ...
4 - 4*x^2
generating function ---------------------
1 - 3*x - x^2 - 5*x^3
=cut
# GP-DEFINE gB2(x) = (4 - 4*x^2) / (1 - 3*x - x^2 - 5*x^3)
# GP-Test Vec(gB2(x) - O(x^9)) == [4,12,36,140,516,1868,6820,24908,90884]
# GP-Test 5*4+1*12+3*36 == 140
# GP-Test 5*12+1*36+3*140 == 516
# not in OEIS: 4, 12, 36, 140, 516, 1868, 6820, 24908
=pod
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.)
=head1 OEIS
Entries in Sloane's Online Encyclopedia of Integer Sequences related to
this path include
=over
L<http://oeis.org/A318438> (etc)
=back
realpart=1 (base i-1, the default)
A318438 X coordinate
A318439 Y coordinate
A318479 norm X^2 + Y^2
A066321 N on X>=0 axis, or N/2 of North-West diagonal
A271472 same in binary
A066323 number of 1 bits
A256441 N on negative X axis, X<=0
A073791 X axis X sorted by N
A320283 Y axis Y sorted by N
A137426 dX/2 at N=2^(k+2)-1
dY at N=2^k-1 (step to next replication level)
A066322 diffs (N at X=16k+4) - (N at X=16k+3)
A340566 permutation N by diagonals +/-
( run in 0.399 second using v1.01-cache-2.11-cpan-df04353d9ac )