[SDOI2014]数表

( ext{Problem})

求

[sum_{i=1}^n sum_{j=1}^m sigma_1(gcd(i,j)) ]

当且仅当 (sigma(gcd(i,j)) le a) 时有贡献
多组询问，答案对 (2^{31}) 取模
(1 le n,m le 10^5,1 le Q le 2 imes 10^4)

( ext{Analysis})

先不考虑 (a) 的限制
然后套路地用莫反推出答案式子

[sum_{T=1}^{min(n,m)} lfloor frac{n}{T} floor lfloor frac{m}{T} floor sum_{d mid T} sigma_1(d) mu(frac{T}{d}) ]

设 (G(T)=sum_{d mid T} sigma_1(d) mu(frac{T}{d}))

加上 (a) 的限制，只有 (sigma_1(d) le a) 时这个 (d) 才会对 (G) 做贡献
考虑离线，把询问按 (a) 排序
把 (sigma_1(d)) 排序
那么我们需要动态把 (d) 的贡献加入 (G)，并且不需要加前面加过的
这样每个 (d) 就只会加一次，枚举它的倍数，这样就能做到 (O(nlog n))
但我们数论分块需要求出 (G) 的前缀和
因为每次 (a) 不同，(G) 值被更新，我们不可以每次重新计算 (G) 的前缀和
所以只要借助能快速动态维护前缀和的工具
树状数组！

修改就做到 (O(n log^2 n)) 了
那总的复杂度就是 (O(q sqrt n log n + n log^2 n))
取模因为是 (2^31)，直接自然溢出即可，最后答案处理一下 Ans = Ans & (2^31-1)
这是一个小 (trick)
但我如果没用，洛谷上能 (A) （因为时限为 (2s)），但 (JZOJ) 上会 (T)

( ext{Code})

#include <cstdio>
#include <algorithm>
#define re register
using namespace std;
typedef unsigned int uint;

const int N = 100000, INF = 1e9 + 1;
int q, n, m, a, lg;
int vis[N + 5], mu[N + 5], pr[N], tp, gd[N + 5], g[N + 5];
uint ans[N + 5];

struct node{int n, m, a, id;}Q[20005];
inline bool cmp(node x, node y){return x.a < y.a;}
inline bool cmpg(int x, int y){return g[x] < g[y];}

void Pre()
{
	vis[1] = mu[1] = 1;
	for(re int i = 1; i <= N; i++)
	{
		if (!vis[i]) pr[++tp] = i, mu[i] = -1;
		for(re int j = 1; j <= tp && pr[j] * i <= N; j++)
		{
			vis[pr[j] * i] = 1;
			if (!(i % pr[j])) break;
			mu[pr[j] * i] = -mu[i];
		}
	}
	for(re int i = 1; i <= N; i++)
	if (g[i] ^ INF)
		for(re int j = i; j <= N; j += i)
		if (g[j] ^ INF)
		{
			g[j] += i;
			if (g[j] >= INF) g[j] = INF;
		}
	for(re int i = 1; i <= N; i++) gd[i] = i;
	sort(gd + 1, gd + N + 1, cmpg), lg = 1;
}

uint f[N + 5];
inline int lowbit(int x){return x & (-x);}
inline void add(int x, uint v){for(; x <= N; x += lowbit(x)) f[x] += v;}
inline uint query(int x)
{
	if (!x) return 0;
	uint res = 0;
	for(; x; x -= lowbit(x)) res += f[x];
	return res;
}

int main()
{
	freopen("table.in", "r", stdin);
	freopen("table.out", "w", stdout);
	scanf("%d", &q), Pre();
	for(re int i = 1; i <= q; i++) 
		scanf("%d%d%d", &Q[i].n, &Q[i].m, &Q[i].a), Q[i].id = i;
	sort(Q + 1, Q + q + 1, cmp);
	for(re int i = 1; i <= q; i++)
	{
		n = Q[i].n, m = Q[i].m, a = Q[i].a;
		while (lg <= N && g[gd[lg]] <= a)
		{
			for(re int j = gd[lg]; j <= N; j += gd[lg])
				add(j, (uint)g[gd[lg]] * mu[j / gd[lg]]);
			++lg;
		}
		uint res = 0;
		for(re int l = 1, r; l <= min(n, m); l = r + 1)
		{
			r = min(n / (n / l), m / (m / l));
			res = res + (uint)(n / l) * (m / l) * (query(r) - query(l - 1));
		}
		ans[Q[i].id] = res;
	}
	for(re int i = 1; i <= q; i++) printf("%u
", ans[i] & (2147483647));
}