~~数论小白，在线挨打~~

很难想象笔者在 48h 之前对这些艰涩的数论名词还是一无所知。然而在两天的探（zhua）索（ba）之后，我也不禁感叹数论的美（e）妙（xin）。

扩 展 欧 几 里 得($\mathtt{exgcd}$)

诸位读者一定对 $\gcd$ 毫不陌生，这是一种递归（或者不递归？）求最大公约数的方法，其核心在于 $\gcd(a,b)=\gcd(b, a\bmod b)$ 而 $\mathtt{exgcd}$ 却与之略有提升，它还是用于求解同余方程 $ax+by=\gcd(a,b)$ 的一种工具。

实现？ ：

  1. 对于我们手上的一个方程，造一个拥有两个参数 $a$ ,$b$ 的函数 $\mathtt{exgcd}$ 层层调用，老千层饼了，递归时参数的改动与 $\gcd$ 相似，都是以 $b$ 成为新 $a$ ，以 $a \bmod b$ 成为新 $b$ 。
  2. 在每一层的递归结束后，还不能结束，我们有一对全局变量 $x,y$ ，它们的存在意义是作为这个方程的最终解。将 $x$ 更改成原来的 $y$ ，将 $y$ 更改成 $x- \left\lfloor\dfrac{a}{b}\right\rfloor * y $
（附）有关 2 中 $x,y$ 变量变化的原因？ ：

不妨设新的 $x$ 为 $ x’$ ，新的 $y$ 为 $ y’$ ；

则有 $ax+by=bx’+(a\bmod b)y’$ ；

也就是 $ax+by=bx’+(a-b * \left\lfloor\dfrac{a}{b}\right\rfloor)y’$ ；

一波运算即得 $ax+by=ay’+b(x’-\left\lfloor\dfrac{a}{b}\right\rfloor * y’) $；

所以 $ x’$ 可以转化成 $y$ ， $ y’$ 可转化为 $x-\left\lfloor\dfrac{a}{b}\right\rfloor * y$。 $\Box$

  1. 我们继续。递归是一定要停下来的，因为只要道路不断延伸…………（雾）。我们要为他设定一个终点，当 $b=0$ 的时候，此时如果再做下去，就会出现 $a \div 0$ 的 RE 情况，观察此时我们手上的方程式 $ax+0 y = \gcd (a,b)$ ，*一定有一个极小解 $x=1,y=0$ ，走到这一步，我们马上改变 $x,y$ ，然后 return 走人。
这里放上我们 $\mathtt{exgcd}$ 的函数代码～：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
int x, y;
void exgcd(int a, int b)
{
     if(b==0)
     {
         x = 1;
         y = 0;
         return;
     }
     exgcd(b, a % b);//递归
     int tmp = x;
     x = y;
     y = tmp - (a / b) * y;//更改x,y的值
}

乘 法 逆 元($\mathtt{inv}$)：

若有 $ ax \equiv1 \bmod b$（ $\equiv $ ，即同余）（要求 $a,b$ 互质) ， 则称 $a$ 关于 1 模 $b$ 的乘法逆元为 $x$ 。也可表示为 $ax≡1 \pmod b$ 。对于一个数关于 1 模 $b$ 的乘法逆元，我们还有另一个通俗的名字“倒数”。

实现？ :

这里给出两种求乘法逆元的方法，分别适用不同的情况：

$\mathtt{exgcd}$法：

适用情形：不管模数是不是质数，较适用于对单个数的查询（言外之意使用面更广，但复杂度不够优）

前文提到，$\mathtt{exgcd}$ 其作用是求解同余方程，它形如 $ax+by=\gcd(a,b)$ ，可我们怎么把求乘法逆元的式子转化成同余方程的形式呢？我们考虑稍加变换。观察这个式子：

我们可以将同余化为等号，用我们的转化带余除法的传统艺能：

再稍加移项～：

啊这，这不是同余方程？那我们对其求解，其答案中的 $x$ 不就是一个满足这个方程的 $a$ 关于 1 模 $f$ 的乘法逆元嘛？

芜湖，起飞！上代码：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
void inv(int a, int b)
{
    if (!b)
    {
        x = 1;
        y = 0;
        return;
    }
    exgcd(b, a % b);
    int tmp = x;
    x = y;
    y = tmp - (a / b) * y;
    return ;
}
int main()
{
    int n, k;
    cin >> n  >> k;
    for (int i = 1; i <= n;i++)
    {
        x = y = 0;
      	inv(i, k);
        cout << (x%k+k)%k << endl;
    }
}

递推法：

适用情形：模数必须是质数，使用于需要求多个数的乘法逆元的情况（换言之就是限制多但是效率快，可以达到线性）

所谓递推，公式求出来，代码就是几行的事，重在公式的推导。来看原式：

我们设这个 $b=k * i+r$ ，也就是设 $k=\left\lfloor\dfrac{b}{i}\right\rfloor$ ， $r=b \bmod i$ 。

左右同乘 $\dfrac{1}{i * r}$ 得：

稍加移项即得：

也就是：

啊这，左边这玩意不就是 $i$ 关于 $1 \bmod b$ 的逆元嘛？有被爽到。而右边那些玩意完全都是已知的。诸位可能疑问的，估计就是最后那个 $(b\bmod i)^{-1}$ 为什么已知，然而 $b\bmod i$ 必然小于 $i$，我们的递推又是从小到大推，所以……

女少 口阿！上代码：
1
2
3
4
5
6
7
8
9
10
11
12
signed main()
{
        int n,p;
        cin >> n >> p;
        inv[1]=1;//inv即为逆元之意
        cout<<1<<endl;
        for (int i = 2; i <= n;i++)
        {
                inv[i] = p - (p / i) * inv[p % i]%p;//递推公式
                printf("%d\n", inv[i]);
        }
}

矩 阵 乘 法

矩阵？乘法？

定义：我们手上有两个矩阵如下:

我们矩阵乘法的规则是“对位相乘”，积矩阵第 $i$ 行第 $j$ 列的值等于第一个乘数矩阵第 $i$ 行第 $k$ 列与第二个乘数矩阵第 $k$ 行第 $j$ 列的乘积的累加（其中$1\le k \le n$）

所以上边的运算结果就是$\begin{pmatrix}a1 a1+a2 b3&a1 b2+a2 b4\\a3 b1+a4 b3&a3 b2+a4 b4\end{pmatrix}$ ！！

学这东西有什么用呐？

一个非常重要的作用就是优化（加速）递推，一般能将线性的递推优化到 $\log(n)$级别，极度实用。

要发挥这一作用，我们可以定义一种名为 $matrix$ 的结构体，里边存的是个二维数组 $num$ ，这就是矩阵里的值了。我们可以为矩阵乘法重载运算符（笔者觉得很香），代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
matrix operator*(matrix x)
        {
                matrix c;
                for (int i = 1; i <= n; i++)
                {
                        for (int j = 1; j <= n; j++)
                        {
                                c.num[i][j] = 0;
                                for (int k = 1; k <= n; k++)
                                {
                                        c.num[i][j] = (c.num[i][j] + (num[i][k] * x.num[k][j]) % MOD) % MOD;
                                }
                        }
                }
                return c;
        }

不妨参照上方定义略加理解，当板子背亦可。

但只有子弹没有枪是不行的，我们需要一个高效的方法帮助递推转移。矩阵优化的基本方略是每次将原来的状态存入矩阵，我们称之为 $ORZ$ ，每次转移时，将 $ORZ$ 乘上一个转移矩阵，我们称之为 $RBQ$ （因为他要被用好多次），使新得的矩阵成为下一个状态，以起到转移状态的目的。如下例是一个快速求斐波那契的矩阵加速：

$\begin{pmatrix}f(n-1)&f(n-2)\\0&0\end{pmatrix}$ 这是 $ORZ$ ，初始矩阵。

$\begin{pmatrix}1&1\\1&0\end{pmatrix}$ 这是 $RBQ$ ，转移矩阵。

$\begin{pmatrix}f(n-1)+f(n-2)&f(n-1)\\0&0\end{pmatrix}$ 这是积，也能作为下一个状态。

因为矩阵乘法是满足结合律的，原因我也不懂 （逃 ，所以如果我们要多次乘以这一矩阵，就可以用快速幂优化累乘。矩阵的快速幂如下代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
matrix quick_pow(matrix b, int k)
{
        if (k == 1)
                return b;
        matrix tmp = quick_pow(b, k / 2);
	//笔者尤其喜欢递归快速幂的写法，毕竟出题人不会恶心到要让人爆栈（？）
        if (k % 2 == 0)
                return tmp * tmp;
        else
        {
                return tmp * tmp * b;
        }
        //因为我们重载了运算符，所以快速幂变得极为方便
}

可是，斐波那契毕竟是数论弃子，递推中的弟弟，要是状态一复杂，推导转移方法将会极考验思维。

笔者在此分享自己推导转移矩阵的方法：

假设我们已知有一种函数 $A(x)$ ，符合 $A(0) = 1 , A(1) = 1 , A(N) = X A(N - 1) + Y A(N - 2)(N \ge 2)$ ，又有函数 $S(x)$ ，符合 $S(N) = A(0)^2 +A(1)^2+……+A(n)^2$ ，给定一个数 $n$，要求 $S(n)$ 。

  1. 先将递推过程中需要的量写在初始矩阵 $ORZ$ 的第一行（依据个人喜好和方便确定是那些量），如下
  2. 对于每一项，按照递推式决定这一项在转移后需要进行怎样的变换，需要哪些“原材料”（这种变换要求只使用我们写下的那行数里的元素）。我们以第一项 $S(n-1)$ 为例，对于 $S(n-1)$ ，他要成为 $S(n)$ ，变化方法即：需要1个 $S(n-1)$ ，$X^2$ 个 $A(N-1)^2$ ， $Y^2$ 个 $A(N-2)^2$ ， $2XY$ 个 $A(n-1)A(n-2)$。

后面几项都可如此推出。

  1. 将每一项的变换方法得知后，对于第 $i$ 项，将每一个元素在造下个状态时所需的个数从上向下排在转移矩阵 $RBQ$ 的第 $i$ 列。为什么可以对于一个元素，只讨论他自己是怎么转移，而不用考虑其他项呢？原因显然，毕竟第二个乘数 $RBQ$ 的各列之间互不影响，这一引论可以手推/肉眼观察法得出。对于那些矩阵里空着的位置，我们用 0 补全。得出的矩阵就像这样：
  1. 来看我们的 $ORZ$ ，刚才第 1 行我们仅仅只是填上了几个字母，但他里面的初始状态究竟是什么，还得咱自己去初始化（具体是什么值题干会给），比如这一题，他的初始矩阵就是:
刚才提到的那道递推题的高清无码代码.avi：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
#include <bits/stdc++.h>
#define int long long
const int M = 10007;
using namespace std;
int n = 4;
struct matrix
{
    int num[17][17];
    matrix operator * (const matrix x) const
    {
        matrix c;
        memset(c.num, 0, sizeof(c.num));
        for (int i = 1; i <= n; i++)
        {
            for (int j = 1; j <= n; j++)
            {
                for (int k = 1; k <= n; k++)
                {
                    c.num[i][j] = (c.num[i][j] + (num[i][k] * x.num[k][j]) % M) % M;
                }
            }
        }
        return c;
    }//重载运算符*
} ORZ, RBQ;
matrix quick_pow(matrix x, int p)
{
    if (p == 1)
        return x;
    matrix tmp = quick_pow(x, p / 2);
    if (p % 2 == 0)
        return tmp * tmp;
    else
        return tmp * tmp * x;
}//矩阵快速幂
signed main()
{
    int N, x, y;
    while (~scanf("%lld%lld%lld", &N, &x, &y))
    {
        memset(RBQ.num, 0, sizeof(RBQ.num));
        ORZ.num[1][1] = 2;
        ORZ.num[1][2] = 1;
        ORZ.num[1][3] = 1;
        ORZ.num[1][4] = 1;
        RBQ.num[1][1] = 1;
        RBQ.num[2][1] = x * x%M;
        RBQ.num[3][1] = y * y%M;
        RBQ.num[4][1] = 2 * x * y%M;
        RBQ.num[2][2] = x * x%M;
        RBQ.num[3][2] = y * y%M;
        RBQ.num[4][2] = 2 * x * y%M;
        RBQ.num[2][3] = 1;
        RBQ.num[2][4] = x%M;
        RBQ.num[4][4] = y%M;
        //初始化，有被爽到                             
        matrix ans = ORZ * quick_pow(RBQ, N-1);
        cout << ans.num[1][1] << endl;
    }
}

欧 拉 函 数（ $\varphi$ ）

定义：

我们定义 $\varphi(d)$ 代表着 $d$ 以内的与 $d$ 互质的数的个数，例如 $\varphi(6)=2$ ，因为 6 之内有 1,5 两个数与它互质。特别的，我们有 $\varphi(1)=1$ 。

这玩意有好多好多奇奇怪怪的性质：

性质0：

对于 $n$ 的每个因数 $i$ ，所有 $\varphi(i)$ 之和等于 $n$
（证明？观察归纳法（雾））

性质0.5：

欧拉函数是积性函数，即有 $\varphi(m\cdot n)=\varphi(m)\cdot \varphi(n)$ 。

性质1：（终于来了个够格的性质）

对于一个质数 $x$ ，$\varphi(x)=x-1$ （毕竟质数以下的所有数都不得不和他互质/kel）。

性质2：

对于一个次方数 $p^k$ （指一个质数 $p$ 的 $k$ 次方的结果，我自己编的QWQ） ，因为他充其量质因子只有 $p$ 一个，所以与他不互质的数只有 $\left\{p,2p,3p,\dots,p^{k-1}p\right\}$ 这里 $p^{k-1}$ 个 $p$ 的倍数，故有 $\varphi(p^k)=p^k-p^{k-1}$ ，提一个 $p^k$ 后稍加变形就有了下面的式子：

（可稍加理解，也可以当结论记）

性质3：

对于任意一个数 $a$ ，都能写作 $a=p_1^{k _ 1}\cdot p _ 2^{k _ 2} \cdot p _ 3^{k _ 3} \cdots $ 的次方数相乘的形式（类比分解质因式的短除法？），也就是 $a= {\textstyle \prod_{i=1}^{n}}p _ i^{k _ i}$ ，又因为性质0.5，欧拉函数是积性函数，所以 $\varphi(a)= {\textstyle \prod_{i=1}^{n}} \varphi(p _ i^{k _ i})$ ，用性质2又可以把 $\varphi(p _i^{k _ i})$ 变形成 $p_i^{k_i}\cdot\dfrac{p_i-1}{p _ i}$ ，那么我们就有：

还记得我们之前 $a$ 是怎么表示的嘛？是了， $a= {\textstyle \prod_{i=1}^{n}}p_i^{k_i}$ 。发现什么，我们可以把 $a$ 代入我们推出的式子！也就是说！！！

怎么求 $\varphi$ :

和 $\mathtt{exgcd}$ 一样，求 $\varphi$ 也有两种方法，下面分别给出：

单点求 $\varphi$:

适用情形：所求数较少，对时间要求高（因为这是一个 $O(\sqrt{n})$ 的做法）。

做法：依据我们已知的性质3可知，一个数 $a$ 不管他是什么，都满足 $\varphi(a)=a\cdot\textstyle\prod_{i=1}^{n}\dfrac{p_i-1}{p_i}$ ，其中 $p_i$ 是 $a$ 的质因子。那我们就根据这一性质，初始化一个 $ans=a$ ,像线性筛素数一样从1枚举到 $\sqrt{a}$ ，找到每一个质因子，累乘上 $\dfrac{p_i-1}{p_i}$ 。

为了防止这个 $a$ 本身是质数，导致枚举过程中一次也没更新 $ans$ ，结果返回了错误答案 $ans=a$ ，我们运用性质1，将其修改为 $ans * \dfrac{p_i-1}{p_i}$ ，最后得到 $ans=a-1$ 。

上代码～：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
int phi(int n)
{
    int ans = n, m = sqrt(n);
    for (int i = 2; i <= m; i++)
    {
        if (n % i == 0)
        {
            ans = ans / i * (i - 1);
            while (n % i == 0)
                n /= i;
 //枚举质因子，之后除去该因子，避免这个因子被多次统计
        }
    }
    if (n >= 2)
        ans = ans / n * (n - 1);
      //特殊处理
    return ans;
}
线性筛 $\varphi$ ：

适用情形：所求数较多，多次 $O(\sqrt{n})$ 复杂度不如单次 $O(n)$ 时。

做法：从2到 $n$ 跑循环，找到未 $vis$ 的一个数 $i$ ，则他是质数，将其加入质数表 $prime$ ，将 $\varphi(i)$ 标为 $i-1$ （性质1）。

再枚举倍数 $j$ ，将这个质数 $i$ 的 $j$ 倍标记为合数，也就是把 $vis$ 标记成1，枚举当前质数表内质数，如果手上这个合数不是一个质数的倍数，那这个合数和这个质数就只能互质。

我们干脆多做点事，求一下他们积的 $\varphi$ 。

已知两个互质数求他们的积的 $\varphi$ ？？？想到了什么，性质0.5！！！ 记得吗？ $\varphi(m\cdot n)=\varphi(m)\cdot \varphi(n)$ ，那么我们就有 $\varphi(i\cdot prime_j)=\varphi(i)\cdot \varphi(prime_j)$ ！！！

那不互质怎么办？我们提到，不互质，就只能是倍数关系。要是倍数关系的话，手上的这个合数根本没有新质因子增添进来。那么我们可以运用性质3，
$a\cdot\textstyle\prod_{i=1}^{n}\dfrac{p_i-1}{p_i}$ 里的 $i$ 取哪些值不会变（因为质因子不变多），只是前面这个 $a$ 乘了 $prime_j$，这个式子在原先 $\varphi(i)$ 的基础上乘上了 $prime_j$ ，这就产生了新的 $\varphi(i\cdot prime_j)=\varphi(i)\cdot prime_j$ ！！！

注意特判 $\varphi(1)=1$ ！！！！！

话不多说，上代码～：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
int phi()
{
    phi[1] = 1;
    for (int i = 2; i <= 1000007; i++)
    {
        if (!vis[i])
        {
            prime[++cntp] = i;
            phi[i] = i - 1;
        }
        for (int j = 1; prime[j] && i * prime[j] <= 1000007; j++)
        {
            vis[i * prime[j]] = 1;
            if (!(i % prime[j])) //case 1：倍数关系
            {
                phi[i * prime[j]] = phi[i] * prime[j];
                break;
            }
            else //case 2:互质关系
                phi[i * prime[j]] = phi[i] * phi[prime[j]];
        }
    }
}

扩 展 欧 拉 定 理

在学扩展之前，这里先提一下欧拉定理来引入主题（工具人）。欧拉定理实则无比简洁，内容如下：

（这里要求 $a,k$ 互质，证明极其麻烦，记结论即可）

还有个叫费马小定理的，写的是 $a^{k-1}≡1 \pmod k $ ，这里要求的 $k$ 是质数。看出点什么了吗，是呐，费马小定理就是欧拉定理的一种特殊情况，在 $k$ 是质数这种情况下， $\varphi(k)$ 也就成了 $k-1$ 。

然而这并不是今天的重点，我们还是得请出今天的主角——扩展欧拉定理！

我们先来看看这个定理长什么样~

其中第一行的条件是 $a,b$ 互质，第二、三行则无须互质，这里的每个算式都是在 $\mod m$ 意义下的。

我们这里尝试推导第三个式子（不感兴趣的读者可自行跳过）：

先来考虑一个质数 $p$ 。

我们设有一个 $m=s\cdot p^r$ ，并满足 $\gcd(s,p)=1$ 。（设，都可以设，只要合理）

就有 $p^{\varphi(s)}≡1\pmod s$ （这是上文提到的欧拉定理），又依据 $m=s\cdot p^r$ 且 $\gcd(s,p)=1$ ，知道了 $m$ 必是 $s$ 的倍数，根据欧拉函数的积性（性质0.5），可知 $\varphi(m)$ 必是 $\varphi(s)$ 的倍数，那让同余方程的左边再乘上若干个 $p^{\varphi(s)}$ ，把 $p$ 的指数补成 $\varphi(m)$ ，右边再乘上若干个1，又有什么影响（有了同余方程的同乘性，说话就是硬气）。那我们手上就变成了：

再来，有了同乘性还有什么做不了的，那两边同乘上 $p^r$ 也是可以的，那就有：

那么我们再来看这个式子：

有性急的同学要开骂了，你这不是显而易见的式子嘛。不急，看下面的变化，我们有了 $p^{\varphi(m)+r}≡p^r\pmod s$ ，不妨把 $p^r$ 代入这个式子，就有:

借助扩展欧拉定理的式子一我们又有了

上面这个式子就可以作为我们的一个半成品了。

相同的规律也可以用在任何一个次方数上，悟一悟下面这个长长的同余式：

其实这一波推导实则通篇只是欧拉定理的推论，结合欧拉定理食用更易消化~

还记得吗，我们之前曾说过任何一个数 $a$ ，他都能表示成 $p_1^{k_1}\cdot p_2^{k_2} \cdot p_3^{k_3} \cdots $ 的次方数相乘的形式。对于每一个 $p^k$ 都满足上面我们推出的公式，那么就有！！！：

$\Box$

高 斯 消 元(Gauss)

这有什么用：

高斯消元，是指一种通过把 $N$ 元一次方程组（下面都以三元一次方程组为例）里的系数，写在一个矩阵里，像下面这样：

其目的是将这个矩阵变换成下面的模样以求解：

这样我们就可以从矩阵最下面一行逐步推出每一个未知数的解啦！（一行一杀，逐行上推）

我们的思想是先找到哪一行的第 $i$ 个数绝对值最大，将这一行全部换到第 $i$ 行的位置 （这里的 $i$ 从1到 $N$ 枚举），如下图的变换(这里我们找到了第一项系数最大值在第三行我们将第三行换到了第一行的位置），我们称之为操作一。

操作一完成以后，再对于第 $i$ 项进行加减消元，也就是要把他下面的每一行的第 $i$ 项系数全部削成0，称作操作二，如下图：

两种操作交替运行，这样最后我们就可以从后往前倒推，求出 $z=1,y=3,x=2$ 辣！

代码来辣：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
for (int i = 1; i <= n; i++)
{    //枚举第i个未知数
    int maxx = i;                    
    //绝对值最大的项所在行
    for (int j = i + 1; j <= n; j++) 
    //从第i行开始往下搜，因为上面的都已经被处理了
    {
        if (fabs(a[j][i]) > fabs(a[maxx][i]))
        {
            maxx = j; //打擂台
        }
    }
    for (int j = 1; j <= n + 1; j++)
        swap(a[i][j], a[maxx][j]); //全体交换
    if (!a[i][i])  //换完以后发现这行要求的未知数前面系数是零
    {
        printf("No Solution\n"); //无解，没得搞了
        return 0;
    }
    for (int j = 1; j <= n; j++)//开消
    {
        if (j == i)
            continue;                   
       //不能拿这一列式子消他自己
        double delta = a[j][i] / a[i][i]; 
        //消元时被削那行需要乘/除上的数
        for (int k = i; k <= n + 1; k++)
            a[j][k] -= a[i][k] * delta; 
       //将第i元系数消为0并对应的把系数消掉
    }
}

这样我们就成功得到了最终我们需要的那个矩阵。在主函数里稍加处理即可。笔者这里采取的是在输出时处理出未知数的值。

1
printf("%.2lf\n", num[i][n + 1] / num[i][i]);

高斯消元一种妙用——解异或方程组：

先咕咕咕~

中 国 剩 余 定 理 （$\mathtt{CRT}$）

他解决什么问题：

举个例子，有一方程组如右： $\begin{cases}x≡2\pmod 3\\x≡3\pmod 5\\x≡2\pmod 7\end{cases}$ 要求求出最小的 $x$ （这里的模数都应该互质）。对于这种问题，我们写成一个通式如下 $\begin{cases}x≡a_1\pmod {m_1}\\x≡a_2\pmod {m_2}\\x≡a_3\pmod {m_3}\end{cases}$

解决的步骤：

  1. 计算出一个 $M=m_1\cdot m_2\cdots m_n$ 。

  2. 开一个 $RBQ$ 数组（老数组名了），对于每个 $RBQ_i$ ，有 $RBQ_i=\dfrac{M}{m_i}$ （即除了 $m_i$ 以外的所有 $m$ 的积），再求出每一个 $RBQ$ 的在 $\bmod m_i$ 意义下的逆元 $RBQ_i^{-1}$ （有时并不用开数组， $RBQ$ 也可即造即用(雾））。

  3. 那么我们所求的 $x$ 就有（不用害怕原数乘逆元等于一，因为一个是在实数域上，一个是在模 $m_i$ 意义下）：

为什么这么做是对的：

对于每一行式子的一个 $m_i$ ，我们有两种情况：

$Case 1(j\neq i)$

因为 $RBQ_j$ 一定是 $m_i$ 的倍数，它本来就是包括 $m_i$ 在内的几个 $m$ 乘起来的结果。

$Case 2(j= i)$

这不必说。

我们都知道的 $RBQ_i \cdot RBQ_i^{-1}\pmod {m_i}$ 他就是1，那这个式子自然成立。

了解了上面两种情况，那么证明就极其简单了：

也就是说让 $x=\sum\limits_{k=1}^n a_k\cdot RBQ_k \cdot RBQ_k^{-1}$ 总能完美契合方程组里每一个同余式，对于里面的每一个同余式都能得到相应的余数。 $\Box$

代码整点康康:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
void exgcd(int a, int b)
{
    if(b==0)
    {
        x = 1;
        y = 0;
        return;
    }
    exgcd(b, a % b);
    int tmp = x;
    x = y;
    y = tmp - (a / b) * y;
}//扩欧求乘法逆元，无需多言 
signed main()
{
    int N;
    cin >> N;
    for (int i = 1; i <= N;i++)
    {
        cin >> m[i] >> a[i];
        M *= m[i];//累乘模数
    }
    for (int i = 1;i<=N;i++)
    {
        RBQ= M / m[i];
        exgcd(RBQ,m[i]);//即造即用RBQ 
        sum_a = (sum_a+a[i] * RBQ * x)%M;//累加a[i]*RBQ*RBQ的逆元 
    }
    cout << sum_a  % M << endl;//具体输出可按题意修改
}

扩 展 中 国 剩 余 定 理 ($\mathtt{EXCRT}$)

我们之前提到 $\mathtt{CRT}$ 解决的问题中方程组的各个模数必须是互质的，可有的问题给出的方程组不一定是互质的，难道这种题就不可做了吗？不慌，我们还有 $\mathtt{CRT}$ 的进阶与扩展， $\mathtt{EXCRT}$ ！

他如何运作：

我们举个例子，有一方程组如下：

我们不妨考虑每次合并第一个方程和第二个方程，而不改变新的方程组求出的解的正确性，最后直至消到只剩一个方程，就可以轻松出结果。那么怎么合并呢，我们挑出两方程来看：

把同余式化成等式的做法我们已经轻车熟路了，我们就有了：

来看右边，我们稍微来一点移项：

我们会觉得这挺像一个同余方程的，想想办法解一下他，要让他可解我们必须让右边变成 $\gcd (m_1,m_2)$ ，没什么可以阻止我们在等号左右同乘上 $\dfrac {\gcd (m_1,m_2)}{a_2-a_1}$ ，由此一来，我们有：

那个大分数也太 $ex$ 了，我们不得不把他换元，分别换成 $K_1,K_2$ 如何：

漂亮的同余方程， $\mathtt{exgcd}$ 一波，推出符合的一组 $K_1,K_2$ ，但这不是我们需要的小 $k$ ，不妨一波返场，对大 $K$ 再乘回 $\dfrac {a_2-a_1}{\gcd (m_1,m_2)}$ ，那没事了，我们这就知道了 $k_1$ 的一个可以的值，把他代回到 $x=m_1\cdot k_1+a_1$ ，我们就暂时知道了一个 $x$ 。

别急，我们还没有做完，这个 $x$ 只能算是局部解，我们应该记得我们此行的目的是把前面两个方程合并到一个方程，怎么搞，我们需要知道新 $a_1$ 和新 $m_1$ ，新 $m_1$ 怎么想都应该知道该是旧 $m_1$ 和旧 $m_2$ 的 $\operatorname{lcm}$，我们用 $\dfrac {m_1\cdot m_2}{\gcd(m_1,m_2)}$ 可以算出最小公倍数（我想这个应该都知道的），那么对于新 $a_1$ ，求之前我们需要知道一个有趣的定理：

（这其实不能叫定理了，笔者对其根据本题情况稍加修改， $x_{common}$ 指的是 $x$ 的通解，我们之前用的是 $K$ ， $x_{special}$ 指的是 $x$ 的特解， $i$ 是随意的一个正整数， $m_1,m_2$ 正如先前所设）

应该会有人还记得，我们要求的 $K$ 的解，也就是这里的 $x$ ，他本该是个什么，是了， $\dfrac {\gcd (m_1,m_2)}{a_2-a_1}k_1$ （是我们嫌他太 $ex$ 才把他换元的） ，那么是时候把他换回来， $x_{special}$ 肯定是不能动了，毕竟我们好不容易才把他求出来，那么 $x_{common}$ 呢，我们尝试变形:

这就是 $k_1$ 的通解了，那么再带回 $x=m_1\cdot k_1+a_1$ ？那么我们又有了:

把 $m_1$ 乘进括号里（分配律），再一波变形如下：

别的显然都不能动了，我们手上的棋子只有 $i$ 这一颗 ，只能说 $x$ 的大小取决于 $\operatorname{lcm}(m_1,m_2)$ 前面的系数，等等，$\operatorname{lcm}(m_1,m_2)$ 是什么，这不是我们定下来的新 $m_1$ 嘛。看来我们已经离胜利只差一步了。什么时候让 $x$ 有最小解，简单， $i=0$ ，那我们就有：

噫，好了，这个形式怎么看怎么爽，是呢，我们之前不管不顾地让 $i=0$
，可 $i$ 的值是不定的，他的背后就是那个模数 $\operatorname{lcm}(m_1,m_2)$ ，那么最终的式子已经呼之欲出了！

这，就是合并了第一第二式的最终式子（呼～）。

不妨对着代码格物致知？：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
int EXCRT()
{
	__int128 A1,A2,MOD1,MOD2;
    //用int128是因为被毒瘤数据爆了long long，千万小心
	A1=A[1],MOD1=MOD[1];
    //a1,m1整来
	for(int i=2;i<=N;i++)
	
		A2=A[i],MOD2=MOD[i];
        //a2,m2整来
		__int128 C=A2-A1,GCD=__gcd(MOD1,MOD2);
        //换个a2-a1,gcd免得式子太ex
		exgcd(MOD1,MOD2);
        //解同余方程
		x=((x*C/GCD)%(MOD2/GCD)+(MOD2/GCD))%(MOD2/GCD);
        //处理一下x，帮助a1出新值
		A1=x*MOD1+A1;
		MOD1=MOD1*MOD2/GCD;
        //a1,m1各更新

	}
	return A1;
}

我NOIp之前再碰数论我就吃奥利给！

完结撒花～

最后更新： 2021年04月19日 11:33

原始链接： https://quest233.github.io/2021/04/17/%E6%B5%85%E8%B0%88%E7%AE%80%E5%8D%95%E6%95%B0%E5%AD%A6/