HDU 4676 Sum Of Gcd（莫队+莫比乌斯反演）

来看答案怎么求？

$$\sum_{i=L}^R\sum_{j=i+1}^Rgcd(a[i],a[j]) =\sum_{i=L}^R\sum_{j=i+1}^R\sum_{d|gcd(a[i],a[j])}\phi(d) =\sum_d\phi(d)\sum_{i=L,d|a[i]}^R\sum_{j=i+1,d|a[j]}^R1$$
后一项就是 $[L,R]$这段区间内有多少个 $d$的倍数。
$cnt[d][x]$代表前x项中d的倍数的个数。
原式可化为
$$\sum_{i=L}^R\sum_{j=i+1}^Rgcd(a[i],a[j]) =\sum_d\phi(d)*(cnt[d][R]-cnt[d][L-1])$$
如果每一次都是直接求，显然不能满足复杂度，也会爆内存，因为只是一道无修改的查询，自然想到某该的莫队！！！
需要维护的就是该区间内 $cnd[d]$的值。
不妨假设 $[L,R]$到 $[L-1,R]$这个区间。用 $cnt[d]$表示 $[L,R]$之间 $d$的倍数。
发现该区间比上一个区间多了一个$\sum_{i=L}^Rgcd(a[L-1],a[i])$。

由上面的反演可得到：

$$\sum_{i=L}^Rgcd(a[L-1],a[i]) =\sum_{d|a[L-1]}\phi(d)*\sum_{i=L,d|a[i]}^R1 =\sum_{d|a[L-1]}\phi(d)*cnt[d]$$

此时就可以做出答案，当然需要找出$a[L-1]$的所有约数。如果$a[i]<=10^6$那么$a[i]$最多有$100$个约数，那么每一次暴力枚举需要$o(100)$，但是这个常数很难达到。
最大的复杂度是$o(n\sqrt{n}*100)$但是因为100很难达到，所有可行。

#include <bits/stdc++.h>
#define LL long long
#define FOR(i,x,y)  for(int i = x;i < y;++ i)
#define IFOR(i,x,y) for(int i = x;i > y;-- i)

using namespace std;

const int maxN = 1000010;
const int maxn = 20020;

int Blocks;
int prime[maxN],phi[maxN];
bool check[maxN];

void Eular(){
    memset(check,false,sizeof(check));
    prime[0] = 0;
    phi[1] = 1;
    FOR(i,2,maxn){
        if(!check[i]){
            prime[++prime[0]] = i;
            phi[i] = i-1;
        }
        FOR(j,1,prime[0]+1){
            if(i*prime[j] >= maxn) break;
            check[i*prime[j]] = true;
            if(i%prime[j]){
                phi[i*prime[j]] = phi[i]*(prime[j]-1);
            }
            else{
                phi[i*prime[j]] = phi[i]*prime[j];
                break;
            }
        }
    }
}

struct Commends{
    int l,r;
    int id;
    bool operator < (const Commends& rhs) const{
        if(l/Blocks == rhs.l/Blocks)    return r < rhs.r;
        return l/Blocks < rhs.l/Blocks;
    }
}cmd[maxn];

int cnt[maxN],n,q;
int pri[maxn],pri_cnt;
vector <int> fac[maxn];

void Get_Fac(int i){
    int len = (1<<pri_cnt);
    FOR(k,0,len){
        int res = 1;
        FOR(j,0,pri_cnt){
            if(k&(1<<j)){
                res *= pri[j];
            }
        }
        fac[i].push_back(res);
    }
    sort(fac[i].begin(),fac[i].end());
    vector <int> :: iterator it = unique(fac[i].begin(),fac[i].end());
    fac[i].erase(it,fac[i].end());
}

LL ans[maxn];

void init(){
    scanf("%d",&n);
    int num;
    memset(cnt,0,sizeof(cnt));
    FOR(i,0,maxn)   fac[i].clear();
    FOR(i,1,n+1){
        scanf("%d",&num);
        pri_cnt = 0;
        FOR(j,1,prime[0]+1){
            if(prime[j]*prime[j] > num) break;
            while(num%prime[j] == 0){
                pri[pri_cnt++] = prime[j];
                num /= prime[j];
            }
        }
        if(num > 1)   pri[pri_cnt++] = num;
        Get_Fac(i);
    }
    scanf("%d",&q);
    FOR(i,0,q){
        scanf("%d%d",&cmd[i].l,&cmd[i].r);
        cmd[i].id = i;
    }
    Blocks = (int)sqrt(n+0.5);
    sort(cmd,cmd+q);
}

void MO(){
    int L=1,R=0;
    LL res = 0;
    FOR(i,0,q){
        while(L > cmd[i].l){
            L --;
            FOR(j,0,(int)fac[L].size()){
                res += cnt[fac[L][j]]*(LL)phi[fac[L][j]];
                cnt[fac[L][j]] ++;
            }
        }
        while(L < cmd[i].l){
            FOR(j,0,(int)fac[L].size()){
                int d = fac[L][j];
                cnt[fac[L][j]] --;
                res -= cnt[fac[L][j]]*(LL)phi[fac[L][j]];
            }
            L ++;
        }
        while(R > cmd[i].r){
            FOR(j,0,(int)fac[R].size()){
                cnt[fac[R][j]] --;
                res -= cnt[fac[R][j]]*(LL)phi[fac[R][j]];
            }
            R --;
        }
        while(R < cmd[i].r){
            R ++;
            FOR(j,0,(int)fac[R].size()){
                res += cnt[fac[R][j]]*(LL)phi[fac[R][j]];
                cnt[fac[R][j]] ++;
            }
        }
        ans[cmd[i].id] = res;
    }
    FOR(i,0,q){
        printf("%I64d\n",ans[i]);
    }
}

int main()
{
    //freopen("test.in","r",stdin);
    int T,tCase = 0;  scanf("%d",&T);
    Eular();
    while(T--){
        printf("Case #%d:\n",++tCase);
        init();
        MO();
    }
    return 0;
}