C++实现大衍求一术：从中国剩余定理到现代代码实践-育师

1. 项目概述：从“物不知数”到现代代码

“大衍求一术”这个名字，听起来就带着一股历史的厚重感。它不是什么新的编程框架，也不是某个时髦的算法，而是源自南宋数学家秦九韶在《数书九章》中记载的一套古老算法，核心是求解一次同余式组，也就是我们常说的“中国剩余定理”问题。最经典的例子就是《孙子算经》里的“物不知数”：有一堆物品，三个三个数剩两个，五个五个数剩三个，七个七个数剩两个，问总数是多少？答案是二十三。这个算法在数论和密码学里有着深远的影响，比如RSA加密、编码理论等领域都能看到它的身影。

今天，我们不谈晦涩的数学史，也不做纯理论推导。我的目标很直接：用现代C++，把这份近八百年前的数学智慧，变成一个清晰、高效、可复用的代码模块。为什么是C++？因为我们需要对整数运算有绝对的控制力，需要追求极致的性能（尤其是在处理大整数时），同时，C++的面向对象特性又能让我们优雅地封装这个算法的各个步骤。无论你是正在学习数论和算法的学生，还是需要在项目中处理模运算、同余方程的开发者，亦或是对古典算法现代化实现感兴趣的爱好者，这篇文章都将带你走通从原理理解到代码落地的完整路径。你会发现，古人的智慧用现代工具重新演绎，既是对经典的致敬，也是一次绝佳的编程实践。

2. 核心原理拆解：大衍术的现代数学翻译

在动手写代码之前，我们必须吃透“大衍求一术”到底在算什么。用现代数学语言重新表述，它的核心是解决如下形式的一次同余式组：

x ≡ r1 (mod m1) x ≡ r2 (mod m2) ... x ≡ rn (mod mn)

我们的目标是找到一个满足所有同余式的最小正整数解x。这里m1, m2, ..., mn被称为“模数”或“问数”，r1, r2, ..., rn是对应的“余数”。

秦九韶的“大衍总数术”是一套完整的解决方案，而“大衍求一术”特指其中最关键的子步骤：求解乘率。什么是乘率？我们一步步来看。

2.1 从“总数术”到“求一术”

整个算法的流程可以概括为以下几个步骤，这也是我们后续代码设计的蓝图：

化约问数为定母：输入的模数mi可能不互质。秦九韶通过“两两连环求等，约奇弗约偶”的方法，将它们化为一组两两互质的“定母”m'i。这一步确保了后续逆元存在的必要条件（模数互质），是整个算法能成立的基础。现代实现中，我们通常直接要求输入模数两两互质，或先调用一个预处理函数来化约。
计算衍母和衍数：
- 衍母 (M)：所有定母的乘积，M = m'1 * m'2 * ... * m'n。最终的解会模M。
- 衍数 (Mi)：对于每个定母m'i，其衍数Mi = M / m'i。即除去自身后，其他所有定母的乘积。
计算奇数：用每个定母m'i去除以其对应的衍数Mi，得到的余数di = Mi % m'i，称为“奇数”。如果di为1，事情就简单了。
核心：大衍求一术求乘率：对于每个di != 1的情况，我们需要求解一个数xi，使得xi * di ≡ 1 (mod m'i)。这个xi就是“乘率”。“求一”之名，正来源于此——求得一个乘数，使得其与奇数之积模定母余一。这本质上就是求解di在模m'i下的模逆元 (Modular Multiplicative Inverse)。
合成最终解：
- 用数：Ni = Mi * xi。
- 各总：Si = ri * Ni。
- 总数：S = Σ Si。
- 所求数：x = S % M。如果x为0，通常取M作为最小正整数解。

可以看到，第4步“求乘率”是算法的核心与难点，也是“大衍求一术”得名的原因。一旦求出所有乘率，剩下的就是简单的乘加运算。

2.2 “求一术”的算法：更相减损的智慧

秦九韶给出的求乘率的方法，是一种基于辗转相除（更相减损）的迭代算法。它比朴素的枚举法高效得多，其思想与扩展欧几里得算法 (Extended Euclidean Algorithm) 异曲同工，都是为了求解方程a * x ≡ 1 (mod m)（其中a是奇数di，m是定母m'i）。

算法步骤（对照秦九韶的“置奇右上，定居右下…”）：

初始化四个变量，我们记为top_left,top_right,bottom_left,bottom_right。
- top_left = 1（“立天元一于左上”）
- top_right = a（奇数，置右上）
- bottom_left = 0
- bottom_right = m（定母，置右下）
循环执行以下操作，直到top_right == 1（“须使右上末后奇一而止”）： a. 计算商q = bottom_right / top_right。 b. 更新右下角：new_bottom_right = bottom_right % top_right。 c. 更新左下角：new_bottom_left = bottom_left + q * top_left。 d. 将右上、右下、左上、左下进行“旋转”赋值，为下一轮迭代做准备： *bottom_right = top_right*bottom_left = top_left*top_right = new_bottom_right*top_left = new_bottom_left（注意：这个“旋转”是理解的关键，它相当于把原来右列作为新的被除数（移到底部），余数作为新的除数（移到顶部））
循环结束时，top_left的值即为所求的乘率x（“乃验左上所得，以为乘率”）。

注意：这个算法描述的是当a和m互质时，最终top_right会变为1。代码实现中必须处理迭代过程，并注意在top_right变为0时的情况（即a与m不互质，此时逆元不存在，应报错）。

为什么这样就能求出逆元？这个迭代过程本质上是在记录扩展欧几里得算法中的系数。当算法终止于top_right == 1时，我们有关系式：top_left * a + bottom_left * m = 1由于bottom_left * m是m的倍数，所以在模m的意义下，上式变为：top_left * a ≡ 1 (mod m)因此，top_left就是a模m的逆元。

理解了这个原理，我们就能用C++将其精确地表达出来。接下来，我们将进入具体的代码设计与实现环节。

3. C++ 类设计与实现

为了将大衍求一术清晰地封装起来，我们设计一个DayanSolver类。这个类负责存储问题参数、执行算法步骤，并最终给出解。使用类的优势在于可以保持状态、复用代码，并且接口清晰。

3.1 数据结构与接口设计

首先，我们需要决定如何表示可能的大整数。对于学习和小型问题，C++内置的int64_t(long long) 通常足够。但为了算法的通用性和致敬古人处理大数的智慧，我们应当考虑使用大整数库，如boost::multiprecision::cpp_int。这里为了代码的纯粹性和可移植性，我们先使用int64_t，但会在架构上留出扩展接口。

#include <vector> #include <stdexcept> #include <cstdint> #include <numeric> // for std::gcd (C++17) #include <iostream> class DayanSolver { private: // 使用 vector 存储同余式组：余数 r_i 和模数 m_i std::vector<int64_t> remainders; std::vector<int64_t> moduli; // 即“问数”，要求输入时两两互质，或后续化约 std::vector<int64_t> reduced_moduli; // 化约后的“定母” int64_t M; // 衍母 std::vector<int64_t> multipliers; // 乘率 x_i int64_t solution; // 最终的最小正整数解 // 内部工具函数 bool arePairwiseCoprime(const std::vector<int64_t>& nums); int64_t computeModularInverse(int64_t a, int64_t m); // 大衍求一术核心 void reduceModuli(); // 化约问数为定母 (简化版) public: // 构造函数：传入余数和模数数组 DayanSolver(const std::vector<int64_t>& r, const std::vector<int64_t>& m); // 主求解函数 void solve(); // 获取结果 int64_t getSolution() const { return solution; } int64_t getDerivedModulus() const { return M; } // 获取衍母 void printSolutionSteps(std::ostream& os = std::cout) const; // 打印求解步骤 };

3.2 核心算法实现：求一术与总数术

接下来是核心算法的实现。我们先实现化约定母的简化版本（假设输入已互质或简单处理），然后重点实现求模逆元的函数。

// 构造函数，进行基础检查 DayanSolver::DayanSolver(const std::vector<int64_t>& r, const std::vector<int64_t>& m) : remainders(r), moduli(m) { if (remainders.size() != moduli.size() || remainders.empty()) { throw std::invalid_argument("Remainders and moduli must have the same non-zero size."); } for (auto mi : moduli) { if (mi <= 0) { throw std::invalid_argument("All moduli must be positive integers."); } } // 简化处理：这里先直接拷贝，假设输入模数已互质。 // 一个健壮的实现应调用 reduceModuli()。 reduced_moduli = moduli; } // 检查一组数是否两两互质 bool DayanSolver::arePairwiseCoprime(const std::vector<int64_t>& nums) { for (size_t i = 0; i < nums.size(); ++i) { for (size_t j = i + 1; j < nums.size(); ++j) { if (std::gcd(nums[i], nums[j]) != 1) { // C++17 return false; } } } return true; } // 大衍求一术：计算 a 在模 m 下的逆元，满足 a*x ≡ 1 (mod m) int64_t DayanSolver::computeModularInverse(int64_t a, int64_t m) { // 确保 a 和 m 互质，否则逆元不存在 if (std::gcd(a, m) != 1) { throw std::runtime_error("Modular inverse does not exist because a and m are not coprime."); } int64_t top_left = 1, top_right = a % m; // 奇数 int64_t bottom_left = 0, bottom_right = m; // 定母 // 秦九韶的迭代过程 while (top_right != 1) { if (top_right == 0) { // 理论上，由于互质，不会走到0。此处为安全。 throw std::runtime_error("Algorithm error: reached zero before one."); } int64_t q = bottom_right / top_right; // 商 int64_t new_bottom_right = bottom_right % top_right; // 新余数 int64_t new_bottom_left = bottom_left + q * top_left; // 更新左下 // 旋转赋值，为下一轮迭代准备 bottom_right = top_right; bottom_left = top_left; top_right = new_bottom_right; top_left = new_bottom_left; } // 此时 top_right == 1, top_left 即为乘率 x // 需要确保结果是正数且在 [0, m) 范围内 int64_t inverse = top_left % m; if (inverse < 0) inverse += m; return inverse; }

现在，实现完整的solve()函数，它串联起“大衍总数术”的所有步骤：

void DayanSolver::solve() { // 步骤1：检查并化约定母（简化版，直接使用输入模数，假设已互质） if (!arePairwiseCoprime(moduli)) { // 在实际完整实现中，这里应调用 reduceModuli() 进行化约。 // 为简化，我们抛出异常，要求输入互质。 throw std::runtime_error("Moduli are not pairwise coprime. A full implementation would reduce them."); } // 本例中，我们使用 reduced_moduli (初始化为 moduli) // 步骤2：计算衍母 M M = 1; for (int64_t mi : reduced_moduli) { // 防止溢出，使用 int64_t，对于大数需用大整数库 if (M > INT64_MAX / mi) { throw std::overflow_error("Product of moduli overflows int64_t. Consider using a big integer library."); } M *= mi; } // 步骤3：计算衍数 Mi 和奇数 di，并求乘率 xi std::vector<int64_t> derived_numbers; // 衍数 Mi multipliers.clear(); multipliers.resize(reduced_moduli.size()); for (size_t i = 0; i < reduced_moduli.size(); ++i) { int64_t Mi = M / reduced_moduli[i]; derived_numbers.push_back(Mi); int64_t di = Mi % reduced_moduli[i]; // 奇数 if (di == 1) { multipliers[i] = 1; // 奇数为一，乘率为一 } else { multipliers[i] = computeModularInverse(di, reduced_moduli[i]); } } // 步骤4：计算用数、各总、总数、最终解 int64_t total_sum = 0; for (size_t i = 0; i < reduced_moduli.size(); ++i) { int64_t Ni = derived_numbers[i] * multipliers[i]; // 用数 int64_t Si = remainders[i] * Ni; // 各总 total_sum += Si; } solution = total_sum % M; // 保证解为正数。如果 solution 为 0，通常问题期望的是 M 本身（最小正整数解）。 if (solution <= 0) { solution += M; } }

最后，实现一个打印详细步骤的函数，便于调试和理解：

void DayanSolver::printSolutionSteps(std::ostream& os) const { os << "=== 大衍求一术求解步骤 ===" << std::endl; os << "同余式组：" << std::endl; for (size_t i = 0; i < remainders.size(); ++i) { os << " x ≡ " << remainders[i] << " (mod " << moduli[i] << ")" << std::endl; } os << "\n1. 定母 (假设输入模数已互质): "; for (auto m : reduced_moduli) os << m << " "; os << "\n2. 衍母 M = " << M << std::endl; os << "\n3. 求乘率过程：" << std::endl; // 这里可以重新计算并打印每个乘率的求解中间过程，为了简洁，我们只打印结果。 os << " 索引 | 定母 m' | 衍数 M_i | 奇数 d_i | 乘率 x_i" << std::endl; os << " ----|---------|----------|----------|---------" << std::endl; for (size_t i = 0; i < reduced_moduli.size(); ++i) { int64_t Mi = M / reduced_moduli[i]; int64_t di = Mi % reduced_moduli[i]; os << " " << i << " | " << reduced_moduli[i] << " | " << Mi << " | " << di << " | " << multipliers[i] << std::endl; } os << "\n4. 合成最终解：" << std::endl; os << " 总数 S = Σ(r_i * M_i * x_i) = "; int64_t total = 0; for (size_t i = 0; i < reduced_moduli.size(); ++i) { int64_t Mi = M / reduced_moduli[i]; int64_t term = remainders[i] * Mi * multipliers[i]; os << (i > 0 ? " + " : "") << "(" << remainders[i] << "*" << Mi << "*" << multipliers[i] << ")"; total += term; } os << " = " << total << std::endl; os << " 最小正整数解 x = S mod M = " << total << " mod " << M << " = " << solution << std::endl; }

4. 完整示例、测试与进阶优化

有了核心类，我们写一个main函数来测试经典的“物不知数”问题。

4.1 基础测试与验证

int main() { // 测试《孙子算经》问题： x ≡ 2 (mod 3), ≡ 3 (mod 5), ≡ 2 (mod 7) std::vector<int64_t> remainders = {2, 3, 2}; std::vector<int64_t> moduli = {3, 5, 7}; try { DayanSolver solver(remainders, moduli); solver.solve(); solver.printSolutionSteps(); std::cout << "\n验证：" << std::endl; int64_t x = solver.getSolution(); for (size_t i = 0; i < moduli.size(); ++i) { std::cout << " " << x << " % " << moduli[i] << " = " << (x % moduli[i]) << ", 期望余数 " << remainders[i]; if (x % moduli[i] == remainders[i]) { std::cout << " ✓" << std::endl; } else { std::cout << " ✗" << std::endl; } } } catch (const std::exception& e) { std::cerr << "求解错误: " << e.what() << std::endl; return 1; } // 测试另一个例子 std::cout << "\n\n=== 测试例子2 ===" << std::endl; std::vector<int64_t> r2 = {1, 2, 3}; std::vector<int64_t> m2 = {3, 5, 7}; // 模数互质 DayanSolver solver2(r2, m2); solver2.solve(); std::cout << "解为: " << solver2.getSolution() << std::endl; // 验证: 52 % 3 = 1, 52 % 5 = 2, 52 % 7 = 3 return 0; }

运行这段代码，你会看到详细的求解步骤和验证结果，最终输出23，与历史答案一致。

4.2 处理非互质模数：化约定母的初步实现

上面的实现假设输入模数两两互质。但秦九韶的原术包含了“约奇弗约偶”等化约步骤来处理非互质情况。这是一个更复杂的部分，其目标是将一组可能不互质的模数{m_i}，转化为一组两两互质的{m'_i}，同时调整余数{r_i}，使得新方程组与原方程组等价。

简化版的思路是：对于每对模数(m_i, m_j)，检查其最大公约数g = gcd(m_i, m_j)。如果g != 1，则需要确保对应的余数r_i和r_j在模g下同余（即r_i % g == r_j % g），否则原方程组无解。然后，可以通过约化模数来简化问题。一个常见的策略是使用模数的质因数分解，但实现起来较为复杂。

这里提供一个非常简化的、针对特定情况的处理思路，仅作示意，完整的通用化约算法需要更细致的数论处理：

void DayanSolver::reduceModuli() { // 这是一个示意性的简化实现，并非完整的“大衍总数术”化约。 // 完整实现需处理“两两连环求等，约奇弗约偶”等复杂规则。 reduced_moduli.clear(); // 简单策略：如果所有模数互质，直接使用。 // 否则，这里可以抛出一个异常，提示用户输入需互质，或引导使用更复杂的库。 if (!arePairwiseCoprime(moduli)) { throw std::runtime_error("Non-coprime moduli detected. The simplified solver requires pairwise coprime moduli. Consider pre-processing the equations."); } reduced_moduli = moduli; }

实操心得：在实际工程中，如果遇到模数不互质的情况，更常见的做法是先检查方程组是否有解（利用所有两两模数的最大公约数检验余数一致性），如果无解则直接报错；如果有解，可以尝试将方程组拆分为多个互质的子方程组分别求解，再用中国剩余定理的推广形式合并。对于学习目的，要求输入互质是合理的简化。

4.3 性能优化与扩展思考

大整数支持：我们的实现基于int64_t，乘积M很容易溢出。对于实际应用，必须使用大整数库，如boost::multiprecision::cpp_int或GMP。只需将代码中的int64_t替换为cpp_int，并移除溢出检查即可。
求逆元算法的选择：我们实现了秦九韶的“求一术”，它与扩展欧几里得算法本质相同。在标准库中，C++17 并未直接提供模逆元函数，但我们可以用std::gcd和扩展欧几里得算法实现一个更常见的版本。两者时间复杂度都是O(log min(a, m))，效率相当。选择“求一术”的实现主要是为了贴合主题。
```
// 使用扩展欧几里得算法求逆元 (替代实现) int64_t modInverseEEA(int64_t a, int64_t m) { int64_t m0 = m, t, q; int64_t x0 = 0, x1 = 1; if (m == 1) return 0; while (a > 1) { q = a / m; t = m; m = a % m; a = t; t = x0; x0 = x1 - q * x0; x1 = t; } if (x1 < 0) x1 += m0; return x1; }
```
批量求解优化：如果需要频繁求解不同余数、相同模数的方程组，可以预先计算所有Mi和multipliers[i]并缓存，求解时只需计算Σ(ri * Mi * xi) mod M，这可以大幅提升性能。
错误处理与鲁棒性：生产代码需要更完善的错误处理，包括输入验证、无解判断、溢出检查等。例如，在computeModularInverse中，我们假设a与m互质，但应在调用前或函数内部明确检查gcd(a, m) == 1。

5. 常见问题与调试技巧

在实现和使用这个大衍求一术求解器时，你可能会遇到以下几个典型问题：

5.1 问题一：程序输出错误或崩溃

可能原因1：整数溢出。
- 排查：检查衍母M的计算。如果模数较大，它们的乘积很可能超出int64_t的范围。在solve()函数中我们加入了简单的溢出检查，但最好使用大整数库。
- 解决：将所有int64_t替换为boost::multiprecision::cpp_int，并链接 Boost 库。
可能原因2：模数不互质，但未正确处理。
- 排查：在构造DayanSolver后、调用solve()前，可以打印或检查模数数组。使用arePairwiseCoprime函数验证。
- 解决：确保输入的模数两两互质。如果问题本身模数不互质，你需要一个更完整的实现来处理化约，或者将问题转化为可解的子问题。
可能原因3：余数ri大于等于模数mi。
- 排查：同余式x ≡ ri (mod mi)中，通常约定0 <= ri < mi。如果ri >= mi，应先取模ri % mi规范化。
- 解决：在构造函数或solve()开始时，对remainders进行规范化处理：ri = ri % mi。

5.2 问题二：求解速度慢

可能原因：模数非常多且很大，导致：
1. 计算衍母M的乘积时，大数乘法开销大。
2. 对每个模数求逆元computeModularInverse，虽然单次是O(log n)，但数量多了也可观。
优化：
- 对于固定模数的问题，使用上述提到的预计算缓存策略。
- 考虑使用更高效的数论库，如GMP或NTL，它们对大数运算和模逆运算有高度优化。
- 检查是否真的需要所有模数。有时问题可以简化。

5.3 问题三：如何验证复杂结果的正确性？

方法1：暴力验证（小范围）：对于解x，直接循环验证x + k * M(k 为小整数) 是否满足所有同余式。因为解在模M意义下唯一，所以只需验证x和x+M等少数几个值。
方法2：使用已知结果或另一个可靠工具交叉验证。例如，可以用 Python 的sympy库的crt(Chinese Remainder Theorem) 函数来验证。
```
# Python 验证代码示例 from sympy.ntheory.modular import crt m = [3, 5, 7] r = [2, 3, 2] x, M = crt(m, r) # 返回 (解, 模) print(x) # 输出 23
```
方法3：分步打印：充分利用我们写的printSolutionSteps()函数，查看每一步的中间结果（定母、衍数、奇数、乘率），与手工计算或理论推导进行比对。

5.4 一个调试案例：乘率计算错误

假设在计算x * 65 ≡ 1 (mod 83)时，我们的computeModularInverse(65, 83)没有返回23。

调试步骤：
1. 在computeModularInverse函数内添加详细打印，输出每一轮迭代的top_left,top_right,bottom_left,bottom_right和商q。
2. 对比秦九韶原术的演算过程（见原理部分表格）。
3. 检查循环终止条件是否为top_right == 1，以及最终top_left的处理（取模确保为正）。
4. 检查输入a=65,m=83是否互质（gcd(65,83)=1）。
5. 最常见错误在于“旋转赋值”的逻辑弄反了。仔细对照代码和算法描述，确保赋值顺序正确。

实现这个近八百年前的算法，不仅是一次编程练习，更是一次与数学史的对话。通过将“大衍求一术”翻译成严谨的C++代码，我们不仅得到了一个解决特定数论问题的工具，更重要的是理解了古人如何系统化、程序化地思考数学问题。这种思想，与现代计算机算法设计的精神是相通的。当你下次再遇到需要处理模运算、同余方程的编程任务时，不妨想想这个来自南宋的智慧，它或许能给你带来不一样的启发。