Line data    Source code

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