zoukankan      html  css  js  c++  java
  • 扩展欧几里得 细讲

    2018/11/8 更新

    求满足 \(ax \mod b = 1\) 的最小正整数 \(x\)。这个 \(x\) 被称为“\(a\) 在模 \(b\) 意义下的逆元”。不少题目中需要求乘法逆元(虽然只是作为计算答案的辅助工具)。解法不止一种,下面用扩展欧几里得求解。

    一步一步来。

    问题转化

    题目问的是满足 \(ax \mod b = 1\) 的最小正整数 \(x\)。(\(a,b\)是正整数)

    但是不能暴力枚举 \(x\),会超时。

    把问题转化一下。观察 \(ax \mod b = 1\),它的实质\(ax + by = 1\):这里 \(y\) 是我们新引入的某个整数,并且似乎是个负数才对。这样表示是为了用扩展欧几里得算法。我们将要努力求出一组 \(x,y\) 来满足这个等式。稍微再等一下——

    问题还需要转化。扩展欧几里得是用来求 \(ax + by = gcd(a,b)\) 中的未知数的,怎么牵扯到等于 \(1\) 呢?

    原理是,方程 \(ax + by = m\) 有解的必要条件是 $m \mod gcd(a,b) = 0 $。


    这个简单证一下。

    由最大公因数的定义,可知 \(a\)\(gcd(a,b)\) 的倍数,且 \(b\)\(gcd(a,b)\) 的倍数,

    \(x,y\) 都是整数,就确定了 \(ax + by\)\(gcd(a,b)\) 的倍数,

    因为 \(m = ax + by\),所以 \(m\) 必须是 \(gcd(a,b)\) 的倍数,

    那么 \(m \mod gcd(a,b) = 0\)


    可得出在这道题中,方程 \(ax + by = 1\) 的有解的必要条件是 \(1 \mod gcd(a,b) = 0\)。可怜的 \(gcd(a,b)\) 只能等于 \(1\) 了。这实际上就是 \(a,b\) 互质。

    扩展欧几里得

    前提是知道了普通欧几里得算法(辗转相除法)。下面字母挺多,希望你耐心地慢慢地读~

    我们拿到了一组 \(a,b\)。目标是求出满足 \(ax + by = gcd(a,b)\) (①) 的整数 \(x\)\(y\)。其中 \(x\) 应当是满足条件的最小正整数,它会成为答案,\(y\) 是辅助答案。

    (注意,虽然刚刚已经证明本题的 \(gcd(a,b)\) 必然等于 \(1\)但是扩展欧几里得算法本身过程求的是 \(ax + by = gcd(a,b)\) 的解。它正好非常适合本题,不过我们要按照它求解的过程去做)

    根据普通欧几里得算法,\(gcd(a,b) = gcd(b, a\mod b)\) (②)。

    如果我们先前已经求出了另一组数 \(x_2, y_2\),它们满足 \(bx_2 + (a \mod b)y_2 = gcd(b, a \mod b)\) (③),(提示一下,这个等式与①的结构一致,系数则是根据②替换的。此处的递归形态已经有所显露,别着急

    结合②③,一定有

    \(bx_2 + (a \mod b)y_2\)

    \(=gcd(b, a \mod b)\)

    \(=gcd(a,b)\)

    在这个“如果”实现的时候,我们的目标变成了

    求出满足 \(ax + by = bx_2 + (a \mod b)y_2\) (④)的 \(x\)\(y\)

    其中\(a,b,x_2,y_2\) 都已知,\(x,y\)待求。未知数比方程更多,没有唯一解。我们先求出一组必然存在的解,最后将在答案处理时转变为使正整数 \(x\) 最小的解。

    取模运算是 \(a \mod b = a - b×(a/b)\),能理解吧?

    等式④实际上是:

    \(ax + by = bx_2 + (a-b×(a/b))y_2\)

    \(ax + by = bx_2 + ay_2 - b × (a/b)y_2\)

    \(ax + by = ay_2 + b(x_2-(a/b)y_2)\)

    看上面这个方程,一组必然存在的解出现了!

    \(x = y_2, y = x_2 - (a/b)y_2\)

    可见我们只要求出 \(x_2,y_2\),就能得出正确的\(x,y\)。问题是 \(x_2,y_2\) 怎么求。

    脑海里抛掉前面的\(x,y\),现在我们手上是 \(b,a \mod b\) 这两个系数,而目标是求出满足 \(bx_2 + (a \mod b)y_2 = gcd(b,a \mod b)\)(③)的 \(x_2\)\(y_2\),以便于待会回馈给⑤。

    等一下。

    比较于③,我再去看看原方程①:

    \(ax + by = gcd(a,b)\) (①)

    原方程中的 \(a\) 变成了③中的 \(b\),原方程中的 \(b\) 变成了③中的 \(a \mod b\) 而已。

    按照上面的一模一样下来(其实都只是推导过程),我们发现,最好有 \(x_3,y_3\) 来支撑 \(x_2, y_2\)

    再一模一样下来,我们又需要 \(x_4,y_4\) 来支撑 \(x_3, y_3\)

    ……

    这个递归中 \(a,b\) 不断变成 \(b, a \mod b\),逼近最后:

    \(b\) 将等于最早的 \(a,b\) 的最大公因数,就像普通 \(gcd\) 的结果。

    但我们再往下一层。此时由于 \(a \mod b = 0\),下一层的 \(b = 0\)

    是时候直接回馈了,我们需要一组 \(x_n,y_n\) 满足 \(a_nx_n + b_ny_n = gcd(a_n,b_n)\)

    然而本层的 \(b_n = 0\),则 \(gcd(a_n,b_n) = gcd(a_n,0) = a_n\)。那么只要等式左边取 \(x_n = 1\),这个等式就妥妥的成立了。

    最后一层回馈的 \(y_n\) 建议用 \(0\),但由于 \(b = 0\),回馈其它数值等式一定也成立。意料之外的是,这样的程序有时候会求错解。如果回馈 \(3\),那么你将收获 \(70\) 分的好成绩。如果回馈 \(-20002\),那么你将收获 \(40\) 分的好成绩。为什么呢?

    原因仅仅是,被回馈(却与 \(0\) 等效)的 \(y_n\) 在回溯时滚雪球式增长,容易数值越界。下面的代码在最后一层令 \(y = 7\),开了long long,没有出问题。

    最后一层结束后,就开始返回,直到最上层。每一层可以轻松地根据下层的 \(x_{k+1},y_{k+1}\) 求出当前层的 \(x_k, y_k\)

    小提示:

    上面的算法中,以辗转相除的方式向下递进的好处(目的)是缩小系数,直到出现一个解可以确定的情况。

    对扩展欧几里得也可以有不同的理解方式,比如可以设 G = gcd(a,b),然后类似地推导,每一层的等式右边都是G。

    答案处理

    一个重大遗留:我们将要求出来的 \(x,y\) 仅仅保证满足 \(ax + by = 1\),而 \(x\) 不一定是最小正整数解。有可能太大,也有可能是负数。这依然可以解决,那就是,\(x\) 批量地减去或加上 \(b\),能保证 \(ax + by = 1\)

    \(ax + by = 1\)

    \(ax + by + k×ba - k×ba = 1\)

    \(a(x+kb) + (y-ka)b = 1\)

    1.显然这并不会把方程中任何整数变成非整数。

    2.“加上或减去 \(b\)”不会使 \(x\) 错过任何解。可以这么理解:


    已经求出一组整数 \(x,y\) 使得 \(ax + by = 1\),也就是 \(\frac{1-ax}{b} = y\)\(y\) 是整数,可见目前 \(1-ax\)\(b\) 的倍数。

    现在想改变 \(x\) 并使得方程仍然成立。已知 \(a,b\) 互质,假若 \(x\)变化量\(\Delta x\))不是 \(b\) 的倍数,则 \(1-ax\) 的变化量(\(-a×\Delta x\))也不是 \(b\) 的倍数,这会使得 \(1-ax\) 不再是 \(b\) 的倍数,则 \(y\) 不是整数了。

    仅当 \(x\) 的变化量是 \(b\) 的倍数时,\(1-ax\) 能保持自己是 \(b\) 的倍数,此时就出现新的解了。


    因此到最后,如果 \(x\) 太小就不断加 \(b\) 直到大于等于 \(0\),太大则一直减 \(b\),直到最小正整数解。直接求模实现更快。

    代码

    推导中的所有 \(x,y\) 共用全局变量 long long x, y,传递也很方便。

    #include<bits/stdc++.h>
    using namespace std;
    
    long long x, y;//目前方程真正的解 
    
    void exgcd(long long a, long long b)
    {
    	//当前目的:求解 ax + by = gcd(a, b) 这么一个方程
    	
    	if(b == 0) //a, b不断改变的过程中,b最终必然会成为0
    	{
    		//在 b = 0 时等式还要成立? 使 x = 1, y = 0 ,必然成立 
    		x = 1;
    		y = 7; //建议返回0。不过y = 7能AC,证明了最后一个等式不受最后一个y影响
    		return;
    	} 
    	
    	exgcd(b, a % b);//把下一层系数传进去(先求下一个方程的解 )
    	
    	//现在我们已经拿到了下一个方程的解x, y
    	long long tx = x;//暂时存一下x,别丢了
    	x = y;
    	y = tx - a / b * y; 
    }
    
    int main()
    {
    	long long a, b;
    	cin >> a >> b;
    	exgcd(a, b);
        
    	while(x < 0)//我们求出来的x必然满足方程,但是不一定是最小正整数解
    		x += b;
    	x %= b;
        //上面三行是“答案处理”,这里用while只是帮助理解,建议写成 x = (x % b + b) % b;
    	printf("%lld\n", x);
    	return 0;
    }
    
    
    
    马上我就弃坑,宝贝。
  • 相关阅读:
    基于OpenVINO的端到端DL网络-A Year in Computer Vision中关于图像增强系列部分
    基于OpenVINO的端到端DL网络-Tesseract5+VS2017+win10源码编译攻略
    基于OpenVINO的端到端DL网络-Intel AI DevCloud 的申请和登陆
    OpenCV和RTSP的综合研究
    基于OpenVINO的端到端DL网络-包含目录、库目录、附加包含目录、附加库目录、附加依赖项之详解【正确的目录添加方法】
    消除临时对象
    Redis数据导入工具优化过程总结
    清楚利弊,用好内联
    STL区间成员函数及区间算法总结
    高效的使用STL
  • 原文地址:https://www.cnblogs.com/cicos/p/10258688.html
Copyright © 2011-2022 走看看