数论分块（一）
数论分块
数论分块是为了解决 \(f(n) = \sum\limits_{i = 1}^{n}\left\lfloor\frac{n}{i}\right\rfloor g(i)\) 形式求值的一种算法，他可以在 \(O(\sqrt{n})\) 的时间内求出上述和式。
我们可以发现 \(\left\lfloor\frac{n}{i}\right\rfloor\) 在一个区间内的取值是一样的，因此我们考虑求出这段区间 \([l, r]\)，并预处理出 \(g(n)\) 的前缀和，每次对一个区间段进行求值，就可以做到 \(O(\sqrt{n})\) 的复杂度。
证明：
-
当 \(0 < x \le \sqrt{n}\) 时，\(\left\lfloor\frac{n}{x}\right\rfloor\) 的取值不超过 \(\sqrt{n}\) 个，这是显然的。
-
当 \(x > \sqrt{n}\) 时，\(\left\lfloor\frac{n}{x}\right\rfloor\) 的取值不超过 \(\sqrt{n}\) 个，因为 \(0 < \left\lfloor\frac{n}{x}\right\rfloor < \sqrt{n}\) 是有限的。
于是我们可以通过迅速求出边界来解决问题，若左边界 \(l\) 已知，则右边界不难求出 \(r = \left\lfloor\frac{n}{\left\lfloor\frac{n}{l}\right\rfloor}\right\rfloor\)，然后令 \(l = r + 1\) 继续遍历。
证明：
-
设 \(k = \left\lfloor\frac{n}{l}\right\rfloor\)，则 \(k \le \frac{n}{l}\)，可知 \(l \le \left\lfloor\frac{n}{k}\right\rfloor\)。
-
\(r = \left\lfloor\frac{n}{\left\lfloor\frac{n}{l}\right\rfloor}\right\rfloor = \left\lfloor\frac{n}{k}\right\rfloor \ge \left\lfloor\frac{n}{\frac{n}{l}}\right\rfloor = l\)，故 \(r\) 为所有满足 \(i \in [l', r'], k = \left\lfloor\frac{n}{i}\right\rfloor\) 的最大值 \(r'\)。
参考代码如下：
for (int l = 1, r; l <= n; l = r + 1)
{
r = n / (n / l);
// 块的和 * 块的权
Sum += (g[r] - g[l - 1]) * (n / l);
}
题目示例
P3935 Calculating
若 \(x\) 分解质因数结果为 \(x=p_1^{k_1}p_2^{k_2}\cdots p_n^{k_n}\)，令\(f(x)=(k_1+1)(k_2+1)\cdots (k_n+1)\)，求 \(\sum_{i=l}^rf(i)\) 对 \(998\,244\,353\) 取模的结果。
\(1\le l \le 10^{14}\), \(1\le r \le 1.6\times 10^{14}\)
考虑任选每个因子，其乘积恰好为因数个数和的式子，因此对 \(f(n)\) 的前缀和求解方式我们可以转化为求 \([1,n]\) 中每个数的倍数有多少个，转化为求下面这个式子：
这就是我们上面所说数论分块的板子，带入即可。
参考代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod = 998244353;
ll calc(ll n)
{
ll res = 0;
for (ll l = 1, r; l <= n; l = r + 1)
{
r = min(n, n / (n / l));
res += (r - l + 1) % mod * (n / l % mod);
res %= mod;
}
return res;
}
int main()
{
ll n, m;
cin >> n >> m;
cout << (calc(m) - calc(n - 1) + mod) % mod;
return 0;
}
P2261 [CQOI2007] 余数求和
给出正整数 \(n\) 和 \(k\)，请计算
其中 \(k\bmod i\) 表示 \(k\) 除以 \(i\) 的余数。
对于 \(100\%\) 的数据，保证 \(1 \leq n, k \leq 10^9\)。
首先可以转化式子的形式 \(k \bmod i = k - \left\lfloor\frac{k}{i}\right\rfloor i\)，对于这个式子我们是可以通过数论分块快速求出的，因此答案为 \(nk - \sum\limits_{i = 1}^n\left\lfloor\frac{k}{i}\right\rfloor i\)，需要注意的是分块上界。
-
\(n > k\) 时，\(i > k\) 的部分为 \(0\)，不用考虑。
-
\(n < k\) 时，\(n < i \le k\) 的部分为 \(0\)，因此也不用考虑。
所以我们取边界为 \(\min(n, k)\) 即可。
参考代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
int main()
{
ll n, k;
cin >> n >> k;
ll sum = n * k;
n = min(n, k);
for (ll l = 1, r; l <= n; l = r + 1)
{
r = min(n, k / (k / l));
sum -= (l + r) * (r - l + 1) / 2 * (k / l);
}
cout << sum;
return 0;
}
CF1485C Floor and Mod
求 \(1\le a\le x,1\le b\le y\) 且 \(\lfloor\frac{a}{b}\rfloor =a\bmod b\) 的\((a,b)\) 个数。
同上可知，\(\lfloor\frac{a}{b}\rfloor = a - \lfloor\frac{a}{b}\rfloor b\)。
化简得：
对每个 \(b\) 枚举合法的 \(a\)，不妨枚举 \(b + 1\) 的倍数 \(k\)，即设 \(a = (b + 1) k\)，则 \(\left\lfloor\frac{a}{b}\right\rfloor = \left\lfloor\frac{(b + 1)k}{b}\right\rfloor = \left\lfloor k + \frac{k}{b}\right\rfloor\)，为满足上式成立，有满足条件：
\(k\) 的取值有限，所以答案为 \(\sum\limits_{b = 1}^{y} \min(b - 1, \left\lfloor\frac{x}{b + 1}\right\rfloor)\)。
由于两边是单调的，不妨记 \(B\) 为两式的分界点，记 \(B\) 为满足 \(\left\lfloor\frac{x}{b + 1}\right\rfloor \ge b - 1\) 的 \(b\) 的取值，很容易将答案化简：
而由 \(\left\lfloor\frac{x}{b + 1}\right\rfloor \ge b - 1\) 容易得到 \(B\) 在 \(\sqrt{x + 1}\) 附近，因此可以首先得到一个估计值然后再进行判断。
式子 \(b < B\) 的部分可以 \(O(1)\) 计算，\(b > B\) 的部分可以通过数论分块实现，令 \(l = b + 1\)，对 \([B, y + 1]\) 做数论分块即可。
总时间复杂度 \(O(T\sqrt{y})\)。
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
void solve()
{
ll x, y;
cin >> x >> y;
ll t = sqrt(x + 1);
while (x / (t + 1) > t - 1) t ++ ;
t = min(t, y + 1);
ll ans = (t - 1) * (t - 2) / 2;
for (ll l = t + 1, r; l <= min(x, y + 1); l = r + 1)
{
r = min(x / (x / l), y + 1);
ans += x / l * (r - l + 1);
}
cout << ans << "\n";
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
int T = 1;
cin >> T;
while (T -- ) solve();
return 0;
}
P3455 [POI2007] ZAP-Queries
\(T\) 组测试数据，每组测试数据给出 \(a,b,d\)，求满足 \(1 \leq x \leq a\)，\(1 \leq y \leq b\)，且 \(\gcd(x,y)=d\) 的二元组 \((x,y)\) 的数量。
对于全部的测试点，保证 \(1 \leq T \leq 5 \times 10^4\)，\(1 \leq d \leq a,b \leq 5 \times 10^4\)。
在推式子之前，你需要了解莫比乌斯反演。
其中 \([P]\) 表示艾弗森括号，当且仅当 \(P\) 为真值是结果为 \(1\)，否则为 \(0\)。
记 \(n = a\)，\(m = b\) 且 \(n \le m\)，原式翻译过来就是：
由于 \(d\) 同时整除 \(i, j\)，不妨枚举 \(i, j\) 为 \(d\) 的倍数，对所有变量同时除以 \(d\) 有：
来一发莫比乌斯反演：
令 \(n' = \left\lfloor\frac{n}{d}\right\rfloor\)，\(m' = \left\lfloor\frac{m}{d}\right\rfloor\) 有：
右边的式子可以做一个二维数论分块，复杂度是根号同阶的。
- 二维数论分块相当于把两组线段集的端点加入同一集合划分为若干连续子段。
时间复杂度 \(O(T\sqrt{n})\)。
参考代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5e4 + 10;
int a, b, c, d, k;
int st[N], primes[N], mu[N], cnt;
void init()
{
mu[1] = 1;
for (int i = 2; i < N; i ++ )
{
if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;
for (int j = 0; primes[j] * i < N; j ++ )
{
st[primes[j] * i] = 1;
if (i % primes[j] == 0) break;
mu[primes[j] * i] = -mu[i];
}
}
for (int i = 1; i < N; i ++ ) mu[i] += mu[i - 1];
}
ll calc(int n, int m)
{
if (n > m) swap(n, m);
ll res = 0;
for (int l = 1, r; l <= n; l = r + 1)
{
r = min(n / (n / l), m / (m / l));
res += 1ll * (mu[r] - mu[l - 1]) * (n / l) * (m / l);
}
return res;
}
void solve()
{
cin >> a >> b >> k;
a = a / k, b = b / k;
cout << calc(a, b) << "\n";
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
init();
int T = 1;
cin >> T;
while (T -- ) solve();
return 0;
}
P2522 [HAOI2011] Problem b
对于给出的 \(T\) 个询问，每次求有多少个数对 \((x,y)\)，满足 \(a \le x \le b\)，\(c \le y \le d\)，且 \(\gcd(x,y) = k\)，\(\gcd(x,y)\) 函数为 \(x\) 和 \(y\) 的最大公约数。
对于 \(100\%\) 的数据满足：\(1 \le T,k \le 5 \times 10^4\)，\(1 \le a \le b \le 5 \times 10^4\)，\(1 \le c \le d \le 5 \times 10^4\)。
双倍经验，如果记上面求解的函数为 \(f(n, m)\)，那么这题可以简单的通过容斥得到答案。
参考代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 5e4 + 10;
int a, b, c, d, k;
int st[N], primes[N], mu[N], cnt;
void init()
{
mu[1] = 1;
for (int i = 2; i < N; i ++ )
{
if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;
for (int j = 0; primes[j] * i < N; j ++ )
{
st[primes[j] * i] = 1;
if (i % primes[j] == 0) break;
mu[primes[j] * i] = -mu[i];
}
}
for (int i = 1; i < N; i ++ ) mu[i] += mu[i - 1];
}
ll calc(int n, int m)
{
if (n > m) swap(n, m);
ll res = 0;
for (int l = 1, r; l <= n; l = r + 1)
{
r = min(n / (n / l), m / (m / l));
res += 1ll * (mu[r] - mu[l - 1]) * (n / l) * (m / l);
}
return res;
}
void solve()
{
cin >> a >> b >> c >> d >> k;
a = (a - 1) / k, b = b / k, c = (c - 1) / k, d = d / k;
cout << calc(b, d) - calc(a, d) - calc(b, c) + calc(a, c) << "\n";
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(0), cout.tie(0);
init();
int T = 1;
cin >> T;
while (T -- ) solve();
return 0;
}
P2568 GCD
给定正整数 \(n\)，求 \(1\le x,y\le n\) 且 \(\gcd(x,y)\) 为素数的数对 \((x,y)\) 有多少对。
- 对于 \(100\%\) 的数据，保证 \(1\le n\le10^7\)。
直接开始推式子：
此时就可以求出答案了，筛出质数后枚举每个质数，求 \(\mu(n)\) 的前缀和后对右边进行数论分块，时间复杂度 \(O(\frac{n}{\ln{n}}\sqrt{n})\)，可以通过。
但我们依旧可以优化，令 \(T = pd\)，按 \(T\) 枚举，化简原式有：
记 \(f(T) = \sum\limits_{p \in \mathbb{P}, p \mid T}\mu\left(\frac{T}{p}\right) = \sum\limits_{p \in \mathbb{P}}\sum\limits_{k = 1}^{\left\lfloor\frac{T}{p}\right\rfloor}\mu(k)\)，容易发现 \(f(T)\) 是可以预处理的，预处理出 \(f(T)\) 的狄利克雷前缀和的时间复杂度是 \(O(n\log\log{n})\)。
每次询问做一次数论分块，时间复杂度是 \(O(n\log\log{n} + \sqrt{n})\)。
参考代码（枚举法）
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e7 + 10;
int st[N], primes[N], mu[N], cnt;
void init()
{
mu[1] = 1;
for (int i = 2; i < N; i ++ )
{
if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;
for (int j = 0; primes[j] * i < N; j ++ )
{
st[primes[j] * i] = 1;
if (i % primes[j] == 0) break;
mu[primes[j] * i] = -mu[i];
}
}
for (int i = 1; i < N; i ++ ) mu[i] += mu[i - 1];
}
int main()
{
init();
int n;
cin >> n;
ll res = 0;
for (int i = 0; primes[i] <= n; i ++ )
{
ll m = n / primes[i];
for (int l = 1, r; l <= m; l = r + 1)
{
r = m / (m / l);
res += 1ll * (mu[r] - mu[l - 1]) * (m / l) * (m / l);
}
}
cout << res;
return 0;
}
参考代码（预处理）
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e7 + 10;
int st[N], primes[N], mu[N], cnt;
ll f[N];
void init()
{
mu[1] = 1;
for (int i = 2; i < N; i ++ )
{
if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;
for (int j = 0; primes[j] * i < N; j ++ )
{
st[primes[j] * i] = 1;
if (i % primes[j] == 0) break;
mu[primes[j] * i] = -mu[i];
}
}
for (int i = 0; i < cnt; i ++ )
for (int j = 1; primes[i] * j < N; j ++ )
f[primes[i] * j] += mu[j];
for (int i = 1; i < N; i ++ ) f[i] += f[i - 1];
}
int main()
{
init();
int n;
cin >> n;
ll res = 0;
for (int l = 1, r; l <= n; l = r + 1)
{
r = n / (n / l);
res += 1ll * (f[r] - f[l - 1]) * (n / l) * (n / l);
}
cout << res;
return 0;
}
P2257 YY的GCD
给定 \(N, M\)，求 \(1 \leq x \leq N\)，\(1 \leq y \leq M\) 且 \(\gcd(x, y)\) 为质数的 \((x, y)\) 有多少对，\(T\) 组测试数据。
- \(T = 10^4\)，\(N, M \leq 10^7\)
这个东西和我们上面求的东西一致，式子要写作这样：
记 \(f(T) = \sum\limits_{p \in \mathbb{P}}\sum\limits_{k = 1}^{\left\lfloor\frac{T}{p}\right\rfloor}\mu(k)\)，求一个前缀和 \(O(n)\)，然后对左边进行二维数论分块。
因此求一次答案的复杂度是根号同阶的，总体时间复杂度 \(O(n\log\log{n} + T\sqrt{n})\)。
参考代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e7 + 10;
int st[N], primes[N], mu[N], cnt;
ll f[N];
void init()
{
mu[1] = 1;
for (int i = 2; i < N; i ++ )
{
if (!st[i]) primes[cnt ++ ] = i, mu[i] = -1;
for (int j = 0; primes[j] * i < N; j ++ )
{
st[primes[j] * i] = 1;
if (i % primes[j] == 0) break;
mu[primes[j] * i] = -mu[i];
}
}
for (int i = 0; i < cnt; i ++ )
for (int j = 1; primes[i] * j < N; j ++ )
f[primes[i] * j] += mu[j];
for (int i = 1; i < N; i ++ ) f[i] += f[i - 1];
}
void solve()
{
int n, m;
cin >> n >> m;
if (n > m) swap(n, m);
ll res = 0;
for (int l = 1, r; l <= n; l = r + 1)
{
r = min(n / (n / l), m / (m / l));
res += 1ll * (f[r] - f[l - 1]) * (n / l) * (m / l);
}
cout << res << "\n";
}
int main()
{
ios::sync_with_stdio(false);
cin.tie(nullptr);
init();
int T;
cin >> T;
while (T -- ) solve();
return 0;
}
