Project Euler problem 312 is one of the harder problems we have worked through. It combines graph theory, fractal geometry, combinatorial counting, and modular arithmetic, and the final answer requires computing C(C(C(10000))) mod 13^8 where C(n) is the number of Hamiltonian cycles on a Sierpiński graph of order n.

We ended up leaning on OEIS to get the recurrence in closed form, which we will confess to below. Full notes: Project Euler/312.

The Problem

A Sierpiński graph of order 1, denoted S_1, is an equilateral triangle - three vertices, three edges.

S_{n+1} is built from three copies of S_n positioned so each pair of copies shares one corner vertex.

Let C(n) be the number of cycles that pass through every vertex of S_n exactly once - the number of Hamiltonian cycles. The problem tells us C(1) = C(2) = 1, and C(5) = 71,328,803,586,048.

The goal: find C(C(C(10000))) mod 13^8.

Two things to notice about the goal. First, C(10000) is astronomically large - you cannot actually compute it and then feed it to C again. Second, mod 13^8 is a big hint: the whole calculation is going to happen in modular arithmetic, and there is going to be some periodicity trick.

Vertices and Edges

Some easy facts to nail down before doing anything clever.

Number of edges of S_n:

Number of vertices of S_n:

So S_10000 has on the order of 3^10000 vertices, which is a number with about 4,700 decimal digits. A Hamiltonian cycle visits every vertex, so the number of Hamiltonian cycles is going to have on the order of 3^10000 factorial-ish in its combinatorics. There is no approach that involves listing anything.

The Recursive Structure

The critical observation is that S_{n+1} is built from three copies of S_n. A Hamiltonian cycle on S_{n+1} has to traverse all three copies, and the only places it can move between copies are at the shared corner vertices. So any Hamiltonian cycle on S_{n+1} consists of:

• A Hamiltonian path through the first copy of S_n, entering and exiting at two of its three outer corners
• A Hamiltonian path through the second copy, similarly
• A Hamiltonian path through the third copy, similarly
• Plus the three edges that connect the copies together at the shared corners

So counting Hamiltonian cycles reduces to counting certain kinds of Hamiltonian paths on smaller graphs.

Let

• C(n) = number of Hamiltonian cycles on S_n
• P(n) = number of Hamiltonian paths on S_n where both endpoints are outer corners
• P̄(n) = number of Hamiltonian paths on S_n where exactly one endpoint is an outer corner

Then

(three copies, each contributing a Hamiltonian path with both endpoints at outer corners, one path per copy, chosen independently).

Similarly, by tracking which paths piece together into which combined paths, we get

These recurrences let us compute P(n), P̄(n), and C(n) in principle. The values grow explosively - by the time n = 5, C(n) is already 7.1 × 10^13.

Closed Form (Confession Time)

At this point we did the thing everyone does: computed C(1) through C(5) from the recurrence, dropped the sequence into OEIS, and immediately hit A246959 - the "number of Hamiltonian cycles in the Sierpiński graph."

The closed form given by OEIS is

which is very compact. Equivalently, the recurrence

Sanity check against the given value C(5) = 71,328,803,586,048. Wolfram Alpha confirms.

First few values:

Note that C(2) = 2 above, but the problem states C(1) = C(2) = 1. The mismatch is a convention issue about what counts as a Hamiltonian cycle on the tiniest graphs (traversal direction, effectively). We used the OEIS convention for anything n ≥ 3, which is where the counting matters.

Nested Evaluation Under a Modulus

Now for the actual goal: C(C(C(10000))) mod 13^8.

The magnitude of C(10000) is not something you can hold in memory. By n = 20 the exponent in 12^((3^(n-2) - 3)/2) is already millions of digits long, and we are being asked about n = 10000. So the whole calculation has to happen without ever materializing an intermediate value in full.

The trick is that we don't need the intermediate values in full. We need C(C(C(10000))) mod 13^8. That means for the outermost C, we only need its argument modulo something that makes the outermost C computable. And that "something" is much smaller than the argument itself.

Three techniques do the work.

Technique 1: Periodicity of C(n) mod M

The sequence C(1), C(2), C(3), ... taken modulo any fixed M is eventually periodic. It has to be - the state at step n in the recurrence C(n) = (3 · C(n-1))^3 mod M is a single residue in Z/MZ, and there are only M possible residues, so the sequence must eventually enter a cycle.

Let P be the period. Then

for X at least as large as the pre-period offset X_0. In practice you compute successive C(n) mod M values until you see a repeat, and that tells you X_0 and P.

This is the workhorse. It lets you reduce "compute C at an astronomically large index" to "compute C at an index in the range [X_0, X_0 + P)."

Technique 2: Choosing the Right Modulus at Each Layer

The three C calls need different moduli.

The outermost C operates modulo 13^8. That is what the problem asks for.

The middle C needs to produce a residue that determines the outermost C(x) mod 13^8. By technique 1, the outermost C is periodic mod 13^8 with some period P_out. So the middle C only needs to compute its result modulo P_out (well, modulo the pre-period plus period - be careful with the offset).

The innermost C(10000) needs to produce a residue that determines the middle C(y) mod P_out. So find the period P_mid of C mod P_out, and compute C(10000) mod P_mid.

So the algorithm is:

1. Find the period of C(n) mod 13^8. Call it P_out.
2. Find the period of C(n) mod P_out. Call it P_mid.
3. Compute C(10000) mod P_mid. Call the result a.
4. Compute C(a) mod P_out. Call the result b.
5. Compute C(b) mod 13^8.

Each step uses either the recurrence directly (b = (3·b_prev)^3 mod M) or, for the innermost step at n = 10000, the recurrence run 10000 times mod the appropriate modulus. Ten thousand modular multiplications is nothing.

Technique 3: Handling the Cube Under a Prime Power

The recurrence C(n) = (3 · C(n-1))^3 involves a cube, so under a prime power modulus you have to be careful about invertibility. Cubing is well-defined mod 13^8 (or mod anything else), so the forward recurrence is not a problem - you just cube and multiply. What can go wrong is if you try to invert the cubing at some point, which you would need if you tried to do matrix exponentiation on the logarithm of the recurrence.

We did not need to invert. The three-layer periodicity approach above uses only forward evaluation of the recurrence, so cubing mod 13^8 is fine.

Chinese Remainder Theorem

CRT is worth mentioning for completeness. If the modulus had been composite with multiple distinct prime factors, we would compute the answer separately modulo each prime power and glue the results together with CRT.

Here the modulus is 13^8, a single prime power, so CRT is not needed - but if the problem had asked for mod (7^5 · 11^3 · 13^8), the approach would be to solve three separate problems and CRT the answers.

What We Learned

Three things worth carrying forward.

Recursive structure beats brute force by a lot. The naive approach of enumerating cycles on S_n is impossible for anything n > 4 or so. The recursive approach reduces the problem to counting paths on S_{n-1} with three specific endpoint types, and that recursion collapses to a single-variable formula. Every time you can find a recursive structure, use it.

OEIS is not cheating. Recognizing a sequence and looking up its closed form is a legitimate research tool, not shortcut. The people who computed and cataloged OEIS A246959 did the work of proving the closed form. Standing on their shoulders is what mathematics has always looked like. The alternative - re-deriving every known result from scratch - is how you never finish anything.

Nested modular evaluation is a real technique. When a problem asks for f(f(f(n))) mod M with a huge n, you rarely need to compute the intermediate values. You work outward-in through periods: the outermost period tells you what modulus the middle layer actually needs to produce, which tells you what modulus the innermost layer needs to produce. Astronomical intermediate values reduce to residues in tractable rings. This shape of problem shows up a lot in number theory and cryptography, and it is worth having the pattern in your head.

References

• Our wiki notes: Project Euler/312
• PE 312 on projecteuler.net
• OEIS A246959 - Hamiltonian cycles in Sierpiński graphs
• Useful paper: https://arxiv.org/pdf/0909.5541 (Hamiltonian properties of Sierpiński graphs)

Project Euler problem 227 - "The Chase" - is a nice example of a problem where the brute force approach is completely reasonable to reach for first, completely infeasible to actually finish, and directly points you at the right mathematical reformulation.

This is what we ended up with. Our full working notes are on the wiki at Project Euler/227.

The Problem

An even number of players sit around a table. Two players sitting directly opposite each other each start with a die. Each round, both players roll:

• Roll a 1 → pass the die to your left neighbor
• Roll a 6 → pass the die to your right neighbor
• Anything else → keep the die

The game ends when one player is holding both dice at the end of a round. That player loses.

The question: with 100 players, what is the expected number of rounds before someone loses?

Reaching For the Brute Force

The first thing we did was write a very small Python simulator. Two integer pointers, one per die. Each round, roll two dice, update each pointer with +1, -1, or nothing, modulo the number of players. Stop when the pointers collide. Average over many trials.

This works. It gives you an estimate that gets more accurate as you run more trials. And then you start measuring how expensive that gets.

The cost has two knobs:

• Trials. Linear. Four decimal places of accuracy needs about 1,000 trials. Ten decimal places needs billions.
• Players. Not linear. Playing with 5, 10, 15, 20, and 25 players and measuring, the runtime looked quadratic in the number of players.

So the total execution time is roughly

Extrapolating our measured constant c from small experiments up to N = 100 players and S = 9 billion trials (to get sufficient accuracy on the answer), we got about 7,000,000 seconds. Call it 80 days.

Project Euler problems are supposed to have solutions that run in about a minute. 80 days is not a minute. Even switching to a compiled language and multithreading, we would need a 10,000× speedup to get into the right ballpark. And the loop body is already about as tight as a loop body gets - two random integers, two pointer updates, a comparison. There is not much to optimize.

So the brute force is out. What is it about the problem that makes brute force so bad, and what does the answer look like if we ask a different question?

Reformulating as a Markov Chain

The key observation is that the game has no memory. The next state only depends on the current state. That is the definition of a Markov chain.

The state is the pair of positions of the two dice, (i, j), where i and j are player indices in {0, 1, ..., N-1}. For 100 players, that's a state space with 10,000 pairs, minus the 100 pairs where i = j (dice have collided - game over). So we have 9,900 transient states and 100 absorbing states.

Note: (i, j) and (j, i) represent the same physical situation. There is a symmetry we could exploit to cut the state count roughly in half, but for the initial derivation we ignored it and paid the extra factor.

Transitions: each die independently keeps, moves left, or moves right, with probabilities 4/6, 1/6, 1/6. So the joint transition from (i, j) to (i', j') has probability equal to the product of the two individual probabilities. Each transient state has at most 9 possible next states (3 outcomes per die × 3 outcomes per die).

The Fundamental Matrix

The clean way to compute expected time to absorption in a finite Markov chain is via the fundamental matrix. It works like this.

Arrange the states so transient states come first and absorbing states come second. Then the transition probability matrix has the block form

where

• Q is t × t, the transition probabilities among transient states
• R is t × r, the transitions from transient to absorbing states
• 0 and I in the bottom row represent the fact that absorbing states don't leave themselves

The fundamental matrix is

Entry N[i, j] gives the expected number of visits to transient state j before absorption, starting from transient state i. That means the expected number of steps until absorption, starting from state i, is the sum of row i of N.

Or equivalently: solve the linear system (I - Q) x = 1, where 1 is a column of ones. Entry x[i] is the expected number of steps starting from state i.

Building Q for 100 Players

For 100 players there are 100 × 99 = 9900 transient states. Q is 9900 × 9900. We only need to fill in the roughly nine nonzero entries per row, so Q is very sparse.

Two things to get right when you build it:

1. State-to-index mapping. Some canonical ordering of the (i, j) pairs with i ≠ j. Any consistent bijection between pairs and row indices works.
2. Iterating over next states. For each (i, j), enumerate the nine (i', j') transitions. For each, check whether i' = j' (in which case that probability mass flows to an absorbing state and doesn't go into Q).

Then invert I - Q, sum the row of N that corresponds to the starting state (0, 50) (opposite sides of the table), and you have your answer.

The Uncomfortable Confession

Ours took 1,687 seconds to run. Twenty-eight minutes. Better than 80 days, but not great.

The slow part is inverting the 9,900 × 9,900 matrix. We tried reformulating as a linear solve (I - Q) x = 1 and using LU decomposition instead of computing the full inverse, on the theory that LU should be faster. It was almost twice as slow. That was a surprise but probably a story about the specific dense linear algebra library we were using rather than a fundamental issue.

The obvious next moves we did not take:

• Sparse matrix representation. Q has at most 9 nonzeros per row out of 9,900. Using a proper sparse solver would probably drop us into the seconds.
• Exploit the symmetry. Treating (i, j) and (j, i) as the same state cuts the matrix size in half. Combined with sparsity, this should be very fast.
• Iterative methods. The system (I - Q) x = 1 is a great candidate for something like BiCGSTAB or GMRES, which don't build the inverse at all.

The Real Solution Times

For fun, we looked at the Project Euler forum posts after we submitted. The first person to solve it did it in 0.1 seconds in C++ in 2009. That is 16,000× faster than ours.

They used a completely different approach, involving grouping the six die outcomes into three groups (+1, -1, keep) and using symmetries to reduce the state space dramatically. There is a shape of insight that turns "big linear algebra problem" into "small closed-form recurrence" that we did not find.

Some of the fun of Project Euler is exactly this: finishing the problem, submitting the answer, and then finding out that half the people who solved it did it 10,000× faster with a technique you did not think of. The correct response is not embarrassment. It is to open the forum, read their write-ups, and steal their trick for next time.

References

• Our full notes: Project Euler/227
• PE 227 on projecteuler.net
• Fundamental matrix: standard reference is Kemeny & Snell, Finite Markov Chains (1960)

Project Euler problem 198 is a number theory problem that turns out to be almost entirely about continued fractions. The problem hides this - the statement talks about "ambiguous" real numbers - but the ambiguity has a clean characterization in terms of continued fraction expansions, and once you see it, the problem gets a lot smaller.

Wiki notes: Project Euler/198.

The Problem

Define a best approximation to a real number x with denominator bound d as a rational r/s in reduced form with s ≤ d, such that any other rational p/q closer to x than r/s has q > d. In other words: it is the rational number with the smallest denominator that gets as close to x as anything with a small enough denominator can.

Usually the best approximation is unique for each denominator bound. But sometimes there are two equally-good best approximations. The problem gives the example of x = 9/40, which has both 1/4 and 1/5 as best approximations for the denominator bound d = 6. Both are tied for closest to 9/40 among rationals with denominator at most 6.

Call a real number x ambiguous if there is at least one denominator bound for which it has two best approximations. Ambiguous numbers are necessarily rational.

The question: how many ambiguous x = p/q with 0 < x < 1/100 and q ≤ 10^8 are there?

Why Continued Fractions

Any real number x has a continued fraction expansion

which we write compactly as [a_0; a_1, a_2, a_3, ...]. For rational numbers, the expansion terminates. For irrationals, it does not.

The convergents of a continued fraction are the rationals you get by truncating the expansion at each step. The convergents of [a_0; a_1, a_2, ..., a_n] are

There is a beautiful classical theorem that says: the convergents of the continued fraction expansion of x are exactly the best approximations to x. Every convergent is a best approximation for some denominator bound, and every best approximation is a convergent or a certain kind of intermediate fraction.

For x = 9/40:

• CF expansion: [0; 4, 2, 4]
• First convergent: [0; 4] = 1/4
• Second convergent: [0; 4, 2] = 1 / (4 + 1/2) = 2/9
• Third (final) convergent: [0; 4, 2, 4] = 9/40

So 1/4 is one best approximation for 9/40 at some bound. Where does the other one, 1/5, come from?

Where Ambiguity Comes From

The answer is that not all best approximations are convergents. There are also intermediate fractions, which sit between successive convergents. Specifically, when a partial quotient a_k is larger than 1, you can form intermediate fractions

for j = 1, 2, ..., a_{k+1} - 1, where p_k / q_k is the k-th convergent. These intermediate fractions are best approximations for certain denominator bounds.

The situation where a real number has two best approximations for some denominator bound is exactly when a convergent and an intermediate fraction happen to sit equidistant from x. This can only happen when the continued fraction expansion has some partial quotient bigger than 1 sitting in the right place.

The clean characterization we ended up with is this:

A rational p/q is non-ambiguous if and only if its continued fraction expansion has the form

Everything else is ambiguous.

This is the crux of the problem. Instead of counting ambiguous fractions directly, count the total number of relevant fractions and subtract the non-ambiguous ones - or, better, characterize the ambiguous ones directly by generating their continued fraction expansions.

Constraints From the Problem

The problem restricts to 0 < p/q < 1/100 and q ≤ 10^8.

Since p ≥ 1 and p/q < 1/100, we have q > 100 p ≥ 100. So q ranges from 101 up to 10^8.

For a given q, valid p values satisfy:

• p ≥ 1
• p ≤ q/100
• gcd(p, q) = 1 (reduced form)

That gives us roughly q/100 candidate numerators per denominator, but we also have to intersect with "coprime to q" and "continued fraction has the right shape to be ambiguous."

The Search

The high-level approach:

1. Iterate over denominators q from 101 to 10^8
2. For each q, enumerate valid p values (coprime to q, in range)
3. Compute the continued fraction expansion of p/q
4. Check whether the expansion has the "ambiguous" shape
5. Count

The bottleneck is that step 1 is 10^8 iterations, and even a very fast inner loop times out. So the real approach is to invert the problem: instead of iterating over fractions and checking whether their CF is ambiguous, iterate over the shapes of ambiguous CF expansions and count the fractions those shapes produce that fall in the range.

This is a nice example of a general technique: when you have a characterization of the objects you want to count in terms of some structural property, generating objects by structure is usually faster than filtering.

The Bigger Point

Problems that look purely computational often turn out to have a clean structural characterization hiding underneath. The PE 198 statement never mentions continued fractions - it's about "best approximations" and "ambiguous real numbers." But the moment you dig into what best approximation actually means for rationals, continued fractions appear naturally, and the entire problem reduces to counting fractions whose CF has a certain shape.

This kind of translation - from a definition in one language to a characterization in another - is where number theory earns its reputation for being satisfying. You start with "count the ambiguous fractions" and you end with "count the fractions whose CF is not of the form [0; a₁, 1, 1, ..., 1]," and the second version is much easier to work with.

References

• Our wiki notes: Project Euler/198
• PE 198 on projecteuler.net
• Classical reference: Hardy & Wright, An Introduction to the Theory of Numbers, chapters on continued fractions