userapps/opensource/sshd/libtommath/bn_mp_gcd.c

   1 /* LibTomMath, multiple-precision integer library -- Tom St Denis
   2  *
   3  * LibTomMath is library that provides for multiple-precision
   4  * integer arithmetic as well as number theoretic functionality.
   5  *
   6  * The library is designed directly after the MPI library by
   7  * Michael Fromberger but has been written from scratch with
   8  * additional optimizations in place.
   9  *
  10  * The library is free for all purposes without any express
  11  * guarantee it works.
  12  *
  13  * Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
  14  */
  15 #include <tommath.h>
  16
  17 /* Greatest Common Divisor using the binary method [Algorithm B, page 338, vol2 of TAOCP]
  18  */
  19 int
  20 mp_gcd (mp_int * a, mp_int * b, mp_int * c)
  21 {
  22   mp_int  u, v, t;
  23   int     k, res, neg;
  24
  25   /* either zero than gcd is the largest */
  26   if (mp_iszero (a) == 1 && mp_iszero (b) == 0) {
  27     return mp_copy (b, c);
  28   }
  29   if (mp_iszero (a) == 0 && mp_iszero (b) == 1) {
  30     return mp_copy (a, c);
  31   }
  32   if (mp_iszero (a) == 1 && mp_iszero (b) == 1) {
  33     mp_set (c, 1);
  34     return MP_OKAY;
  35   }
  36
  37   /* if both are negative they share (-1) as a common divisor */
  38   neg = (a->sign == b->sign) ? a->sign : MP_ZPOS;
  39
  40   if ((res = mp_init_copy (&u, a)) != MP_OKAY) {
  41     return res;
  42   }
  43
  44   if ((res = mp_init_copy (&v, b)) != MP_OKAY) {
  45     goto __U;
  46   }
  47
  48   /* must be positive for the remainder of the algorithm */
  49   u.sign = v.sign = MP_ZPOS;
  50
  51   if ((res = mp_init (&t)) != MP_OKAY) {
  52     goto __V;
  53   }
  54
  55   /* B1.  Find power of two */
  56   k = 0;
  57   while (mp_iseven(&u) == 1 && mp_iseven(&v) == 1) {
  58     ++k;
  59     if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
  60       goto __T;
  61     }
  62     if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
  63       goto __T;
  64     }
  65   }
  66
  67   /* B2.  Initialize */
  68   if (mp_isodd(&u) == 1) {
  69     /* t = -v */
  70     if ((res = mp_copy (&v, &t)) != MP_OKAY) {
  71       goto __T;
  72     }
  73     t.sign = MP_NEG;
  74   } else {
  75     /* t = u */
  76     if ((res = mp_copy (&u, &t)) != MP_OKAY) {
  77       goto __T;
  78     }
  79   }
  80
  81   do {
  82     /* B3 (and B4).  Halve t, if even */
  83     while (t.used != 0 && mp_iseven(&t) == 1) {
  84       if ((res = mp_div_2 (&t, &t)) != MP_OKAY) {
  85         goto __T;
  86       }
  87     }
  88
  89     /* B5.  if t>0 then u=t otherwise v=-t */
  90     if (t.used != 0 && t.sign != MP_NEG) {
  91       if ((res = mp_copy (&t, &u)) != MP_OKAY) {
  92         goto __T;
  93       }
  94     } else {
  95       if ((res = mp_copy (&t, &v)) != MP_OKAY) {
  96         goto __T;
  97       }
  98       v.sign = (v.sign == MP_ZPOS) ? MP_NEG : MP_ZPOS;
  99     }
 100
 101     /* B6.  t = u - v, if t != 0 loop otherwise terminate */
 102     if ((res = mp_sub (&u, &v, &t)) != MP_OKAY) {
 103       goto __T;
 104     }
 105   } while (mp_iszero(&t) == 0);
 106
 107   /* multiply by 2^k which we divided out at the beginning */
 108   if ((res = mp_mul_2d (&u, k, &u)) != MP_OKAY) {
 109     goto __T;
 110   }
 111
 112   mp_exch (&u, c);
 113   c->sign = neg;
 114   res = MP_OKAY;
 115 __T:mp_clear (&t);
 116 __V:mp_clear (&u);
 117 __U:mp_clear (&v);
 118   return res;
 119 }