经典算法（5）- 用二进制方法实现扩展的最大公约数(Extended GCD)

二进制方法中，只需要移位（<<和>>）和加减操作（+和-），不像欧几里德算法中需要乘法和除法运算。虽然算法效率更高，但是程序的可读性和可维护性差一些。

如果设d=gcd(u,v) = u.x + v.y, 本算法涉及到六种操作：

1）已知ext_gcd(u,v)如何求ext_gcd(u,2v)=u'.x' + v'.y'，其中u为奇数，v可奇可偶，d=gcd(u,v)为奇数；

2）已知ext_gcd(u,v)如何求ext_gcd(2u,v)=u'.x' + v'.y'，其中v为奇数，u可奇可偶，d=gcd(u,v)为奇数；

3）已知ext_gcd(u-v,v)如何求ext_gcd(u,v)=u'.x' + v'.y'；

4）已知ext_gcd(u,u-v)如何求ext_gcd(u,v)=u'.x' + v'.y'；

5）已知u/2^c = v = d（c为大于0的整数），如何求ext_gcd(u,v)=u.x' + v.y'；

6）已知u= v/2^c = d（c为大于0的整数），如何求ext_gcd(u,v)=u.x' + v.y'.

其中第1）种操作比较麻烦，第2）种操作只需要在第1）种操作的基础上将u和v交换一下。下面介绍一下第1）种操作的原理：

在第1）操作中，因为d=gcd(u,2v)=gcd(u-v,v)=(u-v)x+v(x+y)，设u1=u-v，x1=x，v1=v，y1=x+y，即问题转化为已知d=u1.x1+u1.y1，求ext_gcd(u,2v)=u.x' + (2v)y'=d；

如果y为偶数，显然可以采用这样一组值：x'=x，y'=y/2；但是如果y为奇数，则需要将d表示为d=u1(x1+v1.k) + v1(y1-u1.k)，其中k为正奇数1，3，5....，找到这样的k，满足x' = x +v.k和y'=(y-u.k)/2都不为0。

第3）种操作比较简单，因为d=gcd(u-v,v) = (u-v)x + vt = ux + v(y-x)，即x'=x，y'=y-x;

第4）种操作类似，因为d=gcd(u,u-v) = ux + (u-v)t = u(x+y) - vy，即x'=x+y，y'=-y;

第5）种操作,可以得到d=v=u + (1-2^c).v，即x'=1，y'=1-2^c；

第6）种操作只需要在第5）种操作的基础上将u和v交换一下。

在下面的实现中，只有上面第2）和6）种操作需要改变v和y的符号。

import java.util.Stack;
/**
 * 
 * @author ljs 2011-5-19
 * 
 * solve extended gcd using Binary Method
 *
 */
public class ExtGCD_Binary {
	/**
	 * caculate x,y to satisfy Bezout Identity: u.x + v.y = gcd(u,v)
	 * @param u
	 * @param v
	 * @return {d,x,y} for: d = u.x + v.y;
	 */
	public static int[] extGCDBinary(int u,int v){
		u=(u<0)?-u:u;
		v=(v<0)?-v:v;
		
		if(u==0)
			return new int[]{v,0,1};
		if(v==0)
			return new int[]{u,1,0};
		
		int k=0;
		while((u & 0x01)==0 && (v & 0x01) == 0){
			u>>=1; //divide by 2
			v>>=1;
			k++;
		}
		//with respect to d=gcd(u,v)=u.s+v.t
		/* the operation type is defined as:
		  1: c, u / 2^c (push two numbers to stack: first push c, then push the op number 1)
		  2: c, v / 2^c  (push two numbers to stack: first push c, then push the op number 2)
		  3: gcd(u-v, v) (only push the op number 3 to stack)
		  4: gcd(u, u-v) (only push the op number 4 to stack)  
		 */
		Stack<Integer> ops= new Stack<Integer>();
		//at this time, there is at least one number is odd between m and n
		int t=-v; //set it negative for later comparison of (t>0)
		if((v & 0x01)==1){
			//if v is odd
			t = u;
		}
		
		//process t as a possible even number
		boolean firstNoOp = true;
		while(t != 0){
			if(firstNoOp) {
				firstNoOp = false;
			}else{
				if(t>0)
					ops.push(3);
				else
					ops.push(4);
			}
			
			int c=0;
			while((t & 0x01)==0){
				//do until t is not even 
				t>>=1;
				c++;
			}
			if(t>0) {//u > v (the max is replaced by |t|)
				u=t; 
				if(c>0) {
					ops.push(c);  
					ops.push(1);		
				}
			}else{ //u<v (the max is replaced by |t|)
				v=-t;
				if(c>0){
					ops.push(c);   
					ops.push(2);
				}
			}
			//now u and v are all odd, then u-v is even
			t = u-v;	
		}
		int d = u;
				
		//the following steps are to caculate x,y in d=u.x + v.y
		
		int x=1,y=0;
		//special processing for the last operation
		if(!ops.isEmpty()){
			//e.g. for input u=3,v=3, then ops is empty, d=3
			int op = ops.pop();
			int c = ops.pop();
			//op number 3 and 4 can not be the last step since: u-v = v, then v is even; or u=u-v, then u=0; 
			//both are not possible		
			int tmp = 1<<c;
			if(op == 1){				
				x = 1;
				y = 1- tmp;
				
				u = v;
				for(int i=0;i<c;i++)
					u <<= 1;				
			}else if(op == 2){
				x = 1- tmp;
				y = -1;
				
				v = -u;
				for(int i=0;i<c;i++)
					v <<= 1;	
			}
		}
		while(!ops.isEmpty()){	
			int op = ops.pop();
			switch(op){
				case 1:{
					int c = ops.pop();
					for(int i=0;i<c;i++){
						if((x & 0x01)==0){
							//if x is even	
							x >>= 1;
						}else{
							y += u;
							int tmp = x - v;
							while(y ==0 || tmp ==0){
								y += (u<<1); //*2
								tmp -= (v<<1); //*2
							};
							x = tmp >>1;
						}
						u <<= 1;
					}
					break;
				}
				case 2:{					
					int c = ops.pop();
					for(int i=0;i<c;i++){
						if((y & 0x01)==0){
							//if y is even								
							y >>= 1;
						}else{
							x += v;
							int tmp = y-u;
							while(x ==0 || tmp ==0){
								x += (v<<1);  //*2
								tmp -= (u<<1);  //*2
							};
							y = tmp >>1;
						}
						v <<= 1;						
					}
					y = -y;
					v = -v;
					break;
				}
				case 3:
					y -= x;
					u += v;
					break;
				case 4:
					x += y;					
					v = u - v;
					
					y = -y;
					break;
			}
		}
		//e.g. gcd(1,4)
		y = (v<0)?-y:y;
		
		for(int i=0;i<k;i++)
			d <<=1;
		return new int[]{d,x,y};
	}
	
	public static void print(int m,int n,int[] extGCDResult){
		m = (m<0)?-m:m;
		n = (n<0)?-n:n;
		System.out.format("extended gcd of %d and %d is: %d = %d.{%d} + %d.{%d}%5s%n",m,n,extGCDResult[0],m,extGCDResult[1],n,extGCDResult[2],(m*extGCDResult[1] + n * extGCDResult[2] == extGCDResult[0])?"OK":"WRONG!!!");		
	}
	
	public static void main(String[] args) {
		int m = -18;
		int n= 12;
		print(m,n,extGCDBinary(m,n));
		
		//co-prime
		m = 15;
		n= 28;
		print(m,n,extGCDBinary(m,n));
				
		m = 6;
		n= 3;
		print(m,n,extGCDBinary(m,n));
		
		m = 6;
		n= 3;
		print(m,n,extGCDBinary(m,n));
		
		m = 6;
		n= 0;
		print(m,n,extGCDBinary(m,n));
		
		m = 0;
		n= 6;
		print(m,n,extGCDBinary(m,n));
		
		m = 0;
		n= 0;
		print(m,n,extGCDBinary(m,n));
		
		m = 1;
		n= 1;
		print(m,n,extGCDBinary(m,n));
		
		m = 3;
		n= 3;
		print(m,n,extGCDBinary(m,n));
		
		m = 2;
		n= 2;
		print(m,n,extGCDBinary(m,n));
		
		m = 1;
		n= 4;
		print(m,n,extGCDBinary(m,n));
		
		m = 4;
		n= 1;
		print(m,n,extGCDBinary(m,n));
		
		m = 10;
		n= 14;
		print(m,n,extGCDBinary(m,n));
		
		m = 14;
		n= 10;
		print(m,n,extGCDBinary(m,n));
		
		m = 10;
		n= 4;
		print(m,n,extGCDBinary(m,n));
		
		
		
		m = 273;
		n= 24;
		print(m,n,extGCDBinary(m,n));
		
		m = 120;
		n= 23;
		print(m,n,extGCDBinary(m,n));
				
		
	}
}