Project Euler Problem 946

Solution

Answer: 585787007

Let

$$\alpha=[2;\underbrace{1,\dots,1}{2},2,\underbrace{1,\dots,1}{3},2, \underbrace{1,\dots,1}_{5},2,\ldots]$$

where the lengths of the blocks of $1$'s are the consecutive primes.

We define

$$\beta=\frac{2\alpha+3}{3\alpha+2}.$$

The task is to compute the sum of the first $10^8$ continued-fraction coefficients of $\beta$.

A standard way to transform continued fractions under a Möbius map

$$x\mapsto \frac{ax+b}{cx+d}$$

is to use the matrix

$$\begin{pmatrix}a&b\ c&d\end{pmatrix}.$$

For each incoming continued-fraction digit $t$, the matrix updates as

$$(a,b,c,d)\mapsto (at+b,\ a,\ ct+d,\ c).$$

Whenever

$$\left\lfloor \frac a c \right\rfloor = \left\lfloor \frac b d \right\rfloor,$$

that common value is the next continued-fraction coefficient of the transformed number.

Here the initial matrix is

$$\begin{pmatrix}2&3\3&2\end{pmatrix},$$

because

$$\beta=\frac{2\alpha+3}{3\alpha+2}.$$

Running this transducer over the continued fraction of $\alpha$ efficiently generates the coefficients of $\beta$. The first ten are indeed

$$[0;1,5,6,16,9,1,10,16,11],$$

whose sum is $75$, matching the statement.

An optimized implementation (constant-size integer state, no big integers needed) computes the first $10^8$ coefficients in under a few seconds.

Python version of the core logic:

from sympy import primerange

N = 10**8

a, b, c, d = 2, 3, 3, 2
count = 0
total = 0

def process(t):
    global a, b, c, d, count, total

    a, b, c, d = a*t + b, a, c*t + d, c

    while c and d and a // c == b // d:
        q = a // c
        total += q
        count += 1

        if count == N:
            return True

        a, b, c, d = c, d, a - q*c, b - q*d

    return False

# initial term of alpha
process(2)

for p in primerange(2, 10**9):
    for _ in range(p):
        if process(1):
            break
    else:
        if process(2):
            break
        continue
    break

print(total)

This yields:

Answer: 585787007