数论

目录

• 数论
• 逆元
• 1.
• 2.
• 3.
• 扩展欧几里得定理
• 问题引入
• 分析
• 类欧几里得算法
• 问题引入
• 拓展
• 推导 $g$
• 推导 $h$
• Reference
• 中国剩余定理（CRT）
• 问题引入
• 分析
• 证明：
• 拓展中国剩余定理 （EXCRT）
• 问题引入
• 分析
• BSGS
• 问题引入
• 分析
• 扩展 BSGS
• 问题引入
• 分析
• 高斯消元
• 问题引入
• 分析
• 一、积性函数
• 1. 积性函数：
• 2. 完全积性函数：
• 3. 证明 $\varphi(n)$ 是积性函数。
• 二、欧拉函数
• 1. 欧拉函数的定义及求法
• 2. 欧拉函数的一些性质
• 拉格朗日插值
• 多项式定义
• 问题引入 1
• 优化 1
• 优化 2
• 一些常识

回到顶部

数论

逆元

若 $b, m$ 互质，并且 $b | a$，则存在一个整数 $x$ 使得

$$\dfrac{a}{b} \equiv a \times x \pmod m$$

则称 $x$ 为 $b$ 的模 $m$ 乘法逆元，记作 $b^{-1}\pmod m$。

例如 $b=31, m=5, a=93$，存在 $x=1$ 满足 $\dfrac{93}{31} \equiv 93 \times 1 \ (\bmod 5)$。

因为乘法满足模运算的分配律，即 $(a \times b) \bmod m = \left( \left(a \bmod m\right) \times \left(b \bmod m\right) \right ) \bmod m$。

而除法不满足，故乘法逆元就可以很好地解决这个问题，即

$$\dfrac{a}{b} \bmod m=((a \bmod m) \times (b^{-1} \bmod m)) \bmod m$$

这里介绍三种求逆元的方法：

1.

因为

$$\dfrac{a}{b}\equiv a \times b^{-1} \equiv \dfrac{a}{b} \times b \times b^{-1} \pmod m$$

所以

$$b \times b^{-1} \equiv 1 \pmod m$$

其实也就是要找到一个 $x$，满足 $b \times x \equiv 1 \pmod m$。

当 $m$ 为质数时，根据费马小定理 $b^m \equiv b \pmod m$ 可得：$b \times b^{m-2} \equiv 1 \pmod m$。

因此，当模数 $m$ 为质数时，$b^{m-2}$ 就是 $b$ 的模 $m$ 乘法逆元。

2.

第一种方法说过，

其实也就是要找到一个 $x$，满足 $b \times x \equiv 1 \pmod m$。

也就是解 $my + bx = 1$ 这个方程，故我们可以使用扩展欧几里得来求解 $x$。

3.

这里介绍一个线性时间求 $1 \sim n$ 在模 $m$ 下的逆元。

我们设 $m = q \times d + r ~ (r < d)$。

可以得到；

1. $q \times d + r \equiv 0 \pmod m$，
2. $q=\left\lfloor\dfrac{m}{d}\right\rfloor$
3. $r = m \bmod d$

对于式子 1. 两边同乘 $d^{-1} \times r^{-1}$，得到

$$q \times r^{-1} + d^{-1} \equiv 0 \pmod m \\ d^{-1} \equiv -q \times r^{-1} \pmod m \\ d^{-1} \equiv -\left\lfloor\dfrac{m}{d}\right\rfloor \times (m \bmod d) ^ {-1} \pmod m$$

故可以根据前面的 $(m \bmod d) ^ {-1}$ 推出后面的 $d^{-1}$。

#include <bits/stdc++.h>
const int N = 3e6 + 10;
using namespace std;

int n, p;
long long ans[N];

int main()
{
    cin >> n >> p;
    ans[1] = 1;
    printf("%lld\n", 1);

    for (int i = 2; i <= n; ++ i)
        ans[i] = ((- (p / i) * ans[p % i]) % p + p) % p,
        printf("%lld\n", ans[i]);

    return 0;
}

扩展欧几里得定理

问题引入

我们要解决一个这样的问题：

$$ax+by=\gcd(a, b)$$

分析

在欧几里得算法的最后一步，即 $b=0$ 时，显然有一对整数 $x=1, y=0$ 满足 $a\times 1+b \times 0=\gcd(a, 0)$。

我们假设有 $x_2, y_2$ 满足 $bx_2+(a \bmod b)y_2=\gcd(b, a\bmod b)$。

由于 $\gcd(a, b)=\gcd(b, a\bmod b)$，故 $ax+by=bx_2+(a \bmod b)y_2$

又因为 $(a \bmod b) = a - \left\lfloor\dfrac{a}{b}\right\rfloor \times b$。

故

$$ax+by=bx_2+(a - \left\lfloor\dfrac{a}{b}\right\rfloor b)y_2 \\ ax+by=bx_2+ay_2-\left\lfloor\dfrac{a}{b}\right\rfloor by_2 \\ ax+by=ay_2+b(x_2-\left\lfloor\dfrac{a}{b}\right\rfloor y_2)$$

由于 $a=a, b=b$，故 $x=y_2, y=x_2-\left\lfloor\dfrac{a}{b}\right\rfloor y_2$，这样从 $x=1, y=0$ 一步一步就可以推到 $ax+by=\gcd(a, b)$ 的解。

对于一般的方程 $ax+by=c$，当且仅当 $\gcd(a,b) | c$ 时有解。

我们可以先求出 $ax+by=\gcd(a, b)$ 的特殊解 $x_0, y_0$，然后令

$$x = x_0 \times \dfrac{c}{\gcd(a, b)}, y=y_0\times \dfrac{c}{\gcd(a, b)}$$

此时 $ax+by=c$，$x, y$ 可能为负数，有的题目会要求 $x, y$ 为正整数，我们可以把原式进行修改，变为

$$ax+by+k\dfrac{ab}{\gcd(a, b)}-k\dfrac{ab}{\gcd(a, b)}=c \\ a(x+k\dfrac{b}{\gcd(a, b)})+b(y-k\dfrac{a}{\gcd(a, b)})=c$$

即

$$x=x+k\dfrac{b}{\gcd(a, b)} \\ y=y-k\dfrac{a}{\gcd(a, b)}$$

其中 $k \in \Z$。

类欧几里得算法

问题引入

设

$$f(a, b, c, n) = \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor$$

其中 $a, b, c, n$ 是常数，需要 $\mathcal O(\log n)$ 的做法。

若 $a \geq c$ 或 $b \geq c$，我们可以将 $a, b$ 对 $c$ 取模以简化问题。

考虑到

$$x = \left \lfloor \frac{x}{c}\right \rfloor c + x \bmod c$$

故

$$\begin{split} f(a, b, c, n) &= \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor \\ &= \sum_{i=0}^n \left\lfloor\frac{(\left \lfloor \frac{a}{c}\right \rfloor c + a \bmod c)i + (\left \lfloor \frac{b}{c}\right \rfloor c + b \bmod c)}{c}\right\rfloor \\ &= \frac{n(n + 1)}{2} \left \lfloor \frac{a}{c}\right \rfloor + (n+1) \left \lfloor \frac{b}{c}\right \rfloor + f(a \bmod c , b\bmod c, c, n) \end{split}$$

此时一定有 $a < c$ 且 $b < c$。

令

$$S(i)={\left\lfloor\frac{ai + b}{c}\right\rfloor} - 1$$

再进行转化

$$\begin{split} \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor &= \sum_{i=0}^n \sum_{j=0}^{S(i)} 1 \\ &= \sum_{j=0}^{S(n)} \sum_{i=0}^n \left[ j \leq S(i) \right] \end{split}$$

考虑到

$$\begin{split} j \leq S(i) &\iff j + 1 \leq {\left\lfloor\frac{ai + b}{c}\right\rfloor} \\ j + 1 \leq {\left\lfloor\frac{ai + b}{c}\right\rfloor} &\iff j + 1 \leq \frac{ai + b}{c} \\ j + 1 \leq \frac{ai + b}{c} &\iff jc + c \leq ai + b \\ jc + c \leq ai + b &\iff jc + c - b \leq ai \\ jc + c - b \leq ai &\iff jc + c - b - 1 < ai \\ jc + c - b - 1 < ai &\iff \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor < i \end{split}$$

故

$$\begin {split} \sum_{j=0}^{S(n)} \sum_{i=0}^n \left[ j \leq S(i) \right] &= \sum_{j=0}^{S(n)} \sum_{i=0}^n \left[ i > \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor\right] \\ &= \sum_{j=0}^{S(n)} \left(n - \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor \right) \\ &= (S(n) + 1)n - \sum_{j=0}^{S(n)} \left \lfloor \frac{c j + (c - b - 1)}{a} \right \rfloor \\ &= (S(n) + 1)n - f(c, c - b - 1, a, S(n)) \end {split}$$

故

$$f(a, b, c, n) = (S(n) + 1)n - f(c, c - b - 1, a, S(n))$$

可以发现，上述式子是一个递归式，我们不断重复上述过程，先取模，后递归，其实就是辗转相除的过程，时间复杂度 $\mathcal O(\log n)$。

拓展

我们再来推导两个变种求和式

$$g(a, b, c, n) = \sum_{i=0}^n i\left\lfloor\frac{ai + b}{c}\right\rfloor \\ h(a, b, c, n) = \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor^2$$

推导 $g$

引理 1.

$$\sum_{i=0}^n i^2 = \frac{n (n+1)(2n+1)}{6}$$

证明如下：

考虑到

$$(n+1)^3=n^3+3n^2+3n+1$$

故

$$\begin {split} (n+1)^3-n^3 &= 3n^2+3n+1 \\ n^3-(n-1)^3 &=3(n-1)^2+3(n-1)+1 \\ &\cdots \\ 2^3 - 1^3&=3 \times (2-1)^2+3 \times(2-1) + 1 \\ \end {split}$$

将这 $n$ 个等式左右两边相加，得到

$$(n+1)^3 - 1 = 3(1^2 + 2^2 + \cdots +n^2) + 3(1 + 2 + \cdots +n) + n$$

即

$$n^3+3n^2+3n = 3(1^2 + 2^2 + \cdots +n^2) + 3\frac{n(1 +n)}{2} + n$$

整理后得

$$\sum_{i=0}^n i^2 = \frac{n (n+1)(2n+1)}{6}$$

首先和 $f$ 一样，对其取模（根据引理 1. 可将其展开）。

$$\begin{split} g(a, b, c, n) &= \sum_{i=0}^n i\left\lfloor\frac{ai + b}{c}\right\rfloor \\ &= \sum_{i=0}^n i\left\lfloor\frac{(\left \lfloor \frac{a}{c}\right \rfloor c + a \bmod c)i + (\left \lfloor \frac{b}{c}\right \rfloor c + b \bmod c)}{c}\right\rfloor \\ &= \frac{n(n + 1)(2n+1)}{6} \left \lfloor \frac{a}{c}\right \rfloor + \frac{n(n + 1)}{2} \left \lfloor \frac{b}{c}\right \rfloor + g(a \bmod c , b\bmod c, c, n) \end{split}$$

其他部分推导与 $f$ 类似。

令

$$S(i)={\left\lfloor\frac{ai + b}{c}\right\rfloor} - 1$$

有

$$\begin {split} g(a, b, c, n) &= \sum_{j=0}^{S(n)} \sum_{i=0}^n i\left[ j \leq S(i) \right] \\ &= \sum_{j=0}^{S(n)} \sum_{i=0}^n i\left[ i > \left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor\right] \\ \end {split}$$

令

$$t(j) =\left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor$$

则有

$$\begin{split} g(a, b, c, n) &= \sum_{j=0}^{S(n)} \left((t(j) + 1) + (t(j) + 2) + \cdots + n\right) \\ &= \sum_{j=0}^{S(n)} \left(\frac{(t(j)+1+n)\times(n-t(j))}{2}\right) \\ &= \frac{1}{2}\left[ (S(n) + 1) n(n+1) - \sum_{j=0}^{S(n)}(t(j))^2 - \sum_{j=0}^{S(n)}t(j)\right] \\ &= \frac{1}{2} \left[(S(n) + 1) n(n+1) - h(c, c - b - 1, a, S(n)) - f(c, c - b - 1, a, S(n))\right] \end {split}$$

推导 $h$

同样套路，先取模

$$\begin{split} h(a, b, c, n) &= \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor^2 \\ &= \sum_{i=0}^n \left\lfloor\frac{(\left \lfloor \frac{a}{c}\right \rfloor c + a \bmod c)i + (\left \lfloor \frac{b}{c}\right \rfloor c + b \bmod c)}{c}\right\rfloor^2 \\ &= \sum_{i=0}^n \left(i \left \lfloor \frac{a}{c}\right \rfloor + \left \lfloor \frac{b}{c}\right \rfloor + \frac{(a \bmod c)i + (b\bmod c)}{c} \right)^2 \\ \end{split}$$

令

$$Q(i)=\frac{(a \bmod c)i + (b\bmod c)}{c}$$

拆开则有

$$\begin {split} h(a, b, c, n) &= \sum_{i=0}^n \left(i^2\left \lfloor \frac{a}{c}\right \rfloor^2 + \left \lfloor \frac{b}{c}\right \rfloor^2 + (Q(i))^2 + 2 i\left \lfloor \frac{a}{c}\right \rfloor \left \lfloor \frac{b}{c}\right \rfloor + 2i\left \lfloor \frac{a}{c}\right \rfloor Q(i) + 2\left \lfloor \frac{b}{c}\right \rfloor Q(i)\right) \\ &= \frac{n(n + 1)(2n + 1)}{6} \left \lfloor \frac{a}{c}\right \rfloor^2 + (n + 1)\left \lfloor \frac{b}{c}\right \rfloor^2 + h(a \bmod c, b \bmod c, c, n) + \\ & n(n+1)\left \lfloor \frac{a}{c}\right \rfloor \left \lfloor \frac{b}{c}\right \rfloor + 2\left \lfloor \frac{a}{c}\right \rfloor g(a \bmod c, b \bmod c, c, n) + 2\left \lfloor \frac{b}{c}\right \rfloor f(a \bmod c, b \bmod c, c, n) \end {split}$$

照样令

$$S(i)={\left\lfloor\frac{ai + b}{c}\right\rfloor} - 1 \\ t(j) =\left \lfloor \frac{jc + c - b - 1}{a} \right \rfloor$$

发现平方不好处理，将平方转化为加法

$$\begin{split} n^2 &= 2\frac{n(n + 1)}{2} - n \\ &= \left(2\sum_{i=0}^n i \right) - n \end{split}$$

则

$$\begin{split} h(a, b, c, n) &= \sum_{i=0}^n \left\lfloor\frac{ai + b}{c}\right\rfloor^2 \\ &= \sum_{i=0}^n \left[\left(2\sum_{j=0}^{S(i) + 1} j \right) - \left\lfloor\frac{ai + b}{c}\right\rfloor\right] \\ &= \left(2 \sum_{i=0}^n \sum_{j=0} ^{S(i) + 1} j\right) - f(a, b, c, n) \end{split}$$

现在只要解决前面的式子即可

$$\begin{split} \sum_{i=0}^n \sum_{j=0} ^{S(i) + 1} j &= \sum_{i=0}^n \sum_{j=0} ^{S(i)} (j+1) \\ &= \sum_{j=0}^{S(n)} (j+1)\sum_{i=0}^n \left[j < \left\lfloor\frac{ai + b}{c}\right\rfloor \right] \\ &= \sum_{j=0}^{S(n)} (j+1) (n - t(j)) \\ &= (S(n) + 1)n + \frac{S(n)(S(n) + 1)}{2} n - g(c, c - b - 1, a, S(n)) - f(c, c - b - 1, a, S(n)) \end{split}$$

综上所述

$$h(a, b, c, n) =(S(n)+1)(S(n) + 2) n - 2g(c, c - b - 1, a, S(n)) - 2f(c, c - b - 1, a, S(n)) - f(a, b, c, n)$$

当 $a=0$ 时，上述三个式子都可以 $\mathcal O(1)$ 计算。

代码实现时，因为三个函数交错递归，考虑三个一起整体递归，同步计算，时间复杂度为 $\mathcal O(\log n)$。

Luogu5170

#include <bits/stdc++.h>
const int mod = 998244353;
int InvTwo, InvSix;
inline int read()
{
    int cnt = 0; char ch = getchar(); bool op = 1;
    for (; ! isdigit(ch); ch = getchar())
        if (ch == '-') op = 0;
    for (; isdigit(ch); ch = getchar())
        cnt = cnt * 10 + ch - 48;
    return op ? cnt : - cnt;
}

inline int quick_pow(int a, int b)
{
    int ret = 1;
    for (; b; b >>= 1)
    {
        if (b & 1) ret = 1ll * ret * a % mod;
        a = 1ll * a * a % mod;
    }
    return ret % mod;
}

struct eu
{
    int f, g, h;
};

inline int S(int a, int b, int c, int i)
{
    return ((1ll * a * i + b) / c) - 1;
}

inline eu solve(int a, int b, int c, int n)
{
    if (a == 0)
    {
        eu now = {0, 0, 0};
        now.f = 1ll * (b / c) * (n + 1) % mod;
        now.h = 1ll * (b / c) * (b / c) % mod * (n + 1) % mod;
        now.g = 1ll * n * (n + 1) % mod * InvTwo % mod * (b / c) % mod;
        return now;
    }
    if (a >= c || b >= c)
    {
        eu now = {0, 0, 0};
        now.f = (now.f + 1ll * n * (n + 1) % mod * InvTwo % mod * (a / c) % mod) % mod;
        now.f = (now.f + 1ll * (n + 1) * (b / c) % mod) % mod;
        now.g = (now.g + 1ll * n * (n + 1) % mod * (2 * n + 1) % mod * InvSix % mod * (a / c) % mod) % mod;
        now.g = (now.g + 1ll * n * (n + 1) % mod * InvTwo % mod * (b / c) % mod) % mod;
        now.h = (now.h + 1ll * n * (n + 1) % mod * (2 * n + 1) % mod * InvSix % mod * (a / c) % mod * (a / c) % mod) % mod;
        now.h = (now.h + 1ll * (n + 1) * (b / c) % mod * (b / c) % mod) % mod;
        now.h = (now.h + 1ll * n * (n + 1) % mod * (a / c) % mod * (b / c) % mod) % mod;
        eu nxt = solve(a % c, b % c, c, n);
        now.f = (now.f + nxt.f) % mod;
        now.g = (now.g + nxt.g) % mod;
        now.h = (now.h + nxt.h) % mod;
        now.h = (now.h + 1ll * 2 * (a / c) % mod * nxt.g % mod) % mod;
        now.h = (now.h + 1ll * 2 * (b / c) % mod * nxt.f % mod) % mod;
        return now;
    }
    else
    {
        eu now = {0, 0, 0};
        eu nxt = solve(c, c - b - 1, a, S(a, b, c, n));
        now.f = 1ll * (S(a, b, c, n) + 1) * n % mod;
        now.f = (now.f - nxt.f) % mod;
        now.f = (now.f + mod) % mod;
        now.g = 1ll * InvTwo * ((((1ll * (S(a, b, c, n) + 1) * n % mod * (n + 1) % mod - nxt.h) % mod) - nxt.f) % mod) % mod;
        now.g = (now.g + mod) % mod;
        now.h = (1ll * (S(a, b, c, n) + 1) * (S(a, b, c, n) + 2) % mod * n % mod) % mod;
        now.h = (now.h - 1ll * 2 * nxt.g) % mod;
        now.h = (now.h - 1ll * 2 * nxt.f) % mod;
        now.h = (now.h - now.f) % mod;
        now.h = (now.h + mod) % mod;
        return now;
    }
}

int main()
{
    int t = read();
    InvTwo = quick_pow(2, mod - 2);
    InvSix = quick_pow(6, mod - 2);
    while (t --)
    {
        int n, a, b, c;
        n = read(), a = read(), b = read(), c = read();
        eu ans = solve(a, b, c, n);
        printf("%d %d %d\n", ans.f, ans.h, ans.g);
    }
    return 0;
}

Reference

oi-wiki

中国剩余定理（CRT）

问题引入

我们需要解决满足下列方程组的一个解 $x$，并且保证所有的 $b$ 都互相互质。

$$\begin{cases}x \equiv a_1 \ (\bmod \ b_1) \\ x \equiv a_2 \ (\bmod \ b_2) \\ ~~~~~~~~~~~~~\vdots \\ x \equiv a_n \ (\bmod \ b_n)\end{cases}$$

分析

我们考虑设 $M=\prod\limits_{i=1}^{n} b_i$，$m_i=\dfrac{M}{b_i}$，$t_i$ 为 $m_ix\equiv 1\pmod {b_i}$ 的一个解，则解为：

$$x=\sum\limits_{i=1}^{n}a_im_it_i$$

我们现在来证明一下 $x$ 是一个满足条件的解。

证明：

对于任意的 $k \in [1, n]$，$a_km_kt_k \equiv a_k \pmod {b_k}$，

而 $a_km_kt_k \equiv 0 \pmod {b_i}(i \in [1, k - 1] \ \bigcup \ [k + 1, n])$，故 $x$ 是一个满足条件的解。

我们知道，$M \equiv 0 \pmod {b_i} (i \in [1, n])$，故 $x + kM \ (k \in \Z)$ 都是满足答案的解，若题目要求最小正整数解，则：

$$x = (x \bmod M + M) \bmod M$$

#include <bits/stdc++.h>
const int N = 10 + 10;
using namespace std;

int n;
long long M = 1, m[N], t[N], ans;
int a[N], b[N];

inline void exgcd(long long a, long long b, long long &x, long long &y)
{
    if (b == 0)
    {
        x = 1, y = 0;
        return;
    }
    exgcd(b, a % b, x, y);
    long long d = x; x = y; y = d - (a / b) * y;
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; ++ i)
        cin >> a[i] >> b[i], M = M * a[i];

    for (int i = 1; i <= n; ++ i)
    {
        long long y;
        m[i] = M / a[i];
        exgcd(m[i], a[i], t[i], y);
        ans = (ans + m[i] * t[i] * b[i]) % M;
    }

    cout << (ans + M) % M << endl;

    return 0;
}

拓展中国剩余定理 （EXCRT）

问题引入

我们要解决的还是这么上文那个问题，其中不保证所有的 $b$ 都互相互质。

分析

我们假设已经有一个解 $x$ 满足前 $k - 1$ 个同余方程，此时我们设 $d=\operatorname{lcm}\{b_1, b_2, \cdots, b_{k-1}\}$，显然，$x + jd \ (j \in \Z)$ 都是满足条件的解。

那么对于第 $k$ 个同余方程，即要满足 $x+jd\equiv a_k \pmod {b_k}$，可转换为 $jd \equiv a_k - x \ (\bmod \ b_{k})$，显然，可以用扩展欧几里得来求解 $j$，若 $\gcd(d, b_k)\nmid a_k-x$ 那么就无解，否则 $x = x+jd, d=\operatorname{lcm}\{b_1, b_2, \cdots, b_k\}$，此时的 $x$ 是满足前 $k$ 个方程的，这样一直做到最后就可以求出答案。

#include <bits/stdc++.h>
const int N = 1e5 + 10;
using namespace std;

int n; long long X, m;
long long a[N], b[N];

inline long long exgcd(long long a, long long b, long long &x, long long &y)
{
    if (b == 0)
    {
        x = 1, y = 0;
        return a;
    }
    long long d = exgcd(b, a % b, x, y);
    long long t = x; x = y; y = t - (a / b) * y;
    return d;
}

inline long long quick_mul(long long a, long long b, long long p)
{
    a %= p, b %= p;
    long long res = 0, op = 1;
    if (b < 0) b = -b, op = -1;
    for (; b; b >>= 1)
    {
        if (b & 1) res = (res + a) % p;
        a = (a + a) % p;
    }
    return res * op;
}

int main()
{
    cin >> n;
    for (int i = 1; i <= n; ++i)
        cin >> a[i] >> b[i];
    X = b[1], m = a[1];

    for (int i = 2; i <= n; ++i)
    {
        long long d, x, y, c = (b[i] - X);
        d = exgcd(m, a[i], x, y);
        x = quick_mul(x, (c / d), (a[i] / d));
        x = (x % (a[i] / d) + (a[i] / d)) % (a[i] / d);
        X = (X + x * m);
        m = m / d * a[i];
        X = (X % m + m) % m;
    }

    cout << X << endl;

    return 0;
}

BSGS

问题引入

我们需要解决 $a^x \equiv b \pmod M$ 的一个解 $x$，其中 $a, M$ 互质。

分析

首先，如果有解，那么解一定在 $[0, M]$ 这个区间内，根据抽屉原理，$a^x \bmod M$ 只会有 $M$ 个取值，故 $[0, M]$ 中至少会存在两个相同余数的数，即一定存在一对 $i, j$ 满足 $a^i\equiv a^j \pmod M$，那么在 $i$ 之后，$a^i \bmod \ M, a^{i + 1} \bmod \ M, \cdots a^{j} \bmod \ M$ 这些数会循环出现，故如果有解，那么解一定在 $[0, M]$ 这个区间内。

我们考虑把 $x$ 进行拆分，$x=qm-p$，其中 $m=\left\lfloor\sqrt{x}\right\rfloor+1$，$q \in [1, m]$，$p \in [0, m]$。

显然，$qm-p$ 可以表示出 $0 \sim m$ 中的所有数。

我们把原式修改一下得到 $a^{qm-p} \equiv b \pmod M$，$a^{qm} \equiv b \times a^{p} \pmod M$（该步后往前推导需要有乘法逆元）

我们发现 $q$ 与 $p$ 没有关系了，考虑先 $\mathcal O(\sqrt{x})$ 枚举 $p$，把 $(b \times a^p) \bmod M$ 放入哈希表（map 也行），然后再 $\mathcal O(\sqrt{x})$ 枚举 $a^{qm}$，找到最小的 $q$ 满足 $a^{qm} \bmod M$ 在哈希表中出现过，那么 $qm+p$ 就是最小的答案。

#include <bits/stdc++.h>

int p, b, n, m;
std::map <int, int> mp;

inline int quick_pow(int x, long long y)
{
    int res = 1; x %= p;
    for (; y; y >>= 1)
    {
        if (y & 1) res = 1ll * res * x % p;
        x = 1ll * x * x % p;
    }
    return res;
}

int main()
{
    std::cin >> p >> b >> n;

    m = sqrt(p) + 1;
    for (int i = 0; i <= m; ++ i)
        mp[1ll * n * quick_pow(b, i) % p] = i;
    for (int i = 1; i <= m; ++ i)
    {
        if (mp[quick_pow(b, 1ll * i * m)])
        {
            std::cout << 1ll * i * m - mp[quick_pow(b, i * m)] << std::endl;
            return 0;
        }
    }
    std::cout << "no solution" << std::endl;

    return 0;
}

扩展 BSGS

问题引入

我们要解决的还是 $a^x\equiv b \pmod p$ 这个问题，只是不满足 $a$ 与 $p$ 互质。

分析

若 $a$ 与 $p$ 互质，那么就会存在逆元，可以将 $x$ 分解为 $qm - p$ 来求解，所以我们要想办法让 $a, p$ 变得互质。

我们考虑把这个式子转化为不定方程，$a^x -yp = b$，首先，我们求出 $d_1=\gcd(a, p)$，然后整个式子都除以 $d_1$，若 $b$ 与 $d_1$ 互质，那么就无解，否则我们继续求出 $d_2=\gcd(a, p)$，一直进行下去，直到 $a$ 与 $p$ 互质。

假设进行了 $k$ 次，最后得到式子就是：

$$\dfrac{a^k}{\prod\limits_{i=1}^{k}d_i}\times a^{x-k}-\dfrac{yp}{\prod\limits_{i=1}^{k}d_i}=\dfrac{b}{\prod\limits_{i=1}^{k}d_i}$$

其中 $a$ 与 $p$ 互质。

我们再化为同余方程的形式：

$$\frac{a^k}{\prod\limits_{i=1}^{k}d_i}\times a^{x-k}\equiv \frac{b}{\prod\limits_{i=1}^{k}d_i} \left( \bmod {\frac{p}{\prod\limits_{i=1}^{k}d_i}}\right)$$

显然，求出

$$\dfrac{a^k}{\prod\limits_{i=1}^{k}d_i}$$

的逆元，然后直接移到右边，直接用 BSGS 来解决，然后最后把答案加上 $k$ 即可。

注意，不保证答案小于 $k$ 的情况，故在消因子前做一下 $O(k)$ 枚举，直接验证 $a^i \equiv b \pmod p$，就能避免这种情况。

code

高斯消元

问题引入

给出 $n$ 组方程，每组方程形如 $\sum\limits_{i=1}^n a_ix_i = b_i$，要求求出 $x_1, x_2 \cdots x_n$ 或告知无解。

分析

我们把这些方程转化为矩阵

$$\begin{vmatrix} & a_{1,1} & a_{1,2} & \cdots & a_{1, n} & b_1 \\ & a_{2,1} & a_{1,2} & \cdots & a_{2, n} & b_2 \\ & \vdots & \vdots & \ddots & \vdots & \vdots & \\ & a_{n,1} & a_{n,2} & \cdots & a_{n, n} & b_n \\ \end{vmatrix}$$

我们希望最后的矩阵形如

$$\begin{vmatrix} & 1 & 0 & 0 & 0 & \cdots & 0 & c_1 \\ & 0 & 1 & 0 & 0 & \cdots & 0 & c_2 \\ & 0 & 0 & 1 & 0 & \cdots & 0 & c_3 \\ & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \\ & 0 & 0 & 0 & \cdots & 1 & 0 & c_{n-1} \\ & 0 & 0 & 0 & \cdots & 0 & 1 & c_n \\ \end{vmatrix}$$

对于第 $i$ 行就能求出 $x_i = c_i$，即可求出解。

流程如下：

1. 我们从上往下做，先将当前系数绝对值最大的移到当前这一行，减小误差。
2. 让当前系数化为 $1$，即让当前行的每一个数都除以当前系数。
3. 再用当前第 $i$ 行的第 $i$ 个系数（即为 $1$）消去 $j \ (j \not= i)$ 的第 $i$ 个系数，把它消为 $0$。

重复做以上的操作，直到做到第 $n$ 行，此时最后的矩阵形如上述。

我们可以在 $\mathcal O(n^3)$ 的时间复杂度完成。

判断解情况：

最后说一下如何判无解和无限解：

1. 无解：当前行系数全为 $0$，但值不为 $0$。

2. 无限解：当前行系数全为 $0$，且值为 $0$。

一、积性函数

1. 积性函数：

对于任意互质的整数 $a$ 和 $b \ (\gcd(a,b)=1)$ 有性质 $f(ab)=f(a)\times f(b)$。

2. 完全积性函数：

对于任意的整数 $a$ 和 $b$ 有性质 $f(ab) = f(a) \times f(b)$。

3. 证明 $\varphi(n)$ 是积性函数。

1. 定理：

当 $\gcd(x, y)=1$ 时，$\varphi(x \times y)=\varphi(x) \times \varphi(y)$。

2. 证明：

我们将 $1 \sim x \times y$ 排列成如下方阵：

$$\begin{matrix} 1 & 2 &\cdots & i & \cdots & y \\ 1 + y & 2 + y & \cdots & i + y & \cdots & 2 \times y \\ 1 + 2 \times y & 2 + 2 \times y & \cdots & i + 2 \times y & \cdots & 3 \times y \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 1 + j \times y & 2 + j \times y & \cdots & i+j \times y & \cdots & (j + 1) \times y \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 1 + (x-1) \times y & 2 + (x-1) \times y & \cdots & i + (x-1) \times y & \cdots & x \times y \\ \end{matrix}$$

由于 $\gcd(x, y)=1$，每一列 $i + j \times y$ 都能构成一个 $\bmod y$ 的完全剩余系 $\{\overline0, \overline1, \cdots, \overline{y-2}， \overline{y-1}\}$，则每一列与 $y$ 互质的数就有 $\varphi(y)$。同理，对于每一行，都能构成一个 $\bmod x$ 的完全剩余系 $\{\overline0, \overline1, \cdots, \overline{x-2}， \overline{x-1}\}$，并且它们都是对齐的，则每一行与 $x$ 互质的数就有 $\varphi(x)$。

若 $a$ 与 $x \times y$ 互质，则 $a$ 与 $x$ 互质且 $a$ 与 $y$ 互质，即 $\varphi(x \times y)=\varphi(x) \times(y)$。

二、欧拉函数

1. 欧拉函数的定义及求法

$\varphi(x)$ 表示有多少个 $i \ (1 \leq i \leq x)$ 满足 $\gcd(i, x) = 1$（即有多少个 $\leq x$ 的 $i$ 与 $x$ 互质）。

显然当 $x$ 为质数时，$\varphi(x) = x - 1$。

把 $x$ 唯一分解：

$$x=\prod\limits_{i=1}^{m}p_i^{c_i}$$

有

$$\varphi(x)=x \times \prod\limits_{i=1}^{m}\dfrac{p_i-1}{p_i}$$

我们来证明一下这个定理：

引理 1：设 $x=p^k \ (p$ 为质数$)$，那么 $\varphi(x)=p^{k-1}\times(p-1)$。

证明 1：此时对于一个 $y$ 与 $x$ 互质，当且仅当 $y \not\equiv 0 \pmod p$，我们考虑将 $x$ 分成许多段，长度为 $p$，有 $\frac{x}{p}=\frac{p^k}{p}=p^{k-1}$ 段。其中每一段都有 $p-1$ 个数与 $p$ 互质，有 $p^{k-1}$ 段，故 $\varphi(x)=p^{k-1}\times(p-1)$。

根据欧拉函数的积性性质，我们可以得到：

$$\varphi(x)=\prod\limits_{i=1}^{m}\varphi(p_i^{c_i})=\prod\limits_{i=1}^mp^{c_i-1}\times(p_i-1)=\prod\limits_{i=1}^m\dfrac{p^{c_i}}{p}\times(p_i-1)=\prod\limits_{i=1}^mp^{c_i}\times\dfrac{p_i-1}{p_i}=x\times\prod\limits_{i=1}^m\dfrac{p_i-1}{p_i}$$

故我们可以使用试除法在 $O(\sqrt{x})$ 的时间内求出 $\varphi(x)$。

2. 欧拉函数的一些性质

定理 1：$\forall \ n > 1$，$1 \sim n$ 之间与 $n$ 互质的数的和为：$\dfrac{n \times \varphi(n)}{2}$。

证明：因为 $\gcd(n, x)=\gcd(n, n-x)$，故与 $n$ 不互质的数是成对出现的，且平均值为 $\dfrac{n}{2}$，又考虑到共有 $\varphi(n)$ 个与其互质的数，故 $\forall \ n > 1$，$1 \sim n$ 之间与 $n$ 互质的数的和为 $\dfrac{n \times \varphi(n)}{2}$。

定理 2：当 $x$ 为奇数时，$\varphi(2 \times x)=\varphi(x)$。

证明 2：把 $x$ 唯一分解：

（其中 $p_i \not=2$），$2\times x=2\times\prod\limits_{i=1}^mp_i^{c_i}$，根据上面的欧拉函数的求法可得：

$$\varphi(2\times x)=2 \times x \times \dfrac{2-1}{2} \times \prod\limits_{i=1}^m\dfrac{p_i - 1}{p_i}=x\times \prod\limits_{i=1}^m\dfrac{p_i-1}{p_i} \\ \varphi(x)=x\times \prod\limits_{i=1}^m\dfrac{p_i-1}{p_i}$$

所以当 $x$ 为奇数时，$\varphi(2\times x)=\varphi(x)$。

定理 3：若 $x > 2$ ，那么 $\varphi(x)$ 是偶数。

证明 3：我们考虑把 $x$ 唯一分解：

$$x=\prod\limits_{i=1}^{m}p_{i}^{c_i} \\ \varphi(x)=x\times\prod\limits_{i=1}^m\dfrac{p_i-1}{p_i}$$

显然，除开 $2$ 这个质数以外，其他质数都是奇数，减去 $1$ 之后就都是偶数了，故不包含 $2$ 这个质因子的数的欧拉函数一定是偶数，而 $2-1=1$，对结果没有影响，并且如果这个数有 $2$ 这个质因子，那么它一定是偶数，所以若 $x > 2$ ，那么 $\varphi(x)$ 是偶数。

定理 4：设 $p$ 为质数，若 $p | n$ 且 $p^2 | n$ 则 $\varphi(n)=\varphi(\dfrac{n}{p})\times p$。

证明 4：由于 $p^2 | n$ 所以 $\dfrac{n}{p}$ 也有 $p$ 这个因子，故 $\varphi(n)$ 和 $\varphi(\dfrac{n}{p})$ 的区别只有 $n$ 和 $\dfrac{n}{p}$ 这个乘积的差别，故可得 $\varphi(n)=\varphi(\dfrac{n}{p})\times p$。

定理 5：设 $p$ 为质数，若 $p | n$ 且 $p^2 \nmid n$ 则 $\varphi(n)=\varphi(\dfrac{n}{p})\times(p-1)$。

证明 5：由于 $p^2 \nmid n$ 所以 $\dfrac{n}{p}$ 与 $p$ 互质，故 $\varphi(n)$ 和 $\varphi(\dfrac{n}{p})$ 的比，多了 $\dfrac{p-1}{p}$ 和 $p$ 这两个乘积，故可得：$\varphi(n)=\varphi(\dfrac{n}{p})\times (p-1)$。

定理 6：

$$x=\sum\limits_{d|x}\varphi(d)$$

证明：设 $f(n)=\sum\limits_{d|n}\varphi(d)$， 根据乘法分配律以及 $\varphi$ 的积性性质，我们可以得到：若 $n, m$ 互质，则

$$f(nm) = \sum\limits_{d | nm} \varphi(d)=\left(\sum\limits_{d|n}\varphi(d)\right)\times \left(\sum\limits_{d|m}\varphi(d) \right)=f(n)\times f(m)$$

故 $f(n)$ 是积性函数。对于单个质因子 ：

$$f(p_i^{c_i})=\sum\limits_{d|p_i^{c_i}}\varphi(d)=\varphi(1)+\varphi(p)+\cdots+\varphi(p_i^{c_i})$$

我们发现这是一个等比数列求和再加 $1$，结果为 $p^{c_i}$。所以

$$f(n)=\prod\limits_{i=1}^mf(p_i^{c_i})=\prod\limits_{i=1}^mp_i^{c_i}=n$$

拉格朗日插值

多项式定义

$n$ 次多项式定义为：

$$f(x) = a_nx^n + a_{n-1}x^{n-1} + a_{n-2}x^{n-2} + \cdots + a_0$$

问题引入 1

给出 $n + 1$ 个点 $(x_i, y_i)$，表示 $f(x_i)=y_i$，要求求出 $f(t)$。

做法 1：

高斯消元。

我们把这 $n + 1$ 个点都写出来：

$$\left\{ \begin{aligned} a_nx_1^n + a_{n-1}x_1^{n-1} + a_{n-2}x_1^{n-2}+\cdots+a_0 = y_1 ~~~~ \\ a_nx_2^n + a_{n-1}x_2^{n-1} + a_{n-2}x_2^{n-2}+\cdots+a_0 = y_2 ~~~~ \\ \vdots ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ \\ a_nx_{n+1}^n + a_{n-1}x_{n+1}^{n-1} + a_{n-2}x_{n+1}^{n-2}+\cdots+a_0 = y_{n+1} \\ \end{aligned} \right.$$

把 $a_n, a_{n-1}, \cdots, a_0$ 看作未知量，把 $x,y$ 看作已知量，这样就可以高斯消元进行求解，解出 $a_n,a_{n-1},\cdots,a_0$。

时间复杂度：$O(n^3)$。

做法 2：

拉格朗日插值。

我们从二次多次项入手，我们给出 $3$ 个点，$(x_1, y_1), (x_2, y_2), (x_3,y_3)$。

可以列出方程组：

$$\left\{ \begin{aligned} a_2x_1^2 + a_1x_1^1+a_0 = y_1 \\ a_2x_2^2 + a_1x_2^1+a_0 = y_2 \\ a_2x_3^2 + a_1x_3^1+a_0 = y_3 \\ \end{aligned} \right.$$

这里我们换一种思路，我们构造三个函数，然后相加得到我们要的函数。

构造函数：

$$f_1(x)=\left\{ \begin{aligned} 1 && (x=x_1) \\ 0 && (x\not=x_1) \\ \end{aligned} \right. , f_2(x)=\left\{ \begin{aligned} 1 && (x=x_2) \\ 0 && (x\not=x_2) \\ \end{aligned} \right. ,f_3(x)=\left\{ \begin{aligned} 1 && (x=x_3) \\ 0 && (x\not=x_3) \\ \end{aligned} \right.$$

那么当 $f(x)=y_1f_1(x)+y_2f_2(x)+y_3f_3(x)$，此时就满足 $x$ 的取值为 $x_1,x_2,x_3$ 时，$f(x)$ 的取值就为 $y_1,y_2,y_3$。

现在考虑 $n$ 次多项式，给出 $n + 1$ 个点：$(x_1, y_1), (x_2, y_2),\cdots,(x_{n + 1}, y_{n + 1})$。

我们考虑构造 $f_1 \sim f_{n+1}$。

设

$$f_i(k)=\prod\limits_{j\not=i}\dfrac{k-x_j}{x_i-x_j}$$

则 $f_i(k)$ 满足在 $k=x_i$ 时函数值为 $1$，在 $k=x_j(j \not= i)$ 时函数值为 $0$。

那么多项式

$$f(k)=\sum\limits_{i=1}^{n+1}y_if_i(k)=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\not=i}\dfrac{k-x_j}{x_i-x_j}$$

就满足当 $k$ 为 $x_1, x_2, \cdots, x_{n+1}$ 时 $f(k)$ 的取值就为 $y_1, y_2, \cdots, y_{n + 1}$。

如果我们只要求解特定函数值，我们可以使用逆元来代替除法。

时间复杂度：$\mathcal O(n^2)$。

优化 1

当 $x_a$ 与 $a$ 有关系时，例如 $x_a=a$ 时，可以使用前缀乘来优化，时间复杂度：$\mathcal O(n)$。

优化 2

重心拉格朗日插值法

当我们要解决两个操作：

1. 往点集里加一个点 $(x, y)$
2. 给出一个 $k$，求出当前点集对应的多项式 $f(k)$ 的值。

考虑把原式子拿去转换。

$$f(k)=\sum\limits_{i=1}^{n+1}y_i\prod\limits_{j\not=i}\dfrac{k-x_j}{x_i-x_j} = \sum\limits_{i=1}^{n + 1} \dfrac{y_i \times \prod\limits_{j\not=i} (k-x_j)}{\prod\limits_{j\not=i} (x_i-x_j)}=\sum\limits_{i=1}^{n+1}\dfrac{y_i \prod\limits_{j=1}^{n + 1}(k-x_j)}{(k-x_i) \prod\limits_{j \not=i}(x_i-x_j)}$$

记

$$g(k) = \prod\limits_{i=1}^{n + 1}(k-x_i)$$

则

$$f(k) = g(k)\times \sum\limits_{i=1}^{n+1}\dfrac{y_i }{(k-x_i) \prod\limits_{j \not=i}(x_i-x_j)}$$

再记

$$t(i) = \dfrac{y_i}{\prod\limits_{j \not= i} (x_i - x_j)}$$

则

$$f(k)=g(k) \times \sum\limits_{i=1}^{n+1}\dfrac{t(i)}{k-x_i}$$

新加点时 $\mathcal O(n \log n)$ 求出 $t(i)$，故每加一个点就能做到 $\mathcal O(n \log n)$ 解决。

一些常识

1. $\sum\limits_{i=1}^n i^k$ 是一个 $k + 1$ 次的多项式，可以插值。