POJ 3904 Sky Code (互素的四元组+容斥原理)

POJ 3904 Sky Code (互素的四元组+容斥原理) ：http://poj.org/problem?id=3904

题目大意：

给定n个数据，求其中互素的四元组（两两可以不互素）的个数。

题目分析：

要想求互素的四元组的个数，我们可以先去求不互素的四元组的个数ans，然后所有可能的取法数C(n,4)（n个里面取4个的组合数公式）减去不互素的四元组的个数。

考虑到不互素的四元组必满足四个数必有除了1以外的大于等于一个的因子（可能是素数，也可能不是素数），则可以对每一个数都进行因子分解，由于这些因子可能会被重复的计算，所以应用容斥原理，对于每一个因子，该因子若含有除了1以外的奇数个因子的话，则加上，否则，含有偶数个因子的话，则减去；举例说明：n个数中，以2为因子的数有a个，以3为因子的数有b个，以6为因子的数有c个，则满足不互素的四元组的个数有C(a,4)+C(b,4)-C(c,4)个。

关于代码实现中的二进制枚举部分，有几个关键的地方：（转）

（1）： 理解cnt[i]表示含有因子i的数据的个数（2,4,6,8,10这组数据中cnt[i]=5） num[i]表示因子i含有的不同素因子的个数（num[6]=2,减）

（2）： 怎样得到cnt[i]、num[i]，或者说怎样实现素因子组合？用到了二进制表示的思想，第j位的1表示第j个素因子参与累乘，0表示不参与累乘，比如含有素因子2*3*5*7，1010表示2*1*5*1，0001表示1*1*1*7，以此类推……这样，只要素因子分解完记录含有素因子的个数tol，那么就有tol位二进制，注意到2*3*5*7*11*13>10000，所以最多有6位二进制；对于任意一个范围内的数i，遍历所有位j（从低位到高位），i&(1<<j)==1，表示i的第j位二进制为1，意思是第j个素因子参与累乘；依据这样的意义，遍历二进制范围内的所有数i和i对应的所有二进制数位j，更新数组cnt[k]和num[k]即可！

代码实现：

#include <iostream>
#include <cstdio>
#include <cstring>

using namespace std;

const int MAXN=10005;

long long  prime[MAXN];  //存放素因子
long long C[MAXN];
long long count[MAXN];
long long num[MAXN];

void zuheshu() //计算组合数，i个里面选4个
{
    memset(C,0,sizeof(C));
    for(long long i=4; i<=10000; i++)
    {
        C[i]=i*(i-1)*(i-2)*(i-3)/24;
    }
}

void divide(int x)
{
    long long tol=0;
    for(long long i=2;i*i<=x;i++)
    {
        if(x%i==0)
        {
            prime[tol++]=i;
        }
        while(x%i==0)
        {
            x/=i;
        }
    }
    if(x!=1)
    {
        prime[tol++]=x;
    } 
    //对x进行素因子分解，tol表示素因子的个数
    for(long long i=1;i<(1<<tol);i++) //利用二进制的每一位0,1对素因子进行组合
    {
        long long k=1;
        long long sum=0;
        for(long long j=0;j<tol;j++) //从低位到高位遍历，找i在哪一位是1
        {
            if(i&(1<<j)) //与运算，j位为1，执行if语句
            {
                k*=prime[j];//累乘素因子
                sum++;
            }
        }
        count[k]++; //记录当前因子的个数
        num[k]=sum; //记录当前因子是由多少个素因子组成
    }
}

int main()
{
    long long  n,m;
    zuheshu();
    while(scanf("%lld",&n)!=EOF)
    {
        memset(count,0,sizeof(count));
        memset(num,0,sizeof(num));
        for(int i=0; i<n; i++)
        {
            scanf("%lld",&m);
            divide(m); //素因子分解，并统计
        }
        long long ans=0;
        for(long long i=0;i<=10000;i++)
        {
            if(count[i])
            {
                if(num[i]%2)
                {
                    ans+=C[count[i]];
                }
                else
                ans-=C[count[i]];
            }
        }
        printf("%lld\n",C[n]-ans);
    }
    return 0;
}