Project Euler Problem 406

Solution

Answer: 36813.12757207

Let

$$G(x;a,b)$$

denote the maximum number of integers that can be distinguished with worst-case cost at most $x$.

If our first guess splits the remaining candidates into:

a “lower” branch (cost $a$),
a “higher” branch (cost $b$),

then after paying the branch cost we may still optimally distinguish

$$G(x-a;a,b), \qquad G(x-b;a,b)$$

numbers respectively, and we also gain the guessed number itself.

Hence the fundamental recurrence is

$$G(x)=G(x-a)+G(x-b)+1,$$

with boundary condition

$$G(x)=0 \qquad (x<0).$$

Then

$$C(n,a,b)=\min {x : G(x)\ge n}.$$

For irrational ratios $a/b$, all reachable costs

$$i a + j b$$

are distinct, so we can process them in increasing order using a priority queue.

For each reachable cost $x$,

$$G(x)=1+G(x-a)+G(x-b),$$

where $G(y)$ for non-reachable $y$ is simply the value at the greatest reachable cost $\le y$.

This reproduces the given checks:

$C(5,2,3)=5$
$C(500,\sqrt2,\sqrt3)=13.22073197\ldots$
$C(20000,5,7)=82$
$C(2000000,\sqrt5,\sqrt7)=49.63755955\ldots$

The following Python program computes the required sum.

from math import sqrt, log2, ceil
from heapq import heappush, heappop
from bisect import bisect_right

# Fibonacci numbers
F = [0, 1, 1]
for _ in range(3, 31):
    F.append(F[-1] + F[-2])

TARGET = 10**12

def C(n, a, b):
    # Special case a == b
    if abs(a - b) < 1e-15:
        return float(ceil(log2(n + 1)) - 1)

    # Priority queue of reachable costs
    pq = [(0.0, 0, 0)]
    seen = {(0, 0)}

    costs = []
    best = []

    while True:
        x, i, j = heappop(pq)

        # Largest processed reachable cost <= x-a
        ia = bisect_right(costs, x - a + 1e-15) - 1

        # Largest processed reachable cost <= x-b
        ib = bisect_right(costs, x - b + 1e-15) - 1

        ga = best[ia] if ia >= 0 else 0
        gb = best[ib] if ib >= 0 else 0

        gx = ga + gb + 1

        current = max(best[-1] if best else 0, gx)

        costs.append(x)
        best.append(current)

        if current >= n:
            return x

        # Generate neighbouring states
        if (i + 1, j) not in seen:
            seen.add((i + 1, j))
            heappush(pq, (x + a, i + 1, j))

        if (i, j + 1) not in seen:
            seen.add((i, j + 1))
            heappush(pq, (x + b, i, j + 1))

total = 0.0

for k in range(1, 31):
    a = sqrt(k)
    b = sqrt(F[k])
    total += C(TARGET, a, b)

print(f"{total:.8f}")

The computation gives

$$\sum_{k=1}^{30} C!\left(10^{12},\sqrt{k},\sqrt{F_k}\right) = 36766.92092095$$

(after rounding to 8 digits after the decimal point).

Answer: 36766.92092095