【算法笔记】第5章数学问题欧几里得算法与扩展欧几里得算法

最新推荐文章于 2024-07-20 09:51:47 发布

Nicole_177

最新推荐文章于 2024-07-20 09:51:47 发布

阅读量383

点赞数 1

分类专栏：算法文章标签：算法

算法专栏收录该内容

4 篇文章 0 订阅

订阅专栏

【算法笔记】入门篇第五章数学问题欧几里得算法与扩展欧几里得算法

欧几里得算法
- - 整数定理与推论
  - 欧几里得算法证明
扩展欧几里得算法
一些题目
Reference

欧几里得算法

	int gcd(int a, int b) {  
	    return !b ? a : gcd(b, a % b);  
	}

整数定理与推论

欧几里得算法证明

扩展欧几里得算法

例子

非递归代码

	int exgcd(int a,int b,int &x,int &y)
	{
	    int xi_1,yi_1,xi_2,yi_2;
	    xi_2=1,yi_2=0;
	    xi_1=0,yi_1=1;
	    x=0,y=1;
	    int r=a%b;
	    int q=a/b;
	    while (r)
	    {
	        x=xi_2-q*xi_1;
	        y=yi_2-q*yi_1;
        
	        xi_2=xi_1;
	        yi_2=yi_1;
	 
	        xi_1=x,yi_1=y;
	        a=b;
	        b=r;
	        r=a%b;
	        q=a/b;
	    }
	    return b;
	}

扩展欧几里得算法；方程 ax + by = gcd(a, b) 的求解

给定2个非零整数a、b，求一组整数解(x, y)，使得ax+ by = gcd(a,b)成立。（通过相关定理知解一定存在）

可以利用欧几里得算法的过程来计算x和y（推倒过程略），代码：

	int exGcd(int a, int b, int &x, int &y) {  
	    if(b == 0) {  
	        x = 1;  //递归到最底层时，给x和y赋值  
	        y = 0;  
	        return a;  
	    }  
	    int g = exGcd(b, a % b, x, y);  
	    int temp = x;  
	    x = y;  
	    y = temp - a / b * y;  
	    return g;  
	}

exGcd函数结束后便得到一对解，可通过下面的下面的式子得到全部解（过程略）：

	x' = x + b / gcd * k            
	y' = y - a / gcd * k           (k为任意整数)

对任意整数来说，（ x % (b / gcd) + (b / gcd)） % (b / gcd)为对应的最小非负整数解。

特殊地，如果gcd == 1，即 ax + by= 1 时，全部解的公式化简为：

	x' = x + b  * k            
	y' = y - a  * k           (k为任意整数)

（ x % b + b） % b 为对应的最小非负整数解。

方程 ax + by = c 的求解

ax + by = c 存在解的充要条件是 c % gcd == 0

假设 ax + by = gcd 的一组解为(x0，y0)，则方程 ax + by = c 的一组解为 (cx0 / gcd，cy0 / gcd)

全部解：

	    x' = x + b / gcd * k = cx0 / gcd + b / gcd * k 
	    y' = y - a / gcd * k = cy0 / gcd - a / gcd * k            (k为任意整数)

对任意整数来说，( x % (b / gcd) + (b / gcd)） % (b / gcd) 为对应的最小非负整数解， x 可取值为 cx0 / gcd

特殊地，如果gcd == 1，全部解的公式化简为：

	x' = x + b  * k = cx0 + bk           
	y' = y - a  * k  = cy0 - ak           (k为任意整数)

( x % b + b) % b为对应的最小非负整数解。

同余式 ax ≡ c(mod m) 的求解

“什么是同余式?”

对整数a,b,m来说，如果m整除a-b（即：(a - b）% m = 0），那么就说a与b模m同余，
对应的同余式为a ≡ b(mod m)，m称为同余式的模。例如：10 ≡ 13(mod 3) 10 ≡ 1(mod 3)
显然每一个整数都各自与[0,m)中的唯一的整数同余。

此处要解决的就是同余式 ax ≡ c(mod m)的求解，根据同余式的定义，有(ax - c)%m = 0成立，
因此存在整数y，使得ax - c = my成立。移项并令y = -y后即得：ax + by = c 我们已知这个方程的解法。

结论：

设 a,c,m是整数，其中 m >= 1，则：

(1)若 c % gcd(a, m) != 0，则同余式方程ax ≡ c(mod m)无解
(2)若c % gcd(a, m) == 0，则同余式方程ax ≡ c(mod m)恰好有gcd(a,m)个模m意义下的不同的解，且解的形式为 x' = x + m / gcd(a, m) * k，k = 0，1，...，gcd(a，m)-1，x是一个解。

逆元的求解以及(b / a)%m 计算

“什么是逆元(乘法逆元)？”

对整数a,b,m来说，m>1，且ab≡1(mod m)成立，那么就说a和b互为模m逆元，
一般也记做a ≡ 1 / b (mod m) 或 b ≡ 1 / a (mod m)，通俗的说，如果两个整数的乘积模m后等于1，就称他们互为模m逆元。

逆元的作用：

对乘法来说，有 (b * a) % m = ( (b % m) * (a % m) ) % m 成立，
但是对除法来讲，(b / a) % m = ( (b % m) / (a % m) ) % m却不成立，(b / a) % m = ( (b % m) / a ) % m也不成立。

此时就需要逆元来计算(b / a) % m。通过找到a模m的逆元x，就有(b / a)% m = (b * x)% m成立（仅考虑整数取模，也即假设 b % a = 0，即b是a的整数倍），于是就把除法取模转换为乘法取模，这对于解决被除数b非常大(即使b已经取过模)的问题来讲是很有用的。

有定义知：求a模m的逆元，就是求解同余式 ax ≡ 1(mod m)，并且在使用中一般把x的最小整数解称为a模m的逆元，下同，显然同余式ax ≡ 1(mod m)是否有解取决于1%gcd(a,m)是否为0，而这等价于gcd(a,m)是否为1：
(1)若gcd(a,m) != 1，那么同余式ax ≡ 1(mod m)无解，a不存在m的逆元。
(2) 若gcd(a,m) == 1，那么同余式ax ≡ 1(mod m)在(0，m)上有唯一解，可以通过求解ax + my = 1得到。
注意：由于gcd(a,m) == 1，因此ax + my = 1 = gcd(a,m)，直接使用扩展欧几里得算法解出x之后就可以用（x%m + m） % m得到（0，m）范围内的解，即所需要的逆元。

下面代码使用扩展欧几里得算法a模m的逆元，使用条件是gcd(a，m) = 1，当然如果m是素数，就肯定成立了

	int inverse(int a, int m) {  
	    int x, y;  
	    int g = exGcd(a, m, x, y);  //求解 ax + my = 1  
	    return (x % m + m) % m;  //a模m的逆元为(x % m + m) % m  
	}

费马小定理

如果m是素数，且a不是m的素数，则还可以直接用费马小定理来得到逆元，这种做法不需要使用扩展欧几里得算法：

a ^ (m-2) % m 就是 a%m的逆元，快速幂求解。

当gcd(a,m) != 1时，扩展欧几里得算法和费马小定理均会失效，（过程略），得：

在a和m有可能不互素的情况下，可以使用公式 (b / a) % m = (b % (am)) / a来计算(b / a) % m的值，但是当a和m的乘积太大时可能溢出…

一些题目

ZOJ 3609 ：http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemId=4712 求最小逆元

	#include <iostream>
	#include <cstdio>
	#include <cstring>
	#include <cmath>
	#include <vector>
	#include <string>
	#include <queue>
	#include <stack>
	#include <algorithm>
	 
	#define INF 0x7fffffff
	#define EPS 1e-12
	#define MOD 1000000007
	#define PI 3.141592653579798
	#define N 100000
	 
	using namespace std;
	 
	typedef long long LL;
	typedef double DB;
	 
	LL e_gcd(LL a,LL b,LL &x,LL &y)
	{
	    if(b==0)
	    {
	        x=1;
	        y=0;
	        return a;
	    }
	    LL ans=e_gcd(b,a%b,x,y);
	    LL temp=x;
	    x=y;
	    y=temp-a/b*y;
	    return ans;
	}
	 
	LL cal(LL a,LL b,LL c)
	{
	    LL x,y;
	    LL gcd=e_gcd(a,b,x,y);
	    if(c%gcd!=0) return -1;
	    x*=c/gcd;
	    b/=gcd;
	    if(b<0) b=-b;
	    LL ans=x%b;
	    if(ans<=0) ans+=b;
	    return ans;
	}
	 
	int main()
	{
	    LL a,b,t;
	    scanf("%lld",&t);
	    while(t--)
	    {
	        scanf("%lld%lld",&a,&b);
	        LL ans=cal(a,b,1);
	        if(ans==-1) printf("Not Exist\n");
	        else printf("%lld\n",ans);
	    }
	    return 0;
	}

ZOJ 3593 http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=3593 求最小的步数，处理特殊一点就过去了

	#include <iostream>
	#include <cstdio>
	#include <cstring>
	#include <cmath>
	#include <string>
	#include <vector>
	#include <stack>
	#include <queue>
	#include <algorithm>
	 
	#define INF 0x7fffffff
	#define EPS 1e-12
	#define MOD 100000007
	#define PI 3.14159265357979823846
	#define N 100005
	 
	using namespace std;
	 
	typedef long long LL;
	 
	LL e_gcd(LL a,LL b,LL &x,LL &y)
	{
	    if(b==0)
	    {
	        x=1;
	        y=0;
	        return a;
	    }
	    LL ans=e_gcd(b,a%b,x,y);
	    LL temp=x;
	    x=y;
	    y=temp-a/b*y;
	    return ans;
	}
	LL cal(LL a,LL b,LL L)
	{
	    LL x,y;
	    LL gcd=e_gcd(a,b,x,y);
	    if(L%gcd!=0) return -1;
	    x*=L/gcd;
	    y*=L/gcd;
	    a/=gcd;
	    b/=gcd;
	    LL ans=((LL)INF)*((LL)INF), f;
	    LL mid=(y-x)/(a+b);
	    for(LL T=mid-1;T<=mid+1;T++)
	    {
	        if(abs(x+b*T)+abs(y-a*T)==abs(x+b*T+y-a*T))
	            f=max(abs(x+b*T),abs(y-a*T));
	        else
	            f=fabs(x-y+(a+b)*T);
	        ans=min(ans,f);
	    }
	    return ans;
	}
	 
	int main()
	{
	    //freopen("in.in","r",stdin);
	    //freopen("out.out","w",stdout);
	    LL A,B,a,b,x,y;
	    int t; scanf("%d",&t);
	    while(t--)
	    {
	        scanf("%lld%lld%lld%lld",&A,&B,&a,&b);
	        LL L=B-A;
	        LL ans=cal(a,b,L);
	        if(ans==-1) printf("-1\n");
	        else printf("%lld\n",ans);
	    }
	    return 0;
	}

POJ 1061 http://poj.org/problem?id=1061 青蛙的约会，裸的扩展欧几里得

	#include <iostream>
	#include <cstdio>
	#include <cstring>
	#include <cmath>
	#include <vector>
	#include <string>
	#include <queue>
	#include <stack>
	#include <algorithm>
	 
	#define INF 0x7fffffff
	#define EPS 1e-12
	#define MOD 1000000007
	#define PI 3.141592653579798
	#define N 100000
	 
	using namespace std;
	 
	typedef long long LL;
	typedef double DB;
	 
	LL e_gcd(LL a,LL b,LL &x,LL &y)
	{
	    if(b==0)
	    {
	        x=1;
	        y=0;
	        return a;
	    }
	    LL ans=e_gcd(b,a%b,x,y);
	    LL temp=x;
	    x=y;
	    y=temp-a/b*y;
	    return ans;
	}
	 
	LL cal(LL a,LL b,LL c)
	{
	    LL x,y;
	    LL gcd=e_gcd(a,b,x,y);
	    if(c%gcd!=0) return -1;
	    x*=c/gcd;
	    b/=gcd;
	    if(b<0) b=-b;
	    LL ans=x%b;
	    if(ans<=0) ans+=b;
	    return ans;
	}
	 
	int main()
	{
	    LL x,y,m,n,L;
	    while(scanf("%lld%lld%lld%lld%lld",&x,&y,&m,&n,&L)!=EOF)
	    {
	        LL ans=cal(m-n,L,y-x);
	        if(ans==-1) printf("Impossible\n");
	        else printf("%lld\n",ans);
	    }
	    return 0;
	}

HDU 1576 http://acm.hdu.edu.cn/showproblem.php?pid=1576 做点处理即可

	#include <iostream>
	#include <cstdio>
	#include <cstring>
	#include <cmath>
	#include <vector>
	#include <string>
	#include <queue>
	#include <stack>
	#include <algorithm>
	 
	#define INF 0x7fffffff
	#define EPS 1e-12
	#define MOD 1000000007
	#define PI 3.141592653579798
	#define N 100000
	 
	using namespace std;
	 
	typedef long long LL;
	typedef double DB;
	 
	LL e_gcd(LL a,LL b,LL &x,LL &y)
	{
	    if(b==0)
	    {
	        x=1;
	        y=0;
	        return a;
	    }
	    LL ans=e_gcd(b,a%b,x,y);
	    LL temp=x;
	    x=y;
	    y=temp-a/b*y;
	    return ans;
	}
	 
	LL cal(LL a,LL b,LL c)
	{
	    LL x,y;
	    LL gcd=e_gcd(a,b,x,y);
	    if(c%gcd!=0) return -1;
	    x*=c/gcd;
	    b/=gcd;
	    if(b<0) b=-b;
	    LL ans=x%b;
	    if(ans<=0) ans+=b;
	    return ans;
	}
	 
	int main()
	{
	    LL n,b,t;
	    scanf("%I64d",&t);
	    while(t--)
	    {
	        scanf("%I64d%I64d",&n,&b);
	        LL ans=cal(b,9973,n);
	        if(ans==-1) printf("Impossible\n");
	        else printf("%lld\n",ans);
	    }
	    return 0;
	}

HDU 2669 http://acm.hdu.edu.cn/showproblem.php?pid=2669 裸的扩展欧几里得

	#include <iostream>
	#include <cstdio>
	#include <cstring>
	#include <cmath>
	#include <vector>
	#include <string>
	#include <queue>
	#include <stack>
	#include <algorithm>
	 
	#define INF 0x7fffffff
	#define EPS 1e-12
	#define MOD 1000000007
	#define PI 3.141592653579798
	#define N 100000
	 
	using namespace std;
	 
	typedef long long LL;
	typedef double DB;
	 
	LL e_gcd(LL a,LL b,LL &x,LL &y)
	{
	    if(b==0)
	    {
	        x=1;
	        y=0;
	        return a;
	    }
	    LL ans=e_gcd(b,a%b,x,y);
	    LL temp=x;
	    x=y;
	    y=temp-a/b*y;
	    return ans;
	}
	 
	LL cal(LL a,LL b,LL c)
	{
	    LL x,y;
	    LL gcd=e_gcd(a,b,x,y);
	    if(c%gcd!=0) return -1;
	    x*=c/gcd;
	    b/=gcd;
	    if(b<0) b=-b;
	    LL ans=x%b;
	    if(ans<=0) ans+=b;
	    return ans;
	}
	 
	int main()
	{
	    LL a,b;
	    while(scanf("%I64d%I64d",&a,&b)!=EOF)
	    {
	        LL ans=cal(a,b,1);
	        if(ans==-1) printf("sorry\n");
	        else printf("%I64d %I64d\n",ans,(1-ans*a)/b);
	    }
	    return 0;
	}

Reference

https://zh.wikipedia.org/zh-hans/扩展欧几里得算法#cite_note-1
https://zhuanlan.zhihu.com/p/34930513
https://blog.csdn.net/u014634338/article/details/40210435
https://blog.csdn.net/zhjchengfeng5/article/details/7786595

Nicole_177

关注

1
点赞
踩
4

收藏

觉得还不错? 一键收藏
0
评论
【算法笔记】第5章数学问题欧几里得算法与扩展欧几里得算法

【算法笔记】入门篇第五章数学问题欧几里得算法与扩展欧几里得算法欧几里得算法整数定理与推论欧几里得算法证明扩展欧几里得算法例子非递归代码扩展欧几里得算法；方程 ax + by = gcd(a, b) 的求解方程 ax + by = c 的求解同余式 ax ≡ c(mod m) 的求解"什么是同余式?"结论：逆元的求解以及(b / a)%m 计算"什么是逆元(乘法逆元)？"逆元的作用：费马小定...
复制链接

扫一扫