[BZOJ3309]DZY Loves Math

试题描述

对于正整数 (n)，定义 (f(n)) 为 (n) 所含质因子的最大幂指数。例如 (f(1960)=f(2^3 imes 5^1 imes 7^2)=3, f(10007)=1, f(1)=0)。

给定正整数 (a,b)，求 (sum_{i=1}^n sum_{j=1}^m f(gcd(i, j)))。

输入

第一行一个数 (T)，表示询问数。

接下来 (T) 行，每行两个数 (a,b)，表示一个询问。

输出

对于每一个询问，输出一行一个非负整数作为回答。

输入示例

4
7558588 9653114
6514903 4451211
7425644 1189442
6335198 4957

输出示例

35793453939901
14225956593420
4332838845846
15400094813

数据规模及约定

(T le 10000)

(1 le a,b le 10^7)

题解

按照套路来，枚举最大公约数（以下有些过程的求和符号的上界应该是 (lfloor frac{min {n, m}}{x} floor)，不失正确性，以下都写成 (lfloor frac{n}{x} floor)）

[sum_{i=1}^n sum_{j=1}^m f(gcd(i, j)) \ = sum_{g=1}^n { f(g) sum_{i=1}^n sum_{j=1}^m [gcd(i, j) = g] } \ = sum_{g=1}^n { f(g) sum_{i=1}^{lfloor frac{n}{g} floor} sum_{j=1}^{lfloor frac{m}{g} floor} sum_{d|i, d|j} mu(d) } \ = sum_{g=1}^n { f(g) sum_{d=1}^{lfloor frac{n}{g} floor} mu(d) lfloor frac{n}{gd} floor lfloor frac{m}{gd} floor } ]

继续按套路来，令 (T = gd)

[= sum_{T=1}^n { lfloor frac{n}{T} floor lfloor frac{m}{T} floor sum_{d|T} mu(d) f left( frac{T}{d} ight) } ]

令 (g(T) = sum_{d|T} mu(d) f left( frac{T}{d} ight))，那么现在的问题就是快速求出 (g(T)) 的前缀和。

[sum_{T=1}^n g(T) \ = sum_{T=1}^n sum_{d|T} mu(d) f left( frac{T}{d} ight) \ = sum_{d=1}^n mu(d) sum_{T=1}^{lfloor frac{n}{d} floor} ]

(f(T), mu(d)) 的前缀和可以线性筛得到，那么我们只需要预处理 (g(T)) 的前 (n^{frac{2}{3}}) 项就可以杜教筛了。这个做法总时间复杂度 (O(T n^{frac{2}{3}}))，在某些 OJ 上能过……

然而 BZOJ 上是不行的，所以需要想办法线性筛出 (g(T))。

接下来是不知道怎么想出来的部分。

令 (T = prod_{i=1}^k p_i^{a_i})，由于 (mu(d)) 不为 (0) 的 (d) 的每个质因子都只有 (1) 次幂，所以有如下结论：

若 (exists i e j, a_i e a_j)，那么 (g(T) = 0)，这个可以用 ((1-1)^k) 的二项式展开证明：把贡献分为 (d) 是否含有所有的指数最大的质因子，两类贡献都是形如 (k') 个数中取偶数个数的方案减去 (k') 个数中取奇数个数的方案，它们都是 (0)；
否则 (g(T) = (-1)^{k+1})，如果有偶数个不同的质因子，会在“全选”的时候少加 (1)，因为 (f) 值变小了 (1)（注意这里是所有 (a_i) 都相等的情况），同理奇数个不同质因子会在“全选”时少减 (1)。

那么就可以线性筛出 (g(T)) 了，复杂度降为 (O(T sqrt n))。

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)

const int BufferSize = 1 << 16;
char buffer[BufferSize], *Head, *Tail;
inline char Getchar() {
	if(Head == Tail) {
		int l = fread(buffer, 1, BufferSize, stdin);
		Tail = (Head = buffer) + l;
	}
	return *Head++;
}
int read() {
	int x = 0, f = 1; char c = Getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = Getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = Getchar(); }
	return x * f;
}

#define maxn 10000010
#define LL long long

bool vis[maxn];
int prime[maxn], cp, mu[maxn], f[maxn], g[maxn], smu[maxn], nos[maxn];
LL sf[maxn], sg[maxn];
void init() {
	int n = maxn - 1;
	mu[1] = smu[1] = 1; f[1] = sf[1] = g[1] = sg[1] = 0;
	rep(i, 2, n) {
		if(!vis[i]) prime[++cp] = i, mu[i] = -1, f[i] = 1, nos[i] = 1, g[i] = 1;
		for(int j = 1; j <= cp && i * prime[j] <= n; j++) {
			vis[i*prime[j]] = 1;
			if(i % prime[j] == 0) {
				mu[i*prime[j]] = 0;
				f[i*prime[j]] = max(f[nos[i]], f[i/nos[i]] + 1);
				if(nos[i] > 1) g[i*prime[j]] = f[nos[i]] == f[i/nos[i]] + 1 ? -g[nos[i]] : 0;
				else g[i*prime[j]] = 1;
				nos[i*prime[j]] = nos[i];
				break;
			}
			mu[i*prime[j]] = -mu[i];
			f[i*prime[j]] = max(f[i], 1);
			g[i*prime[j]] = f[i] == 1 ? -g[i] : 0;
			nos[i*prime[j]] = i;
		}
		smu[i] = smu[i-1] + mu[i];
		sf[i] = sf[i-1] + f[i];
		sg[i] = sg[i-1] + g[i];
	}
	return ;
}

int num[50], cntn;
void putLL(LL x) {
	if(!x) return (void)puts("0");
	cntn = 0;
	while(x) num[cntn++] = x % 10, x /= 10;
	dwn(i, cntn - 1, 0) putchar(num[i] + '0'); putchar('
');
	return ;
}

void work() {
	int n = read(), m = read();
	LL ans = 0;
	for(int T = 1; T <= min(n, m); ) {
		int r = min(min(n / (n / T), m / (m / T)), min(n, m));
		ans += (LL)(n / T) * (m / T) * (sg[r] - sg[T-1]);
		T = r + 1;
	}
	putLL(ans);
	return ;
}

int main() {
	init();
	
	int T = read();
	
	while(T--) work();
	
	return 0;
}

再贴一下杜教筛的暴力

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;
#define rep(i, s, t) for(int i = (s), mi = (t); i <= mi; i++)
#define dwn(i, s, t) for(int i = (s), mi = (t); i >= mi; i--)

const int BufferSize = 1 << 16;
char buffer[BufferSize], *Head, *Tail;
inline char Getchar() {
	if(Head == Tail) {
		int l = fread(buffer, 1, BufferSize, stdin);
		Tail = (Head = buffer) + l;
	}
	return *Head++;
}
int read() {
	int x = 0, f = 1; char c = Getchar();
	while(!isdigit(c)){ if(c == '-') f = -1; c = Getchar(); }
	while(isdigit(c)){ x = x * 10 + c - '0'; c = Getchar(); }
	return x * f;
}

#define maxn 10000010
#define maxm 4000010
#define LL long long

bool vis[maxn];
int prime[maxn], cp, mu[maxn], f[maxn], smu[maxn], nos[maxn];
LL sf[maxn], g[maxm], sg[maxm];
void init() {
	int n = maxn - 1;
	mu[1] = smu[1] = 1; f[1] = sf[1] = 0;
	rep(i, 2, n) {
		if(!vis[i]) prime[++cp] = i, mu[i] = -1, f[i] = 1, nos[i] = 1;
		for(int j = 1; j <= cp && i * prime[j] <= n; j++) {
			vis[i*prime[j]] = 1;
			if(i % prime[j] == 0) {
				mu[i*prime[j]] = 0;
				f[i*prime[j]] = max(f[nos[i]], f[i/nos[i]] + 1);
				nos[i*prime[j]] = nos[i];
				break;
			}
			mu[i*prime[j]] = -mu[i];
			f[i*prime[j]] = max(f[i], 1);
			nos[i*prime[j]] = i;
		}
		smu[i] = smu[i-1] + mu[i];
		sf[i] = sf[i-1] + f[i];
	}
	n = maxm - 1;
	rep(i, 1, n) for(int j = 1; i * j <= n; j++) g[i*j] += mu[i] * f[j];
	rep(i, 1, n) sg[i] = sg[i-1] + g[i];
	return ;
}

const int HMOD = 1000037;
struct Hash {
	int ToT, head[HMOD], nxt[maxn], key[maxn];
	LL val[maxn];
	Hash() { ToT = 0; memset(head, 0, sizeof(head)); }
	LL Find(int x) {
		int u = x % HMOD;
		for(int e = head[u]; e; e = nxt[e]) if(key[e] == x) return val[e];
		return -1;
	}
	void Insert(int x, LL v) {
		int u = x % HMOD;
		key[++ToT] = x; val[ToT] = v; nxt[ToT] = head[u]; head[u] = ToT;
		return ;
	}
} hG;

LL G(int n) {
	if(n < maxm) return sg[n];
	LL getg = hG.Find(n);
	if(getg >= 0) return getg;
	LL ans = 0;
	for(int d = 1; d <= n; ) {
		int r = min(n / (n / d), n);
		ans += sf[n/d] * (smu[r] - smu[d-1]);
		d = r + 1;
	}
	hG.Insert(n, ans);
	return ans;
}

int num[50], cntn;
void putLL(LL x) {
	if(!x) return (void)puts("0");
	cntn = 0;
	while(x) num[cntn++] = x % 10, x /= 10;
	dwn(i, cntn - 1, 0) putchar(num[i] + '0'); putchar('
');
	return ;
}

void work() {
	int n = read(), m = read();
	LL ans = 0;
	for(int T = 1; T <= min(n, m); ) {
		int r = min(min(n / (n / T), m / (m / T)), min(n, m));
		ans += (LL)(n / T) * (m / T) * (G(r) - G(T - 1));
		T = r + 1;
	}
	putLL(ans);
	return ;
}

int main() {
	init();
	
	int T = read();
	
	while(T--) work();
	
	return 0;
}

~~可以看出卡常的痕迹。~~