[SDOI2015]约数个数和

题意

求(sumlimits_i^n sumlimits_j^m d(ij))其中(d(x))为(x)的约数个数

想法

看到这个柿子想到的就是莫反了。
有这样一个结论

( d（xy）= sumlimits_{i|x} sumlimits_{j|y}[gcd(i,j) == 1]　 )

所以变形原柿子(n > m)

(sumlimits_i^n sumlimits_j^m d(ij))

(sumlimits_i^n sumlimits_j^m sumlimits_{x|i} sumlimits_{y|j}[gcd(x,y) == 1]))

用(mu(x)) 的性质：

([n == 1] = sumlimits_{d|n}mu(d))

(sumlimits_i^n sumlimits_j^m sumlimits_{x|i} sumlimits_{y|j}sumlimits_{d|gcd(x,y)}mu(d)))

改为枚举(d)-

(sumlimits_{d = 1}^nsumlimits_{x = 1}^{n/d}sumlimits_{i = 1}^{n/d/x}sumlimits_{y = 1}^{m/d}sumlimits_{j = 1}^{m/d/y}mu(d))

(sumlimits_{d = 1}^nsumlimits_{x = 1}^{n/d}lfloorfrac{n}{d * x} floorsumlimits_{y = 1}^{n/d}lfloorfrac{m}{d * y} floormu(d))

(sumlimits_{d = 1}^nmu(d)sumlimits_{x = 1}^{n/d}lfloorfrac{lfloorfrac{n}{d} floor}{x} floorsumlimits_{y = 1}^{m/d}lfloorfrac{lfloorfrac{m}{d} floor}{y} floor)

后面两个(sum)可以整除分块

复杂度(O(Nsqrt(N)))

代码

#include<bits/stdc++.h>
using namespace std;
#define int long long
int read(){
	char cc = getchar(); int cn = 0, flus = 1;
	while(cc < '0' || cc > '9'){
    	if(cc == '-') flus = -flus;
		cc = getchar();
	}
	while(cc >= '0' && cc <= '9')
	    cn = cn * 10 + cc - '0', cc = getchar();
	return cn * flus;
}
const int N = 50000 + 5;
int u[N], sum[N], pm[N], n, f[N], top;
bool isp[N];
void init(){
	u[1] = 1;
	for(int i = 2; i <= n; i++){
		if(!isp[i]){
			pm[++top] = i;
			u[i] = -1;
		}
		for(int j = 1; j <= top; j++){
			int p = pm[j];
			if(p * i > n) break;
			if(i % p == 0){
				isp[p * i] = 1;
				break;
			}
			isp[p * i] = 1;
			u[p * i] = u[p] * u[i];
		}
	}
	int l, r;
	for(int i = 1; i <= n; i++){
		for(l = 1; l <= i; l = r + 1){
			r = i/(i/l);
			f[i] += (r - l + 1) * (i/l);
		}
		sum[i] = sum[i - 1] + u[i];
	}
}
int solve(int x, int y){
	if(x > y) swap(x, y);
	int ans = 0, l = 1, r;
	for(l = 1; l <= x; l = r + 1){
		r = min((x/(x/l)), (y/(y/l)));
		ans += (sum[r] - sum[l - 1]) * f[x/l] * f[y/l];
	}
	return ans;
}
signed main()
{
	n = 50000;
	init();
	int T = read(), a, b;
	while(T--){
		a = read(); b = read();
		printf("%lld
", solve(a, b));
	}
	return 0;
}