大衍求一术C++ 代码的实现的一点补充

2000年3月傅海伦在《山东师大学报 (自然科学版)》(Vol. 15 No.1 pp.9-15) 发表 〈大衍求一术定数算法的现代计算机程序设计〉论文，他根据李继闵改良的定数算法，以FORTRAN 编写了程序。

今用其定数算法，改写”大衍求一术”C++原程序，并以”积尺寻源 (基广)”为例比较两组定数的结果。

// 大衍求一术
#include <iostream>
#include <vector>
using namespace std;
int gcd(int a, int b);
void reduce(vector<int>& a, int j, int i);

int main() {
    // 输入值
    vector<int> a = {130, 120, 110, 100, 60, 50, 25, 20};
    vector<int> r = {60,   30,  20,  30, 30, 30,  5, 10};

    int n = a.size();           // 问数个数
    cout << "问数 = ";
    for (int num : a) cout << num << " ";
    cout << endl;

    cout << "余数 = ";
    for (int num : r) cout << num << " ";
    cout << endl;

    // 求定数
    for (int i = n - 1; i > 0; --i) {
        for (int j = i - 1; j >= 0; --j) {
            reduce(a, j, i);
        }
    }
    cout << "个数 = " << n << endl;
    cout << "定数 = ";
    for (int num : a) cout << num << " ";
    cout << endl;

    // 创建所需列表
    vector<int> m(n, 0);        // 衍数
    vector<int> p(n, 0);        // 奇数
    vector<int> k(n, 0);        // 乘率
    vector<int> s(n, 0);        // 用数

    // 求衍母
    int M = 1;
    for (int num : a) {
        M *= num;               // M 为衍母
    }
    cout << "衍母 = " << M << endl;

    // 求衍数、奇数、乘率、用数
    for (int i = 0; i < n; ++i) {
        m[i] = M / a[i];        // m[i] 为第 i 个衍数
        p[i] = m[i];
        while (p[i] > a[i]) {
            p[i] -= a[i];       // p[i] 为第 i 个奇数
        }
        int h2 = p[i];          // 置奇右上
        int h1 = a[i];          // 定居右下
        int e2 = 1;             // 立天元一于左上
        int e1 = 0;
        while (h2 > 1) {        // 右上末后奇一而止
            if (h1 > h2) {      // 右行上下以少除多
                h1 -= h2;
                e1 += e2;
            } else if (h2 > h1) {
                h2 -= h1;
                e2 += e1;
            }
        }
        k[i] = e2;              // k[i] 为第 i 个乘率
        s[i] = k[i] * m[i];     // s[i] 为第 i 个用数
    }

    cout << "衍数 = ";
    for (int num : m) cout << num << " ";
    cout << endl;

    cout << "奇数 = ";
    for (int num : p) cout << num << " ";
    cout << endl;

    cout << "乘率 = ";
    for (int num : k) cout << num << " ";
    cout << endl;

    cout << "用数 = ";
    for (int num : s) cout << num << " ";
    cout << endl;

    int N = 0;
    for (int i = 0; i < n; ++i) {
        N += r[i] * s[i];
    }
    cout << "总数 = " << N << endl;
    N = N % M;                  // N 即所求数
    cout << "所求数 = " << N << endl;

    return 0;
}

// 求a, b 最大公约数 (GCD)
int gcd(int a, int b) {
    while (b) {
        int temp = a % b;
        a = b;
        b = temp;
    }
    return a;
}

// 改良定数算法
void reduce(vector<int>& a, int j, int i) {
    int ig = gcd(a[j], a[i]);
    if (ig == 1) {
        return;
    } else {
        a[j] /= ig;              // j 为「奇」，约奇不约偶
    }
    while (ig > 1) {             // 反约
        ig = gcd(a[j], a[i]);
        a[j] *= ig;
        a[i] /= ig;              // i 为「偶」，约偶不约奇
    }
}

运行结果：

问数 = 130 120 110 100 60 50 25 20 
余数 = 60 30 20 30 30 30 5 10 
个数 = 8
定数 = 13 8 11 1 3 1 25 1 
衍母 = 85800
衍数 = 6600 10725 7800 85800 28600 85800 3432 85800 
奇数 = 9 5 1 1 1 1 7 1 
乘率 = 3 5 1 1 1 1 18 1 
用数 = 19800 53625 7800 85800 28600 85800 61776 85800 
总数 = 10125630
所求数 = 1230

原程序求得的定数：

定数 = 13 24 11 25 1 1 1 1

改良版本明显较为优胜, 因为定数的数字较小，令之后计算衍数、用数、总数的工作量降低。当然，两个程序的所求数结果是一样的。

究竟谁更接近作者的原意呢？看一看《数书九章》原著便知道，查《四库全书》电子版：

由此可见，改良版定数跟原著一样。