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

题目大意：求Σ[i|n]i^d

围观题解：http://www.cnblogs.com/jianglangcaijin/p/4033399.html

果然我还是太蒻了- -

此外Σ[1<=i<=n]i^m的零次项注定为0- - 所以常数项不用消了- -

#include <cstdio>#include <cstring>#include <iostream>#include <algorithm>#define M 110#define MOD 1000000007using namespace std;long long d,w;long long f[M][M],a[M];long long h[M],ans;long long Quick_Power(long long x,long long y){long long re=1;(y+=MOD-1)%=(MOD-1);while(y){if(y&1) (re*=x)%=MOD;(x*=x)%=MOD;y>>=1;}return re;}void Gauss_Elimination(){int i,j,k;for(i=1;i<=d+1;i++){long long temp=i;for(j=1;j<=d+1;j++,(temp*=i)%=MOD)f[i][j]=temp;for(j=1,temp=1;j<=d;j++,(temp*=i)%=MOD);(f[i][d+2]=(i?f[i-1][d+2]:0)+temp)%=MOD;}for(i=1;i<=d+1;i++){for(j=i;j<=d+1;j++)if(f[j][i])break;for(k=1;k<=d+2;k++)swap(f[i][k],f[j][k]);for(j=i+1;j<=d+1;j++){long long temp=MOD-f[j][i]*Quick_Power(f[i][i],-1)%MOD;for(k=i;k<=d+2;k++)(f[j][k]+=f[i][k]*temp)%=MOD;}}for(i=d+2;i;i--){for(j=i+1;j<=d+1;j++)(f[i][d+2]-=f[i][j]*a[j])%=MOD;a[i]=(f[i][d+2]*Quick_Power(f[i][i],-1)%MOD+MOD)%MOD;}//for(i=1;i<=d+1;i++)//cout<<a[i]<<endl;}int main(){int i,j;cin>>d>>w;Gauss_Elimination();for(j=1;j<=d+1;j++)h[j]=1;for(i=1;i<=w;i++){int p,a;scanf("%d%d",&p,&a);for(j=1;j<=d+1;j++)( h[j]*=Quick_Power(p,(long long)a*j)%MOD * (1-Quick_Power(p,d-j)%MOD) % MOD )%=MOD;}for(j=1;j<=d+1;j++)(ans+=a[j]*h[j])%=MOD;cout<<(ans%MOD+MOD)%MOD<<endl;return 0;}