\\ 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] ); }