【数论】bzoj3601一个人的数论

题目链接
如果这题叫“一只蛤的神论”…

反正窝只会30分乱搞…
test1:暴力水过 test2:欧拉函数 test3:n/2*欧拉函数…

蒟蒻之前从未见过莫比乌斯反演与高斯消元如此奇妙的结合在一起…
所以来瞎bb一下…

明确目标：$ans=\sum_{i=1}^{n}[gcd(i,n)==1]i^d$，然后就开始漫长的治疗公式恐惧症的历程…

利用莫比乌斯反演最基本的运用将$gcd$展开成莫比乌斯函数。

$$\sum_{i=1}^n[gcd(i,n)==1]i^d=\sum_{i=1}^ni^d\sum_{g|i}\mu(g)$$
利用和式的性质将柿子变一下性。
$$=\sum_{g|n}\mu(g)\sum_{i=1}^n[g|i]*i^d$$
此处 $i$的作用与之前不同，之前的$i$用来枚举 $g$的倍数，此处$gi$用来枚举 $g$的倍数。
$$=\sum_{g|n}\mu(g)\sum_{i=1}^{\frac{n}g}(gi)^d=\sum_{g|n}\mu(g)g^d\sum_{i=1}^{\frac{n}g}i^d$$

然后是窝觉得此题最令人惊叹的地方。
令$s(m)=\sum_{i=1}^mi^d$，我们猜测这个柿子是关于$m$的$d+1$次的多项式。这货只能感性的理解一下…

然后我们就可以用高斯消元来求这个多项式的系数…
当然也可以用拉格朗日插值法，理论复杂度比高斯消元优（然而窝写不来…逃…）
就是暴力求出$1$~$d+1$的$s$值，然后待定系数，然后高斯消元…
然后就知道辣$s(m)=\sum_{i=1}^{d+1}a_im^i$

$$ans=\sum_{g|n}g^d\mu(g)\sum_{i=1}^da_i(\frac{n}g)^i$$
$$=\sum_{i=1}^da_i\sum_{g|n}\mu(g)g^{d-i}*n^i$$
观察与 $g$有关的柿子，$\mu(g)$和 $g^{d-i}$都是积性函数，所以 $h_i(n)=\sum_{g|n}\mu(g)g^{d-i}*n^i$也是积性函数。

加之题目中的$n$是以质因数分解的形式给出，考虑$n$的某一个质因数$p$以及其对应的次幂。

$$h_i(p^a)=\sum_{j=0}^a\mu(p^j)p^{j*(d-i)}p^{ai}$$
根据莫比乌斯函数的定义，在这些柿子中只有$\mu(1)$与$\mu(p)$的值非$0$。
$$=\mu(1)p^{ai}+\mu(p)p^{d-i}p^{ai}$$

$$=p^{ai}-p^{d-i}p^{ai}=p^{ai}(1-p^{ai})$$

至此，问题解决完毕。
大概重新理一下解题过程。
1.用待定系数+高斯消元（拉格朗日插值法）求出某多项式的系数。
2.直接上计算式子$ans==\sum_{i=1}^da_i\prod_ph_i(p^a)$

#include <iostream>
#include <cstdio>
#define mod 1000000007
#define LL long long int
using namespace std;

int n, d, p[1005][2];
int sum[105], mat[155][155], A[155];

int power(int a,LL pos)
{
    int ans=1;
    if(pos<0)pos+=mod-1;
    for(;pos;pos>>=1, a=1ll*a*a%mod)
        if(pos&1)ans=1ll*ans*a%mod;
    return ans;
}

void init()
{
    for(int i=1;i<=d+2;++i)
    {
        sum[i]=(power(i,d)+sum[i-1])%mod;
        mat[i-1][d+2]=sum[i], mat[i-1][0]=1;
        for(int j=1;j<=d+1;++j)mat[i-1][j]=1ll*mat[i-1][j-1]*i%mod;
    }

    for(int i=0, j, k, tmp;i<=d+1;++i)
    {
        for(j=i;j<=d+1;++j)if(mat[j][i])break;
        if(i!=j)for(k=0;k<=d+2;++k)swap(mat[i][k],mat[j][k]);
        for(j=0;j<=d+1;++j)
            if(j!=i&&mat[j][i])
            {
                tmp=1ll*mat[j][i]*power(mat[i][i],mod-2)%mod;
                for(k=0;k<=d+2;++k)mat[j][k]=(mat[j][k]-1ll*tmp*mat[i][k]%mod+mod)%mod;
            }
    }

    for(int i=0;i<=d+1;++i)A[i]=1ll*mat[i][d+2]*power(mat[i][i],mod-2)%mod;
}

int main()
{
    scanf("%d%d",&d,&n);
    init();
    for(int i=1;i<=n;++i)scanf("%d%d",&p[i][0],&p[i][1]);
    int ans=0;
    for(int i=1, tmp;i<=d+1;++i)
    {
        tmp=1;
        for(int j=1;j<=n;++j)
            tmp=1ll*power(p[j][0],1ll*p[j][1]*i%(mod-1))*(1-power(p[j][0],d-i)+mod)%mod*tmp%mod;
        ans=(ans+1ll*A[i]*tmp%mod)%mod;
    }
    printf("%d\n",ans);
    return 0;
}