Detailed mathematical approach

Problem Summary

For a prime \(p\) and an integer \(q\), let \(A_q(p)\) denote the number of \(p\)-element subsets of \(\{1,2,\dots,qp\}\) whose element-sum is divisible by \(p\). Problem 635 asks for

$\sum_{p\le 10^8}\bigl(A_2(p)+A_3(p)\bigr)\pmod{10^9+9}.$

A direct subset search is hopeless: even for a single large prime, the number of candidate subsets is astronomically large. The solution works because, for prime moduli, the counting problem collapses to compact closed forms involving binomial coefficients.

Mathematical Approach

The main task is to count \(p\)-element subsets whose sum is congruent to \(0 \pmod p\). A roots-of-unity filter extracts exactly those sums, and the residue pattern of \(\{1,\dots,qp\}\) makes the nontrivial filter terms identical.

Step 1: Encode \(p\)-element subsets with a generating function

Fix \(q\in\{2,3\}\) and a prime \(p\). Consider the bivariate generating function

$F_q(x,y)=\prod_{n=1}^{qp}(1+y x^n).$

Choosing the term \(y x^n\) means that \(n\) is included in the subset; choosing \(1\) means it is omitted. Therefore the coefficient of \(y^k x^m\) counts \(k\)-element subsets with total sum \(m\).

In particular, \(A_q(p)\) is obtained by taking \(k=p\) and summing only those coefficients whose exponent of \(x\) is a multiple of \(p\).

Step 2: Use a roots-of-unity filter to keep only sums divisible by \(p\)

Let \(\omega\) be a primitive \(p\)-th root of unity. The standard identity

$\frac{1}{p}\sum_{j=0}^{p-1}\omega^{jm}= \begin{cases} 1,& p\mid m,\\ 0,& p\nmid m \end{cases}$

filters out exactly the exponents \(m\) that are divisible by \(p\). Hence

$A_q(p)=\frac{1}{p}\sum_{j=0}^{p-1}[y^p]\prod_{n=1}^{qp}(1+y\omega^{jn}).$

Now group the numbers \(1,2,\dots,qp\) by their residues modulo \(p\). Each residue class appears exactly \(q\) times, so when \(j\neq 0\), multiplication by \(j\) merely permutes the residue classes and we get

$\prod_{n=1}^{qp}(1+y\omega^{jn})=\left(\prod_{r=0}^{p-1}(1+y\omega^r)\right)^q.$

Step 3: Collapse the nonzero filter terms

The cyclotomic product identity

$\prod_{r=0}^{p-1}(1-t\omega^r)=1-t^p$

gives, after substituting \(t=-y\),

$\prod_{r=0}^{p-1}(1+y\omega^r)=1-(-y)^p.$

Therefore every nonzero filter term contributes

$[y^p](1-(-y)^p)^q.$

If \(p\) is odd, then \(1-(-y)^p=1+y^p\), so the coefficient of \(y^p\) is simply \(q\). The \(j=0\) term is different: it is

$[y^p](1+y)^{qp}=\binom{qp}{p}.$

Step 4: Derive the closed forms for odd primes

For every odd prime \(p\), the \(p-1\) nonzero filter terms all contribute \(q\). Combining them with the \(j=0\) term yields

$A_q(p)=\frac{\binom{qp}{p}+q(p-1)}{p}.$

Specializing to the two quantities in the problem gives

$A_2(p)=\frac{\binom{2p}{p}+2(p-1)}{p},$

$A_3(p)=\frac{\binom{3p}{p}+3(p-1)}{p}$

for every odd prime \(p\).

Step 5: Handle the special case \(p=2\)

When \(p=2\), the same identity becomes

$\prod_{r=0}^{1}(1+y\omega^r)=1-y^2,$

so the nonzero filter term contributes \([y^2](1-y^2)^q=-q\) instead of \(+q\). Hence

$A_q(2)=\frac{\binom{2q}{2}-q}{2}.$

In particular,

$A_2(2)=2,\qquad A_3(2)=6,$

which explains why the implementation treats \(p=2\) separately.

Worked Example: \(p=5\)

For \(p=5\), both formulas apply directly:

$A_2(5)=\frac{\binom{10}{5}+2\cdot 4}{5}=\frac{252+8}{5}=52,$

$A_3(5)=\frac{\binom{15}{5}+3\cdot 4}{5}=\frac{3003+12}{5}=603.$

So the prime \(5\) contributes \(52+603=655\) to the final sum. As another useful checkpoint,

$\sum_{p\le 10}A_2(p)=A_2(2)+A_2(3)+A_2(5)+A_2(7)=2+8+52+492=554,$

which matches the internal verification used by the implementation.

How the Code Works

The C++, Python, and Java implementations first sieve all primes up to \(10^8\). They never store complete factorial and inverse-factorial tables up to \(3\cdot 10^8\); that would be unnecessarily large. Instead, they keep only a few checkpoint values for each prime.

A forward scan computes factorials modulo \(10^9+9\) and records the values needed later at \(p-1\), \(2p\), and \(3p\). Then the implementation computes the inverse of \((3\cdot 10^8)!\) once using fast modular exponentiation and performs a backward scan to reconstruct inverse factorials, recording only the checkpoints at \(p\) and \(2p\).

Those saved values are enough to recover

$\frac{1}{p}=\frac{(p-1)!}{p!},\qquad \binom{2p}{p}=\frac{(2p)!}{p!\,p!},\qquad \binom{3p}{p}=\frac{(3p)!}{p!\,(2p)!}$

in constant time for each prime. The modulus is prime and larger than \(3\cdot 10^8\), so every factorial value in that range is invertible modulo \(10^9+9\). After reconstructing the binomial terms, the program applies the \(p=2\) special case or the odd-prime formula, adds the two contributions, and reduces modulo \(10^9+9\) throughout.

Complexity Analysis

Let \(L=10^8\). The prime sieve costs \(O(L\log\log L)\) time. The forward factorial scan and the backward inverse-factorial scan both run to \(3L\), so they add \(O(L)\) time. Once the checkpoint tables are ready, each prime is processed in \(O(1)\) time. The overall running time is therefore \(O(L\log\log L)\).

Memory usage consists of an odd-only sieve plus a constant number of checkpoint arrays indexed by primes. Asymptotically this is \(O(L+\pi(L))\), and the important practical point is that the implementation avoids storing all factorials and inverse factorials up to \(3L\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=635
2. Root-of-unity filter: Wikipedia - Root of unity filter
3. Binomial coefficient: Wikipedia - Binomial coefficient
4. Fermat's little theorem: Wikipedia - Fermat's little theorem
5. Sieve of Eratosthenes: Wikipedia - Sieve of Eratosthenes