关于欧几里得及其扩展算法（C语言实现）

    作为新人Acmer，这两天刚刚学习了欧几里得（扩展算法），为方便以后复习，特地记录一下此算法，作为个人笔记。因水平有限，如有纰漏，日后再完善！
    1.首先我们知道，欧几里得算法是求两个正整数a,b的最大公因数gcd(a,b)，这里不妨设(a>b>0).

先附上代码：int gcd(int a,int b)
{
if(b==0) return a;
else return gcd(b,a%b);
}这是什么意思呢？
    这是一个递归代码，出口条件是b==0,否则则让a=b;  b=a%b;一直递归下去直到满足出口条件.
2.下面给出证明：
设c==gcd(a,b);  r=a%b;
则有 a==x1*c;  b==x2*c;  r==a-x3*b;  (x1,x2,x3均为正整数)
将a,b代入得：r==x1*c-x3*x2*c==(x1-x3*x2)*c;
因此a,b的余数r是gcd(a,b)的倍数；
因此可以得到：c==gcd(b,r);  从而得到gcd(a,b)==gcd(b,r);          
所以当循环一直进行下去直到r==0,所求的数即gcd(a,b);
证毕.
3.欧几里得算法的优点：
    显然欧几里得算法通过降维的方式使得计算最大公因数的时间复杂度大大降低，尤其是针对较大的数字运算时。因为欧几里得算法的特性，人们又给它起了个好听的名字：辗转相除法。

-------------------------------------------------------华丽的分割线-------------------------------------------------------

4.接下来要介绍的是基于欧几里得算法的扩展算法(Euclid v2.0)
      我们知道，对于自然界的正整数...a,b...一定满足ax+by==gcd(a,b);  (解一定存在，根据数论中的相关定理）
下面我们来证明ta:
(1)...ax+by==gcd(a,b);
(2)...bx1+a%by1==gcd(b,a%b);  (运用欧几里算法)
(3)...gcd(a,b)==gcd(b,a%b);  (欧几里得算法)
(4)...ax+by==bx1+a%b*y1;  (在计算机里a%b==(a-a/b*b))
(5)...ax+by==bx1+ay1-a/b*by1;
(6)...ax+by==ay1+b(x1-a/b)y1;  (合并同类项)
(7)...x=y1,y=x1-a/b*y1;  (结论)
    所以我们知道x,y和方程的下一个状态量(x1,y1)有关，而方程最后一个状态（即gcd(a,b)终点，b==0; 还记得吗？）此时：
ax+0==a;  (gcd(a,b)的最终返回值为a)  可以得到一组解{x==1; y==0}，根据此解层层上推即可得到最开始的x,y.
5.下面我们通过代码来实现ta：
我们知道，C语言只有取地址符没有引用，所以我们用指针来实现扩展欧几里得算法（C++用引用实现）：int exgcd(int a,int b,int *x,int *y)
{
if(b==0)
{
*x=1;
*y=0;
return a;
}
int ret=exgcd(b,a%b,x,y);
int t=*x;
*x=*y;
*y=t-a/b*(*y);
return ret;
}注意参数传递方式：（地址传递）

exgcd(a,b,&x,&y);

什么？不怎么会指针？没关系，我们还可以用结构实现ta：struct Euclid
{
int x,y,gcd;
}AC;//全局变量;

int exgcd(int a,int b)
{
if(b==0)
{
AC.x=1;
AC.y=0;
return a;
}
int ret=exgcd(b,a%b);
int t=AC.x;
AC.x=AC.y;
AC.y=t-a/b*AC.y;
return ret;
}引用该模块也很简单：AC.gcd=exgcd(a,b);6.说再多不如水一题    ╮(￣▽￣)╭    ：NYOJ 775 --> 整数性质



这题不多说了，上代码：
#include<stdio.h>
int main()
{
int a,b,s,t,gcd;
int exgcd(int a,int b,int *x,int *y);
while(scanf("%d%d",&a,&b)!=EOF)
{
gcd=exgcd(a,b,&s,&t);
printf("%d %d\n",s,t);
}
return 0;
}
int exgcd(int a,int b,int *x,int *y)//扩展欧几里得算法;
{
if(b==0)
{
*x=1;
*y=0;
return a;
}
int ret=exgcd(b,a%b,x,y);
int t=*x;
*x=*y;
*y=t-a/b*(*y);
return ret;
}


-------------------------------------------------------华丽的分割线-------------------------------------------------------

7.难道欧几里得扩展算法就这样完了吗?--当然不是    (⌒▽⌒)    
我目前遇到的情况就是以下的两个方面：
(1).解不定方程；(2).求乘法逆元；
7.1.下面我们来研究研究不定方程：（￣へ￣）
设方程ax+by==c，对比ax+by==gcd(a,b)，可知当方程有整数解的充要条件是:  (c%gcd(a,b)==0) 即c是gcd(a,b)的整数倍；
为什么呢，因为ax+by==gcd(a,b)一定有整数解，当gcd(a,b)扩大整数倍k时，只需将x,y同时扩大k倍，也就得解了。
上代码：int solve(int a,int b,int *x,int *y)
{
int gcd=exgcd(a,b,x,y);//引用扩展欧几里得算法;
if(c%gcd==0)
{
int k=c/gcd;
*x*=k;//所求的是一组特解;
*y*=k;
return 1;
}
else return 0;//无整数解;
}
    为什么是特解呢？(ﾟДﾟ≡ﾟдﾟ)!? 因为这是不定方程啊！理论上有无数组解，其实这个历史遗留问题还是得回到第4节：我们当时求出的其实是x0，y0（也是一组特解），方程ax+by==gcd(a,b)所有整数解满足：x==x0+b/gcd(a,b)*t;

y==y0-a/gcd(a,b)*t; （t为任意整数）    又是为什么呢？我国有个成语：此消彼长，你喜欢说拆东墙补西墙我也ok  (〜￣△￣)〜 
总之将这些解代入方程，你就会发现加的部分和减的部分刚好抵消了  (°∀°)ﾉ ，得到的还是ax0+by0==gcd(a,b);  显然成立。
那么关于ax+by==c的所有整数解，自然就是：
x==(x0+b/gcd(a,b)*t)*k;

y==(y0-a/gcd(a,b)*t)*k; （t为任意整数,k=c/gcd(a,b））这里就不上例题，大家自行琢磨哈
7.2.关于乘法逆元：
    首先什么叫逆元？广义上说，逆元就是能够抵消一个运算数某种运算的另一个运算数，在数学上，除一个数等于乘以其倒数，倒数被认为是乘法逆元（这只是狭义的），总之概念就是理解一下就行，我们来求同余方程的乘法逆元：
对于一个同余方程ax≡b (mod n)，当且仅当ax≡1 (mod n)时，我们称x为a关于模n的乘法逆元；不难看出，gcd(a,n)==1时a才存在乘法逆元，因为当a,n不互质时，b==0；与  b==1矛盾，所以gcd(a,n)==1为同余方程乘法逆元存在的充要条件。所以我们将同余方程化为ax+by==1;即ax+by==gcd(a,b);
    是不是很熟悉？接下来就是用扩展欧几里得算法求解了，只有一点要注意：逆元里面的求余是只能求余成正数的，范围是（0，n-1），所以一般对于计算机，数论题的x%n我们的写法都是( x % n + n ) % n这样保证能够产生0到n-1的正数，接下来上一个水题： NYOJ 769--> 乘数密码



这里稍微介绍一下乘（法）数密码，可以参考百度百科：乘法密码
注意这题字母表顺序a从0开始而不是1！所以你拿百度的例子运行解密不过！因为百度a从1开始（/TДT)/，关键是掌握解密原理就好了。。。
上代码：
#include<stdio.h>
#include<string.h>
int main()
{
char code[51];
int exgcd(int a,int b,int *x,int *y);
int k,i,_k,y,len,gcd;
while(scanf("%s%d",code,&k)!=EOF)
{
int q=26;
len=strlen(code);
gcd=exgcd(k,q,&_k,&y);
_k=(_k%q+q)%q;//求解逆元_k;
for(i=0;i<len;i++)
code[i]=_k*(code[i]-'A')%q+'A';//逆天了！字母表顺序从0开始啊！
for(i=0;i<len;i++)
printf("%c",code[i]);
printf("\n");
}
return 0;
}
int exgcd(int a,int b,int *x,int *y)//扩展欧几里得算法;
{
if(b==0)
{
*x=1;
*y=0;
return a;
}
int ret=exgcd(b,a%b,x,y);
int t=*x;
*x=*y;
*y=t-a/b*(*y);
return ret;
}
这里再送一个求逆元模板：int reverse(int a,int b,int *x,int *y)
{
int gcd=exgcd(a,b,x,y);
if(gcd==1)
{
*x=(*x%b+b)%b;//保证x为正数;
return 1;
}
else return 0;//逆元不存在;
}

-------------------------------------------------------华丽的分割线-------------------------------------------------------

关于扩展欧几里得算法另一个问题（模线性方程）正在研究中，下次来补上 (⌒▽⌒)