Skip to content
This repository

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Fetching contributors…

Cannot retrieve contributors at this time

file 119 lines (105 sloc) 2.853 kb
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119
/* LibTomMath, multiple-precision integer library -- Tom St Denis
*
* LibTomMath is library that provides for multiple-precision
* integer arithmetic as well as number theoretic functionality.
*
* The library is designed directly after the MPI library by
* Michael Fromberger but has been written from scratch with
* additional optimizations in place.
*
* The library is free for all purposes without any express
* guarantee it works.
*
* Tom St Denis, tomstdenis@iahu.ca, http://math.libtomcrypt.org
*/
#include <tommath.h>

/* Greatest Common Divisor using the binary method [Algorithm B, page 338, vol2 of TAOCP]
*/
int
mp_gcd (mp_int * a, mp_int * b, mp_int * c)
{
  mp_int u, v, t;
  int k, res, neg;

  /* either zero than gcd is the largest */
  if (mp_iszero (a) == 1 && mp_iszero (b) == 0) {
    return mp_copy (b, c);
  }
  if (mp_iszero (a) == 0 && mp_iszero (b) == 1) {
    return mp_copy (a, c);
  }
  if (mp_iszero (a) == 1 && mp_iszero (b) == 1) {
    mp_set (c, 1);
    return MP_OKAY;
  }

  /* if both are negative they share (-1) as a common divisor */
  neg = (a->sign == b->sign) ? a->sign : MP_ZPOS;

  if ((res = mp_init_copy (&u, a)) != MP_OKAY) {
    return res;
  }

  if ((res = mp_init_copy (&v, b)) != MP_OKAY) {
    goto __U;
  }

  /* must be positive for the remainder of the algorithm */
  u.sign = v.sign = MP_ZPOS;

  if ((res = mp_init (&t)) != MP_OKAY) {
    goto __V;
  }

  /* B1. Find power of two */
  k = 0;
  while (mp_iseven(&u) == 1 && mp_iseven(&v) == 1) {
    ++k;
    if ((res = mp_div_2 (&u, &u)) != MP_OKAY) {
      goto __T;
    }
    if ((res = mp_div_2 (&v, &v)) != MP_OKAY) {
      goto __T;
    }
  }

  /* B2. Initialize */
  if (mp_isodd(&u) == 1) {
    /* t = -v */
    if ((res = mp_copy (&v, &t)) != MP_OKAY) {
      goto __T;
    }
    t.sign = MP_NEG;
  } else {
    /* t = u */
    if ((res = mp_copy (&u, &t)) != MP_OKAY) {
      goto __T;
    }
  }

  do {
    /* B3 (and B4). Halve t, if even */
    while (t.used != 0 && mp_iseven(&t) == 1) {
      if ((res = mp_div_2 (&t, &t)) != MP_OKAY) {
        goto __T;
      }
    }

    /* B5. if t>0 then u=t otherwise v=-t */
    if (t.used != 0 && t.sign != MP_NEG) {
      if ((res = mp_copy (&t, &u)) != MP_OKAY) {
        goto __T;
      }
    } else {
      if ((res = mp_copy (&t, &v)) != MP_OKAY) {
        goto __T;
      }
      v.sign = (v.sign == MP_ZPOS) ? MP_NEG : MP_ZPOS;
    }

    /* B6. t = u - v, if t != 0 loop otherwise terminate */
    if ((res = mp_sub (&u, &v, &t)) != MP_OKAY) {
      goto __T;
    }
  } while (mp_iszero(&t) == 0);

  /* multiply by 2^k which we divided out at the beginning */
  if ((res = mp_mul_2d (&u, k, &u)) != MP_OKAY) {
    goto __T;
  }

  mp_exch (&u, c);
  c->sign = neg;
  res = MP_OKAY;
__T:mp_clear (&t);
__V:mp_clear (&u);
__U:mp_clear (&v);
  return res;
}
Something went wrong with that request. Please try again.