[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;
}可以看出卡常的痕迹。
