Luogu

分析

先莫比乌斯反演一下
$$
\begin{aligned}
&\sum_{i=1}^ni^d[\gcd(i,n)=1]\\
=&\sum_{i=1}^ni^d\sum_{x|i,x|n}\mu(x)\\
=&\sum_{x|n}\mu(x)x^d\sum_{i=1}^{n/x}i^d
\end{aligned}
$$
根据常识可以知道 $\sum_{i=1}^{n/x} i^d$ 是关于 $n/x$ 的 $d+1$ 次多项式,所以我们不妨设这个多项式的 $i$ 次项系数为 $f_i$,上式变成
$$
\begin{aligned}
&\sum_{x|n}\mu(x)x^d\sum_{i=0}^{d+1}f_i\left(\frac{n}{x}\right)^i\\
=&\sum_{i=0}^{d+1}f_i\sum_{x|n}\mu(x)x^d\left(\frac{n}{x}\right)^i
\end{aligned}
$$
后面那个东西是 $(\mu\cdot \mathbf{id})*\mathbf{id}$,显然是积性函数,所以我们可以对于每个质因子单独计算后乘起来。而对于每个质因子显然只有 $x=1$ 和 $x=p$ 时 $\mu(x)$ 不为 $0$,所以只要算两项即可。

至于 $f_i$,可以求拉格朗日插值公式中每项系数求出,也可以直接高斯消元。

代码

// ====================================
//   author: M_sea
//   website: https://m-sea-blog.com/
// ====================================
#include <bits/stdc++.h>
#define file(x) freopen(#x".in","r",stdin); freopen(#x".out","w",stdout)
using namespace std;
typedef long long ll;

int read() {
    int X=0,w=1; char c=getchar();
    while (c<'0'||c>'9') { if (c=='-') w=-1; c=getchar(); }
    while (c>='0'&&c<='9') X=X*10+c-'0',c=getchar();
    return X*w;
}

const int N=1000+10,M=100+10;
const int mod=1e9+7;
int qpow(int a,int b) { int c=1;
    for (;b;b>>=1,a=1ll*a*a%mod) if (b&1) c=1ll*c*a%mod;
    return c;
}

int d,w,p[N],c[N];
int a[M][M],f[M];

void Gauss(int n) {
    for (int i=0;i<=n;++i) {
        if (!a[i][i]) {
            for (int j=i+1;j<=n;++j)
                if (a[i][j]) { swap(a[i],a[j]); break; }
        }
        int inv=qpow(a[i][i],mod-2);
        for (int j=0;j<=n;++j) {
            if (i==j) continue;
            int t=1ll*a[j][i]*inv%mod;
            for (int k=0;k<=n+1;++k)
                a[j][k]=(a[j][k]-1ll*a[i][k]*t%mod+mod)%mod;
        }
    }
    for (int i=0;i<=n;++i) f[i]=1ll*a[i][n+1]*qpow(a[i][i],mod-2)%mod;
}

int main() {
    d=read(),w=read();
    for (int i=1;i<=w;++i) p[i]=read(),c[i]=read();
    for (int i=0,s=0;i<=d+1;++i) {
        s=(s+qpow(i+1,d))%mod,a[i][d+2]=s;
        for (int j=0;j<=d+1;++j) a[i][j]=qpow(i+1,j);
    }
    Gauss(d+1);
    int ans=0;
    for (int i=0;i<=d+1;++i) {
        int now=1;
        for (int j=1;j<=w;++j) {
            int n=qpow(p[j],c[j]),ni=qpow(p[j],c[j]-1);
            now=1ll*now*(qpow(n,i)-1ll*qpow(p[j],d)*qpow(ni,i)%mod+mod)%mod;
        }
        ans=(ans+1ll*f[i]*now)%mod;
    }
    printf("%d\n",ans);
    return 0;
}
最后修改:2020 年 08 月 20 日 04 : 09 PM