数论函数相关
数论函数相关
数论函数
定义:
- 数论函数:定义域为正整数的函数称为数论函数。
- 加性函数:满足 \(\forall a, b \in \mathbb{N^*}, a \bot b, f(ab) = f(a) + f(b)\) 的函数。
- 积性函数:满足 \(\forall a, b \in \mathbb{N^*}, a \bot b, f(ab) = f(a) \times f(b)\) 的函数。
- 完全积性函数:满足 \(\forall a, b \in \mathbb{N^*}, f(ab) = f(a) \times f(b)\) 的函数。
数论函数的运算:
- 加法:\((f + g)(x) = f(x) + g(x)\) 。
- 数乘:\((c \cdot f)(x) = c \times f(x)\) 。
- 点乘:\((f \cdot g)(x) = f(x) \times g(x)\) 。
常见数论函数(记 \(n\) 的唯一分解为 \(\prod_{i = 1}^m p_i^{c_i}\) ):
-
单位函数: \(\epsilon(n) = [n - 1]\) ,为完全积性函数。
-
常数函数:\(1(n) = 1\) ,为完全积性函数。
-
幂函数:\(\mathrm{id}_k(n) = n^k\) ,为完全积性函数。
-
除数函数:\(\sigma_k(n) = \sum_{d \mid n} d^k\) ,为积性函数。
-
\(\sigma_0(n) = d(n) = \prod_{i = 1}^m (c_i + 1)\) 表示 \(n\) 的约数个数。
-
\(d(ij) = \sum_{x | i} \sum_{y | j} [\gcd(x, y) = 1]\) 。
证明:考虑一个质因子,若都没被枚举则不统计,否则若被 \(x\) 统计,则贡献为 \(x\) 里对应的次幂,否则贡献为 \(i\) 里对应的次幂加上 \(j\) 里对应的次幂。
-
\(d(ijk) = \sum_{a | i} \sum_{b | j} \sum_{c | k} [\gcd(a, b) = 1] [\gcd(b, c) = 1] [\gcd(a, c) = 1]\) 。
-
\(\sigma_0(n^2) = \sum_{d \mid n} \mu(d) \sigma_0(\frac{n}{d})^2 = \sum_{d \mid n} \mu(d)^2 \sigma_0(\frac{n}{d})\) 。
证明:只需证明 \(n = p^k\) 成立即可,原式等价于 \(2k + 1 = (k + 1)^2 - k^2 = k + (k + 1)\) 。
-
-
\(\sigma_1(n) = \sigma(n) = \prod_{i = 1}^m \frac{p_i^{c_i + 1} - 1}{p_i - 1}\) 表示 \(n\) 的约数和,具有如下性质
-
\(\sum_{i = 1}^n \sigma(i) = \sum_{i = 1}^n \lfloor \frac{n}{i} \rfloor\) 。
证明:考虑每个 \(i\) 的贡献即可。
-
-
\(k \ge 1\) 时 \(\sigma_k(n) = \prod_{i = 1}^n \frac{p_i^{(c_i + 1)k} - 1}{p_i - 1}\) ,这是因为 \(n\) 的所有约数的 \(k\) 次方和可以写作 \(\prod_{i = 1}^m \sum_{j = 0}^{c_i} p_i^{jk}\) ,等比数列求和即可。
-
-
本质不同质因子个数函数:\(\omega(n) = \sum_{p \in \mathbb{P}} [p \mid n]\) ,为加性函数。
线性筛
基础的用法是筛素数:
inline void sieve(int n) {
memset(isp, true, sizeof(isp));
isp[1] = false;
for (int i = 2; i <= n; ++i) {
if (isp[i])
pri[++pcnt] = i;
for (int j = 1; j <= pcnt && i * pri[j] <= n; ++j) {
isp[i * pri[j]] = false;
if (!(i % pri[j]))
break;
}
}
}
线性筛可以 \(O(n)\) 求出某些积性函数在 \(1 \sim n\) 处的所有取值,需要满足 \(f(p^c)\) 的值便于分析。
线性筛时若找到一个质数 \(p\) ,则将所有 \(p\) 的幂处的函数取值都求出。记 \(low_i\) 最小质因子的次幂,则对于合数 \(i\) 有 \(f(i) = f(\frac{i}{low_i}) \times f(low_i)\) 。
inline void sieve(int n) {
memset(isp, true, sizeof(isp));
isp[1] = false, f[1] = 1;
for (int i = 2; i <= n; ++i) {
if (isp[i])
pri[++pcnt] = i, low[i] = i, f[i] = ...;
for (int j = 1; j <= pcnt && i * pri[j] <= n; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j]) {
low[i * pri[j]] = low[i] * pri[j];
if (i == low[i])
f[i * pri[j]] = ...;
else
f[i * pri[j]] = f[i / low[i]] * f[low[i * pri[j]]];
} else {
low[i * pri[j]] = low[j], f[i * pri[j]] = f[i] * f[pri[j]];
break;
}
}
}
}
Dirichlet 卷积
定义狄利克雷卷积:
简记为 \(h = f * g\) ,按定义式计算可以做到 \(O(n \ln n)\) 。
由于 \(\epsilon * f = f\) ,因此 \(\epsilon\) 称为狄利克雷卷积的单位元。
定义 \(f\) 的逆元 \(f^{-1}\) 满足 \(f *f^{-1} = \epsilon\) ,其中 \(f\) 可逆当且仅当 \(f(1) \ne 0\) ,可以计算得出:
性质:
- 交换律:\((f * g)(n) = (g * f)(n)\) 。
- 结合律:\(((f * g) * h)(n) = (f * (g * h))(n)\) 。
- 分配律:\(((f + g) * h)(n) = (f * h + g * h)(n)\) 。
- \(f = g \iff f * h = g * h (h(1) \ne 0)\) 。
- 积性函数的狄利克雷卷积是积性函数。
- 积性函数的逆元是积性函数。
常见的狄利克雷卷积式:
- \(\sigma_0 = 1 * 1\) 。
- \(\mathrm{id} = \varphi * 1\) 。
- \(\sigma = 1 * \mathrm{id} = \sigma_0 * \varphi\) 。
- \(1 = \sigma_0 * \mu\) 。
- \(\epsilon = \mu * 1\) 。
- \(\varphi = \mu * \mathrm{id}\) 。
整除分块
注意到 \(\lfloor \frac{n}{x} \rfloor\) 在 \(x \in [1, n]\) 时的取值不超过 \(2 \sqrt{n}\) 种,对于某些问题,对于一段 \(\lfloor \frac{n}{x} \rfloor\) 均相等的区间,可以直接得出该区间里的所有 \(x\) 对答案的贡献。
如果要实现整块一起统计,则需要求出每一块的块头 \(l\) 和块尾 \(r\),得到答案为:
注意到每一块的 \(l\) 都可以由上一块的 \(r\) 推出,故我们继续讨论如何在已知 \(l\) 的情况下推出 \(r\) 。令 \(t = \lfloor \frac{n}{l} \rfloor\) ,容易得到:
直接计算即可,时间复杂度 \(O(\sqrt{n})\) 。
for (int l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
// Calculate the contribution of [l, r]
}
扩展到二维甚至多维也是类似的,时间复杂度 \(O(k \sqrt{n})\) ,其中 \(k\) 为维度。
二维整除分块:
for (int l = 1, r; l <= min(n, m); l = r + 1) {
r = min(n / (n / l), m / (m / l));
// Calculate the contribution of [l, r]
}
一类题目的整除分块形式为向上取整,只需令 \(r \gets \lfloor \frac{n - 1}{\lceil \frac{n}{l} \rceil - 1} \rfloor\) 即可,这是因为:
注意特判 \(\lceil \frac{n}{l} \rceil = 1\) ,此时 \(r\) 的上界为 \(+ \infty\) ,需要取实际上界。
欧拉函数
欧拉函数 \(\varphi(n)\) 定义为 \(1 \sim n\) 内与 \(n\) 互质的数的个数,具有如下性质:
- 欧拉函数是积性函数,特别的有 \(\varphi(2n) = \varphi(n)\) 。
- \(\sum_{d | n} \varphi(d) = n\) 。
- 证明:设 \(f(x)\) 表示 \(\gcd(k, n) = x (k \in [1, n])\) 的个数,则 \(n = \sum_{d | n} f(d) = \sum_{d | n} \varphi(\frac{n}{d}) = \sum_{d | n} \varphi(d)\) 。
- 当 \(n\) 为质数时, \(\varphi(n) = n - 1\) 。
- 当 \(n = p^k (p \in \mathbb{P})\) 时,\(\varphi(n) = p^k - p^{k - 1}\) 。
- 由唯一分解定理,设 \(n = \prod p_i^{k_i}\) ,则 \(\varphi(n) = n \times \prod \frac{p_i - 1}{p_i}\) 。
- $n > 2 \iff 2 | \varphi(n) $ ,证明考虑对称性即可。
- \(\varphi(mn) \varphi(\gcd(m, n)) = \varphi(m) \varphi(n) \gcd(m, n)\) 。
单个数的求解:
inline int euler_phi(int n) {
int phi = n;
for (int i = 2; i * i <= n; ++i)
if (!(n % i)) {
phi = phi / i * (i - 1);
while (!(n % i))
n /= i;
}
if (n > 1)
phi = phi / n * (n - 1);
return phi;
}
线性筛求解:
inline void sieve(int n) {
memset(isp, true, sizeof(isp));
isp[1] = false, phi[1] = 1;
for (int i = 2; i <= n; ++i) {
if (isp[i])
pri[++pcnt] = i, phi[i] = i - 1;
for (int j = 1; j <= pcnt && i * pri[j] <= n; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j])
phi[i * pri[j]] = phi[pri[j]] * phi[i];
else {
phi[i * pri[j]] = pri[j] * phi[i];
break;
}
}
}
}
利用 \(n = \sum_{d | n} \varphi(d)\) 可以对数论式子做一些变换,如:
由此有结论:
GCDMAT - GCD OF MATRIX
给定 \(n, m\) ,\(T\) 次询问,每次给出 \(i_1, j_1, i_2, j_2\) ,求:
\[\sum_{i = i_1}^{i_2} \sum_{j = j_1}^{j_2} \gcd(i, j) \pmod{10^9 + 7} \]\(n, m \le 5 \times 10^4\)
前缀和差分转化为求:
数论分块可以做到 \(O(n) - O(\sqrt{n})\) 。
#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int Mod = 1e9 + 7;
const int N = 5e4 + 7;
int pri[N], phi[N], sum[N];
bool isp[N];
int T, n, m, pcnt;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline void sieve() {
memset(isp, true, sizeof(isp));
isp[1] = false, phi[1] = 1;
for (int i = 2; i < N; ++i) {
if (isp[i])
pri[++pcnt] = i, phi[i] = i - 1;
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
isp[pri[j] * i] = false;
if (i % pri[j])
phi[i * pri[j]] = phi[i] * phi[pri[j]];
else {
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
}
}
for (int i = 1; i < N; ++i)
sum[i] = add(sum[i - 1], phi[i]);
}
inline int S(int n, int m) {
if (n > m)
swap(n, m);
int ans = 0;
for (int l = 1, r = 0; l <= n; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans = add(ans, 1ll * (sum[r] - sum[l - 1]) * (n / l) % Mod * (m / l) % Mod);
}
return ans;
}
signed main() {
sieve();
scanf("%d%d%d", &T, &n, &m);
while (T--) {
int i1, j1, i2, j2;
scanf("%d%d%d%d", &i1, &j1, &i2, &j2);
printf("%d\n", add(dec(dec(S(i2, j2), S(i1 - 1, j2)), S(i2, j1 - 1)), S(i1 - 1, j1 - 1)));
}
return 0;
}
P5221 Product
求:
\[\prod_{i = 1}^n \prod_{j = 1}^n \dfrac{\operatorname{lcm}(i, j)}{\gcd(i, j)} \pmod{104857601} \]\(n \le 10^6\)
单独考虑下面的部分得到 :
单独考虑指数部分得到:
#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int Mod = 104857601;
const int N = 1e6 + 7;
int pri[N], phi[N], sum[N];
bool isp[N];
int n, pcnt;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline void sieve(int n) {
memset(isp, true, sizeof(isp));
isp[1] = false, phi[1] = 1;
for (int i = 2; i <= n; ++i) {
if (isp[i])
pri[++pcnt] = i, phi[i] = i - 1;
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j])
phi[i * pri[j]] = phi[i] * phi[pri[j]];
else {
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
}
}
for (int i = 1; i <= n; ++i)
sum[i] = (sum[i - 1] + phi[i]) % (Mod - 1);
}
signed main() {
scanf("%d", &n);
sieve(n);
int fac = 1;
for (int i = 1; i <= n; ++i)
fac = 1ll * fac * i % Mod;
int ans = 1ll * mi(fac, n * 2) * fac % Mod * fac % Mod;
for (int i = 1; i <= n; ++i)
ans = 1ll * ans * mi(mi(i, 4ll * sum[n / i] % (Mod - 1)), Mod - 2) % Mod;
printf("%d", ans);
return 0;
}
莫比乌斯函数
定义莫比乌斯函数:
性质:
- \(\varphi(n) = \sum_{d | n} \mu(d) \times \frac{n}{d}\) 。
- \(\sum_{d | n} \mu(d) = [n = 1]\) 。
- \(\mu^2(n) = \sum_{d^2 | n} \mu(d)\) 。
- \(d = 1\) 贡献为 \(1\) ,剩下贡献用二项式定理推即可。
单个数的求解:
inline int Mu(int n) {
int mu = 1;
for (int i = 2; i * i <= n; ++i)
if (!(n % i)) {
if (!(n % (i * i)))
return 0;
mu *= -1, n /= i;
}
return n > 1 ? -mu : mu;
}
线性筛求解:
inline void sieve(int n) {
memset(isp, true, sizeof(isp));
isp[1] = false, mu[1] = 1;
for (int i = 2; i <= n; ++i) {
if (isp[i])
pri[++pcnt] = i, mu[i] = -1;
for (int j = 1; j <= pcnt && i * pri[j] <= n; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j])
mu[i * pri[j]] = -mu[i];
else {
mu[i * pri[j]] = 0;
break;
}
}
}
}
事实上可以通过魔力筛在 \(O(n \log \log n)\) 的时间复杂度求出任意数论函数与 \(\mu\) 的狄利克雷卷积。
设 \(g_{i, n} = \sum_{d \mid n} f(d) \mu(\frac{n}{d})\),其中 \(d\) 只含前 \(i\) 种质因子,则:
需要滚动数组优化。
inline void calc(int n) {
memcpy(g + 1, f + 1, sizeof(int) * n);
for (int i = 1, i <= pcnt; ++i)
for (int j = n / pri[i]; j; --j)
g[j * pri[i]] -= g[j];
}
莫比乌斯反演
设 \(f(n), g(n)\) 为两个数论函数,则:
本质就是对质因子做二项式反演。
常见的应用:
推广:
令 \(d k = T\) ,得到:
P2257 YY的GCD / PGCD - Primes in GCD Table / P2568 GCD
求:
\[\sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) \in \mathbb{P}] \]\(n, m \le 10^7\)
推式子得到答案为:
#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1e7 + 7;
ll f[N], sum[N];
int pri[N], mu[N];
bool isp[N];
int T, n, m, pcnt;
inline void sieve() {
memset(isp, true, sizeof(isp));
mu[1] = 1, isp[1] = false;
for (int i = 2; i < N; ++i) {
if (isp[i])
pri[++pcnt] = i, mu[i] = -1;
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
isp[pri[j] * i] = false;
if (i % pri[j])
mu[i * pri[j]] = -mu[i];
else
break;
}
}
for (int i = 1; i <= pcnt; ++i)
for (int j = 1; pri[i] * j < N; ++j)
f[j * pri[i]] += mu[j];
for (int i = 1; i < N; ++i)
sum[i] = sum[i - 1] + f[i];
}
signed main() {
sieve();
scanf("%d", &T);
while (T--) {
scanf("%d%d", &n, &m);
if (n > m)
swap(n, m);
ll ans = 0;
for (int l = 1, r = 0; l <= n; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans += (sum[r] - sum[l - 1]) * (n / l) * (m / l);
}
printf("%lld\n", ans);
}
return 0;
}
P4449 于神之怒加强版
求:
\[\sum_{i = 1}^n \sum_{j = 1}^m \gcd(i, j)^k \pmod{10^9 + 7} \]\(n, m, k \le 5 \times 10^6\)
答案即为:
线性筛 \(g(T) = \sum_{d | T} d^k \mu(\dfrac{T}{d})\) 即可。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 5e6 + 7;
int pri[N], low[N], g[N], sum[N];
bool isp[N];
int T, n, m, k, pcnt;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline void sieve() {
memset(isp, true, sizeof(isp));
isp[1] = false, low[1] = 1, g[1] = 1;
for (int i = 2; i < N; ++i) {
if (isp[i]) {
pri[++pcnt] = i, low[i] = i, g[i] = dec(mi(i, k), 1);
int x = i;
while (1ll * x * i < N)
x *= i, low[x] = x, g[x] = dec(mi(x, k), mi(x / i, k));
}
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
int x = i * pri[j];
isp[x] = false;
low[x] = i % pri[j] ? pri[j] : low[i] * pri[j];
g[x] = 1ll * g[x / low[x]] * g[low[x]] % Mod;
if (!(i % pri[j]))
break;
}
}
for (int i = 1; i < N; ++i)
sum[i] = add(sum[i - 1], g[i]);
}
signed main() {
scanf("%d%d", &T, &k);
sieve();
while (T--) {
scanf("%d%d", &n, &m);
int ans = 0;
if (n > m)
swap(n, m);
for (int l = 1, r; l <= n; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans = add(ans, 1ll * (n / l) * (m / l) % Mod * dec(sum[r], sum[l - 1]) % Mod);
}
printf("%d\n", ans);
}
return 0;
}
P3327 [SDOI2015] 约数个数和
求:
\[\sum_{i = 1}^n \sum_{j = 1}^m d(ij) \]\(T, n, m \le 5 \times 10^4\)
#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 5e4 + 7;
ll f[N];
int pri[N], mu[N], sum[N];
bool isp[N];
int T, n, m, pcnt;
inline void sieve() {
memset(isp, true, sizeof(isp));
isp[1] = false, mu[1] = 1;
for (int i = 2; i < N; ++i) {
if (isp[i])
pri[++pcnt] = i, mu[i] = -1;
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j])
mu[i * pri[j]] = -mu[i];
else {
mu[i * pri[j]] = 0;
break;
}
}
}
for (int i = 1; i < N; ++i)
sum[i] = sum[i - 1] + mu[i];
for (int i = 1; i < N; ++i)
for (int l = 1, r; l <= i; l = r + 1) {
r = i / (i / l);
f[i] += 1ll * (i / l) * (r - l + 1);
}
}
signed main() {
sieve();
scanf("%d", &T);
while (T--) {
scanf("%d%d", &n, &m);
if (n > m)
swap(n, m);
ll ans = 0;
for (int l = 1, r; l <= n; l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans += f[n / l] * f[m / l] * (sum[r] - sum[l - 1]);
}
printf("%lld\n", ans);
}
return 0;
}
P6055 [RC-02] GCD
求:
\[\sum_{i = 1}^{n} \sum_{j = 1}^n \sum_{p = 1}^{\lfloor \frac{n}{j} \rfloor} \sum_{q = 1}^{\lfloor \frac{n}{j} \rfloor} [\gcd(i, j) = 1] [\gcd(p, q) = 1] \pmod{998244353} \]\(n \le 2 \times 10^9\)
杜教筛处理 \(\mu\) 的前缀和即可。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 998244353;
const int N = 1e7 + 7;
map<int, int> mp;
int pri[N], mu[N], sum[N];
bool isp[N];
int n, pcnt;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline void sieve() {
memset(isp, true, sizeof(isp));
isp[1] = false, mu[1] = 1;
for (int i = 2; i < N; ++i) {
if (isp[i])
pri[++pcnt] = i, mu[i] = Mod - 1;
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j])
mu[i * pri[j]] = Mod - mu[i];
else {
mu[i * pri[j]] = 0;
break;
}
}
}
for (int i = 1; i < N; ++i)
sum[i] = add(sum[i - 1], mu[i]);
}
int Sum(int n) {
if (n < N)
return sum[n];
if (mp.find(n) != mp.end())
return mp[n];
int res = 1;
for (int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
res = dec(res, 1ll * (r - l + 1) * Sum(n / l) % Mod);
}
return mp[n] = res;
}
signed main() {
sieve();
scanf("%d", &n);
int ans = 0;
for (int l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
ans = add(ans, 1ll * (n / l) * (n / l) % Mod * (n / l) % Mod * dec(Sum(r), Sum(l - 1)) % Mod);
}
printf("%d", ans);
return 0;
}
P3768 简单的数学题
求:
\[\sum_{i = 1}^n \sum_{j = 1}^n ij \gcd(i, j) \pmod{p} \]\(n \le 10^{10}\)
其中 \(S(n) = \frac{n(n + 1)}{2}\) ,令 \(T = dk\) 则:
杜教筛处理 \(T^2 \varphi(T)\) 的前缀和即可,具体为构造 \(g(n) = n^2\) 。
#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1e7 + 7;
map<ll, int> mp;
int pri[N], phi[N], sum[N];
bool isp[N];
ll n;
int Mod, pcnt, inv6;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline void sieve() {
memset(isp, true, sizeof(isp));
isp[1] = false, phi[1] = 1;
for (int i = 2; i < N; ++i) {
if (isp[i])
pri[++pcnt] = i, phi[i] = i - 1;
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j])
phi[i * pri[j]] = phi[i] * phi[pri[j]];
else {
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
}
}
for (int i = 1; i < N; ++i)
sum[i] = add(sum[i - 1], 1ll * i * i % Mod * phi[i] % Mod);
}
inline int Sum1(ll n) {
n %= Mod;
return n * (n + 1) / 2 % Mod;
}
inline int Sum2(ll n) {
n %= Mod;
return n * (n + 1) % Mod * (n * 2 + 1) % Mod * inv6 % Mod;
}
int Sum(ll n) {
if (n < N)
return sum[n];
if (mp.find(n) != mp.end())
return mp[n];
int res = 1ll * Sum1(n) * Sum1(n) % Mod;
for (ll l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
res = dec(res, 1ll * dec(Sum2(r), Sum2(l - 1)) * Sum(n / l) % Mod);
}
return mp[n] = res;
}
signed main() {
scanf("%d%lld", &Mod, &n);
sieve();
inv6 = mi(6, Mod - 2);
int ans = 0;
for (ll l = 1, r; l <= n; l = r + 1) {
r = n / (n / l);
ans = add(ans, 1ll * Sum1(n / l) * Sum1(n / l) % Mod * dec(Sum(r), Sum(l - 1)) % Mod);
}
printf("%d", ans);
return 0;
}
P3172 [CQOI2015] 选数
求:
\[\sum_{i_1 = L}^R \sum_{i_2 = L}^R \cdots \sum_{i_n = L}^R [\gcd(i_1, i_2, \cdots, i_n) = k] \pmod{10^9 + 7} \]\(n, k \le 10^9\) ,\(L, R \le 10^9\) ,\(R - L \le 10^5\)
令 \(l = \lceil \dfrac{L}{k} \rceil, r = \lfloor \dfrac{R}{k} \rfloor\) ,则:
#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 1e7 + 7;
map<int, int> mp;
int pri[N], mu[N], sum[N];
bool isp[N];
int n, k, L, R, pcnt;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline void sieve() {
memset(isp, true, sizeof(isp));
isp[1] = false, mu[1] = 1;
for (int i = 2; i < N; ++i) {
if (isp[i])
pri[++pcnt] = i, mu[i] = Mod - 1;
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j])
mu[i * pri[j]] = Mod - mu[i];
else {
mu[i * pri[j]] = 0;
break;
}
}
}
for (int i = 1; i < N; ++i)
sum[i] = add(sum[i - 1], mu[i]);
}
inline int Sum(int n) {
if (n < N)
return sum[n];
if (mp.find(n) != mp.end())
return mp[n];
int res = 1;
for (int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
res = dec(res, 1ll * (r - l + 1) * Sum(n / l) % Mod);
}
return mp[n] = res;
}
signed main() {
sieve();
scanf("%d%d%d%d", &n, &k, &L, &R);
L = (L - 1) / k, R /= k;
int ans = 0;
for (int l = 1, r; l <= R; l = r + 1) {
r = R / (R / l);
if (L / l)
r = min(r, L / (L / l));
ans = add(ans, 1ll * mi(R / l - L / l, n) * dec(Sum(r), Sum(l - 1)) % Mod);
}
printf("%d", ans);
return 0;
}
P4619 [SDOI2018] 旧试题
给定 \(A, B, C\) ,求:
\[\sum_{i = 1}^A \sum_{j = 1}^B \sum_{k = 1}^C d(ijk) \]\(A, B, C \le 2 \times 10^5\)
首先套路地推莫反:
记 \(f(n) = \sum_{i = 1}^n \lfloor \frac{n}{i} \rfloor\) ,这可以 \(O(V \sqrt{V})\) 数论分块预处理,继续推式子得到:
事实上很多 \((x, y, z)\) 的贡献都是 \(0\) ,有贡献的三元组一定满足:
- \(\mu(x) \ne 0\) ,\(\mu(y) \ne 0\) ,\(\mu(z) \ne 0\) 。
- \(\max \{ \mathrm{lcm}(x, y), \mathrm{lcm}(x, z), \mathrm{lcm}(y, z) \} \le \max \{ A, B, C \}\) 。
考虑将可能有贡献的点之间连边,枚举 \(\gcd = i\) ,然后拿出 \(i\) 的倍数 \(j, k\) 满足 \(\mu(j) \ne 0\) 且 \(\mu(k) \ne 0\) 连边。问题转化为三元环计数,可能需要额外处理一下有数字相等的情况。
事实上这张图的边数是很少的 \(m = 821535\) ,因此直接跑三元环计数即可做到 \(O(m \sqrt{m})\) 。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 1e9 + 7;
const int N = 1e5 + 7, S = 1e6 + 7;
vector<int> g[N];
struct Edge {
int u, v, w;
} e[S];
int pri[N], mu[N], deg[N], f[N], tag[N];
bool isp[N];
int A, B, C, n, tot, pcnt;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline void sieve() {
memset(isp, true, sizeof(isp));
isp[1] = false, mu[1] = 1;
for (int i = 2; i < N; ++i) {
if (isp[i])
pri[++pcnt] = i, mu[i] = -1;
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j])
mu[i * pri[j]] = -mu[i];
else {
mu[i * pri[j]] = 0;
break;
}
}
}
}
inline int solve() {
for (int i = 1; i <= n; ++i)
g[i].clear();
for (int i = 1; i <= tot; ++i) {
int &u = e[i].u, &v = e[i].v;
if (!(deg[u] == deg[v] ? u < v : deg[u] < deg[v]))
swap(u, v);
g[u].emplace_back(i);
}
int ans = 0;
for (int i = 1; i <= n; ++i) {
for (int j : g[i])
tag[e[j].v] = e[j].w;
for (int j : g[i])
for (int k : g[e[j].v]) {
if (!tag[e[k].v])
continue;
int u = i, v = e[j].v, w = e[k].v, lcma = tag[e[k].v], lcmb = e[j].w, lcmc = e[k].w, res = 0;
auto calc = [&](int a, int b, int c) {
return 1ll * f[a / lcma] * f[b / lcmb] % Mod * f[c / lcmc] % Mod;
};
if (u == v && v == w)
res = add(res, calc(A, B, C));
else if (u == v || u == w || v == w)
res = add(res, add(calc(A, B, C), add(calc(B, C, A), calc(C, A, B))));
else {
res = add(res, add(calc(A, B, C), calc(A, C, B)));
res = add(res, add(calc(B, A, C), calc(B, C, A)));
res = add(res, add(calc(C, A, B), calc(C, B, A)));
}
if (mu[u] * mu[v] * mu[w] == 1)
ans = add(ans, res);
else
ans = dec(ans, res);
}
for (int j : g[i])
tag[e[j].v] = 0;
}
return ans;
}
signed main() {
sieve();
for (int i = 1; i < N; ++i)
for (int l = 1, r; l <= i; l = r + 1)
r = i / (i / l), f[i] = add(f[i], 1ll * (r - l + 1) * (i / l) % Mod);
int T;
scanf("%d", &T);
while (T--) {
scanf("%d%d%d", &A, &B, &C);
memset(deg + 1, 0, sizeof(int) * (n = max(max(A, B), C))), tot = 0;
for (int i = 1; i <= n; ++i)
for (int j = i; j <= n; j += i) {
if (!mu[j])
continue;
for (int k = j; 1ll * j / i * k <= n; k += i)
if (mu[k] && __gcd(j, k) == i)
e[++tot] = (Edge){j, k, j / i * k}, ++deg[j], ++deg[k];
}
printf("%d\n", solve());
}
return 0;
}
单位根反演
公式:
证明可以等比数列求和,也可以用向量加法。
推论:
LOJ6485. LJJ 学二项式定理
给定 \(n, s, a_{0 \sim 3}\) ,求:
\[\left( \sum_{i = 0}^n \binom{n}{i} \times s^i \times a_{i \bmod 4} \right) \pmod{998244353} \]\(n \le 10^{18}\)
时间复杂度 \(O(\log n)\) 。
#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int Mod = 998244353, buf = 911660635, inv4 = 1ll * ((Mod + 1) / 2) * ((Mod + 2) / 2) % Mod;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
signed main() {
int T;
scanf("%d", &T);
while (T--) {
ll n;
int s, a[4], ans = 0;
scanf("%lld%d%d%d%d%d", &n, &s, a + 0, a + 1, a + 2, a + 3);
for (int i = 0; i <= 3; ++i)
for (int j = 0, pw = 1; j <= 3; ++j, pw = 1ll * pw * buf % Mod)
ans = add(ans, 1ll * a[i] * mi(buf, Mod - 1 - i * j) % Mod *
mi(add(1ll * s * pw % Mod, 1), n % (Mod - 1)) % Mod);
printf("%d\n", 1ll * inv4 * ans % Mod);
}
return 0;
}
P5591 小猪佩奇学数学
给定 \(n, p, k\) ,求:
\[\sum_{i = 0}^n \binom{n}{i} \times p^i \times \lfloor \frac{i}{K} \rfloor \pmod{998244353} \]\(1 \le n, p < 998244353\) ,\(k \in \{ 2^w \mid 0 \le w \le 20 \}\)
用 \(\lfloor \frac{i}{k} \rfloor = \frac{i - i \bmod k}{k}\) 展开得到:
由吸收恒等式得到:
继续考虑后面一项:
考虑后面的 \(\sum_{j = 0}^{k - 1} j \omega_n^{-t j}\) ,设 \(S(n, k) = \sum_{i = 0}^{n - 1} i k^i\) ,则答案为:
下面考虑算 \(S\) 。当 \(k = 1\) 时有 \(S(n, 1) = \sum_{i = 0}^n i = \frac{(n - 1)n}{2}\) ,当 \(k \not = 1\) 时有:
可得 \(S(n, k) = \frac{(n - 1) k^n + \frac{1 - k^n}{k - 1} + 1}{k - 1}\) 。注意到带入的 \(k\) 总是 \(n\) 次单位根,故 \(k^n = 1\) ,则 \(S(n, k) = \frac{n}{k - 1}\) 。
时间复杂度 \(O(k \log n)\) 。
#include <bits/stdc++.h>
using namespace std;
const int Mod = 998244353, rt = 3;
int n, p, k;
inline int add(int x, int y) {
x += y;
if (x >= Mod)
x -= Mod;
return x;
}
inline int dec(int x, int y) {
x -= y;
if (x < 0)
x += Mod;
return x;
}
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline int S(int n, int k) {
return k == 1 ? (1ll * n * (n - 1) / 2) % Mod : 1ll * n * mi(k - 1, Mod - 2) % Mod;
}
signed main() {
scanf("%d%d%d", &n, &p, &k);
vector<int> pw(k + 1);
pw[0] = 1, pw[1] = mi(rt, (Mod - 1) / k);
for (int i = 2; i <= k; ++i)
pw[i] = 1ll * pw[i - 1] * pw[1] % Mod;
int ans = 0;
for (int i = 0; i < k; ++i)
ans = add(ans, 1ll * mi(add(1ll * p * pw[i] % Mod, 1), n) * S(k, pw[k - i]) % Mod);
printf("%d", 1ll * mi(k, Mod - 2) *
dec(1ll * n * p % Mod * mi(p + 1, n - 1) % Mod, 1ll * ans * mi(k, Mod - 2) % Mod) % Mod);
return 0;
}
杜教筛
杜教筛被用于处理一类数论函数的前缀和问题,可以在低于线性时间复杂度内计算某些数论函数的前缀和。
考虑构造一个 \(S(n)\) 关于 \(S(\lfloor \frac{n}{i} \rfloor)\) 的递推式。对于任意数论函数 \(g\) ,必有:
那么可以得到递推式:
假如可以构造恰当的数论函数 \(g\) 使得:
- 可以快速计算 \(\sum_{i = 1}^n (f * g)(i)\) 。
- 可以快速计算 \(g\) 的单点值,以便于用整除分块求解 \(\sum_{i = 2}^n g(i) S(\lfloor \frac{n}{i} \rfloor)\) 。
则可以用杜教筛求出 \(S(n)\) 。预处理 \(S(1 \sim n^{\frac{2}{3}})\) 的值即可做到 \(O(n^{\frac{2}{3}})\) 的复杂度,注意需要记忆化保证复杂度正确性。
P4213【模板】杜教筛(Sum)
求欧拉函数、莫比乌斯函数的前 \(n\) 项和。
\(n < 2^{31}\)
对于欧拉函数前缀和的求解,由 \(\varphi * 1 = id\) ,从而:
对于莫比乌斯函数前缀和的求解,由 \(\mu * 1 = \epsilon\) ,从而:
#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1e7 + 7;
int pri[N];
bool isp[N];
int pcnt;
namespace Phi {
map<ll, ll> mp;
ll s[N];
int f[N];
inline void prework() {
for (int i = 1; i < N; ++i)
s[i] = s[i - 1] + f[i];
}
inline ll Sum(ll n) {
if (n < N)
return s[n];
if (mp.find(n) != mp.end())
return mp[n];
ll res = 1ll * n * (n + 1) / 2;
for (ll l = 2, r; l <= n; l = r + 1)
r = n / (n / l), res -= 1ll * (r - l + 1) * Sum(n / l);
return mp[n] = res;
}
} // namespace Phi
namespace Mu {
map<ll, ll> mp;
ll s[N];
int f[N];
inline void prework() {
for (int i = 1; i < N; ++i)
s[i] = s[i - 1] + f[i];
}
inline ll Sum(ll n) {
if (n < N)
return s[n];
if (mp.find(n) != mp.end())
return mp[n];
ll res = 1;
for (ll l = 2, r; l <= n; l = r + 1)
r = n / (n / l), res -= 1ll * (r - l + 1) * Sum(n / l);
return mp[n] = res;
}
} // namespace Mu
inline void sieve() {
memset(isp, true, sizeof(isp));
isp[1] = false, Phi::f[1] = Mu::f[1] = 1;
for (int i = 2; i < N; ++i) {
if (isp[i])
pri[++pcnt] = i, Phi::f[i] = i - 1, Mu::f[i] = -1;
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j])
Phi::f[i * pri[j]] = Phi::f[i] * (pri[j] - 1), Mu::f[i * pri[j]] = -Mu::f[i];
else {
Phi::f[i * pri[j]] = Phi::f[i] * pri[j], Mu::f[i * pri[j]] = 0;
break;
}
}
}
Phi::prework(), Mu::prework();
}
signed main() {
sieve();
int T;
scanf("%d", &T);
while (T--) {
ll n;
scanf("%lld", &n);
printf("%lld %lld\n", Phi::Sum(n), Mu::Sum(n));
}
return 0;
}
P3768 简单的数学题
给定 \(n, p\) ,求:
\[\left( \sum_{i = 1}^n \sum_{j = 1}^n ij \gcd(i, j) \right) \bmod p \]\(n \le 10^{10}\) ,\(5 \times 10^8 \le p \le 1.1 \times 10^9\) 且 \(p\) 为质数
考虑将 \(\gcd\) 用 \(\varphi\) 反演掉:
直接杜教筛 \(\varphi(d) d^2\) ,为了消掉 \(d^2\) 取 \(g = id^2\) 即可。
#include <bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1e7 + 7;
map<ll, ll> mp;
ll sum[N];
int pri[N], phi[N];
bool isp[N];
ll n, Mod, inv6;
int pcnt;
inline int mi(int a, int b) {
int res = 1;
for (; b; b >>= 1, a = 1ll * a * a % Mod)
if (b & 1)
res = 1ll * res * a % Mod;
return res;
}
inline void sieve() {
memset(isp, true, sizeof(isp));
isp[1] = false, phi[1] = 1;
for (int i = 2; i < N; ++i) {
if (isp[i])
pri[++pcnt] = i, phi[i] = i - 1;
for (int j = 1; j <= pcnt && i * pri[j] < N; ++j) {
isp[i * pri[j]] = false;
if (i % pri[j])
phi[i * pri[j]] = phi[i] * (pri[j] - 1);
else {
phi[i * pri[j]] = phi[i] * pri[j];
break;
}
}
}
for (int i = 1; i < N; ++i)
sum[i] = (sum[i - 1] + 1ll * phi[i] * i % Mod * i % Mod) % Mod;
}
inline ll S1(ll n) {
n %= Mod;
return n * (n + 1) / 2 % Mod;
}
inline ll S2(ll n) {
n %= Mod;
return n * (n + 1) % Mod * (n * 2 + 1) % Mod * inv6 % Mod;
}
inline ll S3(ll n) {
return S1(n) * S1(n) % Mod;
}
ll Sum(ll n) {
if (n < N)
return sum[n];
if (mp.find(n) != mp.end())
return mp[n];
ll res = S3(n);
for (ll l = 2, r; l <= n; l = r + 1)
r = n / (n / l), res = (res - 1ll * (S2(r) - S2(l - 1) + Mod) % Mod * Sum(n / l) % Mod + Mod) % Mod;
return mp[n] = res;
}
signed main() {
scanf("%lld%lld", &Mod, &n);
sieve(), inv6 = mi(6, Mod - 2);
ll ans = 0;
for (ll l = 1, r; l <= n; l = r + 1)
r = n / (n / l), ans = (ans + 1ll * (Sum(r) - Sum(l - 1) + Mod) % Mod * S3(n / l) % Mod) % Mod;
printf("%lld", ans);
return 0;
}