4200 字
21 分钟
竞赛算法(4) 数论
2025-11-24
无标签

数论基础#

整除#

具有如下基础性质:

  1. aba | bab-a | b, aba | -b
  2. aba | b, bcb | caca | c
  3. aba | b, aca | caxb+yca | xb + yc
  4. aba | b, bab | ab=ab = ab=ab = -a
  5. aba | bmamb(m0)ma | mb(m\ne0).
  6. aba | b, b0b\ne0ab|a| \le |b|
  7. aba | bab%aa | b\% a

有上述性质可知整除是整数集上的偏序关系. 若a整除b, 则称a为b的约数, b为a的倍数.

0, 1, -1, b, -b称为平凡约数, 其他约数称为真约数.

b0b\ne0, 当dd遍历bb的全体约数时, b/db/d也遍历bb的全体约数.

带余除法#

b=(b/a)a+b%ab = (b / a)a + b \% a.

相邻的 a 个整数被正整数 a 除后, 恰好取到上述 a 个余数. 特别地, 一定有且仅有一个数被 a 整除.

gcd与lcm#

gcd有如下性质:

  1. gcd(a, b) = gcd(-a, b), 正负不影响结果
  2. gcd(a, b) = gcd(b, a), 可交换
  3. gcd(a, 0) = gcd(a, a) = abs(a)
  4. gcd(bq + r, b) = gcd(r, b)
  5. gcd(a, b, c) = gcd(gcd(a, b), c), 可结合
  6. gcd(a, b, c...) * m = gcd(ma, mb, mc...)
  7. gcd(a^n, b^n) = gcd(a, b)^n
  8. gcd(a, b, c...) = d, 则gcd(a / d, b / d, c / d...) = 1
  9. a | bc, a, b互质, 则a | c.
  10. a | c, b | c, a, b互质, 则ab | c.
  11. a, b互质, 则gcd(a, bc) = gcd(a, c)
  12. A, B数组组间两两互质, 即任意a, b互质, 则mul(A), mul(B)互质.
  13. 若互质的n个数相乘等于某个数vm次方, 则这n个数每个都是完全m次方数.

lcm有如下性质:

  1. 同gcd
  2. 同gcd
  3. 同gcd
  4. a | b, 则lcm(a, b) = b
  5. 可结合
  6. 同gcd
  7. 若数组A每个数a都有a | m, 则lcm(A) | m.
  8. lcm(a, b, c)lcm(ab, bc, ca) = lcm(a, b)lcm(b, c)lcm(a, c)
  9. lcm(a^n, b^n) = lcm(a, b)^n

有一些奇怪的等式:

  1. gcd(a, b)lcm(a, b) = ab
  2. gcd(ab, bc, ca)lcm(a, b, c) = abc
  3. gcd(a, b, c)gcd(a, b, c) / gcd(a, b)gcd(a, c)gcd(b, c) = lcm[a, b, c]lcm[a, b, c] / lcm(a, b)lcm(a, c)lcm(b, c)

质数#

所有大于3的素数都可以表示为6n + 16n + 5的形式.

算术基本定理#

算术基本引理: 若p为素数, p | xy, 则要么p | x, 要么p | y. 算术基本定理/唯一分解定理: a可被唯一的素数数组P表示, a = mul(P) 这两个定理等价.

同余#

m | (a - b)a, bm同余, 记作ab(modm)a\equiv b\pmod m 同余是一种等价关系. 同余的性质:

  1. 自反, 对称, 传递
  2. 线性运算(加减乘)
  3. 若两个整系数多项式的所有系数模m同余, 则这两个多项式模m同余意义下等价. (对于任意s, tm同余, f(s)g(t)m同余)
  4. a, bm同余, 则ak, bkmk同余. 若k整除a, b, m, 则反过来也一样.
  5. d | m, a, bm同余, 则a, bd同余.
  6. a, bm同余, 则gcd(a, m) = gcd(b, m).
  7. a, bm同余, d | md | a, 则d | b.

同余类与剩余系#

模运算算法#

模运算环#

namespace numth {
inline ll norm(ll x, ll mod) { ll tmp = x % mod; return (tmp < 0) ? (tmp + mod) : tmp; }
inline ll mul(ll a, ll b, ll mod) { auto tmp = (i128)a * b % mod; return (tmp < 0) ? (tmp + mod) : tmp; }
inline ll add(ll a, ll b, ll mod) { return norm(a + b, mod); }
inline ll sub(ll a, ll b, ll mod) { return norm(a - b, mod); }
inline ll pow(ll a, ll b, ll mod) {
ll ret = 1;
a = norm(a, mod);
for (; b; b >>= 1, a = mul(a, a, mod)) if (b & 1) ret = mul(ret, a, mod);
return ret;
}
inline ll inv(ll a, ll mod_prime) { return pow(a, mod_prime - 2, mod_prime); }
}

在模数为质数的情况下, 模运算构成一个环. norm: 模. 模运算的结果总是在[0, mod)之间.

模整数类#

代码比模运算环更长, 但健壮性更好(支持合数模数), 且速度更快. 静态u32(最快):

template <int MOD>
struct modint {
static constexpr u32 umod = u32(MOD);
static constexpr int mod = MOD;
static constexpr int get_mod() { return mod; }
static_assert(umod < u32(1) << 31);
u32 val;
// fast construction of raw modint
static modint raw(u32 v) {
modint x;
x.val = v;
return x;
}
// constructors
constexpr modint() : val(0) {}
constexpr modint(u32 x) : val(x % umod) {}
constexpr modint(u64 x) : val(x % umod) {}
constexpr modint(u128 x) : val(x % umod) {}
constexpr modint(int x) : val((x %= mod) < 0 ? x + mod : x) {};
constexpr modint(ll x) : val((x %= mod) < 0 ? x + mod : x) {};
constexpr modint(i128 x) : val((x %= mod) < 0 ? x + mod : x) {};
// operators
bool operator<(const modint &other) const { return val < other.val; }
modint &operator+=(const modint &p) {
if ((val += p.val) >= umod) {
val -= umod;
}
return *this;
}
modint &operator-=(const modint &p) {
if ((val += umod - p.val) >= umod) {
val -= umod;
}
return *this;
}
modint &operator*=(const modint &p) {
val = u64(val) * p.val % umod;
return *this;
}
modint &operator/=(const modint &p) {
// O(log mod)
assert(p.val != 0 && "modint div zero");
*this *= p.inv();
return *this;
}
modint operator-() const { return modint::raw(val ? mod - val : u32(0)); }
modint operator+(const modint &p) const { return modint(*this) += p; }
modint operator-(const modint &p) const { return modint(*this) -= p; }
modint operator*(const modint &p) const { return modint(*this) *= p; }
modint operator/(const modint &p) const { return modint(*this) /= p; } // O(log mod)
bool operator==(const modint &p) const { return val == p.val; }
bool operator!=(const modint &p) const { return val != p.val; }
// quick power of O(log n)
modint qkpow(ll n) const {
if (n < 0) return inv().qkpow(-n);
modint ret(1), mul(val);
while (n > 0) {
if (n & 1) ret *= mul;
mul *= mul;
n >>= 1;
}
return ret;
}
// inv by extgcd of O(log mod)
modint inv() const {
// no inv for zero, no inv for non-coprime
assert(val != 0 && "modint inv zero");
assert((gcd((int)val, mod) == 1) && "modint inv not coprime");
int a = val, b = mod, u = 1, v = 0, t;
while (b > 0) {
t = a / b;
swap(a -= t * b, b), swap(u -= t * v, v);
}
return modint(u);
}
// stream IO
friend std::istream &operator>>(std::istream &is, modint &a) {
ll t;
is >> t;
a = modint(t);
return is;
}
friend std::ostream &operator<<(std::ostream &os, const modint &a) {
return os << a.val;
}
};
using mint107 = modint<1000000007>;
using mint998 = modint<998244353>;

O2有自动Barrett优化, 而且不涉及i128, 常数非常小. 且支持所有常见运算. 无依赖.

Barrett算法#

struct Barrett {
u32 m;
u64 im;
explicit Barrett(u32 m = 1) : m(m), im(u64(-1) / m + 1) {}
u32 umod() const { return m; }
u32 modulo(u64 z) {
if (m == 1) return 0;
u64 x = (u64)(((unsigned __int128)(z)*im) >> 64);
u64 y = x * m;
return (z - y + (z < y ? m : 0));
}
u64 floor(u64 z) {
if (m == 1) return z;
u64 x = (u64)(((unsigned __int128)(z)*im) >> 64);
u64 y = x * m;
return (z < y ? x - 1 : x);
}
pair<u64, u32> divmod(u64 z) {
if (m == 1) return {z, 0};
u64 x = (u64)(((unsigned __int128)(z)*im) >> 64);
u64 y = x * m;
if (z < y) return {x - 1, z - y + m};
return {x, z - y};
}
u32 mul(u32 a, u32 b) { return modulo(u64(a) * b); }
};

模整数类(动态)#

基于手写Barrett算法. 支持u32范围内的模数(更大的模数需要手写Barrett64Dynamic_modint64)

template <int id>
struct Dynamic_Modint {
static constexpr bool is_modint = true;
using mint = Dynamic_Modint;
u32 val;
static Barrett bt;
static u32 umod() { return bt.umod(); }
static int get_mod() { return (int)(bt.umod()); }
static void set_mod(int m) {
assert(1 <= m);
bt = Barrett(m);
}
static Dynamic_Modint raw(u32 v) {
Dynamic_Modint x;
x.val = v;
return x;
}
Dynamic_Modint() : val(0) {}
Dynamic_Modint(u32 x) : val(bt.modulo(x)) {}
Dynamic_Modint(u64 x) : val(bt.modulo(x)) {}
Dynamic_Modint(int x) : val((x %= get_mod()) < 0 ? x + get_mod() : x) {}
Dynamic_Modint(ll x) : val((x %= get_mod()) < 0 ? x + get_mod() : x) {}
Dynamic_Modint(i128 x) : val((x %= get_mod()) < 0 ? x + get_mod() : x){};
mint& operator+=(const mint& rhs) {
val = (val += rhs.val) < umod() ? val : val - umod();
return *this;
}
mint& operator-=(const mint& rhs) {
val = (val += umod() - rhs.val) < umod() ? val : val - umod();
return *this;
}
mint& operator*=(const mint& rhs) {
val = bt.mul(val, rhs.val);
return *this;
}
mint& operator/=(const mint& rhs) { return *this = *this * rhs.inv(); }
mint operator-() const { return mint() - *this; }
mint pow(ll n) const {
assert(0 <= n);
mint x = *this, r = 1;
while (n) {
if (n & 1) r *= x;
x *= x, n >>= 1;
}
return r;
}
mint inv() const {
int x = val, mod = get_mod();
int a = x, b = mod, u = 1, v = 0, t;
while (b > 0) {
t = a / b;
swap(a -= t * b, b), swap(u -= t * v, v);
}
if (u < 0) u += mod;
return u;
}
friend mint operator+(const mint& lhs, const mint& rhs) { return mint(lhs) += rhs; }
friend mint operator-(const mint& lhs, const mint& rhs) { return mint(lhs) -= rhs; }
friend mint operator*(const mint& lhs, const mint& rhs) { return mint(lhs) *= rhs; }
friend mint operator/(const mint& lhs, const mint& rhs) { return mint(lhs) /= rhs; }
friend bool operator==(const mint& lhs, const mint& rhs) { return lhs.val == rhs.val; }
friend bool operator!=(const mint& lhs, const mint& rhs) { return lhs.val != rhs.val; }
};
using dmint = Dynamic_Modint<-1>;
template <int id>
Barrett Dynamic_Modint<id>::bt;

数论函数#

正整数到整数(在广义数论中可以是复数)的函数. 可以被看作一个数列.

积性函数#

  1. f(1)=1f(1) = 1.
  2. 对于互质的xx, yy, 有f(xy)=f(x)f(y)f(xy) = f(x)f(y).
完全积性函数
  1. f(1)=1f(1) = 1.
  2. 对于任意的xxyy, 有f(xy)=f(x)f(y)f(xy) = f(x)f(y).
积性函数的性质

f(x),g(x)f(x), g(x)都是积性函数:

  1. f(xp)f(x^p)也是积性函数.
  2. fp(x)f^p(x)也是积性函数.
  3. f(x)g(x)f(x)g(x)也是积性函数.
  4. dxf(d)g(xd)\displaystyle \sum _{d|x}f(d)g(\frac{x}{d})也是积性函数.
  5. 容易注意到正整数在唯一分解后, 积性函数等于其每一个质因数幂次的函数值相乘, 完全积性函数还等于其每一个质因数的函数值相乘.
常见的积性函数
  1. 单位函数ϵ(n)=[n=1].\epsilon (n) = [n = 1].n=1n = 1的指示函数. 为完全积性函数.
  2. 恒等函数idk(n)=nk,id(n)=nid_k(n) = n^k, id(n) = n. 完全积性.
  3. 常数函数1(n)=11(n) = 1. 完全积性.
  4. 除数函数σk(n)=dndk\displaystyle \sigma _k(n) = \sum _{d|n}d^k. 其中σ0(n)\sigma _0(n)常记作d(n)d(n)τ(n)\tau(n), σ1(n)\sigma_1(n)常记作σ(n)\sigma(n).
  5. 欧拉函数φ(n)=i=1n[gcd(i,n)=1]\varphi(n) = \displaystyle\sum_{i = 1}^{n}[gcd(i, n) = 1].
  6. 莫比乌斯函数μ(n)\mu (n). nn11时值为11, nn为完全平方数的整数倍时值为00, 其他时候值为(1)k(-1)^k, 其中kknn的不相等质因子个数.

欧拉函数#

一种数论函数, 记作φ(n)\varphi(n), 是一种积性函数, 其含义为小于等于n的正整数中与n互质的数的个数.

欧拉函数的性质
  1. 欧拉函数是积性函数. 特别地, 当nn为奇数时, φ(2n)=φ(n)\varphi(2n) = \varphi(n).
  2. n=dnφ(d)n = \displaystyle\sum_{d|n}\varphi(d). 即恒等函数idid是常数函数11和欧拉函数φ\varphi的狄利克雷卷积.
  3. n=pkn = p^k, 则φ(n)=pkpk1\varphi(n) = p^k - p^{k - 1}.
  4. φ(n)\varphi(n)等于nn除以所有不同的质因数, 再乘以所有不同质因数-1的结果.

莫比乌斯函数#

一种数论函数, 记作μ(n)\mu(n):

  1. n=1n = 1μ(n)=1\mu(n) = 1.
  2. nn的质因数分解形如p1p2p3pkp_1p_2p_3\cdots p_k, 其中pip_i之间互不相等时, μ(n)=(1)k\mu(n) = (-1)^k.
  3. 其他时候, μ(n)=0\mu(n) = 0.

莫比乌斯函数是一种积性函数.

莫比乌斯函数可以使用质因数分解在O(n)O(\sqrt{n})时间内求出. 也可使用线性筛在O(n)O(n)时间内预处理出1n1\sim n内所有的莫比乌斯函数值.

欧拉筛#

一种O(n)O(n)的筛法, 也叫线性筛.

:::[欧拉筛处理质数的过程]

  1. 维护质数表primes(要求支持O(1)O(1)尾插), 最小质因数表minp.
  2. 循环, 到i时(先判后筛):
    1. 判: 若minp[i] == 0, 说明没被筛过, 为质数. minp[i] = i, primes尾插i.
    2. 筛: 遍历所有小于等于minp[i]的质数pj. 我们考虑i * pj, 可以知道pj一定是i * pj的最小质因数(由唯一分解定理). 于是记录minp[i * pj] = pj, 表示i * pj已经被筛过. :::

这个过程中, minp可以换成一个bool数组vis以节省空间, 这样的话筛的时候条件改为”从小到大遍历所有i % pj != 0的质数和第一个i % pj == 0的质数”. 常数时间接近, 空间节省, 但缺少了副产物minp表.

欧拉筛用于预处理φ(n)\varphi(n), μ(n)\mu(n), d(n)d(n), σ(n)\sigma(n)等特殊积性函数时, 利用其特殊性质可以避开处理一般积性函数时较差的时空复杂度.

using i8 = int8_t;
using ll = long long;
// 全局
int minp[N]; // minp[i]反映i的最小质因数
/*
int phi[N], d[N];
i8 mu[N];
int times[N]; // 求d用, 表示最小质因子的次数.
ll sigma[N];
int g[N]; // 求sigma用, 表示最小质因子的幂次贡献的约数和.
*/
// 局部
vector<int> primes; // 质数表, primes[0] = 2
auto euler = [&](int n) -> void {
primes.clear();
memset(minp, 0, sizeof(int) * (n + 3));
/*
memset(phi, 0, sizeof(int) * (n + 3));
memset ... 以下省略
phi[1] = mu[1] = d[1] = sigma[1] = 1;
*/
for (int i = 2; i <= n; i++) {
if (minp[i] == 0) {
primes.push_back(i);
minp[i] = i;
/*
phi[i] = i - 1;
mu[i] = -1;
d[i] = 2;
times[i] = 1;
sigma[i] = i + 1;
g[i] = i + 1;
*/
}
for (const auto& pj : primes) {
if (pj > minp[i] || pj > n / i) break;
/*
if (pj == minp[i]) {
phi[i * pj] = phi[i] * pj;
mu[i * pj] = 0;
d[i * pj] = d[i] / (times[i] + 1) * (times[i] + 2);
times[i * pj] = times[i] + 1;
g[i * pj] = g[i] * pj + 1;
sigma[i * pj] = sigma[i] / g[i] * g[i * pj];
}
else {
phi[i * pj] = phi[i] * (pj - 1);
mu[i * pj] = -mu[i];
d[i * pj] = d[i] * 2;
times[i * pj] = 1;
sigma[i * pj] = sigma[i] * sigma[pj];
g[i * pj] = pj + 1;
}
*/
minp[i * pj] = pj;
}
}
};
// 判定质数
auto isprime = [&](int x) -> bool {
return x == minp[x];
};

欧拉筛求一般积性函数#

若可以在关于kk的低次多项式时间(在算法竞赛中, 可以忽略”低次”, 因为kk可以看作常数)内得到积性函数ffpkp^k处的值(其中pp为质数), 便可在O(n)O(n)内使用线性筛预处理出1n1\sim n范围内所有的f(n)f(n)值.

欧拉筛预处理积性函数的过程
  1. 维护质数表primes(O(1)O(1)尾插). <int>minp<bool>vis用于标记. low用于记录i的最小质因子的完整幂. f用于储存答案. 定义函数getf(int p, int pk)用于快速求f(pk)f(p^k). 初始化f[1] = 1, low[1] = 1.
  2. 循环, 到i时(先判后筛):
    1. 判: 若为质数, 更新4项(包括primes).
    2. 筛, 其实就是3个分支, (注意pj > n / i时要提前break):
      1. pj < minp[i], 更新3项.
      2. pj == minp[i], 更新low[tgt]minp[tgt], 判断, 最后break:
        1. i == low[i], 说明i是质数整幂, 调用getf;
        2. i != low[i], 直接计算f.
// 全局
int minp[N], low[N];
ll f[N];
// 局部
vector<int> primes;
auto euler = [&](int n) -> void {
primes.clear();
memset(minp, 0, sizeof(int) * (n + 3));
memset(low, 0, sizeof(int) * (n + 3));
memset(f, 0, sizeof(int) * (n + 3));
f[1] = 1, low[1] = 1;
for (int i = 2; i <= n; i++) {
if (minp[i] == 0) {
primes.push_back(i);
low[i] = i;
minp[i] = i;
f[i] = getf(i, i);
}
for (const auto& pj : primes) {
if (pj > n / i) break;
int tgt = pj * i;
minp[tgt] = pj;
if (pj != minp[i]) {
low[tgt] = pj;
f[tgt] = f[pj] * f[i];
}
else {
low[tgt] = pj * low[i];
if (i == low[i]) {
f[tgt] = getf(pj, low[tgt]);
}
else {
f[tgt] = f[tgt / low[tgt]] * f[low[tgt]];
}
break;
}
}
}
}

预处理后快速分解质因数#

在预处理出minp数组后, 可以在O(logn)O(\log n)的时间内完成质因数的单次询问. 不断除以minp[i]即可.

vector<int> getfac(int x) {
vector<int> ans;
while (x > 1) {
ans.push_back(minp[x]);
x /= minp[x];
}
return ans;
}

整除分块(数论分块)#

O(1)O(1)得到fpref_{pre}gg的前提下, 以O(n)O(\sqrt{n})的复杂度计算i=1limitf(i)g(ni)\displaystyle\sum_{i=1}^{limit}f(i)g(\lfloor \frac{n}{i}\rfloor).

template <typename F>
void block_floor(ll N, F f) {
ll sq = sqrtl(N);
ll n = (sq * sq + sq <= N ? sq : sq - 1);
for (ll l = 1; l <= sq; ++l) {
f(N / l, l, l + 1);
}
for (ll q = n; q >= 1; --q) {
f(q, N / (q + 1) + 1, N / q + 1);
}
}

狄利克雷卷积#

数论函数f(n)f(n)g(n)g(n)的狄利克雷卷积记作fgf*g, 定义为knf(k)g(nk)=kl=nf(k)g(l)\displaystyle\sum_{k|n}f(k)g(\frac{n}{k}) = \sum_{kl=n}f(k)g(l).

常见的狄利克雷卷积
  1. 单位函数ε\varepsilon是默比乌斯函数μ\mu和常数函数11的狄利克雷卷积(莫比乌斯反演的核心):ε=μ1    ε(n)=dnμ(d)\varepsilon = \mu * 1 \iff \varepsilon(n) = \displaystyle\sum_{d|n}\mu(d)
  2. 除数函数dd是常数函数1111的狄利克雷卷积:d(n)=dn1d(n) = \displaystyle\sum_{d|n}1
  3. 除数和函数σ\sigma是恒等函数idid11的狄利克雷卷积:σ(n)=dnd\sigma(n) = \displaystyle\sum_{d|n}d.
  4. 欧拉函数φ\varphi是恒等函数idid和莫比乌斯函数μ\mu的狄利克雷卷积:φ(n)=dndμ(nd)\varphi(n) = \displaystyle\sum_{d|n}d\mu(\frac{n}{d})
  5. gcd和函数P(n)=i=1ngcd(i,n)P(n) = \displaystyle\sum_{i=1}^ngcd(i, n)是恒等函数idid与欧拉函数φ\varphi的狄利克雷卷积(所谓欧拉反演的核心).
狄利克雷卷积的性质
  1. 交换律
  2. 结合律
  3. 对加法的分配律
  4. 存在单位元, 为单位函数ε\varepsilon
  5. f(1)0f(1)\ne 0时存在逆元, 称为狄利克雷逆元. 逆元满足递推公式g(n)=ε(n)kl=n,k1f(k)g(l)f(1)\displaystyle g(n) = \frac{\varepsilon(n) - \sum_{kl = n, k \ne 1}f(k)g(l)}{f(1)}

如上性质表明全体数论函数在加法运算和狄利克雷卷积运算下构成交换环, 称为狄利克雷环.

莫比乌斯反演#

由狄利克雷卷积可以知道μ1=ε\mu * 1 = \varepsilon, 即在狄利克雷环上, μ\mu11互为逆元. 我们可以利用这个性质得到: 若F=f1F = f * 1, 则f=Fμf = F * \mu. 更加正式地: F(n)=dnf(d),f(n)=dnμ(d)F(nd)\text{若}\displaystyle F(n) = \sum_{d|n}f(d), \text{则}f(n) = \sum_{d|n}\mu(d)F(\frac{n}{d}).

这个过程称作莫比乌斯反演. 这提醒我们: 对于某些函数f(n)f(n), 如果我们很难直接求出它的值, 但却容易求出其倍数和或约数和g(n)g(n)的值, 那么我们可以通过莫比乌斯反演简化运算, 从而求得f(n)f(n)的值.

竞赛算法(4) 数论
https://fuwari.vercel.app/posts/竞赛算法4-数论/
作者
ykindred
发布于
2025-11-24
许可协议
CC BY-NC-SA 4.0