userapps/opensource/sshd/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@iahu.ca, http://math.libtomcrypt.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) == 1 && mp_iszero (b) == 0) {
  26     return mp_abs (b, c);
  27   }
  28   if (mp_iszero (a) == 0 && mp_iszero (b) == 1) {
  29     return mp_abs (a, c);
  30   }
  31
  32   /* optimized.  At this point if a == 0 then
  33    * b must equal zero too
  34    */
  35   if (mp_iszero (a) == 1) {
  36     mp_zero(c);
  37     return MP_OKAY;
  38   }
  39
  40   /* get copies of a and b we can modify */
  41   if ((res = mp_init_copy (&u, a)) != MP_OKAY) {
  42     return res;
  43   }
  44
  45   if ((res = mp_init_copy (&v, b)) != MP_OKAY) {
  46     goto LBL_U;
  47   }
  48
  49   /* must be positive for the remainder of the algorithm */
  50   u.sign = v.sign = MP_ZPOS;
  51
  52   /* B1.  Find the common power of two for u and v */
  53   u_lsb = mp_cnt_lsb(&u);
  54   v_lsb = mp_cnt_lsb(&v);
  55   k     = MIN(u_lsb, v_lsb);
  56
  57   if (k > 0) {
  58      /* divide the power of two out */
  59      if ((res = mp_div_2d(&u, k, &u, NULL)) != MP_OKAY) {
  60         goto LBL_V;
  61      }
  62
  63      if ((res = mp_div_2d(&v, k, &v, NULL)) != MP_OKAY) {
  64         goto LBL_V;
  65      }
  66   }
  67
  68   /* divide any remaining factors of two out */
  69   if (u_lsb != k) {
  70      if ((res = mp_div_2d(&u, u_lsb - k, &u, NULL)) != MP_OKAY) {
  71         goto LBL_V;
  72      }
  73   }
  74
  75   if (v_lsb != k) {
  76      if ((res = mp_div_2d(&v, v_lsb - k, &v, NULL)) != MP_OKAY) {
  77         goto LBL_V;
  78      }
  79   }
  80
  81   while (mp_iszero(&v) == 0) {
  82      /* make sure v is the largest */
  83      if (mp_cmp_mag(&u, &v) == MP_GT) {
  84         /* swap u and v to make sure v is >= u */
  85         mp_exch(&u, &v);
  86      }
  87
  88      /* subtract smallest from largest */
  89      if ((res = s_mp_sub(&v, &u, &v)) != MP_OKAY) {
  90         goto LBL_V;
  91      }
  92
  93      /* Divide out all factors of two */
  94      if ((res = mp_div_2d(&v, mp_cnt_lsb(&v), &v, NULL)) != MP_OKAY) {
  95         goto LBL_V;
  96      }
  97   }
  98
  99   /* multiply by 2**k which we divided out at the beginning */
 100   if ((res = mp_mul_2d (&u, k, c)) != MP_OKAY) {
 101      goto LBL_V;
 102   }
 103   c->sign = MP_ZPOS;
 104   res = MP_OKAY;
 105 LBL_V:mp_clear (&u);
 106 LBL_U:mp_clear (&v);
 107   return res;
 108 }
 109 #endif