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

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