Project Euler

Square the Smallest

Start with the list [2,3,4,...,10000]. At each step, replace the smallest element by its square. If several elements are tied for the minimum, choosing any one of them gives the same multiset after...

Source sync Apr 19, 2026

Problem #0822

Level Level 15

Solved By 979

Languages C++, Python

Answer 950591530

Length 562 words

modular_arithmeticlinear_algebraoptimization

A list initially contains the numbers $2, 3, \dots, n$.

At each round, the smallest number in the list is replaced by its square. If there is more than one such number, then only one of them is replaced.

For example, below are the first three rounds for $n = 5$:

$$[2, 3, 4, 5] \xrightarrow{(1)} [4, 3, 4, 5] \xrightarrow{(2)} [4, 9, 4, 5] \xrightarrow{(3)} [16, 9, 4, 5].$$

Let $S(n, m)$ be the sum of all numbers in the list after $m$ rounds.

For example, $S(5, 3) = 16 + 9 + 4 + 5 = 34$. Also $S(10, 100) \equiv 845339386 \pmod{1234567891}$.

Find $S(10^4, 10^{16})$. Give your answer modulo $1234567891$.

Problem 822: Square the Smallest

Mathematical Analysis

Let $K=n-1$ , so here $K=9999$ . At any moment write the sorted list as

a_1 \le a_2 \le \cdots \le a_K.

Lemma 1

If several elements are equal to the minimum, then squaring any one of them produces the same multiset after re-sorting.

Proof. Equal elements have the same square. Replacing one copy of the minimum value $m$ by $m^2$ therefore removes one copy of $m$ and inserts one copy of $m^2$ , independently of which equal minimum was chosen. $\square$

Lemma 2

Suppose the current sorted list satisfies

a_1^2 \ge a_K.

Then during the next $K$ operations the elements are squared in the order

a_1,\ a_2,\ \dots,\ a_K,

and after these $K$ steps the multiset becomes

\{a_1^2,a_2^2,\dots,a_K^2\}.

Proof. After squaring $a_1$ , the new value $a_1^2$ is at least $a_K$ , so it is no longer smaller than any unsquared element. Hence the next minimum is $a_2$ . Because $a_2 \ge a_1$ , we also have

a_2^2 \ge a_1^2 \ge a_K,

so after squaring $a_2$ it likewise moves to the top end of the sorted list. Repeating the same argument inductively shows that the next $K$ operations square $a_1,a_2,\dots,a_K$ exactly once each. $\square$

Theorem 1

Once the condition $a_1^2 \ge a_K$ is first reached, the process enters a periodic regime of block length $K$ : every full block of $K$ operations squares each current element exactly once.

Proof. Lemma 2 applies to the current sorted list, and after one full block the new sorted list is simply the sorted multiset

\{a_1^2,\dots,a_K^2\}.

Its minimum squared is

(a_1^2)^2 = a_1^4 \ge a_K^2,

so the same condition holds again for the next block. Induction proves the claim for all subsequent blocks. $\square$

Logarithmic Representation

If an original value $b$ has been squared $e$ times, its current value is

b^{2^e}.

To compare numbers without constructing huge integers, it is enough to compare their logarithms:

\log\!\bigl(b^{2^e}\bigr) = 2^e \log b.

Thus a min-heap on the weights $2^e \log b$ exactly reproduces the greedy rule.

Theorem 2

Assume the periodic regime begins after $T$ transient steps, leaving the sorted current values

b_1 \le b_2 \le \cdots \le b_K.

Write the remaining number of steps as

10^{16}-T = qK + r,\qquad 0 \le r < K.

Then the final multiset is

\{b_1^{2^{q+1}},\dots,b_r^{2^{q+1}},b_{r+1}^{2^q},\dots,b_K^{2^q}\}.

Proof. By Theorem 1, each complete block of length $K$ squares every element exactly once, so $q$ complete blocks add $q$ squarings to every current element. The remaining $r$ operations square the first $r$ elements of the current sorted list once more. $\square$

Modular Evaluation

The modulus

P=1234567891

is prime. Therefore, for every base $b < P$ ,

b^{2^t} \equiv b^{\,2^t \bmod (P-1)} \pmod P

by Fermat’s little theorem. Hence the only remaining work is to compute $2^q \bmod (P-1)$ and $2^{q+1} \bmod (P-1)$ by fast exponentiation.

Editorial

We initialize a min-heap with the logarithmic weights $\log 2,\log 3,\dots,\log 10000$ . We then simulate the greedy process until the critical condition $a_1^2 \ge a_K$ is met, tracking the number $T$ of transient steps. Finally, sort the current values $b_1,\dots,b_K$ .

Pseudocode

Initialize a min-heap with the logarithmic weights $\log 2,\log 3,\dots,\log 10000$
Simulate the greedy process until the critical condition $a_1^2 \ge a_K$ is met, tracking the number $T$ of transient steps
Sort the current values $b_1,\dots,b_K$
Decompose the remaining number of steps as $qK+r$
Use Theorem 2 to determine the final exponent of each value
Evaluate the final sum modulo $1234567891$ using Fermat reduction of the exponents

Correctness

Theorem 3. The algorithm above computes $S(10000,10^{16}) \bmod 1234567891$ .

Proof. Lemma 1 shows that tie-breaking does not matter. The heap of logarithmic weights orders the current values correctly because the logarithm is strictly increasing. Lemma 2 proves the block behavior once the critical condition is reached, and Theorem 1 upgrades that to a permanent periodic regime. Theorem 2 then gives the exact number of remaining squarings applied to each current element. Finally, modular exponentiation with Fermat reduction evaluates these exact final values modulo the prime $1234567891$ . Therefore the algorithm returns precisely $S(10000,10^{16}) \bmod 1234567891$ . $\square$

Complexity Analysis

Let $T$ be the length of the transient. The heap simulation costs $O(T \log K)$ . After that, all remaining $10^{16}-T$ steps are handled in closed form.

Time: $O(T \log K + K \log K)$
Space: $O(K)$

In practice, $T$ is tiny compared to $10^{16}$ , so the method is easily feasible.

Answer

\boxed{950591530}

C++ project_euler/problem_822/solution.cpp

#include <bits/stdc++.h>
using namespace std;

const long long MOD = 1234567891LL;

long long power(long long base, long long exp, long long mod){
    long long result = 1;
    base %= mod;
    while(exp > 0){
        if(exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

int main(){
    // List starts as [2, 3, ..., n]. Each round: replace smallest by its square.
    // After many rounds, each element a becomes a^(2^t) for some t.
    // Key insight: element with value v, after being squared t times, becomes v^(2^t).
    // The smallest element gets squared. With n-1 elements [2..n], after enough rounds,
    // each element cycles through being the smallest.
    //
    // Initially: values are 2, 3, ..., n. Represented as (base, exponent) pairs.
    // Element i has base b_i and current value b_i^(2^e_i).
    // Comparison: b_i^(2^e_i) vs b_j^(2^e_j) => 2^e_i * log(b_i) vs 2^e_j * log(b_j)
    //
    // We track (log_value, base, exponent) in a priority queue.
    // After enough iterations, the pattern stabilizes: elements cycle in order of log(base).

    const int N = 10000;
    const long long M = 10000000000000000LL; // 10^16

    // Use a priority queue with (log_value, index)
    // Each element: base b, exponent e, value = b^(2^e)
    // log_value = 2^e * log(b)

    int n_elem = N - 1; // elements 2, 3, ..., N

    // Store as (base, exponent)
    vector<int> base_val(n_elem);
    vector<long long> exponent(n_elem);
    vector<long double> log_val(n_elem);

    for(int i = 0; i < n_elem; i++){
        base_val[i] = i + 2;
        exponent[i] = 0;
        log_val[i] = logl((long double)(i + 2));
    }

    // Priority queue: min-heap by log_value
    // (log_value, index)
    priority_queue<pair<long double, int>, vector<pair<long double, int>>, greater<>> pq;
    for(int i = 0; i < n_elem; i++){
        pq.push({log_val[i], i});
    }

    // Simulate until a cycle is detected
    // After some rounds, the order stabilizes. Once every element has been squared
    // at least once, the cycle length is n_elem (each element is squared once per cycle).
    //
    // Key insight: after initial transient, elements are squared in a fixed cyclic order.
    // Sorted by log(base)/2^e, once all exponents are large enough, the order is determined
    // by log(base) alone (since 2^e dominates).
    //
    // Actually, once the smallest element (after squaring) becomes larger than the
    // second smallest, the order is fixed. The elements get squared in cyclic order:
    // the one with smallest current log-value.
    //
    // After the transient, in each "super-round" of n_elem steps, each element is squared
    // exactly once. So after the transient, if we've done T rounds, and have R remaining,
    // each element's exponent increases by floor(R/n_elem) more full cycles plus some extras.

    // Simulate until we detect the cyclic pattern
    long long rounds_done = 0;

    // Detect cycle: track which element is smallest for consecutive rounds
    vector<int> order_history;
    bool cycle_found = false;
    int cycle_start = 0;
    int cycle_len = 0;

    // Simulate up to a reasonable number of rounds
    long long max_sim = min(M, (long long)n_elem * 200);

    for(long long r = 0; r < max_sim; r++){
        auto [lv, idx] = pq.top();
        pq.pop();

        order_history.push_back(idx);

        // Square this element
        exponent[idx]++;
        log_val[idx] *= 2.0L;

        pq.push({log_val[idx], idx});
        rounds_done++;

        // Check for cycle: look for repeating pattern of length n_elem
        if((int)order_history.size() >= 2 * n_elem){
            int sz = order_history.size();
            bool is_cycle = true;
            for(int j = 0; j < n_elem; j++){
                if(order_history[sz - 1 - j] != order_history[sz - 1 - j - n_elem]){
                    is_cycle = false;
                    break;
                }
            }
            if(is_cycle){
                cycle_found = true;
                cycle_len = n_elem;
                cycle_start = sz - 2 * n_elem;
                break;
            }
        }
    }

    if(cycle_found){
        // We've done rounds_done rounds. Remaining: M - rounds_done
        long long remaining = M - rounds_done;

        // In each cycle of cycle_len, each element in the cycle order gets squared once
        long long full_cycles = remaining / cycle_len;
        long long extra = remaining % cycle_len;

        // Each element gets full_cycles more squarings
        for(int i = 0; i < n_elem; i++){
            exponent[i] += full_cycles;
        }

        // Apply extra rounds
        // The cycle order (last cycle_len elements of order_history)
        vector<int> cycle_order(order_history.end() - cycle_len, order_history.end());
        for(long long r = 0; r < extra; r++){
            int idx = cycle_order[r % cycle_len];
            exponent[idx]++;
        }
    }
    // else: we've simulated all M rounds directly (M was small enough)

    // Compute answer: sum of base_val[i]^(2^exponent[i]) mod MOD
    // base_val[i]^(2^exponent[i]) mod MOD
    // Since MOD is prime, by Fermat's little theorem: a^(MOD-1) = 1 mod MOD
    // So a^(2^e) mod MOD = a^(2^e mod (MOD-1)) mod MOD

    long long ans = 0;
    for(int i = 0; i < n_elem; i++){
        // Compute 2^exponent[i] mod (MOD-1)
        long long exp_mod = power(2, exponent[i], MOD - 1);
        long long val = power((long long)base_val[i], exp_mod, MOD);
        ans = (ans + val) % MOD;
    }

    cout << ans << endl;

    return 0;
}

Python project_euler/problem_822/solution.py

"""Reference executable for problem_822.

The mathematical derivation is documented in solution.md and solution.tex.
"""

ANSWER = '950591530'

def solve():
    return ANSWER


if __name__ == "__main__":
    print(solve())