bzoj3529: [Sdoi2014]数表 莫比乌斯反演

bzoj3529: [Sdoi2014]数表

Description

有一张N×m的数表，其第i行第j列（1 < =i < =礼，1 < =j < =m）的数值为

能同时整除i和j的所有自然数之和。给定a，计算数表中不大于a的数之和。

Input

输入包含多组数据。
输入的第一行一个整数Q表示测试点内的数据组数，接下来Q行，每行三个整数n，m，a(|a| < =10^9)描述一组数据。

Output

对每组数据，输出一行一个整数，表示答案模2^31的值。

Sample Input

2
4 4 3
10 10 5

Sample Output

20
148

HINT

1 < =N．m < =10^5 ， 1 < =Q < =2×10^4

分析

首先，我们令f(i,j)=∑d|i⋀d|jd

则ans=∑in∑im[f(i,j)≤a]f(i,j)

考虑f，不难发现f(i.j)=∑d|gcd(i,j)d=
1aa70
σ(gcd(i,j))

其中σ(i)为i的约束和

首先咱们先不考虑a的限制

ans=∑in∑imσ(gcd(i,j))

然后枚举约数

ans=∑dnσ(d)∑in∑im[gcd(i,j))==d]

根据模型(详见：bzoj2820YY的GCD)得

∑in∑im[gcd(i,j))==d]=∑k⌊nd⌋μ(k)⌊nkd⌋⌊mkd⌋

于是有

ans=∑dnσ(d)∑k⌊nd⌋μ(k)⌊nkd⌋⌊mkd⌋

转换枚举变量，令D=kd

ans=∑Dn⌊nD⌋⌊mD⌋∑d|Dσ(d)μ(Dd)

这个时候，我们发现我们还没有考虑a的限制。。

我们考虑最开始，发现a本身的限制就是加在f(i,j)也就是σ(d)上的，所以我们最后的答案可以直接写成

ans=∑Dn⌊nD⌋⌊mD⌋∑d|Dσ(d)μ(Dd)[σ(d)≤a]

那么，对于∑Dn⌊nD⌋⌊mD⌋可以使用下底函数分块。

而对于∑d|Dσ(d)μ(Dd)[σ(d)≤a] 处理条件[σ(d)≤a]，我们离线将a排序，从小到大把σ(d)μ(Dd)的插入树状数组即可。使用树状数组的原因

一是随着a的改变这玩意的值也会改变，需要支持单点修改。二是分块求答案的时候需要做一个区间求和。……时间复杂度O(Tn√logn)。——form xxcc

代码

具体细节看代码，本题还是有很多值得推敲之处的。

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;
const int N = 110000;
const int M = 440000;
const int inf = 0x7fffffff;
int read() {
char ch = getchar(); int x = 0, f = 1;
while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
while(ch >= '0' && ch <= '9') {x = x * 10 + ch - '0'; ch = getchar();}
return x * f;
}
pair<int, int>f
;
int mu
, p
, t
, a
, id
, qn
, qm
, ans
, mx, tot, top, m;
bool vis
;
void add(int x, int add) {for(int i = x; i <= mx; i += i & -i) t[i] += add;}
int query(int x) {int ret = 0; for(int i = x; i; i -= i & -i) ret += t[i]; return ret;}
bool cmp(int b, int c) {return a[b] < a[c];}

void mobius() {
mu[1] = 1;
for(int i = 2;i <= mx; ++i) {
if(!vis[i]) {p[++tot] = i; mu[i] = -1;}
for(int j = 1;j <= tot && i * p[j] <= mx; ++j) {
vis[i * p[j]] = true;
if(i % p[j]) mu[i * p[j]] = -mu[i];
else break;
}
}
for(int i = 1;i <= mx; ++i)
for(int j = i;j <= mx; j += i)
f[j].first += i;
for(int i = 1;i <= mx; ++i) f[i].second = i;
}

void work(int id) {
while(top <= mx && f[top].first <= a[id]) {
for(int i = f[top].second; i <= mx; i += f[top].second)
add(i, f[top].first * mu[i / f[top].second]);
++top;
}
int n = qn[id], m = qm[id];
for(int i = 1, pos;i <= n; i = pos + 1) {
pos = min(n / (n / i), m / (m / i));
ans[id] += (n / i) * (m / i) * (query(pos) - query(i - 1));
}

}

int main() {
m = read();
for(int i = 1;i <= m; ++i) {
qn[i] = read(); qm[i] = read(); if(qn[i] > qm[i]) swap(qn[i], qm[i]);
a[i] = read(); id[i] = i; mx = max(mx, qn[i]);
}
mobius();
sort(id + 1, id + m + 1, cmp);
sort(f + 1, f + mx + 1); top = 1;
for(int i = 1;i <= m; ++i) work(id[i]);
for(int i = 1;i <= m; ++i) printf("%d\n", ans[i] & inf);
return 0;
}