Problem 498: Remainder of Polynomial Division

View on Project Euler

Project Euler Problem 498 Solution

EulerSolve provides an optimized solution for Project Euler Problem 498, Remainder of Polynomial Division, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Let \(R_{n,m}(x)\) be the remainder when \(x^n\) is divided by \((x-1)^m\). For integers \(0 \le d \lt m \le n\), the quantity of interest is the absolute value of the coefficient of \(x^d\) in that remainder. The real input is enormous, so the solution cannot expand polynomials term by term; it must turn the coefficient into a direct combinatorial formula and then evaluate that formula modulo the prime \(p=999999937\). Mathematical Approach The key observation is that the remainder modulo \((x-1)^m\) is exactly the Taylor expansion of \(x^n\) around \(x=1\), truncated after degree \(m-1\). Step 1: Write the remainder as a truncated expansion around \(x=1\) Since \(x=1+(x-1)\), the binomial theorem gives $$x^n=(1+(x-1))^n=\sum_{k=0}^{n}\binom{n}{k}(x-1)^k.$$ When we divide by \((x-1)^m\), every term with \(k \ge m\) is a multiple of \((x-1)^m\). Therefore the unique remainder of degree less than \(m\) is $$R_{n,m}(x)=\sum_{k=0}^{m-1}\binom{n}{k}(x-1)^k.$$ This already explains why the code never performs polynomial long division: the remainder has a closed form immediately....

Detailed mathematical approach

Problem Summary

Let \(R_{n,m}(x)\) be the remainder when \(x^n\) is divided by \((x-1)^m\). For integers \(0 \le d \lt m \le n\), the quantity of interest is the absolute value of the coefficient of \(x^d\) in that remainder. The real input is enormous, so the solution cannot expand polynomials term by term; it must turn the coefficient into a direct combinatorial formula and then evaluate that formula modulo the prime \(p=999999937\).

Mathematical Approach

The key observation is that the remainder modulo \((x-1)^m\) is exactly the Taylor expansion of \(x^n\) around \(x=1\), truncated after degree \(m-1\).

Step 1: Write the remainder as a truncated expansion around \(x=1\)

Since \(x=1+(x-1)\), the binomial theorem gives

$$x^n=(1+(x-1))^n=\sum_{k=0}^{n}\binom{n}{k}(x-1)^k.$$

When we divide by \((x-1)^m\), every term with \(k \ge m\) is a multiple of \((x-1)^m\). Therefore the unique remainder of degree less than \(m\) is

$$R_{n,m}(x)=\sum_{k=0}^{m-1}\binom{n}{k}(x-1)^k.$$

This already explains why the code never performs polynomial long division: the remainder has a closed form immediately.

Step 2: Extract the coefficient of \(x^d\)

Inside one term \((x-1)^k\), the coefficient of \(x^d\) is

$$[x^d](x-1)^k=\binom{k}{d}(-1)^{k-d},\qquad k \ge d.$$

So if \(a_d\) denotes the coefficient of \(x^d\) in the remainder, then

$$a_d=\sum_{k=d}^{m-1}\binom{n}{k}\binom{k}{d}(-1)^{k-d}.$$

The problem asks for \(|a_d|\), not the signed coefficient itself.

Step 3: Separate the fixed part from the alternating sum

Use the standard identity

$$\binom{n}{k}\binom{k}{d}=\binom{n}{d}\binom{n-d}{k-d}.$$

After substituting \(t=k-d\), the coefficient becomes

$$a_d=\binom{n}{d}\sum_{t=0}^{m-d-1}(-1)^t\binom{n-d}{t}.$$

Now the whole problem is reduced to one alternating partial sum of binomial coefficients.

Step 4: Collapse the alternating binomial sum

For \(0 \le r \lt N\), the identity

$$\sum_{t=0}^{r}(-1)^t\binom{N}{t}=(-1)^r\binom{N-1}{r}$$

applies. With \(N=n-d\) and \(r=m-d-1\), we obtain

$$a_d=(-1)^{m-d-1}\binom{n}{d}\binom{n-d-1}{m-d-1}.$$

Hence the absolute value is

$$\boxed{|a_d|=\binom{n}{d}\binom{n-d-1}{m-d-1}.}$$

This closed formula is the central mathematical fact used by the implementation.

Step 5: Evaluate the formula modulo the prime \(p\)

The exact coefficient is far too large to build directly, so the computation is done in \(\mathbb{F}_p\) with

$$p=999999937.$$

For a prime modulus, Lucas's theorem says that if

$$N=\sum_i N_i p^i,\qquad K=\sum_i K_i p^i,$$

then

$$\binom{N}{K}\equiv \prod_i \binom{N_i}{K_i}\pmod{p}.$$

Each digit-level binomial is small enough to evaluate with the multiplicative formula

$$\binom{u}{v}=\frac{u(u-1)\cdots(u-v+1)}{v!},$$

and the division by \(v!\) is performed modulo \(p\) using Fermat's little theorem, namely

$$q^{-1}\equiv q^{p-2}\pmod{p} \qquad (q \not\equiv 0 \pmod{p}).$$

So the final answer is

$$|a_d| \bmod p \equiv \binom{n}{d}\binom{n-d-1}{m-d-1}\pmod{p}.$$

Worked Example: \((n,m,d)=(6,3,1)\)

The truncated expansion is

$$R_{6,3}(x)=\binom{6}{0}+\binom{6}{1}(x-1)+\binom{6}{2}(x-1)^2.$$

Substituting the values gives

$$R_{6,3}(x)=1+6(x-1)+15(x-1)^2=15x^2-24x+10.$$

The coefficient of \(x^1\) is \(-24\), so the required absolute value is \(24\).

The closed formula gives exactly the same result:

$$\binom{6}{1}\binom{6-1-1}{3-1-1}=\binom{6}{1}\binom{4}{1}=6\cdot 4=24.$$

How the Code Works

The implementation first handles the invalid range: if \(d \ge m\) or \(m > n\), the desired coefficient is zero. Otherwise it computes the two binomial factors from the closed formula separately modulo \(p\), and then multiplies them together modulo \(p\).

To evaluate a binomial coefficient modulo the prime, the C++, Python, and Java implementations apply Lucas's theorem digit by digit in base \(p\). For each digit pair they form the numerator and denominator multiplicatively, replace division by a modular inverse obtained from fast exponentiation, and accumulate the digit contributions into the full Lucas product.

The C++ implementation also checks the formula on small cases by comparing it with direct coefficient expansion of the truncated remainder, which confirms both the combinatorial identity and the modular computation.

Complexity Analysis

Let \(L=\lfloor \log_p n \rfloor + 1\), the number of base-\(p\) digits. Lucas's theorem reduces each large binomial to \(L\) digit-level binomials. If a digit pair is \((u,v)\), the multiplicative evaluation costs \(O(\min(v,u-v))\) modular multiplications, plus \(O(\log p)\) time for the modular inverse by fast exponentiation. Therefore the overall running time is

$$O\left(\sum_{i=0}^{L-1}\min(K_i,N_i-K_i)+L\log p\right)$$

for each large binomial, with only \(O(1)\) extra memory. In the actual problem, \(p\) is close to \(10^9\), so the number of base-\(p\) digits is extremely small.

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=498
  2. Binomial theorem: Wikipedia — Binomial theorem
  3. Binomial coefficient: Wikipedia — Binomial coefficient
  4. Lucas's theorem: Wikipedia — Lucas's theorem
  5. Fermat's little theorem: Wikipedia — Fermat's little theorem

Problem 498 source code

C++

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

namespace {

using u64 = std::uint64_t;
using u128 = __uint128_t;

constexpr u64 kPrime = 999'999'937ULL;

u64 mod_pow(u64 base, u64 exp, const u64 mod) {
    u64 result = 1ULL % mod;
    u64 cur = base % mod;
    u64 e = exp;
    while (e > 0ULL) {
        if (e & 1ULL) {
            result = static_cast<u64>((static_cast<u128>(result) * cur) % mod);
        }
        cur = static_cast<u64>((static_cast<u128>(cur) * cur) % mod);
        e >>= 1ULL;
    }
    return result;
}

u64 binom_mod_small_prime_digit(u64 n, u64 k, const u64 prime) {
    if (k > n) {
        return 0ULL;
    }
    k = std::min(k, n - k);
    if (k == 0ULL) {
        return 1ULL;
    }

    u64 numerator = 1ULL;
    u64 denominator = 1ULL;
    for (u64 i = 1ULL; i <= k; ++i) {
        numerator = static_cast<u64>((static_cast<u128>(numerator) * (n - k + i)) % prime);
        denominator = static_cast<u64>((static_cast<u128>(denominator) * i) % prime);
    }
    const u64 inv_denominator = mod_pow(denominator, prime - 2ULL, prime);
    return static_cast<u64>((static_cast<u128>(numerator) * inv_denominator) % prime);
}

u64 binom_mod_prime_lucas(u64 n, u64 k, const u64 prime) {
    if (k > n) {
        return 0ULL;
    }
    u64 result = 1ULL;
    u64 nn = n;
    u64 kk = k;
    while (nn > 0ULL || kk > 0ULL) {
        const u64 ni = nn % prime;
        const u64 ki = kk % prime;
        if (ki > ni) {
            return 0ULL;
        }
        const u64 digit = binom_mod_small_prime_digit(ni, ki, prime);
        result = static_cast<u64>((static_cast<u128>(result) * digit) % prime);
        nn /= prime;
        kk /= prime;
    }
    return result;
}

u64 coefficient_mod(const u64 n, const u64 m, const u64 d, const u64 prime) {
    if (d >= m || m > n) {
        return 0ULL;
    }
    const u64 first = binom_mod_prime_lucas(n, d, prime);
    const u64 second = binom_mod_prime_lucas(n - d - 1ULL, m - d - 1ULL, prime);
    return static_cast<u64>((static_cast<u128>(first) * second) % prime);
}

u64 binom_exact_small(const u64 n, const u64 k) {
    if (k > n) {
        return 0ULL;
    }
    const u64 kk = std::min(k, n - k);
    u128 result = 1;
    for (u64 i = 1ULL; i <= kk; ++i) {
        result = (result * static_cast<u128>(n - kk + i)) / static_cast<u128>(i);
    }
    return static_cast<u64>(result);
}

u64 coefficient_abs_exact_small(const u64 n, const u64 m, const u64 d) {
    if (d >= m || m > n) {
        return 0ULL;
    }
    const u64 first = binom_exact_small(n, d);
    const u64 second = binom_exact_small(n - d - 1ULL, m - d - 1ULL);
    return first * second;
}

u64 coefficient_abs_bruteforce_small(const u64 n, const u64 m, const u64 d) {
    if (d >= m || m > n) {
        return 0ULL;
    }

    std::vector<long long> coeffs(static_cast<std::size_t>(m), 0LL);
    for (u64 k = 0ULL; k < m; ++k) {
        const u64 choose_n_k = binom_exact_small(n, k);
        for (u64 j = 0ULL; j <= k; ++j) {
            const long long choose_k_j = static_cast<long long>(binom_exact_small(k, j));
            const long long sign = ((k - j) % 2ULL == 0ULL) ? 1LL : -1LL;
            coeffs[static_cast<std::size_t>(j)] +=
                static_cast<long long>(choose_n_k) * choose_k_j * sign;
        }
    }
    const long long value = coeffs[static_cast<std::size_t>(d)];
    return static_cast<u64>(value >= 0LL ? value : -value);
}

bool run_checkpoints() {
    if (coefficient_abs_exact_small(6ULL, 3ULL, 1ULL) != 24ULL) {
        std::cerr << "Checkpoint failed: C(6,3,1)\n";
        return false;
    }
    if (coefficient_abs_exact_small(100ULL, 10ULL, 4ULL) != 227'197'811'615'775ULL) {
        std::cerr << "Checkpoint failed: C(100,10,4)\n";
        return false;
    }

    for (u64 n = 1ULL; n <= 12ULL; ++n) {
        for (u64 m = 1ULL; m <= n; ++m) {
            for (u64 d = 0ULL; d < m; ++d) {
                const u64 closed_form = coefficient_abs_exact_small(n, m, d);
                const u64 brute = coefficient_abs_bruteforce_small(n, m, d);
                if (closed_form != brute) {
                    std::cerr << "Bruteforce mismatch at n=" << n << ", m=" << m << ", d=" << d
                              << '\n';
                    return false;
                }
                const u64 via_mod = coefficient_mod(n, m, d, kPrime);
                if (via_mod != (closed_form % kPrime)) {
                    std::cerr << "Modulo mismatch at n=" << n << ", m=" << m << ", d=" << d
                              << '\n';
                    return false;
                }
            }
        }
    }
    return true;
}

}  // namespace

int main() {
    if (!run_checkpoints()) {
        return 1;
    }

    constexpr u64 n = 10'000'000'000'000ULL;
    constexpr u64 m = 1'000'000'000'000ULL;
    constexpr u64 d = 10'000ULL;

    const u64 answer = coefficient_mod(n, m, d, kPrime);
    std::cout << answer << '\n';
    return 0;
}

Python

def mod_pow(base, exp, mod):
    return pow(base, exp, mod)

def binom_mod_small_prime_digit(n, k, prime):
    if k > n:
        return 0
    k = min(k, n - k)
    if k == 0:
        return 1

    numerator = 1
    denominator = 1
    for i in range(1, k + 1):
        numerator = (numerator * (n - k + i)) % prime
        denominator = (denominator * i) % prime
        
    inv_denominator = mod_pow(denominator, prime - 2, prime)
    return (numerator * inv_denominator) % prime

def binom_mod_prime_lucas(n, k, prime):
    if k > n:
        return 0
    result = 1
    nn = n
    kk = k
    while nn > 0 or kk > 0:
        ni = nn % prime
        ki = kk % prime
        if ki > ni:
            return 0
        digit = binom_mod_small_prime_digit(ni, ki, prime)
        result = (result * digit) % prime
        nn //= prime
        kk //= prime
    return result

def coefficient_mod(n, m, d, prime):
    if d >= m or m > n:
        return 0
    first = binom_mod_prime_lucas(n, d, prime)
    second = binom_mod_prime_lucas(n - d - 1, m - d - 1, prime)
    return (first * second) % prime

def solve():
    kPrime = 999999937
    n = 10000000000000
    m = 1000000000000
    d = 10000
    ans = coefficient_mod(n, m, d, kPrime)
    return str(ans)

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

Java

public class Euler498 {

    private static final long kPrime = 999999937L;

    private static long modPow(long base, long exp, long mod) {
        long result = 1L % mod;
        long cur = base % mod;
        long e = exp;
        while (e > 0L) {
            if ((e & 1L) != 0L) {
                result = (result * cur) % mod;
            }
            cur = (cur * cur) % mod;
            e >>= 1L;
        }
        return result;
    }

    private static long binomModSmallPrimeDigit(long n, long k, long prime) {
        if (k > n) {
            return 0L;
        }
        k = Math.min(k, n - k);
        if (k == 0L) {
            return 1L;
        }

        long numerator = 1L;
        long denominator = 1L;
        for (long i = 1L; i <= k; ++i) {
            numerator = (numerator * (n - k + i)) % prime;
            denominator = (denominator * i) % prime;
        }
        long invDenominator = modPow(denominator, prime - 2L, prime);
        return (numerator * invDenominator) % prime;
    }

    private static long binomModPrimeLucas(long n, long k, long prime) {
        if (k > n) {
            return 0L;
        }
        long result = 1L;
        long nn = n;
        long kk = k;
        while (nn > 0L || kk > 0L) {
            long ni = nn % prime;
            long ki = kk % prime;
            if (ki > ni) {
                return 0L;
            }
            long digit = binomModSmallPrimeDigit(ni, ki, prime);
            result = (result * digit) % prime;
            nn /= prime;
            kk /= prime;
        }
        return result;
    }

    private static long coefficientMod(long n, long m, long d, long prime) {
        if (d >= m || m > n) {
            return 0L;
        }
        long first = binomModPrimeLucas(n, d, prime);
        long second = binomModPrimeLucas(n - d - 1L, m - d - 1L, prime);
        return (first * second) % prime;
    }

    public static void main(String[] args) {
        long n = 10000000000000L;
        long m = 1000000000000L;
        long d = 10000L;

        long answer = coefficientMod(n, m, d, kPrime);
        System.out.println(answer);
    }
}