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