Sum Of Gcd(欧拉函数+莫队算法详解)

Sum Of Gcd

 Given you a sequence of number a 1, a 2, ..., a n, which is a permutation of 1...n.
You need to answer some queries, each with the following format:
Give you two numbers L, R, you should calculate sum of gcd(a[i], a[j]) for every L <= i < j <= R. 

Input
First line contains a number T(T <= 10),denote the number of test cases.
Then follow T test cases.
For each test cases,the first line contains a number n(1<=n<= 20000).
The second line contains n number a 1,a 2,…,a n.
The third line contains a number Q(1<=Q<=20000) denoting the number of queries.
Then Q lines follows,each lines contains two integer L,R(1<=L<=R<=n),denote a query.
Output
For each case, first you should print “Case #x:”, where x indicates the case number between 1 and T.
Then for each query print the answer in one line.
Sample Input

1
5
3 2 5 4 1
3
1 5
2 4
3 3

Sample Output

Case #1:
11
4
0

题意：

给定一个$1到n的全排列A，若干个询问，每次询问给出一个区间[l,r]$，要求

$$\sum_{l\le i < j \le r} gcd(A_i,A_j)$$

分析：

对于这道题目，我们比较容易统计的是各个数对公约数的出现次数，而不容易统计最大公约数的出现次数

为什么呢？

举个例子假设我有数列$\{4,5,8\}$

4的约数有1,2,4；
5的约数有1,5；
8的约数有1,2,3,4,8；

我们设$cnt(d)$ 表示约数d在数列中出现的次数

根据上面枚举我们知道约数1出现了三次$cnt(1)=3$，约数2出现了两次$cnt(2) = 2$，约数5出现了1次$cnt(5)=1$，约数4出现了2次$cnt(4) = 2$，约数8出现了1次$cnt(8)=1$

那么既然是两个数的公约数d,对于每个约数，只要从约数个数中任取其中2个即为两个数的公约数，所以很明显通过组合数我们知道公约数为d可能的情况为

$$C_{cnt(d)}^2 = \frac{cnt(d)\times (cnt(d)-1)}{2}$$

所以和为公约数乘其出现次数即

$$C_{cnt(d)}^2 \times d$$

但是！

很明显这并不是我们想要求的，我们原本是想枚举d作为最大公约数，这样我们求得了最大公约的数对的个数然后乘d，最后合起来不就行了。但是我们发现其中会有大量重复，比如上面例子中对于4和8来说就很难办，因为它会把其公约数1,2,4全部统计一遍加进去，而我们只想要4的情况

因此这种想法以失败告终。。

但是其中思路和想法仍然对于我们进一步分析有重大意义

---

既然我们只能统计公因数，而不能统计最大公因数，那么我们就将计就计，我们就统计公约数

但是题目要求最大公约数的和

因此我们想要找到两个数最大公约数和公约数的关系：

1）首先我们要知道的第一个常识是两个数的公约数一定是这两个数的最大公约数的约数

2）欧拉函数求和公式：
设$d_1,d_2,\cdots d_r是n的因数$，则

$$\phi(d_1)+\phi(d_2)+\cdots +\phi(d_r) = n$$
这里不进行证明了

因此根据上面的性质个关系我们得到

$$gcd(A_i,A_j) = \sum_{d|gcd(A_i,A_j)}\phi(d)$$

（d是gcd的因子）

因此原式

$$\sum_{l\le i < j \le r} gcd(A_i,A_j)$$

$$=\sum_{l\le i < j \le r} \sum_{d|gcd(A_i,A_j)}\phi(d)$$

这样变成了求公约数d的欧拉函数$\phi(d)$之和

$\phi(d)$完全可以预处理出来，我们只要知道他有多少个就行了，一乘就是和

因此继续化简

$$原式 = \sum \phi(d)\sum_{l\le i < j \le r} [d|A_i]\cdot [d|A_j]$$

显然式子中$\sum_{l\le i < j \le r} [d|A_i]\cdot [d|A_j]$代表的就是d作为公约数时的所有可能情况

只要我们统计出数d是区间内多少个数的公约数即一开始我们设的$cnt(d)$那么

$$\sum_{l\le i < j \le r} [d|A_i]\cdot [d|A_j] =C_{cnt(d)}^2$$

---

现在问题的核心部分已经解决了我们已经将原式转化乘求某个数d作为公约数的可能数乘其欧拉函数，最后求和。每个数的约数和欧拉函数预处理出来。

因为题目给的是区间查询，而且是无修改的，因此想到可以使用莫队算法进行离线查询

对于假设区间$[L,R]$的答案我们已经知道了

1） 我们想要算$[L-1,R]或[L,R+1]$怎么办？

其实求这两个区间属于同一类都是扩展，增加了新的数n

对于新增加的数n，我们枚举其因子d

此时原来的$\sum_{d|n}\phi(d)\times C_{cnt(d)}^2变成\sum_{d|n}\phi(d)\times C_{cnt(d)+1}^2$

根据组合数的性质$C_n^k = C_{n-1}^k + C_{n-1}^{k-1}$

那么我们要求$\sum_{d|n}\phi(d)\times C_{cnt(d)+1}^2$只需要求出$C_{cnt(d)+1}^2$

而对于当前我们有$C_{cnt(d)}^2$只需要加上$C_{cnt(d)}^1 = cnt(d)这个cnt(d)是未更新前的cnt(d)$

所以对于答案而言我们只需加上$\phi(d) \times cnt(d)$即完成更新

然后在更新$cnt(d)$

2）如果想要求$[L+1,R]或[L,R-1]$我们发现他们都是需要除去一个数 n

对于要除去的n，这个数肯定是原来区间内的数（边界上的数）

同理枚举其约数d

仍然根据公式

$$C_{cnt(d)}^2 = C_{cnt(d)-1}^2 + C_{cnt(d)-1}^1$$

我们现在知道的是$C_{cnt(d)}^2$需要得到$C_{cnt(d)-1}^2$

所以只需减去$C_{cnt(d)-1}^1 = cnt(d)-1$即可

然后更新$cnt(d)$

当然这个题应该是可以用莫比乌斯反演做，但是看了网上的一篇代码写快300行就没看

code:

#include <stdio.h>
#include <string.h>
#include <iostream>
#include <algorithm>
#include <vector>
#include <queue>
#include <set>
#include <map>
#include <string>
#include <math.h>
#include <stdlib.h>
#include <time.h>
using namespace std;
typedef long long ll;
const int maxn = 20010;
int T,n,m,blocklen,block[maxn],a[maxn],fac[maxn][210];
ll phi[maxn],p[maxn],cnt[maxn],ans[maxn],sum;
bool isprime[maxn];
struct Query{
    int id,l,r;
    bool operator < (const Query &a)const{   /分块排序
        if(block[l] == block[a.l])
            return r < a.r;
        return block[l] < block[a.l];
    }
}q[maxn];

void Euler(){//预处理处欧拉函数
    phi[1] = 1;
    for(int i = 2; i < maxn; i++){
        phi[i] = i;
    }
    for(int i = 2; i < maxn; i++){
        if(phi[i] == i){
            for(int j = i; j < maxn; j += i){
                phi[j] = phi[j] / i * (i - 1);
            }
        }
    }
}

void init(){//记录下每个数的因数
    for(int i = 1; i < maxn; i++){
        fac[i][0] = 0;//记录因数个数
        for(int j = 1; j * j <= i; j++){
            if(i % j == 0){
                fac[i][++fac[i][0]] = j;
                if(j * j != i) fac[i][++fac[i][0]] = i / j;
            }
        }
    }
}

void expend(int x,ll add){
    for(int i = 1; i <= fac[a[x]][0]; i++){
        int j = fac[a[x]][i];//暂存a[x]的一个因子
        //更新
        if(add == 1) sum += cnt[j] * phi[j]; //这里add=1时乘的是cnt[j]
        else sum -= (cnt[j] - 1) * phi[j];//减的时候乘的是cnt[j]-1，原因看上面解释
        cnt[j] += add;//最后更新这个因子的个数（有多少个数含有这个因子）
    }
}

void Mo(){
    int l = 1,r = 0;
    memset(cnt,0,sizeof(cnt));
    sum = 0;
    for(int i = 1; i <= m; i++){
        while(q[i].l < l) expend(--l,1);
        while(q[i].r > r) expend(++r,1);//前两个是扩展所以要用前缀加减，这样处理加上新的数
        while(q[i].l > l) expend(l++,-1);
        while(q[i].r < r) expend(r--,-1);//后两个是减少所以要把原来的边界减去，所以处理的是原来的，传进去的原来边界值
        ans[q[i].id] = sum;
    }
}
int main(){
    Euler();
    init();
    scanf("%d",&T);
    int cas = 0;
    while(T--){
        scanf("%d",&n);
        for(int i = 1; i <= n; i++) scanf("%d",&a[i]);
        scanf("%d",&m);
        blocklen = (int)sqrt(n);
        for(int i = 1; i <= n; i++){
            block[i] = i / blocklen;
        }
        for(int i = 1; i <= m; i++){
            scanf("%d%d",&q[i].l,&q[i].r);
            q[i].id = i;
        }
        sort(q+1,q+1+m);
        Mo();
        printf("Case #%d:\n",++cas);
        for(int i = 1; i <= m; i++)
            printf("%lld\n",ans[i]);
    }
    return 0;
}