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
mpeuclid(mpint * a,mpint * b,mpint * d,mpint * x,mpint * y)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 }
86