[莫比乌斯反演 高斯消元 数学技巧] BZOJ 3601 一个人的数论

推导过程直接拉过来不是很好 

看这位神犇的吧

http://www.cnblogs.com/jianglangcaijin/p/4033399.html

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

const ll P=1e9+7;

inline ll Pow(ll a,ll b){
  ll ret=1;
  for (;b;b>>=1,a=a*a%P)
    if (b&1)
      ret=ret*a%P;
  return ret;
}

inline ll Inv(ll a){
  return Pow(a,P-2);
}

const int N=105;

int w,d;
ll ip[1005],ia[1005];
ll a[N][N],A[N];

int main(){
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  scanf("%d%d",&d,&w);
  for (int i=1;i<=w;i++) scanf("%lld%lld",ip+i,ia+i);
  ll isum=0;
  for (int i=1;i<=d+2;i++){
    ll tem=1;
    for (int j=1;j<=d+2;j++) a[i][j]=tem,tem=tem*i%P;
    a[i][d+3]=a[i-1][d+3]+Pow(i,d);
  }
  int m=d+2;
  for (int i=1;i<=m;i++){
    int k=i;
    while (!a[k][i]) k++;
    for (int j=1;j<=m+1;j++) swap(a[i][j],a[k][j]);
    ll inv=Inv(a[i][i]);
    for (int j=1;j<=m;j++){
      if (j==i) continue;
      ll t=a[j][i]*inv%P;
      for (int k=1;k<=m+1;k++)
	(a[j][k]+=P-a[i][k]*t%P)%=P;
    }
  }
  for (int i=0;i<=d+1;i++)
    A[i]=a[i+1][m+1]*Inv(a[i+1][i+1])%P;
  ll Ans=0;
  for (int i=0;i<=d+1;i++){
    ll hi=1;
    for (int j=1;j<=w;j++)
      (hi*=(Pow(ip[j],ia[j]*i%(P-1))+P-Pow(ip[j],(d+(ia[j]-1)*i)%(P-1)))%P)%=P;
    (Ans+=A[i]*hi%P)%=P;
  }
  printf("%lld\n",Ans);
  return 0;
}