zoukankan      html  css  js  c++  java
  • 扩展欧几里得算法

    【转载】http://blog.csdn.net/qq_34494458/article/details/52637193

    一:欧几里得算法(辗转相除法)

                           基本算法:设a=qb+r,其中a,b,q,r都是整数,则gcd(a,b)=gcd(b,r),即gcd(a,b)=gcd(b,a%b)。

    证明:

           a可以表示成a = kb + r,则r = a mod b

      假设d是a,b的一个公约数,则有

      d|a, d|b,而r = a - kb,因此d|r

      因此d是(b,a mod b)的公约数

      假设d 是(b,a mod b)的公约数,则

      d | b , d |r ,但是a = kb +r

      因此d也是(a,b)的公约数

      因此(a,b)和(b,a mod b)的公约数是一样的,其最大公约数也必然相等,得证

    算法实现:

    int gcd(int a,int b) {             //递归算法  
        return b ? gcd(b, a%b) : a;  
    }  
      
    int gcd(int a, int b) {      //迭代算法  
        while(b) {  
            int c = a%b;  
            a = b;  
            b = c;  
        }  
        return a;  
    }  

    二   扩展欧几里得算法:

                   基本算法:对于不完全为 0 的非负整数 a,b,gcd(a,b)表示 a,b 的最大公约数,必然存在整数对 x,y ,使得 gcd(a,b)=ax+by。

       证明:设 a>b。

      1,显然当 b=0,gcd(a,b)=a。此时 x=1,y=0;

      2,ab!=0 时

      设 ax1+by1=gcd(a,b);

      bx2+(a mod b)y2=gcd(b,a mod b);

      根据朴素的欧几里德原理有 gcd(a,b)=gcd(b,a mod b);

      则:ax1+by1=bx2+(a mod b)y2;

      即:ax1+by1=bx2+(a-(a/b)*b)y2=ay2+bx2-(a/b)*by2;

      根据恒等定理得:x1=y2; y1=x2-(a/b)*y2;

          这样我们就得到了求解 x1,y1 的方法:x1,y1 的值基于 x2,y2.

       上面的思想是以递归定义的,因为 gcd 不断的递归求解一定会有个时候 b=0,所以递归可以结束。

    //递归代码  
    int exgcd(int a, int b, int&x, int&y) {  
       if (b == 0) {  
           x = 1;  
           y = 0;  
           return a;  
       }  
       int r = exgcd(b, a%b, y, x);  
       int t = x;  
       y = y - a/b*t;  
       return r;  
      }
       // 非递归算法  
    int exgcd(int m, int n, int &x, int &y) {  
       int x1, y1, x0, y0;  
       x0 = 1; y0 = 0;  
       x1 = 0; y1 = 1;  
       x = 0; y = 1;  
       int r = m%n;  
       int q = (m-r)/n;  
       while(r) {  
           x = x0 - q*x1;  
           y = y0 - q*y1;  
           x0 = x1; y0 = y1;  
           x1 = x; y1 = y;  
           m = n; n = r; r = m%n;  
           q = (m-r)/n;  
       }  
       return n;  
    }

    刘汝佳的:

    void gcd(int a, int b, int& d, int& x, int& y) {  
        if (!b) {d = a; x = 1; y = 0; }  
        else {gcd(b, a%b, d, y, x); y -= x*(a/b); }  
    }  

    上面求出的 x(当a>b时)都是最接近0的正整数。(不知道为什么)

    扩展欧几里德算法的应用主要有以下两个方面:

    (1)求解不定方程;

    (2)求解模线性方程(线性同余方程)与逆元;

    (1)求解不定方程;

      1.对于不定整数方程ax+by=m,若 m mod Gcd(a, b)=0,则该方程存在整数解,否则不存在整数解。

         证明:

       首先扩展欧几里德主要是用来与求解线性方程相关的问题,所以我们从一个线性方程开始分析。现在假设这个线性方程为a*x+b*y=m,如果这个线性方程有解,那么一定有gcd(a,b) | m,即a,b的最大公约数能够整除m(m%gcd(a,b)==0)。证明很简单,由于a%gcd(a,b)==b%gcd(a,b)==0,所以a*x+b*y肯定能够整除gcd(a,b),如果线性方程成立,那么就可以用m代替a*x+b*y,从而得到上面的结论,利用上面的结论就可以用来判断一个线性方程是否有解。
      2.ax+by=Gcd(a, b) 两边同时乘以 m/Gcd(a, b)得a(x*m/Gcd(a, b))+b(y*m/Gcd(a, b))=m;

    上面已经列出找一个整数解的方法,在找到 a*x+  b*y = Gcd(a, b)的一组解x0,y0后,不定方程ax+by=m的一组解满足 x = m(x0)/gcd , y = m*(y0)/gcd;
    所以通解为  x = m*(x0)/gcd + k*lcm/a , y = m*(y0)/gcd + k*lcm/b (其中k为任意整数);
    证明:

          令a1=a/gcd(a,b),b1=b/gcd(a,b),m1=m/gcd(a,b)。那么x,y的一组解就是x0*m1,y0*m1,但是由于满足方程的解无穷多个,在实际的解题中一般都会去求解x或是y的最小正数的值。以求x为例,又该如何求解呢?还是从方程入手,现在的x,y已经满足a*x+b*y=m,那么a*(x+n*b)+b*(y-n*a)=m显然也是成立的。可以得出x+n*b(n=…,-2,-1,0,1,2,…)就是方程的所有x解的集合,由于每一个x都肯定有一个y和其对应,所以在求解x的时候可以不考虑y的取值。取k使得x+k*b>0,x的最小正数值就应该是(x+k*b)%b,但是这个值真的是最小的吗??如果我们将方程最有两边同时除以gcd(a,b),则方程变为a1*x+b1*y=m1,同上面的分析可知,此时的最小值应该为(x+k*b1)%b1,由于b1<=b,所以这个值一定会小于等于之前的值。在实际的求解过程中一般都是用while(x<0)x+=b1来使得为正的条件满足,为了更快的退出循环,可以将b1改为b(b是b1的倍数),并将b乘以一个倍数后再加到x上。

    代码:

        //不定方程  
    void linear_equation(int a, int b, int c, int &x, int &y) {  
        int d = exgcd(a, b, x, y);  
        if (c%d)  
            return ;  
        int k = c/d;  
        x *= k; y *= k; // 一组解  
        return ;  
    }  

    相关题目:poj2115 poj2142 poj1061。

    (2)求解模线性方程(线性同余方程)与逆元;

            同余方程 ax≡b (mod n)对于未知数 x 有解,当且仅当 gcd(a,n) | b。且方程有解时,方程有 gcd(a,n) 个解。

               求解方程 ax≡b (mod n) 相当于求解方程 ax+ (-ny)= b, (x, y为整数)。

       //ax = b (mod n) 即 (ax) mod n = b mod n

     算法导论上有两个定理:

      定理一:设d = gcd(a, n), 假定对整数x', y', 有d = ax' + ny', 如果d | b, 则方程ax = b(mod n)有一个解的值为x0, 满足:、

          x0 = x'(b / d)(mod n)

      定理二:假设方程ax = b(mod n)有解, x0是方程的任意一个解, 则方程对模n恰有d个不同的解, 分别为: xi = x0 + i * (n / d), 其中 i = 1,2,3......d - 1

      有了这两个定理, 解方程就不难了。

    代码:

       // 模线性方程  
    bool modular_linear_equation(int a, int b, int n) {  
        int x, y, x0, i;  
        int d = exgcd(a, n, x, y);  
        if (b%d)  
            return false;  
        x0 = x*(b/d)%n; //特解  
        for(i = 0; i < d; i++)  
            printf("%d
    ", (x0 + i*(n/d))%n);  
        return true;  
    }  

    相关题目  hdu1576.

    同余方程ax≡b (mod n),如果 gcd(a,n)== 1,则方程只有唯一解。
          在这种情况下,如果 b== 1,同余方程就是 ax=1 (mod n ),gcd(a,n)= 1。
          这时称求出的 x 为 a 的对模 n 乘法的逆元。
          对于同余方程 ax= 1(mod n ), gcd(a,n)= 1 的求解就是求解方程
          ax+ ny= 1,x, y 为整数。这个可用扩展欧几里德算法求出,原同余方程的唯一解就是用扩展欧几里德算法得出的 x 。

    代码:

     
    ll inv(ll a, ll n) {  
        ll d, x, y;  
        d = exgcd(a, n, x, y);  
        return (d == 1) ? (x%n + n)%n : -1;  
    }  

    相关题目: hdu2115.

  • 相关阅读:
    Anagram
    HDU 1205 吃糖果(鸽巢原理)
    Codeforces 1243D 0-1 MST(补图的连通图数量)
    Codeforces 1243C Tile Painting(素数)
    Codeforces 1243B2 Character Swap (Hard Version)
    Codeforces 1243B1 Character Swap (Easy Version)
    Codeforces 1243A Maximum Square
    Codeforces 1272E Nearest Opposite Parity(BFS)
    Codeforces 1272D Remove One Element
    Codeforces 1272C Yet Another Broken Keyboard
  • 原文地址:https://www.cnblogs.com/pach/p/5911163.html
Copyright © 2011-2022 走看看