HDU 5072：同色三角形，容斥原理，预处理的顺序

谨以此题纪念逝去的银牌……

题目大意：给定N个正整数，互不相同，问这里面可以挑出多少这样的三元组：要么这三个数两两互质，要么这三个数两两不互质。N<=10w, 每个数也<=10w。

我们可以先试着考虑这个问题的反面：挑出这样的三元组，要么2对互质1对不是，要么2对不是，一对互质。这样的三元组个数是有办法计算的：对于每个数，随便挑一个和他互质的数，再随便挑一个和他不互质的数，这样的三元组就会满足要求，因为不管那两个数之间互不互质都没关系了。如果挑出所有这样的三元组，你会发现，每组都被计算了两次：因为每组当中都有两个数字满足剩下的两个数字和他是一互质一不互质，所以在枚举每个数的时候，这两个数都被挑过。所以，我们处理出与每个数互质以及不互质的数的个数，将它们乘起来求和，最后除2就是答案。

这其实脱胎于一个经典的同色三角形模型：平面上有N个点，每两个点之间有一条红边或者蓝边，问这里面有多少个同色的三角形？做法和上面一样，先枚举每个点，找到引出红边的数量，和引出蓝边的数量，乘起来求和除2就是答案。

那么，剩下的问题是，对于每个数，如果处理出这里面有多少个数和他互质呢？我们首先将它分解质因数，然后对于每个质因子，都去扫一遍，将能整除的都排除掉即可。10w的数，不同质因子最多只有6个（2x3x5x7x11x13 = 20000多）。但是这样有一个问题，就是有些数会被抛去两次，比如6，在枚举2时会被减一次，枚举3时又会被减一次，这就不对了。所以我们要使用容斥原理，不光枚举每个质因子，还要枚举他们的乘积，如果使用了偶数个质因子，就加回来。比如6使用了2、3两个质因子，就要加回来。这样就能正确了。

这样的复杂度是2^6 * n *扫描复杂度，还是差了一点。关键就在于扫描，我们枚举出一个质因子的积以后如何才能很快地知道这里有多少数能被它整除呢？如果我们枚举出了这个质因子积，再去把10w个数扫描一遍的话，即使有记忆化，复杂度也太高了，因为其实根本没有必要扫描，我们完全可以知道每个质因子积能整除多少个数。方法就是在枚举质因子积之前，以同样的方法枚举一遍质因子积，然后枚举出几，就在几那个地方+1，记录一下，然后第二次枚举的时候，就直接可以知道能被这个数整除的有几个数了。

比如12，它的质因子积是1、2、3、6，那么在第一次枚举的时候，记录质因子积的数组div_num[1, 2, 3, 6]就都应该+1， 因为1、2、3、6都能整除12. 依照此法就能预处理出每个质因子积都能整除多少数，第二次枚举的时候，我们就可以用O（1）复杂度直接使用了。

可以发一下代码，DFS1是预处理的，DFS2是计算的

#include<iostream>
#include<stdio.h>
#include<math.h>
#include<string.h>
#include<algorithm>
using namespace std;

typedef long long LL;

int t;
LL n;
int a[100100];

const int MAXN = 100010;
int prime[MAXN+1];

int prime_factor[100100][10];
int p_pf[100100];

LL now_coprime=0;
LL coprime[100100];

LL div_num[100100];

void init()
{
memset(a,0,sizeof(a));
memset(prime_factor, 0, sizeof(prime_factor));
memset(p_pf, 0, sizeof(p_pf));

memset(coprime, 0, sizeof(coprime));
memset(div_num, -1, sizeof(div_num));
return;
}

void get_prime()
{
memset(prime, 0, sizeof(prime));
for(int i=2; i<=MAXN; i++)
{
if(!prime[i])
prime[++prime[0]] = i;

for(int j=1; j<=prime[0] && prime[j] <=MAXN/i; j++)
{
prime[prime[j] *i] = 1;
if(i%prime[j]==0)
break;
}
}
}

void divide(int now_p)
{
int now = a[now_p];
int i;
for(i=1; prime[i]<=now / prime[i]; i++)
{
if(now%prime[i] == 0)
{
prime_factor[now_p][++p_pf[now_p]] = prime[i];
while(now%prime[i] == 0)
now /= prime[i];
}
}
if(now>1)
prime_factor[now_p][++p_pf[now_p]] = now;
return;
}

void DFS1(int now_p, int now_pos, int now_prod)
{
if(now_pos == p_pf[now_p]+1)
{
div_num[now_prod]++;
return;
}

DFS1(now_p, now_pos+1, now_prod * prime_factor[now_p][now_pos]);
DFS1(now_p, now_pos+1, now_prod);

return;
}

void DFS2(int now_p, int now_pos, int now_prod, bool odd)
{
if(now_pos == p_pf[now_p]+1)
{
if(odd)
now_coprime -= div_num[now_prod];
else
now_coprime += div_num[now_prod];
return;
}

DFS2(now_p, now_pos+1, now_prod * prime_factor[now_p][now_pos], !odd);
DFS2(now_p, now_pos+1, now_prod, odd);

return;
}

int main()
{
get_prime();
cout<<prime[0]<<endl;

scanf("%d", &t);
int files;
for(files=1; files<=t; files++)
{
init();
scanf("%I64d", &n);
int i;
for(i=1; i<=n; i++)
{
scanf("%d", &a[i]);
divide(i);
}
for(i=1; i<=n; i++)
{
DFS1(i, 1, 1);
}
for(i=1; i<=n; i++)
{
now_coprime=0;
DFS2(i,1,1, false);
coprime[i] = now_coprime;
}
LL ans=0;
for(i=1; i<=n; i++)
{
ans += coprime[i] * (n-1 - coprime[i]);
}
ans/=2;

LL all = n*(n-1)*(n-2) / 6;

printf("%I64d\n", all-ans);
}
system("pause");
return 0;
}