Project Euler Problem 811

Solution

Answer: 327287526

A key observation is that the recursion depends only on the pattern of zero-runs between the 1-bits of $n$.

For

$$n=(2^t+1)^r,$$

when $t$ is sufficiently large (and here $t=10^{14}+31$), the binary blocks coming from the binomial expansion

$$(1+2^t)^r=\sum_{k=0}^r \binom{r}{k}2^{kt}$$

do not overlap. Thus the 1-bit positions of $n$ are obtained directly from the binary expansions of the binomial coefficients.

Using the resulting closed-form evaluation of $A(n)$, one can compute

$$H(10^{14}+31,62)\pmod{1,000,062,031}.$$

The following Python program evaluates it efficiently:

MOD = 1_000_062_031

def b(n):
    return n & -n

def max_binom_bitlen(r):
    c = 1
    mx = 1
    for k in range(r + 1):
        mx = max(mx, c.bit_length())
        if k < r:
            c = c * (r - k) // (k + 1)
    return mx

def one_positions_via_binom(t, r):
    pos = []
    c = 1
    for k in range(r + 1):
        x = c
        while x:
            lsb = x & -x
            bit = lsb.bit_length() - 1
            pos.append(k * t + bit)
            x -= lsb
        if k < r:
            c = c * (r - k) // (k + 1)
    return pos

def A_from_positions(pos, mod):
    m = len(pos)

    v = [0] * m
    v[0] = 1
    for k in range(1, m):
        v[k] = (5 * v[k - 1] + 3) % mod

    ans = 1
    desc = pos[::-1]

    for i in range(m - 1):
        gap = desc[i] - desc[i + 1] - 1
        if gap > 0:
            ans = ans * pow(v[i + 1], gap, mod) % mod

    return ans

def H(t, r, mod):
    pos = one_positions_via_binom(t, r)
    return A_from_positions(pos, mod)

print(H(10**14 + 31, 62, MOD))

Verification against the sample:

$$H(3,2)=A(81)=636056,$$

which the program reproduces correctly.

Therefore,

Answer: 327287526