xref: /plan9/sys/src/cmd/unix/drawterm/libmp/mpextendedgcd.c (revision 8ccd4a6360d974db7bd7bbd4f37e7018419ea908)
1*8ccd4a63SDavid du Colombier #include "os.h"
2*8ccd4a63SDavid du Colombier #include <mp.h>
3*8ccd4a63SDavid du Colombier 
4*8ccd4a63SDavid du Colombier #define iseven(a)	(((a)->p[0] & 1) == 0)
5*8ccd4a63SDavid du Colombier 
6*8ccd4a63SDavid du Colombier // extended binary gcd
7*8ccd4a63SDavid du Colombier //
8*8ccd4a63SDavid du Colombier // For a anv b it solves, v = gcd(a,b) and finds x and y s.t.
9*8ccd4a63SDavid du Colombier // ax + by = v
10*8ccd4a63SDavid du Colombier //
11*8ccd4a63SDavid du Colombier // Handbook of Applied Cryptography, Menezes et al, 1997, pg 608.
12*8ccd4a63SDavid du Colombier void
mpextendedgcd(mpint * a,mpint * b,mpint * v,mpint * x,mpint * y)13*8ccd4a63SDavid du Colombier mpextendedgcd(mpint *a, mpint *b, mpint *v, mpint *x, mpint *y)
14*8ccd4a63SDavid du Colombier {
15*8ccd4a63SDavid du Colombier 	mpint *u, *A, *B, *C, *D;
16*8ccd4a63SDavid du Colombier 	int g;
17*8ccd4a63SDavid du Colombier 
18*8ccd4a63SDavid du Colombier 	if(a->top == 0){
19*8ccd4a63SDavid du Colombier 		mpassign(b, v);
20*8ccd4a63SDavid du Colombier 		mpassign(mpone, y);
21*8ccd4a63SDavid du Colombier 		mpassign(mpzero, x);
22*8ccd4a63SDavid du Colombier 		return;
23*8ccd4a63SDavid du Colombier 	}
24*8ccd4a63SDavid du Colombier 	if(b->top == 0){
25*8ccd4a63SDavid du Colombier 		mpassign(a, v);
26*8ccd4a63SDavid du Colombier 		mpassign(mpone, x);
27*8ccd4a63SDavid du Colombier 		mpassign(mpzero, y);
28*8ccd4a63SDavid du Colombier 		return;
29*8ccd4a63SDavid du Colombier 	}
30*8ccd4a63SDavid du Colombier 
31*8ccd4a63SDavid du Colombier 	g = 0;
32*8ccd4a63SDavid du Colombier 	a = mpcopy(a);
33*8ccd4a63SDavid du Colombier 	b = mpcopy(b);
34*8ccd4a63SDavid du Colombier 
35*8ccd4a63SDavid du Colombier 	while(iseven(a) && iseven(b)){
36*8ccd4a63SDavid du Colombier 		mpright(a, 1, a);
37*8ccd4a63SDavid du Colombier 		mpright(b, 1, b);
38*8ccd4a63SDavid du Colombier 		g++;
39*8ccd4a63SDavid du Colombier 	}
40*8ccd4a63SDavid du Colombier 
41*8ccd4a63SDavid du Colombier 	u = mpcopy(a);
42*8ccd4a63SDavid du Colombier 	mpassign(b, v);
43*8ccd4a63SDavid du Colombier 	A = mpcopy(mpone);
44*8ccd4a63SDavid du Colombier 	B = mpcopy(mpzero);
45*8ccd4a63SDavid du Colombier 	C = mpcopy(mpzero);
46*8ccd4a63SDavid du Colombier 	D = mpcopy(mpone);
47*8ccd4a63SDavid du Colombier 
48*8ccd4a63SDavid du Colombier 	for(;;) {
49*8ccd4a63SDavid du Colombier //		print("%B %B %B %B %B %B\n", u, v, A, B, C, D);
50*8ccd4a63SDavid du Colombier 		while(iseven(u)){
51*8ccd4a63SDavid du Colombier 			mpright(u, 1, u);
52*8ccd4a63SDavid du Colombier 			if(!iseven(A) || !iseven(B)) {
53*8ccd4a63SDavid du Colombier 				mpadd(A, b, A);
54*8ccd4a63SDavid du Colombier 				mpsub(B, a, B);
55*8ccd4a63SDavid du Colombier 			}
56*8ccd4a63SDavid du Colombier 			mpright(A, 1, A);
57*8ccd4a63SDavid du Colombier 			mpright(B, 1, B);
58*8ccd4a63SDavid du Colombier 		}
59*8ccd4a63SDavid du Colombier 
60*8ccd4a63SDavid du Colombier //		print("%B %B %B %B %B %B\n", u, v, A, B, C, D);
61*8ccd4a63SDavid du Colombier 		while(iseven(v)){
62*8ccd4a63SDavid du Colombier 			mpright(v, 1, v);
63*8ccd4a63SDavid du Colombier 			if(!iseven(C) || !iseven(D)) {
64*8ccd4a63SDavid du Colombier 				mpadd(C, b, C);
65*8ccd4a63SDavid du Colombier 				mpsub(D, a, D);
66*8ccd4a63SDavid du Colombier 			}
67*8ccd4a63SDavid du Colombier 			mpright(C, 1, C);
68*8ccd4a63SDavid du Colombier 			mpright(D, 1, D);
69*8ccd4a63SDavid du Colombier 		}
70*8ccd4a63SDavid du Colombier 
71*8ccd4a63SDavid du Colombier //		print("%B %B %B %B %B %B\n", u, v, A, B, C, D);
72*8ccd4a63SDavid du Colombier 		if(mpcmp(u, v) >= 0){
73*8ccd4a63SDavid du Colombier 			mpsub(u, v, u);
74*8ccd4a63SDavid du Colombier 			mpsub(A, C, A);
75*8ccd4a63SDavid du Colombier 			mpsub(B, D, B);
76*8ccd4a63SDavid du Colombier 		} else {
77*8ccd4a63SDavid du Colombier 			mpsub(v, u, v);
78*8ccd4a63SDavid du Colombier 			mpsub(C, A, C);
79*8ccd4a63SDavid du Colombier 			mpsub(D, B, D);
80*8ccd4a63SDavid du Colombier 		}
81*8ccd4a63SDavid du Colombier 
82*8ccd4a63SDavid du Colombier 		if(u->top == 0)
83*8ccd4a63SDavid du Colombier 			break;
84*8ccd4a63SDavid du Colombier 
85*8ccd4a63SDavid du Colombier 	}
86*8ccd4a63SDavid du Colombier 	mpassign(C, x);
87*8ccd4a63SDavid du Colombier 	mpassign(D, y);
88*8ccd4a63SDavid du Colombier 	mpleft(v, g, v);
89*8ccd4a63SDavid du Colombier 
90*8ccd4a63SDavid du Colombier 	mpfree(A);
91*8ccd4a63SDavid du Colombier 	mpfree(B);
92*8ccd4a63SDavid du Colombier 	mpfree(C);
93*8ccd4a63SDavid du Colombier 	mpfree(D);
94*8ccd4a63SDavid du Colombier 	mpfree(u);
95*8ccd4a63SDavid du Colombier 	mpfree(a);
96*8ccd4a63SDavid du Colombier 	mpfree(b);
97*8ccd4a63SDavid du Colombier }
98