crypto/external/bsd/heimdal/dist/lib/hcrypto/libtommath/bn_mp_gcd.c

   1 /*      $NetBSD: bn_mp_gcd.c,v 1.1.1.2 2014/04/24 12:45:31 pettai 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@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  */