欧几里得算法--二进制欧几里得算法--拓展欧几里得算法

首先说一下，欧几里得算法：求两个数的最大公约数

教材的写法是：

//输入要求a>=b>=0 
unsigned euclid(unsigned a,unsigned b){
    if(b==0) return a;  
    return euclid(b,a%b);
}

当然也可以改成非递归的迭代版本：

int Gcd(int a, int b)
{
    while(b != 0)
    {
    　　int r = b;
    　　b = a % b;
    　　a = r;
    }
    return a;
}

上个世纪60年代研究人员对这个古老的算法进行了改进，称为二进制欧几里得算法：

基于的原理是：

1. 若 N 和 M 是偶数，gcd(N, M) = 2 gcd(N/2, M/2) ，
2. 若 N 是偶数而 M 是奇数，那么 gcd(N, M) = gcd(N/2, M) ，
3. 若 N 和 M 都是奇数，那么 (由于 N-M 是偶数) |N-M| < max(N,M) 。用 |N-M| 替换两者中的较大者   
   
   该算法称为 二进制 因为，与原来的不同，它不使用通常的除法而只使用除以2的。由于计算机数字总是用二进制系统表示，你不会对有一种特殊的 机器指令 以高效的方式实现 除以 2而感到惊讶。这称为 右移位，其中最右一位被舍弃，其他的位向右移动一位，而最左一位被设为0。
   

   另一个顺手拈来的操作是按位与 & 。对任意整数 N 来说，N & 1 是 1 或 0 。当且仅当 N 是奇数时它为 1 。按位与也在硬件中实现，例如，作为 机器指令 

递归版本：

unsigned int gcd(unsigned int u, unsigned int v)
{
    // simple cases (termination)
    if (u == v) return u;
    if (u == 0) return v;
    if (v == 0) return u;

    // look for factors of 2
    if (~u & 1) // u is even
    {
        if (v & 1)  return gcd(u >> 1, v);             // u is even,v is odd
        else        return gcd(u >> 1, v >> 1) << 1;   // both u and v are even
            
    }
    if (~v & 1)    return gcd(u, v >> 1);              // u is odd, v is even
    
    
    // both u and v are odd,reduce larger argument
    if (u > v)     return gcd((u - v) >> 1, v);

    return gcd((v - u) >> 1, u);
}

当然也可以改成迭代版本：

unsigned int gcd2(unsigned int u, unsigned int v)
{
  int shift;

  /* GCD(0,v) == v; GCD(u,0) == u, GCD(0,0) == 0 */
  if (u == 0) return v;
  if (v == 0) return u;
 
  /* Let shift := lg K, where K is the greatest power of 2
        dividing both u and v. */
  for (shift = 0; ((u | v) & 1) == 0; ++shift) {
         u >>= 1;
         v >>= 1;
  }
 
  while ((u & 1) == 0)
    u >>= 1;
 
  /* From here on, u is always odd. */
  do {
       /* remove all factors of 2 in v -- they are not common */
       /*   note: v is not zero, so while will terminate */
       while ((v & 1) == 0)  /* Loop X */
           v >>= 1;

       /* Now u and v are both odd. Swap if necessary so u <= v,
          then set v = v - u (which is even). For bignums, the
          swapping is just pointer movement, and the subtraction
          can be done in-place. */
       if (u > v) {
         unsigned int t = v; v = u; u = t;}  // Swap u and v.
       v = v - u;                       // Here v >= u.
     } while (v != 0);

  /* restore common factors of 2 */
  return u << shift;
}

拓展欧几里得的算法是：

        欧几里得定理告诉我们，只要能找到整数x和y，使得ax+by=d,则d就是a和b的最大公因子。而且如果d是a和b的最大公因子，则d一定可以表示为ax+by的形式。    

        运用下列函数的前提是gcd(a,b)=d,即d是a和b的最大公约数，则该函数可以计算得到ax+by=d中的x值。

//要求是gcd(a,b)=d      
int e_gcd(int a,int b,int &x,int &y){
    if(b==0){
        x=1;
        y=0;
        return a;
    }
    int ans=e_gcd(b,a%b,x,y);
    int temp=x;
    x=y;
    y=temp-a/b*y;
    return ans;
} 

拓展欧几里得算法的一个运用是计算模n逆元：如果有ax=1(mod N),那么我们称x是a的模N逆元。

//计算得到 “a模m逆元” 
int cal(int a,int m){   
    int x,y;
    int gcd=e_gcd(a,m,x,y);
    if(1%gcd!=0)  return -1;
    x*=1/gcd;
    m=abs(m);
    int ans=x%m;
    if(ans<=0) ans+=m;
    return ans;
}

其实，“a的模N逆元”，主要在RSA中计算解密密钥时用的较多。关于RSA的细节，这里不赘述，详情参加另外一篇博客。在计算RSA中加密时，需要用到快速幂取模算法，详见九月份的一篇博客。

贴一个应用：

ural  Idempotents

Time Limit:1000MS     Memory Limit:65536KB     64bit IO Format:%I64d & %I64u

Description

The number x is called an idempotent modulo n if

x* x = x (mod n)

Write the program to find all idempotents modulo n, where n is a product of two distinct primes p and q.

Input

First line contains the number k of test cases to consider (1 ≤ k ≤ 1000). Each of the following k lines contains one number n < 10 ⁹.

Output

Write on the i-th line all idempotents of i-th test case in increasing order. Only nonnegative solutions bounded by n should be printed.

Sample Input

input	output
3 6 15 910186311	0 1 3 4 0 1 6 10 0 1 303395437 606790875

AC代码：

#include <iostream>
using namespace std;


int getp(int n){   //得到大于1的最小的因数 
    int p=2;
    while(n%p!=0){
        p++;
    }
    return p;
}

void e_gcd(int a,int b,int &x,int &y){
    if(b==0){
        x=1;
        y=0;
        return ;
    }
    e_gcd(b,a%b,x,y);
    int temp=x;
    x=y;
    y=temp-a/b*y;
    return ;
} 

int main()
{
    int k;
    cin>>k;
    while(k--)
    {
        int n;
        cin>>n;
        int p=getp(n),q=n/p;
        int x,y;
        e_gcd(p,q,x,y) ;
        int x1=x*p;
        if(x1<0) x1+=n;
        
        e_gcd(q,p,x,y) ; 
        int x2=x*q;
        if(x2<0) x2+=n;
        if(x2<x1){
            int temp=x2;
            x2=x1;
            x1=temp;
        }
        cout<<0<<" "<<1<<" "<<x1<<" "<<x2<<endl;
    } 
    
}