1. 程式人生 > 實用技巧 >[SDOI2015]約數個數和

[SDOI2015]約數個數和

題意

\(\sum\limits_i^n\ \sum\limits_j^m\ d(ij)\)其中\(d(x)\)\(x\)的約數個數

想法

看到這個柿子想到的就是莫反了。
有這樣一個結論

\( d(xy)= \sum\limits_{i|x} \ \sum\limits_{j|y}[gcd(i,j) == 1]  \)

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

\(\sum\limits_i^n\ \sum\limits_j^m\ d(ij)\)

\(\sum\limits_i^n\ \sum\limits_j^m\ \sum\limits_{x|i} \ \sum\limits_{y|j}[gcd(x,y) == 1])\)

\(\mu(x)\) 的性質:

\([n == 1] = \sum\limits_{d|n}\mu(d)\)

\(\sum\limits_i^n\ \sum\limits_j^m\ \sum\limits_{x|i} \ \sum\limits_{y|j}\sum\limits_{d|gcd(x,y)}\mu(d))\)

改為列舉\(d\)-

\(\sum\limits_{d = 1}^n\sum\limits_{x = 1}^{n/d}\sum\limits_{i = 1}^{n/d/x}\sum\limits_{y = 1}^{m/d}\sum\limits_{j = 1}^{m/d/y}\mu(d)\)

\(\sum\limits_{d = 1}^n\sum\limits_{x = 1}^{n/d}\lfloor\frac{n}{d * x}\rfloor\sum\limits_{y = 1}^{n/d}\lfloor\frac{m}{d * y}\rfloor\mu(d)\)

\(\sum\limits_{d = 1}^n\mu(d)\sum\limits_{x = 1}^{n/d}\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{x}\rfloor\sum\limits_{y = 1}^{m/d}\lfloor\frac{\lfloor\frac{m}{d}\rfloor}{y}\rfloor\)

後面兩個\(\sum\)可以整除分塊

複雜度\(O(N\sqrt(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\n", solve(a, b));
	}
	return 0;
}