Solution

Answer: 1966666166408794329

Let

$$S_n=\sum_{k=1}^n k a_k.$$

The recurrence is

$$a_n=S_{n-1}\pmod n.$$

We want

$$f(N,M)=#{(p,q):1\le p\le q\le N,\ \sum_{i=p}^q a_i\equiv0\pmod M}.$$

The standard prefix-sum trick converts this into a counting problem.

Define

$$B_n=\sum_{i=1}^n a_i,\qquad B_0=0.$$

Then

$$\sum_{i=p}^q a_i = B_q-B_{p-1},$$

so the condition is

$$B_q\equiv B_{p-1}\pmod M.$$

Hence if a residue $r$ occurs $c_r$ times among

$$B_0,B_1,\dots,B_N \pmod M,$$

then it contributes $\binom{c_r}{2}$ valid pairs. Therefore

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

So the entire problem reduces to understanding the sequence $a_n$.

1. Deriving the structure of $a_n$

We compute a few terms:

$$1,1,0,3,0,3,5,4,1,9,1,6,9,7,2,15,\dots$$

A strong $6$-term pattern appears.

We claim:

$$\boxed{ \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}}$$

This can be proved by induction from the defining recurrence.

Computing the prefix sums $B_n$

Summing one full block:

$$(4m+1)+(3m+1)+m+(6m+3)+m+(3m+3)=18m+8.$$

Thus

$$B_{6m} =\sum_{j=0}^{m-1}(18j+8) =9m^2-m.$$

Within a block:

$$\boxed{ \begin{aligned} B_{6m} &= 9m^2-m,\ B_{6m+1} &= 9m^2+3m+1,\ B_{6m+2} &= 9m^2+6m+2,\ B_{6m+3} &= 9m^2+7m+2,\ B_{6m+4} &= 9m^2+13m+5,\ B_{6m+5} &= 9m^2+14m+5. \end{aligned}}$$

These formulas are exact.

2. Periodicity modulo $M$

Reduce modulo $M$.

For $m\mapsto m+M$,

$$9(m+M)^2-(m+M)\equiv 9m^2-m \pmod M,$$

and similarly for all six residue classes.

Therefore:

$$B_{n+6M}\equiv B_n\pmod M.$$

So the sequence $B_n\bmod M$ has period $6M$.

For this problem:

$$M=10^6,\qquad P=6M=6\cdot10^6.$$

3. Counting residues

Let

$$N=qP+r.$$

Then the first $q$ full periods contribute equally, and we only need one extra partial block.

For each residue $x\in[0,M-1]$, let $c_x$ be its total frequency among

$$B_0,B_1,\dots,B_N.$$

Finally:

$$f(N,M)=\sum_x \frac{c_x(c_x-1)}2.$$

4. Python implementation

from collections import defaultdict

def F(N, M):
    P = 6 * M

    # Count residues in one full period
    period_count = [0] * M

    for m in range(M):
        base = (9 * m * m - m) % M

        vals = [
            base,
            (base + 4 * m + 1) % M,
            (base + 7 * m + 2) % M,
            (base + 8 * m + 2) % M,
            (base + 14 * m + 5) % M,
            (base + 15 * m + 5) % M,
        ]

        for v in vals:
            period_count[v] += 1

    q, r = divmod(N, P)

    # Total counts
    total = [q * c for c in period_count]

    # Add the remaining prefix terms
    full_blocks = r // 6
    rem = r % 6

    for m in range(full_blocks + 1):
        base = (9 * m * m - m) % M

        vals = [
            base,
            (base + 4 * m + 1) % M,
            (base + 7 * m + 2) % M,
            (base + 8 * m + 2) % M,
            (base + 14 * m + 5) % M,
            (base + 15 * m + 5) % M,
        ]

        limit = 6 if m < full_blocks else rem + 1

        for i in range(limit):
            total[vals[i]] += 1

    # Count equal pairs
    ans = 0
    for c in total:
        ans += c * (c - 1) // 2

    return ans

print(F(10, 10))         # 4
print(F(10**4, 10**3))   # 97158
print(F(10**12, 10**6))

5. Code walkthrough

Step 1

We use the closed forms for the six values of $B_n$ inside each block of length $6$.

Step 2

We count how many times each residue modulo $M$ appears during one full period of length $6M$.

Step 3

Since the sequence is periodic, we multiply those counts by the number of complete periods.

Step 4

We add the remaining incomplete tail.

Step 5

If a residue appears $c$ times, then choosing any two occurrences gives a valid interval sum divisible by $M$. Hence contribution:

$$\binom c2.$$

6. Verification

The implementation reproduces both values given in the problem:

$$f(10,10)=4,$$

and

$$f(10^4,10^3)=97158.$$

Therefore the derivation is consistent.

The computation for

$$f(10^{12},10^6)$$

gives:

$$1966666166408794329.$$

Answer: 1966666166408794329

Project Euler Problem 326