Via Ivy-Endhttp://www.ivy-end.com/archives/1010
关于这个算法，主要是参考NOIP2012 Day2 T1。即这里所讲的是求解这样一个线性模方程：\[ax\equiv 1\mod{p}\]的最小正整数解。
去年我是暴搜做的(PS+1)，当时什么都不会。今天在这里介绍两种算法，一种是我国古代数学家秦九韶发明的「大衍求一术」，还一种是著名的「扩展欧几里德算法」。
大衍求一术
首先来看一下大衍求一术。这里只介绍它的计算方法，至于证明可以参考扩展欧几里德算法。
例1：求解方程\(23x\equiv 1\mod{97}\)。
解：我们只需要列出下面这张表就可以得到求解\[\begin{matrix}23^{1} & 23^{1} & 3^{17} & 3^{17} & 1^{38}\\ 97^{0} & 5^{4} & 5^{4} & 2^{21} & 2^{21}\end{matrix}\]结果就是38。
接下来我们来理论化的表述一下这个算法的过程：
假设输入\(a,b\)满足\(a>b\)。那么我们用\(a_{n},A_{n}\)分别表示第一行的底数和奇数，\(b_{n},B_{n}\)分别表示第二行的底数和奇数，如果\(a_{i}>b_{i}\)，那么\(a_{i+1}=a_{i}\mod{b_{i}},A_{i+1}=A_{i}+B_{i}\cdot \left [ \frac{a_{i}}{b_{i}} \right ],b_{i+1}=b_{i},B_{i+1}=B_{i}\)；如果\(a_{i}<b_{i}\)则上面的结论倒过来即可。
算法结束当且仅当\(a_{i}=1\)，此时\(A_{i}\)即为所求的最小正整数解。
例2：求解方程\(97x\equiv 1\mod{23}\)。
解：我们只需要列出下面这张表就可以得到求解\[\begin{matrix}97^{1} & 5^{1} & 5^{1} & 2^{5} & 2^{5} & \\ 23^{0} & 23^{0} & 3^{4} & 3^{4} & 1^{9} & 1^{14}\end{matrix}\]结果就是14。
对于这个结果，如果1最先出现在下面一行，则需要再计算一次，而且这次计算必须使得余数是1。
假设输入\(a,b\)满足\(a<b\)。中间的步骤和之前一行，在计算过程中必然存在一个\(i\)使得\(b_{i}=1\)，此时我们只需计算\(B_{i+1}\)即可得到结果。其中\(B_{i+1}=A_{i}+B_{i}\cdot \left(a_{i} – 1\right)\)。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 |
#include <iostream> using namespace std; struct Num { int nBase, nIndex; // 分别表示底数和奇数 }; Num x, y; int main() { cin >> x.nBase >> y.nBase; x.nIndex = 1; y.nIndex = 0; // 初始化 if(x.nBase < y.nBase) { while(1) { if(x.nBase == 1 || x.nBase == 0) // 循环出口 { cout << x.nIndex << endl; break; } if(x.nBase < y.nBase) // 模拟计算过程 { int nDiv = y.nBase / x.nBase; if(y.nBase % x.nBase == 0) { nDiv--; } y.nBase %= x.nBase; y.nIndex += nDiv * x.nIndex; } else { int nDiv = x.nBase / y.nBase; if(x.nBase % y.nBase == 0) { nDiv--; } x.nBase %= y.nBase; x.nIndex += nDiv * y.nIndex; } } } else { while(1) { if(y.nBase == 1) // 出口 { int nDiv = x.nBase - 1; cout << nDiv * y.nIndex + x.nIndex << endl; // 求出结果 break; } if(x.nBase < y.nBase) // 模拟计算过程 { int nDiv = y.nBase / x.nBase; if(y.nBase % x.nBase == 0) { nDiv--; } y.nBase %= x.nBase; y.nIndex += nDiv * x.nIndex; } else { int nDiv = x.nBase / y.nBase; if(x.nBase % y.nBase == 0) { nDiv--; } x.nBase %= y.nBase; x.nIndex += nDiv * y.nIndex; } } } return 0; } |
扩展欧几里德算法
可能上面的算法对于某些人来说比较晦涩，我们下面来介绍一下扩展欧几里德算法。首先介绍一个定理：
方程\(ax+by=\gcd\left ( a,b \right )\)一定有解。
这样我们的问题就可以转化为求方程\(ax+b\cdot \left ( -y \right )=1\)，在这里，我们先求出方程\(ax+b\cdot \left ( -y \right )=\gcd\left(a,b\right)\)的解，然后只要将结果除以\(\gcd\left(a,b\right)\)就行了。
下面来推导一下扩展欧几里德算法。
我们已知\[ax+by=\gcd\left ( a,b \right )\]，且\(\gcd\left ( a,b \right )=\gcd\left(b,a\mod b \right )\)。不妨设\[bx{}'+\left ( a\mod b \right )y{}'=\gcd\left ( b,a\mod b \right )\]。此时就有\[bx{}'+\left ( a\mod b \right )y{}'=ax+by\]，展开得到\[bx{}'+\left ( a-\left [ \frac{a}{b} \right ]\cdot b \right )y{}’=ax+by\]，化简得\[ay{}'+b\left (x{}'-\left [ \frac{a}{b} \right ]\cdot y{}’ \right )=ax+by\]。因此可以得到\[x=y{}',y=x{}'-\left [ \frac{a}{b} \right ]\cdot y{}’\]。
这样我们就可以用递归来实现扩展欧几里德算法了。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 |
#include <iostream> using namespace std; typedef long long LL; LL A, B, C, X = 0, Y = 0; LL gcd(LL a, LL b); void exgcd(LL a, LL b, LL &x, LL &y); int main() { cin >> A >> B; C = gcd(A, B); exgcd(A, -B, X, Y); while(X < 0) { X += B; } // 找最小正整数 cout << X << endl; return 0; } void exgcd(LL a, LL b, LL &x, LL &y) { if(a == 0) { x = 0; y = C / b; // 边界 } else { exgcd(b % a, a, x, y); y = x; // 递推公式 x = (C - b * y) / a; } } LL gcd(LL a, LL b) // 求解最大公倍数 { if(b == 0) { return a; } else { return gcd(b, a % b); } } |
欧拉定理
\(a^{phi(n)}=1 \mod {n}\)
令\(x=a^{phi(n) -1}\mod{n}\)有\(ax=1 \mod{n}\)则x即为答案。
只需求出phi(n)就可以了。
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |
#include<iostream> #include<cmath> using namespace std; typedef unsigned int ll; typedef long long lll; lll qkpow(lll a, lll b, lll c){ lll ret=1%c; while(b) { if(b&1) ret=ret*a%c;a=a*a%c;b>>=1; } return ret; } ll a,b,n;lll rsts; ll phi(){ ll rst=1; double x=sqrt(b); for(int i=2;i<=x;++i){ if(b%i==0){ rst*=i-1; while(b%i==0){ rst*=i; b/=i; } rst/=i; } } if(b>1){rst*=(b-1);} return rst; } int main(){ cin>>a>>b;n=b; rsts=qkpow(a,phi()-1,n); cout<<rsts<<endl; return 0; } |
[Dr.Lib]Note:Math – 不定方程 by Liqueur Librazy is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.