SPOJ 4191. Sky Code

http://www.spoj.com/problems/MSKYCODE/

容斥原理

在我不断的优化过后，终于过了

代码：

#include <iostream>
#include <cstdio>
#include <map>
#include <algorithm>
#include <cstring>
#include <string>
#include <vector>

#define ll long long

using namespace std;

const int N=10005;
const long long LINF=(long long)(1e18);
int a1[N],a2[N],r1[N],r2[N];
int p[N];
int b[N];
int num[N],k[N];
ll c[N][5];
bool prime[N];
int main()
{
    //freopen("data.in","r",stdin);
    for(int i=0;i<N;++i)
    for(int j=0;j<=i&&j<=4;++j)
    if(i==j||j==0) c[i][j]=1;
    else c[i][j]=c[i-1][j-1]+c[i-1][j];
    int n;
    while(scanf("%d",&n)!=EOF)
    {
        int M=0;
        memset(num,0,sizeof(num));
        for(int i=0;i<n;++i)
        {
            scanf("%d",&b[i]);
            ++num[b[i]];
            M=max(M,b[i]);
        }
        memset(k,0,sizeof(k));
        for(int i=2;i<=M;++i)
        for(int j=i;j<=M;j=j+i)
        k[i]+=num[j];
        int pn=0;
        memset(prime,true,sizeof(prime));
        for(int i=2;i<=M;++i)
        if(prime[i])
        {
            p[pn++]=i;
            for(int j=i+i;j<=M;j=j+i)
            prime[j]=false;
        }
        ll sum=c[n][4];
        for(int i=0;i<pn;++i)
        {
            a1[i]=p[i];
            r1[i]=i;
        }
        int len=pn;
        int f=1;
        while(len>0)
        {
            f=-f;
            int m=0;
            for(int i=0;i<len;++i)
            {
                if(k[a1[i]]>=4)
                {
                    r2[m]=r1[i];
                    a2[m++]=a1[i];
                    sum+=(f*c[k[a1[i]]][4]);
                }
            }
            len=0;
            for(int i=0;i<m;++i)
            {
                for(int j=r2[i]+1;j<pn;++j)
                {
                    if(a2[i]%p[j]&&a2[i]*p[j]<=M)
                    {
                        r1[len]=j;
                        a1[len++]=a2[i]*p[j];
                    }
                }
            }
        }
        cout<<sum<<endl;
    }
    return 0;
}