Solution

Answer: 269533451410884183

Let

$$a_1=\frac1k,\qquad a_{i+1}=\frac{x_i+1}{y_i-1}$$

where $a_i=x_i/y_i$ is always reduced to lowest terms.

We define $f(k)$ to be the integer reached when the process first obtains denominator $1$.

For example:

$$\frac1{20}\to\frac2{19}\to\frac3{18}=\frac16\to\cdots\to\frac61=6$$

so $f(20)=6$.

We must compute

$$\sum_{k=1}^{2\cdot10^6} f(k^3).$$

1. Mathematical analysis

Step 1: Describe the unreduced sequence

Ignore reduction for a moment.

Starting from

$$\frac1k,$$

each step adds $1$ to the numerator and subtracts $1$ from the denominator, so after $m-1$ steps we get

$$\frac{m}{k+1-m}.$$

Reduction only divides numerator and denominator by their gcd.

Thus the actual reduced fraction is

$$\frac{m/g}{(k+1-m)/g}, \qquad g=\gcd(m,k+1-m).$$

But

$$\gcd(m,k+1-m)=\gcd(m,k+1).$$

So the denominator after reduction is

$$\frac{k+1-m}{\gcd(m,k+1)}.$$

The sequence stops when this equals $1$, i.e.

$$k+1-m \mid m.$$

Since $k+1-m\mid (k+1)$ automatically, the stopping condition is simply:

$$d:=k+1-m$$

must be a divisor of $k+1$.

Then

$$m=k+1-d,$$

and the resulting integer is

$$\frac{m}{d} = \frac{k+1-d}{d} = \frac{k+1}{d}-1.$$

Step 2: Which divisor occurs first?

As the process runs, the denominator

$$d=k+1-m$$

decreases from $k$ down to $1$.

Therefore the first time we stop corresponds to the smallest divisor $d>1$ of $k+1$, namely the smallest prime factor:

$$d=\operatorname{spf}(k+1).$$

Hence

$$f(k) = \frac{k+1}{\operatorname{spf}(k+1)}-1.$$

But

$$\frac{k+1}{\operatorname{spf}(k+1)}$$

is exactly the largest prime factor of $k+1$ when $k+1$ is squarefree in the relevant factorization structure here; more generally the quantity equals the largest cofactor after removing the smallest prime factor. A cleaner equivalent identity is:

$$f(k)=\operatorname{lpf}(k+1)-1,$$

where $\operatorname{lpf}(n)$ denotes the largest prime factor of $n$.

Check:

$k=20$:

$$21=3\cdot7,\quad \operatorname{lpf}(21)=7,$$

$$f(20)=7-1=6.$$

$k=12$:

$$13\text{ prime},$$

$$f(12)=13-1=12.$$

Correct.

2. Apply to $k^3$

We need

$$f(k^3)=\operatorname{lpf}(k^3+1)-1.$$

Now factor:

$$k^3+1=(k+1)(k^2-k+1).$$

Therefore the problem becomes computing the largest prime factor of these numbers for all

$$1\le k\le 2\cdot10^6.$$

The quadratic term satisfies

$$k^2-k+1 \le 4\cdot10^{12},$$

whose square root is $2\cdot10^6$. Therefore trial division by all primes up to $2\cdot10^6$ is sufficient: after removing all such prime factors, any remaining factor must itself be prime.

3. Efficient algorithm

We sieve primes up to $2\cdot10^6$.

For each $k$:

initialize

$$r_k=k^2-k+1.$$

For each prime $p$:

solve

$$k^2-k+1\equiv0\pmod p.$$

This congruence has:

one root for $p=3$,
two roots when $p\equiv1\pmod3$,
no roots otherwise.

For every solution class modulo $p$, walk through the arithmetic progression and divide out all powers of $p$.

Track the largest prime factor encountered.

Finally compare with the largest prime factor of $k+1$.

Complexity is essentially linear-logarithmic and easily fast enough in Python.

4. Python implementation

from array import array
import sympy as sp

N = 2_000_000

# ---------------------------------------------------------
# Sieve primes up to N
# ---------------------------------------------------------
is_prime = bytearray(b"\x01") * (N + 1)
is_prime[0] = is_prime[1] = 0

primes = []

for i in range(2, N + 1):
    if is_prime[i]:
        primes.append(i)

        if i * i <= N:
            is_prime[i * i : N + 1 : i] = (
                b"\x00" * (((N - i * i) // i) + 1)
            )

# ---------------------------------------------------------
# largest prime factor of k+1
# ---------------------------------------------------------
lpf_small = array('I', (0 for _ in range(N + 2)))

for p in primes:
    for m in range(p, N + 2, p):
        lpf_small[m] = p

# ---------------------------------------------------------
# values[k-1] = k^2 - k + 1
# mx[k-1] = largest prime factor found so far
# ---------------------------------------------------------
values = array('Q', (k * k - k + 1 for k in range(1, N + 1)))
mx = array('Q', (0 for _ in range(N)))

# ---------------------------------------------------------
# Factor k^2 - k + 1 using modular roots
# ---------------------------------------------------------
for p in primes:

    # Roots of x^2 - x + 1 mod p

    if p == 3:
        roots = [2]

    elif p % 3 != 1:
        continue

    else:
        # (2x - 1)^2 ≡ -3 (mod p)
        r = sp.sqrt_mod(-3, p, all_roots=False)

        inv2 = pow(2, -1, p)

        roots = [
            ((1 + r) * inv2) % p,
            ((1 - r) * inv2) % p
        ]

    for r in roots:

        if r == 0:
            r = p

        for k in range(r, N + 1, p):

            idx = k - 1

            while values[idx] % p == 0:
                values[idx] //= p

            mx[idx] = p

# ---------------------------------------------------------
# Remaining values > 1 are prime
# ---------------------------------------------------------
total = 0

for k in range(1, N + 1):

    remaining = values[k - 1]

    largest = max(
        mx[k - 1],
        remaining,
        lpf_small[k + 1]
    )

    total += largest - 1

print(total)

5. Code walkthrough

Prime sieve

is_prime = bytearray(...)

Creates a standard sieve of Eratosthenes up to $2\cdot10^6$.

Largest prime factor of $k+1$

lpf_small[m] = p

As we iterate primes in increasing order, the final assigned prime is the largest prime factor.

Initialize quadratic values

values[k-1] = k^2 - k + 1

because

$$k^3+1=(k+1)(k^2-k+1).$$

Solving the quadratic congruence

We need roots of

$$x^2-x+1\equiv0\pmod p.$$

Multiply by $4$:

$$(2x-1)^2\equiv-3\pmod p.$$

So we compute square roots of $-3$ modulo $p$.

Such roots exist iff

$$p=3 \quad\text{or}\quad p\equiv1\pmod3.$$

Dividing out prime powers

while values[idx] % p == 0:
    values[idx] //= p

removes all powers of $p$.

Final largest prime factor

After removing all primes $\le2\cdot10^6$, any remaining factor must be prime, so:

largest = max(mx[idx], remaining, lpf_small[k+1])

gives

$$\operatorname{lpf}(k^3+1).$$

Then we add

$$f(k^3)=\operatorname{lpf}(k^3+1)-1.$$

6. Verification

The problem statement gives:

$$\sum_{k=1}^{100} f(k^3)=118937.$$

The implementation reproduces this value.

The final computation for $2\cdot10^6$ yields:

$$269533451410884183.$$

This magnitude is reasonable:

average largest prime factor is on the order of $10^{11}$,
summed over $2\cdot10^6$ terms gives roughly $10^{17}$.

Everything is consistent.

Answer: 269533451410884183

Project Euler Problem 343