擴展歐幾里得算法+乘法逆元詳解


困在這個算法快一個禮拜了,在經過不斷的百度查找博客學習中終於弄懂了這個算法,並找到一個寫的非常好的大牛的博客,故特意保留下來以便以后復習

本博客轉載自:http://blog.csdn.net/zhjchengfeng5/article/details/7786595

 擴展歐幾里德算法

    誰是歐幾里德?自己百度去

    先介紹什么叫做歐幾里德算法

    有兩個數 a b,現在,我們要求 a b 的最大公約數,怎么求?枚舉他們的因子?不現實,當 a b 很大的時候,枚舉顯得那么的naïve ,那怎么做?

    歐幾里德有個十分又用的定理: gcd(a, b) = gcd(b , a%b) ,這樣,我們就可以在幾乎是 log 的時間復雜度里求解出來 a 和 b 的最大公約數了,這就是歐幾里德算法,用 C++ 語言描述如下:

    

 

    由於是用遞歸寫的,所以看起來很簡潔,也很好記憶。那么什么是擴展歐幾里德呢?

    現在我們知道了 a 和 b 的最大公約數是 gcd ,那么,我們一定能夠找到這樣的 x 和 y ,使得: a*x + b*y = gcd 這是一個不定方程(其實是一種丟番圖方程),有多解是一定的,但是只要我們找到一組特殊的解 x0 和 y0 那么,我們就可以用 x0 和 y0 表示出整個不定方程的通解:

        x = x0 + (b/gcd)*t

        y = y0 – (a/gcd)*t

    為什么不是:

        x = x0 + b*t

        y = y0 – a*t

    這個問題也是在今天早上想通的,想通之后忍不住噴了自己一句弱逼。那是因為:

    b/gcd 是 b 的因子, a/gcd 是 a 的因子是吧?那么,由於 t的取值范圍是整數,你說 (b/gcd)*t 取到的值多還是 b*t 取到的值多?同理,(a/gcd)*t 取到的值多還是 a*gcd 取到的值多?那肯定又要問了,那為什么不是更小的數,非得是 b/gcd 和a/gcd ?

    注意到:我們令 B = b/gcd , A = a、gcd , 那么,A 和 B 一定是互素的吧?這不就證明了 最小的系數就是 A 和 B 了嗎?要是實在還有什么不明白的,看看《基礎數論》(哈爾濱工業大學出版社),這本書把關於不定方程的通解講的很清楚

    現在,我們知道了一定存在 x 和 y 使得 : a*x + b*y = gcd , 那么,怎么求出這個特解 x 和 y 呢?只需要在歐幾里德算法的基礎上加點改動就行了。

    我們觀察到:歐幾里德算法停止的狀態是: a= gcd , b = 0 ,那么,這是否能給我們求解 x y 提供一種思路呢?因為,這時候,只要 a = gcd 的系數是 1 ,那么只要 b 的系數是 0 或者其他值(無所謂是多少,反正任何數乘以 0 都等於 0 但是a 的系數一定要是 1),這時,我們就會有: a*1 + b*0 = gcd

    當然這是最終狀態,但是我們是否可以從最終狀態反推到最初的狀態呢?

    假設當前我們要處理的是求出 a 和 b的最大公約數,並求出 x 和 y 使得 a*x + b*y= gcd ,而我們已經求出了下一個狀態:b 和 a%b 的最大公約數,並且求出了一組x1 和y1 使得: b*x1 + (a%b)*y1 = gcd , 那么這兩個相鄰的狀態之間是否存在一種關系呢?

    我們知道: a%b = a - (a/b)*b(這里的 “/” 指的是整除,例如 5/2=2 , 1/3=0),那么,我們可以進一步得到:

        gcd = b*x1 + (a-(a/b)*b)*y1

            = b*x1 + a*y1 – (a/b)*b*y1

            = a*y1 + b*(x1 – a/b*y1)

    對比之前我們的狀態:求一組 x 和 y 使得:a*x + b*y = gcd ,是否發現了什么?

    這里:

        x = y1

        y = x1 – a/b*y1

    以上就是擴展歐幾里德算法的全部過程,依然用遞歸寫:

    

 

    依然很簡短,相比歐幾里德算法,只是多加了幾個語句而已。

    這就是理論部分,歐幾里德算法部分我們好像只能用來求解最大公約數,但是擴展歐幾里德算法就不同了,我們既可以求出最大公約數,還可以順帶求解出使得: a*x + b*y = gcd 的通解 x 和 y

    擴展歐幾里德有什么用處呢?

    求解形如 a*x +b*y = c 的通解,但是一般沒有誰會無聊到讓你寫出一串通解出來,都是讓你在通解中選出一些特殊的解,比如一個數對於另一個數的乘法逆元

    什么叫乘法逆元?

    

    這里,我們稱 x 是 a 關於 m 的乘法逆元

    這怎么求?可以等價於這樣的表達式: a*x + m*y = 1

    看出什么來了嗎?沒錯,當gcd(a , m) != 1 的時候是沒有解的這也是 a*x + b*y = c 有解的充要條件: c % gcd(a , b) == 0

    接着乘法逆元講,一般,我們能夠找到無數組解滿足條件,但是一般是讓你求解出最小的那組解,怎么做?我們求解出來了一個特殊的解 x0 那么,我們用 x0 % m其實就得到了最小的解了。為什么?

可以這樣思考:

    x 的通解不是 x0 + m*t 嗎?

    那么,也就是說, a 關於 m 的逆元是一個關於 m 同余的,那么根據最小整數原理,一定存在一個最小的正整數,它是 a 關於m 的逆元,而最小的肯定是在(0 , m)之間的,而且只有一個,這就好解釋了。

    可能有人注意到了,這里,我寫通解的時候並不是 x0 + (m/gcd)*t ,但是想想一下就明白了,gcd = 1,所以寫了跟沒寫是一樣的,但是,由於問題的特殊性,有時候我們得到的特解 x0 是一個負數,還有的時候我們的 m 也是一個負數這怎么辦?

    當 m 是負數的時候,我們取 m 的絕對值就行了,當 x0 是負數的時候,他模上 m 的結果仍然是負數(在計算機計算的結果上是這樣的,雖然定義的時候不是這樣的),這時候,我們仍然讓 x0 對abs(m) 取模,然后結果再加上abs(m) 就行了,於是,我們不難寫出下面的代碼求解一個數 a 對於另一個數 m 的乘法逆元:

    

 

    還有最小整數解之類的問題,但都是大同小異,只要細心的推一推就出來了,這里就不一一介紹了,下面給一些題目還有AC代碼,僅供參考

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

    

[cpp]  view plain  copy
 
  1. #include <iostream>  
  2. #include <cstdio>  
  3. #include <cstring>  
  4. #include <cmath>  
  5. #include <vector>  
  6. #include <string>  
  7. #include <queue>  
  8. #include <stack>  
  9. #include <algorithm>  
  10.   
  11. #define INF 0x7fffffff  
  12. #define EPS 1e-12  
  13. #define MOD 1000000007  
  14. #define PI 3.141592653579798  
  15. #define N 100000  
  16.   
  17. using namespace std;  
  18.   
  19. typedef long long LL;  
  20. typedef double DB;  
  21.   
  22. LL e_gcd(LL a,LL b,LL &x,LL &y)  
  23. {  
  24.     if(b==0)  
  25.     {  
  26.         x=1;  
  27.         y=0;  
  28.         return a;  
  29.     }  
  30.     LL ans=e_gcd(b,a%b,x,y);  
  31.     LL temp=x;  
  32.     x=y;  
  33.     y=temp-a/b*y;  
  34.     return ans;  
  35. }  
  36.   
  37. LL cal(LL a,LL b,LL c)  
  38. {  
  39.     LL x,y;  
  40.     LL gcd=e_gcd(a,b,x,y);  
  41.     if(c%gcd!=0) return -1;  
  42.     x*=c/gcd;  
  43.     b/=gcd;  
  44.     if(b<0) b=-b;  
  45.     LL ans=x%b;  
  46.     if(ans<=0) ans+=b;  
  47.     return ans;  
  48. }  
  49.   
  50. int main()  
  51. {  
  52.     LL a,b,t;  
  53.     scanf("%lld",&t);  
  54.     while(t--)  
  55.     {  
  56.         scanf("%lld%lld",&a,&b);  
  57.         LL ans=cal(a,b,1);  
  58.         if(ans==-1) printf("Not Exist\n");  
  59.         else printf("%lld\n",ans);  
  60.     }  
  61.     return 0;  
  62. }  

ZOJ 3593 http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=3593 求最小的步數,處理特殊一點就過去了

 

 

[cpp]  view plain  copy
 
  1. #include <iostream>  
  2. #include <cstdio>  
  3. #include <cstring>  
  4. #include <cmath>  
  5. #include <string>  
  6. #include <vector>  
  7. #include <stack>  
  8. #include <queue>  
  9. #include <algorithm>  
  10.   
  11. #define INF 0x7fffffff  
  12. #define EPS 1e-12  
  13. #define MOD 100000007  
  14. #define PI 3.14159265357979823846  
  15. #define N 100005  
  16.   
  17. using namespace std;  
  18.   
  19. typedef long long LL;  
  20.   
  21. LL e_gcd(LL a,LL b,LL &x,LL &y)  
  22. {  
  23.     if(b==0)  
  24.     {  
  25.         x=1;  
  26.         y=0;  
  27.         return a;  
  28.     }  
  29.     LL ans=e_gcd(b,a%b,x,y);  
  30.     LL temp=x;  
  31.     x=y;  
  32.     y=temp-a/b*y;  
  33.     return ans;  
  34. }  
  35. LL cal(LL a,LL b,LL L)  
  36. {  
  37.     LL x,y;  
  38.     LL gcd=e_gcd(a,b,x,y);  
  39.     if(L%gcd!=0) return -1;  
  40.     x*=L/gcd;  
  41.     y*=L/gcd;  
  42.     a/=gcd;  
  43.     b/=gcd;  
  44.     LL ans=((LL)INF)*((LL)INF), f;  
  45.     LL mid=(y-x)/(a+b);  
  46.     for(LL T=mid-1;T<=mid+1;T++)  
  47.     {  
  48.         if(abs(x+b*T)+abs(y-a*T)==abs(x+b*T+y-a*T))  
  49.             f=max(abs(x+b*T),abs(y-a*T));  
  50.         else  
  51.             f=fabs(x-y+(a+b)*T);  
  52.         ans=min(ans,f);  
  53.     }  
  54.     return ans;  
  55. }  
  56.   
  57. int main()  
  58. {  
  59.     //freopen("in.in","r",stdin);  
  60.     //freopen("out.out","w",stdout);  
  61.     LL A,B,a,b,x,y;  
  62.     int t; scanf("%d",&t);  
  63.     while(t--)  
  64.     {  
  65.         scanf("%lld%lld%lld%lld",&A,&B,&a,&b);  
  66.         LL L=B-A;  
  67.         LL ans=cal(a,b,L);  
  68.         if(ans==-1) printf("-1\n");  
  69.         else printf("%lld\n",ans);  
  70.     }  
  71.     return 0;  
  72. }  


POJ 1061 http://poj.org/problem?id=1061 青蛙的約會,裸的擴展歐幾里得

 

 

[cpp]  view plain  copy
 
  1. #include <iostream>  
  2. #include <cstdio>  
  3. #include <cstring>  
  4. #include <cmath>  
  5. #include <vector>  
  6. #include <string>  
  7. #include <queue>  
  8. #include <stack>  
  9. #include <algorithm>  
  10.   
  11. #define INF 0x7fffffff  
  12. #define EPS 1e-12  
  13. #define MOD 1000000007  
  14. #define PI 3.141592653579798  
  15. #define N 100000  
  16.   
  17. using namespace std;  
  18.   
  19. typedef long long LL;  
  20. typedef double DB;  
  21.   
  22. LL e_gcd(LL a,LL b,LL &x,LL &y)  
  23. {  
  24.     if(b==0)  
  25.     {  
  26.         x=1;  
  27.         y=0;  
  28.         return a;  
  29.     }  
  30.     LL ans=e_gcd(b,a%b,x,y);  
  31.     LL temp=x;  
  32.     x=y;  
  33.     y=temp-a/b*y;  
  34.     return ans;  
  35. }  
  36.   
  37. LL cal(LL a,LL b,LL c)  
  38. {  
  39.     LL x,y;  
  40.     LL gcd=e_gcd(a,b,x,y);  
  41.     if(c%gcd!=0) return -1;  
  42.     x*=c/gcd;  
  43.     b/=gcd;  
  44.     if(b<0) b=-b;  
  45.     LL ans=x%b;  
  46.     if(ans<=0) ans+=b;  
  47.     return ans;  
  48. }  
  49.   
  50. int main()  
  51. {  
  52.     LL x,y,m,n,L;  
  53.     while(scanf("%lld%lld%lld%lld%lld",&x,&y,&m,&n,&L)!=EOF)  
  54.     {  
  55.         LL ans=cal(m-n,L,y-x);  
  56.         if(ans==-1) printf("Impossible\n");  
  57.         else printf("%lld\n",ans);  
  58.     }  
  59.     return 0;  
  60. }  


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

 

 

[cpp]  view plain  copy
 
  1. #include <iostream>  
  2. #include <cstdio>  
  3. #include <cstring>  
  4. #include <cmath>  
  5. #include <vector>  
  6. #include <string>  
  7. #include <queue>  
  8. #include <stack>  
  9. #include <algorithm>  
  10.   
  11. #define INF 0x7fffffff  
  12. #define EPS 1e-12  
  13. #define MOD 1000000007  
  14. #define PI 3.141592653579798  
  15. #define N 100000  
  16.   
  17. using namespace std;  
  18.   
  19. typedef long long LL;  
  20. typedef double DB;  
  21.   
  22. LL e_gcd(LL a,LL b,LL &x,LL &y)  
  23. {  
  24.     if(b==0)  
  25.     {  
  26.         x=1;  
  27.         y=0;  
  28.         return a;  
  29.     }  
  30.     LL ans=e_gcd(b,a%b,x,y);  
  31.     LL temp=x;  
  32.     x=y;  
  33.     y=temp-a/b*y;  
  34.     return ans;  
  35. }  
  36.   
  37. LL cal(LL a,LL b,LL c)  
  38. {  
  39.     LL x,y;  
  40.     LL gcd=e_gcd(a,b,x,y);  
  41.     if(c%gcd!=0) return -1;  
  42.     x*=c/gcd;  
  43.     b/=gcd;  
  44.     if(b<0) b=-b;  
  45.     LL ans=x%b;  
  46.     if(ans<=0) ans+=b;  
  47.     return ans;  
  48. }  
  49.   
  50. int main()  
  51. {  
  52.     LL n,b,t;  
  53.     scanf("%I64d",&t);  
  54.     while(t--)  
  55.     {  
  56.         scanf("%I64d%I64d",&n,&b);  
  57.         LL ans=cal(b,9973,n);  
  58.         if(ans==-1) printf("Impossible\n");  
  59.         else printf("%lld\n",ans);  
  60.     }  
  61.     return 0;  
  62. }  


HDU 2669 http://acm.hdu.edu.cn/showproblem.php?pid=2669 裸的擴展歐幾里得

 

 

[cpp]  view plain  copy
 
  1. #include <iostream>  
  2. #include <cstdio>  
  3. #include <cstring>  
  4. #include <cmath>  
  5. #include <vector>  
  6. #include <string>  
  7. #include <queue>  
  8. #include <stack>  
  9. #include <algorithm>  
  10.   
  11. #define INF 0x7fffffff  
  12. #define EPS 1e-12  
  13. #define MOD 1000000007  
  14. #define PI 3.141592653579798  
  15. #define N 100000  
  16.   
  17. using namespace std;  
  18.   
  19. typedef long long LL;  
  20. typedef double DB;  
  21.   
  22. LL e_gcd(LL a,LL b,LL &x,LL &y)  
  23. {  
  24.     if(b==0)  
  25.     {  
  26.         x=1;  
  27.         y=0;  
  28.         return a;  
  29.     }  
  30.     LL ans=e_gcd(b,a%b,x,y);  
  31.     LL temp=x;  
  32.     x=y;  
  33.     y=temp-a/b*y;  
  34.     return ans;  
  35. }  
  36.   
  37. LL cal(LL a,LL b,LL c)  
  38. {  
  39.     LL x,y;  
  40.     LL gcd=e_gcd(a,b,x,y);  
  41.     if(c%gcd!=0) return -1;  
  42.     x*=c/gcd;  
  43.     b/=gcd;  
  44.     if(b<0) b=-b;  
  45.     LL ans=x%b;  
  46.     if(ans<=0) ans+=b;  
  47.     return ans;  
  48. }  
  49.   
  50. int main()  
  51. {  
  52.     LL a,b;  
  53.     while(scanf("%I64d%I64d",&a,&b)!=EOF)  
  54.     {  
  55.         LL ans=cal(a,b,1);  
  56.         if(ans==-1) printf("sorry\n");  
  57.         else printf("%I64d %I64d\n",ans,(1-ans*a)/b);  
  58.     }  
  59.     return 0;  
  60. }  


暫時就這么多了吧

 

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM