拓展欧几里德

写这个也是为了自己更好的学习数学专题

1.欧几里德算法

作用：求两个数a和b的最大公约数。

gcd(a,b)=gcd(b,a%b)

我们可以看到,上式是将原来函数的第二个参数放到了第一个参数的位置，然后第二个参数的位置放上 原来的第一个参数%原来第二参数。

高中时学过一个辗转相除法。

我们可以举个例子。

求24和10的最大公约数：

24  ÷ 10  =  2 ……4 

10  ÷  4   =  2 ……2      下面一个式子的被除数是上面一个式子的除数，除数是上一个式子的余数。

 4  ÷   2  =   2 ……0 

 2 ÷  0  (除数为0)→结果是 2

我们可以根据这个写一个程序：

#include <iostream>#include <cstdio>using namespace std;int gcd(int a,int b){    if(b==0)  return a;    return gcd(b,a%b);}int main(){    int a,b;    while(~scanf("%d%d",&a,&b))    {        printf("%d\n",gcd(a,b));    }    return 0;}

我们来测试测试

我给出的例子是被除数(a)比除数(b)大的情况，这个很好理解。（我并没给出严格的证明，我们只要能理解、能运用就好）

如果a<b呢？

举个例子：10和24

10 ÷ 24   =  0 ……10

24  ÷ 10  =  2 ……4 

.......

看到了如果 a<b那么gcd(a,b)的结果相当于 swap(a,b),gcd(a,b);

会把两个参数a，b交换位置，再调用一次。所以gcd(10,24)的结果gcd(24,10)相同，

（ps:这不是显然的吗?a、b的最大公约数是个定值)

我来考虑下这个函数的边界情况：

             if(b==0 && a),那么直接返回a了，结果=a，

              if(b &&  a==0),那么a<b,交换a，b的位置再调用一次,结果=b

              if(a==0&&b==0) 结果返回0

注意：实际上我们不会去求0和另一个数的最大公约数。

2.证明两整数a,b互质的充要条件是：存在两个整数x,y满足ax+by=1

充分性：设gcd(a,b)=h;

              所以  h|a ,  h|b;

               所以  h|a*x  ,h| b*x;

               所以   h| (a*x+b*y );

              又    a*x+b*y=1;

              所以h|1;

               所以h=1;

必要性：

        已知条件：(a,b)=1;

          设集合S={s|s=ax+by   },设d是集合s中最小的正整数,

        ①我们先证明一个事实：  集合S中任意一个元素都是d的倍数。（当然S中的元素可以是负数，但并不影响。）

        设t是S中一个任意元素,那么  t可以表示为这种形式   ：t=  p*d+r (0<=r<|d|，均为整数),

         因为t∈S，p*d∈S，故有r∈S（很显然，这是基于s=a*x+b*y这种形式）

         有因为0<=r<|d|，且d是集合s中最小的正整数

         故r=0;

        所以t∈S，故S中任意元素均为d的倍数。

        

       ②设 d=a*x0+b*y0;

              令t=a*(x0+1)+b*y0;

            又 d|t  ,d|a*x0+b*y0;

            所以 d|a;

             同理d|b;

           所以d|gcd(a,b);

            又gcd(a,b)=1;

           所以d|1;所以d=1;

          因为d是集合s中最小的正整数

          所以存在两个整数x,y满足ax+by=1

证毕。

通过这个结论我们很容易知道：

整数方程  a*X+b*Y=c的有解的充要条件是  gcd(a,b)|c   ;

(先对是a和b约分得到 a'X+b'Y=c/gcd(a,b)  ，然后就得到了附加条件gcd(a',b')=1，利用2， 得证）

证： 

设a' *gcd(a,b)=a,   b'* gcd(a,b) =b;

左边=gcd(a,b) (a'*X +b'*Y)=  gdc(a,b) * k  * (1)   (k属于正整数)

 右边=c,

若有解，必有   gcd(a,b)*k=c 

即 gcd(a,b)|c

3.拓展欧几里德算法

 

 我们回顾一下gcd(a,b)的求解方法，对于b!=0，结果=gcd(b,a%b)，它的意义就是a和b的最大公约数与b,a%b相等。

  对于递归调用的一系列函数中的a,b,他们的最大公约数都=gcd;

拓展欧几里德算法用于找到整数方程  a*X+b*Y=c ①的一组解。

首先我们已经知道，有解 等价于  gcd(a,b)|c  ,

那么我们可以先考虑求解 a*X+b*Y=gcd(a,b)            ②;

只需将②的解 乘以c/gcd(a,b) 即可得到原方程的一组解。（至少可以的到一组解，更多的解稍加处理就能得到）

现在考虑怎么求a*X+b*Y=gcd(a,b)            ②的解。

在求gcd(a,b)时，递归终点是b==0，这时我们可以写一个非常确定的等式：

递归终点方程组：

a*X+b*Y=gcd   ;

X=1；

这里X就=1，然而Y可以任意取值。

gcd就是要求的答案，在整个递归调用过程中都未曾发生改变。

我们有了一个起点，可以做的就是推理啊，

思想：

     gcd(a，b)调用的每一步都可以对应一个方程，其形式为a*X+b*Y=gcd   ，我们可以知道a，b在改变，gcd始终不变。

     现在我们知道了最后一步的X，Y值，我们要做的就是倒推，找出第一次调用函数时 0的X，Y值，那个就是原方程的一组解。

 

    设a*X+b*Y=gcd,递归调用得到子方程：b*X1+a%b*Y1=gcd;

    那么  b*X1+(a-a/b*b)*Y1=gcd;

           即 a*Y1+b*(X1-a/b*Y1)=gcd;

         对比a*X +b*Y                =gcd;

可得：

     X=Y1;

     Y=X1-a/b*Y1;

得解

程序:

#include <iostream>#include <cstdio>using namespace std;int exgcd(int a,int b,int &X,int &Y)//采用引用&，从递归终点不断修改X，Y值{    if(!b)    {        X=1,Y=0;//Y可取任意值，这里姑且取0，这样得到的结果abs  更小        return a;    }    int ans=exgcd(b,a%b,X,Y);    int X1=X,Y1=Y;    X=Y1;    Y=X1-a/b*Y1;    return ans;}int main(){    int a,b,X,Y;    while(~scanf("%d%d",&a,&b))    {        exgcd(a,b,X,Y);        printf("结果:\n%d  %d\n",X,Y);    }    return 0;}

说明：

 1.我们看到了代码中 终点的Y赋值为0，其实Y可以任意赋值，但是得到相应的解会不同，但都是原方程的解。

 2. exgcd的返回值是a,b的最大公约数，同时能得到a*X+b*Y=gcd(a,b) 的一组解X，Y。

X，Y不用赋初值。

3.解出的X、Y可能是负数，输入的a、b也可以是负数。

a*X+b*Y=gcd(a,b)的其它解：

我们已经得到了一组解x0,y0;

其它解可以表示为:

X=X0+   t* b/gcd(a,b);

Y=Y0- t*a/gcd(a,b); (t是整数)

为什么？：

我们这样写是为了得到连续的解，就是在某区内不漏掉一些解。

很好理解一点：X=X0+t*b;

                         Y=Y0-t*a  一定是原方程的解。

因为gcd(a,b)|b  并且gcd(a,b)|a 那么对于

X=X0+   t* b/gcd(a,b);

Y=Y0-   t*a/gcd(a,b); 

一定满足

a*X+b*Y=gcd(a,b)

这个时候X，Y的改变跨度一定比

X=X0+t*b;

Y=Y0-t*a 

小

而且关于X、Y的改变已经达到最小跨度，无法更小。

这是因为gcd（a,b）表示的是a和b的最大公约数，那么/gcd(a,b)就是相对应最小的，无法更小。

4.乘法逆元，拓展欧几里德求逆元

1.一般逆元

乘法逆元：如果a*x≡1(mod m) ，那么称x为a关于m的乘法逆元

给定a和m，拓展欧几里德可以求出a的关于m的乘法逆元。

解法：a*x=1(mod m) 等价于a*X+ m*Y=1; 

那么只用调用函数exgcd(a,m,X,Y)，然后X的值即是 a关于m的逆元。

如果看了前面的1、2、3，很容易看出逆元存在的条件：(a,m)=1  

2.最小正乘法逆元

我们通过1.得到的一般逆元 X0 进行调整 来得到最小正逆元，就是最小的正数X，满足方程a*X+m*Y=1;

首先我们由    3.拓展欧几里德算法

可以知道：

X=X0+ t*m/gcd;

Y=Y0- t*a/gcd;

现在我们只关心X， 因为gcd=1，所以

X=X0+t*m;

所以求最小的X，只需X%m即可，等等，还要考虑负数的情况，

经过我试验对于%运算,

a%b,不论分子分母的正负号情况,结果无外乎就是这两者之一的情况,(非负数)  m或者（负数） m-b;

所以只需要   b=abs(b); int  ans=a%b,if(ans<=0) ans+=b; 即为正解。

我看到网上的代码，保证了分母非负，这样也好，因为编译器不同，结果可能不一样。

代码如下：

#include <iostream>#include <cstdio>#include <cmath>using namespace std;int exgcd(int a,int b,int &X,int &Y)//采用引用&，从递归终点不断修改X，Y值{    if(!b)    {        X=1,Y=0;        return a;    }    int ans=exgcd(b,a%b,X,Y);    int X1=X,Y1=Y;    X=Y1;    Y=X1-a/b*Y1;    return ans;}int cal(int a,int m)//a*X=1(mod m);{    int X,Y;    int gcd=exgcd(a,m,X,Y);    if(1%gcd)  return -1;//此时方程无解    m=abs(m);//保证分母非负    int ans=X%m;   if(ans<=0)  ans+=m;   return ans;   }int main(){    int a,m;    while(~scanf("%d%d",&a,&m))    {        printf("结果:\n%d\n",cal(a,m));    }    return 0;}

3.a*X≡c(mod m)型

等价于求解   a*X+m*Y=c,利用exgcd(a,m,X,Y)求出一组解，然后

X=X0+ t*m/gcd;(注：m%gcd!=0时，无解)

所以最小的正数X=X0%(m/gcd),(这时没考虑分子分母为负数的情况)

代码中考虑了：

代码：

//a*X=c(mod m)int cal(int a,int m,int c){    int X,Y;    int gcd=exgcd(a,m,X,Y);    if(c%gcd) return -1;    X*=c/gcd;    m/=gcd;    m=abs(m);    return (X%m+m)%m;}

   

简单的题目:

1.ZOJ Problem Set - 3609

Modular Inverse

---

Time Limit: 2 Seconds      Memory Limit: 65536 KB

---

The modular modular multiplicative inverse of an integer a modulo m is an integer x such that a-1≡x (modm). This is equivalent toax≡1 (modm).

Input

There are multiple test cases. The first line of input is an integer T ≈ 2000 indicating the number of test cases.

Each test case contains two integers 0 < a ≤ 1000 and 0 < m ≤ 1000.

Output

For each test case, output the smallest positive x. If such x doesn't exist, output "Not Exist".

Sample Input

33 114 125 13

Sample Output

4Not Exist8

References

• http://en.wikipedia.org/wiki/Modular_inverse

---

#include<cstdio>#include<string>#include<cstring>#include<iostream>#include<cmath>#include<algorithm>#include<climits>#include<queue>#include<vector>#include<map>#include<sstream>#include<set>#include<stack>#include<cctype>#include<utility>#pragma comment(linker, "/STACK:102400000,102400000")#define PI 3.1415926535897932384626#define eps 1e-10#define sqr(x) ((x)*(x))#define FOR0(i,n)  for(int i=0 ;i<(n) ;i++)#define FOR1(i,n)  for(int i=1 ;i<=(n) ;i++)#define FORD(i,n)  for(int i=(n) ;i>=0 ;i--)#define  lson   num<<1,le,mid#define rson    num<<1|1,mid+1,ri#define MID   int mid=(le+ri)>>1#define zero(x)((x>0? x:-x)<1e-15)#define mk    make_pair#define _f     first#define _s     second/*const double  h=6.63*pow(10,-34);const double   c=3*pow(10,8);const double u=1.660552*1e-27;*/using namespace std;//const int INF=    ;typedef long long ll;//const ll inf =1000000000000000;//1e15;//ifstream fin("input.txt");//ofstream fout("output.txt");//fin.close();//fout.close();//freopen("a.in","r",stdin);//freopen("a.out","w",stdout);const int INF =0x3f3f3f3f;//const int maxn=    ;//const int maxm=    ;//by yskysker123ll exgcd(ll a,ll b,ll& X,ll &Y){    if(!b)    {        X=1,Y=0;        return a;    }    ll ans=exgcd(b,a%b,X,Y);    ll X1=X,Y1=Y;    X=Y1;    Y=X1-a/b*Y1;    return ans;}ll cal(ll a,ll m){    ll X,Y;    ll gcd=exgcd(a,m,X,Y);    if(1%gcd)  return 0;    m=abs(m);    ll ans=X%m;    if(ans<=0)  ans+=m;    return ans;}int main(){    ll a,m;    int T;scanf("%d",&T);    ll tmp;    while(T--)    {        scanf("%lld%lld",&a,&m);        if(!(tmp=cal(a,m) ) )  puts("Not Exist");        else printf("%lld\n",tmp);    }    return 0;}

2.ZOJ Problem Set - 3593

One Person Game

---

Time Limit: 2 Seconds      Memory Limit: 65536 KB

---

There is an interesting and simple one person game. Suppose there is a number axis under your feet. You are at pointA at first and your aim is pointB. There are 6 kinds of operations you can perform in one step. That is to go left or right bya,b andc, here c always equals to a+b.

You must arrive B as soon as possible. Please calculate the minimum number of steps.

Input

There are multiple test cases. The first line of input is an integer T(0 <T ≤ 1000) indicates the number of test cases. ThenT test cases follow. Each test case is represented by a line containing four integers 4 integersA,B,a and b, separated by spaces. (-2³¹ ≤A,B < 2³¹, 0 <a,b < 2³¹)

Output

For each test case, output the minimum number of steps. If it's impossible to reach pointB, output "-1" instead.

Sample Input

20 1 1 20 1 2 4

Sample Output

1-1

#include<cstdio>#include<string>#include<cstring>#include<iostream>#include<cmath>#include<algorithm>#include<climits>#include<queue>#include<vector>#include<map>#include<sstream>#include<set>#include<stack>#include<cctype>#include<utility>#pragma comment(linker, "/STACK:102400000,102400000")#define PI 3.1415926535897932384626#define eps 1e-10#define sqr(x) ((x)*(x))#define FOR0(i,n)  for(int i=0 ;i<(n) ;i++)#define FOR1(i,n)  for(int i=1 ;i<=(n) ;i++)#define FORD(i,n)  for(int i=(n) ;i>=0 ;i--)#define  lson   num<<1,le,mid#define rson    num<<1|1,mid+1,ri#define MID   int mid=(le+ri)>>1#define zero(x)((x>0? x:-x)<1e-15)#define mk    make_pair#define _f     first#define _s     second/*const double  h=6.63*pow(10,-34);const double   c=3*pow(10,8);const double u=1.660552*1e-27;*/using namespace std;//const int INF=    ;typedef long long ll;const ll inf =1000000000000000;//1e15;//ifstream fin("input.txt");//ofstream fout("output.txt");//fin.close();//fout.close();//freopen("a.in","r",stdin);//freopen("a.out","w",stdout);const int INF =0x3f3f3f3f;//const int maxn=    ;//const int maxm=    ;//by yskysker123ll exgcd(ll a,ll b,ll &X ,ll &Y ){    if(!b)    {       X=1,Y=0;       return a;    }    ll ans=exgcd(b,a%b,X,Y);    ll X1=X,Y1=Y;    X=Y1;    Y=X1-a/b*Y1;    return ans;}ll cal(ll a,ll b,ll c){    ll X,Y;    ll gcd=exgcd(a,b,X,Y);    if(c%gcd)  return -1;    X*=c/gcd;    Y*=c/gcd;    a/= gcd;//    a=abs(a);    b/= gcd;//    b=abs(b);    ll t=(Y-X)/(b+a);//可以认为y在以a的速度减小，X在以b的速度增加。注意a和b都是/l了gcd的，也可里认为是相遇问题。    ll ans=1e18;//之前ans初始化 值太小了，结果wa了    for(ll i=t-1;i<=t+1;i++)    {        ll tx=X+b*i;//可以这样理解tx和ty的正负代表了方向，长度为a的步子要么左走，要么右走，b同。        ll ty=Y-a*i;        ll tmp;        if(abs(ty)+abs(tx)==abs(tx+ty))//tx和ty同号（对于tx+ty=c情况下的判断方法），长度为a和b中的部分步子可以合并为1长度为c的步子            //且三者方向相同。          tmp=max( abs(ty),abs(tx) );        else  //tx和ty异号          tmp=abs(ty-tx)  ;//长度为a和b中的部分步子可以合并为1长度为c的步子，        //只不过a或b中1个和c 步子方向相同，另外一个与之方向相反。          ans=min(ans,tmp);    }    return ans;}int main(){  int T;  scanf("%d",&T);  ll A,B,a,b;  while(T--)  {      scanf("%lld%lld%lld%lld",&A,&B,&a,&b);      ll C=B-A;      ll tmp;      if((tmp=cal( a,b,C) )!=-1 ) printf("%lld\n",tmp);      else puts("-1");  }  return 0;}

3.hdu 2669  Romantic

Romantic

Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 3932    Accepted Submission(s): 1630

Problem Description

The Sky is Sprite.
The Birds is Fly in the Sky.
The Wind is Wonderful.
Blew Throw the Trees
Trees are Shaking, Leaves are Falling.
Lovers Walk passing, and so are You. 
................................Write in English class by yifenfei

 

Girls are clever and bright. In HDU every girl like math. Every girl like to solve math problem!
Now tell you two nonnegative integer a and b. Find the nonnegative integer X and integer Y to satisfy X*a + Y*b = 1. If no such answer print "sorry" instead.

 

Input

The input contains multiple test cases.
Each case two nonnegative integer a,b (0<a, b<=2^31)

 

Output

output nonnegative integer X and integer Y, if there are more answers than the X smaller one will be choosed. If no answer put "sorry" instead. 

 

Sample Input

77 5110 4434 79

 

Sample Output

2 -3sorry7 -3

 

Author

yifenfei

 

Source

HDU女生专场公开赛——谁说女子不如男

 

Recommend

lcy   |   We have carefully selected several similar problems for you:  2668 2674 2671 2670 2672 

 

#include<cstdio>#include<string>#include<cstring>#include<iostream>#include<cmath>#include<algorithm>#include<climits>#include<queue>#include<vector>#include<map>#include<sstream>#include<set>#include<stack>#include<cctype>#include<utility>#pragma comment(linker, "/STACK:102400000,102400000")#define PI 3.1415926535897932384626#define eps 1e-10#define sqr(x) ((x)*(x))#define FOR0(i,n)  for(int i=0 ;i<(n) ;i++)#define FOR1(i,n)  for(int i=1 ;i<=(n) ;i++)#define FORD(i,n)  for(int i=(n) ;i>=0 ;i--)#define  lson   num<<1,le,mid#define rson    num<<1|1,mid+1,ri#define MID   int mid=(le+ri)>>1#define zero(x)((x>0? x:-x)<1e-15)#define mk    make_pair#define _f     first#define _s     second/*const double  h=6.63*pow(10,-34);const double   c=3*pow(10,8);const double u=1.660552*1e-27;*/using namespace std;//const int INF=    ;typedef long long ll;//const ll inf =1000000000000000;//1e15;//ifstream fin("input.txt");//ofstream fout("output.txt");//fin.close();//fout.close();//freopen("a.in","r",stdin);//freopen("a.out","w",stdout);const int INF =0x3f3f3f3f;//const int maxn=    ;//const int maxm=    ;//by yskysker123ll exgcd(ll a,ll b,ll &X,ll &Y){    if(!b)    {        X=1,Y=0;        return a;    }    ll gcd=exgcd( b,a%b,X,Y );    ll X1=X,Y1=Y;    X=Y1;    Y=X1-a/b*Y1;    return gcd;}bool work(ll a,ll b,ll &X,ll &Y){   ll gcd= exgcd(a,b,X,Y);   if( 1%gcd ) return false;   X=(X%b+b)%b;   Y=  (1-X*a)/b;   return true;}int main(){  ll a,b, X,Y;  while(~scanf("%lld%lld",&a,&b))  {      if(!work(a,b,X,Y)) puts("sorry");      else printf("%lld %lld\n",X,Y);  }    return 0;}

4.  hdu 1576   A/B

A/B

Time Limit: 1000/1000 MS (Java/Others)    Memory Limit: 32768/32768 K (Java/Others)
Total Submission(s): 3306    Accepted Submission(s): 2507

Problem Description

要求(A/B)%9973，但由于A很大，我们只给出n(n=A%9973)(我们给定的A必能被B整除，且gcd(B,9973) = 1)。

 

Input

数据的第一行是一个T，表示有T组数据。
每组数据有两个数n(0 <= n < 9973)和B(1 <= B <= 10^9)。

 

Output

对应每组数据输出(A/B)%9973。

 

Sample Input

21000 5387 123456789

 

Sample Output

79226060

 

Author

xhd

 

Source

HDU 2007-1 Programming Contest

 

Recommend

linle

 

若A*X=1(mod m)，则称X为A关于m的逆元，A为X关于m的逆元。

那么C*A%m=C/X%m ;

       C/A%m= C*X%m ;

#include<cstdio>#include<string>#include<cstring>#include<iostream>#include<cmath>#include<algorithm>#include<climits>#include<queue>#include<vector>#include<map>#include<sstream>#include<set>#include<stack>#include<cctype>#include<utility>#pragma comment(linker, "/STACK:102400000,102400000")#define PI 3.1415926535897932384626#define eps 1e-10#define sqr(x) ((x)*(x))#define FOR0(i,n)  for(int i=0 ;i<(n) ;i++)#define FOR1(i,n)  for(int i=1 ;i<=(n) ;i++)#define FORD(i,n)  for(int i=(n) ;i>=0 ;i--)#define  lson   num<<1,le,mid#define rson    num<<1|1,mid+1,ri#define MID   int mid=(le+ri)>>1#define zero(x)((x>0? x:-x)<1e-15)#define mk    make_pair#define _f     first#define _s     second/*const double  h=6.63*pow(10,-34);const double   c=3*pow(10,8);const double u=1.660552*1e-27;*/using namespace std;//const int INF=    ;typedef long long ll;//const ll inf =1000000000000000;//1e15;//ifstream fin("input.txt");//ofstream fout("output.txt");//fin.close();//fout.close();//freopen("a.in","r",stdin);//freopen("a.out","w",stdout);const int INF =0x3f3f3f3f;//const int maxn=    ;//const int maxm=    ;//by yskysker123const ll P=9973;ll exgcd(ll a,ll b,ll &X,ll &Y){    if(!b)    {        X=1,Y=0;        return a;    }    ll ans=exgcd(b,a%b,X,Y);    ll X1=X,Y1=Y;    X=Y1;    Y=X1-a/b*Y1;    return ans;}ll cal(ll a,ll m){    ll X,Y;    ll gcd=exgcd(a,m,X,Y);    m=abs(m);    X=(X%m+m)%m;    return X;}int main(){    int T;    scanf("%d",&T);//已知A%9973,要求A/B%9973,等价于A*(B的逆元) %9973,只需求 B的逆元；    ll A,B;    while(T--)    {        scanf("%lld%lld",&A,&B);        ll t=cal( B,P);         printf("%lld\n",A*t%P  );    }    return 0;}