数学进阶中的数论与组合魔法：从同余到高斯消元

大家好！在编程学习的路上，我们常常会遇到一些有趣的数学问题。比如，你知道今天星期三，那么100天后是星期几吗？其实只要看100除以7的余数（100 mod 7 = 2），星期三加两天就是星期五。这种“只看余数”的玩法就是同余的核心思想。掌握同余以及相关的数论知识，不但能轻松解决生活中的周期问题，还能在密码学、组合计数、算法优化中发挥巨大作用。这篇文章将带你从同余出发，一路走过模逆元、扩展欧几里得、中国剩余定理、欧拉定理、容斥原理、卡特兰数，最后到达矩阵快速幂和高斯消元。每个知识点都会用生活例子帮你理解，用C++代码教会你实现。

1. 同余式与模逆元：星期几的数学秘密

什么是同余？

同余就像“除以某个数后只看余数”。如果两个数a和b除以m的余数相同，我们就说 a ≡ b (mod m)。比如：

今天星期三，再过10天，10 mod 7 = 3，所以相当于星期三+3天=星期六。用同余式写就是 10 ≡ 3 (mod 7)。
你的零花钱有15元，朋友有8元，你想知道两人钱数模5是不是一样？15 mod 5 = 0，8 mod 5 = 3，不一样。

生活中： 同余可以用来算周期。比如学校每隔7天放一次假，今天放假，那么第19天后放假吗？19 mod 7 = 5，所以不会，要再过5天才放假。

什么是模逆元？

在普通算术里，一个数a的倒数（乘法逆元）是满足 a × x = 1 的数x。在模的世界里，我们想找这样一个数x，使得 a × x ≡ 1 (mod m)。这个x就叫 a模m的逆元，记作 a⁻¹ mod m。

例子： 在模7的世界里，2的逆元是4，因为2×4=8≡1 (mod 7)。就像数字的“倒数”一样，有了逆元，我们就可以做“模除法”：计算 (5 ÷ 2) mod 7，就是5乘以2的逆元4，得到20 mod 7 = 6。

怎么找逆元？

最简单的方法是暴力尝试，从1到m-1一个一个试，看哪个数乘a后模m等于1。但这种方法很慢（尤其m很大时）。C++里可以用扩展欧几里得算法快速求出（后面会学到）。

下面是用暴力法求逆元的代码（适合m较小且保证逆元存在的情况）：

#include <iostream>
using namespace std;

// 求 a 在模 m 下的逆元（暴力法，m较小且保证存在）
int modInverse(int a, int m) {
    for (int x = 1; x < m; x++) {          // x从1试到m-1
        if ((a * x) % m == 1)              // 如果乘起来模m等于1
            return x;                      // x就是逆元
    }
    return -1; // 不存在逆元
}

int main() {
    int a = 2, m = 7;
    cout << "2 的模 7 逆元是: " << modInverse(a, m) << endl; // 输出4
    return 0;
}

常见错误与注意

逆元存在的条件： a和m必须互质（最大公约数为1）。如果不互质，逆元不存在。例如2在模4下就没有逆元，因为2和4不互质（gcd=2），找遍0~3都没有数使得2×x≡1 mod 4。
暴力法只适合小模数： 当m很大（比如10⁹）时，循环m次会超时，必须用扩展欧几里得。
结果可能是负数？ 逆元通常取0到m-1之间的正整数，如果扩展欧几里得返回了负数，要加m调整为正。

2. 扩展欧几里得算法：水壶倒水的智慧

它用来干什么？

欧几里得算法（辗转相除法）可以快速求两个数的最大公约数gcd(a,b)。但有时候我们不仅想知道gcd，还想解方程 ax + by = gcd(a,b) 的整数解(x,y)。这就是扩展欧几里得算法的本领。

生活例子： 你有两个水壶，容量分别为a升和b升，你想倒出恰好等于它们最大公约数升的水。每次可以装满、倒空或互相倒水。扩展欧几里得就能告诉你最少需要多少次操作，以及每次怎么倒（实际上它给出了x和y，代表a的倍数加上b的倍数等于gcd）。

算法核心

递归求解：已知 gcd(b, a % b) 可以求出它们的解，然后反推回原方程。公式推导如下：

递归基：当 b=0 时，gcd(a,0)=a，此时方程变为 a·x + 0·y = a，显然 x=1, y=0。
递归过程：设 g = exgcd(b, a % b, x, y) 返回了gcd，并得到了方程 b·x + (a%b)·y = g 的一组解。注意 a%b = a - (a/b)*b，代入整理可得 a·y + b·(x - (a/b)*y) = g。所以新的x' = y, y' = x - (a/b)*y。

代码实现如下：

#include <iostream>
using namespace std;

// 扩展欧几里得，返回gcd，并设置x,y使 ax + by = gcd
int exgcd(int a, int b, int &x, int &y) {
    if (b == 0) {            // 递归结束
        x = 1;
        y = 0;
        return a;
    }
    int g = exgcd(b, a % b, x, y);  // 递归求下一层
    int temp = x;                   // 保存原来的x
    x = y;                          // 新x = 旧y
    y = temp - (a / b) * y;         // 新y = 旧x - (a/b)*旧y
    return g;
}

int main() {
    int x, y;
    int g = exgcd(18, 12, x, y);
    cout << "gcd = " << g << ", x = " << x << ", y = " << y << endl; // 得 gcd=6, x=1, y=-1
    // 验证：18*1 + 12*(-1) = 6
    return 0;
}

求模逆元

还记得之前要解 a * x ≡ 1 (mod m) 吗？这等价于 a·x + m·y = 1（因为a·x - 1是m的倍数，设为 m·(-y)）。用exgcd(a, m, x, y)求出x后，调整x为正数就是a的逆元。例如求2模7的逆元：exgcd(2,7,x,y)得到x=-3，y=1，因为-3 mod 7 = 4，所以逆元是4。

// 用 exgcd 求逆元（更快）
int modInverseExgcd(int a, int m) {
    int x, y;
    int g = exgcd(a, m, x, y);
    if (g != 1) return -1;      // 不互质，无逆元
    return (x % m + m) % m;     // 调整为正数
}

常见错误

忘记调整x为0~m-1：exgcd可能返回负数，须加上模数再取模。
调用时参数顺序：exgcd(a, m, x, y) 求解 a·x + m·y = gcd(a,m)，注意是a和m，不要搞反。
当gcd不为1时逆元不存在，要检查返回值。

3. 中国剩余定理（CRT）：多个闹钟同时响的时刻

故事背景

《孙子算经》有个经典问题：“一个数除以3余2，除以5余3，除以7余2，问这个数最小是多少？”答案23。中国剩余定理（CRT）就是解决这类“多个模条件”的通用方法。

生活例子： 你有很多个闹钟，闹钟1每3小时响一次，第一次在2小时后；闹钟2每5小时一次，第一次在3小时后；闹钟3每7小时一次，第一次在2小时后。下次所有闹钟同时响起是什么时候？CRT告诉你最早时间是23小时后。

算法步骤

计算总模数M：把所有模数乘起来，M = m1 × m2 × ... × mk。
对每个条件i：
- 计算 Mi = M / mi（去掉当前模数）。
- 求 Mi 在模 mi 下的逆元 ti（即解 Mi·ti ≡ 1 (mod mi)）。
组合答案：答案 ans = sum( ai × Mi × ti ) mod M，并取正数。

代码实现（用之前写的exgcd求逆元）

#include <iostream>
using namespace std;

// 扩展欧几里得求逆元（省略，直接用之前函数）
int exgcd(int a, int b, int &x, int &y) {
    if (b == 0) { x = 1; y = 0; return a; }
    int g = exgcd(b, a % b, x, y);
    int temp = x;
    x = y;
    y = temp - (a / b) * y;
    return g;
}

// 求逆元（模m，a与m互质）
int inv(int a, int m) {
    int x, y;
    exgcd(a, m, x, y);
    return (x % m + m) % m;
}

// 中国剩余定理，返回最小非负整数解
int crt(int n, int a[], int m[]) {
    int M = 1, ans = 0;
    for (int i = 0; i < n; i++) M *= m[i];            // 总模数
    for (int i = 0; i < n; i++) {
        int Mi = M / m[i];                            // 去掉当前模数
        int ti = inv(Mi % m[i], m[i]);                // Mi模mi的逆元
        ans = (ans + (long long)a[i] * Mi % M * ti) % M;
    }
    return ans;
}

int main() {
    int a[] = {2, 3, 2}; // 余数
    int m[] = {3, 5, 7}; // 模数
    cout << "答案: " << crt(3, a, m) << endl; // 输出23
    return 0;
}

常见错误与注意

模数必须两两互质：如果模数不互质，上面的简单CRT不成立（结果可能错误）。需要更复杂的合并方法（例如用扩展欧几里得合并两个同余式）。
答案需要取模：最后ans是M的剩余类，取模后就是最小非负解。
注意整数溢出：a[i]*Mi可能很大，先用long long计算再取模。

4. 欧拉函数与欧拉定理：数出“朋友”个数，简化大幂运算

欧拉函数 φ(n) 是什么？

欧拉函数 φ(n) 表示1到n中与n互质的整数个数。例如 φ(6) = 2，因为1和5与6互质（2、3、4与6不互质）。怎么算呢？先对n分解质因数：n = p1^e1 × p2^e2 × ...，则 φ(n) = n × (1-1/p1) × (1-1/p2) × ...。

生活例子： 班级里有n个同学，其中有多少人和你（编号为某个数）没有共同的班级编号因子？比如n=6，你想找到和你“互质”的同学，只有2个人。

欧拉定理

如果a和n互质，那么 a^φ(n) ≡ 1 (mod n)。这个定理超级强大！比如计算 7^100 mod 10，因为φ(10)=4，7和10互质，所以7^4 ≡ 1 (mod 10)，那么7^100 = (7^4)^25 ≡ 1^25 = 1。不用真的算7的100次方。

代码实现欧拉函数（试除法）

#include <iostream>
using namespace std;

int phi(int n) {
    int res = n;                          // 初始化为n
    for (int i = 2; i * i <= n; i++) {    // 试除到√n
        if (n % i == 0) {                 // i是质因子
            res = res / i * (i - 1);      // 应用公式 n*(1-1/i)
            while (n % i == 0) n /= i;    // 除掉所有i因子
        }
    }
    if (n > 1) res = res / n * (n - 1);   // 如果还有大质因子
    return res;
}

int main() {
    cout << "phi(6) = " << phi(6) << endl;  // 2
    cout << "phi(10) = " << phi(10) << endl; // 4
    return 0;
}

用欧拉定理做快速幂

假设我们要计算 a^b mod n，其中a与n互质，可以用 b = q·φ(n) + r 来简化：a^b ≡ a^r (mod n)。不过更通用的方法还是快速幂。

常见错误

忘记检查互质性：欧拉定理要求a和n互质，如果gcd(a,n)≠1，不能直接用。
φ(n)计算要完整：分解质因数时，要确保把所有因子都除掉，最后剩余的>1也要处理。

5. 费马小定理与威尔逊定理：质数的魔法性质

费马小定理

如果p是质数，且a不是p的倍数（即gcd(a,p)=1），那么 a^(p-1) ≡ 1 (mod p)。这是欧拉定理的特例，因为对于质数p，φ(p)=p-1。

应用1：快速计算大幂模p
例如计算 2^100 mod 7。因为7是质数，2^6 ≡ 1 (mod 7)，所以2^100 = 2^(6×16+4) = (2^6)^16 × 2^4 ≡ 1^16 × 16 mod 7 ≡ 2 (mod 7)。

应用2：质数判定（不靠谱）
如果对于某个a，a^(p-1) ≡ 1 mod p不成立，则p一定不是质数。但反过来不成立（有一些合数欺骗性很强，称为伪素数），比如2^340 ≡ 1 mod 341（341=11×31），所以费马小定理不能完全用来判定质数，但可以快速筛掉大多数合数。

威尔逊定理

p是质数当且仅当 (p-1)! ≡ -1 (mod p)。例如p=5，(4!)=24 ≡ 4 ≡ -1 (mod 5)。这个定理理论意义大，但直接算阶乘太大，不能用于实际判质。

代码：快速幂验证费马小定理

#include <iostream>
using namespace std;

// 快速幂计算 a^b mod p
long long quickPow(long long a, long long b, long long p) {
    long long res = 1;
    while (b) {
        if (b & 1) res = res * a % p;   // 如果最低位是1
        a = a * a % p;                  // a自乘
        b >>= 1;                        // 右移
    }
    return res;
}

int main() {
    long long p = 7;
    cout << "2^100 mod 7 = " << quickPow(2, 100, p) << endl; // 2
    // 验证费马小定理：2^6 mod 7 = 1
    cout << "2^6 mod 7 = " << quickPow(2, 6, p) << endl; // 1
    return 0;
}

常见错误

误用费马小定理：p必须是质数且a不被p整除。如果p不是质数，结果不一定成立。
快速幂乘法溢出：两个long long相乘可能超过范围，用 (a*a)%p 已经用模运算保证了，但中间乘积在C++中可能先溢出，建议用 128位或写乘法取模函数（如快速乘）。

6. 容斥原理：数人数时去掉重复的

基本思想

容斥原理用来计算多个集合并集的元素个数。最简单的公式：|A∪B| = |A| + |B| - |A∩B|。比如班里学数学的有10人，学语文的有8人，两科都学有3人，那么总人数是10+8-3=15。

对于三个集合，公式是： |A∪B∪C| = |A|+|B|+|C| - |A∩B| - |A∩C| - |B∩C| + |A∩B∩C|。

一般地，奇数个交集相加，偶数个交集相减。

生活例子

求1到100中能被2、3或5整除的数的个数。直接数很麻烦，用容斥：

能被2整除：100/2=50个
能被3整除：33个
能被5整除：20个
能被2和3（即6）整除：16个
能被2和5（即10）整除：10个
能被3和5（即15）整除：6个
能被2、3、5（即30）整除：3个总数 = 50+33+20 -16-10-6 +3 = 74。

代码实现：枚举子集（二进制）

#include <iostream>
using namespace std;

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

int main() {
    int n = 100;
    int primes[] = {2, 3, 5};
    int k = 3;
    int cnt = 0;
    // 枚举所有非空子集（用二进制位表示）
    for (int mask = 1; mask < (1<<k); mask++) {
        int l = 1, bits = 0;
        for (int i = 0; i < k; i++) {
            if (mask & (1<<i)) {
                l = lcm(l, primes[i]);   // 计算子集元素的最小公倍数
                bits++;
            }
        }
        int num = n / l;                  // 能被这个最小公倍数整除的个数
        if (bits % 2 == 1) cnt += num;    // 奇数个元素加
        else cnt -= num;                  // 偶数个元素减
    }
    cout << "1~100中能被2,3,5之一整除的数有: " << cnt << "个" << endl; // 74
    return 0;
}

常见错误

最小公倍数可能溢出：当数字较大时，l = lcm(l, prime) 可能超过int范围，应使用long long或适时判断。
重复减加顺序：注意符号，奇数加偶数减。
忘记包含空子集：mask从1开始，不包括空集。

7. 卡特兰数：排队不插队的计数规则

它是什么？

卡特兰数是组合数学里一系列整数，第n项记为C_n。它有着各种各样的几何和组合意义。最著名的定义是：n对括号的合法匹配数。例如n=3时，有5种：
((())) (()()) (())() ()(()) ()()()

其他等价问题：栈的出栈序列数、n个节点能形成的二叉树个数、凸多边形三角形划分方法数等。

生活例子： 你有一串票，要按顺序进游乐场，但排队时想用栈来中转，有多少种不同的排队顺序？也是卡特兰数。

计算公式

递推式：C0=1, Cn = (2*(2n-1)/(n+1)) * C(n-1)
封闭公式：Cn = (1/(n+1)) * C(2n, n)

代码实现（递推求前10项）

#include <iostream>
using namespace std;

int main() {
    long long catalan[20] = {1};               // C0=1
    for (int i = 1; i <= 10; i++) {
        catalan[i] = catalan[i-1] * 2 * (2*i - 1) / (i + 1);
        cout << "C(" << i << ") = " << catalan[i] << endl;
    }
    return 0;
}

常见错误

整数除法的顺序：递推式中除法必须保证整数。先乘后除，由于分子包含 (i+1) 的因子，结果一定是整数。
数字增长很快：n=20时卡特兰数约6.56e9，已经超过int，要用long long或高精度。

8. 矩阵快速幂：飞快地算递推数列

背景

普通快速幂可以快速计算 a^b。如果把a换成矩阵，就能在O(logn)时间内计算矩阵的n次幂。这在解递推数列时特别有用。例如斐波那契数列F0=0, F1=1, Fn=Fn-1+Fn-2。可以写成矩阵乘法：

[ Fn ] = [1 1] ^ (n-1) * [F1] [ Fn-1 ] [1 0] [F0]

只要计算出矩阵的 (n-1) 次幂，再乘以初始向量，就能得到Fn。矩阵乘法有结合律，所以可以用快速幂。

代码实现（求斐波那契第n项模1e9+7）

#include <iostream>
#include <vector>
using namespace std;

typedef vector<vector<long long>> Matrix;   // 矩阵类型
const long long MOD = 1e9 + 7;

// 矩阵乘法
Matrix mul(Matrix a, Matrix b) {
    int n = a.size();
    Matrix c(n, vector<long long>(n, 0));
    for (int i = 0; i < n; i++)
        for (int k = 0; k < n; k++)           // 优化：交换k和j循环可提高缓存友好性
            for (int j = 0; j < n; j++)
                c[i][j] = (c[i][j] + a[i][k] * b[k][j]) % MOD;
    return c;
}

// 矩阵快速幂
Matrix matrixPow(Matrix a, long long p) {
    int n = a.size();
    Matrix res(n, vector<long long>(n, 0));
    for (int i = 0; i < n; i++) res[i][i] = 1; // 单位矩阵
    while (p) {
        if (p & 1) res = mul(res, a);
        a = mul(a, a);
        p >>= 1;
    }
    return res;
}

// 求斐波那契第n项（F0=0, F1=1）
long long fib(long long n) {
    if (n == 0) return 0;
    Matrix base = {{1, 1}, {1, 0}};
    Matrix m = matrixPow(base, n - 1);
    return m[0][0]; // F_n
}

int main() {
    cout << "F(10) = " << fib(10) << endl; // 55
    return 0;
}

常见错误

矩阵乘法维度：必须是方阵，且大小相同。
快速幂的循环：指数p要反复除以2，直到0。
初始化单位矩阵：对角线为1，其他地方0。
乘法取模：每个加法后取模，避免溢出。

9. 高斯消元法：解方程组的系统方法

它是什么？

解方程组时，我们常常通过“消去未知数”来逐步化简。高斯消元法用矩阵（增广矩阵）的形式，通过行变换（交换行、倍乘行、一行加若干倍另一行）将矩阵化为行阶梯形，然后回代求解。

生活例子： 你有两种铅笔，一次买1支A和2支B花了5元，一次买3支A和4支B花了11元。问A和B各多少钱？解方程组： 1·x + 2·y = 5 3·x + 4·y = 11
高斯消元能帮你找到x=1, y=2。

算法步骤

建立增广矩阵：把系数和常数写在一起。
对于第col列（从0开始），在当前行row开始找该列绝对值最大的元素（选主元），交换到当前行。
用主元行消去下面所有行的第col列（把下面的该列变为0）。
重复下一列，直到处理完所有列或行。
回代：从最后一行开始，逐行求解每个未知数。

代码实现（浮点数）

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;

const double EPS = 1e-9;   // 误差阈值

vector<double> gauss(vector<vector<double>> a) {
    int n = a.size();      // 方程个数
    for (int col = 0, row = 0; col < n && row < n; col++) {
        int sel = row;
        for (int i = row; i < n; i++)               // 选主元
            if (fabs(a[i][col]) > fabs(a[sel][col]))
                sel = i;
        if (fabs(a[sel][col]) < EPS) continue;      // 该列全为0，跳过
        swap(a[row], a[sel]);                       // 交换到当前行
        for (int i = row + 1; i < n; i++) {         // 消去下面行的该列
            double f = a[i][col] / a[row][col];
            for (int j = col; j <= n; j++)
                a[i][j] -= f * a[row][j];
        }
        row++;
    }
    // 回代
    vector<double> x(n, 0);
    for (int i = n - 1; i >= 0; i--) {
        double sum = a[i][n];
        for (int j = i + 1; j < n; j++)
            sum -= a[i][j] * x[j];
        if (fabs(a[i][i]) < EPS) {                  // 无解或无穷多解
            // 这里简单返回空向量；实际需要额外处理
            return vector<double>();
        }
        x[i] = sum / a[i][i];
    }
    return x;
}

int main() {
    vector<vector<double>> a = {
        {1, 2, 5},
        {3, 4, 11}
    };
    vector<double> x = gauss(a);
    if (x.empty()) {
        cout << "无解或无穷多解" << endl;
    } else {
        cout << "x = " << x[0] << ", y = " << x[1] << endl; // x=1, y=2
    }
    return 0;
}

常见错误

浮点精度问题：判断是否为0时用EPS，不要直接比较 ==0。
主元选择：如果选到绝对值很小的数，会导致数值不稳定，所以通常取绝对值最大的。
无解或无穷多解：当回代时某行系数全0但常数不为0，则无解；全0常数也为0则有无穷多解。实际程序中要判断。

完整综合示例：用数论知识解决“星期几”扩展问题

假设今天星期一，问经过 a 天后是星期几？如果 a 很大（比如10^100），直接计算很麻烦。我们可以利用费马小定理或欧拉定理简化。例如模7的世界里，10^100 mod 7 = ？因为7是质数，10 mod 7 = 3，3^6 ≡ 1 mod 7，所以3^100 = 3^(6*16+4) ≡ 3^4 = 81 ≡ 4 mod 7，所以从星期一加4天是星期五。这比直接算10^100快多了！

#include <iostream>
using namespace std;

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

int main() {
    // 题目：今天星期一，问 10^100 天后是星期几？
    // 星期天为0, 星期一为1, ..., 星期六为6
    // 返回 (1 + 10^100) mod 7
    int today = 1; // 星期一
    long long base = 10, exp = 100, mod = 7;
    long long days = quickPow(base % mod, exp, mod);
    int result = (today + days) % mod;
    cout << "10^100 天后是星期" << result << "（0~6分别表示星期日~星期六）" << endl;
    return 0;
}

例题精讲