Home | History | Annotate | Line # | Download | only in libtommath 
      1 /*	$NetBSD: bn_mp_gcd.c,v 1.2 2017/01/28 21:31:47 christos 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 (at) 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