2024.12.11 数论 1 讲课笔记
记号说明
\(f(n)=O(n^\varepsilon)\iff \forall \varepsilon>0,f(n)=O(n^\varepsilon)\)
\(\gamma\)：欧拉常数
\(f\ast g\)：狄利克雷卷积，有交换律和结合律
\(\varepsilon(n)\)：单位函数
\(\operatorname{1}(n)\)：常数函数
\(\operatorname{id}(n)\)：恒等函数
\(a//b\)：\(a\) 除以 \(b\) 结果下取整
一些式子的渐进估计
\(d(n)=O(n^{\varepsilon})\)
\(\sum_{i=1}^n\frac 1i=\ln n+\gamma+O(n^{-1})\)
欧几里得辗转相除法
int gcd(int a, int b){return b? gcd(b, a % b) : a;}
扩展欧几里得算法
递归：
int exgcd(int a, int b, int &x, int &y){
if (b == 0){x = 1, y = 0;return a;}
int res = exgcd(b, a % b, x, y);
int X = x;x = y;y = X - a / b * y;return res;
}
非递归：
int exgcd(int a, int b, int &x1, int &x2){
x1 = 1, x2 = 0; int x3 = 0, x4 = 1;
while (b){int c = a / b;tie(x1, x2, x3, x4, a, b) = make_tuple(x3, x4, x1 - x3 * c, x2 - x4 * c, b, a - b * c);}
return a;
}
质数判断
朴素实现
bool is_prime(int v){
if (v < 2)return 0;
for (int i = 2; i * i <= v; ++i)if (v % i == 0)return 0;
return 1;
}
二次探测定理
Rabin-Miller 算法
模板：
const int prl[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
long long mul(long long a, long long b, long long M){return __int128(a) * b % M;}
long long gcd(long long a, long long b){return __gcd(a, b);}
long long add(long long a, long long b, long long M){a += b; while (a >= M)a -= M; return a;}
long long fpw(long long a, long long b, long long M){long long ret = 1;for (; b; b >>= 1, a = mul(a, a, M))if (b & 1)ret = mul(ret, a, M);return ret;}
bool check(int a, long long v){
long long d = v - 1, G = fpw(a, d, v);if (G != 1)return 1;
while ((d & 1) == 0){d >>= 1;if ((G = fpw(a, d, v)) == v - 1)return 0;else if (G != 1)return 1;}
return 0;
}
bool is_prime(long long v){
if (v <= 40){for (int a : prl)if (v == a)return 1;return 0;}
if (v < (1ll << 32)){if (check(2, v) || check(7, v) || check(61, v))return 0;return 1;}
for (int a : prl)if (check(a, v))return 0;return 1;
}
分解质因数
Pollard-Rho 算法
模板：
const int prl[] = {2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37};
long long mul(long long a, long long b, long long M){return __int128(a) * b % M;}
long long gcd(long long a, long long b){return __gcd(a, b);}
long long add(long long a, long long b, long long M){a += b; while (a >= M)a -= M; return a;}
long long fpw(long long a, long long b, long long M){long long ret = 1;for (; b; b >>= 1, a = mul(a, a, M))if (b & 1)ret = mul(ret, a, M);return ret;}
bool check(int a, long long v){
long long d = v - 1, G = fpw(a, d, v);if (G != 1)return 1;
while ((d & 1) == 0){d >>= 1;if ((G = fpw(a, d, v)) == v - 1)return 0;else if (G != 1)return 1;}
return 0;
}
bool prm(long long v){
if (v <= 40){for (int a : prl)if (v == a)return 1;return 0;}
if (v < (1ll << 32)){if (check(2, v) || check(7, v) || check(61, v))return 0;return 1;}
for (int a : prl)if (check(a, v))return 0;return 1;
}
long long f(long long x, long long c, long long M){return add(mul(x, x, M), c, M);}
mt19937_64 rnd(time(nullptr));
long long PR(long long v){
long long c = rnd() % (v - 1) + 1;
long long t = f(0, c, v), r = f(f(0, c, v), c, v);
while (t != r){
long long d = gcd(abs(t - r), v);if (d > 1)return d;
t = f(t, c, v), r = f(f(r, c, v), c, v);
}
return v;
}
long long mxv;
void dfs(long long v){
if (v < mxv)return;
if (v < 100){
for (long long i = 2; i * i <= v; ++i)
if (v % i == 0){mxv = max(mxv, i);while (v % i == 0)v /= i;}
if (v != 1)mxv = max(mxv, v);
return;
}
long long d;while ((d = PR(v)) == v);if (prm(d))mxv = max(mxv, d); else dfs(d);
d = v / d;if (prm(d))mxv = max(mxv, d); else dfs(d);
}
long long get(long long v){mxv = 0; dfs(v); return mxv;}
void work(){
long long n, v;cin >> n;
if (prm(n))cout << "Prime" << endl;
else cout << get(n) << endl;//输出最小质因子
}
扩展中国剩余定理
模板：
#define int long long
int n, m[100010], a[100010];
bool res = 1;
void mul(int &a, int b, int c){
int ret = 0;
if (b < 0)b = -b, a = -a;
for (; b; b >>= 1, a = (a + a) % c)
if (b & 1)ret = (ret + a) % c;
a = ret;
}
int _exgcd(int a, int b, int &x, int &y){
if (b == 0){x = 1, y = 0;return a;}
int r = _exgcd(b, a % b, x, y);y = x - a / b * (x = y);return r;
}
int exgcd(int a, int b, int c, int &x){int y, t = _exgcd(a, b, x, y);mul(x, c / t, b / t);return t;}
void un(int a1, int a2, int m1, int m2, int &A, int &M){
int k1, t = exgcd(m1, m2, a1 - a2, k1);
M = m1 / t * m2;A = (a1 % M - k1 * m1 % M + M) % M;
}
int excrt(){
int A = a[1], M = m[1];
for (int i = 2; i <= n; ++i)un(A, a[i], M, m[i], A, M);
return A;
}
Lucas 定理
模板：
int fc[200010], inv[200010];
int C(int n, int m, int p){
if (m < n)return 0;
if (n >= p || m >= p)return 1ll * C(n / p, m / p, p) * C(n % p, m % p, p) % p;
return 1ll * fc[m] * inv[n] % p * inv[m - n] % p;
}
int fp(int a, int b, int p){
int ret = 1;
for (; b; b >>= 1, a = 1ll * a * a % p)
if (b & 1)ret = 1ll * ret * a % p;
return ret;
}
//预处理
fc[0] = 1;for (int i = 1; i < p; ++i)fc[i] = 1ll * fc[i - 1] * i % p;
inv[p - 1] = fp(fc[p - 1], p - 2, p);for (int i = p - 2; i >= 0; --i)inv[i] = 1ll * inv[i + 1] * (i + 1) % p;
扩展 Lucas
模板：
#include <bits/stdc++.h>
#define int long long
using namespace std;
int mul(int a, int b, int M){return __int128(a) * b % M;}
int mulc(int M, int a, int b, int c){return mul(mul(a, b, M), c, M);}
int fpw(int a, int b){int ret = 1;for (; b; b >>= 1, a *= a)if (b & 1)ret *= a;return ret;}
int fpw(int a, int b, int m){int ret = 1;for (; b; b >>= 1, a = mul(a, a, m))if (b & 1)ret = mul(ret, a, m);return ret;}
int fct(int n, int p, int M){//n!/p^x mod M
vector<int> f(p); f[0] = 1;
for (int i = 1; i < p; ++i)f[i] = f[i - 1] * i % M;
int res = 1, tg = 0;
while (n > 1){tg ^= (n / p) & 1;res = res * f[n % p] % M;n /= p;}
if (tg)res = M - res;return res;
}
int G(int n, int p){return n < p? 0 : G(n / p, p) + n / p;}
void exgcd(int a, int b, int &x, int &y){
if (!b)return x = 1, y = 0, (void)0;
exgcd(b, a % b, x, y); y = x - a / b * (x = y);
}
int inv(int a, int p){int x, y;exgcd(a, p, x, y);return (x + p) % p;}
pair<int, int> C(int n, int m, int a, int r){//C(n,m)%a^r, a^r
int p = fpw(a, r);
int N = fct(n, a, p), M = fct(m, a, p), Nm = fct(n - m, a, p);
return {mulc(p, fpw(a, G(n, a) - G(m, a) - G(n - m, a), p), N, inv(M * Nm % p, p)), p};
}
int crt(vector<pair<int, int> > ve, int n){
int ans = 0;
for (int i = 0; i < ve.size(); ++i){
auto [a, r] = ve[i];int m = n / r, b, y;
exgcd(m, r, b, y);ans = (ans + a * m * b % n) % n;
}
return (ans % n + n) % n;
}
signed main(){
int n, m, p, t;
cin >> n >> m >> p, t = p;
vector<pair<int, int> > fc;
for (int i = 2; i * i <= p; ++i)
if (p % i == 0){
int cnt = 0;while (p % i == 0)p /= i, ++cnt;
fc.emplace_back(i, cnt);
}
if (p > 1)fc.emplace_back(p, 1);
vector<pair<int, int> > crtv;
for (auto [a, r] : fc)
crtv.emplace_back(C(n, m, a, r));
cout << crt(crtv, t) << endl;
return 0;
}
阶 和 原根
阶
若 \(b\) 为满足 \(a^b\equiv 1\pmod p\) 的最小值，则 \(b\) 为 \(a\) 模 \(p\) 的阶，记为 \(\operatorname{ord}_p(a)\) 或 \(\delta_p(a)\)，阶一定存在
性质 1
\([a^i\bmod p\mid 1\le i\le\delta_p(a)]\) 两两不同
性质 2
若 \(a^n\equiv 1\pmod p\)，则 \(\delta_p(a)\mid p\)
推论
若 \(a^x\equiv a^y\pmod p\)，则 \(x\equiv y\pmod{\delta_p(a)}\)
性质 3
若 \(p\in \Bbb N^+,\,a,b\in\Bbb Z,\,a\bot p,\,b\bot p\)，则 \(\delta_p(ab)=\delta_p(a)\delta_p(b)\iff \delta_p(a)\bot\delta_p(b)\)
性质 4
若 \(k\in\Bbb N,p\in\Bbb N^+,a\in\Bbb Z,a\bot p\)，则 \(\delta_p(a^k)=\dfrac{\delta_p(a)}{\gcd(\delta_p(a),k)}\)
原根定义
若 \(a\bot p\)，且 \(\delta_p(a)=\varphi(p)\)，则 \(a\) 为模 \(p\) 意义下的原根
原根判定定理
若 \(a\bot p\)，则 \(a\) 为模 \(p\) 意义下的原根当且仅当 \(\forall (q\mid \varphi(p),q\in\{Prime\}),a^{\large\frac{\varphi(p)}q}\not\equiv 1\pmod p\)
原根存在定理
\(p\) 存在原根当且仅当 \(p\in\{2,4,q^k,2q^k\}\;(q\in\{Prime\},q\ne 2,k\in\Bbb N^+)\)
原根个数
若 \(p\) 有原根，则其原根数量为 \(\varphi(\varphi(p))\)
最小原根的范围估计
素数 \(p\) 的最小原根 \(g_p=O(p^{0.25+\varepsilon})\;(\varepsilon>0)\)，且 \(g_p=\Omega(\log p)\)
找出一个数的所有原根
先判断是否存在，再暴力找到最小的原根 \(g_p\)，然后暴力枚举 \(k\)，若 \(\varphi(p)\bot k\)，则 \(g_p^k\) 为原根
模板题：P6091 【模板】原根
模板：
int fpw(int a, int b, int M){int ret = 1;for (; b; b >>= 1, a = 1ll * a * a % M)if (b & 1)ret = 1ll * ret * a % M;return ret;}
int phi(int a){
int ret = a;
for (int i = 2; i * i <= a; ++i)
if (a % i == 0){while (a % i == 0)a /= i;ret = ret / i * (i - 1); }
if (a > 1)ret = ret / a * (a - 1);return ret;
}
void sol(vector<int> &ret, int n){//ret 保存了 n 的所有原根
if (n == 2){ret = {1};return ;}
int ph = phi(n), t = ph;
vector<int> P;
for (int i = 2; i * i <= ph; ++i)
if (ph % i == 0){
P.emplace_back(t / i);
while (ph % i == 0)ph /= i;
}
if (ph > 1)P.emplace_back(t / ph);
int G = 0;
for (int g = 1; g < n; ++g){
if (__gcd(g, n) != 1)continue;
bool fl = 1;
for (int p : P)
if (fpw(g, p, n) == 1){fl = 0;break;}
if (fl){G = g;break;}
}
if (G)
for (int i = 1, T = 1; i < t; ++i){
T = 1ll * T * G % n;
if (__gcd(i, t) == 1)ret.emplace_back(T);
}
sort(ret.begin(), ret.end());
}
例：[ABC212G] Power Pair
给定 \(p\)，求 \(\sum_{x=0}^{p-1}\sum_{y=0}^{p-1}[\exists n\in\Bbb N,x^n\equiv y\pmod p]\) 取模，\(p\le10^{12},p\in\{Prime\}\)
若 \(x=0\)，则 \(y=0\)；反之若 \(y=0\)，则 \(x=0\)
因此答案为 \(1+\sum_{x=1}^{p-1}\sum_{y=1}^{p-1}[\exists n\in\Bbb N,x^n\equiv y\pmod p]\)
令 \(p\) 的某一原根为 \(g\)，则 \(x\) 可以表示为 \(g^a\)，\(y\) 可以表示为 \(g^b\)（其中 \(1\le a,b\le p-1\)），则 \(x^n\equiv y\pmod p\iff na\equiv b\pmod{(p-1)}\)
根据裴蜀定理，\([\exists n\in\Bbb N,na\equiv b\pmod{(p-1)}]=[\gcd(p-1,a)\mid b]\)
因此答案为：
直接计算即可
时间复杂度为 \(O(n^{\varepsilon}\sqrt n)=O(n^{0.5+\varepsilon})\)
离散对数
给定 \(a,b,p\)，求出最小的 \(x\)，使 \(a^x\equiv b\pmod p,p\le10^9\)
a 和 p 互质的情况
令 \(s=\lceil\sqrt p\rceil,x=rs-c\)
则 \(a^{rs-c}\equiv b\pmod p\)
即 \((a^s)^r\equiv b^c\pmod p\)
将所有 \(b^c\bmod p\;(0\le c<s)\) 保存到 map 或 unordered_map 中，枚举 \(r\) 查询即可
时间复杂度 \(O(s\log s)=O(\sqrt p\log p)\) 或 \(O(s)=O(\sqrt p)\)
一般情况
方法 1
令 \(d_1=\gcd(a,p)\)，若 \(d_1\nmid b\) 则无解，否则方程变为 \(\frac{a}{d_1}a^{x-1}\equiv\frac b{d_1}\pmod {\frac p{d_1}}\)
重复上述操作直到 \(a\) 和 \(\frac pd\) 互质（其中 \(d\) 为每次的 \(d\) 的积）
此时方程变为 \(\frac{a^k}da^{x-k}\equiv \frac bd\pmod {\frac p{d}}\)
转化为 \(a^{x-k}\equiv \frac bd{\left(\frac{a^k}d\right)}^{-1}\pmod {\frac p{d}}\)，用 \(a\) 和 \(p\) 互质时的方法求解即可
需要特判 \(x\le k\) 的情况
方法 2
可证答案上界为 \(2\varphi(p)\)，因此类似一般情况根号平衡，令 \(t=\left\lceil\sqrt{2\varphi(p)}\right\rceil\)
将 \(a^{pt}\bmod p\;(0\le p\le t)\) 存入哈希表，值相同的保留 \(p\) 最小的两个，枚举 \(ba^q\) 并检验两值，若合法则 \(pt-q\) 为一合法解
可证这样一定不会漏最小解
时间复杂度 \(O(\sqrt p)\)
二次剩余
即模意义下的开根
定义
对于 \(a,p\)，若存在 \(x\) 使得 \(x^2\equiv a\pmod p\)，则 \(a\) 为模 \(p\) 下的二次剩余，否则为非二次剩余
数量
假定 \(p\) 为奇质数
则二次剩余有 \(\frac {p+1}2\) 个，非二次剩余有 \(\frac{p-1}2\) 个
勒让德记号
定义勒让德记号 \(\left(\frac ap\right)\)，当 \(p\mid a\) 时值为 \(0\)，若 \(a\) 为模 \(p\) 的二次剩余则值为 \(1\)，否则为 \(-1\)
根据欧拉判别准则，有 \(\left(\frac ap\right)\equiv a^{\frac{p-1}2}\bmod p\)
Cipolla 算法
求出一个 \(x\) 使得 \(x^2\equiv a\pmod p\)
首先判断是否是二次剩余
然后找到一个 \(w\) 满足 \(w^2-a\) 为模 \(p\) 下的非二次剩余（每次随机选一个，约 \(\frac 12\) 的概率非二次剩余，期望选 \(O(1)\) 次）
定义虚数单位 \(i\) 为满足 \(i^2\equiv w^2-a\pmod p\) 的数（类比负数开根得到的虚数），则 \((w+i)^{\frac{p-1}2}\bmod p\) 即为 \(x^2\equiv a\pmod p\) 的一个解（计算过程类似复数运算，实部和虚部分别对 \(p\) 取模，可证结果的虚部一定为 \(0\)）
莫比乌斯反演
莫比乌斯函数
\(\mu(n)=\begin{cases} 1&n=1\\ 0&\exists k>1,k^2\mid n\\ (-1)^{x}&n\ne 1,\nexists k>1,k^2\mid n,n=\prod_{i=1}^x p_i,\forall p_i,p_i\in\{Prime\} \end{cases}\)
显然为积性函数
一种简单的 \(O(n\log n)\) 写法：
mu[1] = 1;
for (int i = 1; i <= n; ++i)
for (int j = i + i; j <= n; j += i)mu[j] -= mu[i];
性质
即在狄利克雷卷积中，莫比乌斯函数为常数函数的逆元
常用扩展：
常用公式（本质为莫反的特殊形式）
莫比乌斯变换
形式 \(1\)：
另一种写法：\(f=\operatorname{1}\ast \,g\iff g=\mu\ast f\)（左式两边同时卷 \(\mu\) 即可得右式）
这种形式下，\(f\) 为 \(g\) 的莫比乌斯变换，\(g\) 为 \(n\) 的莫比乌斯逆变换（反演）
形式 \(2\)：
扩展
设 \(t\) 为完全积性函数，则
例 1：P2522 [HAOI2011] Problem b
\(t\) 组询问，给定 \(a,b,c,d,k\)，查询 \(\sum_{x=a}^b\sum_{y=c}^d[\gcd(x,y)=k]\)，\(t,a,b,c,d,k\le5\times10^4\)
先拆为四个
这等于
令 \(N=\lfloor\frac nk\rfloor,M=\lfloor\frac mk\rfloor\)
则上式转化为
显然可以数论分块
时间复杂度 \(O\left(\sum\sqrt{\normalsize\frac{\min(n,m)}k}\right)\)
例 2：P1829 [国家集训队] Crash的数字表格 / JZPTAB
给定 \(n,m\)，求 \(\sum_{i=1}^n\sum_{j=1}^m\operatorname{lcm}(n,m)\)，\(n,m\le10^7\)
化式子：
令 \(f(x)=\frac{x(x+1)}2\)，\(F(n,m)=\sum_{D=1}^{\min(n,m)}\mu(D)D^2 f(n//D)f(m//D)\)
则上式等于 \(\sum_{d=1}^{\min(n,m)}d F(\lfloor\frac nd\rfloor,\lfloor\frac md\rfloor)\)
显然可以预处理加数论分块套数论分块解决
时间复杂度 \(O(n+m)\)
例 3：P3327 [SDOI2015] 约数个数和
给定 \(T\) 组 \(n,m\)，求出 \(\sum_{i=1}^n\sum_{j=1}^m d(ij)\)，其中 \(d(x)\) 为 \(x\) 的因子数量，\(T,n,m\le5\times10^4\)
可证 \(d(ij)=\sum_{x\mid i}\sum_{y\mid j}\varepsilon(\gcd(x,y))\)
因此
预处理出所有 \(s_x=\sum_{i=1}^x(x//i)\)，则上式容易数论分块解决
时间复杂度 \(O(\max n+\sum\sqrt n)\)（假设 \(n,m\) 同阶）
例 4：CF2039F2 Shohag Loves Counting (Hard Version)
对于 \(a_{1\sim n}\)，定义 \(f(k)=\gcd_{i=1}^{n=k+1}(\sum_{j=i}^{j+k-1}a_j)\)。定义 \(a_{1\sim n}\) 是好的当且仅当 \(\forall 1\le i<j\le n,f(i)\ne f(j)\)。给定 \(m\)，求出仅包含数字 \(1\sim m\) 的好的序列的数量取模。\(m\le10^6\)，多测 \(T\le3\times10^5\)
可证 \(a\) 是好的当且仅当 \(a\) 单谷 且 将其从大到小排序后其任意前缀 \(\gcd\) 互不相同
先考虑其弱化版，\(\sum m\le10^5\) 的情况（即 CF2039F1 Shohag Loves Counting (Easy Version)），对于每组 \(m\) 依次求解
对于一个从大到小排序了的 \(a_{1\sim n}\)，其排序前有 \(2^{n-1}\) 种可能（每个 \(a_{1\sim {n-1}}\) 都可以在 \(a_n\) 左侧或右侧）
令 \(f_{i,j}\) 为排序后的 \(a\) 前缀 \(\gcd\) 等于 \(i\)，长为 \(j\) 的方案数（显然 \(j\) 为 \(O(\log n)\) 的）
由于要求序列递减，因此从 \(n\) 到 \(1\) 加入最后一个数 \(x\)，考虑其对整个 \(f\) 数组的贡献
显然只会影响 \(i\mid x\) 的 \(f_{i,j}\)，其余不变
令 \(f'\) 为修改后的 \(f\)，则有：
答案为 \(\sum_i\sum_j2^{j-1}f_{i,j}\)
直接做 为 \(O(m^2+\sum m^{3/2}\log^2 n)\) 的（预处理了 \(\gcd\)）
发现最后 \(2^{j-1}\) 的系数可以移到转移中：
答案为 \(\sum_i\sum_jf_{i,j}\)
此时可以省略第二维，新的 \(f_i\) 存储 \(\sum_j f_{i,j}\)，转移为
答案为 \(\sum f_i\)
这样 时间复杂度为 \(O(\sum n^{3/2}\log^2 n)\)
令 \(s_i=\sum_{i\mid j}f_j\)，则第二种转移变为：
先根据旧的 \(f\) 求出新的 \(f_i\) 相对旧 \(f_i\) 的增量，记为 \(tmp_i\)，然后用 \(tmp_i\) 统一更新 \(f\) 和 \(s\)
预处理所有数的因数列表，则时间复杂度为 \(O(\sum n^{1+\varepsilon})\)（因为 \(d(n)=n^{\varepsilon}\)）
这样足够 通过简单版 了
考虑如何解决原问题
令 \(\vec v_x\) 为加入 \(x\) 后 \(f\)，\(tmp\)，\(s\) 三个数组和 \(1\) 依次拼接起来得到的向量
则 \(\vec v_x\) 到 \(\vec v_{x-1}\) 的转移可以写为一个矩阵
矩阵的转置有性质 \(AB=(B^TA^T)^T\)
因此若将初始向量和每个转移矩阵转置，则原本 \(n\sim 1\) 的递推可以变为 \(1\sim n\) 的正向递推，这样可以一次递推求出 \(m=1\sim 10^6\) 的答案
具体地，定义一个变量 \(one\) 表示原本 \(1\) 的位置，初始 \(s_1=1\)。在之前的算法的基础上，将外层循环顺序改为 \(1\sim n\)，循环内各子语句顺序倒置，原本形如 \(x\leftarrow ay\) 的语句变为 \(y\leftarrow ax\)。\(m=x\) 的答案即为该次循环结束后 \(one\) 的值
时间复杂度 \(O(m\log^2 m+T)\)
例 5：P5518 [MtOI2019] 幽灵乐团 / 莫比乌斯反演基础练习题
给定 \(T\) 组 \(A,B,C\)，每组对 \(t=0/1/2\) 分别求出 \(\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C\left(\dfrac{\operatorname{lcm}(i,j)}{\gcd(i,k)}\right)^{f(t)}\)，其中 \(f(0)=1\)，\(f(1)=ijk\)，\(f(2)=\gcd(i,j,k)\)，答案对 \(p\) 取模，\(T=70,A,B,C\le10^5\)，\(p\) 为足够大的质数
由于
因此转化为两个子问题，求出 \(\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^Ci^{f(t)}\) 和 \(\prod_{i=1}^A\prod_{j=1}^B\prod_{k=1}^C\gcd(i,j)^{f(t)}\)
令 \(Adn(x)=\frac{x(x+1)}2\)
接下来对于 \(t=0/1/2\) 分别解决这两个问题
t = 0
子问题 1
显然
预处理阶乘即可（注意指数要对 \(p-1\) 取模，下同）
子问题 2
可以 \(O(n\log n)\) 预处理 \(Tp\) 数组，\({Tp}_T=\prod_{d\mid T}d^{^{\normalsize\mu\left(\frac Td\right)}}\)，每次询问数论分块
t = 1
子问题 1
显然
预处理 \(i^i\) 的前缀积即可
子问题 2
显然
而
令 \(Tmp(A,B)=\sum_{D=1}^{\min(A,B)}\mu(D)D^2Adn(A//D)Adn(B//D)\)（容易通过预处理加数论分块单次 \(O(\sqrt A+\sqrt B)\) 计算）
则原式 \(\;=\prod_{d=1}^{\min(A,B)}\left(d^{d^2}\right)^{Tmp(A//d,B//d)}\)，可以再次数论分块实现
t = 2
子问题 1
预处理后直接分块即可
子问题 2
其中前一项
显然可以预处理加数论分块
后一项
之前预处理了 \({Tp}_T=\prod_{d\mid T}d^{^{\large\mu\left(\frac Td\right)}}\)，带入可得
令 \(Ftp(n,m)=\prod_F{Tp}_F^{^{\normalsize(n//F)(m//F)}}\)，则上式等于 \(\prod_T\left(Ftp(A//T,B//T)\right)^{^{\normalsize\varphi(T) (C//T)}}\)，容易通过数论分块套数论分块解决
时间复杂度
假定 \(A,B,C\) 同阶，\(n\) 为各组数据中 \(A,B,C\) 的最大值，则总时间复杂度为：
其中 \(\sum_{i=1}^{\sqrt A}\sqrt i\) 显然不超过 \(O(A^{3/4})\)
而
因此总时间复杂度为 \(O(n\log n+TA^{3/4}\log p)\)
例 6：P4449 于神之怒加强版
给定 \(n,m,k\)，求出 \(\sum_{i=1}^n\sum_{j=1}^m\gcd(i,j)^k\)，多次 \(T\le2000\)，\(n,m,k\le5\times10^6\)
根据一般套路，原式可化为 \(\sum_{T=1}^{\min(n,m)}(n//T)(m//T)\sum_{d\mid T}d^k\mu(\frac Td)\)
令 \(g(T)=\sum_{d\mid T}d^k\mu(\frac Td)\)，则 \(g=\operatorname{id}_k\ast \mu\)
由于 \(\operatorname{id}_k\) 和 \(\mu\) 都是积性函数，根据狄利克雷卷积的性质，\(g\) 也是积性函数，因此可以线性筛
例 7：hydro #H1060. 绝世傻逼题
给定 \(t\) 组 \(A,B,C,D,E\)，其中 \(t\le2000,1\le A,B,C,D,E\le5\times10^7\)，分别求出以下式子的值
其中 \([x,y,z]\) 表示取遍 \(\{a,b,c,d,e\}\) 的大小为三的子集，\([x,y]\) 表示取遍大小为二的子集
考虑化简求和中的分式
实际上其等于 \(\gcd(a,b,c,d,e)^6\)，证明需要分别考虑每个质因子，将 \(\gcd\) 和 \(\operatorname{lcm}\) 转化为 \(\min\) 和 \(\max\)，然后假定五个数中该质因子幂次的顺序，化简后即可证
于是转化为求
数论分块即可，时间复杂度 \(O(V+t\sqrt V)\)
vjudge 练习
Number Theory 1，pw: edgeedgeweloveyou
浙公网安备 33010602011771号