Project Euler

Sequences with Nice Divisibility Properties

Let S(n, k, m) denote the number of n -tuples (a_1, a_2,..., a_n) of positive integers satisfying: 1. 1 <= a_i <= m for all 1 <= i <= n. 2. a_1 + a_2 +... + a_n equiv 0 (mod k). Compute S(n, k, m)...

Source sync Apr 19, 2026

Problem #0511

Level Level 23

Solved By 504

Languages C++, Python

Answer 935247012

Length 381 words

modular_arithmeticalgebranumber_theory

Let $Seq(n,k)$ be the number of positive-integer sequences $\{a_i\}_{1 \le i \le n}$ of length $n$ such that:

$n$ is divisible by $a_i$ for $1 \le i \le n$, and
$n + a_1 + a_2 + \cdots + a_n$ is divisible by $k$.

Example:

$Seq(3,4) = 4$, and the $4$ sequences are: \begin{align*} \{1, 1, 3\} \\ \{1, 3, 1\} \\ \{3, 1, 1\} \\ \{3, 3, 3\} \end{align*}

$Seq(4,11) = 8$, and the $8$ sequences are: \begin{align*} \{1, 1, 1, 4\} \\ \{1, 1, 4, 1\} \\ \{1, 4, 1, 1\} \\ \{4, 1, 1, 1\} \\ \{2, 2, 2, 1\} \\ \{2, 2, 1, 2\} \\ \{2, 1, 2, 2\} \\ \{1, 2, 2, 2\} \end{align*}

The last nine digits of $Seq(1111,24)$ are $840643584$.

Find the last nine digits of $Seq(1234567898765,4321)$.

Problem 511: Sequences with Nice Divisibility Properties

Mathematical Foundation

Definition 1. Let $\omega = e^{2\pi i / k}$ denote a fixed primitive $k$ -th root of unity. For $0 \leq j \leq k-1$ , define the character sum

g_j \;=\; \sum_{a=1}^{m} \omega^{ja}.

Theorem 1 (Roots of Unity Filter). For any integer $s$ ,

[k \mid s] \;=\; \frac{1}{k}\sum_{j=0}^{k-1} \omega^{js},

where $[P]$ denotes the Iverson bracket (1 if $P$ is true, 0 otherwise).

Proof. We distinguish two cases.

Case 1: $k \mid s$ . Then $\omega^s = 1$ , so $\omega^{js} = 1$ for every $j$ . The sum equals $k$ , and $k/k = 1$ .

Case 2: $k \nmid s$ . Then $\omega^s \neq 1$ . The sum is a geometric series with ratio $\omega^s$ :

\sum_{j=0}^{k-1} \omega^{js} = \frac{(\omega^s)^k - 1}{\omega^s - 1} = \frac{(\omega^k)^s - 1}{\omega^s - 1} = \frac{1 - 1}{\omega^s - 1} = 0.

In both cases the identity holds. $\square$

Theorem 2 (Counting Formula). With the notation above,

S(n, k, m) = \frac{1}{k}\sum_{j=0}^{k-1} g_j^{\,n}.

Proof. By definition,

S(n, k, m) = \sum_{\substack{(a_1,\ldots,a_n) \\ 1 \leq a_i \leq m}} [k \mid a_1 + \cdots + a_n].

Substituting the identity from Theorem 1:

S(n, k, m) = \sum_{\substack{(a_1,\ldots,a_n) \\ 1 \leq a_i \leq m}} \frac{1}{k}\sum_{j=0}^{k-1} \omega^{j(a_1+\cdots+a_n)}.

Since all sums are finite, we may exchange the order of summation (Fubini’s theorem for finite sums):

S(n, k, m) = \frac{1}{k}\sum_{j=0}^{k-1} \sum_{\substack{(a_1,\ldots,a_n) \\ 1 \leq a_i \leq m}} \prod_{i=1}^{n} \omega^{j a_i} = \frac{1}{k}\sum_{j=0}^{k-1} \prod_{i=1}^{n}\left(\sum_{a=1}^{m}\omega^{ja}\right) = \frac{1}{k}\sum_{j=0}^{k-1} g_j^{\,n}.

The factorization into a product holds because the exponent $j(a_1 + \cdots + a_n)$ separates multiplicatively, and the constraint on each $a_i$ is identical and independent. $\square$

Lemma 1 (Closed Form for $g_j$ ). For $0 \leq j \leq k - 1$ :

g_j = \begin{cases} m & \text{if } j = 0, \\ \displaystyle \omega^j \cdot \frac{1 - \omega^{jm}}{1 - \omega^j} & \text{if } 1 \leq j \leq k - 1. \end{cases}

Proof. When $j = 0$ , every term $\omega^{0 \cdot a} = 1$ , so $g_0 = \sum_{a=1}^{m} 1 = m$ .

When $1 \leq j \leq k - 1$ , we have $\omega^j \neq 1$ (since $\omega$ is a primitive $k$ -th root). The sum $g_j = \omega^j + \omega^{2j} + \cdots + \omega^{mj}$ is a finite geometric series with first term $\omega^j$ , common ratio $\omega^j$ , and $m$ terms. By the standard geometric series formula:

g_j = \omega^j \cdot \frac{(\omega^j)^m - 1}{\omega^j - 1} = \omega^j \cdot \frac{1 - \omega^{jm}}{1 - \omega^j}. \qquad \square

Corollary 1 (Special Case $m = k$ ). When $m = k$ , we have $S(n, k, k) = k^{n-1}$ .

Proof. For $j = 0$ : $g_0 = k$ . For $1 \leq j \leq k - 1$ : $\sum_{a=1}^{k}\omega^{ja} = \omega^j \sum_{a=0}^{k-1}\omega^{ja} - \omega^{j \cdot 0} + \omega^{jk} = \omega^j \cdot 0 - 1 + 1$ . More directly, $\sum_{a=0}^{k-1}\omega^{ja} = 0$ for $j \not\equiv 0 \pmod{k}$ (by Theorem 1 with $s = j$ ), so $\sum_{a=1}^{k}\omega^{ja} = \sum_{a=0}^{k-1}\omega^{ja} + \omega^{jk} - \omega^{j \cdot 0} = 0 + 1 - 1 = 0$ . Therefore $S(n,k,k) = k^n / k = k^{n-1}$ .

A combinatorial verification: choose $a_1, \ldots, a_{n-1} \in \{1,\ldots,k\}$ freely ( $k^{n-1}$ ways), then $a_n$ is uniquely determined as the element of $\{1,\ldots,k\}$ satisfying $a_n \equiv -(a_1 + \cdots + a_{n-1}) \pmod{k}$ . $\square$

Proposition 1 (Modular Arithmetic Realization). Suppose $p$ is a prime with $k \mid (p-1)$ . Let $g$ be a primitive root modulo $p$ and set $\hat{\omega} = g^{(p-1)/k} \bmod p$ . Then $\hat{\omega}$ is a primitive $k$ -th root of unity in $\mathbb{F}_p$ , and the formula from Theorem 2 can be evaluated entirely in $\mathbb{F}_p$ :

S(n, k, m) \equiv k^{-1} \sum_{j=0}^{k-1} \hat{g}_j^{\,n} \pmod{p},

where $\hat{g}_j$ is the modular analogue of $g_j$ , computed using $\hat{\omega}$ in place of $\omega$ .

Proof. Since $g$ has order $p - 1$ in $\mathbb{F}_p^*$ , the element $\hat{\omega} = g^{(p-1)/k}$ has order exactly $k$ , hence is a primitive $k$ -th root of unity in $\mathbb{F}_p$ . The algebraic identities in Theorems 1, 2, and Lemma 1 hold over any field containing such a root, in particular over $\mathbb{F}_p$ . Division by $k$ is well-defined since $p \nmid k$ (as $k \mid (p-1)$ and $p$ is prime). $\square$

Editorial

Count n-tuples (a_1,…,a_n) with 1 <= a_i <= m and sum divisible by k. Uses the roots-of-unity filter (discrete Fourier transform technique). We precondition: p prime, k | (p - 1). Finally, else.

Pseudocode

Precondition: p prime, k | (p - 1)
else

Complexity Analysis

Time: $O(k \log n)$ . The loop iterates $k$ times; within each iteration, computing $g_j^n \bmod p$ costs $O(\log n)$ via binary exponentiation.
Space: $O(1)$ auxiliary beyond the input.

Answer

\boxed{935247012}

C++ project_euler/problem_511/solution.cpp

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

ll mod_pow(ll base, ll exp, ll mod) {
    ll result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = (__int128)result * base % mod;
        base = (__int128)base * base % mod;
        exp >>= 1;
    }
    return result;
}

ll mod_inv(ll a, ll mod) {
    return mod_pow(a, mod - 2, mod);
}

// Find primitive root of prime p
ll primitive_root(ll p) {
    ll phi = p - 1;
    // Factor phi
    vector<ll> factors;
    ll n = phi;
    for (ll d = 2; d * d <= n; d++) {
        if (n % d == 0) {
            factors.push_back(d);
            while (n % d == 0) n /= d;
        }
    }
    if (n > 1) factors.push_back(n);

    for (ll g = 2; g < p; g++) {
        bool ok = true;
        for (ll f : factors) {
            if (mod_pow(g, phi / f, p) == 1) {
                ok = false;
                break;
            }
        }
        if (ok) return g;
    }
    return -1;
}

// S(n, k, m) mod p using roots of unity filter
// Requires k | (p-1)
ll S_mod(ll n, ll k, ll m, ll p) {
    assert((p - 1) % k == 0);

    ll g = primitive_root(p);
    ll omega = mod_pow(g, (p - 1) / k, p);
    ll inv_k = mod_inv(k, p);

    ll total = 0;
    for (ll j = 0; j < k; j++) {
        ll gj;
        if (j == 0) {
            gj = m % p;
        } else {
            ll w = mod_pow(omega, j, p);
            ll wm = mod_pow(w, m, p);
            // gj = w * (1 - w^m) / (1 - w)
            ll num = w % p * ((1 - wm % p + p) % p) % p;
            ll den = (1 - w % p + p) % p;
            gj = num % p * mod_inv(den, p) % p;
        }
        total = (total + mod_pow(gj, n, p)) % p;
    }

    return total % p * inv_k % p;
}

int main() {
    // Simple case: S(n, k) = k^(n-1) when m = k
    cout << "Simple case verification:" << endl;
    for (int n = 1; n <= 5; n++) {
        for (int k : {2, 3, 5}) {
            ll result = 1;
            for (int i = 0; i < n - 1; i++) result *= k;
            cout << "  S(" << n << ", " << k << ") = " << result << endl;
        }
    }

    // General case with modular arithmetic
    ll p = 1000000007;  // 10^9 + 7
    ll k = 7;
    // Check k | (p-1): (10^9 + 6) % 7 = ?
    if ((p - 1) % k == 0) {
        ll n = 100, m = 20;
        ll result = S_mod(n, k, m, p);
        cout << "\nS(" << n << ", " << k << ", " << m << ") mod " << p
             << " = " << result << endl;
    }

    return 0;
}

Python project_euler/problem_511/solution.py

"""
Problem 511: Sequences with Nice Divisibility Properties
Count n-tuples (a_1,...,a_n) with 1 <= a_i <= m and sum divisible by k.
Uses the roots-of-unity filter (discrete Fourier transform technique).
"""

import numpy as np
from math import gcd

def S_simple(n: int, k: int):
    """
    Count n-tuples from {1,...,k} with sum divisible by k.
    Answer: k^(n-1) (since choosing n-1 values freely determines the last).
    """
    return pow(k, n - 1)

def S_general_float(n: int, k: int, m: int):
    """
    Count n-tuples from {1,...,m} with sum divisible by k.
    Uses roots-of-unity filter with floating-point arithmetic.
    For exact results, use modular arithmetic version.
    """
    omega = np.exp(2j * np.pi / k)
    total = 0.0

    for j in range(k):
        if j == 0:
            gj = m
        else:
            w = omega ** j
            # g_j = sum_{a=1}^{m} w^a = w * (1 - w^m) / (1 - w)
            gj = w * (1 - w ** m) / (1 - w)
        total += gj ** n

    result = total / k
    return int(round(result.real))

def S_general_mod(n: int, k: int, m: int, p: int):
    """
    Count n-tuples from {1,...,m} with sum divisible by k, modulo prime p.
    Requires k | (p - 1) so that a k-th root of unity exists mod p.
    """
    # Find a primitive k-th root of unity mod p
    # First find a primitive root mod p, then take its ((p-1)/k)-th power
    assert (p - 1) % k == 0, f"k={k} does not divide p-1={p-1}"

    # Find primitive root of p
    def primitive_root(p):
        if p == 2:
            return 1
        phi = p - 1
        # Factor phi
        factors = set()
        n = phi
        for d in range(2, int(n**0.5) + 1):
            while n % d == 0:
                factors.add(d)
                n //= d
            if n == 1:
                break
        if n > 1:
            factors.add(n)

        for g in range(2, p):
            ok = True
            for f in factors:
                if pow(g, phi // f, p) == 1:
                    ok = False
                    break
            if ok:
                return g
        return -1

    g = primitive_root(p)
    omega = pow(g, (p - 1) // k, p)

    total = 0
    inv_k = pow(k, p - 2, p)  # modular inverse of k

    for j in range(k):
        if j == 0:
            gj = m % p
        else:
            w = pow(omega, j, p)
            wm = pow(w, m, p)
            # gj = w * (1 - w^m) / (1 - w) mod p
            num = w * (1 - wm) % p
            den = (1 - w) % p
            gj = num * pow(den, p - 2, p) % p

        total = (total + pow(gj, n, p)) % p

    return total * inv_k % p

def brute_force(n: int, k: int, m: int):
    """Brute-force for small n, m: enumerate all tuples."""
    if n == 0:
        return 1 if 0 % k == 0 else 0

    count = 0

    def recurse(depth, running_sum):
        nonlocal count
        if depth == n:
            if running_sum % k == 0:
                count += 1
            return
        for a in range(1, m + 1):
            recurse(depth + 1, running_sum + a)

    recurse(0, 0)
    return count

# Verify simple case
print("Simple case S(n, k) = k^(n-1):")
for n in range(1, 6):
    for k in [2, 3, 5]:
        exact = brute_force(n, k, k)
        formula = S_simple(n, k)
        status = "OK" if exact == formula else "MISMATCH"
        print(f"  S({n}, {k}) = {formula} (brute={exact}) [{status}]")

# Verify general case
print("\nGeneral case S(n, k, m):")
for n, k, m in [(2, 3, 5), (3, 4, 6), (2, 5, 7), (3, 3, 4)]:
    exact = brute_force(n, k, m)
    approx = S_general_float(n, k, m)
    status = "OK" if exact == approx else f"MISMATCH(got {approx})"
    print(f"  S({n}, {k}, {m}) = {exact} (float={approx}) [{status}]")

# Modular computation for larger values
# Need p such that k | (p-1)
# For k=1234567, we need p with (p-1) % 1234567 == 0
# Example with moderate k
k_test = 7
n_test = 100
m_test = 20
p_mod = 10**9 + 7  # (10^9+6) % 7 = (10^9+6) mod 7

# Check if 7 | (p-1)
if (p_mod - 1) % k_test == 0:
    result_mod = S_general_mod(n_test, k_test, m_test, p_mod)
    result_float = S_general_float(n_test, k_test, m_test)
    print(f"\nS({n_test}, {k_test}, {m_test}) mod {p_mod} = {result_mod}")
else:
    print(f"\nCannot use p={p_mod} for k={k_test} (k does not divide p-1)")