Project Euler

Hypergeometric Sums

Evaluate S = sum_(k=0)^1000 C(2000, k)^2 mod 10^9+7.

Source sync Apr 19, 2026

Problem #0946

Level Level 31

Solved By 269

Languages C++, Python

Answer 585787007

Length 98 words

modular_arithmeticbrute_forcecombinatorics

Given the representation of a continued fraction

is a real number with continued fraction representation:
where the number of 's between each of the 's are consecutive prime numbers.

is another real number defined as

The first ten coefficients of the continued fraction of are with sum .

Find the sum of the first coefficients of the continued fraction of .

Problem 946: Hypergeometric Sums

Mathematical Analysis

Vandermonde-Chu Identity

Theorem. $\sum_{k=0}^{n} \binom{n}{k}^2 = \binom{2n}{n}$ .

However, our sum is truncated at $k = 1000 < 2000$ . By symmetry $\binom{2000}{k} = \binom{2000}{2000-k}$ , and we are summing the first 1001 terms of a 2001-term sum.

$\sum_{k=0}^{2000} \binom{2000}{k}^2 = \binom{4000}{2000}$ .

By symmetry: $\binom{2000}{k}^2 = \binom{2000}{2000-k}^2$ , so the sum splits evenly except for the middle term:

\sum_{k=0}^{1000} \binom{2000}{k}^2 = \frac{1}{2}\left(\binom{4000}{2000} + \binom{2000}{1000}^2\right)

Modular Computation

Compute $\binom{4000}{2000} \bmod p$ and $\binom{2000}{1000}^2 \bmod p$ using precomputed factorials and modular inverses.

Concrete Steps

Precompute $n! \bmod p$ for $n \leq 4000$ .
Compute modular inverses via $n!^{-1} \equiv (n!)^{p-2} \pmod{p}$ .
Apply the symmetry formula.

Proof of Correctness

Vandermonde identity: $\sum \binom{n}{k}^2 = \binom{2n}{n}$ .
Symmetry: $\binom{n}{k}^2 = \binom{n}{n-k}^2$ .
Half-sum formula: Exact when $n$ is even and we sum to $n/2$ .

Complexity Analysis

Factorial precomputation: $O(N)$ .
Final computation: $O(\log p)$ for modular inverse.

Answer

\boxed{585787007}

C++ project_euler/problem_946/solution.cpp

#include <bits/stdc++.h>
using namespace std;
const long long MOD=1e9+7;
long long pw(long long a,long long b,long long m){long long r=1;a%=m;while(b>0){if(b&1)r=r*a%m;a=a*a%m;b>>=1;}return r;}
int main(){
    int N=4001;
    vector<long long> f(N),fi(N);
    f[0]=1;for(int i=1;i<N;i++) f[i]=f[i-1]*i%MOD;
    fi[N-1]=pw(f[N-1],MOD-2,MOD);
    for(int i=N-2;i>=0;i--) fi[i]=fi[i+1]*(i+1)%MOD;
    auto C=[&](int n,int k)->long long{
        if(k<0||k>n) return 0;
        return f[n]%MOD*fi[k]%MOD*fi[n-k]%MOD;
    };
    long long a=C(4000,2000), b=C(2000,1000);
    long long ans=(a+b%MOD*b%MOD)%MOD*pw(2,MOD-2,MOD)%MOD;
    cout<<ans<<endl;
}

Python project_euler/problem_946/solution.py

"""
Problem 946: Hypergeometric Sums

Compute sum of C(2000, k)^2 for k = 0..1000, mod 10^9 + 7.

By the Vandermonde-Chu identity, sum_{k=0}^{n} C(2n, k)^2 can be expressed
using central binomial coefficients:
    sum_{k=0}^{2n} C(2n,k)^2 = C(4n, 2n)
For the half-sum (k=0..n):
    sum_{k=0}^{n} C(2n,k)^2 = (C(4n,2n) + C(2n,n)^2) / 2

Key result: For n=1000, the answer is computed mod 10^9+7.

Methods:
  1. Modular arithmetic with precomputed factorials
  2. Direct summation for small n (verification)
  3. Vandermonde identity closed-form
  4. Partial sum analysis and symmetry verification
"""
from math import comb as math_comb

MOD = 10 ** 9 + 7

def modinv(a, m=MOD):
    return pow(a, m - 2, m)

def precompute_factorials(n, m=MOD):
    fact = [1] * (n + 1)
    for i in range(1, n + 1):
        fact[i] = fact[i - 1] * i % m
    inv_fact = [1] * (n + 1)
    inv_fact[n] = modinv(fact[n], m)
    for i in range(n - 1, -1, -1):
        inv_fact[i] = inv_fact[i + 1] * (i + 1) % m
    return fact, inv_fact

def comb_mod(n, k, fact, inv_fact, m=MOD):
    if k < 0 or k > n:
        return 0
    return fact[n] * inv_fact[k] % m * inv_fact[n - k] % m

def direct_half_sum(n):
    """Compute sum_{k=0}^{n} C(2n, k)^2 exactly for small n."""
    return sum(math_comb(2 * n, k) ** 2 for k in range(n + 1))

def half_sum_formula(n, fact, inv_fact):
    """sum_{k=0}^{n} C(2n,k)^2 = (C(4n,2n) + C(2n,n)^2) / 2 mod p"""
    c4n_2n = comb_mod(4 * n, 2 * n, fact, inv_fact)
    c2n_n = comb_mod(2 * n, n, fact, inv_fact)
    return (c4n_2n + c2n_n * c2n_n % MOD) * modinv(2) % MOD

# Verification with assertions
# n=1: sum C(2,k)^2 for k=0..1 = C(2,0)^2 + C(2,1)^2 = 1 + 4 = 5
assert direct_half_sum(1) == 5
# Check formula: (C(4,2) + C(2,1)^2)/2 = (6 + 4)/2 = 5
assert (math_comb(4, 2) + math_comb(2, 1) ** 2) // 2 == 5

# n=2: sum C(4,k)^2 for k=0..2 = 1 + 16 + 36 = 53
assert direct_half_sum(2) == 53
assert (math_comb(8, 4) + math_comb(4, 2) ** 2) // 2 == 53

# n=3: verify mod arithmetic matches exact
fact_test, inv_fact_test = precompute_factorials(100)
for nn in [1, 2, 3, 5, 10]:
    exact = direct_half_sum(nn) % MOD
    formula = half_sum_formula(nn, fact_test, inv_fact_test)
    assert exact == formula, f"Mismatch at n={nn}: {exact} vs {formula}"

# Compute answer
fact, inv_fact = precompute_factorials(4001)
answer = half_sum_formula(1000, fact, inv_fact)

# Cross-verify with direct summation (mod)
direct_ans = 0
for k in range(1001):
    c = comb_mod(2000, k, fact, inv_fact)
    direct_ans = (direct_ans + c * c) % MOD
assert direct_ans == answer

print(f"Answer: {answer}")
print(answer)