Problem 823: Factor Shuffle

Project Euler Problem 823 Solution

EulerSolve provides an optimized solution for Project Euler Problem 823, Factor Shuffle, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For every integer $j\in\{2,\dots,n\}$, write its prime factors in nondecreasing order and treat that ordered list as one pile. One round removes the first factor from every nonempty pile, sorts the extracted factors, and appends them as a new pile. After $m$ rounds, $S(n,m)$ is the sum of the products of the factors still remaining in the active piles, taken modulo $1234567891$. The challenge is that the target round count is $10^{16}$, so the solution must replace long simulation by a structural description of the stabilized shuffle. Mathematical Approach The implementations work with prime-factor lists directly. The key is to separate two layers of the process: the evolution of pile lengths, and the evolution of the sorted factors extracted at each round. Step 1: Model the piles and track the conserved quantity Let the initial pile for $j$ be the ascending prime-factor list of $j$. If $\Omega(j)$ denotes the total number of prime factors of $j$, counted with multiplicity, then the total number of stored factors is $$T=\sum_{j=2}^{n}\Omega(j).$$ This quantity never changes. In every round we remove exactly one factor from each active pile, then we create one new pile whose length is exactly the number of active piles in that round. Therefore only the distribution of lengths changes, not the total factor count....

Detailed mathematical approach

Problem Summary

For every integer $j\in\{2,\dots,n\}$, write its prime factors in nondecreasing order and treat that ordered list as one pile. One round removes the first factor from every nonempty pile, sorts the extracted factors, and appends them as a new pile. After $m$ rounds, $S(n,m)$ is the sum of the products of the factors still remaining in the active piles, taken modulo $1234567891$. The challenge is that the target round count is $10^{16}$, so the solution must replace long simulation by a structural description of the stabilized shuffle.

Mathematical Approach

The implementations work with prime-factor lists directly. The key is to separate two layers of the process: the evolution of pile lengths, and the evolution of the sorted factors extracted at each round.

Step 1: Model the piles and track the conserved quantity

Let the initial pile for $j$ be the ascending prime-factor list of $j$. If $\Omega(j)$ denotes the total number of prime factors of $j$, counted with multiplicity, then the total number of stored factors is

$$T=\sum_{j=2}^{n}\Omega(j).$$

This quantity never changes. In every round we remove exactly one factor from each active pile, then we create one new pile whose length is exactly the number of active piles in that round. Therefore only the distribution of lengths changes, not the total factor count.

If the current pile lengths are $\lambda_1,\dots,\lambda_k$, then the next round transforms them into the positive parts of

$$(\lambda_1-1,\dots,\lambda_k-1,k).$$

This is the same length update that appears in Bulgarian solitaire. It explains why the number of relevant pile positions is only on the order of $\sqrt{T}$: after the transient, the lengths cluster around a staircase shape rather than staying spread across all $T$ factors.

Step 2: Encode each round by a sorted extracted row

At round $t$, let the extracted factors after sorting be

$$E_t=\bigl(a_1(t),a_2(t),\dots,a_{k_t}(t)\bigr),\qquad a_1(t)\le a_2(t)\le \dots \le a_{k_t}(t),$$

where $k_t$ is the number of active piles at that round. For bookkeeping, extend the definition by setting

$$a_k(t)=1,\quad k>k_t.$$

This padding is harmless because multiplying by $1$ does not change a pile product.

The new pile created at round $t$ is exactly the list $E_t$. One round later its first entry $a_1(t)$ is removed; two rounds later its second entry $a_2(t)$ is removed; in general, its $k$-th entry is extracted exactly $k$ rounds after the pile is born.

Step 3: Why periodic rows appear

Once the length dynamics has settled, the shuffle becomes shift-invariant: the same geometric position inside a pile is revisited after a fixed number of rounds. Because the $k$-th entry of a newborn pile survives for exactly $k$ rounds, the $k$-th extracted row naturally repeats with period $k$.

So after a transient time $t_0$, the implementations use the periodic description

$$a_k(t+k)=a_k(t)\qquad (t\ge t_0),$$

or equivalently

$$a_k(t)=p_k(t\bmod k),$$

where $p_k$ is the stored cycle for row $k$. The code does not assume these cycles in advance; it simulates until the last $k$ values of row $k$ repeat consistently, then records that row as periodic.

Step 4: Reconstruct the remaining pile products along diagonals

Consider the state after round $m$. A pile created $d$ rounds earlier was born at round $m-1-d$. By then it has already lost its first $d$ entries, so its remaining product is

$$\prod_{k=d+1}^{K} a_k(m-1-d),$$

where $K$ is the largest row that is ever nontrivial in the stabilized regime. The pile exists only if it originally had at least $d+1$ entries, which is equivalent to

$$a_{d+1}(m-1-d)\ne 1.$$

Therefore the total sum is a sum of diagonal suffix-products through the row table:

$$S(n,m)=\sum_{d=0}^{K-1}\mathbf{1}\!\bigl(a_{d+1}(m-1-d)\ne 1\bigr)\prod_{k=d+1}^{K} a_k(m-1-d)\pmod{1234567891}.$$

This is the decisive transformation: instead of advancing one round at a time, we can jump directly to the large target round by evaluating finitely many periodic rows.

Step 5: Reduce a huge round number to modular row indices

After the transient, only residue classes inside each row period matter. If $r_0=m-t_0-1$, then the large-round formula becomes

$$S(n,m)=\sum_{d=0}^{K-1}\mathbf{1}\!\bigl(p_{d+1}((r_0-d)\bmod(d+1))\ne 1\bigr)\prod_{k=d+1}^{K} p_k((r_0-d)\bmod k)\pmod{1234567891}.$$

All dependence on the enormous value of $m$ is now reduced to modular indexing inside short stored cycles.

Worked Example: $S(5,3)=21$

Start from the four piles coming from $2,3,4,5$:

$$[2],\ [3],\ [2,2],\ [5].$$

Round 1 extracts $[2,3,2,5]$, which sorts to $[2,2,3,5]$. The piles after the round are $[2]$ and $[2,2,3,5]$, so the sum of residual products is

$$2+(2\cdot 2\cdot 3\cdot 5)=62.$$

Round 2 extracts $[2,2]$. The piles become $[2,3,5]$ and $[2,2]$, giving

$$2\cdot 3\cdot 5+2\cdot 2=30+4=34.$$

Round 3 again extracts $[2,2]$. The remaining piles are $[3,5]$, $[2]$, and $[2,2]$, so

$$S(5,3)=3\cdot 5+2+2\cdot 2=15+2+4=21.$$

This is the first small checkpoint used by the implementations.

How the Code Works

The C++, Python, and Java implementations first build a smallest-prime-factor sieve up to $n$, then factor every integer $2,3,\dots,n$ into an ascending prime list. That produces the initial piles efficiently.

They then simulate the shuffle while monitoring the $k$-th extracted value of each round. For every row $k$, the implementation stores the last $k$ observed values and checks whether the current value matches the one from $k$ rounds earlier; missing entries are treated as $1$. Once all monitored rows have repeated long enough, the row cycles are recorded.

If the requested round still lies inside the transient, the answer is obtained by direct simulation. Otherwise the implementation evaluates the diagonal formula above, using only modular indexing into the stored cycles. Small checkpoints such as $S(5,3)=21$ and $S(10,100)=257$ verify the logic before the final target value is produced.

Complexity Analysis

Let $T=\sum_{j=2}^{n}\Omega(j)$, and let $K$ be the highest periodic row that contains a value other than $1$. Building the smallest-prime-factor sieve takes $O(n\log\log n)$ time and $O(n)$ memory, while constructing the initial piles takes $O(T)$ total factor extractions. If periodicity is detected after $R$ simulated rounds and $k_t$ piles are active at round $t$, the transient costs $O\!\left(\sum_{t=1}^{R} k_t\log k_t\right)$ because each round sorts the extracted factors. After that, a huge-round query is answered in $O(K^2)$ time using the diagonal product formula, with $O(K^2)$ memory for the stored cycles.

Footnotes and References

Problem page: https://projecteuler.net/problem=823
Prime factorization: Wikipedia — Prime factor
Bulgarian solitaire: Wikipedia — Bulgarian solitaire
Integer partitions: Wikipedia — Partition (number theory)
Modular arithmetic: Wikipedia — Modular arithmetic

Problem 823 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <deque>
#include <iostream>
#include <vector>

using i64 = std::int64_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

static constexpr i64 kMod = 1'234'567'891;

struct SimResult {
    i64 end_t;
    std::vector<std::vector<int>> patterns;
    int kmax;
};

static std::vector<int> sieve_spf(int n) {
    std::vector<int> spf(n + 1);
    for (int i = 0; i <= n; ++i) {
        spf[i] = i;
    }
    int lim = static_cast<int>(std::sqrt(static_cast<long double>(n)));
    for (int i = 2; i <= lim; ++i) {
        if (spf[i] == i) {
            for (int j = i * i; j <= n; j += i) {
                if (spf[j] == j) {
                    spf[j] = i;
                }
            }
        }
    }
    return spf;
}

static std::vector<int> factor_list(int x, const std::vector<int>& spf) {
    std::vector<int> out;
    while (x > 1) {
        int p = spf[x];
        out.push_back(p);
        x /= p;
    }
    return out;
}

static i64 direct_sum_mod(int n, i64 m, i64 mod) {
    const auto spf = sieve_spf(n);

    std::vector<std::vector<int>> piles;
    std::vector<int> pos;
    piles.reserve(n);
    pos.reserve(n);
    for (int i = 2; i <= n; ++i) {
        piles.push_back(factor_list(i, spf));
        pos.push_back(0);
    }

    for (i64 t = 0; t < m; ++t) {
        const int k = static_cast<int>(piles.size());
        std::vector<int> extracted(k);
        std::vector<std::vector<int>> next_piles;
        std::vector<int> next_pos;
        next_piles.reserve(k + 1);
        next_pos.reserve(k + 1);

        for (int i = 0; i < k; ++i) {
            extracted[i] = piles[i][pos[i]];
            ++pos[i];
            if (pos[i] < static_cast<int>(piles[i].size())) {
                next_piles.push_back(std::move(piles[i]));
                next_pos.push_back(pos[i]);
            }
        }

        std::sort(extracted.begin(), extracted.end());
        next_piles.push_back(std::move(extracted));
        next_pos.push_back(0);

        piles.swap(next_piles);
        pos.swap(next_pos);
    }

    i64 total = 0;
    for (int i = 0; i < static_cast<int>(piles.size()); ++i) {
        i64 prod = 1;
        for (int j = pos[i]; j < static_cast<int>(piles[i].size()); ++j) {
            prod = static_cast<i64>((static_cast<u128>(prod) * static_cast<u64>(piles[i][j])) % mod);
        }
        total += prod;
        if (total >= mod) {
            total -= mod;
        }
    }
    return total;
}

static SimResult simulate_until_periodic(
    int n,
    int k_extra = 10,
    int streak_needed = 2000,
    int max_rounds = 200000
) {
    const auto spf = sieve_spf(n);

    std::vector<std::vector<int>> piles;
    std::vector<int> pos;
    piles.reserve(n);
    pos.reserve(n);

    i64 total_factors = 0;
    for (int i = 2; i <= n; ++i) {
        auto f = factor_list(i, spf);
        total_factors += static_cast<i64>(f.size());
        piles.push_back(std::move(f));
        pos.push_back(0);
    }

    const int k_lim = static_cast<int>(std::sqrt(2.0L * static_cast<long double>(total_factors))) + k_extra;

    std::vector<std::deque<int>> bufs(k_lim + 1);

    int stable_streak = 0;

    for (int t = 1; t <= max_rounds; ++t) {
        const int k = static_cast<int>(piles.size());
        std::vector<int> extracted(k);

        std::vector<std::vector<int>> next_piles;
        std::vector<int> next_pos;
        next_piles.reserve(k + 1);
        next_pos.reserve(k + 1);

        for (int i = 0; i < k; ++i) {
            extracted[i] = piles[i][pos[i]];
            ++pos[i];
            if (pos[i] < static_cast<int>(piles[i].size())) {
                next_piles.push_back(std::move(piles[i]));
                next_pos.push_back(pos[i]);
            }
        }

        std::sort(extracted.begin(), extracted.end());
        next_piles.push_back(extracted);
        next_pos.push_back(0);

        piles.swap(next_piles);
        pos.swap(next_pos);

        int mm = std::min(static_cast<int>(extracted.size()), k_lim);

        if (t <= k_lim) {
            for (int kk = 1; kk <= mm; ++kk) {
                auto& b = bufs[kk];
                b.push_back(extracted[kk - 1]);
                if (static_cast<int>(b.size()) > kk) {
                    b.pop_front();
                }
            }
            for (int kk = mm + 1; kk <= k_lim; ++kk) {
                auto& b = bufs[kk];
                b.push_back(1);
                if (static_cast<int>(b.size()) > kk) {
                    b.pop_front();
                }
            }
            stable_streak = 0;
            continue;
        }

        bool all_ok = true;

        for (int kk = 1; kk <= mm; ++kk) {
            auto& b = bufs[kk];
            const int v = extracted[kk - 1];
            if (b.front() != v) {
                all_ok = false;
            }
            b.push_back(v);
            if (static_cast<int>(b.size()) > kk) {
                b.pop_front();
            }
        }

        for (int kk = mm + 1; kk <= k_lim; ++kk) {
            auto& b = bufs[kk];
            if (b.front() != 1) {
                all_ok = false;
            }
            b.push_back(1);
            if (static_cast<int>(b.size()) > kk) {
                b.pop_front();
            }
        }

        if (all_ok) {
            ++stable_streak;
            if (stable_streak >= streak_needed) {
                std::vector<std::vector<int>> patterns(k_lim + 1);
                int kmax = 0;
                for (int kk = 1; kk <= k_lim; ++kk) {
                    patterns[kk] = std::vector<int>(bufs[kk].begin(), bufs[kk].end());
                    bool has_non_one = false;
                    for (int x : patterns[kk]) {
                        if (x != 1) {
                            has_non_one = true;
                            break;
                        }
                    }
                    if (has_non_one) {
                        kmax = kk;
                    }
                }
                return {t, std::move(patterns), kmax};
            }
        } else {
            stable_streak = 0;
        }
    }

    throw std::runtime_error("periodicity not detected");
}

static i64 sum_at_round_mod(int n, i64 m, i64 mod) {
    const SimResult sim = simulate_until_periodic(n);

    if (m <= sim.end_t) {
        return direct_sum_mod(n, m, mod);
    }

    const i64 r0 = m - sim.end_t - 1;
    i64 total = 0;

    for (int d = 0; d < sim.kmax; ++d) {
        const i64 r = r0 - d;
        if (sim.patterns[d + 1][static_cast<int>(r % (d + 1))] == 1) {
            continue;
        }

        i64 prod = 1;
        for (int k = sim.kmax; k > d; --k) {
            const int v = sim.patterns[k][static_cast<int>(r % k)];
            if (v != 1) {
                prod = static_cast<i64>((static_cast<u128>(prod) * static_cast<u64>(v)) % mod);
            }
        }

        total += prod;
        if (total >= mod) {
            total -= mod;
        }
    }

    return total;
}

int main() {
    if (direct_sum_mod(5, 3, (1LL << 62)) != 21) {
        std::cerr << "Validation failed for S(5,3).\n";
        return 1;
    }
    if (direct_sum_mod(10, 100, (1LL << 62)) != 257) {
        std::cerr << "Validation failed for S(10,100).\n";
        return 1;
    }

    std::cout << sum_at_round_mod(10'000, 10'000'000'000'000'000LL, kMod) << '\n';
    return 0;
}

Python

import math
from collections import deque

kMod = 1234567891

def sieve_spf(n):
    spf = list(range(n + 1))
    lim = math.isqrt(n)
    for i in range(2, lim + 1):
        if spf[i] == i:
            for j in range(i * i, n + 1, i):
                if spf[j] == j:
                    spf[j] = i
    return spf

def factor_list(x, spf):
    out = []
    while x > 1:
        p = spf[x]
        out.append(p)
        x //= p
    return out

def direct_sum_mod(n, m, mod):
    spf = sieve_spf(n)
    piles = []
    pos = []
    for i in range(2, n + 1):
        piles.append(factor_list(i, spf))
        pos.append(0)
        
    for _ in range(m):
        k = len(piles)
        extracted = [0] * k
        next_piles = []
        next_pos = []
        
        for i in range(k):
            extracted[i] = piles[i][pos[i]]
            pos[i] += 1
            if pos[i] < len(piles[i]):
                next_piles.append(piles[i])
                next_pos.append(pos[i])
                
        extracted.sort()
        next_piles.append(extracted)
        next_pos.append(0)
        
        piles = next_piles
        pos = next_pos
        
    total = 0
    for i in range(len(piles)):
        prod = 1
        for j in range(pos[i], len(piles[i])):
            prod = (prod * piles[i][j]) % mod
        total = (total + prod) % mod
        
    return total

def simulate_until_periodic(n, k_extra=10, streak_needed=2000, max_rounds=200000):
    spf = sieve_spf(n)
    piles = []
    pos = []
    total_factors = 0
    
    for i in range(2, n + 1):
        f = factor_list(i, spf)
        total_factors += len(f)
        piles.append(f)
        pos.append(0)
        
    k_lim = int(math.sqrt(2.0 * total_factors)) + k_extra
    bufs = [deque() for _ in range(k_lim + 1)]
    stable_streak = 0
    
    for t in range(1, max_rounds + 1):
        k = len(piles)
        extracted = [0] * k
        next_piles = []
        next_pos = []
        
        for i in range(k):
            extracted[i] = piles[i][pos[i]]
            pos[i] += 1
            if pos[i] < len(piles[i]):
                next_piles.append(piles[i])
                next_pos.append(pos[i])
                
        extracted.sort()
        next_piles.append(extracted)
        next_pos.append(0)
        
        piles = next_piles
        pos = next_pos
        
        mm = min(len(extracted), k_lim)
        
        if t <= k_lim:
            for kk in range(1, mm + 1):
                b = bufs[kk]
                b.append(extracted[kk - 1])
                if len(b) > kk:
                    b.popleft()
            for kk in range(mm + 1, k_lim + 1):
                b = bufs[kk]
                b.append(1)
                if len(b) > kk:
                    b.popleft()
            stable_streak = 0
            continue
            
        all_ok = True
        for kk in range(1, mm + 1):
            b = bufs[kk]
            v = extracted[kk - 1]
            if b[0] != v:
                all_ok = False
            b.append(v)
            if len(b) > kk:
                b.popleft()
                
        for kk in range(mm + 1, k_lim + 1):
            b = bufs[kk]
            if b[0] != 1:
                all_ok = False
            b.append(1)
            if len(b) > kk:
                b.popleft()
                
        if all_ok:
            stable_streak += 1
            if stable_streak >= streak_needed:
                patterns = [[] for _ in range(k_lim + 1)]
                kmax = 0
                for kk in range(1, k_lim + 1):
                    patterns[kk] = list(bufs[kk])
                    if any(x != 1 for x in patterns[kk]):
                        kmax = kk
                return t, patterns, kmax
        else:
            stable_streak = 0
            
    raise RuntimeError("periodicity not detected")

def sum_at_round_mod(n, m, mod):
    end_t, patterns, kmax = simulate_until_periodic(n)
    
    if m <= end_t:
        return direct_sum_mod(n, m, mod)
        
    r0 = m - end_t - 1
    total = 0
    
    for d in range(kmax):
        r = r0 - d
        if patterns[d + 1][r % (d + 1)] == 1:
            continue
            
        prod = 1
        for k in range(kmax, d, -1):
            v = patterns[k][r % k]
            if v != 1:
                prod = (prod * v) % mod
                
        total = (total + prod) % mod
        
    return total

def solve():
    ans = sum_at_round_mod(10000, 10000000000000000, kMod)
    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.Collections;
import java.util.Deque;
import java.util.LinkedList;

public class Euler823 {

    static final long kMod = 1234567891L;

    static int[] sieveSpf(int n) {
        int[] spf = new int[n + 1];
        for (int i = 0; i <= n; ++i) {
            spf[i] = i;
        }
        int lim = (int) Math.sqrt(n);
        for (int i = 2; i <= lim; ++i) {
            if (spf[i] == i) {
                for (int j = i * i; j <= n; j += i) {
                    if (spf[j] == j) {
                        spf[j] = i;
                    }
                }
            }
        }
        return spf;
    }

    static ArrayList<Integer> factorList(int x, int[] spf) {
        ArrayList<Integer> out = new ArrayList<>();
        while (x > 1) {
            int p = spf[x];
            out.add(p);
            x /= p;
        }
        return out;
    }

    static class SimResult {
        long endT;
        ArrayList<ArrayList<Integer>> patterns;
        int kmax;

        SimResult(long endT, ArrayList<ArrayList<Integer>> patterns, int kmax) {
            this.endT = endT;
            this.patterns = patterns;
            this.kmax = kmax;
        }
    }

    static long directSumMod(int n, long m, long mod) {
        int[] spf = sieveSpf(n);
        ArrayList<ArrayList<Integer>> piles = new ArrayList<>(n);
        ArrayList<Integer> pos = new ArrayList<>(n);

        for (int i = 2; i <= n; ++i) {
            piles.add(factorList(i, spf));
            pos.add(0);
        }

        for (long t = 0; t < m; ++t) {
            int k = piles.size();
            ArrayList<Integer> extracted = new ArrayList<>(k);
            ArrayList<ArrayList<Integer>> nextPiles = new ArrayList<>(k + 1);
            ArrayList<Integer> nextPos = new ArrayList<>(k + 1);

            for (int i = 0; i < k; ++i) {
                int pIdx = pos.get(i);
                extracted.add(piles.get(i).get(pIdx));
                pos.set(i, pIdx + 1);
                if (pos.get(i) < piles.get(i).size()) {
                    nextPiles.add(piles.get(i));
                    nextPos.add(pos.get(i));
                }
            }

            Collections.sort(extracted);
            nextPiles.add(extracted);
            nextPos.add(0);

            piles = nextPiles;
            pos = nextPos;
        }

        long total = 0;
        for (int i = 0; i < piles.size(); ++i) {
            long prod = 1;
            for (int j = pos.get(i); j < piles.get(i).size(); ++j) {
                prod = (prod * piles.get(i).get(j)) % mod;
            }
            total += prod;
            if (total >= mod) {
                total -= mod;
            }
        }
        return total;
    }

    static SimResult simulateUntilPeriodic(int n, int kExtra, int streakNeeded, int maxRounds) {
        int[] spf = sieveSpf(n);

        ArrayList<ArrayList<Integer>> piles = new ArrayList<>(n);
        ArrayList<Integer> pos = new ArrayList<>(n);

        long totalFactors = 0;
        for (int i = 2; i <= n; ++i) {
            ArrayList<Integer> f = factorList(i, spf);
            totalFactors += f.size();
            piles.add(f);
            pos.add(0);
        }

        int kLim = (int) Math.sqrt(2.0 * totalFactors) + kExtra;
        ArrayList<Deque<Integer>> bufs = new ArrayList<>(kLim + 1);
        for (int i = 0; i <= kLim; i++) {
            bufs.add(new LinkedList<>());
        }

        int stableStreak = 0;

        for (int t = 1; t <= maxRounds; ++t) {
            int k = piles.size();
            ArrayList<Integer> extracted = new ArrayList<>(k);
            ArrayList<ArrayList<Integer>> nextPiles = new ArrayList<>(k + 1);
            ArrayList<Integer> nextPos = new ArrayList<>(k + 1);

            for (int i = 0; i < k; ++i) {
                int pIdx = pos.get(i);
                extracted.add(piles.get(i).get(pIdx));
                pos.set(i, pIdx + 1);
                if (pos.get(i) < piles.get(i).size()) {
                    nextPiles.add(piles.get(i));
                    nextPos.add(pos.get(i));
                }
            }

            Collections.sort(extracted);
            nextPiles.add(extracted);
            nextPos.add(0);

            piles = nextPiles;
            pos = nextPos;

            int mm = Math.min(extracted.size(), kLim);

            if (t <= kLim) {
                for (int kk = 1; kk <= mm; ++kk) {
                    Deque<Integer> b = bufs.get(kk);
                    b.addLast(extracted.get(kk - 1));
                    if (b.size() > kk)
                        b.removeFirst();
                }
                for (int kk = mm + 1; kk <= kLim; ++kk) {
                    Deque<Integer> b = bufs.get(kk);
                    b.addLast(1);
                    if (b.size() > kk)
                        b.removeFirst();
                }
                stableStreak = 0;
                continue;
            }

            boolean allOk = true;

            for (int kk = 1; kk <= mm; ++kk) {
                Deque<Integer> b = bufs.get(kk);
                int v = extracted.get(kk - 1);
                if (b.peekFirst() != v) {
                    allOk = false;
                }
                b.addLast(v);
                if (b.size() > kk)
                    b.removeFirst();
            }

            for (int kk = mm + 1; kk <= kLim; ++kk) {
                Deque<Integer> b = bufs.get(kk);
                if (b.peekFirst() != 1) {
                    allOk = false;
                }
                b.addLast(1);
                if (b.size() > kk)
                    b.removeFirst();
            }

            if (allOk) {
                stableStreak++;
                if (stableStreak >= streakNeeded) {
                    ArrayList<ArrayList<Integer>> patterns = new ArrayList<>(kLim + 1);
                    for (int i = 0; i <= kLim; i++)
                        patterns.add(new ArrayList<>());

                    int kmax = 0;
                    for (int kk = 1; kk <= kLim; ++kk) {
                        ArrayList<Integer> pat = new ArrayList<>(bufs.get(kk));
                        patterns.set(kk, pat);
                        boolean hasNonOne = false;
                        for (int x : pat) {
                            if (x != 1) {
                                hasNonOne = true;
                                break;
                            }
                        }
                        if (hasNonOne) {
                            kmax = kk;
                        }
                    }
                    return new SimResult(t, patterns, kmax);
                }
            } else {
                stableStreak = 0;
            }
        }

        throw new RuntimeException("periodicity not detected");
    }

    static long sumAtRoundMod(int n, long m, long mod) {
        SimResult sim = simulateUntilPeriodic(n, 10, 2000, 200000);

        if (m <= sim.endT) {
            return directSumMod(n, m, mod);
        }

        long r0 = m - sim.endT - 1;
        long total = 0;

        for (int d = 0; d < sim.kmax; ++d) {
            long r = r0 - d;
            int patIndex = (int) (r % (d + 1));
            if (patIndex < 0)
                patIndex += (d + 1); // just in case

            if (sim.patterns.get(d + 1).get(patIndex) == 1) {
                continue;
            }

            long prod = 1;
            for (int k = sim.kmax; k > d; --k) {
                int innerIndex = (int) (r % k);
                if (innerIndex < 0)
                    innerIndex += k;
                int v = sim.patterns.get(k).get(innerIndex);
                if (v != 1) {
                    prod = (prod * v) % mod;
                }
            }

            total += prod;
            if (total >= mod) {
                total -= mod;
            }
        }

        return total;
    }

    public static String solve() {
        return Long.toString(sumAtRoundMod(10000, 10000000000000000L, kMod));
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}