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 */
mp_gcd(mp_int * a,mp_int * b,mp_int * c)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 */
108