大衍求一术的计算方法

Posted 2022-06-10 mooleetzi

引言

大衍求一术是秦九韶发明的一种求特殊一次同余式的方法，以下摘自百度百科

秦九韶，字道古，生活于南宋时期，自幼喜好数学，经过长期积累和苦心钻研，于公元1247年写成《数书九章》。这部中世纪的数学杰作，在许多方面都有创造，其中求解一次同余组的“大衍求一术”和求高次方程数值解的“正负开方术”，更是具有世界意义的成就。

---

 

进入正题

对于形如ax ≡ 1 (mod n)，(a,n) = 1的同余式，秦九韶的大衍求一术给出了一种优秀的解法。

注意a,n必须互素

第一步我们令k₀ = 0，k₁ = 1, r₀ = n, r₁ = a，（k叫乘率，也是最后答案）

第二步我们让 r₀ 和 r₁ 作带余除法，即 r₀ = q₂r₁ + r_2，

第三步我们用 r₁ 和 r₂ 作带余除法，即r₁ = q₃r₂ + r₃

...

直到出现r_n = 1， 这时我们带回同余方程

定义k_n = k_n-2 - q_nk_n-1，依次回代求得k_n即是同余方程的解

---

 

举个栗子

37x ≡ 1 (mod 73)

(37,73)=1,所以满足使用大衍求一术的条件

首先我们令k₀ = 0，k₁ = 1，r₀= 73，r₁=37

r₀/r₁: 73=37*1+36 =>q₂=1,r₂=36

r₁/r₂: 37=36*1+1 =>q₃=1,r₃=1

此时r₃=1，停止除法，代回求解

k₂=k₀-q₂k₁，

k₃=k₁-q₃k₂

=>k₃=2

因此37 x≡ 1 (mod 73)的解为

x ≡ 2 (mod 73)

 

---

 

代码实现

 1 int dyqy(int a, int p)
 2 
 3     if (__gcd(a, p) != 1)
 4         return -1; //不能用大衍求一术
 5     vector<int> k, r, q;
 6     k.emplace_back(0), k.emplace_back(1); //k0=0,k1=1
 7     r.emplace_back(p), r.emplace_back(a); //r0=p,r1=a
 8     q.emplace_back(-1), q.emplace_back(-1);
 9     int i = 0;
10     while (1)
11     
12         int ra = r[i], rb = r[i + 1];
13         int qtmp = ra / rb;
14         int rtmp = ra % rb;
15         q.emplace_back(qtmp);
16         r.emplace_back(rtmp);
17         ++i;
18         if (rtmp == 1)
19             break;
20     
21     int sz = r.size();
22     for (int i = 2; i < sz; i++)
23     
24         int ktmp = k[i - 2] - k[i - 1] * q[i];
25         k.emplace_back(ktmp);
26     
27     return k.back();
28 

大衍求一术