BZOJ 3601: 一个人的数论

小清新数学题，全是套路真是舒适无比的说

我们先化一下式子：

\[ans=\sum_{i=1}^n [\gcd(i,n)==1] i^d \]

\[=\sum_{i=1}^n \sum_{t|i,t|n} \mu(t) \times i^d \]

\[=\sum_{t|n}\mu(t) \sum_{i=1}^{\frac{n}{t}} (it)^d \]

\[=\sum_{t|n}\mu(t)\times t^d\times \sum_{i=1}^{\frac{n}{t}} i^d \]

后面的那个\(\sum_{i=1}^{\frac{n}{t}} i^d\)我们都不要太熟悉的说，自然数幂和

然后我们开始冷静分析一下，\(n\)很大必然是不能直接求的，也不能直接枚举\(n\)的约数因为实在太多了

看起来GG了，我们反过来审视一下题目，发现它给我们的是\(n\)的质因数分解，看到这个我们再回头看一下式子

\(\mu(t)\)是积性函数，\(t^d\)是完全积性函数，那么如果后面的自然数幂和换成一个积性函数的话我们就可以对每一个\(p_i\)单独计算出答案然后乘起来即可

在拉格朗日插值中我们得知\(\sum_{i=1}^n i^d\)是一个关于\(n\)的\(d+1\)次多项式，如果我们假设这个多项式的系数是\(a\)的话，则有：

\[\sum_{i=1}^n i^d=\sum_{i=1}^{d+1} a_i\times n^i \]

那么假设我们知道了\(a\)，原式子就变成了：

\[=\sum_{i=1}^{d+1} a_i \times \sum_{t|n}\mu(t)\times t^d\times (\frac{n}{t})^i \]

后面的那个是什么，狄利克雷卷积啊！我们知道它也是个积性函数，然后就可以算了

等下你说就算是\(p_i^{\alpha_i}\)的约数也有\(\alpha_i\)个啊，不能一个个枚举。想一下\(\mu(t)\)的定义吧，当\(t\)有平方质因子是\(\mu(t)=0\)

那么也就意味着我们只用考虑\(t=1\)和\(t=p_i\)两种情况了，顺利解决

最后对于那个\(a\)的计算显然可以用拉格朗日插值的插系数的\(O(n^2)\)写法（不过需要手写多项式除法，比较麻烦），我们看到数据范围\(d\le 100\)直接暴力上高斯消元就可以了

#include<cstdio>
#include<iostream>
#define RI register int
#define CI const int&
using namespace std;
const int N=1005,mod=1e9+7;
int w,d,p[N],alp[N],a[N][N],val[N],ans;
inline void inc(int& x,CI y)
{
    if ((x+=y)>=mod) x-=mod;
}
inline void dec(int& x,CI y)
{
    if ((x-=y)<0) x+=mod;
}
inline int quick_pow(int x,int p=mod-2,int mul=1)
{
    for (;p;p>>=1,x=1LL*x*x%mod) if (p&1) mul=1LL*mul*x%mod; return mul;
}
inline void Gauss(CI n)
{
    RI i,j,k; for (i=1;i<=n;++i)
    {
        for (j=i;!a[j][i];++j); swap(a[i],a[j]); swap(val[i],val[j]);
        for (j=i+1;j<=n;++j)
        {
            int tp=1LL*a[j][i]*quick_pow(a[i][i])%mod; dec(val[j],1LL*val[i]*tp%mod);
            for (k=i;k<=n;++k) dec(a[j][k],1LL*a[i][k]*tp%mod);
        }
    }
    for (i=n;i;--i)
    {
        for (j=i+1;j<=n;++j) dec(val[i],1LL*val[j]*a[i][j]%mod);
        val[i]=1LL*val[i]*quick_pow(a[i][i])%mod;
    }
}
inline int calc(CI x,int ret=1)
{
    for (RI i=1;i<=w;++i)
    {
        int tp=quick_pow(p[i],1LL*alp[i]*x%(mod-1));
        if (alp[i]) dec(tp,1LL*quick_pow(p[i],d)*quick_pow(p[i],1LL*(alp[i]-1)*x%(mod-1))%mod);
        ret=1LL*ret*tp%mod;
    }
    return ret;
}
int main()
{
    RI i,j; for (scanf("%d%d",&d,&w),i=1;i<=w;++i) scanf("%d%d",&p[i],&alp[i]);
    for (i=1;i<=d+1;++i) for (j=1;j<=d+1;++j)
    a[i][j]=quick_pow(i,j),j<=i&&(inc(val[i],quick_pow(j,d)),0);
    for (Gauss(d+1),i=1;i<=d+1;++i) inc(ans,1LL*val[i]*calc(i)%mod);
    return printf("%d",ans),0;
}