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

标签：bzoj   bzoj3601   莫比乌斯反演   高斯消元   

题目大意：求Σ[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 1000000007
using 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;
}

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

标签：bzoj   bzoj3601   莫比乌斯反演   高斯消元   

原文地址：http://blog.csdn.net/popoqqq/article/details/42738567