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可能不互质。秦九韶通过“两两连环求等约奇弗约偶”的方法将它们化为一组两两互质的“定母”mi。这一步确保了后续逆元存在的必要条件模数互质是整个算法能成立的基础。现代实现中我们通常直接要求输入模数两两互质或先调用一个预处理函数来化约。计算衍母和衍数衍母 (M)所有定母的乘积M m1 * m2 * ... * mn。最终的解会模M。衍数 (Mi)对于每个定母mi其衍数Mi M / mi。即除去自身后其他所有定母的乘积。计算奇数用每个定母mi去除以其对应的衍数Mi得到的余数di Mi % mi称为“奇数”。如果di为1事情就简单了。核心大衍求一术求乘率对于每个di ! 1的情况我们需要求解一个数xi使得xi * di ≡ 1 (mod mi)。这个xi就是“乘率”。“求一”之名正来源于此——求得一个乘数使得其与奇数之积模定母余一。这本质上就是求解di在模mi下的模逆元 (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是奇数dim是定母mi。算法步骤对照秦九韶的“置奇右上定居右下…”初始化四个变量我们记为top_left,top_right,bottom_left,bottom_right。top_left 1“立天元一于左上”top_right a奇数置右上bottom_left 0bottom_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 (C17) #include iostream class DayanSolver { private: // 使用 vector 存储同余式组余数 r_i 和模数 m_i std::vectorint64_t remainders; std::vectorint64_t moduli; // 即“问数”要求输入时两两互质或后续化约 std::vectorint64_t reduced_moduli; // 化约后的“定母” int64_t M; // 衍母 std::vectorint64_t multipliers; // 乘率 x_i int64_t solution; // 最终的最小正整数解 // 内部工具函数 bool arePairwiseCoprime(const std::vectorint64_t nums); int64_t computeModularInverse(int64_t a, int64_t m); // 大衍求一术核心 void reduceModuli(); // 化约问数为定母 (简化版) public: // 构造函数传入余数和模数数组 DayanSolver(const std::vectorint64_t r, const std::vectorint64_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::vectorint64_t r, const std::vectorint64_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::vectorint64_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) { // C17 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::vectorint64_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::vectorint64_t remainders {2, 3, 2}; std::vectorint64_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::vectorint64_t r2 {1, 2, 3}; std::vectorint64_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并移除溢出检查即可。求逆元算法的选择我们实现了秦九韶的“求一术”它与扩展欧几里得算法本质相同。在标准库中C17 并未直接提供模逆元函数但我们可以用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 问题二求解速度慢可能原因模数非常多且很大导致计算衍母M的乘积时大数乘法开销大。对每个模数求逆元computeModularInverse虽然单次是O(log n)但数量多了也可观。优化对于固定模数的问题使用上述提到的预计算缓存策略。考虑使用更高效的数论库如GMP或NTL它们对大数运算和模逆运算有高度优化。检查是否真的需要所有模数。有时问题可以简化。5.3 问题三如何验证复杂结果的正确性方法1暴力验证小范围对于解x直接循环验证x k * M(k 为小整数) 是否满足所有同余式。因为解在模M意义下唯一所以只需验证x和xM等少数几个值。方法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。调试步骤在computeModularInverse函数内添加详细打印输出每一轮迭代的top_left,top_right,bottom_left,bottom_right和商q。对比秦九韶原术的演算过程见原理部分表格。检查循环终止条件是否为top_right 1以及最终top_left的处理取模确保为正。检查输入a65,m83是否互质gcd(65,83)1。最常见错误在于“旋转赋值”的逻辑弄反了。仔细对照代码和算法描述确保赋值顺序正确。实现这个近八百年前的算法不仅是一次编程练习更是一次与数学史的对话。通过将“大衍求一术”翻译成严谨的C代码我们不仅得到了一个解决特定数论问题的工具更重要的是理解了古人如何系统化、程序化地思考数学问题。这种思想与现代计算机算法设计的精神是相通的。当你下次再遇到需要处理模运算、同余方程的编程任务时不妨想想这个来自南宋的智慧它或许能给你带来不一样的启发。
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可能不互质。秦九韶通过“两两连环求等约奇弗约偶”的方法将它们化为一组两两互质的“定母”mi。这一步确保了后续逆元存在的必要条件模数互质是整个算法能成立的基础。现代实现中我们通常直接要求输入模数两两互质或先调用一个预处理函数来化约。计算衍母和衍数衍母 (M)所有定母的乘积M m1 * m2 * ... * mn。最终的解会模M。衍数 (Mi)对于每个定母mi其衍数Mi M / mi。即除去自身后其他所有定母的乘积。计算奇数用每个定母mi去除以其对应的衍数Mi得到的余数di Mi % mi称为“奇数”。如果di为1事情就简单了。核心大衍求一术求乘率对于每个di ! 1的情况我们需要求解一个数xi使得xi * di ≡ 1 (mod mi)。这个xi就是“乘率”。“求一”之名正来源于此——求得一个乘数使得其与奇数之积模定母余一。这本质上就是求解di在模mi下的模逆元 (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是奇数dim是定母mi。算法步骤对照秦九韶的“置奇右上定居右下…”初始化四个变量我们记为top_left,top_right,bottom_left,bottom_right。top_left 1“立天元一于左上”top_right a奇数置右上bottom_left 0bottom_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 (C17) #include iostream class DayanSolver { private: // 使用 vector 存储同余式组余数 r_i 和模数 m_i std::vectorint64_t remainders; std::vectorint64_t moduli; // 即“问数”要求输入时两两互质或后续化约 std::vectorint64_t reduced_moduli; // 化约后的“定母” int64_t M; // 衍母 std::vectorint64_t multipliers; // 乘率 x_i int64_t solution; // 最终的最小正整数解 // 内部工具函数 bool arePairwiseCoprime(const std::vectorint64_t nums); int64_t computeModularInverse(int64_t a, int64_t m); // 大衍求一术核心 void reduceModuli(); // 化约问数为定母 (简化版) public: // 构造函数传入余数和模数数组 DayanSolver(const std::vectorint64_t r, const std::vectorint64_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::vectorint64_t r, const std::vectorint64_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::vectorint64_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) { // C17 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::vectorint64_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::vectorint64_t remainders {2, 3, 2}; std::vectorint64_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::vectorint64_t r2 {1, 2, 3}; std::vectorint64_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并移除溢出检查即可。求逆元算法的选择我们实现了秦九韶的“求一术”它与扩展欧几里得算法本质相同。在标准库中C17 并未直接提供模逆元函数但我们可以用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 问题二求解速度慢可能原因模数非常多且很大导致计算衍母M的乘积时大数乘法开销大。对每个模数求逆元computeModularInverse虽然单次是O(log n)但数量多了也可观。优化对于固定模数的问题使用上述提到的预计算缓存策略。考虑使用更高效的数论库如GMP或NTL它们对大数运算和模逆运算有高度优化。检查是否真的需要所有模数。有时问题可以简化。5.3 问题三如何验证复杂结果的正确性方法1暴力验证小范围对于解x直接循环验证x k * M(k 为小整数) 是否满足所有同余式。因为解在模M意义下唯一所以只需验证x和xM等少数几个值。方法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。调试步骤在computeModularInverse函数内添加详细打印输出每一轮迭代的top_left,top_right,bottom_left,bottom_right和商q。对比秦九韶原术的演算过程见原理部分表格。检查循环终止条件是否为top_right 1以及最终top_left的处理取模确保为正。检查输入a65,m83是否互质gcd(65,83)1。最常见错误在于“旋转赋值”的逻辑弄反了。仔细对照代码和算法描述确保赋值顺序正确。实现这个近八百年前的算法不仅是一次编程练习更是一次与数学史的对话。通过将“大衍求一术”翻译成严谨的C代码我们不仅得到了一个解决特定数论问题的工具更重要的是理解了古人如何系统化、程序化地思考数学问题。这种思想与现代计算机算法设计的精神是相通的。当你下次再遇到需要处理模运算、同余方程的编程任务时不妨想想这个来自南宋的智慧它或许能给你带来不一样的启发。
