最大公约数与欧几里得算法(带扩展)

有时我们会有求最大公约数的需求，一般我们需要使用欧几里得算法。值得注意的是，在数学里面最大公约数是没有负的定义的，所以负数是不谈最大公约数的。但是在讲解欧几里得算法之前，我们先讲解一下其他方式求解最大公约数。

思路1：试除法

设数A、B，以及最大公约数C，那么必有A%C==0 && B%C==0。而基于整除的性质，我们可以直接枚举更小的那个数的因子，时间复杂度为O( √max( A , B ) )。

int x = min( A , B );
if( A % x == 0 && B % x == 0 ) return x;
for(int i=sqrt(x);i>=1;--i){
    if(A % i ==0 && B % i ==0 )return i;
}

思路2：辗转相减法

又称更相减损术，基本算法流程如下：

初始化时，令 fac = 1

A.待判定数 A B 如果均为偶数，则转而求 A=A/2 与 B=B/2 的最大公约数，fac *= 2

B.将A与B做差   C=|A-B|，如果 C = min( A , B )，那么算法结束，返回 C * fac   ，否则重复上述流程

证明：设A>B，且最大公约数为X，那么有 A = a * X，B = b * X，且 (a > b，gcd(a,b) = 1)

1. ，那么X必然是C的一个因数。

2.由于 gcd ( a , b ) = 1 并且 a > b，可设 a = r * b + p ，gcd( b ,p ) = 1，那么 C = ( ( r - 1 ) * b + r )  * X。因为 p 与 b互质，所以我们得到( ( r - 1 ) * b + p ) 与 b 互质，即 gcd ( B , C) = X = gcd ( A , B )。不过值得注意的是我们并不知道B与C的大小关系，所以编码的时候需要注意下，防止出现负数。

int gcd( int a , int b ) {
    int fac = 1;
    while( a != b ){
        if( (a & 1) == 0 ) && ( b & 1 ) == 0 ){
            a >>= 1;
            b >>= 1;
            fac <<= 1;    
        }
        if( a > b ) {
            a -= b;        
        }else{
            b -= a;
        }
    }
    return fac * a;
}

思路3：Stein算法

一般情况下欧几里得虽然复杂度可以达到O(logn)级别，但是在面对大整数时，需要手动设计除法和取模运算，这十分消耗时间。而Stein算法在辗转相减的基础上做了改进，满足了在大整数下求解最大公约数的需要。首先注意到两个结论：

1.        

2. 

Stein算法的流程如下：

1. 待判定数 A B 如果均为偶数，则转而求 A=A/2 与 B=B/2 的最大公约数，fac *= 2

2. 如果只有一个为奇数，那么直接将他变为一半，因为一奇一偶两个数不会以2作为公约数

3.若两个都是奇数，当 a = b时，返回a，否则令 a = | a - b |，b = min( a , b )，即 gcd( a , b ) = gcd ( | a - b | , min ( a , b ) )

int SteinGcd(int a,int b){
    if ( a == b ) return a;
    if( (a & 1) == 0 ) && ( b & 1 ) == 0 ){
        return SteinGcd ( a >> 1 , b >> 1 ) << 1;    
    }if( ( a & 1 ) == 0 ) ) {
        return SteinGcd ( a >> 1 , b );
    }if( ( b & 1 ) == 0 ) { 
        return SteinGcd ( a , b >> 1 );
    }else{
        if ( a > b ){
            return SteinGcd( a - b , b );    
        }else{
            return SteinGcd( a , b - a );        
        }
    }
}

思路4：辗转相除法

欧几里得算法：辗转相除法，用来求两个数的最大公约数。

它是个递归算法，递推公式为 gcd( a , b ) = gcd( b , a%b )。

gcd算法递归到最后，就会变成 a = gcd(a,b)，b = 0。

证明：设 a = x * b + r ( r < b ) 即 a % b = r ，设 d 为 a，b 的最大公约数，由 d | a，d | b，可得d | r，即我们知道了 a 与 b 的最大公约数，也是 b 与 r 的最大公约数，也就是 gcd(a,b) = gcd(b,a%b)。

具体代码：

递归方式：

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

非递归方式：用抑或做交换

int gcd(int a,int b){
    while(b^=a^=b^=a%=b);
    return a;
}

扩展欧几里得算法

贝祖公式告诉我们，对于任意一对正整数 a , b ，一定存在整数 x , y ，有。

那么当我们要求 a * m + b * n = q 时，如果 q % gcd(a,b) != 0，那么此方程无解，反之则会有无数组解。因为对于二元一次方程来说，要么无解，要么有无数组解。一个方程约束不了两个未知数。

我们来一遍推导：

首先呢，由欧几里得算法，我们有 gcd( a , b ) = gcd( b , a%b )。

由贝祖公式我们可以设 gcd( a , b ) = a * x1 + b * y1 ，gcd( b , a%b ) = b * x2 + ( a % b ) * y2 = b * x2 + ( a - a / b * b ) * y2【此处为整除】，按照不变的量a和b整理一下得到， a * x1 = a * y2 ，b * y1 = b * ( x2 - a / b * y2 )。故x1 = y2，y1 = x2 - a / b * y2。我们可以通过递归过程和引用将结果返回。

int extgcd(int a,int b,int &x,int &y){
    if(!b){
        x=1;y=0; //一般我们预设最底层的系数 x y 为 1 0
        return a;
    }
    int gcd=extgcd(b,a%b,x,y);
    int temp=x;
    x=y;
    y=temp-a/b*y;
    return gcd;
}

关于扩展欧几里得，还有个比较重要的东西，我们可以通过这个简单的得到一组x,y使得 a * x + b * y = gcd( a , b )或者是 a * x + b * y = c【可以判断有无解或者得到一组特解】。

但是刚刚也说了，二元一次方程如果有解是无数组的，那么如何通过这一组特解得到他们的通解呢。

设特解为 x0 , y0 ,通解为 x0 + m , y0 + n。

1.对于特解而言，我们有恒等式 a * x0 + b * y0 = gcd( a , b )。

2.也就是说我们需要求得一组 m , n 使得 a * m + b * n=0。a*b/gcd(a,b)=b*a/gcd(a,b)，所以我们可以令m=b/gcd(a,b)*t，n=-a/gcd(a,b)*t。于是通解可以表示为x=x0+b/gcd(a,b)*t，y=y0-a/gcd(a,b)*t。在计算机中，负数mod正数，如果负数可以整除这个正数，返回0，否则就反悔满足同余的情况下的最大负数。所以说如果要求最小的非负x：

如果x>=0，那么直接x%(b/gcd(a,b))。

如果x<0，那么应该用x%(b/gcd(a,b))，得到满足同余的最大负数后，再加上b/gcd(a,b)得到最小正数解。也就是x%(b/gcd(a,b))+(b/gcd(a,b))

下面证明一下上述假设：对于a * m + b * n = 0，由于等式右边为0，因此无论给两边同乘任何因子，结果都应为0(前提是a * m + b * n = 0本身是成立的)。因此我们所要求解的m0，n0应满足gcd(m0,n0) = 1，在这个基础上，我们可以得到通解为 m = factor * m0， n = factor * n0。对于a * m0 + b * n0 = 0变形可得 n0 = - (a * m0) / b，要使得(a * m0) / b最小，只能的情况是(a * m0) = gcm(a,b) = a * b / gcd(a, b)，即为a * m0是 a 和 b 的最小公倍数，于是可以得到 m0 = b / gcd(a,b)，同理也可得到 n0 = - a / gcd(a,b)。m0和n0是m和n取最小绝对值的情况。

于是整个通解便可以表示为 x = x0 + k * m，y = y0 + k * n。

#include<iostream>
using namespace std;
 
int exgcd(int A,int B,int& x1,int& y1){//x1、y1示当前层的系数   x0、y0示下层的系数
	if(B==0){
		x1=1;
		y1=0;
		return A;
	}
	int d=exgcd(B,A%B,x1,y1);
	int x0=y1,y0=x1-A/B*y1;
	x1=x0;
	y1=y0;
	return d;
}

int main(){
	int A=35,B=14,x,y;
	int d=exgcd(A,B,x,y);
	cout<<"公因数:"<<d<<endl;
	cout<<A<<"*"<<x<<"+"<<B<<"*"<<y<<"="<<d<<endl;
	//解X的最小正整数 
	cout<<"X最小正整数解:"<<(x%(B/d)+B/d)%(B/d)<<endl;
	//解Y的最小正整数 
	cout<<"Y最小正整数解:"<<(y%(A/d)+A/d)%(A/d)<<endl;
} 

 

洛谷 P1082 同余方程

提交36.47k

通过20.88k

时间限制1.00s

内存限制125.00MB

题目描述

求关于xx的同余方程 a x \equiv 1 \pmod {b}ax≡1(modb) 的最小正整数解。

输入格式

一行，包含两个正整数 a,ba,b，用一个空格隔开。

输出格式

一个正整数 x_0x0​，即最小正整数解。输入数据保证一定有解。

输入输出样例

输入 #1复制

3 10

输出 #1复制

7

说明/提示

【数据范围】

对于 40%的数据，2 ≤b≤ 1,0002≤b≤1,000；

对于 60%的数据，2 ≤b≤ 50,000,0002≤b≤50,000,000；

对于 100%的数据，2 ≤a, b≤ 2,000,000,0002≤a,b≤2,000,000,000。

#include<bits/stdc++.h>
using namespace std;

int exgcd(int A,int B,int& x1,int& y1){
	if(B==0){
		x1=1;
		y1=0;
		return A;
	}
	int d=exgcd(B,A%B,x1,y1);
	int x0=y1;
	int y0=x1-A/B*y1;
	x1=x0;
	y1=y0;
	return d;
}

int main(){
	int A,B,x,y,d;
	cin>>A>>B;
	d=exgcd(A,B,x,y);
	cout<<(x%B+B)%B;
	return 0;
}