【题解】Sum Of Gcd 莫队+数学 hdu4676 ——quhengyi11

传送门

这道题比较神仙所以我来写篇 blog 记录下 qwq

首先它让我们求的是

$$\sum_{L\leq i<j\leq R}gcd(a[i],a[j])$$

我刚开始是想到 $gcd$是具有差分不变的性质，但是这似乎只能转化区间呀所以没啥用

后来想记录因数来统计，可是感觉要容斥搞半天就自闭了

好了现在我们来讲正解

首先我们可以化一化式子

$$\sum_{L\leq i<j\leq R}\sum_{d|gcd(a[i],a[j])}\phi(d)$$

这样我们就可以通过统计因子个数来计算答案了呢

记因子 $d$出现的次数为 $cnt[d]$，显然上式等价于

$$\sum_{d=1}^n \sum_{L\leq i<j\leq R} \phi(d)C_{cnt[d]}^2$$

经过打表我们可以知道 $2e4$以内平均因数个数约为 $10$个，且最多的有 $80$个因子

于是发现这挺棒的呀我们只需要能快速维护区间每个因子的个数就行了

事实上用线段树啥的维护的话，有个问题就是当你询问的时候你不能确定哪些因子在这个区间内，所以你就不可避免地要枚举因子，然后就会 $TLE$

我们考虑莫队

增加一个数，要将它的所有因子 $d$的 $cnt$全部 $+1$，由 $pascal$恒等式，一个因子对答案贡献了 $C_{cnt[d]}^1=cnt[d]$，枚举一遍因子就行了，减少一个数同理。时间就是一个 $10$的常数。

代码

#include<bits/stdc++.h>
#define fo(i, a, b) for (Re int i = (a); i <= (b); ++i)
#define fd(i, a, b) for (Re int i = (a); i >= (b); --i)
#define edge(i, u) for (Re int i = head[u], v = e[i].v; i; i = e[i].nxt, v = e[i].v)
#define Re register
#define pb push_back
#define F first
#define S second
#define ll unsigned long long
#define lowbit(x) (x & -x)
#define mod 1000000007
#define N 20005
int T, n, a[N], m, block[N], cnt[N], l, r;
ll anss[N], ans;
std::vector<int> d[N];
struct ques{
    int x, y, id;
    friend bool operator < (ques x, ques y)
    {
        return block[x.x] == block[y.x] ? x.y < y.y: block[x.x] < block[y.x];
    }
}q[N];
std::bitset<N> vis;
int p[N], tot, phi[N];
inline void init ()
{
    int n = N - 3;
    phi[1] = 1;
    fo (i, 2, n)
    {
        if (!vis[i])
        {
            p[++tot] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; j <= tot && i * p[j] <= n; ++j)
        {
            vis[i * p[j]] = 1;
            if (!(i % p[j]))
            {
                phi[i * p[j]] = phi[i] * p[j];
                break;
            }
            phi[i * p[j]] = phi[i] * phi[p[j]];
        }
    }
}
inline void add (int x, bool opt)
{
    if (opt)
    {
        int sz = d[x].size();
        for (int i = 0; i < sz; ++i)
        {
            int now = d[x][i];
            ans += cnt[now] * phi[now];
            ++cnt[now];
        }
    }
    else
    {
        int sz = d[x].size();
        for (int i = 0; i < sz; ++i)
        {
            int now = d[x][i];
            --cnt[now];
            ans -= cnt[now] * phi[now];
        }
    }
}
int main ()
{
    init();
    int up = 20000;
    fo (i, 1, up)
    {
        int sq = (int) sqrt(i - 0.5);
        fo (j, 1, sq)
            if (!(i % j))
            {
                d[i].pb(j);
                d[i].pb(i / j);
            }
        ++sq;
        if (sq * sq == i) d[i].pb(sq);
    }
    scanf("%d", &T);
    fo (qwq, 1, T)
    {
        memset(cnt, 0, sizeof cnt);
        scanf("%d", &n);
        fo (i, 1, n) scanf("%d", &a[i]);
        int sq = (int) sqrt(n);
        fo (i, 1, n) block[i] = i / sq;
        scanf("%d", &m);
        fo (i, 1, m)
        {
            scanf("%d %d", &q[i].x, &q[i].y); 
            q[i].id = i;
        }
        std::sort(q + 1, q + m + 1);
        l = 1; r = 0; ans = 0;
        fo (i, 1, m)
        {
            while (l > q[i].x) add(a[--l], 1);
            while (r < q[i].y) add(a[++r], 1);
            while (l < q[i].x) add(a[l++], 0);
            while (r > q[i].y) add(a[r--], 0);
            anss[q[i].id] = ans;
        }
        printf("Case #%d:\n", qwq);
        fo (i, 1, m)
            printf("%lld\n", anss[i]);
    }
    return 0;
}

总结：

1. 在求最大公约数 $n$对答案的贡献的时候，我们可以将它转化为 $\sum_{d|n}\phi(d)$，这样就把贡献转移到了因子的上面。

2. 每个数的因子个数似乎也是很小的呢，$2\times 10^6$内的数平均下来才 $14$个因子。然而它的个数似乎还是找不到公式呢，反正最大的应该是小于 $sqrtn$的 (臆断)，可以现场打表检验呢反正也不麻烦。