zoukankan      html  css  js  c++  java
  • 拓展欧几里德

    例题 同余方程

    题目描述

    求关于 (x) 的同余方程 (a x equiv 1 pmod {b}) 的最小正整数解。
    输入输出格式
    输入格式:

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

    输出格式:

    一个正整数 (x_0)​,即最小正整数解。输入数据保证一定有解。

    输入输出样例
    输入样例#1:

    3 10
    输出样例#1:

    7
    说明

    【数据范围】
    对于 100%的数据,(2 ≤a, b≤ 2,000,000,000)

    NOIP 2012 提高组 第二天 第一题

    问题转化

    • 题目问的是满足 (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) 。设 (G = gcd(a, b)) 。那么目标是求出满足 (ax + by = G) (①) 的整数 (x)(y) 。其中 (x) 应当是满足条件的最小正整数,即答案,而 (y) 是辅助答案。

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

    • 如果我们先前已经求出了另一组数 (x_2, y_2) ,它们满足这么一个式子: (bx_2 + (a mod b)y_2 = G) (②),则此时结合①②一定有:

    (ax + by = bx_2 + (a mod b)y_2)​ (③)

    • 可见,在这个“如果”实现的时候,我们的目标变成了“求出满足上式的 (x)(y) ”。

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

    • 怎么求呢?取模运算是 (a mod b = a - b×(a/b)) ,所以方程③实际上是:

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

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

    (Rightarrow 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) 怎么求。

    现在我们手上是 (b,a mod b) 这两个系数,而目标是求出 (x_2)​ 和 (y_2) 满足:

    (bx_2 + (a mod b)y_2 = G)(②)

    把①和②对比一下:

    (ax + by = G) (①)

    原方程中的系数 (a) 变成了②中的系数 (b) ,原方程中的系数 (b) 变成了②中的 (a mod b)而已。

    • 所以,把新的方程也看作 (ax + by = G) 的形式(只是系数 (a)(b) 的具体数值改变了)。然后按照上面的一模一样下来(其实都只是推导过程),我们发现,最好有 (x_3,y_3) 来支撑 (x_2, y_2)

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

    ……

    这个递归中 (a, b) 不断被替代为 (b, a mod b),这个替换方式与普通欧几里得是一样的,所以最后会出现 (a_n = G, b_n = 0)

    这时要直接返回了,我们需要一组 (x_n,y_n) 满足 (a_nx_n + b_ny_n = G)(⑤)。

    然而该层的 (a_n = G, b_n = 0) 。所以只要⑤左边取 (x_n = 1) ,这个方程就妥妥的成立了。

    (最后一层的 (y_n) 建议取 0 。然而由于 (b = 0) ,就算返回其它数值,方程也一定成立。但这样的程序容易出错,因为 (y_n)​ 在回溯时滚雪球式增长,容易数值越界。

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

    整个过程就是:以辗转相除的方式向下递进,不断缩小系数,保证会出现有确定解的最后一层。
    答案处理

    一个遗留问题:我们将要求出来的 (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),也就是 ((1−ax)/b=y)(y) 是整数,可见目前 (1−ax)(b) 的倍数。

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

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

    因此到最后,如果 (x) 太小就不断加 (b) 直到大于等于 (0)$ ,太大则一直减 (b) ,直到最小正整数解。也就是这么写:

    x = (x % b + b) % b;//括号中取模再加,可以处理负数
    

    代码

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

    #include <bits/stdc++.h>
    #define int long long //全局变量 long long
    using namespace std;
    int x,y;////目前方程真正的解 
    void exgcd(int a,int b){////当前目的:求解 ax + by = gcd(a, b) 这么一个方程
    	if(b==0){ //在 b = 0 时方程还要成立, 使 x = 1, y = 0 ,必然成立 
    		x=1;y=0;//建议返回0。
    		return;
    	}
    	exgcd(b,a%b);//先求下一个方程的解 
    	//现在我们已经拿到了下一个方程的解x, y
    	int res=x;//暂时存一下x
    	x=y;
    	y=res-(a/b)*y;
    }
    signed main(){
    	int a,b;
    	scanf("%lld%lld",&a,&b);
    	exgcd(a,b);
    	printf("%lld
    ",(x%b+b)%b);//我们求出来的x必然满足方程,但不一定是最小正整数解,所以要进行答案处理
    	return 0;
    }
    
  • 相关阅读:
    apollo实现c#与android消息推送(三)
    apollo实现c#与android消息推送(二)
    apollo实现c#与android消息推送(一)
    成为架构师,需要哪些技能
    Centos 7.x Smokeping部署安装及使用
    ISIS实验配置
    Netty网编程实战:四种解决粘包方式切换、两种生产级双向监听模式并行、高效编解码、多处理器协同作战
    STP+基于LACP的portchannel 实验分享
    Java基础之反射
    IntelliJ IDEA 可以使用中文了
  • 原文地址:https://www.cnblogs.com/Lour688/p/exgcd.html
Copyright © 2011-2022 走看看