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

来源:互联网 发布:淘宝上怎么买汽枪 编辑:程序博客网 时间:2024/06/05 14:49

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

看这位神犇的吧

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;}


0 0
原创粉丝点击