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