\\ These pari programs use the implementation of arithmetic
\\ operations in pari to compute GCD and the associated
\\ solution of the linear form using different techniques.

euclidgcd(a,b)=
{
	local(t);

	while(b != 0, \\ I.e. while b is non-zero
		t=a%b; \\ This is a mod b written in pari
		a=b;
		b=t;
	);
	a; \\ The answer.
}

twoval(n)=
{
	local(v);

	v=0;while(n%2 == 0,
		n=n\2;
		v++;\\ increment
	);
	[v,n];
}

binarygcd(a,b)=
{
	local(k,m,n);

	if(b == 0, return(a)); \\ Don't remove infinitely many
		\\ powers of 2 from 0!
	if(a == 0, return(b)); \\ Same for a
	m = twoval(a); n = twoval(b);
	
	k = min(m[1],n[1]); a = m[2]; b = n[2];
	\\ a and b are now odd 
	c = a - b;
	while(c != 0,
		c = twoval(c)[2];
		if( c < 0,
			b = -c,
			a = c
		);
		c = a - b;
	);
	2^k*a;
}

\\ For the Lehmer GCD we need a "wordsize" in bits
wsize=16;

\\ And a way of getting the length of an integer in bits
wlen(n)=
{
	length(binary(n));
}

\\ And a way of getting top word of an integer of bitsize p
tword(n,p)=
{
	if(p > wsize,
		n\( 2^(p-wsize) ),
		n
	);
}

\\ The transformation matrix a la extended Euclid
\\ note u,v are digits
findmat(u,v)=
{
	local(M,x,y,t);

	M=[[1,0],[0,1]];
	while(1, 
		if(v + M[2][2] == 0, return(M));
		if(v + M[1][2] == 0, return(M));
		x=(u+M[2][1])\(v+M[2][2]);y=(u+M[1][1])\(v+M[1][2]);
		if(x != y, return(M));
		t = u-x*v; u = v; v = t;
		t = M[1][1]-x*M[1][2]; M[1][1]=M[1][2]; M[1][2]=t;		
		t = M[2][1]-x*M[2][2]; M[2][1]=M[2][2]; M[2][2]=t;		
	);
}

lehmergcd(a,b)=
{
	local(M,w,p,t,u,v);
	w=2^wsize;

	if(a < b, t = b; b = a; a = t);
	while(b > 0,
		if(b < w, return(euclidgcd(a,b)));
		p=wlen(a);
		u=tword(a,p);v=tword(b,p);
		M=findmat(u,v);
		if(M[2][1] == 0,
			t = a%b; a = b; b = t,
			t = M[1][1]*a + M[2][1]*b; b = M[1][2]*a + M[2][2]*b; a = t
		);
	);
	a;
}

\\ Extending the Euclid GCD is easy
euclidegcd(a,b)=
{
	local(t,q,A,B);

	A=[1,a];B=[0,b];
	while(B[2] != 0, 
		q=A[2]\B[2];
		t=A-q*B;
		A=B;
		B=t;
	);
	[A[1],(A[2]-a*A[1])/b]; \\ The answer.
}

\\ We need this subroutine for Lehmer extended
subeuclidegcd(A,B)=
{
	local(t,q);

	while(B[2] != 0, 
		q=A[2]\B[2];
		t=A-q*B;
		A=B;
		B=t;
	);
	[A[1],(A[2]-a*A[1])/b]; \\ The answer.
}

\\ Lehmer extends quite easily as well
lehmeregcd(a,b)=
{
	local(M,w,p,t,A,B);
	w=2^wsize;
	
	if(a < b,
		A=[0,b];B=[1,a],
		A=[1,a];B=[0,b]
	);	

	while( B[2] != 0,
		if(B[2] < w, return(subeuclidegcd(A,B)));
		p=wlen(A[2]);
		M=findmat( tword(A[2],p), tword(B[2],p));
		if(M[2][1] == 0,
			M = [[0,1],[1,-A[2]\B[2]]];
		);
		t = M[1][1]*A + M[2][1]*B; B = M[1][2]*A + M[2][2]*B; A = t
	);
	[A[1],(A[2]-a*A[1])/b]; \\ The answer.
}

\\ The slightly more tricky Extended Binary GCD
\\ We have to ensure that both components of the vector
\\ are divisible by two.
binaryegcd(a,b)=
{
	local(f,A,B,c,d);

	if(b == 0, return([1,0])); \\ Don't remove infinitely many
		\\ powers of 2 from 0!
	if(a == 0, return([0,1])); \\ Same for a

	while( (a%2 + b%2) == 0,
		a = a/2; b = b/2;
	);
	\\ One of a and b is now odd 

	if( b%2 == 0,
		d = a; c = b; f = 1,
		c = a; d = b; f = 0
	);
	\\ Now d is certainly odd

	B=[0,d]; A = [1,c];
	if( c%2 == 0,
		C = A,
		C = A - B
	);
	\\ C[2] is definitely even
	while( C[2] != 0,
		until( C[2]%2 == 1,
			if( C[1]%2 == 0,
				C = [ C[1]\2, C[2]\2 ],
				if ( C[1] > 0,
					 C = [ (C[1]-d)\2, C[2]\2 ],
					 C = [ (C[1]+d)\2, C[2]\2 ])
			); \\ C[2]-C[1]*c remains divisible by d
		);
		\\ Now C[2] is odd
		if( C[2] < 0, \\ So A[2] < B[2] and c was not even
			B = -C, \\ The first term is to get positivity of B[1]
			A = C
		);
		C = A - B;
	);
	if( f == 1, \\ Interchange case
		[ (A[2]-c*A[1])/d, A[1] ],
		[ A[1], (A[2]-c*A[1])/d]
	);
}