一个人的数论：莫比乌斯反演，伯努利数

在做多项式的时候发现自己不会伯努利数，被推荐回来做这道题。

然后就发现这真的是个大神题，而且有一些思路很重要。

以后不能再漏题丢知识点了不然被迫填坑是真的难受。

 

题意式子：

$ \sum\limits_{i=1}^{N} \left[ gcd \left( i,N \right) =1 \right] i^d $

经典反演：

$=\sum\limits_{i=1}^{N} i^d \sum\limits_{j|i \ and \ j|N} \mu(j)$

$=\sum\limits_{j|N} \sum\limits_{i=1}^{\frac{N}{j}} (ij)^d \mu(j)$

$=\sum\limits_{j|N} \mu(j) j^d \sum\limits_{i=1}^{\frac{N}{j}}i^d$

后面是自然数幂和的公式。如果原来没有听说过伯努利数，那么你可以凭空猜测下列结论。

如果听说过伯努利数，那么你就知道对于确定的d，$\sum\limits_{i=1}^{n} i^d =\sum\limits_{i=0}^{d+1} a_i \times n^{i}$

其中$a_i$表示一个特定的系数，不随n变化而变化。事实证明这是正确的，我不会理论证明。

这样的话就得到了一个关于n的d+1次多项式，把求和的复杂度由与n相关变为了与d相关。（伯努利数的应用）

伯努利数是有公式的：$\sum\limits_{i=1}^{n} i^d = \frac{1}{d+1} \sum\limits_{i=0}^{d+1} C_{d+1}^{i} B_i n^{d+1-i} $

这样的话只要我们求解出伯努力数就可以得到上面的所有$a_i$值。

但是首先我们什么也不会，其次这题的数据范围很小，所以我们可以选择将d+1个x值带入上面的多项式进行瓜丝消元，就可以求解出所有$a_i$值。

那么自然数幂和的问题（假如）我们已经解决了。继续化简原式：

$=\sum\limits_{j|N} \mu(j) j^d \sum\limits_{i=0}^{d+1}(\frac{N}{j})^i \times a_i$

$=\sum\limits_{i=0}^{d+1} a_i \sum\limits_{j|N} \mu(j) \times j^d \times  (\frac{N}{j})^i$

$=\sum\limits_{i=0}^{d+1} a_i \times N^i \sum\limits_{j|N} \mu(j) \times j^{d-i}$

观察后面那个和式，它是两个积性函数相乘，依然是积性函数，而题目已经给出了质因数分解，这很方便。

再从含义出发，依次考虑每个质因子，只有在质因子的次数为0或1时对答案有贡献，其余时候$\mu(j)$均为0。

所以只用考虑1与p的和式值，然后因为是积性函数，所以对于每个质因子把结果相乘。

式子的最终形态：

$=\sum\limits_{i=0}^{d+1} a_i \times N^i \prod\limits_{j=1}^{w} (1-p_j^{d-i})$

确实是痛苦AC。感谢miku大神的指导。

#include<cstdio>
#define int long long
#define mod 1000000007
int d,w,fac[105],inv[105],N=1,p[1005],t[1005],mx[102][102],ans;
int pow(int b,int t,int a=1){for(;t;t>>=1,b=b*b%mod)if(t&1)a=a*b%mod;return a;}
void Gauss(){
    for(int i=1;i<=d+1;++i){
        int x=pow(mx[i][i],mod-2);
        for(int j=i;j<=d+2;++j)mx[i][j]=mx[i][j]*x%mod;
        for(int j=i+1;j<=d+1;++j)for(int k=d+2;k>=i;--k)mx[j][k]=(mx[j][k]-mx[i][k]*mx[j][i]%mod+mod)%mod;
    }
    for(int i=d+1;i;--i)for(int j=i-1;j;--j)mx[j][d+2]=(mx[j][d+2]-mx[j][i]*mx[i][d+2]%mod+mod)%mod;
}
signed main(){
    scanf("%lld%lld",&d,&w);fac[0]=1;
    for(int i=1;i<=d;++i)fac[i]=fac[i-1]*i%mod;
    inv[d]=pow(fac[d],mod-2);
    for(int i=d-1;~i;--i)inv[i]=inv[i+1]*(i+1)%mod;
    for(int i=1;i<=w;++i)scanf("%lld%lld",&p[i],&t[i]),N=N*pow(p[i],t[i])%mod;
    for(int i=1;i<=d+1;++i)mx[i][d+2]=(mx[i-1][d+2]+pow(i,d))%mod,mx[i][1]=i;
    for(int i=1;i<=d+1;++i)for(int j=2;j<=d+1;++j)mx[i][j]=mx[i][j-1]*i%mod;
    Gauss();
    for(int i=1;i<=d+1;++i){
        int base=mx[i][d+2]*pow(N,i)%mod;
        for(int j=1;j<=w;++j)base=base*(mod+1-pow(p[j],d-i+mod-1))%mod;
        ans=(ans+base)%mod;
    }printf("%lld\n",ans);
}