bzoj3529/SDOI2014 数表（莫比乌斯函数+离散化+树状数组

题目描述

有一张N*m的数表，其第i行第j列（1 < =i < =n，1 < =j < =m）的数值为能同时整除i和j的所有自然数之和。给定a，计算数表中不大于a的数之和。

输入格式

输入包含多组数据。

输入的第一行一个整数Q表示测试点内的数据组数

接下来Q行，每行三个整数n，m，a(|a| < =10^9)描述一组数据。

输出格式

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

输入输出样例

输入 #1

2
4 4 3
10 10 5

输出 #1

20
148
说明/提示
1 < =N．m < =10^5 ， 1 < =Q < =2*10^4

分析

这题有点难度啊。。。
记 $f(i)$ 表示 $i$ 的约数和
题目要求的是
$ans=\sum _{i=1}^n\sum _{j=1}^m\sum _{d|i,d|j}d\ast (\sum _{d|i,d|j}d\le a)=\sum _{i=1}^n\sum _{j=1}^{m}f(gcd(i,j))\ast (f(gcd(i,j))\le a)$
先不考虑 $a$ 的限制
我们先求求 $\sum _{i=1}^n\sum _{j=1}^{m}f(gcd(i,j))$
$\sum _{i=1}^n\sum _{j=1}^{m}f(gcd(i,j))=\sum _{d=1}^{n}f(d)\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)==d]=\sum _{d=1}^{n}f(d)\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }\sum _{x|i,x|j}\mu (x)=\sum _{d=1}^{n}f(d)\sum _{x=1}^{\lfloor \frac{n}{d}\rfloor }\mu (x)\lfloor \frac{n}{dx}\rfloor \lfloor \frac{m}{dx}\rfloor$
然后我们枚举 $d$ 与 $x$ 的乘积
那么原式等价于求
$\sum _{T=1}^n\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor \sum _{d|T}f(d)\mu (\frac{T}{d})$
记 $h(T)=\sum _{d|T}f(d)\ast \mu (\frac{T}{d})$
那么我们要求的就是
$\sum _{T=1}^n\lfloor \frac{n}{T}\rfloor \lfloor \frac{m}{T}\rfloor h(T)$
现实总是很残酷的，给了个 $a$ 的限制。也就是我们对于每个 $a$ 只有小于 $a$ 的 $f(d)$ 可以对 $h(T)$ 产生贡献。
我们每次修改都是再修改 $h(T)$ ，我们每次询问，都是要求 $h(T)$ 的前缀和，而 $h(T)$ 是在变化的，于是我们想到用树状数组来维护 $h(T)$。
所以我们先将 $a$ 排序，每次将小于 $a$ 的 $f(d)$ 加入贡献，然后询问的时候除法分块就好了。
最后，由于模数是 $2^{31}$，自然溢出就好了。

代码如下

#include <bits/stdc++.h>
#define N 100005
#define LL long long
#define inf (1 << 31)
using namespace std;
int tree[N], ans[N], x[N], p[N], f[N], mu[N], g[N], pp[N], cnt, maxn = 100000;
LL z = 1, t;
struct node{
	int a, i, n, m;
	bool operator < (const node & A) const{
		return a < A.a;
	}
}d[N], e[N];
void add(int i, int x){
	for(; i <= maxn; i += i&-i) tree[i] += x;
}
int query(int i){
	int s = 0;
	for(; i; i -= i&-i) s += tree[i];
	return s;
}
int main(){
	int i, j, b, k, n, m, q, a, l, r;
	mu[1] = 1;
	for(i = 2; i <= maxn; i++){
		if(!x[i]) x[i] = i, p[++cnt] = i, mu[i] = -1;
		for(j = 1; j <= cnt; j++){
			t = z * i * p[j];
			if(t > maxn) break;
			x[t] = p[j];
			if(i % p[j] == 0) break;
			mu[t] = -mu[i];
		}
	}
	for(i = 1; i <= maxn; i++){
		for(j = 1; z * i * j <= maxn; j++){
			f[i * j] += i;
		}
		d[i].a = f[i]; d[i].i = i;
	}
	sort(d + 1, d + i);
	scanf("%d", &q);
	for(i = 1; i <= q; i++){
		scanf("%d%d%d", &e[i].n, &e[i].m, &e[i].a);
		e[i].i = i;
	}
	sort(e + 1, e + i);
	for(i = 1; i <= q; i++) pp[e[i].i] = i;
	j = 1;
	for(i = 1; i <= q; i++){
		while(j <= maxn && d[j].a <= e[i].a){
			a = d[j].a;
			b = d[j].i;
			for(k = 1; z * b * k <= maxn; k++){
				add(b * k, a * mu[k]);
			}
			j++;
		}
		n = e[i].n; m = e[i].m;
		if(n > m) swap(n, m);
		for(l = 1; l <= n; l = r + 1){
			r = min(n / (n / l), m / (m / l));
			ans[i] += (n / l) * (m / l) * (query(r) - query(l - 1));
		}
		if(ans[i] < 0) ans[i] += inf;
	}
	for(i = 1; i <= q; i++) printf("%d\n", ans[pp[i]]);
	return 0;
}