BZOJ3601：一个人的数论（莫比乌斯反演+伯努利数）

题面

题意：给出d和n（n以分解质因数给出），问所有与n互质的数的d次幂之和，即

$$\sum_xx^d[gcd(x,n)==1]$$

套路推♂倒

$$=\sum_{i|n}\mu(i)*i^d*\sum_{x=1}^{\frac{n}{i}}x^d$$
fa♂现最右边是个幂和，设为 $h(\frac{n}{i})$，整个就是一个狄利克雷卷积。

根据题目n以分解质因数输入的套路，应该是找到积性函数，然后一个个质因数乘起来
$\mu(i)*i^d$是积性的，但是幂和应该不是积性的。

dyh的ppt告诉我们，幂和是一个多项式，而多项式的每一项显然是积性的
所以我们分开每一项来做，最后求和。

对于$x^i$，系数为$\frac{1}{d+1}*B_{d-i+1}*C_{d+1}^{d-i+1}$，设为$a_i$
（伯努利数递推求出）

则原式为

$$\sum_{i=1}^{d+1}a_i\sum_{c|n}\mu(c)*c^d*(\frac{n}{c})^i$$ 设 $g[n]=\sum_{c|n}\mu(c)*c^d*(\frac{n}{c})^i$，由于积性，仅考虑 $g[p^q]$的值

$$g[p^q]=\sum_{j=0}^q\mu(p^j)*p^{jd}*p^{qi-ji}$$
仅当j取0和1时有贡献
则 $g[p^q]=p^{qi}-p^{qi+d-i}$

我是看到高斯消元才点开这题的，AC了也没想懂这怎么是高斯消元。
刚才下去买零食看到L指导
才想起他曾讲过
幂和如果不记得
就打高斯消元

#include <iostream>
#include <fstream>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <cstdio>
#include <cstdlib>
#include <cstring>

using namespace std;
#define mmst(a, b) memset(a, b, sizeof(a))
#define mmcp(a, b) memcpy(a, b, sizeof(b))

typedef long long LL;

const LL mo=1e9+7;
const LL N=100010;

LL d,w;
LL a[N],p[N];
LL Ijc[N],I[N],jc[N];
LL B[N],A[N];
LL ans;

LL cheng(LL a,LL b)
{
    LL res=1;
    for(;b;b>>=1,a=a*a%mo)
    if(b&1)
    res=res*a%mo;
    return res;
}

int main()
{
    cin>>d>>w;
    for(int i=1;i<=w;i++)
    scanf("%lld%lld",&p[i],&a[i]);

    I[1]=Ijc[0]=jc[0]=1;

    for(int i=2;i<N;i++)
    I[i]=I[mo%i]*(mo-mo/i)%mo;

    for(int i=1;i<N;i++)
    jc[i]=jc[i-1]*i%mo,Ijc[i]=Ijc[i-1]*I[i]%mo;

    B[0]=1;
    for(int i=1;i<=2000;i++)
    {
        for(int j=0;j<i;j++)
        B[i]=(B[i]+jc[i+1]*Ijc[j]%mo*Ijc[i+1-j]%mo*B[j]%mo)%mo;
        B[i]=(mo-B[i])%mo;
        B[i]=B[i]*I[i+1]%mo;
    }

    for(int i=0;i<=d;i++)
    A[d+1-i]=I[d+1]*jc[d+1]%mo*Ijc[i]%mo*Ijc[d+1-i]%mo*B[i]%mo;

    for(int i=1;i<=d+1;i++)
    {
        LL g=1;
        for(int j=1;j<=w;j++)
        {
            LL tu=cheng(cheng(p[j],a[j]),i);
            LL hy=(tu-tu*cheng(p[j],d)%mo*cheng(cheng(p[j],i),mo-2)%mo+mo)%mo;
            g=g*hy%mo;
        }
        ans=(ans+g*A[i])%mo;
    }

    cout<<ans<<endl;

    return 0;
}