Mercurial > dropbear

comparison bn_mp_gcd.c @ 2:86e0b50a9b58 libtommath-orig ltm-0.30-orig

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
ltm 0.30 orig import
author Matt Johnston <matt@ucc.asn.au>
date Mon, 31 May 2004 18:25:22 +0000
parents
children d29b64170cf0
comparison
equal deleted inserted replaced
-1:000000000000 2:86e0b50a9b58
1 /* LibTomMath, multiple-precision integer library -- Tom St Denis
2 *
3 * LibTomMath is a library that provides multiple-precision
4 * integer arithmetic as well as number theoretic functionality.
5 *
6 * The library was designed directly after the MPI library by
7 * Michael Fromberger but has been written from scratch with
8 * additional optimizations in place.
9 *
10 * The library is free for all purposes without any express
11 * guarantee it works.
12 *
13 * Tom St Denis, [email protected], http://math.libtomcrypt.org
14 */
15 #include <tommath.h>
16
17 /* Greatest Common Divisor using the binary method */
18 int mp_gcd (mp_int * a, mp_int * b, mp_int * c)
19 {
20 mp_int u, v;
21 int k, u_lsb, v_lsb, res;
22
23 /* either zero than gcd is the largest */
24 if (mp_iszero (a) == 1 && mp_iszero (b) == 0) {
25 return mp_abs (b, c);
26 }
27 if (mp_iszero (a) == 0 && mp_iszero (b) == 1) {
28 return mp_abs (a, c);
29 }
30
31 /* optimized. At this point if a == 0 then
32 * b must equal zero too
33 */
34 if (mp_iszero (a) == 1) {
35 mp_zero(c);
36 return MP_OKAY;
37 }
38
39 /* get copies of a and b we can modify */
40 if ((res = mp_init_copy (&u, a)) != MP_OKAY) {
41 return res;
42 }
43
44 if ((res = mp_init_copy (&v, b)) != MP_OKAY) {
45 goto __U;
46 }
47
48 /* must be positive for the remainder of the algorithm */
49 u.sign = v.sign = MP_ZPOS;
50
51 /* B1. Find the common power of two for u and v */
52 u_lsb = mp_cnt_lsb(&u);
53 v_lsb = mp_cnt_lsb(&v);
54 k = MIN(u_lsb, v_lsb);
55
56 if (k > 0) {
57 /* divide the power of two out */
58 if ((res = mp_div_2d(&u, k, &u, NULL)) != MP_OKAY) {
59 goto __V;
60 }
61
62 if ((res = mp_div_2d(&v, k, &v, NULL)) != MP_OKAY) {
63 goto __V;
64 }
65 }
66
67 /* divide any remaining factors of two out */
68 if (u_lsb != k) {
69 if ((res = mp_div_2d(&u, u_lsb - k, &u, NULL)) != MP_OKAY) {
70 goto __V;
71 }
72 }
73
74 if (v_lsb != k) {
75 if ((res = mp_div_2d(&v, v_lsb - k, &v, NULL)) != MP_OKAY) {
76 goto __V;
77 }
78 }
79
80 while (mp_iszero(&v) == 0) {
81 /* make sure v is the largest */
82 if (mp_cmp_mag(&u, &v) == MP_GT) {
83 /* swap u and v to make sure v is >= u */
84 mp_exch(&u, &v);
85 }
86
87 /* subtract smallest from largest */
88 if ((res = s_mp_sub(&v, &u, &v)) != MP_OKAY) {
89 goto __V;
90 }
91
92 /* Divide out all factors of two */
93 if ((res = mp_div_2d(&v, mp_cnt_lsb(&v), &v, NULL)) != MP_OKAY) {
94 goto __V;
95 }
96 }
97
98 /* multiply by 2**k which we divided out at the beginning */
99 if ((res = mp_mul_2d (&u, k, c)) != MP_OKAY) {
100 goto __V;
101 }
102 c->sign = MP_ZPOS;
103 res = MP_OKAY;
104 __V:mp_clear (&u);
105 __U:mp_clear (&v);
106 return res;
107 }