1#include <tommath.h>
2#ifdef BN_MP_GCD_C
3/* LibTomMath, multiple-precision integer library -- Tom St Denis
4 *
5 * LibTomMath is a library that provides multiple-precision
6 * integer arithmetic as well as number theoretic functionality.
7 *
8 * The library was designed directly after the MPI library by
9 * Michael Fromberger but has been written from scratch with
10 * additional optimizations in place.
11 *
12 * The library is free for all purposes without any express
13 * guarantee it works.
14 *
15 * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
16 */
17
18/* Greatest Common Divisor using the binary method */
19int mp_gcd (mp_int * a, mp_int * b, mp_int * c)
20{
21  mp_int  u, v;
22  int     k, u_lsb, v_lsb, res;
23
24  /* either zero than gcd is the largest */
25  if (mp_iszero (a) == MP_YES) {
26    return mp_abs (b, c);
27  }
28  if (mp_iszero (b) == MP_YES) {
29    return mp_abs (a, c);
30  }
31
32  /* get copies of a and b we can modify */
33  if ((res = mp_init_copy (&u, a)) != MP_OKAY) {
34    return res;
35  }
36
37  if ((res = mp_init_copy (&v, b)) != MP_OKAY) {
38    goto LBL_U;
39  }
40
41  /* must be positive for the remainder of the algorithm */
42  u.sign = v.sign = MP_ZPOS;
43
44  /* B1.  Find the common power of two for u and v */
45  u_lsb = mp_cnt_lsb(&u);
46  v_lsb = mp_cnt_lsb(&v);
47  k     = MIN(u_lsb, v_lsb);
48
49  if (k > 0) {
50     /* divide the power of two out */
51     if ((res = mp_div_2d(&u, k, &u, NULL)) != MP_OKAY) {
52        goto LBL_V;
53     }
54
55     if ((res = mp_div_2d(&v, k, &v, NULL)) != MP_OKAY) {
56        goto LBL_V;
57     }
58  }
59
60  /* divide any remaining factors of two out */
61  if (u_lsb != k) {
62     if ((res = mp_div_2d(&u, u_lsb - k, &u, NULL)) != MP_OKAY) {
63        goto LBL_V;
64     }
65  }
66
67  if (v_lsb != k) {
68     if ((res = mp_div_2d(&v, v_lsb - k, &v, NULL)) != MP_OKAY) {
69        goto LBL_V;
70     }
71  }
72
73  while (mp_iszero(&v) == 0) {
74     /* make sure v is the largest */
75     if (mp_cmp_mag(&u, &v) == MP_GT) {
76        /* swap u and v to make sure v is >= u */
77        mp_exch(&u, &v);
78     }
79
80     /* subtract smallest from largest */
81     if ((res = s_mp_sub(&v, &u, &v)) != MP_OKAY) {
82        goto LBL_V;
83     }
84
85     /* Divide out all factors of two */
86     if ((res = mp_div_2d(&v, mp_cnt_lsb(&v), &v, NULL)) != MP_OKAY) {
87        goto LBL_V;
88     }
89  }
90
91  /* multiply by 2**k which we divided out at the beginning */
92  if ((res = mp_mul_2d (&u, k, c)) != MP_OKAY) {
93     goto LBL_V;
94  }
95  c->sign = MP_ZPOS;
96  res = MP_OKAY;
97LBL_V:mp_clear (&u);
98LBL_U:mp_clear (&v);
99  return res;
100}
101#endif
102
103/* $Source$ */
104/* $Revision$ */
105/* $Date$ */
106