Problem 326: Modulo Summations

Project Euler Problem 326 Solution

EulerSolve provides an optimized solution for Project Euler Problem 326, Modulo Summations, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The sequence is defined by $$a_1=1,\qquad a_n\equiv \sum_{k=1}^{n-1}k\,a_k\pmod n\qquad(n\ge 2).$$ For given $N$ and $M$, we must count how many subarrays of $$a_1,a_2,\dots,a_N$$ have sum divisible by $M$. The final instance is $$f(10^{12},10^6).$$ A direct $O(N)$ pass is impossible, so the key is to expose the hidden periodic structure. Mathematical Approach 1) Prefix sums turn divisibility into equal residues. Define prefix sums modulo $M$: $$P_t\equiv \sum_{i=1}^{t}a_i\pmod M,\qquad P_0=0.$$ Then for any $0\le i \lt j\le N$, $$\sum_{k=i+1}^{j}a_k\equiv 0\pmod M \iff P_i\equiv P_j\pmod M.$$ So if residue $r$ appears $c_r$ times among $$P_0,P_1,\dots,P_N,$$ then it contributes $$\binom{c_r}{2}$$ valid subarrays. Therefore $$f(N,M)=\sum_{r=0}^{M-1}\binom{c_r}{2}.$$ 2) Introduce the weighted running sum. Let $$W_n=\sum_{k=1}^{n}k\,a_k.$$ Then the recurrence is simply $$a_n=W_{n-1}\bmod n,\qquad W_n=W_{n-1}+n\,a_n.$$ This is exactly the pair of variables tracked in the C++ code. 3) The first terms reveal a 6-step block pattern. The code checks the first ten values $$1,1,0,3,0,3,5,4,1,9,$$ and if we continue a little further, the sequence groups naturally into blocks of length 6: $$[1,1,0,3,0,3],$$ $$[5,4,1,9,1,6],$$ $$[9,7,2,15,2,9],\dots$$ This is not an accident....

Detailed mathematical approach

Problem Summary

A sequence is built out of its own weighted running total. It starts with $a_1=1$, and for $n\ge2$

$$a_n=\Bigl(\sum_{k=1}^{n-1}k\,a_k\Bigr)\bmod n,$$

so its first ten terms are $1,1,0,3,0,3,5,4,1,9$. For positive integers $N$ and $M$, let $f(N,M)$ be the number of pairs $(p,q)$ with $1\le p\le q\le N$ whose block of terms has a sum divisible by $M$, that is

$$\Bigl(\sum_{i=p}^{q}a_i\Bigr)\bmod M=0.$$

The statement gives $f(10,10)=4$, realised by the four pairs $(3,3)$, $(5,5)$, $(7,9)$ and $(9,10)$, and $f(10^{4},10^{3})=97158$. We must find $f(10^{12},10^{6})$. ^[1]

The whole solution in one sentence: a block sum is divisible by $M$ exactly when two running totals leave the same remainder modulo $M$, and those remainders turn out to repeat every $6M$ steps, so one short pass records how often each remainder occurs and the pairs are then counted group by group.

Mathematical Approach

1. From blocks to running totals

Let $S_0=0$ and $S_n=a_1+a_2+\cdots+a_n$, the running totals of the sequence. ^[6] Subtracting two of them cancels everything outside a block, so for $1\le p\le q$

$$\sum_{i=p}^{q}a_i=S_q-S_{p-1}.$$

A difference is a multiple of $M$ exactly when its two parts are congruent modulo $M$, ^[4] so the condition on $(p,q)$ is

$$S_q\equiv S_{p-1}\pmod M.$$

The map $(p,q)\mapsto(p-1,q)$ sends the pairs with $1\le p\le q\le N$ one to one onto the pairs $(i,j)$ with $0\le i\lt j\le N$. So $f(N,M)$ counts the two-element subsets $\{i,j\}$ of $\{0,1,\dots,N\}$ whose running totals are congruent. Group those $N+1$ indices by that remainder: for $0\le r\lt M$ put

$$c_r=\#\{\,n\le N:S_n\bmod M=r\,\},$$

where $n$ runs over $0,1,\dots,N$. A group of size $c$ contains $\binom{c}{2}$ two-element subsets, and those are exactly the subsets we want. ^[7] Hence

$$f(N,M)=\sum_{r=0}^{M-1}\binom{c_r}{2}.$$

Nothing about this particular sequence has been used yet — the identity holds for any integer sequence — and we no longer need the blocks, only how often each remainder appears.

2. The checkpoint $f(10,10)$ by hand

With the ten given terms, the running totals $S_0,\dots,S_{10}$ are

$$0,1,2,2,5,5,8,13,17,18,27,$$

and their remainders modulo $10$ are

$$0,1,2,2,5,5,8,3,7,8,7.$$

Four remainders occur twice: $2$ at the indices $n=2,3$; $5$ at $n=4,5$; $8$ at $n=6,9$; and $7$ at $n=8,10$. Three more — $0$, $1$ and $3$ — occur once, and $4$, $6$ and $9$ not at all, so only the four groups of size $2$ contribute, $\binom{2}{2}=1$ each, and $f(10,10)=4$. Translating each index pair $\{i,j\}$ back with $p=i+1$, $q=j$ gives $(3,3)$, $(5,5)$, $(7,9)$, $(9,10)$ — the four pairs named in the statement.

3. A weighted running total

The definition of $a_n$ sums all earlier terms, which looks quadratic in $n$; one extra variable removes that. Let

$$W_n=\sum_{k=1}^{n}k\,a_k,\qquad W_0=0.$$

Starting from $a_1=1$, so $W_1=1$, the definition reads $a_n=W_{n-1}\bmod n$ for $n\ge2$, and $W$ is updated by $W_n=W_{n-1}+n\,a_n$. Each term now costs a constant number of operations, so the first $n$ terms cost $O(n)$ in total — still hopeless at $n=10^{12}$. We need structure, not speed.

4. The terms fall into blocks of six

Writing the first thirty terms in rows of six makes the structure visible.

$1,\ 1,\ 0,\ 3,\ 0,\ 3$
$5,\ 4,\ 1,\ 9,\ 1,\ 6$
$9,\ 7,\ 2,\ 15,\ 2,\ 9$
$13,\ 10,\ 3,\ 21,\ 3,\ 12$
$17,\ 13,\ 4,\ 27,\ 4,\ 15$

Read down the columns instead of along the rows. From each row to the next the entries grow by $4$, $3$, $1$, $6$, $1$ and $3$ — the same six steps every time, so every column is an arithmetic progression. That suggests the following claim, proved below for every $m\ge0$:

$$\begin{aligned}a_{6m+1}&=4m+1,\\ a_{6m+2}&=3m+1,\\ a_{6m+3}&=m,\\ a_{6m+4}&=6m+3,\\ a_{6m+5}&=m,\\ a_{6m+6}&=3m+3.\end{aligned}$$

A sequence given by one polynomial on each residue class of its index is called quasi-polynomial; the polynomials are its constituents and the number of classes its period. ^[3] All six constituents here have degree $1$. The sequence is catalogued as OEIS A372865, which records this closed form and a rational generating function. ^[2]

5. Proof of the closed form

We argue by induction on the block index $m$, carrying the weighted total at the end of each block as the invariant:

$$W_{6m}=36m^{3}-3m.$$

Base. At $m=0$ we have $W_0=0$ and $a_1=1$, the given initial value, which matches $4m+1$. The recurrence then generates the first block: $W_1=1$ gives $a_2=1$, $W_2=3$ gives $a_3=0$, $W_3=3$ gives $a_4=3$, $W_4=15$ gives $a_5=0$, $W_5=15$ gives $a_6=3$, and $W_6=33=36-3$.

Step. Assume $W_{6m}=36m^{3}-3m$ with $m\ge1$. We walk through the block in order: reduce the total we already have to obtain the next term, then add that term times its index to the total, and repeat. The totals that arise along the way are

$$\begin{aligned}W_{6m+1}&=36m^{3}+24m^{2}+7m+1,\\ W_{6m+2}&=36m^{3}+42m^{2}+19m+3,\\ W_{6m+3}&=36m^{3}+48m^{2}+22m+3,\\ W_{6m+4}&=36m^{3}+84m^{2}+64m+15,\\ W_{6m+5}&=36m^{3}+90m^{2}+69m+15.\end{aligned}$$

Each one is the previous total plus $(6m+j)\,a_{6m+j}$ with the term just established. Since $a_{6m+j+1}=W_{6m+j}\bmod(6m+j+1)$, every stage asks for a reduction modulo $6m+j+1$, and three rules do all the work. Modulo $6m+j+1$ we have $6m\equiv-(j+1)$; squaring gives $36m^{2}\equiv(j+1)^{2}$, hence $36m^{3}\equiv(j+1)^{2}m$; and a middle term $6\alpha m^{2}$, being $\alpha m$ times $6m$, satisfies $6\alpha m^{2}\equiv-\alpha(j+1)m$. Every coefficient of $m^{2}$ above is a multiple of $6$, so all three rules apply. Carrying them out:

$j=0$: $m-3m=-2m$, so $W_{6m}\equiv-2m$; adding $6m+1$ gives $4m+1$, which lies in $\{0,\dots,6m\}$. So $a_{6m+1}=4m+1$.
$j=1$: $4m-8m+7m=3m$, so $W_{6m+1}\equiv 3m+1$, already below $6m+2$. So $a_{6m+2}=3m+1$.
$j=2$: $9m-21m+19m=7m$, so $W_{6m+2}\equiv 7m+3$; subtracting $6m+3$ leaves $m$. So $a_{6m+3}=m$.
$j=3$: $16m-32m+22m=6m$, so $W_{6m+3}\equiv 6m+3$, already below $6m+4$. So $a_{6m+4}=6m+3$.
$j=4$: $25m-70m+64m=19m$, so $W_{6m+4}\equiv 19m+15$; subtracting $3(6m+5)$ leaves $m$. So $a_{6m+5}=m$.
$j=5$: $36m-90m+69m=15m$, so $W_{6m+5}\equiv 15m+15$; subtracting $2(6m+6)$ leaves $3m+3$, which is below $6m+6$. So $a_{6m+6}=3m+3$.

In each line we produced a number congruent to $W$ that is non-negative and strictly smaller than the modulus. Division with remainder has a unique answer in that range, so each of these numbers is the remainder. ^[5] Finally, adding $(6m+6)(3m+3)$, which is $18m^{2}+36m+18$, to $W_{6m+5}$ gives

$$\begin{aligned}W_{6m+6}&=36m^{3}+108m^{2}\\ &\qquad+105m+33\\ &=36(m+1)^{3}-3(m+1),\end{aligned}$$

so the invariant holds for $m+1$ and the induction is complete. The restriction $m\ge1$ is needed only at $j=0$: at $m=0$ the modulus $6m+1$ is $1$, every remainder is $0$, and $a_1=1$ comes from the initial condition instead.

6. The remainders repeat every $6M$ steps

Fix a modulus $M\ge1$; two consequences of the closed form give the periodicity.

First, term by term: replacing $m$ by $m+M$ in the six formulas adds $4M$, $3M$, $M$, $6M$, $M$ and $3M$, all multiples of $M$. Hence

$$a_{n+6M}\equiv a_n\pmod M\qquad(n\ge1).$$

Second, one whole stretch of $6M$ terms starting at the beginning sums to a multiple of $M$. Adding the six formulas of block $m$, the coefficients of $m$ give $4+3+1+6+1+3=18$ and the constants give $1+1+0+3+0+3=8$, so the block contributes $18m+8$ and

$$S_{6M}=\sum_{m=0}^{M-1}(18m+8)=M(9M-1).$$

Combining the two, for every $n\ge0$,

$$\begin{aligned}S_{n+6M}&=S_{6M}+\sum_{i=1}^{n}a_{i+6M}\\ &\equiv S_n\pmod M.\end{aligned}$$

So the remainders $S_n\bmod M$ are periodic with period $6M$, starting already at $n=0$. The factor $6$ is not obviously slack: a direct search over $1\le M\le80$ finds the smallest period to be $1$ at $M=1$, where every remainder is $0$, and exactly $6M$ for every other modulus in that range. The method below needs only a valid period, which the proof above supplies for every $M\ge1$.

The running totals have a closed form too, used by a second computation below. Writing $n=6m+r$ with $0\le r\lt6$ and summing the first $r$ formulas of block $m$ on top of $S_{6m}=9m^{2}-m$ gives

$$S_{6m+r}=9m^{2}-m+\lambda_r m+\mu_r,$$

with $(\lambda_r,\mu_r)$ equal to $(0,0)$, $(4,1)$, $(7,2)$, $(8,2)$, $(14,5)$ and $(15,5)$ for $r=0,\dots,5$. Replacing $m$ by $m+M$ changes this degree-$2$ expression by $M(18m+9M-1+\lambda_r)$, a multiple of $M$ — the periodicity again, in one line.

7. One period and a tail

Put $P=6M$ and write $N=QP+R$ with $Q=\lfloor N/P\rfloor$ and $0\le R\lt P$. Since $S_n\bmod M$ depends only on $n\bmod P$, it is enough to know the remainders of $S_0,\dots,S_{P-1}$ and how many indices of $\{0,1,\dots,N\}$ fall into each class modulo $P$. For $0\le j\lt P$ that class consists of the numbers $j+kP$ with $k\ge0$ and $j+kP\le QP+R$. Every $k$ from $0$ to $Q-1$ qualifies, because $j\lt P$ gives $j+(Q-1)P\lt QP$; when $Q=0$ that list is simply empty. The one further choice $k=Q$ qualifies exactly when $j\le R$, and any larger $k$ exceeds $N$. So the class of $j$ has $Q+1$ members when $j\le R$ and $Q$ members otherwise. With the two histograms

$$h_r=\#\{\,0\le j\lt P:S_j\bmod M=r\,\},$$

$$u_r=\#\{\,0\le j\le R:S_j\bmod M=r\,\},$$

we therefore get

$$c_r=Q\,h_r+u_r,$$

$$f(N,M)=\sum_{r=0}^{M-1}\binom{c_r}{2}.$$

As a check on the bookkeeping, $\sum_r h_r=P$ and $\sum_r u_r=R+1$, so $\sum_r c_r=QP+R+1=N+1$: every index is counted once. The huge bound $N$ has disappeared into one division, and a pass of length $6M$ is all that is left.

8. A worked example from end to end

Take $N=20$ and $M=2$, so $P=12$. The running totals $S_0,\dots,S_{12}$ are

$$0,1,2,2,5,5,8,13,17,18,27,28,34,$$

whose parities for $n=0,\dots,11$ are

$$0,1,0,0,1,1,0,1,1,0,1,0.$$

Note that $S_{12}=34=M(9M-1)$ at $M=2$ is even, so $S_{12}$ and $S_0$ have the same parity, exactly as Section 6 predicts. Over the full period both parities occur six times, so $h_0=h_1=6$. Here $Q=\lfloor20/12\rfloor=1$ and $R=8$; among $S_0,\dots,S_8$ the even ones sit at $n=0,2,3,6$ and the odd ones at $n=1,4,5,7,8$, so $u_0=4$ and $u_1=5$. Then

$$\begin{aligned}c_0&=1\cdot6+4=10,\\ c_1&=1\cdot6+5=11,\end{aligned}$$

whose total is $21=N+1$, and

$$\begin{aligned}f(20,2)&=\binom{10}{2}+\binom{11}{2}\\ &=45+55=100.\end{aligned}$$

A direct double loop over all $\binom{21}{2}=210$ pairs $(p,q)$ also returns $100$, so the regrouping is correct on this instance.

How the Code Works

All three programs put Sections 1, 6 and 7 into a single function; the C++ version ^[10] is described here. count_pairs(N, M) starts with the three quantities of Section 7: period is $6M$, full_cycles is $Q$ and remainder is $R$. Two vectors of $M$ counters, freq_in_period and freq_in_tail, hold the histograms $h_r$ and $u_r$. Two scalars carry the state of the recurrence: weighted_sum is $W_n$ of Section 3 and prefix_mod is $S_n\bmod M$.

The main loop runs n from $0$ to period - 1 and, crucially, counts before it advances. It first increments freq_in_period[prefix_mod], and also freq_in_tail[prefix_mod] while n <= remainder; only then does it form idx = n + 1, take a = weighted_sum % idx (or $1$ when idx is $1$), add idx * a to weighted_sum and move prefix_mod on by a % M. Because the counting comes first, iteration n records $S_n$, so the loop covers $S_0,\dots,S_{P-1}$, including the empty prefix $S_0=0$ that the index $i=p-1=0$ of Section 1 needs. The tail condition is n <= remainder and not n < remainder for the same reason: $u_r$ counts the $R+1$ indices $0,\dots,R$.

A second loop over the $M$ remainders forms count = full_cycles * freq_in_period[r] + freq_in_tail[r], which is $c_r$, and adds count * (count - 1) / 2. One of the two consecutive factors is even, so the halving is exact. A remainder that never occurs is harmless even though count - 1 then wraps around: unsigned wraparound is well defined, and the product with a zero factor is still zero. to_string_u128 finally extracts the decimal digits by repeated division, reverses them and returns a string, which main writes to std::cout; the standard stream has no overload for unsigned __int128.

Only the C++ version ships self-checks. run_checkpoints regenerates the first ten terms and compares them with $1,1,0,3,0,3,5,4,1,9$, then verifies count_pairs(10, 10) == 4 and count_pairs(10000, 1000) == 97158, the two checkpoints of the statement; main runs them unless --skip-checkpoints is passed and exits with status $2$ on failure. The Python and Java files compute the target directly, with no checkpoint routine.

Integer sizes in the three languages

The weighted total is the one genuinely large intermediate value. By Section 5, $W_{6M}=36M^{3}-3M$, about $3.6\cdot10^{19}$ at $M=10^{6}$ and so nearly twice $2^{64}\approx1.8\cdot10^{19}$. A 64-bit accumulator would wrap during the last blocks of the pass, so C++ keeps weighted_sum in unsigned __int128, a GCC and Clang extension. ^[8]

Java has no such type and emulates one with two long bit patterns. A long is signed, so let $H$ and $L$ be the unsigned interpretations of sumHigh and sumLow; at the start of iteration n they satisfy $W_n=2^{64}H+L$. The next term is $W_n$ reduced modulo idx, which the program assembles from $r_{64}=2^{64}\bmod\mathtt{idx}$, obtained as (Long.remainderUnsigned(-1L, idx) + 1) % idx, and the two reduced half-words, as (hMod * r64 + lMod) % idx. Adding idx * a to the low word may overflow it, which Long.compareUnsigned(sumLow, val) detects as a carry into the high word. Both idx * a and hMod * r64 + lMod stay below $3.6\cdot10^{13}$ here, so each fits in a signed long. ^[9]

For the final accumulation only $\sum_r c_r=N+1$ is known in advance, so the a priori worst case is a single group of size $N+1$, giving $\binom{N+1}{2}\approx5\cdot10^{23}$ — far beyond the signed 64-bit range. C++ stays in 128 bits, Java switches to BigInteger for this loop, and Python's integers are unbounded anyway. ^[9] Python otherwise mirrors the C++ function, with n_limit and m_val as its parameters and two lists of $M$ counters.

One design decision is worth naming. Once Section 4 is proved, each $a_n$ could be read off from $n$ alone and the weighted total, the reason for the wide accumulator, would disappear. The programs instead generate the first period from the original recurrence and use the proved periodicity only to extend its histogram to all $N+1$ indices, so none of the six term formulas has to be implemented. The final pair count needs its own size analysis either way.

Complexity and Verification

The main loop makes $6M$ iterations and the final loop $M$, so the method performs $O(M)$ arithmetic operations and stores $O(M)$ counters. No loop bound depends on $N$: it enters only as $Q$ and $R$, which fix the multiplicities and the length of the tail. With fixed-width arithmetic the running time is therefore $O(M)$; in Python the cost of an operation also grows with the size of its operands. C++ and Java hold two arrays of $M$ 64-bit counters, a payload of $16$ MB at $M=10^{6}$ before allocator overhead, whereas Python holds lists of references to arbitrary-precision integers, for which that calculation says nothing; measured peak resident memory was about $18$ MB for C++ and $28$ MB for Python. At the target parameters the pass is $6\,000\,000$ iterations long, with $Q=166\,666$ and $R=4\,000\,000$, and on one desktop machine the three programs took about $0.07$ s (C++), $0.15$ s (Java) and $0.95$ s (Python).

Beyond the two packaged checkpoints, the following was computed for this article.

The closed form of Section 4 agrees with the definition for every $n\le60\,000$, and the invariants $W_{6m}=36m^{3}-3m$ and $S_{6m}=9m^{2}-m$ hold for every $m\le10\,000$.
The five intermediate polynomials for $W_{6m+j}$ were checked for every $m\le5\,000$, and the six reductions of Section 5 gave the claimed remainders for every $m$ from $1$ to $5\,000$. At $m=0$ the five reductions with $j\ge1$ hold as well, while the one with $j=0$ does not — exactly the exception noted in Section 5.
The quasi-polynomial for $S_{6m+r}$ with the coefficients $(\lambda_r,\mu_r)$ of Section 6 matches the running totals for every $n\le60\,000$, and $a_{n+6M}\equiv a_n\pmod M$ was confirmed for all $M\le40$ and $n\le2\,000$.
The published Python function agrees with a direct double loop over all pairs for every $N\le120$ and $M\le12$, which is $1440$ instances, and reproduces $f(10,10)=4$ and $f(10^{4},10^{3})=97158$.
At the target parameters the C++, Python and Java programs return the same value, and so does a fourth program written only for this check. It uses the closed form instead of the recurrence: for each $r$ it runs $m$ over $0,\dots,M-1$, evaluates $9m^{2}-m+\lambda_r m+\mu_r$ modulo $M$ and adds that class's multiplicity $\lfloor(L_r-m)/M\rfloor+1$, where $L_r=\lfloor(N-r)/6\rfloor$. It reaches the prefix remainders by a different route and shares no code, but it does share the closed form, the periodicity and the pair-counting identity, so its agreement cross-checks two implementations rather than proving those three facts. It reproduces the small instances and both checkpoints too.

No finite run can test the periodicity the programs rely on at every index, and the checkpoints do not establish it: that is the work of the proofs. Section 5 settles the closed form, Section 6 derives the period, and Sections 1 and 7 justify the pair count; the computations above only confirm that the programs follow those arguments.

Footnotes and References

The references below supply the background facts used above. The reduction to remainder histograms, the induction proof of the closed form, the period $6M$ and the counting formula are derived explicitly in this article.

Project Euler 326 — Modulo Summations. The official statement defines the sequence by a(1) = 1 and a(n) = (sum of k·a(k) for k = 1 to n−1) mod n, lists its first ten terms as 1, 1, 0, 3, 0, 3, 5, 4, 1, 9, defines f(N, M) as the number of pairs (p, q) with 1 ≤ p ≤ q ≤ N whose block sum is divisible by M, and asks for f(10¹², 10⁶). Its two given values, f(10, 10) = 4 with the pairs (3,3), (5,5), (7,9), (9,10) and f(10⁴, 10³) = 97158, serve here as checkpoints.
OEIS A372865. The entry is named "a(n) = (Sum_{k = 1..n-1} k*a(k)) (mod n) with a(1) = 1" and lists the first thirty terms, which match the five rows of six in Section 4. Its formula section calls the sequence quasipolynomial in n, proved by induction, and gives a(6n) = 3n for n ≥ 1 together with the five formulas for a(6n+1) … a(6n+5) proved in Section 5, as well as a rational generating function.
Wikipedia — Quasi-polynomial. A function q on the integers is a quasi-polynomial when there are a positive integer s and polynomials p₀, …, p₍ₛ₋₁₎ with q(n) = pᵢ(n) whenever n ≡ i (mod s); these are its constituents, and the smallest such s is its minimal period. Sections 4 and 6 use this language.
Wikipedia — Modular arithmetic. Two integers a and b are congruent modulo m when their difference a − b is an integer multiple of m, and congruence is compatible with addition, subtraction and multiplication. The first fact turns the divisibility condition of Section 1 into an equality of remainders, the second licenses the polynomial reductions of Section 5.
Wikipedia — Euclidean division. For integers a and b with b nonzero there exist unique integers q and r with a = bq + r and 0 ≤ r < |b|. The uniqueness is what allows Section 5 to identify a remainder by exhibiting any congruent value that lies in the right range.
Wikipedia — Prefix sum. The prefix sum, cumulative sum or scan of a sequence is defined there as the sequence of its running totals. It does not state the block identity used here, which Section 1 derives from that definition.
Wikipedia — Binomial coefficient. The article defines C(n, k) as the number of k-element subsets of an n-element set and gives the factorial formula for 0 ≤ k ≤ n. Section 1 uses the case k = 2 to count the index pairs inside one remainder group.
GCC manual — 128-bit integers. The page documents the __int128 and unsigned __int128 extensions, supported on targets with an integer mode wide enough to hold 128 bits. The C++ solution holds the weighted total and the final sum in unsigned __int128.
Java — class Long documents remainderUnsigned, which divides with both arguments and the result read as unsigned 64-bit values, and compareUnsigned, which compares two long values as unsigned; the Java solution builds its wide accumulator on them. Java — BigInteger provides immutable arbitrary-precision integers, used there for the final sum, and Python — numeric types states that Python integers have unlimited precision.
C++ — Euler326.cpp. The linked immutable revision contains count_pairs, run_checkpoints, to_string_u128 and the --skip-checkpoints option described above, and is the implementation record for this article. The fourth program of the verification section is not part of it.

Problem 326 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

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

std::string to_string_u128(u128 value) {
    if (value == 0) {
        return "0";
    }
    std::string s;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        s.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    std::reverse(s.begin(), s.end());
    return s;
}

u128 count_pairs(const u64 N, const int M) {
    const u64 period = 6ULL * static_cast<u64>(M);
    const u64 full_cycles = N / period;
    const u64 remainder = N % period;

    std::vector<u64> freq_in_period(static_cast<std::size_t>(M), 0ULL);
    std::vector<u64> freq_in_tail(static_cast<std::size_t>(M), 0ULL);

    u128 weighted_sum = 0;  // sum_{k=1}^{n} k * a_k
    int prefix_mod = 0;     // sum_{i=1}^{n} a_i (mod M)

    for (u64 n = 0; n < period; ++n) {
        ++freq_in_period[static_cast<std::size_t>(prefix_mod)];
        if (n <= remainder) {
            ++freq_in_tail[static_cast<std::size_t>(prefix_mod)];
        }

        const u64 idx = n + 1ULL;
        const u64 a = (idx == 1ULL) ? 1ULL : static_cast<u64>(weighted_sum % idx);
        weighted_sum += static_cast<u128>(idx) * static_cast<u128>(a);

        prefix_mod += static_cast<int>(a % static_cast<u64>(M));
        if (prefix_mod >= M) {
            prefix_mod %= M;
        }
    }

    u128 answer = 0;
    for (int r = 0; r < M; ++r) {
        const u128 count = static_cast<u128>(full_cycles) * static_cast<u128>(freq_in_period[static_cast<std::size_t>(r)]) +
                           static_cast<u128>(freq_in_tail[static_cast<std::size_t>(r)]);
        answer += count * (count - 1) / 2;
    }

    return answer;
}

bool run_checkpoints() {
    {
        // Statement: first 10 values of a_n.
        std::array<u64, 10> expected = {1, 1, 0, 3, 0, 3, 5, 4, 1, 9};
        std::array<u64, 10> got{};

        u128 weighted_sum = 0;
        for (u64 n = 1; n <= 10; ++n) {
            const u64 a = (n == 1ULL) ? 1ULL : static_cast<u64>(weighted_sum % n);
            got[static_cast<std::size_t>(n - 1)] = a;
            weighted_sum += static_cast<u128>(n) * static_cast<u128>(a);
        }

        if (got != expected) {
            std::cerr << "Checkpoint failed: first 10 sequence elements mismatch\n";
            return false;
        }
    }

    if (count_pairs(10ULL, 10) != static_cast<u128>(4ULL)) {
        std::cerr << "Checkpoint failed: f(10,10)\n";
        return false;
    }

    if (count_pairs(10000ULL, 1000) != static_cast<u128>(97158ULL)) {
        std::cerr << "Checkpoint failed: f(10^4,10^3)\n";
        return false;
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    bool skip_checkpoints = false;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            skip_checkpoints = true;
        } else {
            std::cerr << "Unknown argument: " << arg << '\n';
            return 1;
        }
    }

    if (!skip_checkpoints && !run_checkpoints()) {
        return 2;
    }

    const u128 answer = count_pairs(1000000000000ULL, 1000000);
    std::cout << to_string_u128(answer) << '\n';
    return 0;
}

Python

def count_pairs(n_limit=1000000000000, m_val=1000000):
    period = 6 * m_val
    full_cycles = n_limit // period
    remainder = n_limit % period
    
    freq_in_period = [0] * m_val
    freq_in_tail = [0] * m_val
    
    weighted_sum = 0
    prefix_mod = 0
    
    for n in range(period):
        freq_in_period[prefix_mod] += 1
        if n <= remainder:
            freq_in_tail[prefix_mod] += 1
            
        idx = n + 1
        a = 1 if idx == 1 else weighted_sum % idx
            
        weighted_sum += idx * a
        prefix_mod = (prefix_mod + a % m_val) % m_val
        
    answer = 0
    for r in range(m_val):
        count = full_cycles * freq_in_period[r] + freq_in_tail[r]
        answer += count * (count - 1) // 2
        
    return str(answer)

def solve():
    return count_pairs()

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

Java

import java.util.*;

public class Euler326 {
    public static String solve() {
        long limitN = 1000000000000L;
        int mVal = 1000000;

        long period = 6L * mVal;
        long fullCycles = limitN / period;
        long remainder = limitN % period;

        long[] freqInPeriod = new long[mVal];
        long[] freqInTail = new long[mVal];

        long sumLow = 0;
        long sumHigh = 0;
        int prefixMod = 0;

        for (long n = 0; n < period; n++) {
            freqInPeriod[prefixMod]++;
            if (n <= remainder) {
                freqInTail[prefixMod]++;
            }

            long idx = n + 1;
            long a;
            if (idx == 1) {
                a = 1;
            } else {
                long r64 = (Long.remainderUnsigned(-1L, idx) + 1) % idx;
                long hMod = Long.remainderUnsigned(sumHigh, idx);
                long lMod = Long.remainderUnsigned(sumLow, idx);
                a = (hMod * r64 + lMod) % idx;
            }

            long val = idx * a;
            sumLow += val;
            if (Long.compareUnsigned(sumLow, val) < 0) {
                sumHigh++;
            }

            prefixMod = (int) ((prefixMod + a % mVal) % mVal);
        }

        // Java 128 bit answer... wait!
        // answer can be (6,000,000)^2 / 2 ~ 1.8 * 10^13 ? No!
        // max value of answer:
        // sum_r (count_r * count_r / 2)
        // count_r is full_cycles * freq_in_period + freq_in_tail
        // fullCycles = 10^12 / 6M = 166,666.
        // sum(freq) = 6,000,000.
        // So sum(count) = 10^12.
        // answer is at worst 10^12 * 10^12 / 2 = 5 * 10^{23}.
        // This EXCEEDS Long.MAX_VALUE ~ 9 * 10^{18}.
        // We MUST use BigInteger for the final answer accumulation!

        java.math.BigInteger answer = java.math.BigInteger.ZERO;

        java.math.BigInteger fcBig = java.math.BigInteger.valueOf(fullCycles);

        for (int r = 0; r < mVal; r++) {
            java.math.BigInteger count = fcBig.multiply(java.math.BigInteger.valueOf(freqInPeriod[r]))
                    .add(java.math.BigInteger.valueOf(freqInTail[r]));
            java.math.BigInteger terms = count.multiply(count.subtract(java.math.BigInteger.ONE))
                    .divide(java.math.BigInteger.valueOf(2));
            answer = answer.add(terms);
        }

        return answer.toString();
    }

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