poj1091 容斥原理的应用

　　这个题的意思是给你n+1个数， a1 a2 ... an an+1 其中 1<=ai<=M i!=n+1 an+1 = M, 求解存在x1..xn+1使x1*a1+x2*a2+x3*a3+..+xn+1*an+1 = 1的a序列的个数， 由数论部分知识我们可以知道上式要满足就得使（a1, a2, a3, .. an, an+1) = 1, 我们还能知道存在一些a序列他们的最大公约数不等于1， 然而我们要统计其中最大公约数为1的个数， 就可以使用容斥原理， 稍加思考我们可以发现(a1, a2 .. an+1) = 1的反面是（a1, a2, ... an+1）= pi的倍数， pi是M = pi^ai次方中的pi, 因此我们定义Ai为a序列最大公约数为pi的倍数的方案数根据容斥原理有 res = s - (A1+A2+...Ak) + (Ai并Aj i!=j) - ...... 代码如下：

#include <cstdio>
#include <algorithm>
#include <cstring>
#include <iostream>
#include <cmath>

using namespace std;
typedef long long LL;
const int maxn = 100000 + 10;
LL N, M;
LL pi[maxn], npi;

bool vis[maxn];
int prime[maxn], num;
void shai(int n)
{
memset(vis, 0, sizeof(vis));
int m = sqrt(n+0.5);
for(int i=2; i<=m; i++) if(!vis[i])
for(int j=i*i; j<=n; j+=i) vis[j] = 1;
num = 0;
for(int i=2; i<=n; i++) if(vis[i] == 0)
prime[num++] = i;
}

LL pow(LL A, LL B)  // A^B
{
LL res = 1;
while(B > 0)
{
if(B%2 == 1)
res = res * A;
A = A*A;
B /= 2;
}
return res;
}

int main()
{
shai(20010);
while( cin>>N>>M )
{
LL tp = M;
npi = 0;
for(int i=0; i<num; i++)
{
bool flog = false;
while(tp%prime[i] == 0)
{
tp /= prime[i];
pi[npi] = prime[i];
flog = true;
}
if(flog) npi++;
if(tp == 1) break;
}
if(tp!=1) pi[npi++] = tp;              //这里出完所有素数后有可能不为1
LL res = 0;
for(int i=1; i<(1<<npi); i++)
{
LL ge=0, d=1;
for(int j=0; j<npi; j++)
if(((i>>j)&1) == 1) ge++, d*=pi[j];
if(ge%2==1) res -= pow(M/d, N);
else res += pow(M/d, N);
}
cout<<res+pow(M, N)<<endl;
}
return 0;
}