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

Start with $D_0=\mathtt{Fa}$ and rewrite it again and again: in every round each $\mathtt a$ is replaced by $\mathtt{aRbFR}$ and, at the same time, each $\mathtt b$ by $\mathtt{LFaLb}$, while $\mathtt F$, $\mathtt L$ and $\mathtt R$ stay as they are. After $n$ rounds we have $D_n$; thus $D_1=\mathtt{FaRbFR}$ and $D_2=\mathtt{FaRbFRRLFaLbFR}$. The word drives a cursor that starts at $(0,0)$ pointing towards $(0,1)$: $\mathtt F$ moves it forward one unit, $\mathtt L$ and $\mathtt R$ turn it $90^\circ$ to the left or right, and $\mathtt a$ and $\mathtt b$ are ignored; the drawing of $D_n$ is the Heighway dragon of order $n$. Counting each $\mathtt F$ as one step, the cursor is at $(18,16)$ after $500$ steps of $D_{10}$, and we must find its position after $10^{12}$ steps of $D_{50}$, in the form x,y with no spaces. ^[1]

The whole solution in one sentence: each block of instructions that grows out of one letter $\mathtt a$ or $\mathtt b$ is replaced by its net effect on the cursor (a displacement, a turn and a number of steps), which the two rules give for all orders, and the $10^{12}$ steps are then walked by jumping over whole blocks.

Mathematical Approach

We draw the first dragons by hand, cut $D_n$ into blocks that grow out of single letters, describe each block by its net effect and walk through just enough blocks.

1. The first dragons

The letters $\mathtt a$ and $\mathtt b$ never move the cursor, but they drive the rewriting, so we keep them through all rounds and delete them only for drawing; from now on $D_n$ means this word of commands. Thus $D_0$, $D_1$ and $D_2$ become $\mathtt F$, $\mathtt{FRFR}$ and $\mathtt{FRFRRLFLFR}$. Let us follow $D_2$ from $(0,0)$, facing north. $\mathtt F$ leads to $(0,1)$; $\mathtt R$ turns the cursor east and $\mathtt F$ leads to $(1,1)$; $\mathtt{RRL}$ is one right turn in total, so the cursor faces south and $\mathtt F$ leads to $(1,0)$; $\mathtt L$ turns it east again, $\mathtt F$ ends at $(2,0)$, and the last $\mathtt R$ only turns it. In the same way $D_3$, which has $22$ instructions, makes $8$ steps through $(0,1)$, $(1,1)$, $(1,0)$, $(2,0)$, $(2,-1)$, $(1,-1)$, $(1,-2)$ and $(2,-2)$.

Three things stand out: the number of steps doubles from one order to the next; every dragon begins with the whole previous one; and each of the end points $(0,1)$, $(1,1)$, $(2,0)$, $(2,-2)$ lies $\sqrt2$ times as far from the origin as the one before, turned by $45^\circ$ clockwise. Sections 2, 3 and 11 explain all three.

2. Blocks that grow out of one letter

Rewriting acts on every letter separately, as in an L-system, where all letters are rewritten at once and letters without a rule, the constants, stay unchanged. ^[3] Let $A_k$ and $B_k$ be the words into which $\mathtt a$ and $\mathtt b$ grow in $k$ rounds, with $\mathtt a$ and $\mathtt b$ deleted; $A_0$ and $B_0$ are empty. In $k+1$ rounds, $\mathtt a$ first becomes $\mathtt{aRbFR}$, and then each of these five letters grows for $k$ more rounds, and similarly for $\mathtt b$. Hence

$$\begin{aligned}A_{k+1}&=A_k\,\mathtt R\,B_k\,\mathtt F\,\mathtt R,\\ B_{k+1}&=\mathtt L\,\mathtt F\,A_k\,\mathtt L\,B_k.\end{aligned}$$

For example $A_1=\mathtt{RFR}$, $B_1=\mathtt{LFL}$ and $A_2=\mathtt{RFRRLFLFR}$. As $D_n$ is $\mathtt F$ followed by what $\mathtt a$ grows into, $D_n=\mathtt FA_n$, which for $n=2$ is the word of Section 1. Writing $A_n$ by its rule gives $D_n=D_{n-1}\,\mathtt R\,B_{n-1}\,\mathtt{FR}$, so every dragon begins with the previous one, and hence with all earlier ones: the position after $N$ steps is the same in every $D_n$ that has at least $N$ steps.

3. Counting steps

Let $s_k$ be the number of letters $\mathtt F$ in $A_k$. Each rule consists of one block $A_k$, one block $B_k$, one $\mathtt F$ and two turns, so $A_k$ and $B_k$ have equally many steps, and

$$s_0=0,\qquad s_{k+1}=2s_k+1.$$

By induction $s_k=2^k-1$, since $2(2^k-1)+1=2^{k+1}-1$; in the same way a block of order $k$ has $2(2^k-1)$ turns. So $D_n=\mathtt FA_n$ makes exactly $2^n$ steps, and $D_{50}$ makes $2^{50}=1\,125\,899\,906\,842\,624$ steps. Since $2^{40}=1\,099\,511\,627\,776$ already exceeds $10^{12}$, Section 2 even allows us to replace $D_{50}$ by $D_{40}$. Executing the steps one by one would still mean about $3\cdot10^{12}$ instructions.

4. Position and heading

Let $\rho(x,y)=(y,-x)$ be the clockwise quarter turn of the plane; its matrix is that of the anticlockwise rotation by $270^\circ$. ^[5] We record the heading of the cursor as a number $h\in\{0,1,2,3\}$, meaning that the cursor faces $\rho^h(1,0)$: east, south, west or north. Then $\mathtt F$ adds $\rho^h(1,0)$ to the position, $\mathtt R$ adds $1$ to $h$ and $\mathtt L$ adds $3$, always modulo $4$. The cursor of the problem starts at $(0,0)$ with $h=3$. As in turtle graphics, every instruction acts relative to the current position and heading. ^[4]

5. The net effect of a word

For a word $W$ of instructions let $\Sigma(W)=(v,d,s)$: a cursor that starts at $(0,0)$ facing east and executes $W$ ends at $v$ with heading $d$, and $W$ contains $s$ steps. Thus $\Sigma(\mathtt F)=((1,0),0,1)$, $\Sigma(\mathtt R)=((0,0),1,0)$ and $\Sigma(\mathtt L)=((0,0),3,0)$; we call these triples $f$, $r$ and $l$.

A cursor at $p$ with heading $h$ that executes $W$ ends at $p+\rho^h(v)$ with heading $h+d$, after $s$ more steps. Indeed, the map $q\mapsto p+\rho^h(q)$, a rotation by $h$ quarter turns followed by a shift by $p$, takes the origin facing east to $p$ facing $\rho^h(1,0)$; as the instructions act relative to the cursor, it maps the run of $W$ from the origin onto the run of $W$ from $p$. For $W_2$ executed after $W_1$, with $(p,h)=(v_1,d_1)$, this gives the rule for joining two words:

$$\begin{aligned}\Sigma(W_1W_2)=\bigl(&v_1+\rho^{d_1}(v_2),\\ &d_1+d_2,\,s_1+s_2\bigr).\end{aligned}$$

We write the right-hand side as $\Sigma(W_1)\star\Sigma(W_2)$, with headings added modulo $4$. Both sides describe one run of the joined word, so a longer word may be cut into pieces in any way. The state of the cursor (position, heading and steps so far) is itself such a triple, and executing $W$ turns the state $S$ into $S\star\Sigma(W)$. The start of the problem is $((0,0),3,0)=l$: facing north is facing east after one left turn.

6. The effects of all blocks

Since $\Sigma$ turns joining into $\star$, the rules of Section 2 become

$$\begin{aligned}\Sigma(A_{k+1})&=\Sigma(A_k)\star r\star\Sigma(B_k)\\ &\quad{}\star f\star r,\\ \Sigma(B_{k+1})&=l\star f\star\Sigma(A_k)\\ &\quad{}\star l\star\Sigma(B_k),\end{aligned}$$

starting from $\Sigma(A_0)=\Sigma(B_0)=((0,0),0,0)$, so $2n$ evaluations give the effects of all blocks up to order $n$. By hand: in $A_1=\mathtt{RFR}$ the cursor turns south, steps to $(0,-1)$ and turns west, so $\Sigma(A_1)=((0,-1),2,1)$, and likewise $\Sigma(B_1)=((0,1),2,1)$. Then $\Sigma(A_1)\star r=((0,-1),3,1)$; joining $\Sigma(B_1)$ adds $\rho^3(0,1)=(-1,0)$ and gives $((-1,-1),1,2)$; joining $f$ adds $\rho(1,0)=(0,-1)$ and gives $((-1,-2),1,3)$; and the last $r$ gives $\Sigma(A_2)=((-1,-2),2,3)$. For $D_2=\mathtt FA_2$ from the start $l$, the first $\mathtt F$ leads to $(0,1)$ and $A_2$ adds $\rho^3(-1,-2)=(2,-1)$, so $D_2$ ends at $(2,0)$, as in Section 1. Continuing, $\Sigma(A_3)=((-3,-2),2,7)$, $\Sigma(B_3)=((3,2),2,7)$, $\Sigma(A_4)=((-5,0),2,15)$ and $\Sigma(B_4)=((5,0),2,15)$. All of them have $d=2$, and by induction so does every nonempty block, because the headings in the rules add up to $2+1+2+0+1$ and $3+0+2+3+2$, both $\equiv2\pmod4$: each block turns the cursor around.

7. Walking a prefix

To find the position after $N$ steps of $D_n$, where $1\le N\le2^n$, we keep the state $S$ of the cursor and a budget $m$, the number of steps still to be made. We execute the first $\mathtt F$, so that $m=N-1$, and then walk the block $A_n$. A block $X$ of order $k$, which is $A_k$ or $B_k$, is walked as follows.

If $m=0$ or $k=0$, nothing happens.
If $m\ge2^k-1$, the whole block fits: $S$ becomes $S\star\Sigma(X)$, and $m$ decreases by $2^k-1$.
Otherwise the block is opened, and the five parts of its rule are visited from left to right: a turn is joined to $S$, an explicit $\mathtt F$ is joined and decreases $m$ by $1$, a block of order $k-1$ is walked by this same rule, and the walk stops as soon as $m=0$.

8. Why the walk is right

If a block $X$ of order $k$ is walked from the state $S$ with budget $m$, then for $m\ge2^k-1$ the walk ends in the state $S\star\Sigma(X)$ with budget $m-2^k+1$, and for $m\lt2^k-1$ it ends with budget $0$ and with the cursor where a cursor in state $S$ stands after the first $m$ steps of $X$. We use induction on $k$. For $k=0$ the block is empty, and the first case is Section 5. In the second case the parts are visited in order: turns and explicit steps are executed as written, and by the induction hypothesis every sub-block is either executed completely or stopped with budget $0$ at the right step. A sub-block that is executed completely leaves the full state, heading included, for the next part. As $X$ contains more than $m$ steps, the budget runs out inside $X$ and the walk stops there; after that no step is taken, so the turns it skips cannot change the position. For $D_n=\mathtt FA_n$ the walk thus ends at the position after $N$ steps.

9. Why the walk is short

A block is opened only when the budget runs out inside it, and inside an opened block the walk stops in the first part in which the budget runs out. So at most one block of each order is opened, and inside it at most two blocks of the next lower order are walked. The walk therefore handles at most $2n+1$ blocks, each with a fixed number of joins, instead of executing about $3N$ instructions.

10. Worked example: the checkpoint

For $D_{10}$ and $N=500$ the budget after the first $\mathtt F$ is $499$. The blocks $A_{10}$ and $A_9$ have $1023$ and $511$ steps, more than $499$, so both are opened. The first part of $A_9$ is $A_8$, whose $255$ steps fit, and after $256$ steps the cursor stands at $(0,16)$, the end of $D_8$. Next come the $\mathtt R$ of $A_9$ and its block $B_8=\mathtt{LF}A_7\mathtt LB_7$, which is opened: $\mathtt L$, $\mathtt F$, the whole $A_7$ and $\mathtt L$ are executed, and $B_7$ is opened in the same way, and so on down to $B_4$, whose part $A_3$ is opened and ends the walk with its first part $A_2$. Apart from the turns, the walk executes $\mathtt FA_8$, $\mathtt FA_7$, $\mathtt FA_6$, $\mathtt FA_5$, $\mathtt FA_4$ and $\mathtt FA_2$. After $256$, $384$, $448$, $480$, $496$ and $500$ steps the cursor stands at $(0,16)$, $(8,8)$, $(16,8)$, $(20,12)$, $(20,16)$ and $(18,16)$, as in the statement. Each piece $\mathtt FA_e$ makes $2^e$ steps, so here the pieces simply follow the binary digits of $500=2^8+2^7+2^6+2^5+2^4+2^2$.

11. Folding paper

The Heighway dragon was found by folding a strip of paper in half again and again and opening every fold to a right angle; its turns form the regular paperfolding sequence. ^[2] Let $\overline W$ be the word $W$ read backwards with $\mathtt L$ and $\mathtt R$ exchanged. Then $B_k=\overline{A_k}$ for all $k$: both are empty for $k=0$, and if $B_k=\overline{A_k}$, then also $A_k=\overline{B_k}$, so

$$\begin{aligned}\overline{A_{k+1}}&=\mathtt{LF}\,\overline{B_k}\,\mathtt L\,\overline{A_k}\\ &=\mathtt{LF}A_k\mathtt LB_k=B_{k+1}.\end{aligned}$$

Hence $B_n\mathtt F=\overline{D_n}$ and $D_{n+1}=D_n\,\mathtt R\,\overline{D_n}\,\mathtt R$; read as a sequence of net turns between steps, this is the construction described in ^[2]. A path walked backwards meets its turns in reverse order with left and right exchanged. If the last step of $D_n$ has heading $h$, the turns after it are undone at the start of $\overline{D_n}$, so the second half of $D_{n+1}$ starts with heading $h+1$, a quarter turn anticlockwise from the heading $h+2$ of the backward walk. So the second half is $D_n$ walked backwards and turned by $90^\circ$ anticlockwise about the end point $e_n$ of $D_n$: for $2^n\lt N\le2^{n+1}$, the point after $N$ steps is the image of the point after $2^{n+1}-N$ steps under this quarter turn. For $N=2^{n+1}$ the origin is carried to $e_{n+1}=e_n+\rho(e_n)$, so the end point moves $\sqrt2$ times farther out and turns by $45^\circ$ clockwise from one order to the next: this is the spiral of Section 1. The turns also follow the binary digits: if $j=2^tq$ with $q$ odd, the net turn after the $j$-th step is $\mathtt R$ when $q\equiv1\pmod4$ and $\mathtt L$ when $q\equiv3\pmod4$; so step $6=2\cdot3$ is followed by $\mathtt L$ and step $10=2\cdot5$ by $\mathtt R$. ^[6]

How the Code Works

The C++ program ^[7] stores an effect $(v,d,s)$ as a Segment with the fields x, y, dir and steps, using signed 64-bit coordinates and an unsigned 64-bit step count. combine(lhs, rhs) is the join $\star$: it turns the displacement of rhs by lhs.dir quarter turns $(x,y)\mapsto(y,-x)$, adds displacements and step counts, and adds the headings with & 3, that is, modulo $4$. The constants kF, kL and kR are $f$, $l$ and $r$, and the options --order=, --steps= and --skip-checkpoints override $n=50$ and $N=10^{12}$ or skip the check.

The constructor of Walker fills the vectors a and b from the back: entry $i$ holds $\Sigma(A_{n-i})$ or $\Sigma(B_{n-i})$, so index $n$ holds the empty blocks, and each pass evaluates the formulas of Section 6. The state starts as $((0,0),3,0)$, and remaining is the budget $m$. run executes the first $\mathtt F$ and calls execute_a(0). execute_a and execute_b walk the block of order $n$ minus level as in Section 7: they return when the budget is $0$ or the level is $n$, apply the stored effect when remaining is at least its step count, and otherwise go through $A_{k-1}\mathtt RB_{k-1}\mathtt{FR}$ or $\mathtt{LF}A_{k-1}\mathtt LB_{k-1}$, checking the budget after every part. run_checkpoints requires solve(10, 500) to return $(18,16)$, or the program exits with code $2$; then main prints the position. If --steps exceeds $2^n$, the walk executes all of $D_n$ and the program prints $e_n$ without a warning. The parsers check neither for overflow nor for orders above $64$, whose step counts no longer fit in 64 bits; the defaults are far from these limits.

The Python and Java versions follow the same plan, with the order and the number of steps fixed and without the checkpoint. Python keeps effects as tuples and walks with nested functions that update the state and the budget through nonlocal; its sys.setrecursionlimit(200) lowers the default limit but leaves ample room for a recursion about $n$ levels deep. Java stores an effect as four long values and the state in the static fields sx, sy, sd and remaining, and walks with exeA and exeB. The largest step count, $2^{50}-1$, needs more than 32 bits, and no coordinate can exceed the number of steps, so 64-bit integers suffice.

Complexity and Verification

The tables take $2n$ evaluations of Section 6, each with four joins of at most three quarter turns, and by Section 9 the walk handles at most $2n+1$ blocks, so time and memory are $O(n)$, counting an operation on a 64-bit integer as one step, with a recursion at most $n+1$ levels deep. For the target, the walk opens the eleven blocks $A_{50},\dots,A_{40}$, which are longer than the budget, and in all executes $13$ whole blocks, each right after an explicit $\mathtt F$ and one for each binary digit $1$ of $10^{12}$; it stops after $51$ calls of execute_a and execute_b and $50$ joins. Each program runs in a few hundredths of a second, mostly start-up time.

Unless --skip-checkpoints is given, the C++ program checks the checkpoint before the main computation; the Python and Java versions contain no check. For this article we also verified:

for every $n\le14$ and every $N$ with $0\le N\le2^n$, $32\,782$ cases, a transcription of the walk agrees with rewriting $D_n$ letter by letter and drawing it;
the tables agree with the words $A_k$ and $B_k$ executed instruction by instruction for $k\le10$, and $d=2$ and the end-point rule of Section 11 hold up to order $60$;
with --order and --steps, the C++ program agrees with a separate program that expands the rules recursively and executes every instruction, for $N=10^6$, $10^9$, $2^{32}-1$ and $12\,345\,678\,901$;
the three programs print the same position, and so do the folding recursion of Section 11, which shares no code with them, and the instruction-by-instruction program of the previous item, which took about $36$ minutes for all $10^{12}$ steps of $D_{50}$.

Footnotes and References

The references below supply the background facts; the effects of the blocks, the walk with its proof and cost, and the folding symmetry of the words are derived in this article.

Project Euler 220 — Heighway Dragon. The official statement defines D₀ = Fa, the rules a → aRbFR and b → LFaLb and the meaning of the letters, starts the cursor at (0, 0) pointing towards (0, 1), notes the position (18, 16) after 500 steps of D₁₀ and asks for the position after 10¹² steps of D₅₀.
Wikipedia — Dragon curve. The Heighway dragon, first investigated by the NASA physicists John Heighway, Bruce Banks and William Harter, can be constructed by folding a strip of paper; the turn sequence of each iteration is the previous one, then R, then the previous one in reverse order with L and R swapped.
Wikipedia — L-system. A parallel rewriting system: as many rules as possible are applied simultaneously in each iteration, and a symbol without a rule is a constant that stays unchanged.
Wikipedia — Turtle graphics. The turtle has a location, an orientation and a pen, and its commands are relative to its own position.
Wikipedia — Rotation matrix. The anticlockwise rotation through 270 degrees has the matrix with the rows (0, 1) and (−1, 0), which sends (x, y) to (y, −x).
OEIS A014577 — The regular paper-folding sequence. It gives the turns of the Heighway dragon, 1 for right and 0 for left, and a(n) is the complement of the bit to the left of the least significant 1 in the binary expansion of n + 1; for j = n + 1 this is the rule on q in Section 11.
C++ — Euler220.cpp. The linked immutable C++ revision contains combine, the Walker with execute_a and execute_b, and run_checkpoints. It is the implementation record for this article.

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] };
    }
}