Project Euler 650

题目

Divisors of Binomial Product

Let $B(n) = \displaystyle \prod_{k=0}^n {n \choose k}$, a product of binomial coefficients.

For example, $B(5) = {5 \choose 0} \times {5 \choose 1} \times {5 \choose 2} \times {5 \choose 3} \times {5 \choose 4} \times {5 \choose 5} = 1 \times 5 \times 10 \times 10 \times 5 \times 1 = 2500$.

Let $D(n) = \displaystyle \sum_{d\mid B(n)} d$, the sum of the divisors of $B(n)$.

For example, the divisors of $B(5)$ are $1, 2, 4, 5, 10, 20, 25, 50, 100, 125, 250, 500, 625, 1250$ and $2500$,

so $D(5) = 1 + 2 + 4 + 5 + 10 + 20 + 25 + 50 + 100 + 125 + 250 + 500 + 625 + 1250 + 2500 = 5467$.

Let $S(n) = \displaystyle \sum_{k=1}^n D(k)$.

You are given $S(5) = 5736$, $S(10) = 141740594713218418$ and $S(100) \bmod 1\,000\,000\,007 = 332792866$.

Find $S(20\,000) \bmod 1\,000\,000\,007$.

解决方案

根据组合数的定义，不难计算出

$$B(n)=\dfrac{(n!)^{n+1}}{\prod_{i=0}^n (i!)^2}$$

那么，不难将$B(n)$可以写成一个递推式：

$$B(n)=B(n-1)\cdot \dfrac{n^n}{n!}$$

直接维护$B(n)$的分解质因数$B(n)=p_1^{e_1}\cdot p_2^{e_2}\cdot p_1^{e_3}\dots$中$e_1,e_2,e_3\dots$的值。每一次求好之后利用因数和定理$\sigma(B(n))=\prod\dfrac{p_i^{e_i+1}-1}{p_i-1}$和快速幂算法进行计算。

代码

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=20000;
bool vis[N+4];
int p[24],e[24],o=0;
int pr[N/10+1000],m=0;
ll mod=1e9+7;

void fact(int n){
    o=0;
    for(int i=1;i<=m&&pr[i]*pr[i]<=n;i++)
        if(n%pr[i]==0){
            p[++o]=pr[i];e[o]=0;
            for(;n%pr[i]==0;n/=pr[i],++e[o]);
        }
    if(n!=1){
        p[++o]=n;e[o]=1;
    }
}
ll qpow(ll n,ll m){
    ll a=1;
    for(;m;m>>=1){
        if(m&1) a=a*n%mod;
        n=n*n%mod;
    }
    return a;
}
int fac_fact[N+4],now_fact[N+4];
ll inv_p_1[N+4];
int main(){
    for(int i=2;i<=N;i++){
        if(vis[i]) continue;
        pr[++m]=i;
        for(ll j=1ll*i*i;j<=N;j+=i)
            vis[j]=true;
    }
    for(int i=1;i<=m;i++)
        inv_p_1[pr[i]]=qpow(pr[i]-1,mod-2);
    ll ans=0;
    for(int i=1;i<=N;i++){
        fact(i);
        for(int j=1;j<=o;j++) {
            fac_fact[p[j]]+=e[j];
            now_fact[p[j]]+=e[j]*i;
        }
        ll s=1;
        for(int j=1;j<=m&&pr[j]<=i;j++){
            now_fact[pr[j]]-=fac_fact[pr[j]];
            s=s*(qpow(pr[j],now_fact[pr[j]]+1)-1)%mod*inv_p_1[pr[j]]%mod;
        }
        ans=(ans+s)%mod;
    }
    printf("%lld\n",ans);
}