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