hdu 4059 容斥原理 通项公式

其实这道题目的难点不在容斥原理,在于大数越界、取余等等的细节

首先求出通项公式为

n(2n+1)(n+1)(3n² +3n-1)/30
求(a/b)%mod有两种方法
1、原式=a%(b*n)/b
2、原式=a*b^(phi(mod)-1)%mod;(其中b与mod互质)
本题中,由于mod过大,采用第一种方法会出错,于是用了第二种方法。
容斥过程中比如要加上2的四次,4的四次,6的四次。。。。可把2的四次提取出来,另外一项就变成了四次方和了
View Code
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
using namespace std;
const __int64 mod = 1000000007;
__int64 n;
__int64 power(__int64 a,__int64 b,__int64 c)
{
__int64 ret=1;
while(b>0)
{
if(b&1)
ret=ret*a%c;
a=a*a%c;
b>>=1;
}
return ret;
}
__int64 calc(__int64 n)
{
__int64 a=n,b=(2*n+1),c=(n+1),d=(3*n*n+3*n-1)%mod;
__int64 sum=1;
sum*=a;sum%=mod;
sum*=b;sum%=mod;
sum*=c;sum%=mod;
sum*=d;sum%=mod;
__int64 t=power(30,mod-2,mod);
return sum*t%mod;
}
__int64 sici(__int64 n)
{
__int64 sum=1;
for(int i=1;i<=4;i++)
{
sum*=n;
if(sum>=mod) sum%=mod;
}
return sum;
}
__int64 solve(__int64 r,__int64 n)
{
vector<__int64> p;
__int64 i;
for(i=2;i*i<=r;i++)
{
if(r%i==0)
{
p.push_back(i);
while(r%i==0) r/=i;
}
}
if(r>1) p.push_back(r);
__int64 sum=0;
for(__int64 num=1;num<(1<<p.size());num++)
{
__int64 mult=1,ones=0;
for(i=0;i<p.size();i++)
{
if(num&(1<<i))
{
ones++;
mult*=p[i];
if(mult>=mod) mult%=mod;
}
}
if(ones%2) sum+=sici(mult)*calc(n/mult);
else sum-=sici(mult)*calc(n/mult);
sum%=mod;
}
return sum;
}
int main()
{
__int64 t,i,j;
scanf("%I64d",&t);
while(t--)
{
scanf("%I64d",&n);
printf("%I64d\n",(calc(n)-solve(n,n)%mod+mod)%mod);
}
return 0;
}



posted @ 2012-02-26 19:21  Because Of You  Views(827)  Comments(0Edit  收藏  举报