#######################################
# Euclid.py                           #
#######################################
# Author:  Frank.Schuhmacher@segrids.com
# Date:    2015-03-03
# Licence: GPL
#
# gcd (Euclidian algorithm) and extended gcd
#

import numpy as np

def gcd(X,Y):
	if X < Y:
		(X,Y) = (Y,X)
	while Y != 0:
		(X,Y) = (Y, X%Y)
	return X

# Denote (X_i, Y_i) the values of (X,Y) in step i of the Euclidian gcd, then
#
# 	X_0 = X
# 	Y_0 = Y
# and
# 	X_i = A_i[0,0] * X_0 + A_i[0,1] * Y_0     for i in 0...n
# 	Y_i = A_i[1,0] * X_0 + A_i[1,1] * Y_0     for i in 0...n
#
# for an adequate 2*2 matrix A. In the last step, we get
#
# 	X_n = gcd(X,Y)
# 	Y_n = 0
#
# Using  
# 	X_{i+1} = Y_i
# 	Y_{i+1} = X_i % Y_i = X_i - (X_i//Y_i)*Y_i 
#
# derive
#	A_{i+1}[0,0] = A_i[1,0]
#	A_{i+1}[0,1] = A_i[1,1]
#	A_{i+1}[1,0] = A_i[0,0] - (X_i//Y_i)*A_i[1,0]
#	A_{i+1}[1,1] = A_i[0,1] - (X_i//Y_i)*A_i[1,1],
# i.e.
#	A_0 = Id,
# and                
#	A_{i+1} = [[0 ,  1       ],
#	           [1,  -X_i//Y_i]] * A_i
#
check = True

def extended_gcd(intX,intY):
	X = intX
	Y = intY
	if X < Y:
		(X,Y) = (Y,X)
		swap = True
	else:
		swap = False
	A = np.matrix([[1,0],[0,1]], dtype='int')
	while Y != 0:
		A = np.matrix([[0,1],[1,-(X//Y)]], dtype='int') * A
		(X,Y) = (Y, X%Y)
	if swap:
		(Xdash, Ydash) = (A[0,1], A[0,0])
	else:
		(Xdash, Ydash) = (A[0,0], A[0,1])
	if check:
		print('gcd(%d' %intX + ',%d)' %intY + ' = %d' %X)
		print('%d' %Xdash + '*%d + ' %intX +\
                      '%d' %Ydash + '*%d' %intY +\
                      ' = %d' %(intX * Xdash + intY * Ydash))
	return (Xdash, Ydash)