1 /* mpz/gcd.c: Calculate the greatest common divisor of two integers.
3 Copyright 1991, 1993, 1994, 1996, 2000, 2001, 2002, 2005 Free Software
6 This file is part of the GNU MP Library.
8 The GNU MP Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
13 The GNU MP Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
31 mpz_gcd (mpz_ptr g, mpz_srcptr u, mpz_srcptr v)
32 #else /* BERKELEY_MP */
33 gcd (mpz_srcptr u, mpz_srcptr v, mpz_ptr g)
34 #endif /* BERKELEY_MP */
36 unsigned long int g_zero_bits, u_zero_bits, v_zero_bits;
37 mp_size_t g_zero_limbs, u_zero_limbs, v_zero_limbs;
40 mp_size_t usize = ABS (u->_mp_size);
42 mp_size_t vsize = ABS (v->_mp_size);
52 if (g->_mp_alloc < vsize)
53 _mpz_realloc (g, vsize);
54 MPN_COPY (g->_mp_d, vp, vsize);
64 if (g->_mp_alloc < usize)
65 _mpz_realloc (g, usize);
66 MPN_COPY (g->_mp_d, up, usize);
73 g->_mp_d[0] = mpn_gcd_1 (vp, vsize, up[0]);
80 g->_mp_d[0] = mpn_gcd_1 (up, usize, vp[0]);
86 /* Eliminate low zero bits from U and V and move to temporary storage. */
89 u_zero_limbs = up - u->_mp_d;
90 usize -= u_zero_limbs;
91 count_trailing_zeros (u_zero_bits, *up);
93 up = TMP_ALLOC_LIMBS (usize);
96 mpn_rshift (up, tp, usize, u_zero_bits);
97 usize -= up[usize - 1] == 0;
100 MPN_COPY (up, tp, usize);
104 v_zero_limbs = vp - v->_mp_d;
105 vsize -= v_zero_limbs;
106 count_trailing_zeros (v_zero_bits, *vp);
108 vp = TMP_ALLOC_LIMBS (vsize);
109 if (v_zero_bits != 0)
111 mpn_rshift (vp, tp, vsize, v_zero_bits);
112 vsize -= vp[vsize - 1] == 0;
115 MPN_COPY (vp, tp, vsize);
117 if (u_zero_limbs > v_zero_limbs)
119 g_zero_limbs = v_zero_limbs;
120 g_zero_bits = v_zero_bits;
122 else if (u_zero_limbs < v_zero_limbs)
124 g_zero_limbs = u_zero_limbs;
125 g_zero_bits = u_zero_bits;
129 g_zero_limbs = u_zero_limbs;
130 g_zero_bits = MIN (u_zero_bits, v_zero_bits);
133 /* Call mpn_gcd. The 2nd argument must not have more bits than the 1st. */
134 vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1]))
135 ? mpn_gcd (vp, vp, vsize, up, usize)
136 : mpn_gcd (vp, up, usize, vp, vsize);
138 /* Here G <-- V << (g_zero_limbs*GMP_LIMB_BITS + g_zero_bits). */
139 gsize = vsize + g_zero_limbs;
140 if (g_zero_bits != 0)
143 gsize += (vp[vsize - 1] >> (GMP_NUMB_BITS - g_zero_bits)) != 0;
144 if (g->_mp_alloc < gsize)
145 _mpz_realloc (g, gsize);
146 MPN_ZERO (g->_mp_d, g_zero_limbs);
148 tp = g->_mp_d + g_zero_limbs;
149 cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
155 if (g->_mp_alloc < gsize)
156 _mpz_realloc (g, gsize);
157 MPN_ZERO (g->_mp_d, g_zero_limbs);
158 MPN_COPY (g->_mp_d + g_zero_limbs, vp, vsize);