Problem 987: Cyclotomic Polynomial Evaluation

Mathematical Analysis

Definition and Fundamental Identity

The cyclotomic polynomial is:

Its degree is $\varphi(n)$ (Euler’s totient). The fundamental factorization of $x^n - 1$ over $\mathbb{Z}[x]$ is:

Mobius Inversion Formula

Applying Mobius inversion to the logarithm of (1):

where $\mu$ is the Mobius function. This is the explicit product formula and is both theoretically elegant and computationally useful.

Evaluating at $x = 2$

Substituting $x = 2$ into (2):

Alternatively, from (1) directly:

This gives the divisor recurrence:

Key Properties of $\Phi_n(2)$

1. Mersenne connection: For prime $p$, $\Phi_p(2) = \frac{2^p - 1}{2 - 1} = 2^p - 1 = M_p$, the $p$-th Mersenne number.

2. Aurifeuillean factorizations: For certain $n$, $\Phi_n(2)$ admits algebraic factorizations beyond the cyclotomic structure (e.g., $\Phi_{105}(2)$ has a non-trivial factorization discovered by Aurifeuille).

3. Divisibility: $\Phi_n(2) \mid 2^n - 1$ for all $n \ge 1$, and the multiplicative order of 2 modulo any prime $p \mid \Phi_n(2)$ is exactly $n$.

4. Growth rate: $\Phi_n(2) \approx 2^{\varphi(n)}$ as $n \to \infty$, since the roots of $\Phi_n$ lie on the unit circle and $|2 - e^{i\theta}| \approx 2$ for most $\theta$.

Concrete Examples

Note the surprising result $\Phi_6(2) = 3 = \Phi_2(2)$. In general, small values of $\Phi_n(2)$ occur when $n$ has many small prime factors.

Derivation

Algorithm 1: Divisor Recurrence (used in implementation)

Compute $\Phi_n(2) \bmod p$ for $n = 1, 2, \ldots, N$ using (5):

1. For each $n$, compute $2^n - 1 \bmod p$ via modular exponentiation.
2. Compute $D = \prod_{d \mid n, d < n} \Phi_d(2) \bmod p$ using previously computed values.
3. Set $\Phi_n(2) \equiv (2^n - 1) \cdot D^{-1} \pmod{p}$, where $D^{-1}$ is the modular inverse via Fermat’s little theorem ($p$ is prime).

Algorithm 2: Mobius Product (alternative)

Using (3) directly:

1. Precompute $\mu(k)$ for $k = 1, \ldots, N$ via a sieve.
2. Precompute $2^d - 1 \bmod p$ for $d = 1, \ldots, N$.
3. For each $n$, compute $\Phi_n(2) = \prod_{d \mid n} (2^d - 1)^{\mu(n/d)} \bmod p$.

Factors with $\mu(n/d) = 0$ contribute 1 (skip). Factors with $\mu(n/d) = -1$ contribute a modular inverse.

Verification

Both algorithms yield identical results. Cross-check: $\Phi_1(2) + \Phi_2(2) + \cdots + \Phi_6(2) = 1 + 3 + 7 + 5 + 31 + 3 = 50$.

From (4): $\sum_{d \mid 6} \Phi_d(2) = \Phi_1(2) \cdot \Phi_2(2) \cdot \Phi_3(2) \cdot \Phi_6(2) = 1 \cdot 3 \cdot 7 \cdot 3 = 63 = 2^6 - 1$. Correct.

Proof of Correctness

Theorem. The identity $x^n - 1 = \prod_{d \mid n} \Phi_d(x)$ holds in $\mathbb{Z}[x]$.

Proof. The roots of $x^n - 1$ are all $n$-th roots of unity $\{e^{2\pi i k/n} : 0 \le k < n\}$. Each such root is a primitive $d$-th root of unity for exactly one $d \mid n$ (namely $d = n / \gcd(k, n)$). Since $\Phi_d(x)$ is the minimal polynomial of the primitive $d$-th roots, the factorization follows from unique factorization in $\mathbb{Z}[x]$. $\square$

Corollary. The Mobius inversion formula (2) follows by applying $\mu$-inversion to $\log(x^n - 1) = \sum_{d \mid n} \log \Phi_d(x)$.

Lemma. For prime $p = 10^9 + 7$, the modular computation is exact: $\Phi_n(2) \not\equiv 0 \pmod{p}$ for all $1 \le n \le 500$.

Proof. If $p \mid \Phi_n(2)$, then $\text{ord}_p(2) = n$, i.e., $n \mid p - 1 = 10^9 + 6$. One verifies that $10^9 + 6 = 2 \times 500000003$ and $500000003$ is prime. So the only $n \le 500$ with $n \mid p - 1$ are $n \in \{1, 2\}$. But $\Phi_1(2) = 1$ and $\Phi_2(2) = 3$, neither divisible by $p$. $\square$

Complexity Analysis

• Divisor recurrence: $O(N \sqrt{N})$ for enumerating divisors + $O(N \log N)$ for modular exponentiations = $O(N^{3/2})$ total.
• Mobius product: Same complexity with a different constant factor.
• Space: $O(N)$ to store the computed values.

For $N = 500$, both algorithms run in microseconds.

Answer

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

/*
 * Problem 987: Cyclotomic Polynomial Evaluation
 *
 * Compute S = sum_{n=1}^{500} Phi_n(2) mod (10^9 + 7).
 *
 * Two methods implemented:
 *   1. Divisor recurrence:  Phi_n(2) = (2^n - 1) / prod_{d|n, d<n} Phi_d(2)
 *   2. Mobius product:      Phi_n(2) = prod_{d|n} (2^d - 1)^{mu(n/d)}
 *
 * Both yield the same result; we use method 1 as primary and method 2
 * as a cross-check.
 */

const long long MOD = 1e9 + 7;
const int N = 500;

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

long long modinv(long long a, long long mod = MOD) {
    return power(a, mod - 2, mod);
}

vector<int> get_divisors(int n) {
    vector<int> divs;
    for (int d = 1; (long long)d * d <= n; d++) {
        if (n % d == 0) {
            divs.push_back(d);
            if (d != n / d) divs.push_back(n / d);
        }
    }
    sort(divs.begin(), divs.end());
    return divs;
}

// Method 1: Divisor recurrence
// Phi_n(2) = (2^n - 1) / prod_{d|n, d<n} Phi_d(2)
long long solve_recurrence() {
    vector<long long> phi(N + 1, 0);
    long long total = 0;

    for (int n = 1; n <= N; n++) {
        long long num = (power(2, n, MOD) - 1 + MOD) % MOD;
        long long den = 1;
        for (int d : get_divisors(n)) {
            if (d < n) {
                den = den * phi[d] % MOD;
            }
        }
        phi[n] = num % MOD * modinv(den) % MOD;
        total = (total + phi[n]) % MOD;
    }
    return total;
}

// Method 2: Mobius product formula
// Phi_n(2) = prod_{d|n} (2^d - 1)^{mu(n/d)}
long long solve_mobius() {
    // Sieve Mobius function
    vector<int> mu(N + 1, 0);
    mu[1] = 1;
    vector<bool> is_prime(N + 1, true);
    vector<int> primes;

    for (int i = 2; i <= N; i++) {
        if (is_prime[i]) {
            primes.push_back(i);
            mu[i] = -1;
        }
        for (int p : primes) {
            if ((long long)i * p > N) break;
            is_prime[i * p] = false;
            if (i % p == 0) {
                mu[i * p] = 0;
                break;
            }
            mu[i * p] = -mu[i];
        }
    }

    // Precompute 2^d - 1 mod p
    vector<long long> pow2m1(N + 1);
    for (int d = 1; d <= N; d++) {
        pow2m1[d] = (power(2, d, MOD) - 1 + MOD) % MOD;
    }

    long long total = 0;
    for (int n = 1; n <= N; n++) {
        long long phi_n = 1;
        for (int d : get_divisors(n)) {
            int m = mu[n / d];
            if (m == 1) {
                phi_n = phi_n * pow2m1[d] % MOD;
            } else if (m == -1) {
                phi_n = phi_n * modinv(pow2m1[d]) % MOD;
            }
            // m == 0: skip (contributes 1)
        }
        total = (total + phi_n) % MOD;
    }
    return total;
}

int main() {
    long long ans1 = solve_recurrence();
    long long ans2 = solve_mobius();

    // Cross-check both methods
    assert(ans1 == ans2);

    cout << ans1 << endl;
    return 0;
}

"""
Problem 987: Cyclotomic Polynomial Evaluation

The n-th cyclotomic polynomial Phi_n(x) is the minimal polynomial over Q
whose roots are the primitive n-th roots of unity. We compute:

    S = sum_{n=1}^{500} Phi_n(2)  mod (10^9 + 7)

Key identity:  x^n - 1 = prod_{d|n} Phi_d(x)
Mobius form:   Phi_n(x) = prod_{d|n} (x^d - 1)^{mu(n/d)}
Recurrence:    Phi_n(x) = (x^n - 1) / prod_{d|n, d<n} Phi_d(x)
"""

from math import gcd

MOD = 10**9 + 7

# --- Mobius function sieve ---
def mobius_sieve(n: int):
    """Compute mu(k) for k = 0, 1, ..., n via linear sieve."""
    mu = [0] * (n + 1)
    mu[1] = 1
    is_prime = [True] * (n + 1)
    primes = []
    for i in range(2, n + 1):
        if is_prime[i]:
            primes.append(i)
            mu[i] = -1  # prime => mu = -1
        for p in primes:
            if i * p > n:
                break
            is_prime[i * p] = False
            if i % p == 0:
                mu[i * p] = 0  # p^2 divides i*p
                break
            else:
                mu[i * p] = -mu[i]
    return mu

def divisors(n: int):
    """Return sorted list of divisors of n."""
    divs = []
    d = 1
    while d * d <= n:
        if n % d == 0:
            divs.append(d)
            if d != n // d:
                divs.append(n // d)
        d += 1
    return sorted(divs)

# --- Method 1: Divisor recurrence ---
def solve_recurrence(N: int, mod: int):
    """Compute Phi_n(2) mod p for n=1..N using the divisor recurrence.

    Phi_n(2) = (2^n - 1) / prod_{d|n, d<n} Phi_d(2)  mod p
    """
    phi = [0] * (N + 1)
    for n in range(1, N + 1):
        num = (pow(2, n, mod) - 1) % mod
        den = 1
        for d in divisors(n):
            if d < n:
                den = den * phi[d] % mod
        phi[n] = num * pow(den, mod - 2, mod) % mod
    return phi, sum(phi[1:]) % mod

# --- Method 2: Mobius product formula ---
def solve_mobius(N: int, mod: int):
    """Compute Phi_n(2) mod p using the Mobius product formula.

    Phi_n(2) = prod_{d|n} (2^d - 1)^{mu(n/d)}  mod p
    """
    mu = mobius_sieve(N)
    pow2m1 = [0] * (N + 1)  # 2^d - 1 mod p
    for d in range(1, N + 1):
        pow2m1[d] = (pow(2, d, mod) - 1) % mod

    phi = [0] * (N + 1)
    for n in range(1, N + 1):
        result = 1
        for d in divisors(n):
            m = mu[n // d]
            if m == 1:
                result = result * pow2m1[d] % mod
            elif m == -1:
                result = result * pow(pow2m1[d], mod - 2, mod) % mod
            # m == 0 contributes factor 1
        phi[n] = result
    return phi, sum(phi[1:]) % mod

# --- Method 3: Exact computation (small n, for verification) ---
def cyclotomic_exact(N: int) -> list:
    """Compute exact Phi_n(2) for n = 1..N (no modular reduction)."""
    phi = [0] * (N + 1)
    for n in range(1, N + 1):
        num = 2**n - 1
        den = 1
        for d in divisors(n):
            if d < n:
                den *= phi[d]
        phi[n] = num // den
    return phi

# --- Compute and verify ---
N = 500

phi_rec, ans_rec = solve_recurrence(N, MOD)
phi_mob, ans_mob = solve_mobius(N, MOD)

assert ans_rec == ans_mob, f"Mismatch: {ans_rec} != {ans_mob}"

# Verify small cases with exact computation
phi_exact = cyclotomic_exact(50)
expected = {1: 1, 2: 3, 3: 7, 4: 5, 5: 31, 6: 3, 7: 127, 8: 17, 10: 11, 12: 13}
for n, val in expected.items():
    assert phi_exact[n] == val, f"Phi_{n}(2) = {phi_exact[n]}, expected {val}"

# Verify fundamental identity: prod_{d|n} Phi_d(2) = 2^n - 1
for n in range(1, 51):
    product = 1
    for d in divisors(n):
        product *= phi_exact[d]
    assert product == 2**n - 1, f"Identity failed for n={n}"

# Verify Mobius formula against exact for small n
for n in range(1, 51):
    assert phi_rec[n] == phi_exact[n] % MOD

print(ans_rec)