组合数学——容斥原理及其应用
容斥原理是计数中常用的一种方法。在讨论容斥原理的过程中,要用到以下集合论的基本性质。
德摩根(De Morgan)定理
若A和B是集合U的子集,则
例题
题意:给一个数n(n=10^8),求X1^4+X2^4+..Xk^4的和模1,000,000,007,其中1<=Xi<n,且Xi与n互素。
解:
四次方求和公式为:n*(n+1)*(2*n+1)*(3*n*n+3*n-1)/30
将n分解素因子为p1^e1*p2^e2*..*pk^ek,设集合Ai为pi的整数倍的数的集合。则答案为U-|A1UA2U...UAk|。
求|A1UA2U...UAk|即可用容斥原理求解。
注意先求30在模1,000,000,007下的逆元。
时间复杂度分析:
n为10^8,素因子分解最多有10个左右,容斥原理求解计算次数大约为2^10,复杂度已经很低了。
#include <iostream>
#include <cstdio>
#include <cstring>
#include <vector>
#include <cmath>
using namespace std;
const int Max = 10005;
const int MOD = 1000000007;
int prime[Max],num=0;
bool isPrime[Max]={0};
int x,y,rev;
vector<int> factors;
void getPrime(int Max)
{
int i,j;
int t;
for(i = 2; i < Max; i++)
{
if(!isPrime[i])
{
prime[num++] = i;
for(j = 2; (t=i*j) < Max; j++)
isPrime[t] = 1;
}
}
}
void extend_Euclid(int a, int b)
{
if(b==0)
{
x = 1;
y = 0;
return;
}
extend_Euclid(b, a%b);
int t = x;
x = y;
y = t - a/b*y;
}
void init()
{
extend_Euclid(30,MOD);
rev = (x%MOD+MOD)%MOD;
getPrime(Max);
}
void decompose(int n)
{
int t = (int)sqrt(n*1.0);
for(int i = 0; i<num && prime[i]<=t; i++)
{
if(n%prime[i]==0)
{
factors.push_back(prime[i]);
while(n%prime[i]==0)
n = n/prime[i];
}
}
if(n>1)
factors.push_back(n);
}
__int64 getSum(__int64 n)
{
__int64 result = n;
result = result*(n+1)%MOD;
result = result*(2*n+1)%MOD;
result = result*(((n*n%MOD*3+3*n%MOD-1)%MOD+MOD)%MOD)%MOD;
result = result*rev%MOD;
return result;
}
__int64 getPow(__int64 n)
{
return (n*n)%MOD*n%MOD*n%MOD;
}
__int64 in_out_principle(int start, __int64 n)
{
int i;
__int64 result = 0;
for(i = start; i < factors.size(); i++)
{
int tmpt = factors[i];
result = (result + getSum(n/tmpt)*getPow(tmpt)%MOD)%MOD;
result = ((result - in_out_principle(i+1, n/tmpt)*getPow(tmpt)%MOD) + MOD)%MOD;
}
return result;
}
int main()
{
int t;
init();
int n;
scanf("%d",&t);
while(t--)
{
scanf("%d",&n);
factors.clear();
decompose(n);
printf("%I64d\n",((getSum(n)-in_out_principle(0,n))%MOD+MOD)%MOD);
}
return 0;
}