1/* mpz/gcd.c:   Calculate the greatest common divisor of two integers.
2
3Copyright 1991, 1993, 1994, 1996, 2000, 2001, 2002, 2005 Free Software
4Foundation, Inc.
5
6This file is part of the GNU MP Library.
7
8The GNU MP Library is free software; you can redistribute it and/or modify
9it under the terms of the GNU Lesser General Public License as published by
10the Free Software Foundation; either version 3 of the License, or (at your
11option) any later version.
12
13The GNU MP Library is distributed in the hope that it will be useful, but
14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
16License for more details.
17
18You should have received a copy of the GNU Lesser General Public License
19along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
20
21#include "gmp.h"
22#include "gmp-impl.h"
23#include "longlong.h"
24#ifdef BERKELEY_MP
25#include "mp.h"
26#endif
27
28
29void
30#ifndef BERKELEY_MP
31mpz_gcd (mpz_ptr g, mpz_srcptr u, mpz_srcptr v)
32#else /* BERKELEY_MP */
33gcd (mpz_srcptr u, mpz_srcptr v, mpz_ptr g)
34#endif /* BERKELEY_MP */
35{
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;
38  mp_ptr tp;
39  mp_ptr up = u->_mp_d;
40  mp_size_t usize = ABS (u->_mp_size);
41  mp_ptr vp = v->_mp_d;
42  mp_size_t vsize = ABS (v->_mp_size);
43  mp_size_t gsize;
44  TMP_DECL;
45
46  /* GCD(0, V) == V.  */
47  if (usize == 0)
48    {
49      g->_mp_size = vsize;
50      if (g == v)
51	return;
52      if (g->_mp_alloc < vsize)
53	_mpz_realloc (g, vsize);
54      MPN_COPY (g->_mp_d, vp, vsize);
55      return;
56    }
57
58  /* GCD(U, 0) == U.  */
59  if (vsize == 0)
60    {
61      g->_mp_size = usize;
62      if (g == u)
63	return;
64      if (g->_mp_alloc < usize)
65	_mpz_realloc (g, usize);
66      MPN_COPY (g->_mp_d, up, usize);
67      return;
68    }
69
70  if (usize == 1)
71    {
72      g->_mp_size = 1;
73      g->_mp_d[0] = mpn_gcd_1 (vp, vsize, up[0]);
74      return;
75    }
76
77  if (vsize == 1)
78    {
79      g->_mp_size = 1;
80      g->_mp_d[0] = mpn_gcd_1 (up, usize, vp[0]);
81      return;
82    }
83
84  TMP_MARK;
85
86  /*  Eliminate low zero bits from U and V and move to temporary storage.  */
87  while (*up == 0)
88    up++;
89  u_zero_limbs = up - u->_mp_d;
90  usize -= u_zero_limbs;
91  count_trailing_zeros (u_zero_bits, *up);
92  tp = up;
93  up = TMP_ALLOC_LIMBS (usize);
94  if (u_zero_bits != 0)
95    {
96      mpn_rshift (up, tp, usize, u_zero_bits);
97      usize -= up[usize - 1] == 0;
98    }
99  else
100    MPN_COPY (up, tp, usize);
101
102  while (*vp == 0)
103    vp++;
104  v_zero_limbs = vp - v->_mp_d;
105  vsize -= v_zero_limbs;
106  count_trailing_zeros (v_zero_bits, *vp);
107  tp = vp;
108  vp = TMP_ALLOC_LIMBS (vsize);
109  if (v_zero_bits != 0)
110    {
111      mpn_rshift (vp, tp, vsize, v_zero_bits);
112      vsize -= vp[vsize - 1] == 0;
113    }
114  else
115    MPN_COPY (vp, tp, vsize);
116
117  if (u_zero_limbs > v_zero_limbs)
118    {
119      g_zero_limbs = v_zero_limbs;
120      g_zero_bits = v_zero_bits;
121    }
122  else if (u_zero_limbs < v_zero_limbs)
123    {
124      g_zero_limbs = u_zero_limbs;
125      g_zero_bits = u_zero_bits;
126    }
127  else  /*  Equal.  */
128    {
129      g_zero_limbs = u_zero_limbs;
130      g_zero_bits = MIN (u_zero_bits, v_zero_bits);
131    }
132
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);
137
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)
141    {
142      mp_limb_t cy_limb;
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);
147
148      tp = g->_mp_d + g_zero_limbs;
149      cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
150      if (cy_limb != 0)
151	tp[vsize] = cy_limb;
152    }
153  else
154    {
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);
159    }
160
161  g->_mp_size = gsize;
162  TMP_FREE;
163}
164