1 /* mpz/gcd.c: Calculate the greatest common divisor of two integers.
3 Copyright (C) 1991, 1993, 1994, 1996, 2000 Free Software Foundation, Inc.
5 This file is part of the GNU MP Library.
7 The GNU MP Library is free software; you can redistribute it and/or modify
8 it under the terms of the GNU Lesser General Public License as published by
9 the Free Software Foundation; either version 2.1 of the License, or (at your
10 option) any later version.
12 The GNU MP Library is distributed in the hope that it will be useful, but
13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15 License for more details.
17 You should have received a copy of the GNU Lesser General Public License
18 along with the GNU MP Library; see the file COPYING.LIB. If not, write to
19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
20 MA 02111-1307, USA. */
33 mpz_gcd (mpz_ptr g, mpz_srcptr u, mpz_srcptr v)
40 #else /* BERKELEY_MP */
43 gcd (mpz_srcptr u, mpz_srcptr v, mpz_ptr g)
50 #endif /* BERKELEY_MP */
53 unsigned long int g_zero_bits, u_zero_bits, v_zero_bits;
54 mp_size_t g_zero_limbs, u_zero_limbs, v_zero_limbs;
57 mp_size_t usize = ABS (u->_mp_size);
59 mp_size_t vsize = ABS (v->_mp_size);
69 if (g->_mp_alloc < vsize)
70 _mpz_realloc (g, vsize);
71 MPN_COPY (g->_mp_d, vp, vsize);
81 if (g->_mp_alloc < usize)
82 _mpz_realloc (g, usize);
83 MPN_COPY (g->_mp_d, up, usize);
90 g->_mp_d[0] = mpn_gcd_1 (vp, vsize, up[0]);
97 g->_mp_d[0] = mpn_gcd_1 (up, usize, vp[0]);
103 /* Eliminate low zero bits from U and V and move to temporary storage. */
106 u_zero_limbs = up - u->_mp_d;
107 usize -= u_zero_limbs;
108 count_trailing_zeros (u_zero_bits, *up);
110 up = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
111 if (u_zero_bits != 0)
113 mpn_rshift (up, tp, usize, u_zero_bits);
114 usize -= up[usize - 1] == 0;
117 MPN_COPY (up, tp, usize);
121 v_zero_limbs = vp - v->_mp_d;
122 vsize -= v_zero_limbs;
123 count_trailing_zeros (v_zero_bits, *vp);
125 vp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
126 if (v_zero_bits != 0)
128 mpn_rshift (vp, tp, vsize, v_zero_bits);
129 vsize -= vp[vsize - 1] == 0;
132 MPN_COPY (vp, tp, vsize);
134 if (u_zero_limbs > v_zero_limbs)
136 g_zero_limbs = v_zero_limbs;
137 g_zero_bits = v_zero_bits;
139 else if (u_zero_limbs < v_zero_limbs)
141 g_zero_limbs = u_zero_limbs;
142 g_zero_bits = u_zero_bits;
146 g_zero_limbs = u_zero_limbs;
147 g_zero_bits = MIN (u_zero_bits, v_zero_bits);
150 /* Call mpn_gcd. The 2nd argument must not have more bits than the 1st. */
151 vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1]))
152 ? mpn_gcd (vp, vp, vsize, up, usize)
153 : mpn_gcd (vp, up, usize, vp, vsize);
155 /* Here G <-- V << (g_zero_limbs*BITS_PER_MP_LIMB + g_zero_bits). */
156 gsize = vsize + g_zero_limbs;
157 if (g_zero_bits != 0)
160 gsize += (vp[vsize - 1] >> (BITS_PER_MP_LIMB - g_zero_bits)) != 0;
161 if (g->_mp_alloc < gsize)
162 _mpz_realloc (g, gsize);
163 MPN_ZERO (g->_mp_d, g_zero_limbs);
165 tp = g->_mp_d + g_zero_limbs;
166 cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
172 if (g->_mp_alloc < gsize)
173 _mpz_realloc (g, gsize);
174 MPN_ZERO (g->_mp_d, g_zero_limbs);
175 MPN_COPY (g->_mp_d + g_zero_limbs, vp, vsize);