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

来看答案怎么求？

∑i=LR∑j=i+1Rgcd(a[i],a[j])=∑i=LR∑j=i+1R∑d|gcd(a[i],a[j])ϕ(d)=∑dϕ(d)∑i=L,d|a[i]R∑j=i+1,d|a[j]R1\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][L,R]这段区间内有多少个dd的倍数。

cnt[d][x]cnt[d][x]代表前x项中d的倍数的个数。

原式可化为

∑i=LR∑j=i+1Rgcd(a[i],a[j])=∑dϕ(d)∗(cnt[d][R]−cnt[d][L−1])\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]cnd[d]的值。

不妨假设[L,R][L,R]到[L−1,R][L-1,R]这个区间。用cnt[d]cnt[d]表示[L,R][L,R]之间dd的倍数。

发现该区间比上一个区间多了一个∑Ri=Lgcd(a[L−1],a[i])\sum_{i=L}^Rgcd(a[L-1],a[i])。

由上面的反演可得到：

∑i=LRgcd(a[L−1],a[i])=∑d|a[L−1]ϕ(d)∗∑i=L,d|a[i]R1=∑d|a[L−1]ϕ(d)∗cnt[d]\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[L-1]的所有约数。如果a[i]<=106a[i]<=10^6那么a[i]a[i]最多有100100个约数，那么每一次暴力枚举需要o(100)o(100)，但是这个常数很难达到。

最大的复杂度是o(nn−−√∗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;
}