扩展欧几里得算法

        给予二整数 a 与 b, 必存在有整数 x 与 y 使得  ax + by = gcd(a,b)  。

       有两个数a,b，对它们进行辗转相除法，可得它们的最大公约数——这是众所周知的。然后，收集辗转相除法中产生的式子，倒回去，可以得到ax+by=gcd(a,b)的整数解。

1、算法原理

        给定整数a、b，若gcd(a,b)=1。则存在c，满足ac=1 mod b，c即为a模b的乘法逆元。

       利用扩展的欧几里得算法求得满足条件的c：先做辗转相除，当a、b互素时，最后一步得到的余数为1，再从1出发，对前面得到的所有除法算式进行变形，将余数用除数和被除数表示，最终便可将1表示为a与b的一种线性组合，即

       从而x就是a模b的乘法逆元。因此寻找乘法逆元的过程就是求x和y的过程。 

       这里我们先看一下人工计算的过程：

       具体实现时使用三组变量：x1，x2，x3；y1，y2，y3；t1，t2，t3。初始化时，给x1，x2，x3别赋值1、0、a，则有

1 × a + 0 × b = a

       类似地，给y1，y2，y3分别赋值0、1、b，有

0 × a + 1 × b = b

       接下来开始迭代，保证在每次迭代中，有

a × t1 + b × t2 = t3

成立。为此只需在每一步中为 ti 分别赋值  (i=1,2,3) ，其中q为x3除以y3的商。

       当最后t3迭代至计算结果为1时，相应的 t1 就是 a模b的乘法速元，而 t2 是b模a的乘法逆元。

2、代码部分

#include <stdio.h>
int main(int argc, char *argv[])
{
	int temp,q,t1,t2,t3;
	int a,b,swap=0;
	int x1,x2,x3,y1,y2,y3;
	printf("Input two integers: ");
	scanf("%d",&a);
	scanf("%d",&b);
	if (a<b) {   //保证a>b
		swap = 1;
		temp = a;
		a = b;
		b = temp;		
	}
	x1=1;	x2=0;	x3=a;    // 初始
	y1=0;	y2=1;	y3=b;
	while (y3!=0)
	{
		q=x3/y3; //商
		t1=x1-q*y1;	t2=x2-q*y2;	t3=x3-q*y3;
		x1=y1;		x2=y2;		x3=y3;
		y1=t1;		y2=t2;		y3=t3;
		printf("\nt1,t2,t3\t%d\t%d\t%d ;",t1,t2,t3);	
	}
	printf("\n");
	if(x3 == 1)   // 这里使用x3，组后一轮循环中，x3保存了上一步中值为1的t3
	{
		if(swap==1)
		{
			printf("\ninverse of %d mod %d is:%d",b,a,x2);
			printf("\ninverse of %d mod %d is:%d",a,b,x1);		
		}
		else{
			printf("\ninverse of %d mod %d is:%d",a,b,x2);
			printf("\ninverse of %d mod %d is:%d",b,a,x1);	
		}			
	}
	else  printf("no inverse");
	printf("\n\n\n");
	return 0;
}

注释：

       关于判断条件是  x3==1，而不是 t3==1 的问题，while循环中本质使用的是辗转相除法，对于 ti 的赋值表达式，先看t3，t3=x3-q*y3 ，q=x3/y3 ，我理解的是 t3 就是 x3除以y3的余数，加上开始时,给 x3=a，y3=b ，所以对t3的迭代可以看成对 （a,b）做辗转相除法。

      再看 t1 和 t2，为了保持线性关系ax + by = gcd(a,b)，相应的x和y的值也需要变化，这就是 t1 和 t2 的作用。

其实，一直就想写这篇博客的，因为上次在写RSA程序，我看了网上的资料都不够简洁，有的提到了方法，但是代码却不够简便，例如欧几里得扩展算法的矩阵形式，以及商的倒序求逆，方法众多，但是代码实现却感觉有点复杂。况且我也看了前人的博客，需要开辟空间对商进行存储。这次我就把我目前感觉最简单的求模逆方法分享给大家~有错误请大家多多指出！

问题：

求e关于模p的逆元d，即要求出整数d，使e * d mod p = 1 （或 ed+px=1），这里要求e与p互素。

方法：

这里是引用初等数论（严士健 第三版）第27页的辗转相除法的表格形式的辗转相除法。

首先我们先要引入概念：
对于二元一次不定方程的求解：
即 ax+by=1的求解方法，

ax+by=1,(a,b)=1
a=bq₁ + r ₁ , 0<r₁ <b,
b=r₁q₂ + r ₂ , 0<r₂ <r₁ ,
……
r_n-2 =r_n-1 q_n + r_n , 0<r_n <r_n-1 ,
r_n-1 =r_n q_n+1 .

因为（a,b）=1,故r_n =1
于是有 a[(-1)^n-1Q_n ]+b[(-1)ⁿP_n ]=1
P₀ =1 , P₁ =q₁ ,P_k =q_k P_k-1 +P_k-2 ,
Q₀ =0 , Q₁ =1,Q_k =q_k Q_k-1 +Q_k-2,
k=2,…,n.

由于我时间不是很多，我就不举具体例子了 ，防止大家视疲劳我还是手写吧。

直接上代码 ，代码思想就是我所列举的辗转相除法思想，具体证明过程有兴趣的朋友可以参考：《初等数论》(第三版）作者：闵嗣鹤 严士健
第一章ξ3定理1的证明过程 以及 第二章ξ1定理2的证明过程
（不过看到好像这么多观看量，我还是列出来吧~~）

数学证明：

定理1：若a,b是任意两个正整数，则
Q_k a - P_k b = (-1)^k-1 r_k ,k=1,…，n;
其中
P₀ =1 , P₁ =q₁ ,P_k =q_k P_k-1 +P_k-2 ,
Q₀ =0 , Q₁ =1,Q_k =q_k Q_k-1 +Q_k-2 ,
k=2,…,n.

证：

当k=1时，上式显然成立，当k=2时，
r₂ = -[aq₂ - b(1 +q₁ q₂ )].
但
1 +q₁ q₂ =q₂ P₁ +P₀ ,
q₂ =q₂ *1 +0 = q₂ Q₁ +Q₀ ,
故
Q₂ a - P₂ b = (-1)^2-1 r₂ ,
P₂ =q₂ P₁ +P₀ ,
Q₂ =q₂ Q₁ +Q₀ …
假定上述式子对于不超过k>=2的正整数都成立，则
(-1)^k r_k+1
= (-1)^k (r_k-1 - q_k+1 r_k )
= (Q_k-1 a - P_k-1 b) + q_k+1 (Q_k a - P_k b)
= (q_k+1 Q_k + Q_k-1 )a - (q_k+1 P_k + P_k-1 )b.
故
Q_k a - P_k b = (-1)^k-1 r_k ,
其中
P_k+1 =q_k+1 P_k +P_k-1 , Q_k+1 =q_k+1 Q_k +Q_k-1 ,
由归纳法，定理为真

代码：

我们先判断素数，素数的判断方法我所知道的有三种，纳筛法，miller-素性检验，暴力求解，miller-素性检验和纳筛法由于我的实验报告被我清理文件的时候删了，我就用我同学的暴力求解法（一般小素数我们为了节约时间优先考虑纳筛法，但是对于大素数推荐采用miller-素性检验，但是miller素性检验是概率算法，1/4可能性会产生伪素数，需要多次验证，验证k次，产生伪素数的可能性就减少为(1/4)^k ）

//判断n是否为素数，此处使用最暴力的方法，将n取余小于n的所有值
void searchprime(int n)
{
	   int i;
	   for(i=2;i<=n-1;i++)
	   {
		   if(n%i==0)
			   break;
	   }//将小于n的数依次除n，若可以整除，则不是素数，退出循环
       if(i<n-1)
		  printf("%d不是素数\n",n);
	   else
		  printf("%d是素数\n",n);
}

然后我们进行求模逆操作

//模逆运算，用模逆算法，此处编写程序时候u,v为实际计算过程中v的各个值
int mod_invese(int d,int n)//求d模n的逆，n 为正整数
{
    int a,b,q,r,u=0,v=1,t;//a为被除数，b为除数，q为商
    a=n;
    b=(d>=0)?(d%n):-(d%n);//规定正负数mod（n）都取正；
    while(b!=0)
	{
        q=(int)a/b;
        r=a-b*q;
        a=b;
        b=r;
        t=v;
        v=u-q*v;
        u=t;//把每次没有进行运算的v的值给u
	}//辗转相除，直到余数b=0，跳出循环
    if(a!=1) return 0;
    return((u<0)?u+n:u);//返回的值都是mod(n)后的值，使得最后的值是正数；
}