libs/libtommath/bn_mp_gcd.c

   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 */
  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.5 $ */
 105 /* $Date: 2006/12/28 01:25:13 $ */