Solution

Answer: 950591530

Let

$$x_i(k)=i^{2^k}$$

be the value obtained from the initial number $i$ after it has been squared $k$ times.

At every step we choose the smallest current value and square it once more.

We must compute

$$S(10^4,10^{16}) \pmod{1234567891}.$$

1. Mathematical analysis

For each $i$, define

$$L_i(k)=\log_2(\log x_i(k)).$$

Since

$$x_i(k)=i^{2^k},$$

we have

$$\log x_i(k)=2^k\log i,$$

and therefore

$$L_i(k) =\log_2(2^k\log i) = k+\log_2(\log i).$$

Define

$$c_i=\log_2(\log i).$$

Then after $k_i$ squarings, the “priority value” of $i$ is simply

$$k_i+c_i.$$

Squaring the current minimum increases its priority by exactly $1$.

So the whole process becomes:

Start with numbers $c_i$.
Each operation picks the smallest current value.
Replace it by value $+1$.

This is a classic balancing/water-filling process.

Structure of the final counts

Suppose after all operations, number $i$ was squared $k_i$ times.

Then the values used from sequence $i$ are:

$$c_i,\ c_i+1,\ c_i+2,\ \dots,\ c_i+k_i-1.$$

Globally, the algorithm selects the $m$ smallest elements among all these arithmetic progressions.

Write

$$c_i=q_i+f_i,$$

where

$$q_i=\lfloor c_i\rfloor,\qquad 0\le f_i<1.$$

For sufficiently large $Y$,

$$k_i = Y-q_i$$

for everyone, plus possibly one extra operation for some indices.

Indeed,

$$\sum_i (Y-q_i) = (n-1)Y-\sum_i q_i.$$

Let

$$B(Y)=(n-1)Y-\sum_i q_i.$$

Choose the largest $Y$ such that

$$B(Y)\le m.$$

Then define

$$d=m-B(Y).$$

Exactly $d$ indices receive one additional squaring, namely those with the smallest fractional parts $f_i$.

Thus:

$$k_i= Y-q_i+\varepsilon_i,$$

where $\varepsilon_i\in{0,1}$, and $\varepsilon_i=1$ for the $d$ smallest $f_i$.

2. Modular evaluation

We need

$$\sum_{i=2}^{10000} i^{2^{k_i}} \pmod{1234567891}.$$

Let

$$M=1234567891.$$

Since $M$ is prime and $i<M$,

Fermat gives

$$i^{a\bmod (M-1)} \equiv i^a \pmod M.$$

So we compute:

$$e_i = 2^{k_i}\bmod (M-1),$$

then

$$i^{e_i}\bmod M.$$

3. Python implementation

from math import log, log2, floor

MOD = 1234567891

def solve(n, m):
    data = []
    sum_q = 0

    # Compute c_i = log2(log(i))
    for i in range(2, n + 1):
        c = log2(log(i))

        q = floor(c)
        f = c - q

        data.append((f, i, q))
        sum_q += q

    N = n - 1

    # Largest Y with base_count <= m
    Y = (m + sum_q) // N

    base_count = N * Y - sum_q
    extra = m - base_count

    # Smallest fractional parts get one extra operation
    data.sort()

    extra_indices = set(i for _, i, _ in data[:extra])

    ans = 0

    for f, i, q in data:
        k = Y - q

        if i in extra_indices:
            k += 1

        # exponent = 2^k mod (MOD-1)
        exponent = pow(2, k, MOD - 1)

        ans = (ans + pow(i, exponent, MOD)) % MOD

    return ans

# Checks from the statement
print(solve(5, 3))      # 34
print(solve(10, 100))   # 845339386

# Final problem
print(solve(10_000, 10**16))

4. Code walkthrough

Step 1

For every initial value $i$:

c = log2(log(i))
q = floor(c)
f = c - q

We decompose

$$c_i=q_i+f_i.$$

Step 2

Find the balancing level $Y$:

Y = (m + sum_q) // N

because

$$B(Y)=NY-\sum q_i.$$

Step 3

The remaining

$$d=m-B(Y)$$

operations are assigned to the smallest fractional parts.

data.sort()
extra_indices = set(i for _, i, _ in data[:extra])

Step 4

Compute final squaring count:

k = Y - q
if i in extra_indices:
    k += 1

Step 5

Compute

$$i^{2^{k_i}} \pmod{MOD}$$

using Fermat reduction:

exponent = pow(2, k, MOD - 1)
pow(i, exponent, MOD)

5. Verification

The program reproduces both values from the statement:

$$S(5,3)=34$$

and

$$S(10,100)\equiv 845339386 \pmod{1234567891}.$$

So the derivation and implementation are consistent.

Therefore:

Answer: 950591530

Project Euler Problem 822