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 all the way around become big if norm=r*r+1 is big.

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.

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.

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.

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

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

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->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, ...
                                 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.)