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