1 #include "os.h"
2 #include <mp.h>
3 
4 /* extended euclid */
5 /* */
6 /* For a and b it solves, d = gcd(a,b) and finds x and y s.t. */
7 /* ax + by = d */
8 /* */
9 /* Handbook of Applied Cryptography, Menezes et al, 1997, pg 67 */
10 
11 void
12 mpeuclid(mpint *a, mpint *b, mpint *d, mpint *x, mpint *y)
13 {
14 	mpint *tmp, *x0, *x1, *x2, *y0, *y1, *y2, *q, *r;
15 
16 	if(a->sign<0 || b->sign<0)
17 		sysfatal("mpeuclid: negative arg");
18 
19 	if(mpcmp(a, b) < 0){
20 		tmp = a;
21 		a = b;
22 		b = tmp;
23 		tmp = x;
24 		x = y;
25 		y = tmp;
26 	}
27 
28 	if(b->top == 0){
29 		mpassign(a, d);
30 		mpassign(mpone, x);
31 		mpassign(mpzero, y);
32 		return;
33 	}
34 
35 	a = mpcopy(a);
36 	b = mpcopy(b);
37 	x0 = mpnew(0);
38 	x1 = mpcopy(mpzero);
39 	x2 = mpcopy(mpone);
40 	y0 = mpnew(0);
41 	y1 = mpcopy(mpone);
42 	y2 = mpcopy(mpzero);
43 	q = mpnew(0);
44 	r = mpnew(0);
45 
46 	while(b->top != 0 && b->sign > 0){
47 		/* q = a/b */
48 		/* r = a mod b */
49 		mpdiv(a, b, q, r);
50 		/* x0 = x2 - qx1 */
51 		mpmul(q, x1, x0);
52 		mpsub(x2, x0, x0);
53 		/* y0 = y2 - qy1 */
54 		mpmul(q, y1, y0);
55 		mpsub(y2, y0, y0);
56 		/* rotate values */
57 		tmp = a;
58 		a = b;
59 		b = r;
60 		r = tmp;
61 		tmp = x2;
62 		x2 = x1;
63 		x1 = x0;
64 		x0 = tmp;
65 		tmp = y2;
66 		y2 = y1;
67 		y1 = y0;
68 		y0 = tmp;
69 	}
70 
71 	mpassign(a, d);
72 	mpassign(x2, x);
73 	mpassign(y2, y);
74 
75 	mpfree(x0);
76 	mpfree(x1);
77 	mpfree(x2);
78 	mpfree(y0);
79 	mpfree(y1);
80 	mpfree(y2);
81 	mpfree(q);
82 	mpfree(r);
83 	mpfree(a);
84 	mpfree(b);
85 }