P6271 [湖北省队互测2014]一个人的数论

为了方便计算，这里用变量名 ${k}$ 表示题目中的 ${d}$

求解：

${\sum_{i=1}^{n}{[gcd(i,n)=1]i^k}}$

则原式：

$={\sum_{i=1}^{n}{\varepsilon (gcd(i,n)=1)i^k}}$

$={\sum_{i=1}^{n}\sum_{d|i,d|n}{\mu(d)*i^k}}$

$={\sum_{d|n}d^k*\mu(d)*\sum_{i=1}^{⌊\frac{n}{d}⌋}i^k}$

自然数幂和是一个 ${k+1}$ 次多项式，可以用伯努利数算出系数，不妨记为 ${f_i}$

由伯努利公式得： ${f_{k+1-i}=\frac{\binom{k+1}{i}B_i}{k+1}}$

则原式：

$={\sum_{d|n}d^k*\mu(d)*\sum_{i=1}^{k+1}f_i*(\frac{n}{d})^i}$

$=\sum_{i=1}^{k+1}f_i*n^i*{\sum_{d|n}\mu(d)*d^{k-i}}$

令 ${g(x)=\sum_{d|n}\mu(d)*d^{k-c}}$ ，易知 ${g}$ 是积性函数，令 $y=p_i^{a_i}$ ，则有

${g(y)=u(p_i^0)*(p_i^0)^{k-c}+u(p_i^1)*(p_i^1)^{k-c}+u(p_i^2)*(p_i^2)^{k-c}+...+u(p_i^{a_i})*(p_i^{a_i})^{k-c}}$

${=1-p_i^{k-c}+0+...+0}$

则 ${g(x)=g(p_1^{a_1})*g(p_2^{a_2})*...**g(p_w^{a_w})=\prod_{i=1}^{w}(1-p_i^{k-c})}$

Code:

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn=107,mod=1e9+7;
ll qpow(ll a,ll b) {
    ll res=1;
    while(b) {
        if(b&1) res=res*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return res;
}
ll inv[maxn],fac[maxn];
int B[maxn],f[maxn];
ll C(int n,int m) {
    return fac[n]*inv[m]%mod*inv[n-m]%mod;
}
void Init_Bernoulli(const int k) {
    B[0]=1;
    for(int i=1;i<=k;++i) {
        ll sum(0);
        for(int j=0;j<i;++j) sum=(sum+C(i+1,j)*B[j])%mod;
        B[i]=(-sum*qpow(C(i+1,i),mod-2)%mod+mod)%mod;
    }
    for(int i=0;i<=k;++i) f[k+1-i]=C(k+1,i)*B[i]%mod*qpow(k+1,mod-2)%mod;
}
void init(const int k) {
    fac[ 0 ] = 1;
    for( int i = 1 ; i <= k ; i ++ ) fac[ i ] = fac[ i - 1 ]*i%mod;
    inv[ k ] = qpow(fac[k],mod-2);
    for( int i = k ; i >= 1 ; i -- ) inv[ i - 1 ] = inv[ i ]*i%mod;
}
int k,w,p[1005];
int sum(int pw) {
    if(pw<0) pw+=mod-1;
    ll res=1;
    for(int i=1;i<=w;++i) res=res*(1-qpow(p[i],pw))%mod;
    return (res+mod)%mod;
}
signed main() {
    cin.sync_with_stdio(false), cin.tie(nullptr),cout.tie(nullptr);
    cin>>k>>w;
    ll ans(0),n=1,x=1;
    for(int i=1,y;i<=w;++i) {
        cin>>p[i]>>y;
        n=n*qpow(p[i],y)%mod;
    }
    init(k+1);
    Init_Bernoulli(k);
    for(int i=1;i<=k+1;++i) {
        x=x*n%mod;
        ans=(ans+f[i]*x%mod*sum(k-i))%mod;
//        cout<<f[i]<<' '<<x<<' '<<sum(k-i)<<'\n';
    }
    cout<<ans<<'\n';
    return 0;
}