[BZOJ3309]DZY Loves Math

[BZOJ3309]DZY Loves Math

试题描述

对于正整数 \(n\)，定义 \(f(n)\) 为 \(n\) 所含质因子的最大幂指数。例如 \(f(1960)=f(2^3 \times 5^1 \times 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} \rfloor\)，不失正确性，以下都写成 \(\lfloor \frac{n}{x} \rfloor\)）

\[\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} \rfloor} \sum_{j=1}^{\lfloor \frac{m}{g} \rfloor} \sum_{d|i, d|j} \mu(d) } \\ = \sum_{g=1}^n { f(g) \sum_{d=1}^{\lfloor \frac{n}{g} \rfloor} \mu(d) \lfloor \frac{n}{gd} \rfloor \lfloor \frac{m}{gd} \rfloor } \]

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

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

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

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

\(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 \ne j, a_i \ne 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('\n');
	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('\n');
	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;
}

可以看出卡常的痕迹。