Project Euler

Heighway Dragon

The Heighway Dragon curve is defined by the Lindenmayer system: Axiom: F_a Rules: F_a -> F_a R F_b, F_b -> F_a L F_b Here F_a and F_b mean "draw forward", R means "turn right 90 degrees", and L mea...

Source sync May 21, 2026

Problem #0220

Level Level 09

Solved By 2,458

Languages C++, Python

Answer 139776,963904

Length 377 words

modular_arithmeticrecursionbrute_force

Let $D_0$ be the two-letter string "Fa". For $n \geq 1$, derive Dn from $D_{n - 1}$ by the string-rewriting rules:

"a" â†’ "aRbFR"
"b" â†’ "LFaLb"

Thus, $D_0 =$ "Fa", $D_1 = $ "FaRbFR", $D_2 = $ "FaRbFRRLFaLbFR", and so on.

These strings can be interpreted as instructions to a computer graphics program, with "F" meaning "draw forward one unit", "L" meaning "turn left $90$ degrees", "R" meaning "turn right 90 degrees", and "a" and "b" being ignored. The initial position of the computer cursor is $(0,0)$, pointing up towards $(0,1)$.

Then Dn is an exotic drawing known as the Heighway Dragon of order n. For example, D10 is shown below; counting each "F" as one step, the highlighted spot at $(18,16)$ is the position reached after $500$ steps.

What is the position of the cursor after $1012$ steps in $D_{50}$ ?

Give your answer in the form $x,y$ with no spaces.

Problem 220: Heighway Dragon

Mathematical Development

String Expansion

Let $F_a^{(n)}$ and $F_b^{(n)}$ denote the level- $n$ expansions of the two symbols. Their recursive structure is

F_a^{(n)} = F_a^{(n-1)} \; R \; F_b^{(n-1)}, \qquad F_b^{(n)} = F_a^{(n-1)} \; L \; F_b^{(n-1)}.

Each level doubles the number of forward moves, so both expansions contain exactly

2^n

draw steps.

Precomputed Displacements

For each level $n$ , each symbol $s \in \{a,b\}$ , and each starting direction $d \in \{0,1,2,3\}$ corresponding to north, east, south, west, define:

$\Delta x(n,s,d)$ and $\Delta y(n,s,d)$ : the net displacement after executing the whole expansion
$\operatorname{rot}(n,s)$ : the net change of direction modulo 4

At level 0, both symbols are just one forward step, so

\operatorname{rot}(0,a) = \operatorname{rot}(0,b) = 0

and the displacement is one unit in the current direction.

From the recursive rules,

\operatorname{rot}(n,a) = \operatorname{rot}(n-1,a) + 1 + \operatorname{rot}(n-1,b) \pmod 4,

\operatorname{rot}(n,b) = \operatorname{rot}(n-1,a) - 1 + \operatorname{rot}(n-1,b) \pmod 4.

The displacement recurrences follow by concatenating the displacement of the first half with the appropriately rotated displacement of the second half.

Recursive Query

To find the position after only the first $k$ drawing steps of $F_s^{(n)}$ , observe that the first half always consists of $F_a^{(n-1)}$ and contains exactly

2^{n-1}

steps.

If $k \leq 2^{n-1}$ , recurse into the first half.
If $k > 2^{n-1}$ , take the full displacement of the first half, apply the intervening turn ( $R$ for $a$ , $L$ for $b$ ), and recurse into the second half $F_b^{(n-1)}$ with the remaining step count.

Because each recursive call drops the level by 1, the query depth is at most 50.

Editorial

The full dragon string is astronomically large, but we never need to build it. The recursive production rules already tell us how a level- $n$ block is assembled from two level- $(n-1)$ blocks and one turn, so the useful information is just the net effect of a block: where it ends and which way it is facing.

Once those net effects are precomputed, answering the $10^{12}$ -step query becomes a binary decomposition problem. Each level splits into a first half of size $2^{n-1}$ and a second half of the same size, so we only need to decide which half contains the requested step count. The recursion walks down at most 50 levels and accumulates the displacement along the branches that are taken in full.

Pseudocode

For each direction, record the unit step vectors.

Initialize level-0 data:
    each symbol moves one step in the current direction
    and produces zero net rotation.

For n from 1 to 50:
    Compute rot(n, a) and rot(n, b) from the two recursive formulas.
    For each starting direction d:
        Build the displacement of F_a^(n) by:
            taking the full displacement of F_a^(n-1),
            turning right,
            then adding the displacement of F_b^(n-1) in the new direction.
        Build the displacement of F_b^(n) similarly, but turn left.

Define Query(n, symbol, direction, steps):
    If steps = 0:
        return (0, 0)
    If n = 0:
        return one unit step in the current direction

    half = 2^(n - 1)
    If steps <= half:
        return Query(n - 1, a, direction, steps)

    first = full displacement of F_a^(n-1) in this direction
    direction_after_first = direction rotated by rot(n - 1, a)
    direction_for_second = apply R or L depending on the symbol

    second = Query(n - 1, b, direction_for_second, steps - half)
    return first + second

Return Query(50, a, north, 10^12)

Complexity Analysis

Time: $O(50 \cdot 4)$ for the precomputation tables and $O(50)$ for the recursive query.
Space: $O(50 \cdot 4)$ for the displacement tables.

Answer

\boxed{139776,963904}

C++ project_euler/problem_220/solution.cpp

#include <bits/stdc++.h>
using namespace std;

/*
 * Problem 220: Heighway Dragon
 *
 * L-system: Fa -> Fa R Fb, Fb -> Fa L Fb
 * Start at origin facing up (north). Find position after 10^12 steps
 * at expansion level 50.
 *
 * Approach: Recursive decomposition with precomputed displacement tables.
 *
 * Answer: 139776,963904
 */

const int MAXN = 51;

// Directions: 0=N(+y), 1=E(+x), 2=S(-y), 3=W(-x)
int dx_dir[4] = {0, 1, 0, -1};
int dy_dir[4] = {1, 0, -1, 0};

// For each level n, symbol s (0=a, 1=b), direction d:
// disp_x[n][s][d], disp_y[n][s][d] = net displacement
// rot[n][s] = net rotation (mod 4)

long long disp_x[MAXN][2][4];
long long disp_y[MAXN][2][4];
int rot[MAXN][2];

void precompute() {
    // Base case: level 0
    for (int s = 0; s < 2; s++) {
        rot[0][s] = 0;
        for (int d = 0; d < 4; d++) {
            disp_x[0][s][d] = dx_dir[d];
            disp_y[0][s][d] = dy_dir[d];
        }
    }

    // Recursive case
    for (int n = 1; n < MAXN; n++) {
        // Fa^(n) = Fa^(n-1) R Fb^(n-1)
        // rot(n, a) = rot(n-1, a) + 1 + rot(n-1, b)  (mod 4)
        rot[n][0] = (rot[n-1][0] + 1 + rot[n-1][1]) % 4;

        // Fb^(n) = Fa^(n-1) L Fb^(n-1)
        // rot(n, b) = rot(n-1, a) - 1 + rot(n-1, b)  (mod 4)
        rot[n][1] = (rot[n-1][0] - 1 + rot[n-1][1] + 4) % 4;

        for (int d = 0; d < 4; d++) {
            // Fa^(n) starting in direction d:
            // First half: Fa^(n-1) in direction d
            long long x1 = disp_x[n-1][0][d];
            long long y1 = disp_y[n-1][0][d];
            int d_after_first_a = (d + rot[n-1][0]) % 4;
            // Turn right
            int d_second_a = (d_after_first_a + 1) % 4;
            // Second half: Fb^(n-1) in direction d_second_a
            long long x2 = disp_x[n-1][1][d_second_a];
            long long y2 = disp_y[n-1][1][d_second_a];
            disp_x[n][0][d] = x1 + x2;
            disp_y[n][0][d] = y1 + y2;

            // Fb^(n) starting in direction d:
            // First half: Fa^(n-1) in direction d
            x1 = disp_x[n-1][0][d];
            y1 = disp_y[n-1][0][d];
            int d_after_first_b = (d + rot[n-1][0]) % 4;
            // Turn left
            int d_second_b = (d_after_first_b + 3) % 4; // +3 = -1 mod 4
            // Second half: Fb^(n-1) in direction d_second_b
            x2 = disp_x[n-1][1][d_second_b];
            y2 = disp_y[n-1][1][d_second_b];
            disp_x[n][1][d] = x1 + x2;
            disp_y[n][1][d] = y1 + y2;
        }
    }
}

// Find position after k steps of symbol s at level n, starting at direction d
// Returns (x, y)
pair<long long, long long> query(int n, int s, int d, long long k) {
    if (k == 0) return {0, 0};

    if (n == 0) {
        // Single step
        return {dx_dir[d], dy_dir[d]};
    }

    long long half = 1LL << (n - 1);

    // First half is always Fa^(n-1)
    // Turn depends on symbol: R for a, L for b
    // Second half is always Fb^(n-1)

    if (k <= half) {
        // Only in first half (Fa^(n-1))
        return query(n - 1, 0, d, k);
    } else {
        // Full first half + turn + partial second half
        long long x1 = disp_x[n-1][0][d];
        long long y1 = disp_y[n-1][0][d];
        int d_after = (d + rot[n-1][0]) % 4;
        int turn = (s == 0) ? 1 : 3; // R=+1 for a, L=-1=+3 for b
        int d_second = (d_after + turn) % 4;

        auto [x2, y2] = query(n - 1, 1, d_second, k - half);
        return {x1 + x2, y1 + y2};
    }
}

int main() {
    precompute();

    long long steps = 1000000000000LL; // 10^12
    int level = 50;

    // Start at origin, facing north (direction 0), symbol Fa at level 50
    auto [x, y] = query(level, 0, 0, steps);

    cout << x << "," << y << endl;
    return 0;
}

Python project_euler/problem_220/solution.py

"""
Problem 220: Heighway Dragon

L-system: Fa -> Fa R Fb, Fb -> Fa L Fb
Start at origin (0,0) facing north (up).
Find position after 10^12 drawing steps at expansion level 50.

Approach: Precompute net displacement and rotation for each level/symbol/direction,
then recursively decompose the query.

Answer: 139776,963904
"""

import sys
sys.setrecursionlimit(100)

MAXN = 51

# Directions: 0=N(+y), 1=E(+x), 2=S(-y), 3=W(-x)
DX = [0, 1, 0, -1]
DY = [1, 0, -1, 0]

def solve():
    # Precompute displacement and rotation tables
    # disp_x[n][s][d], disp_y[n][s][d], rot[n][s]
    disp_x = [[[0]*4 for _ in range(2)] for _ in range(MAXN)]
    disp_y = [[[0]*4 for _ in range(2)] for _ in range(MAXN)]
    rot = [[0]*2 for _ in range(MAXN)]

    # Base case: level 0
    for s in range(2):
        rot[0][s] = 0
        for d in range(4):
            disp_x[0][s][d] = DX[d]
            disp_y[0][s][d] = DY[d]

    # Build up
    for n in range(1, MAXN):
        # Fa^(n) = Fa^(n-1) R Fb^(n-1)
        rot[n][0] = (rot[n-1][0] + 1 + rot[n-1][1]) % 4
        # Fb^(n) = Fa^(n-1) L Fb^(n-1)
        rot[n][1] = (rot[n-1][0] - 1 + rot[n-1][1]) % 4

        for d in range(4):
            # Fa^(n)
            x1 = disp_x[n-1][0][d]
            y1 = disp_y[n-1][0][d]
            d_after = (d + rot[n-1][0]) % 4
            d_second = (d_after + 1) % 4  # R = +1
            x2 = disp_x[n-1][1][d_second]
            y2 = disp_y[n-1][1][d_second]
            disp_x[n][0][d] = x1 + x2
            disp_y[n][0][d] = y1 + y2

            # Fb^(n)
            x1 = disp_x[n-1][0][d]
            y1 = disp_y[n-1][0][d]
            d_after = (d + rot[n-1][0]) % 4
            d_second = (d_after - 1) % 4  # L = -1
            x2 = disp_x[n-1][1][d_second]
            y2 = disp_y[n-1][1][d_second]
            disp_x[n][1][d] = x1 + x2
            disp_y[n][1][d] = y1 + y2

    def query(n, s, d, k):
        """Position after k steps of symbol s at level n, starting direction d."""
        if k == 0:
            return (0, 0)
        if n == 0:
            return (DX[d], DY[d])

        half = 1 << (n - 1)

        if k <= half:
            return query(n - 1, 0, d, k)
        else:
            x1 = disp_x[n-1][0][d]
            y1 = disp_y[n-1][0][d]
            d_after = (d + rot[n-1][0]) % 4
            turn = 1 if s == 0 else -1  # R for a, L for b
            d_second = (d_after + turn) % 4
            x2, y2 = query(n - 1, 1, d_second, k - half)
            return (x1 + x2, y1 + y2)

    steps = 10**12
    level = 50
    x, y = query(level, 0, 0, steps)
    print(f"{x},{y}")

if __name__ == "__main__":
    solve()