1 /* mpz_gcd -- Calculate the greatest common divisior of two integers.
3 Copyright (C) 1991 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 General Public License as published by
9 the Free Software Foundation; either version 2, or (at your option)
12 The GNU MP Library is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
17 You should have received a copy of the GNU General Public License
18 along with the GNU MP Library; see the file COPYING. If not, write to
19 the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA. */
28 mpz_gcd (MP_INT *w, const MP_INT *u, const MP_INT *v)
35 #else /* BERKELEY_MP */
38 gcd (const MP_INT *u, const MP_INT *v, MP_INT *w)
45 #endif /* BERKELEY_MP */
47 mp_size usize, vsize, wsize;
57 usize = ABS (u->size);
58 vsize = ABS (v->size);
64 _mpz_realloc (w, vsize);
67 MPN_COPY (w->d, v->d, vsize);
75 _mpz_realloc (w, usize);
78 MPN_COPY (w->d, u->d, usize);
82 /* Make U odd by shifting it down as many bit positions as there
83 are zero bits. Put the result in temporary space. */
84 up = (mp_ptr) alloca (usize * BYTES_PER_MP_LIMB);
86 for (i = 0; (d = up_in[i]) == 0; i++)
88 count_leading_zeros (bcnt, d & -d);
89 bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
90 usize = mpn_rshift (up, up_in + i, usize - i, bcnt);
92 bcnt += i * BITS_PER_MP_LIMB;
95 /* Make V odd by shifting it down as many bit positions as there
96 are zero bits. Put the result in temporary space. */
97 vp = (mp_ptr) alloca (vsize * BYTES_PER_MP_LIMB);
99 for (i = 0; (d = vp_in[i]) == 0; i++)
101 count_leading_zeros (bcnt, d & -d);
102 bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
103 vsize = mpn_rshift (vp, vp_in + i, vsize - i, bcnt);
105 /* W_BCNT is set to the minimum of the number of zero bits in U and V.
106 Thus it represents the number of common 2 factors. */
107 bcnt += i * BITS_PER_MP_LIMB;
115 cmp = usize - vsize != 0 ? usize - vsize : mpn_cmp (up, vp, usize);
117 /* If U and V have become equal, we have found the GCD. */
123 /* Replace U by (U - V) >> cnt, with cnt being the least value
124 making U odd again. */
126 usize += mpn_sub (up, up, usize, vp, vsize);
127 for (i = 0; (d = up[i]) == 0; i++)
129 count_leading_zeros (bcnt, d & -d);
130 bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
131 usize = mpn_rshift (up, up + i, usize - i, bcnt);
135 /* Replace V by (V - U) >> cnt, with cnt being the least value
136 making V odd again. */
138 vsize += mpn_sub (vp, vp, vsize, up, usize);
139 for (i = 0; (d = vp[i]) == 0; i++)
141 count_leading_zeros (bcnt, d & -d);
142 bcnt = BITS_PER_MP_LIMB - 1 - bcnt;
143 vsize = mpn_rshift (vp, vp + i, vsize - i, bcnt);
147 /* GCD(U_IN, V_IN) now is U * 2**W_BCNT. */
149 wsize = usize + w_bcnt / BITS_PER_MP_LIMB + 1;
150 if (w->alloc < wsize)
151 _mpz_realloc (w, wsize);
155 MPN_ZERO (wp, w_bcnt / BITS_PER_MP_LIMB);
157 cy_digit = mpn_lshift (wp + w_bcnt / BITS_PER_MP_LIMB, up, usize,
158 w_bcnt % BITS_PER_MP_LIMB);
159 wsize = usize + w_bcnt / BITS_PER_MP_LIMB;
162 wp[wsize] = cy_digit;