Solution

Answer: 605857431263981935

Let

$$N(p,q)=\sum_{n=0}^{q} T_n p^n,$$

where

$$S_0=290797,\qquad S_{n+1}=S_n^2 \bmod 50515093,\qquad T_n=S_n\bmod p.$$

We want

$$\operatorname{NF}(p,q)=v_p(N(p,q)!),$$

the exponent of the prime $p$ in $N(p,q)!$.

For this problem:

$$p=61,\qquad q=10^7,$$

and we need

$$\operatorname{NF}(61,10^7)\bmod 61^{10}.$$

1. Mathematical analysis

The key fact is Legendre’s formula:

$$v_p(n!)=\sum_{k\ge1}\left\lfloor \frac{n}{p^k}\right\rfloor.$$

A standard identity derived from base-$p$ digit sums is

$$v_p(n!)=\frac{n-s_p(n)}{p-1},$$

where $s_p(n)$ is the sum of the digits of $n$ in base $p$.

Step 1: Base-$p$ representation of $N(p,q)$

Because each

$$0\le T_n < p,$$

the number

$$N(p,q)=\sum_{n=0}^{q}T_n p^n$$

already is the base-$p$ expansion of $N$.

Therefore,

$$s_p(N)=\sum_{n=0}^{q}T_n.$$

Hence

$$\operatorname{NF}(p,q) = \frac{N-\sum T_n}{p-1}.$$

Substitute the definition of $N$:

$$\operatorname{NF}(p,q) = \frac{\sum_{n=0}^{q}T_n p^n-\sum_{n=0}^{q}T_n}{p-1}.$$

Factor term-by-term:

$$\operatorname{NF}(p,q) = \sum_{n=0}^{q} T_n\frac{p^n-1}{p-1}.$$

Define

$$R_n=\frac{p^n-1}{p-1}=1+p+p^2+\cdots+p^{n-1}.$$

$$\boxed{ \operatorname{NF}(p,q)=\sum_{n=0}^{q}T_nR_n }.$$

Step 2: Reduce modulo $p^{10}$

We only need modulo $61^{10}$.

Observe:

$$R_n = 1+p+\cdots+p^{n-1}.$$

For any $n\ge 10$,

$$R_n \equiv 1+p+\cdots+p^9 \pmod{p^{10}},$$

because all higher powers are divisible by $p^{10}$.

Let

$$R=1+p+\cdots+p^9.$$

Then

$$\operatorname{NF}(p,q) \equiv \sum_{n=0}^{9}T_nR_n + R\sum_{n=10}^{q}T_n \pmod{p^{10}}.$$

This is the crucial simplification: we never need huge powers of $p$, only the first ten.

2. Python implementation

def solve():
    p = 61
    q = 10_000_000
    mod = p ** 10

    # Random generator seed
    S = 290797

    # R = 1 + p + p^2 + ... + p^9  (mod p^10)
    R = sum(p ** i for i in range(10)) % mod

    # Store first 10 T_n values separately
    first_terms = []

    # Sum of all T_n for n >= 10
    tail_sum = 0

    # Generate T_n values
    for n in range(q + 1):
        T = S % p

        if n < 10:
            first_terms.append(T)
        else:
            tail_sum += T

        S = (S * S) % 50515093

    # Compute contribution of first 10 terms exactly
    ans = 0

    for n in range(10):
        # R_n = 1 + p + ... + p^(n-1)
        Rn = sum(p ** i for i in range(n))
        ans += first_terms[n] * Rn

    # Contribution from all remaining terms
    ans += R * tail_sum

    ans %= mod

    return ans

print(solve())

3. Code walkthrough

Initialization

p = 61
q = 10_000_000
mod = p ** 10

We work modulo $61^{10}$.

Compute the stabilized repunit

R = sum(p ** i for i in range(10)) % mod

For every $n \ge 10$,

$$R_n \equiv R \pmod{61^{10}}.$$

Generate the pseudo-random sequence

S = 290797

Then repeatedly:

T = S % p
S = (S * S) % 50515093

exactly as required by the problem statement.

Separate first 10 terms from the tail

if n < 10:
    first_terms.append(T)
else:
    tail_sum += T

The first ten coefficients have distinct weights $R_n$.

All later coefficients share the same weight $R$.

Final accumulation

For the first terms:

Rn = sum(p ** i for i in range(n))
ans += first_terms[n] * Rn

For the rest:

ans += R * tail_sum

Finally reduce modulo $61^{10}$.

4. Verification on the given example

Using the same method for $p=3$, $q=10000$:

624955285

which matches the problem statement:

$$\operatorname{NF}(3,10000)\bmod 3^{20}=624955285.$$

So the derivation is correct.

5. Final computation

Running the algorithm for

$$(p,q)=(61,10^7)$$

gives:

$$605857431263981935.$$

Answer: 605857431263981935

Project Euler Problem 288