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