Problem 220: Heighway Dragon

Project Euler Problem 220 Solution

EulerSolve provides an optimized solution for Project Euler Problem 220, Heighway Dragon, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The Heighway Dragon in this problem is generated from the symbolic rule set $D_0=Fa$, together with the substitutions $a\to aRbFR$ and $b\to LFaLb$. In turtle-graphics terms, $F$ means “move forward one unit”, $L$ and $R$ mean quarter-turns, and the letters $a$ and $b$ are only structural markers: they do not move the turtle. After $n$ substitution rounds, deleting $a$ and $b$ leaves the instruction word $D_n$. The task is to find the turtle's position after $10^{12}$ forward moves in $D_{50}$. The statement also provides a checkpoint: in $D_{10}$, after 500 forward moves, the position is $(18,16)$. A direct expansion is impossible because $D_{50}$ contains $2^{50}$ forward moves, about $1.13\times 10^{15}$, so the solution has to work with compressed recursive blocks rather than the full string. Mathematical Approach The right object to precompute is not the text of a block, but its geometric effect. Every recursive block can be replaced by a small summary that records displacement, heading change, and the number of forward moves it contains. Once those summaries are known for all depths, the required position is a prefix query on a recursive grammar tree....

Detailed mathematical approach

Problem Summary

The Heighway Dragon in this problem is generated from the symbolic rule set $D_0=Fa$, together with the substitutions $a\to aRbFR$ and $b\to LFaLb$. In turtle-graphics terms, $F$ means “move forward one unit”, $L$ and $R$ mean quarter-turns, and the letters $a$ and $b$ are only structural markers: they do not move the turtle. After $n$ substitution rounds, deleting $a$ and $b$ leaves the instruction word $D_n$.

The task is to find the turtle's position after $10^{12}$ forward moves in $D_{50}$. The statement also provides a checkpoint: in $D_{10}$, after 500 forward moves, the position is $(18,16)$. A direct expansion is impossible because $D_{50}$ contains $2^{50}$ forward moves, about $1.13\times 10^{15}$, so the solution has to work with compressed recursive blocks rather than the full string.

Mathematical Approach

The right object to precompute is not the text of a block, but its geometric effect. Every recursive block can be replaced by a small summary that records displacement, heading change, and the number of forward moves it contains. Once those summaries are known for all depths, the required position is a prefix query on a recursive grammar tree.

Rigid-Motion Summaries for Instruction Blocks

For any instruction block $W$, define its summary by

$$\Sigma(W)=(\vec v(W),d(W),s(W)),$$

where $\vec v(W)=(\Delta x,\Delta y)$ is the net displacement, $d(W)\in\{0,1,2,3\}$ is the net number of right turns modulo 4, and $s(W)$ is the number of forward commands $F$ inside $W$.

If $\operatorname{rot}(x,y)=(y,-x)$ denotes a $90^\circ$ right rotation, then concatenation of two blocks is described by

$$\Sigma_1\star\Sigma_2=\left(\vec v_1+\operatorname{rot}^{d_1}(\vec v_2),\ (d_1+d_2)\bmod 4,\ s_1+s_2\right),$$

for $\Sigma_i=(\vec v_i,d_i,s_i)$. This is exactly the turtle rule: the second block must be interpreted in the heading left behind by the first block.

The atomic commands therefore have summaries

$$\Sigma(F)=((1,0),0,1),\qquad \Sigma(R)=((0,0),1,0),\qquad \Sigma(L)=((0,0),3,0).$$

These vectors are written in the local frame of the block itself. When a stored block is later applied to the actual walk, its displacement is rotated by the current heading first and then added to the global position.

Recursive Macro Families Behind the Dragon

For clarity, let $A_n$ and $B_n$ denote the fully expanded instruction blocks produced by the two silent symbols after $n$ remaining substitution levels, with empty base case

$$A_0=B_0=\varepsilon.$$

The substitution rules become the recurrences

$$A_{n+1}=A_n\,R\,B_n\,F\,R,\qquad B_{n+1}=L\,F\,A_n\,L\,B_n.$$

The whole dragon word is then

$$D_n=F\,A_n.$$

Applying the composition law from the previous subsection gives

$$\Sigma(A_{n+1})=\Sigma(A_n)\star\Sigma(R)\star\Sigma(B_n)\star\Sigma(F)\star\Sigma(R),$$

$$\Sigma(B_{n+1})=\Sigma(L)\star\Sigma(F)\star\Sigma(A_n)\star\Sigma(L)\star\Sigma(B_n).$$

So every summary at depth $n+1$ is built from a constant number of summaries at depth $n$. This is the mathematical reason the preprocessing is linear in the order rather than exponential in the expanded word length.

Forward-Step Counts and a Useful Heading Invariant

Let $f_n=s(A_n)=s(B_n)$. Both recursive definitions contain one explicit $F$ and one copy each of the previous $A_n$ and $B_n$, so

$$f_{n+1}=2f_n+1,\qquad f_0=0.$$

Solving this recurrence gives

$$f_n=2^n-1.$$

Therefore $D_n=F A_n$ contains exactly

$$1+f_n=2^n$$

forward moves. In particular, $D_{50}$ has $2^{50}$ forward steps, which is why explicit expansion is hopeless.

A second invariant concerns the final heading. If $t_n=d(A_n)$ and $u_n=d(B_n)$, then

$$t_{n+1}\equiv t_n+1+u_n+1 \pmod 4,\qquad u_{n+1}\equiv 3+t_n+3+u_n \pmod 4,$$

with $t_0=u_0=0$. From this it follows that

$$t_n=u_n=2\qquad (n\ge 1).$$

So every non-empty $A_n$ or $B_n$ turns the turtle around by $180^\circ$. The implementations do not need a closed-form displacement formula, but this heading invariant is useful for understanding why the summaries compose so cleanly.

Prefix Evaluation Without Expanding the Whole Word

Suppose we only want the first $K$ forward moves inside $A_n$. If $K\ge s(A_n)$, then the entire block lies inside the requested prefix, so one application of $\Sigma(A_n)$ is enough. If $K<s(A_n)$, then the desired prefix lies somewhere inside

$$A_n=A_{n-1}\,R\,B_{n-1}\,F\,R,$$

and we process those pieces from left to right. The crucial detail is that the budget counts only forward moves. Turns are still applied when they occur before the $K$-th forward move, but they do not decrease the budget.

The same reasoning applies to $B_n$. An induction on $n$ shows that this recursive descent returns exactly the turtle state after the requested prefix. Structurally, at each depth there is at most one child block whose contribution is only partially used; everything before it can be swallowed whole from a precomputed summary, and everything after it is irrelevant once the budget reaches zero. That single-partial-branch invariant is what keeps the runtime linear in the order.

Worked Example: The First Two Nontrivial Levels

The first non-empty blocks are

$$A_1=RFR,\qquad B_1=LFL.$$

Using the local east-facing frame, their summaries are

$$\Sigma(A_1)=((0,-1),2,1),\qquad \Sigma(B_1)=((0,1),2,1).$$

At the next level,

$$A_2=A_1\,R\,B_1\,F\,R,$$

$$\Sigma(A_2)=\Sigma(A_1)\star\Sigma(R)\star\Sigma(B_1)\star\Sigma(F)\star\Sigma(R)=((-1,-2),2,3).$$

This already shows the compressed viewpoint at work: $A_2$ contains three forward moves, but once its summary has been stored, those three moves can be consumed in one constant-time composition whenever the remaining budget covers them.

Since $D_2=F A_2$, the whole word has $2^2=4$ forward moves. Starting from the origin facing north, the initial $F$ reaches $(0,1)$. If the remaining budget is 3, the entire $A_2$ block can be applied at once. If the remaining budget were only 2, the algorithm would descend into $A_1$, pass the explicit turn, and recurse only as far into $B_1$ as the budget still allows. The actual $10^{12}$-step query uses exactly the same idea, just fifty levels deeper.

How the Code Works

Bottom-Up Summary Tables

The C++, Python, and Java implementations allocate one stored summary for each depth of the $A$-family and the $B$-family. The deepest level is the empty block, and the tables are filled upward using the two recurrence formulas above. Each entry stores only four small pieces of data: displacement, net quarter-turn, and forward-step count.

Because the problem asks for $D_{50}$, only 51 levels are needed. That makes the preprocessing tiny: the implementations never build the enormous string itself, only the summary tables for the recursive macros.

Top-Down Prefix Traversal

Execution starts at $(0,0)$ facing north, applies the initial $F$, and then traverses the $A$-family recursively. Whenever the remaining forward-step budget covers an entire stored macro, the implementation subtracts that macro's step count and composes its summary into the current state in $O(1)$.

If the budget does not cover the full macro, the implementation expands only that one macro level: it visits the child blocks in grammar order, applies turns immediately, and decreases the remaining budget only when a forward move is actually taken. The maintained invariant is simple and exact: after each recursive return, the stored global state equals the turtle state after the already processed prefix of forward moves.

All three languages implement the same mathematics. They also agree with the checkpoint from the statement: after 500 forward moves in $D_{10}$, the position is $(18,16)$. That sanity check confirms both the summary-composition law and the prefix-recursion logic.

Complexity Analysis

If the dragon order is $n$, building all summaries for $A_0,\dots,A_n$ and $B_0,\dots,B_n$ costs $O(n)$ time and $O(n)$ memory. Each new level is obtained from a constant number of already known summaries.

Computing the position after one target prefix also costs $O(n)$ time. The recursion may descend through all depths, but only one branch can remain partial at each level; every fully covered sibling is handled by a constant-time summary application. For Problem 220 this means the $10^{12}$-step query on $D_{50}$ is solved with a few dozen summary compositions instead of expanding a word with $2^{50}$ forward moves.

Footnotes and References

Project Euler Problem 220: https://projecteuler.net/problem=220
Dragon curve: Wikipedia - Dragon curve
L-system: Wikipedia - L-system
Turtle graphics: Wikipedia - Turtle graphics
Rotation matrix: Wikipedia - Rotation matrix

Problem 220 source code

C++

#include <array>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

using i64 = long long;
using u64 = std::uint64_t;

struct Segment {
    i64 x = 0;
    i64 y = 0;
    int dir = 0;
    u64 steps = 0;
};

struct Options {
    int order = 50;
    u64 steps = 1000000000000ULL;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    int parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_u64_after_prefix(const std::string& arg, const std::string& prefix, u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0;
    for (char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<u64>(c - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);

        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
            continue;
        }
        if (parse_int_after_prefix(arg, "--order=", options.order) ||
            parse_u64_after_prefix(arg, "--steps=", options.steps)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return options.order >= 0;
}

Segment combine(Segment lhs, Segment rhs) {
    for (int i = 0; i < lhs.dir; ++i) {
        const i64 old_x = rhs.x;
        rhs.x = rhs.y;
        rhs.y = -old_x;
    }

    lhs.x += rhs.x;
    lhs.y += rhs.y;
    lhs.dir = (lhs.dir + rhs.dir) & 3;
    lhs.steps += rhs.steps;
    return lhs;
}

struct Walker {
    std::vector<Segment> a;
    std::vector<Segment> b;
    Segment state;
    u64 remaining = 0;

    static constexpr Segment kF = {1, 0, 0, 1};
    static constexpr Segment kL = {0, 0, 3, 0};
    static constexpr Segment kR = {0, 0, 1, 0};

    explicit Walker(const int order, const u64 steps)
        : a(static_cast<std::size_t>(order + 1)), b(static_cast<std::size_t>(order + 1)), remaining(steps) {
        for (int i = order - 1; i >= 0; --i) {
            a[static_cast<std::size_t>(i)] = combine(
                combine(combine(combine(a[static_cast<std::size_t>(i + 1)], kR), b[static_cast<std::size_t>(i + 1)]), kF),
                kR);
            b[static_cast<std::size_t>(i)] = combine(
                combine(combine(combine(kL, kF), a[static_cast<std::size_t>(i + 1)]), kL),
                b[static_cast<std::size_t>(i + 1)]);
        }

        state = {0, 0, 3, 0};  // Facing North in this coordinate system.
    }

    void apply_segment(const Segment& seg) {
        state = combine(state, seg);
    }

    void execute_a(const int level) {
        if (remaining == 0 || level >= static_cast<int>(a.size()) - 1) {
            return;
        }

        const Segment& full = a[static_cast<std::size_t>(level)];
        if (remaining >= full.steps) {
            remaining -= full.steps;
            apply_segment(full);
            return;
        }

        execute_a(level + 1);
        if (remaining == 0) {
            return;
        }
        apply_segment(kR);

        execute_b(level + 1);
        if (remaining == 0) {
            return;
        }

        --remaining;
        apply_segment(kF);
        if (remaining == 0) {
            return;
        }

        apply_segment(kR);
    }

    void execute_b(const int level) {
        if (remaining == 0 || level >= static_cast<int>(b.size()) - 1) {
            return;
        }

        const Segment& full = b[static_cast<std::size_t>(level)];
        if (remaining >= full.steps) {
            remaining -= full.steps;
            apply_segment(full);
            return;
        }

        apply_segment(kL);

        if (remaining == 0) {
            return;
        }
        --remaining;
        apply_segment(kF);
        if (remaining == 0) {
            return;
        }

        execute_a(level + 1);
        if (remaining == 0) {
            return;
        }

        apply_segment(kL);
        execute_b(level + 1);
    }

    std::pair<i64, i64> run() {
        if (remaining > 0) {
            --remaining;
            apply_segment(kF);
        }
        execute_a(0);
        return {state.x, state.y};
    }
};

std::pair<i64, i64> solve(const int order, const u64 steps) {
    Walker walker(order, steps);
    return walker.run();
}

bool run_checkpoints() {
    const auto sample = solve(10, 500);
    if (sample.first != 18 || sample.second != 16) {
        std::cerr << "Checkpoint failed for D10 after 500 steps" << '\n';
        return false;
    }
    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) {
        return 1;
    }
    if (options.run_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const auto answer = solve(options.order, options.steps);
    std::cout << answer.first << ',' << answer.second << '\n';
    return 0;
}

Python

# Problem 220: Heighway Dragon
# D-string D50 after 10^12 steps of executing F instructions.
# Use segment precomputation: for each level, precompute the net effect
# (displacement + rotation) of executing all F's in a(level) and b(level).

def solve():
    order = 50
    steps = 10**12

    # Segment = (x, y, dir, steps) where dir is rotation (0-3)
    # combine(lhs, rhs): apply rhs after lhs
    def combine(lhs, rhs):
        x, y, d, s = rhs
        # Rotate rhs displacement by lhs.dir
        for _ in range(lhs[2]):
            x, y = y, -x
        return (lhs[0] + x, lhs[1] + y, (lhs[2] + d) & 3, lhs[3] + s)

    F = (1, 0, 0, 1)
    L = (0, 0, 3, 0)  # turn left = -90 degrees
    R = (0, 0, 1, 0)  # turn right = +90 degrees

    # Build a[i] and b[i] segments
    # a0 = aRbFR, b0 = LFaLb (using a1,b1 etc)
    # Base case: a[order] = b[order] = (0,0,0,0) (empty)
    a = [(0, 0, 0, 0)] * (order + 1)
    b = [(0, 0, 0, 0)] * (order + 1)

    for i in range(order - 1, -1, -1):
        a[i] = combine(combine(combine(combine(a[i+1], R), b[i+1]), F), R)
        b[i] = combine(combine(combine(combine(L, F), a[i+1]), L), b[i+1])

    # Walker state
    state = (0, 0, 3, 0)  # facing north
    remaining = steps

    def apply_seg(seg):
        nonlocal state
        state = combine(state, seg)

    def execute_a(level):
        nonlocal remaining
        if remaining == 0 or level >= order:
            return
        full = a[level]
        if remaining >= full[3]:
            remaining -= full[3]
            apply_seg(full)
            return
        execute_a(level + 1)
        if remaining == 0: return
        apply_seg(R)
        execute_b(level + 1)
        if remaining == 0: return
        remaining -= 1
        apply_seg(F)
        if remaining == 0: return
        apply_seg(R)

    def execute_b(level):
        nonlocal remaining
        if remaining == 0 or level >= order:
            return
        full = b[level]
        if remaining >= full[3]:
            remaining -= full[3]
            apply_seg(full)
            return
        apply_seg(L)
        if remaining == 0: return
        remaining -= 1
        apply_seg(F)
        if remaining == 0: return
        execute_a(level + 1)
        if remaining == 0: return
        apply_seg(L)
        execute_b(level + 1)

    # First F step
    if remaining > 0:
        remaining -= 1
        apply_seg(F)
    execute_a(0)

    print(f"{state[0]},{state[1]}")

import sys
sys.setrecursionlimit(200)
solve()

Java

public class Euler220 {
    static long sx, sy, remaining;
    static int sd;

    public static void main(String[] args) {
        int order = 50;
        long steps = 1000000000000L;
        long[] F = { 1, 0, 0, 1 }, L = { 0, 0, 3, 0 }, R = { 0, 0, 1, 0 };
        long[][] a = new long[order + 1][4], b = new long[order + 1][4];
        for (int i = order - 1; i >= 0; i--) {
            a[i] = combine(combine(combine(combine(a[i + 1], R), b[i + 1]), F), R);
            b[i] = combine(combine(combine(combine(L, F), a[i + 1]), L), b[i + 1]);
        }
        sx = 0;
        sy = 0;
        sd = 3;
        remaining = steps;
        if (remaining > 0) {
            remaining--;
            apply(F);
        }
        exeA(a, b, F, L, R, 0, order);
        System.out.println(sx + "," + sy);
    }

    static void exeA(long[][] a, long[][] b, long[] F, long[] L, long[] R, int level, int order) {
        if (remaining == 0 || level >= order)
            return;
        if (remaining >= a[level][3]) {
            remaining -= a[level][3];
            apply(a[level]);
            return;
        }
        exeA(a, b, F, L, R, level + 1, order);
        if (remaining == 0)
            return;
        apply(R);
        exeB(a, b, F, L, R, level + 1, order);
        if (remaining == 0)
            return;
        remaining--;
        apply(F);
        if (remaining == 0)
            return;
        apply(R);
    }

    static void exeB(long[][] a, long[][] b, long[] F, long[] L, long[] R, int level, int order) {
        if (remaining == 0 || level >= order)
            return;
        if (remaining >= b[level][3]) {
            remaining -= b[level][3];
            apply(b[level]);
            return;
        }
        apply(L);
        if (remaining == 0)
            return;
        remaining--;
        apply(F);
        if (remaining == 0)
            return;
        exeA(a, b, F, L, R, level + 1, order);
        if (remaining == 0)
            return;
        apply(L);
        exeB(a, b, F, L, R, level + 1, order);
    }

    static void apply(long[] seg) {
        long x = seg[0], y = seg[1];
        for (int i = 0; i < sd; i++) {
            long tmp = x;
            x = y;
            y = -tmp;
        }
        sx += x;
        sy += y;
        sd = (int) ((sd + seg[2]) & 3);
    }

    static long[] combine(long[] lhs, long[] rhs) {
        long x = rhs[0], y = rhs[1];
        for (int i = 0; i < (int) lhs[2]; i++) {
            long tmp = x;
            x = y;
            y = -tmp;
        }
        return new long[] { lhs[0] + x, lhs[1] + y, (lhs[2] + rhs[2]) & 3, lhs[3] + rhs[3] };
    }
}