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