xref: /plan9/sys/src/cmd/unix/drawterm/libmp/mpextendedgcd.c (revision 8ccd4a6360d974db7bd7bbd4f37e7018419ea908)
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
mpextendedgcd(mpint * a,mpint * b,mpint * v,mpint * x,mpint * y)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->top == 0){
19 		mpassign(b, v);
20 		mpassign(mpone, y);
21 		mpassign(mpzero, x);
22 		return;
23 	}
24 	if(b->top == 0){
25 		mpassign(a, v);
26 		mpassign(mpone, x);
27 		mpassign(mpzero, y);
28 		return;
29 	}
30 
31 	g = 0;
32 	a = mpcopy(a);
33 	b = mpcopy(b);
34 
35 	while(iseven(a) && iseven(b)){
36 		mpright(a, 1, a);
37 		mpright(b, 1, b);
38 		g++;
39 	}
40 
41 	u = mpcopy(a);
42 	mpassign(b, v);
43 	A = mpcopy(mpone);
44 	B = mpcopy(mpzero);
45 	C = mpcopy(mpzero);
46 	D = mpcopy(mpone);
47 
48 	for(;;) {
49 //		print("%B %B %B %B %B %B\n", u, v, A, B, C, D);
50 		while(iseven(u)){
51 			mpright(u, 1, u);
52 			if(!iseven(A) || !iseven(B)) {
53 				mpadd(A, b, A);
54 				mpsub(B, a, B);
55 			}
56 			mpright(A, 1, A);
57 			mpright(B, 1, B);
58 		}
59 
60 //		print("%B %B %B %B %B %B\n", u, v, A, B, C, D);
61 		while(iseven(v)){
62 			mpright(v, 1, v);
63 			if(!iseven(C) || !iseven(D)) {
64 				mpadd(C, b, C);
65 				mpsub(D, a, D);
66 			}
67 			mpright(C, 1, C);
68 			mpright(D, 1, D);
69 		}
70 
71 //		print("%B %B %B %B %B %B\n", u, v, A, B, C, D);
72 		if(mpcmp(u, v) >= 0){
73 			mpsub(u, v, u);
74 			mpsub(A, C, A);
75 			mpsub(B, D, B);
76 		} else {
77 			mpsub(v, u, v);
78 			mpsub(C, A, C);
79 			mpsub(D, B, D);
80 		}
81 
82 		if(u->top == 0)
83 			break;
84 
85 	}
86 	mpassign(C, x);
87 	mpassign(D, y);
88 	mpleft(v, g, v);
89 
90 	mpfree(A);
91 	mpfree(B);
92 	mpfree(C);
93 	mpfree(D);
94 	mpfree(u);
95 	mpfree(a);
96 	mpfree(b);
97 }
98