SP4191 天空代码 分析

题目概述

问有多少个 \((a,b,c,d)\)，在 \(n\) 个数的 \(x\) 满足 \(\gcd\{x_a,x_b,x_c,x_d\}=1\).

其中，\(n,\max x\leq 10^4\)。

分析

套路经典题目，记录一下。

设 \(f(d)\) 表示选 \(4\) 个数，其最大公约数为 \(d\) 的倍数的个数。

设 \(F(d)\) 表示选 \(4\) 个数，其最大公约数为 \(d\) 的个数。

那么我们 \(f(d)\) 是好求的：

\[f(d)=C_{cnt_d}^{4} \]

我们考虑怎么通过 \(f\) 求 \(F\)，只需要考虑容斥即可。

因为我要 \(1\times d\) 的 \(f\)，再减去 \(2\times d\) 的 \(f\)，再减去 \(3\times d\) 的 \(f\)，再加上 \(6\times d\) 的 \(f\) 等等即可。

我们发现这其实就是其分解质因数的个数，容易想到莫比乌斯函数——一个天然的容斥系数。

于是：

\[F(d)=\sum_{d\mid k}\mu(\frac{k}{d})f(k) \]

所以说我们求的 \(F(1)=\sum_{k=1}^{mx}mu(k)f(k)\)。

然后直接做就行了哦。

代码

时间复杂度 \(\mathcal{O}(\sum n\sqrt n)\)。

#include <iostream>
#include <cstdio>
#include <cstring>
#include <stdlib.h>
#include <algorithm>
#include <vector>
#include <cmath>
#define int long long
#define N 10005
using namespace std;
int gcd(int a,int b) {
    return b ? gcd(b,a % b) : a;
}
int calc(int x) {
    if (x <= 3) return 0;
    return x * (x - 1) * (x - 2) * (x - 3) / 24;
}
int prime[N];
bool vis[N];
int mu[N];
void init(int n) {
    mu[1] = 1;
    int cnt = 0;
    for (int i = 2;i <= n;i ++) {
        if (!vis[i]) prime[++cnt] = i,mu[i] = -1;
        for (int j = 1;j <= cnt && prime[j] * i <= n;j ++) {
            vis[i * prime[j]] = 1;
            if (i % prime[j] == 0) break;
            mu[i * prime[j]] = -mu[i];
        }
    }
}
int a[N],cnt[N],f[N];
signed main(){
    init(1e4);
    int n;
    int mx = 0;
    while(~scanf("%lld",&n)) {
        memset(f,0,sizeof f),memset(cnt,0,sizeof cnt);
        // for (int i = 0;i <= mx;i ++) f[i] = cnt[i] = 0;
        mx = 0;
        for (int i = 1;i <= n;i ++) scanf("%lld",&a[i]),mx = max(mx,a[i]);
        for (int i = 1;i <= n;i ++) {
            int x = a[i];
            int t = sqrt(x);
            for (int j = 1;j <= t;j ++)
                if (x % j == 0) {
                    cnt[j] ++;
                    if (j != x / j) cnt[x / j] ++;
                }
        }
        for (int i = 1;i <= mx;i ++) f[i] = calc(cnt[i]);
        int ans = 0;
        for (int i = 1;i <= mx;i ++) ans += mu[i] * f[i];
        printf("%lld\n",ans);
        // int ans = 0;
        // for (int i = 1;i <= n;i ++)
        //     for (int j = i + 1;j <= n;j ++)
        //             for (int k = j + 1;k <= n;k ++)
        //                     for (int l = k + 1;l <= n;l ++)
        //                         if (gcd(a[i],gcd(a[j],gcd(a[k],a[l]))) == 1)
        //                             ans ++;
        // printf("%lld\n",ans);
    }
    return 0;
}