Codeforces Round #146 E. Number Challenge 莫比烏斯反演

Description

定義 d ( i ) d(i) d(i)為 i i i的約數個數，求：

∑ i = 1 a ∑ j = 1 b ∑ k = 1 c d ( i j k ) \sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^cd(ijk) i=1∑a​j=1∑b​k=1∑c​d(ijk)

Sample Input

2 2 2

Sample Output

20

首先有一個性質：

∑ i = 1 a ∑ j = 1 b d ( i j ) = ∑ i = 1 a ∑ j = 1 b ∑ x ∣ i ∑ y ∣ j [ g c d ( x , y ) = = 1 ] \sum_{i=1}^a\sum_{j=1}^bd(ij)=\sum_{i=1}^a\sum_{j=1}^b\sum_{x|i}\sum_{y|j}[gcd(x,y)==1] i=1∑a​j=1∑b​d(ij)=i=1∑a​j=1∑b​x∣i∑​y∣j∑​[gcd(x,y)==1]

根據這個式子，可以推導：

d ( i j k ) = ∑ x ∣ i ∑ y ∣ j k [ g c d ( x , y ) = = 1 ] = ∑ x ∣ i ∑ y ∣ j ∑ z ∣ k [ g c d ( x , y z ) = = 1 ] [ g c d ( y , z ) = = 1 ] d(ijk)=\sum_{x|i}\sum_{y|jk}[gcd(x,y)==1]=\sum_{x|i}\sum_{y|j}\sum_{z|k}[gcd(x,yz)==1][gcd(y,z)==1] d(ijk)=x∣i∑​y∣jk∑​[gcd(x,y)==1]=x∣i∑​y∣j∑​z∣k∑​[gcd(x,yz)==1][gcd(y,z)==1]

= ∑ x ∣ i ∑ y ∣ j ∑ z ∣ k [ g c d ( x , y ) = = 1 ] [ g c d ( x , z ) = = 1 ] [ g c d ( y , z ) = = 1 ] =\sum_{x|i}\sum_{y|j}\sum_{z|k}[gcd(x,y)==1][gcd(x,z)==1][gcd(y,z)==1] =x∣i∑​y∣j∑​z∣k∑​[gcd(x,y)==1][gcd(x,z)==1][gcd(y,z)==1]

那麼，我們将化出來的帶進式子中：

∑ i = 1 a ∑ j = 1 b ∑ k = 1 c ∑ x ∣ i ∑ y ∣ j ∑ z ∣ k [ g c d ( x , y ) = = 1 ] [ g c d ( x , z ) = = 1 ] [ g c d ( y , z ) = = 1 ] \sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c\sum_{x|i}\sum_{y|j}\sum_{z|k}[gcd(x,y)==1][gcd(x,z)==1][gcd(y,z)==1] i=1∑a​j=1∑b​k=1∑c​x∣i∑​y∣j∑​z∣k∑​[gcd(x,y)==1][gcd(x,z)==1][gcd(y,z)==1]

= ∑ i = 1 a ⌊ a i ⌋ ∑ d μ [ d ] ∑ d ∣ j , g c d ( i , j ) = = 1 ⌊ b j ⌋ ∑ d ∣ k , g c d ( i , k ) = = 1 ⌊ c k ⌋ =\sum_{i=1}^a\lfloor \frac ai\rfloor\sum_d\mu[d]\sum_{d|j,gcd(i,j)==1}\lfloor \frac bj\rfloor\sum_{d|k,gcd(i,k)==1}\lfloor \frac ck\rfloor =i=1∑a​⌊ia​⌋d∑​μ[d]d∣j,gcd(i,j)==1∑​⌊jb​⌋d∣k,gcd(i,k)==1∑​⌊kc​⌋

= ∑ i = 1 a ∑ j = 1 b ∑ k = 1 c ∑ x ∣ i ∑ y ∣ j ∑ z ∣ k [ g c d ( x , y ) = = 1 ] [ g c d ( x , z ) = = 1 ] [ g c d ( y , z ) = = 1 ] =\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c\sum_{x|i}\sum_{y|j}\sum_{z|k}[gcd(x,y)==1][gcd(x,z)==1][gcd(y,z)==1] =i=1∑a​j=1∑b​k=1∑c​x∣i∑​y∣j∑​z∣k∑​[gcd(x,y)==1][gcd(x,z)==1][gcd(y,z)==1]

= ∑ i = 1 a ⌊ a i ⌋ ∑ d μ [ d ] ∑ g c d ( i , d j ) = = 1 ⌊ b d j ⌋ ∑ g c d ( i , d k ) = = 1 ⌊ c d k ⌋ =\sum_{i=1}^a\lfloor \frac ai\rfloor\sum_d\mu[d]\sum_{gcd(i,dj)==1}\lfloor \frac b{dj}\rfloor\sum_{gcd(i,dk)==1}\lfloor \frac c{dk}\rfloor =i=1∑a​⌊ia​⌋d∑​μ[d]gcd(i,dj)==1∑​⌊djb​⌋gcd(i,dk)==1∑​⌊dkc​⌋

以 g c d ( i , d j ) = = 1 gcd(i,dj)==1 gcd(i,dj)==1為例，然後你考慮如何滿足。

枚舉 i i i，再枚舉 d d d，首先使 g c d ( i , d ) = = 1 gcd(i,d)==1 gcd(i,d)==1，

對于 j j j從 1 1 1枚舉到 b d \frac bd db​，使 g c d ( i , j ) = = 1 gcd(i,j)==1 gcd(i,j)==1。

根據調和計數，這樣的複雜度是 O ( n 2 l o g n ) O(n^2log n) O(n2logn)的。

#include <ctime>
#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>

using namespace std;
typedef long long LL;
const int N = 2001;
const int mod = 1073741824;
int _max(int x, int y) {return x > y ? x : y;}
int _min(int x, int y) {return x < y ? x : y;}
int read() {
	int s = 0, f = 1; char ch = getchar();
	while(ch < '0' || ch > '9') {if(ch == '-') f = -1; ch = getchar();}
	while(ch >= '0' && ch <= '9') s = s * 10 + ch - '0', ch = getchar();
	return s * f;
}

int plen, mu[N], p[N], gcd[N][N];
bool v[N];

int Gcd(int a, int b) {
	if(gcd[a][b]) return gcd[a][b];
	return b == 0 ? a : Gcd(b, a % b);
}

void get_p() {
	mu[1] = 1;
	for(int i = 2; i < N; i++) {
		if(!v[i]) p[++plen] = i, mu[i] = -1;
		for(int j = 1; j <= plen && (LL)p[j] * i < N; j++) {
			v[i * p[j]] = 1;
			if(i % p[j] == 0) {mu[i * p[j]] = 0; break;}
			mu[i * p[j]] = -mu[i];
		}
	}
}

LL calc(int n, int x) {
	LL sum = 0;
	for(int i = 1; i <= n; i++) if(gcd[i][x] == 1) (sum += n / i) %= mod;
	return sum;
}

int main() {
	get_p();
	int a = read(), b = read(), c = read();
	int u = _max(a, _max(b, c));
	for(int i = 1; i <= u; i++) {
		for(int j = i; j <= u; j++) {
			gcd[i][j] = Gcd(i, j);
			gcd[j][i] = gcd[i][j];
		}
	}
	LL ans = 0;
	for(int i = 1; i <= a; i++) {
		u = _min(b, c);
		LL sum = 0;
		for(int d = 1; d <= u; d++) if(gcd[i][d] == 1){
			(sum += (LL)mu[d] * calc(b / d, i) % mod * calc(c / d, i) % mod) %= mod;
		} (ans += sum * (a / i) % mod) %= mod;
	} printf("%lld\n", (ans + mod) % mod);
	return 0;
}