mp_gcd.c

   1 #include "tommath_private.h"
   2 #ifdef MP_GCD_C
   3 /* LibTomMath, multiple-precision integer library -- Tom St Denis */
   4 /* SPDX-License-Identifier: Unlicense */
   5
   6 /* Greatest Common Divisor using the binary method */
   7 mp_err mp_gcd(const mp_int *a, const mp_int *b, mp_int *c)
   8 {
   9    mp_int  u, v;
  10    int     k, u_lsb, v_lsb;
  11    mp_err err;
  12
  13    /* either zero than gcd is the largest */
  14    if (mp_iszero(a)) {
  15       return mp_abs(b, c);
  16    }
  17    if (mp_iszero(b)) {
  18       return mp_abs(a, c);
  19    }
  20
  21    /* get copies of a and b we can modify */
  22    if ((err = mp_init_copy(&u, a)) != MP_OKAY) {
  23       return err;
  24    }
  25
  26    if ((err = mp_init_copy(&v, b)) != MP_OKAY) {
  27       goto LBL_U;
  28    }
  29
  30    /* must be positive for the remainder of the algorithm */
  31    u.sign = v.sign = MP_ZPOS;
  32
  33    /* B1.  Find the common power of two for u and v */
  34    u_lsb = mp_cnt_lsb(&u);
  35    v_lsb = mp_cnt_lsb(&v);
  36    k     = MP_MIN(u_lsb, v_lsb);
  37
  38    if (k > 0) {
  39       /* divide the power of two out */
  40       if ((err = mp_div_2d(&u, k, &u, NULL)) != MP_OKAY) {
  41          goto LBL_V;
  42       }
  43
  44       if ((err = mp_div_2d(&v, k, &v, NULL)) != MP_OKAY) {
  45          goto LBL_V;
  46       }
  47    }
  48
  49    /* divide any remaining factors of two out */
  50    if (u_lsb != k) {
  51       if ((err = mp_div_2d(&u, u_lsb - k, &u, NULL)) != MP_OKAY) {
  52          goto LBL_V;
  53       }
  54    }
  55
  56    if (v_lsb != k) {
  57       if ((err = mp_div_2d(&v, v_lsb - k, &v, NULL)) != MP_OKAY) {
  58          goto LBL_V;
  59       }
  60    }
  61
  62    while (!mp_iszero(&v)) {
  63       /* make sure v is the largest */
  64       if (mp_cmp_mag(&u, &v) == MP_GT) {
  65          /* swap u and v to make sure v is >= u */
  66          mp_exch(&u, &v);
  67       }
  68
  69       /* subtract smallest from largest */
  70       if ((err = s_mp_sub(&v, &u, &v)) != MP_OKAY) {
  71          goto LBL_V;
  72       }
  73
  74       /* Divide out all factors of two */
  75       if ((err = mp_div_2d(&v, mp_cnt_lsb(&v), &v, NULL)) != MP_OKAY) {
  76          goto LBL_V;
  77       }
  78    }
  79
  80    /* multiply by 2**k which we divided out at the beginning */
  81    if ((err = mp_mul_2d(&u, k, c)) != MP_OKAY) {
  82       goto LBL_V;
  83    }
  84    c->sign = MP_ZPOS;
  85    err = MP_OKAY;
  86 LBL_V:
  87    mp_clear(&u);
  88 LBL_U:
  89    mp_clear(&v);
  90    return err;
  91 }
  92 #endif