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@gmail.com, http://libtom.org
18 */
19
20/* Greatest Common Divisor using the binary method */
21int 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;
99LBL_V:mp_clear (&u);
100LBL_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