解题思路

核心思想

利用容斥原理计算不满足条件的四元组数目（即 GCD 大于 1 的四元组），再用总四元组数目减去该值，得到符合条件的数目。

步骤解析

1. 统计每个数的出现次数：使用数组 cnt[d] 记录 ID 能被 d 整除的数的个数。
2. 容斥原理求 GCD 为 d 的数的个数：
   • 从大到小遍历所有可能的 d（1 到 10000）。
   • 对于每个 d，若 cnt[d] 表示能被 d 整除的数的个数，则 至少能被 d 整除的数的个数为 f[d] = sum_{k是d的倍数} cnt[k]。
   • 根据容斥，恰好 GCD 为 d 的数的个数为 g[d] = f[d] - sum_{k是d的倍数且k>d} g[k]。
3. 计算四元组数目：
   • 对于每个 d > 1，若 g[d] >= 4，则贡献的四元组数目为 C(g[d], 4)（组合数）。
   • 总不满足条件的数目为所有 d > 1 的贡献之和。
   • 最终结果 = 总四元组数目（若 N ≥ 4，为 C(N, 4)，否则为 0） - 不满足条件的数目。

关键优化

• 预处理倍数关系：倒序遍历 d（从 10000 到 1），利用倍数关系快速计算 f[d] 和 g[d]，避免重复计算。
• 组合数计算：提前预处理组合数或使用公式 C(n, 4) = n*(n-1)*(n-2)*(n-3)/24，注意处理 n < 4 时结果为 0。

代码实现

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<string>
#include<cstring>
#include<cmath>
#include<ctime>
#include<algorithm>
#include<utility>
#include<stack>
#include<queue>
#include<vector>
#include<set>
#include<map>
#define EPS 1e-9
#define PI acos(-1.0)
#define INF 0x3f3f3f3f
#define LL long long
const int MOD = 1E9+7;
const int N = 10000+5;
const int dx[] = {0,0,-1,1,-1,-1,1,1};
const int dy[] = {-1,1,0,0,-1,1,-1,1};
using namespace std;
 
LL n;
LL mu[N],prime[N];
bool bprime[N];
LL a[N],num[N];
void getMu(LL n){//线性筛求莫比乌斯函数
 
    mu[1]=1;//根据定义，μ(1)=1
 
    LL cnt=0;
    memset(bprime,false,sizeof(bprime));
    for(LL i=2;i<=n;i++){//求2~n的莫比乌斯函数
        if(!bprime[i]){
            prime[cnt++]=i;//存储质数
            mu[i]=-1;//i为质数时，μ(1)=-1
        }
        for(LL j=0;j<cnt;j++){//枚举i之前的素数个数
            LL k=i*prime[j];//素数的乘积
            if(k>n)//剪枝
                break;
            bprime[k]=true;//不是质数
            if(i%prime[j])//i不是primes[j]的整数倍时，i*prime[j]就不会包含相同质因子
                mu[k]=-mu[i];//mu[k]=mu[i]*mu[prime[j]]，因为prime[j]是质数，mu值为-1
            else{
                mu[k]=0;
                break;//留到后面再筛
            }
        }
    }
}
LL Cn4(LL m){
    if(m==0)
        return 0;
    return m*(m-1)*(m-2)*(m-3)/24;
}
void getM(){
    memset(num,0,sizeof(num));
    for(LL i=0;i<n;i++){
        LL x=a[i];
        LL m=sqrt(x);
        for(LL j=1;j<=m;j++){
            if(x%j==0){
                num[j]++;
                num[x/j]++;
            }
        }
        if(m*m==x)
            num[m]--;
    }
}
int main(){
    getMu(N);
    while(scanf("%lld",&n)!=EOF&&n){
        for(LL i=0;i<n;i++)
            scanf("%lld",&a[i]);
        getM();
 
        LL res=0;
        for(LL i=1;i<N;i++)
            res+=mu[i]*Cn4(num[i]);
        printf("%lld\n",res);
    }
    return 0;
}