bzoj 3601: 一个人的数论 莫比乌斯反演+高斯消元

题目大意：
给定$n=p_1^{{\alpha }_1}\ast p_2^{{\alpha }_1}\ast ...\ast p_w^{{\alpha }_w}$和$d$，求
$$\sum _{i=1}^n{i}^d[gcd(i,n)==1]$$
其中，$d\le 100$，$w\le 1000$，$p_i$为质数，$p_i,q_i\le 10^9$。

分析：
根据题意推式子得到
$$ans=\sum _{p|n}\mu (p)\ast p^d\ast \sum _{i=1}^{n/p}{i}^d$$
对于最后面的自然数幂和，可以通过高斯消元求出多项式，设多项式为$a_{d+1}\ast x^{d+1}+a_d\ast x^d+...+a_0$.
显然$a_0=0$，这一位就不用消了（避免出错）。
原式变成
$$ans=\sum _{p|n}\mu (p)\ast p^d\ast \sum _{i=1}^{d+1}{a}_i\ast (n/p{)}^i$$
把$i$移到前面，得到
$$ans=\sum _{i=1}^{d+1}{a}_i\sum _{p|n}\mu (p)\ast p^d\ast (n/p{)}^i$$
我们可以枚举$i$，当$i$一定时，后面这些是一个积性函数，我们可以分别计算出每一个$p_i^{{\alpha }_i}$的贡献，乘起来即可。而根据$\mu$的性质，只有$p=1$和$p=p_i$时有贡献，这个可以$O(1)$算出。总复杂度是$O(wd)$的，另外高斯消元是$O(n^3)$。

代码：

/**************************************************************
    Problem: 3601
    User: ypxrain
    Language: C++
    Result: Accepted
    Time:812 ms
    Memory:9240 kb
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long
 
const int maxn=1007;
const LL mod=1e9+7;
 
using namespace std;
 
int n,k;
int p[maxn],q[maxn];
LL a[maxn][maxn],b[maxn],r[maxn];
LL sum,ans,num;
 
LL power(LL x,LL y)
{
    if (y==1) return x;
    if (y==0) return 1;
    LL c=power(x,y/2);
    c=(c*c)%mod;
    if (y&1) c=(c*x)%mod;
    return c;
}
 
void guass(int n)
{
    for (int i=1;i<=n;i++)
    {
        for (int j=i+1;j<=n;j++)
        {
            if (a[j][i]>a[i][i])
            {
                for (int k=i;k<=n;k++) swap(a[i][k],a[j][k]);
                swap(b[i],b[j]);
            }
        }
        LL inv=power(a[i][i],mod-2);
        for (int j=1;j<=n;j++)
        {
            if (i==j) continue;
            LL rate=a[j][i]*inv%mod;
            for (int k=i;k<=n;k++) a[j][k]=(a[j][k]+mod-a[i][k]*rate%mod)%mod;
            b[j]=(b[j]+mod-b[i]*rate%mod)%mod;
        }
    }
    for (int i=1;i<=n;i++) r[i]=b[i]*power(a[i][i],mod-2)%mod;
}
 
int main()
{
    scanf("%d%d",&k,&n);
    for (int i=1;i<=n;i++)
    {
        scanf("%d%d",&p[i],&q[i]);
        q[i]=power(p[i],q[i]);
    }   
    for (int i=1;i<=k+1;i++) a[i][1]=i;
    for (int j=2;j<=k+1;j++)
    {
        for (int i=1;i<=k+1;i++) a[i][j]=(a[i][j-1]*(LL)i)%mod;
    }
    for (int i=1;i<=k+1;i++)
    {
        sum=(sum+power(i,k))%mod;
        b[i]=sum;
    }
    guass(k+1);
    for (int i=1;i<=k+1;i++)
    {
        num=1;
        for (int j=1;j<=n;j++)
        {
            LL e=q[j]*power(p[j],mod-2)%mod;
            num=num*((power(q[j],i)+mod-power(p[j],k)*power(e,i)%mod)%mod)%mod;
        }
        ans=(ans+r[i]*num%mod)%mod;
    }
    printf("%lld\n",ans);
}