我与数论不共戴天

将当时高中学的稀烂的东西重新整合一下。

希望对自己有所帮助。

快速幂

oi-wiki 并没有将快速幂分到数论中。

但是毕竟很典,而且很重要。

模板链接

快速幂,用于解决形如 \(a^b\bmod p\) 的问题。

朴素的想法是,直接按照指数 \(b\)for

int pow(int a, int b, int p)
{
    int res = 1;
    for (int i = 1; i <= b; i++)
        res = res * a % p;
    return res;
}

但是时间复杂度为 \(O(b)\),常见题目中 \(b\) 一般为 \(1e9\) 范围。

此时考虑二进制优化的快速幂。

将指数 \(b\) 写为二进制的形式,并且求出此时对应的 \(a\) 的幂次。

\(3^{10}\) 可以写成 \(3^{1010_2}=3^{10_2}\times 3^{1000_2}=3^2\times 3^8\)

此时便可通过判断 \(b\) 的二进制位上是否为 \(1\) 来确定 \(a\) 的该幂次是否在答案中。

int ksm(int a, int b, int p)
{
    int res = 1;
    while (b)
    {
        if (b & 1)
            res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
}

习题

矩阵乘与矩阵快速幂

同样,这块应该属于线性代数而不是数论()

矩阵乘

关于矩阵乘,只有当前一个矩阵的行数和后一个矩阵的列数相等时才有意义。

对于两矩阵 \(A_{n\cdot m}\)\(B_{m\cdot p}\) 相乘的公式为:

\[c_{i,j}=\displaystyle\sum_{k=1}^m a_{i,k}\times b_{k,j} \]

仔细想想,在实数中,是不是有 \(n\times 1=n\)

那么,在矩阵中,是不是也有类似于 \(1\) 的“ \(1\) 矩阵”呢?

答案是有的,但它不叫“ \(1\) 矩阵”,而是叫 单位矩阵 ,经常用 \(I\) 来表示。

单位矩阵满足主对角线为 1 ,其他值均为 0

大概长这个样子:

\[I= \begin{pmatrix} 1 & 0 & \cdots & 0 & 0\\ 0 & 1 & \ddots & 0 & 0\\ \vdots & \ddots & \ddots & \ddots & \vdots\\ 0 & 0 & \ddots & 1 & 0\\ 0 & 0 & \cdots & 0 & 1 \end{pmatrix}\]

根据上边的公式手推就能得到任何矩阵乘单位矩阵都还是原矩阵的原因。

一般来说,计算机中的矩阵是用二维数组模拟的。

两矩阵 \(A_{n\cdot m}\)\(B_{m\cdot p}\) 相乘的代码为:

for (int i = 1; i <= n; i++)
    for (int j = 1; j <= n; j++)
        for (int k = 1; k <= n; k++)
            c[i][j] += a[i][k] * b[k][j];

矩阵乘法一般都是配合着快速幂,叫做 矩阵快速幂

而矩阵快速幂最常见的应用就是加速递推或 dp,将 \(O(n)\) 的线性递推优化成 \(O(\log n)\)

矩阵快速幂

矩阵快速幂,用于快速求解一个矩阵 \(n\)\(k\) 次幂。

模板

首先可以确定: 矩阵乘+快速幂=矩阵快速幂

然后。

全文完。

其实就是根据快速幂原理,将其中的数乘换成矩阵乘。

为方便,将矩阵封装在结构体中

struct Matrix
{
    int n, m; // 行、列
    i64 mode; // 模数
    vector<vector<i64>> A;

    Matrix(int n_, int m_, i64 mod_, bool mat_1 = 0)
        : n(n_), m(m_), mode(mod_), A(n_, vector<i64>(m_, 0))
    {
        if (mat_1)
        {
            for (int i = 0; i < min(n_, m_); i++)
                A[i][i] = 1;
        }
    }

    static Matrix mat_1(int n, i64 mode)
    { // 单位矩阵
        return Matrix(n, n, mode, 1);
    }

    vector<i64> &operator[](int i) { return A[i]; }
    const vector<i64> &operator[](int i) const { return A[i]; }

    Matrix operator*(const Matrix &b) const
    { // 矩阵乘法
        assert(m == b.n && mode == b.mode);
        Matrix c(n, b.m, mode);
        for (int i = 0; i < n; i++)
            for (int k = 0; k < m; k++)
                for (int j = 0; j < b.m; j++)
                    c[i][j] = (c[i][j] + A[i][k] * b[k][j]) % mode;
        return c;
    }

    Matrix ksm(i64 b) const
    { // 矩阵快速幂
        assert(n == m);
        Matrix ret = mat_1(n, mode);
        Matrix a = *this;
        while (b > 0)
        {
            if (b & 1)
                ret = ret * a;
            a = a * a;
            b >>= 1;
        }
        return ret;
    }

    void print() const
    {
        for (int i = 0; i < n; i++)
        {
            for (int j = 0; j < m; j++)
                cerr << A[i][j] << ' ';
            cerr << '\n';
        }
    }
};

结合代码来看,Matrix operator* 将矩阵乘重载,Matrix.ksm 是用矩阵快速幂实现自身 \(b\) 次幂。

由此,可以得到矩阵加速的两种操作:

void ak()
{// 求 A ^ k
    int n;
    i64 k;
    cin >> n >> k;
    Matrix a(n, n, mod);
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            cin >> a[i][j];

    a = a.ksm(k);

    int ans = 0;
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            ans = (ans + a[i][j]) % mod;

    cout << ans << '\n';
}
void fib()
{ // 求斐波那契
    i64 n;
    cin >> n;
    Matrix a(2, 1, mod);
    Matrix f(2, 2, mod);

    a.A = {{1}, {1}};
    f.A = {{1, 1}, {1, 0}};

    if (n > 2)
    {
        f = f.ksm(n - 2);
        a = f * a;
    }

    cout << a[0][0] << '\n';
}

可爱的 练习题们:

  1. 矩阵快速幂
  2. 【模板】矩阵加速
  3. 斐波那契数列
  4. 广义斐波那契数列
  5. 随机数生成器
  6. 数列
  7. 刷题比赛

分解质因数

唯一分解定理:

对于任意的整数 \(n\) ,有 \(n=\displaystyle\prod_{i=1}^k p_i^{e_i}\) ,其中 \(p_i\) 均为质数, \(e_i\) 均为正整数。

对于 \(n\),在分解质因数的时候只需找到符合条件的 \(p_i\),然后将其中的 \(e_i\) 一步一步的除至 \(0\)。这样的时间复杂度为 \(O(\sqrt n)\)

void factorization(int n)
{
    for (int i = 2; i * i <= n; i++)
    {
        if (n % i == 0)
        {
            p.push_back(i), e.push_back(0);
            while (n % i == 0)
                n /= i, e.back()++;
        }
    }
    if (n > 1)
        p.push_back(n), e.push_back(1);
}

质数相关

试除法

根据定义 for 一遍。由于 \(i\)\(\dfrac{n}{i}\)是对称的,所以只需要枚举到 \(\sqrt{n}\) 即可,时间复杂度 \(O(n\sqrt{n})\)

bool is_prime(int n)
{
    if (n < 2)
        return 0;
    for (int i = 2; i * i <= n; i++)
        if (n % i == 0)
            return 0;
    return 1;
}

埃氏筛

原理:质数的倍数绝对不是质数,时间复杂度 \(O(n\log \log n)\)

/*
 * using v = is_not_prime;
 * using p=prime
 */
v[0] = v[1] = 1;
for (int i = 2; i * i <= n; i++)
{
    if (!v[i])
    {
        p.push_back(i);
        for (int j = i << 1; j <= n; j += i)
            v[j] = 1;
    }
}

优点:便于理解

弊端:某些数会筛除多次。

欧拉筛

原理:根据唯一分解定理,确保每个数都只被它的最小质因数筛掉,避免埃氏筛中某些数重复筛除,提高效率。由于每个数只会被其最小质因数筛掉,因此时间复杂度为 \(O(n)\)

/*
 * using v = is_not_prime;
 * using p=prime
 */
void ol_shai()
{
    v[0] = v[1] = 1;
    for (int i = 2; i <= n; i++)
    {
        if (!v[i])
            p.push_back(i);
        for (int j = 0; j < p.size() && i * p[j] <= n; j++)
        {
            v[i * p[j]] = 1;
            if (i % p[j] == 0)
                break;
        }
    }
}

对于新手而言,可能会难以理解这个 break

i % p[j] == 0 成立时,说明 \(i\) 的最小质因数是 \(p_j\)

若令 \(i=k\times p_j\),若继续找到 \(p_{j+1}\),则有 \(i \times p_{j+1} = k \times p_j \times p_{j+1}\)

此时便可发现 \(p_{j+1}\) 并非是 \(i \times p_{j+1}\) 的最小质因数。

因此需要在 i % p[j] == 0break

习题:

P5736 质数筛

P3383 线性筛素数

应用:线性复杂度求一些函数的函数值。

同余简介

\((a-b)\bmod p=0\) ,则 \(a\)\(b\) 在模 \(p\) 意义下同余,记作 \(a\equiv b\pmod p\)(并非碳碳三键)。

说人话就是 \(a\bmod p=b\)

同余有很多性质,但事实上和普通的加减乘并无多大区别。

性质
  1. \(a\equiv a\pmod p\)
  2. \(a\equiv b\pmod p\) ,则 \(b\equiv a\pmod p\)
  3. \(a\equiv b\pmod p\) ,且 \(b\equiv c\pmod p\) ,则 \(a\equiv c\pmod p\)
  4. \(a\equiv b\pmod p\) ,且 \(c\equiv d\pmod p\) ,则 \(a\pm c\equiv b\pm d\pmod p\)
  5. \(a\equiv b\pmod p\) ,且 \(c\equiv d\pmod p\) ,则 \(ac\equiv bd\pmod p\)
  6. 如果 \(ac\equiv bc\pmod p\) ,且 \(c\)\(p\) 互质,则有 \(a\equiv b\pmod p\)
  7. \(a^c\equiv b^c\pmod p\) ,则 \(a\equiv b\pmod p\)
  8. \(a\equiv b\pmod p\) ,且 \(m|p\) ,则 \(a\equiv b\pmod m\)
  9. 若 $a\equiv b\pmod {p_i},i\in\mathbf N_* $,则 \(a\equiv b\pmod {lcm(p_i)}\)

除此之外,同余还有一些让人看不懂的东西:

同余系,完全剩余系,缩剩余系
  • 同余系(剩余系):

    对于所有模 n 余 r 的整数,我们可以将其分为 n 类,

    那么 \(\overline{r}_ n=\{m\in \mathbf{Z}\mid m\cdot n+r\}\) 就为模 n 余 r 的同余系。

    举个例子: \(\overline{4}_ {10}=\{\cdots,-16,-6,4,14,24,\cdots\}\) 就是模 10 余 4 的同余系。

  • 完全剩余系(完系):

    若从 \(\overline{0}_ n,\overline{1}_ n,\overline{2}_ n,\cdots,\overline{(n-1)}_ n\) 中各挑出一个数,便组成了模 n 的完系 \(R(n)\)

    \(R(n)=\{r_0,r_1,r_2,\cdots,r_{n-1}\}\)

    其中, \(r_0\in \overline{0}_ n,r_1\in \overline{1}_ n,r_2\in \overline{2}_ n,\cdots,r_{n-1}\in \overline{(n-1)}_ n\)

    也就是说, n 个模 n 结果不同余的整数便组成了模 n 的完系。而 \(R(n)=\{0,1,2,\cdots,n-1\}\) 则称为模 n 的最小非负完系。

  • 缩剩余系(缩系):

    对于模 n 的完系,取所有与 n 互质的数,即为模 n 的缩系 \(\Phi_n\)

    \(\Phi_n=\{c_1,c_2,\cdots,c_{\varphi(n)}\}\)

    若缩系 \(\Phi_n\) 满足 \(c_1,c_2,\cdots,c_{\varphi(n)}\in[1,n-1]\) ,那么就称为模 n 的最小正缩系。

这些东西是我小时候稀里糊涂写的,当时应该也没理解,现在更看不懂了。

逆元

数论中的逆元一般为乘法逆元。

乘法逆元的定义为:

\(a\cdot x\equiv 1\pmod p\),则 \(x\)\(a\) 的乘法逆元,记作 \(x=a^{-1}\)

\(\gcd(a,p)=1\) 互质为 \(a\) 存在模 \(p\) 逆元的充要条件。

求逆元的方式主要有

  1. 费马小定理(特殊情况,最常用);
  2. 拓展欧几里得,即 exgcd(一般情况);
  3. 线性递推(\(O(n)\)\(1\sim n\) 的逆元)。

费马小定理

费马小定理有 \(a^{p-1}\equiv 1\pmod p\)

因此 \(a\times a^{p-2}\equiv 1\pmod p\)

所以 \(a^{p-2}\) 满足刚刚的 \(x\),故 \(a^{p-2}\)\(a\) 的逆元。

使用快速幂即可求解。

证明详见 oiwiki。

P2613 有理数取余

由于题目数据很大,需要在使用快读输入的时候取模,或者以 string 形式读入。

#include <bits/stdc++.h>
#define int long long
using namespace std;
const int mod = 19260817;
int re()
{
    int s = 0, f = 1;
    char ch = getchar();
    while (ch > '9' || ch < '0')
    {
        if (ch == '-')
            f = -1;
        ch = getchar();
    }
    while (ch >= '0' && ch <= '9')
        s = (s * 10 + ch - 48) % mod, ch = getchar();
    return s * f;
}
void wr(int s)
{
    if (s < 0)
        putchar('-'), s = -s;
    if (s > 9)
        wr(s / 10);
    putchar(s % 10 + 48);
}
void wr(int s, char ch) { wr(s), putchar(ch); }
int ksm(int a, int b, int p)
{
    int ans = 1;
    while (b)
    {
        if (b & 1)
            ans = (ans * a) % p;
        a = (a * a) % p;
        b >>= 1;
    }
    return ans % p;
}
signed main()
{
    int a = re(), b = re();
    wr(a * ksm(b, mod - 2, mod) % mod, '\n');
    return 0;
}

exgcd

见后文。

线性递推

先说结论:inv[i] = p - (p / i) * inv[p % i] % p;。即:

\[\text{inv}_i=-\dfrac{p}{i}\cdot \text{inv}_{p\bmod i}\bmod p \]

推导过程

显然 \(p\equiv 0\pmod p\)

\(k\times i + r = p\)

替换得 \(k\times i + r\equiv 0\pmod p\)

其中, \(k=\left\lfloor \dfrac{p}{i}\right\rfloor, r=p\bmod i\)

上式左右同乘 \(i^{-1}\cdot r^{-1}\)

\(k\times r^{-1} + i^{-1}\equiv 0\pmod p\)

移项得 \(i^{-1}\equiv -k\cdot r^{-1}\pmod p\)

\(k,r\) 代入得

\(i^{-1}\equiv -\left\lfloor\dfrac{p}{i}\right\rfloor \times(p\bmod i)^{-1}\pmod p\)

显然, \(p\bmod i<i\)\((p\bmod i)^{-1}\) 会在 \(i^{-1}\) 之前求得

用数组 \(\text{inv}_i\) 表示 \(i\) 的逆元,初始化 \(inv_1=1\)

则递推式为 \(\text{inv}_i=-\dfrac{p}{i}\cdot \text{inv}_{p\bmod i}\bmod p\)

这是我小时候写的推导过程,算是没有跳步吧。

P3811 乘法逆元

cin >> n >> p;
inv[1] = 1;
for (int i = 2; i <= n; i++)
    inv[i] = p - (p / i) * inv[p % i] % p;
for (int i = 1; i <= n; i++)
    cout << inv[i] << '\n';

欧几里得相关

欧几里得算法

欧几里得算法即辗转相除法,用于求两数的 \(\gcd\)

对于任意整数 \(a,b\),不妨设 \(a>b\),有 \(\gcd(a,b)=\gcd(b,a\bmod b)\)

递归代码实现:

int gcd(int a, int b)
{
    return b ? gcd(b, a % b) : a;
}

迭代代码实现:

int gcd(int a, int b)
{
    while (b ^= a ^= b ^= a %= b);
    return a;
}

拓展欧几里得

直接看

前置知识:裴蜀定理

对于任意整数 \(a\)\(b\),一定存在整数 \(x, y\) 满足 \(ax+by=\gcd(a,b)\)

裴蜀定理指出,若 \(ax+by=c\) 有整数解,当且仅当 \(\gcd(a,b)\mid c\)

有解时,可以通过拓展欧几里得算法求出一组特解。

代码如下:

int exgcd(int a, int b, int &x, int &y)
{
    if (!b)
    {
        x = 1, y = 0;
        return a;
    }
    int g = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return g;
}
证明

\(a x_1+b y_1=\gcd(a,b)\)

\(b x_2+(a\bmod b) y_2 =\gcd(b,a\bmod b)\)

由欧几里得算法可知

\(\gcd(a,b) = \gcd(b, a\bmod b)\)

\(a x_1 + b y_1 = b x_2 + (a\bmod b) \cdot y_2\)

\(a\bmod b\) 拆为 \(a - \left \lfloor \dfrac{a}{b} \right \rfloor \cdot b\),带入上式,化简得

\[\begin{array}{rl} b x_2+(a\bmod b) y_2 &= a y_2 + b (x_2-\left\lfloor\dfrac{a}{b}\right\rfloor\cdot y_2)\\ & = a x_1 + b y_1 \end{array} \]

\(x_1 = y_2, y_1 = x_2 - \left\lfloor\dfrac{a}{b}\right\rfloor \cdot y_2\)

对于这个题,还需要考虑输出格式。代码:

a = re(), b = re(), c = re();
int gcd = exgcd(a, b, x_1, y_1);
if (c % gcd)
{ // 无解
    puts("-1");
    return;
}

// 一组特解
x_1 = x_1 * c / gcd;
y_1 = y_1 * c / gcd;

// 步长
dx = b / gcd;
dy = a / gcd;

/**
 * x = x_1 + k * dx >= 1
 * y = y_1 - k * dy >= 1
 * 求出 k 的范围 [le, ri]
 */
int le = ceil((double)(1 - x_1) / dx),
    ri = floor((double)(y_1 - 1) / dy);

if (le <= ri)
{ // 有正整数解
    wr(ri - le + 1, ' ');
    x_min = x_1 + le * dx, y_min = y_1 - ri * dy;
    x_max = x_1 + ri * dx, y_max = y_1 - le * dy;
    wr(x_min, ' '), wr(y_min, ' '), wr(x_max, ' '), wr(y_max, '\n');
}
else
{ // 无正整数解
    x_min = x_1 + le * dx, y_min = y_1 - ri * dy;
    wr(x_min, ' '), wr(y_min, '\n');
}

exgcd 求逆元

上文提到,乘法逆元同样可以用 exgcd 求。

且并无类似费马小定理那样的限制。

乘法逆元即为求 \(x\) 满足 \(a x \equiv 1\pmod p\)

改写一下即为 \(ax+py=1\)

上文提到,乘法逆元存在的充要条件是 \(\gcd(a,p)=1\)

直接套用 exgcd 即可。

int exgcd(int a, int b, int &x, int &y)
{
    if (!b)
    {
        x = 1, y = 0;
        return a;
    }
    int ret = exgcd(b, a % b, y, x);
    y -= a / b * x;
    return ret;
}
int main()
{
    int a, b, x, y;
    cin >> a >> b;
    exgcd(a, b, x, y);
    x = (x % b + b) % b;
    cout << x << '\n';
    return 0;
}

积性函数相关

欧拉函数

欧拉函数 \(\varphi(n)\) 表示在 \(1\sim n\) 中与 \(n\) 互质的数的个数。

显然,对于质数 \(p\) 而言,\(\varphi(p)=p-1\)

对于一般数,有:

\[\varphi(n) = n\displaystyle\prod_{p_i\mid n}\dfrac{p_i-1}{p_i} \]

性质:

  1. 欧拉函数是积性函数。

    即当 \(\gcd(a, b)=1\) 时,\(\varphi(a\cdot b)=\varphi(a)\cdot \varphi (b)\)

  2. 对于任意质数 \(p\) ,有 \(\varphi(p^k)=p^{k-1}\cdot (p-1)\)

对于单个数的欧拉函数值,可以直接根据公式使用质因数分解求得。

而线性求欧拉函数值需要欧拉筛。

void Phi()
{
    v[0] = v[1] = 1;
    phi[1] = 1;
    for (int i = 2; i <= n; i++)
    {
        if (!v[i])
        {
            p.push_back(i);
            phi[i] = i - 1;
        }
        for (int j = 0; j < p.size() && i * p[j] <= n; j++)
        {
            v[i * p[j]] = 1;
            if (i % p[j])
                phi[i * p[j]] = phi[i] * phi[p[j]];
            else
            {
                phi[i * p[j]] = phi[i] * p[j];
                break;
            }
        }
    }
}
  • i % p[j] != 0,说明此时 \(p_j\)\(i\cdot p_j\) 的新质因数,\(\gcd(i,p_j)=1\)

    由积性函数定义,\(\varphi(i\cdot p_j)=\varphi(i)\cdot \varphi(p_j)\)

  • i % p[j] == 0,说明 \(p_j\)\(i\) 的质因数,令 \(i=t\cdot p_j^k\)。那么有:

    \[\begin{array}{rl} \varphi(i\cdot p_j)&=\varphi(t\cdot p_j^{k+1})\\ &=\varphi(t) \cdot \varphi(p_j^{k+1})\\ &=\varphi(t) \cdot \varphi(p_j^{k})\cdot p_j\\ &=\varphi(i)\cdot p_j \end{array} \]

莫比乌斯函数

这里仅简单介绍,详细见后文莫比乌斯反演(会更的)。

莫比乌斯函数 \(\mu(n)\),简单来说与 \(n\) 的质因子的数目有关。

\[\mu(n)= \begin{cases} 1&n=1\\ (-1)^k&n=\displaystyle\prod_{i=1}^{k}p_i\\ 0&\text{otherwise} \end{cases} \]

直接看表达式可能不太好理解。

文字语言描述就是:

  • \(n=1\) 时,\(\mu(n)=1\)
  • 对于 \(n=\displaystyle\prod_{i=1}^k p_i^{e_i}\)(唯一分解定理),若 \(\forall e_i=1\),则 \(\mu(n)=(-1)^k\)。换句话说,若 \(n\) 有奇数个不同的质因子,则 \(\mu(n)=-1\),否则 \(\mu(n)=1\)
  • \(n\) 有平方因子时,即 \(\exists e_i>1\),则 \(\mu(n)=0\)

莫比乌斯函数也有很多性质,不过这里不过多介绍。

同为积性函数,也有线性筛求莫比乌斯函数。

void mu_(int n)
{
    mu[1] = 1;
    for (int i = 2; i <= n; i++)
    {
        if (!v[i])
        {
            p.push_back(i);
            mu[i] = -1;
        }
        for (int j = 0; j < p.size() && i * p[j] <= n; j++)
        {
            v[i * p[j]] = 1;
            if (i % p[j])
                mu[i * p[j]] = -mu[i];
            else
            {
                mu[i * p[j]] = 0;
                break;
            }
        }
    }
}

欧拉定理

欧拉定理,也可以叫欧拉降幂定理。

\(a\)\(p\) 互质时,有\(a^{\varphi(p)} \equiv 1 \pmod p\)

这里有一个特例和一个推论。

  • 特例:\(p\) 为质数。

    此时 \(\varphi(p)=p-1\),则 \(a^{p-1}\equiv 1\pmod p\),即费马小定理。

  • 推论:同样,当 \(a\)\(p\) 互质时,有 \(a^{b}\equiv a^{b\bmod \varphi(p)}\pmod p\),所谓的欧拉降幂定理就是这里来的。

    证明:

    \[\begin{array}{rl} a^{b}\bmod p & = a^{\varphi(p) \cdot \left\lfloor\frac{b}{\varphi(p)}\right\rfloor + b \bmod \varphi(p)}\bmod p \\ \\ & = (a^{\varphi(p) \cdot \left\lfloor\frac{b}{\varphi(p)}\right\rfloor}\bmod p) \cdot (a^{b \bmod \varphi(p)}\bmod p) \\ \\ & = 1 \cdot a^{b\bmod \varphi(p)} \bmod p \\ \\ & = a ^ {b \bmod \varphi(p)} \bmod p \end{array}\]

限制:仅当 \(a\)\(p\) 互质时可以使用。

拓展欧拉定理

实际做时,很难保证 \(a\)\(p\) 互质。当不满足互质条件是,需要拓展欧拉定理:

\[a^b\equiv \begin{cases} a ^ {b \bmod \varphi(p)} & \pmod p & \gcd(a, p) = 1 \\ a ^ b & \pmod p & \gcd(a,p) \ne 1 \ \text{AND} \ b < \varphi(p) \\ a ^ {(b \bmod \varphi(p)) + \varphi(p)} & \pmod p & \gcd(a,p) \ne 1 \ \text{AND} \ b \ge \varphi(p) \\ \end{cases} \]

简化版为:

\[a^b\equiv \begin{cases} a ^ b & \pmod p & b < \varphi(p) \\ a ^ {(b \bmod \varphi(p)) + \varphi(p)} & \pmod p & b \ge \varphi(p) \\ \end{cases} \]

int a, p;
cin >> a >> p;

auto Phi = [](int n) -> int
{
    int res = n;
    for (int i = 2; i <= n; i++)
    {
        if (n % i == 0)
        {
            res = res / i * (i - 1);
            while (n % i == 0)
                n /= i;
        }
    }
    if (n > 1)
        res = res / n * (n - 1);
    return res;
};
int phi = Phi(p);

auto re = [&]() -> int
{
    int s = 0;
    bool pd = 0;
    char ch = getchar();
    while (ch > '9' || ch < '0')
        ch = getchar();
    while (ch >= '0' && ch <= '9')
    {
        s = s * 10 + ch - 48, ch = getchar();
        if (s >= phi)
            s %= phi, pd = 1;
        //此处使用为简化版
    }
    return (pd) ? (s + phi) : (s);
};
int b = re();

auto ksm = [](int a, int b, int p) -> int
{
    int res = 1;
    while (b)
    {
        if (b & 1)
            res = (res * a) % p;
        a = (a * a) % p;
        b >>= 1;
    }
    return res;
};
cout << ksm(a, b, p) << '\n';

练习

代码
int n;
cin >> n;
vector<int> a(n + 7);
for (int i = 1; i <= n; i++)
    cin >> a[i];

auto Phi = [](int n) -> int
{
    int res = n;
    for (int i = 2; i * i <= n; i++)
    {
        if (n % i == 0)
        {
            res = res / i * (i - 1);
            while (n % i == 0)
                n /= i;
        }
    }
    if (n > 1)
        res = res / n * (n - 1);
    return res;
};
auto ksm = [](int a, int b, int p) -> int
{
    int res = 1;
    while (b)
    {
        if (b & 1)
            res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
};
auto gcd = [&](auto &gcd, int a, int b) -> int
{
    return b ? gcd(gcd, b, a % b) : a;
};

auto calc = [&](auto &calc, int i, int p) -> int
{//此处使用原三段式
    if (i == n)
        return a[n] %= p;
    int phi = Phi(p), nex = calc(calc, i + 1, phi);
    if (gcd(gcd, a[i], p) == 1)
        return a[i] = ksm(a[i], nex, p);
    else
    {
        if (a[i + 1] < phi)
            return a[i] = ksm(a[i], a[i + 1], p);
        return a[i] = ksm(a[i], nex + phi, p);
    }
};

cout << calc(calc, 1, mod) << '\n';

中国剩余定理

感觉很没用。

模板

设正整数 \(p_1,p_2,p_3, \cdots,p_n\) 两两互质,则同余方程组:

\[\begin{aligned} x \equiv & a_1 \pmod{p_1} \\ x \equiv & a_2 \pmod{p_2} \\ & \vdots{} \\ x \equiv & a_n \pmod{p_n} \end{aligned} \]

有整数解。

并且在模 \(P = \displaystyle \prod_{i = 1} ^ n p_i\) 下的解是唯一的,解为

\[x = \left( \displaystyle \sum_{i = 1} ^ k a_i \cdot P_i \cdot P_i ^ {-1} \right) \bmod P \]

其中, \(P_i=\dfrac{P}{p_i}\)\(P_i^{-1}\)\(P_i\) 在模 \(p_i\) 意义下的逆元。

无法保证 \(p_i\) 为质数,所以 \(P_i^{-1}\) 通过 exgcd 求出。

拓展

不会

组合数相关

高中知识,此处仅说明,不过多赘述。

加法原理、乘法原理;

排列数 \(A_n^m=n(n-1)(n-2)\cdots (n-m+1)=\dfrac{n!}{(n-m)!}\)

全排列 \(A_n^n=n!\)

组合数 \(C_n^m=\dfrac{A_n^m}{m!}=\dfrac{n!}{m!(n-m!)}\),记作 \(\binom{n}{m}\)

二项式定理 \((a+b)^n = \displaystyle \sum_{i=0}^n \binom{n}{i} a^{n-i}b^i\)

组合数性质:

  1. \(\displaystyle \binom{n}{m} = \binom{n}{n-m}\)
  2. \(\displaystyle \binom{n}{k} = \dfrac{n}{k} \binom{n-1}{k-1}\)
  3. \(\displaystyle \binom{n}{m} = \binom{n-1}{m} + \binom{n-1}{m-1}\)
  4. \(\displaystyle \sum _{i=0} ^n \binom{n}{i}=2^n\)

其余详见 oiwiki

卢卡斯定理

规模较小的组合数,可以直接用上述的性质 \(3\) 递推求得,复杂度为 \(O(nk)\);也可直接根据定义,分别计算分子分母的阶乘,复杂度为 \(O(n)\)

规模较大(\(n\sim 10^{18}\))但模数较小时(\(p\sim 10^6\)),便考虑使用 lucas 定理。

\(p\) 为质数时,有

\[\displaystyle \binom{n}{k} \equiv \binom{\lfloor \frac{n}{p} \rfloor}{\lfloor \frac{k}{p} \rfloor} \binom{ n \bmod p }{ k \bmod p } \pmod p\]

显然 \(n \bmod p, k \bmod p < p\),可以用阶乘求得。

\(\displaystyle \binom{\lfloor \frac{n}{p} \rfloor}{\lfloor \frac{k}{p} \rfloor}\) 继续带入 lucas 定理递归求得。

代码

int n, m, p;
cin >> n >> m >> p;
vector<int> f(p);
f[0] = 1;
for (int i = 1; i < p; i++)
    f[i] = f[i - 1] * i % p;

auto ksm = [](int a, int b, int p) -> int
{
    int res = 1;
    while (b)
    {
        if (b & 1)
            res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
};
auto inv = [&](int x) -> int
{
    return ksm(x, p - 2, p);
};

auto C = [&](int n, int m, int p) -> int
{
    if (m < 0 || m > n)
        return 0;
    return f[n] * inv(f[m]) % p * inv(f[n - m]) % p;
};

auto Lucas = [&](auto &Lucas, int n, int m, int p) -> int
{
    if (m == 0)
        return 1;
    return Lucas(Lucas, n / p, m / p, p) * C(n % p, m % p, p) % p;
};

cout << Lucas(Lucas, n + m, n, p) << '\n';

数论分块

数论分块可以解决形如 \(\displaystyle \sum _{i=1} ^n f(i) g(\left\lfloor\frac{n}{i}\right\rfloor)\) 的问题。

其思想和根号分治有点像()

最简单的数论分块问题为,求 \(\displaystyle \sum _{i=1} ^n \left\lfloor\frac{n}{i}\right\rfloor\),等价于求所有 \(1\)\(n\) 的整数的约数个数之和。

直接枚举的话,复杂度为 \(O(n)\)

观察可知,当 \(i\) 逐渐增大时,\(\left\lfloor\dfrac{n}{i}\right\rfloor\) 的变化会越来越慢。

\(\left\lfloor\dfrac{n}{i}\right\rfloor\) 的值出现的分段现象。

设该值为 \(v\),分段长度为 \(l\),便可直接计算整段的贡献为 \(v\times l\)

令当前 \(i\)\(l\),有 \(v=\left\lfloor \dfrac{n}{l} \right\rfloor\),可求得 \(r=\left\lfloor \dfrac{n}{v} \right\rfloor\),即 \(r=\left\lfloor \dfrac{n}{\lfloor \frac{n}{l} \rfloor} \right\rfloor\)

注意到,这样的取值共有 \(\sqrt n\) 个,故复杂度为 \(O(\sqrt n)\)

int D(int n)
{
    int res = 0, l = 1, r = 0;
    while (l <= n)
    {
        r = n / (n / l);
        res += (r - l + 1) * (n / l);
        l = r + 1;
    }
    return res;
}

推广到 \(\displaystyle \sum _{i=1} ^n f(i) g(\left\lfloor\frac{n}{i}\right\rfloor)\) 后,代码变为

int FG(int n)
{
    int res = 0, l = 1, r = 0;
    while (l <= n)
    {
        r = n / (n / l);
        res += (f(r) - f(l - 1)) * g(n / l);
        l = r + 1;
    }
    return res;
}

其中 \(f\) 为可以简单计算前缀和的函数。

莫比乌斯反演

莫比乌斯反演,在已知 \(g(n)\)\(g(n)=\displaystyle\sum _{d\mid n}f(d)\) 时,求出 \(f(n)\)

先说结论:\(f(n)=\displaystyle\sum_{d \mid n}\mu(d) \cdot g(\dfrac{n}{d})\)

\(\mu\) 即为上文提到过莫比乌斯函数:

\[\mu(n)= \begin{cases} 1&n=1\\ (-1)^k&n=\displaystyle\prod_{i=1}^{k}p_i\\ 0&\text{otherwise} \end{cases} \]

莫比乌斯函数有一个非常重要的性质,俗称性质 \(0\)

\[\displaystyle \sum _{d\mid n} \mu(d) = [ n = 1 ] \]

可以通过构造 \(f(n)=[n=1], g(n)=1\),带入上述公式来证明。

注,在数论中,\([条件]\) 代表

\[[条件]= \begin{cases} 1, &\text{满足条件} \\ 0, &\text{否则} \end{cases} \]

先看一道很经典的莫反题

经过题意转化,要求的即为 \(S=\displaystyle \sum _{i=1} ^ {n-1} \sum _{j=1} ^ {n-1} [\gcd(i,j)=1]\)

根据上述性质 \(0\),将 \(\gcd(i,j)\) 替换 \(n\),可转化为:

\[\begin{array}{rl} S &= \displaystyle \sum _{i=1} ^ {n-1} \sum _{j=1} ^ {n-1} [\gcd(i,j)=1] \\ & = \displaystyle \sum _{i=1} ^ {n-1} \sum _{j=1} ^ {n-1} \sum _{d\mid \gcd(i,j)} \mu(d) \\ \end{array} \]

注意到 \(d\mid \gcd(i,j)\)\([d\mid i] [d\mid j]\),因此

\[\begin{array}{rl} S & = \displaystyle \sum _{i=1} ^ {n-1} \sum _{j=1} ^ {n-1} \sum _{d=1} ^{n-1} \mu(d) [d\mid i][d\mid j] \\ & = \displaystyle \sum _{d=1} ^{n-1} \mu(d) \sum _{i=1} ^ {n-1} \sum _{j=1} ^ {n-1} [d\mid i][d\mid j] \\ \end{array} \]

\([d\mid i][d\mid j]\) 代表 \(i, j\) 均为 \(d\) 的倍数,在 \(n-1\) 范围内的 \(d\) 的倍数有 \(\left\lfloor\dfrac{n-1}{d}\right\rfloor\) 个,故原式可整理得

\[\begin{array}{rl} S = \displaystyle \sum _{d=1} ^{n-1} \mu(d) \left\lfloor \dfrac{n-1}{d}\right\rfloor ^2 \\ \end{array} \]

\(\mu(n)\) 可通过欧拉筛求得,复杂度 \(O(n)\)

int n;
cin >> n;
if (n == 1)
    return cout << 0 << '\n', void();
vector<int> v(n), p;
vector<int> mu(n);
v[0] = v[1] = 1;
mu[1] = 1;
for (int i = 2; i < n; i++)
{
    if (!v[i])
    {
        p.push_back(i);
        mu[i] = -1;
    }
    for (int j = 0; j < p.size() && i * p[j] < n; j++)
    {
        v[i * p[j]] = 1;
        if (i % p[j])
            mu[i * p[j]] = -mu[i];
        else
        {
            mu[i * p[j]] = 0;
            break;
        }
    }
}
int ans = 0;
auto pf = [](int x) -> int { return x * x; };
for (int i = 1; i < n; i++)
    ans += mu[i] * pf((n - 1) / i);
cout << ans + 2 << '\n';

莫比乌斯函数和欧拉函数也存在一定得联系。

\[\varphi (n) = \displaystyle \sum _{d\mid n} \mu(d) \dfrac{n}{d} \]

学不懂,鸽了。

posted @ 2026-03-16 18:50  Zvelig1205  阅读(24)  评论(0)    收藏  举报