Home | History | Annotate | Download | only in libtommath 
      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 (at) gmail.com, http://math.libtomcrypt.com
     16  */
     17 
     18 /* Greatest Common Divisor using the binary method */
     19 int 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;
     97 LBL_V:mp_clear (&u);
     98 LBL_U:mp_clear (&v);
     99   return res;
    100 }
    101 #endif
    102 
    103 /* $Source: /cvs/libtom/libtommath/bn_mp_gcd.c,v $ */
    104 /* $Revision: 1.4 $ */
    105 /* $Date: 2006/03/31 14:18:44 $ */
    106