Annotation of OpenXM_contrib/gmp/mpz/gcd.c, Revision 1.1.1.3
1.1 maekawa 1: /* mpz/gcd.c: Calculate the greatest common divisor of two integers.
2:
1.1.1.3 ! ohara 3: Copyright 1991, 1993, 1994, 1996, 2000, 2001, 2002 Free Software Foundation,
! 4: Inc.
1.1 maekawa 5:
6: This file is part of the GNU MP Library.
7:
8: The GNU MP Library is free software; you can redistribute it and/or modify
1.1.1.2 maekawa 9: it under the terms of the GNU Lesser General Public License as published by
10: the Free Software Foundation; either version 2.1 of the License, or (at your
1.1 maekawa 11: option) any later version.
12:
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
1.1.1.2 maekawa 15: or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
1.1 maekawa 16: License for more details.
17:
1.1.1.2 maekawa 18: You should have received a copy of the GNU Lesser General Public License
1.1 maekawa 19: along with the GNU MP Library; see the file COPYING.LIB. If not, write to
20: the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
21: MA 02111-1307, USA. */
22:
23: #include "gmp.h"
24: #include "gmp-impl.h"
25: #include "longlong.h"
1.1.1.2 maekawa 26: #ifdef BERKELEY_MP
27: #include "mp.h"
28: #endif
1.1 maekawa 29:
30:
31: void
1.1.1.3 ! ohara 32: #ifndef BERKELEY_MP
1.1 maekawa 33: mpz_gcd (mpz_ptr g, mpz_srcptr u, mpz_srcptr v)
34: #else /* BERKELEY_MP */
35: gcd (mpz_srcptr u, mpz_srcptr v, mpz_ptr g)
36: #endif /* BERKELEY_MP */
37: {
38: unsigned long int g_zero_bits, u_zero_bits, v_zero_bits;
39: mp_size_t g_zero_limbs, u_zero_limbs, v_zero_limbs;
40: mp_ptr tp;
41: mp_ptr up = u->_mp_d;
42: mp_size_t usize = ABS (u->_mp_size);
43: mp_ptr vp = v->_mp_d;
44: mp_size_t vsize = ABS (v->_mp_size);
45: mp_size_t gsize;
46: TMP_DECL (marker);
47:
48: /* GCD(0, V) == V. */
49: if (usize == 0)
50: {
51: g->_mp_size = vsize;
52: if (g == v)
53: return;
54: if (g->_mp_alloc < vsize)
55: _mpz_realloc (g, vsize);
56: MPN_COPY (g->_mp_d, vp, vsize);
57: return;
58: }
59:
60: /* GCD(U, 0) == U. */
61: if (vsize == 0)
62: {
63: g->_mp_size = usize;
64: if (g == u)
65: return;
66: if (g->_mp_alloc < usize)
67: _mpz_realloc (g, usize);
68: MPN_COPY (g->_mp_d, up, usize);
69: return;
70: }
71:
72: if (usize == 1)
73: {
74: g->_mp_size = 1;
75: g->_mp_d[0] = mpn_gcd_1 (vp, vsize, up[0]);
76: return;
77: }
78:
79: if (vsize == 1)
80: {
81: g->_mp_size = 1;
82: g->_mp_d[0] = mpn_gcd_1 (up, usize, vp[0]);
83: return;
84: }
85:
86: TMP_MARK (marker);
87:
88: /* Eliminate low zero bits from U and V and move to temporary storage. */
89: while (*up == 0)
90: up++;
91: u_zero_limbs = up - u->_mp_d;
92: usize -= u_zero_limbs;
93: count_trailing_zeros (u_zero_bits, *up);
94: tp = up;
95: up = (mp_ptr) TMP_ALLOC (usize * BYTES_PER_MP_LIMB);
96: if (u_zero_bits != 0)
97: {
98: mpn_rshift (up, tp, usize, u_zero_bits);
99: usize -= up[usize - 1] == 0;
100: }
101: else
102: MPN_COPY (up, tp, usize);
103:
104: while (*vp == 0)
105: vp++;
106: v_zero_limbs = vp - v->_mp_d;
107: vsize -= v_zero_limbs;
108: count_trailing_zeros (v_zero_bits, *vp);
109: tp = vp;
110: vp = (mp_ptr) TMP_ALLOC (vsize * BYTES_PER_MP_LIMB);
111: if (v_zero_bits != 0)
112: {
113: mpn_rshift (vp, tp, vsize, v_zero_bits);
114: vsize -= vp[vsize - 1] == 0;
115: }
116: else
117: MPN_COPY (vp, tp, vsize);
118:
119: if (u_zero_limbs > v_zero_limbs)
120: {
121: g_zero_limbs = v_zero_limbs;
122: g_zero_bits = v_zero_bits;
123: }
124: else if (u_zero_limbs < v_zero_limbs)
125: {
126: g_zero_limbs = u_zero_limbs;
127: g_zero_bits = u_zero_bits;
128: }
129: else /* Equal. */
130: {
131: g_zero_limbs = u_zero_limbs;
132: g_zero_bits = MIN (u_zero_bits, v_zero_bits);
133: }
134:
1.1.1.2 maekawa 135: /* Call mpn_gcd. The 2nd argument must not have more bits than the 1st. */
1.1 maekawa 136: vsize = (usize < vsize || (usize == vsize && up[usize-1] < vp[vsize-1]))
1.1.1.2 maekawa 137: ? mpn_gcd (vp, vp, vsize, up, usize)
138: : mpn_gcd (vp, up, usize, vp, vsize);
1.1 maekawa 139:
140: /* Here G <-- V << (g_zero_limbs*BITS_PER_MP_LIMB + g_zero_bits). */
141: gsize = vsize + g_zero_limbs;
142: if (g_zero_bits != 0)
143: {
144: mp_limb_t cy_limb;
1.1.1.3 ! ohara 145: gsize += (vp[vsize - 1] >> (GMP_NUMB_BITS - g_zero_bits)) != 0;
1.1 maekawa 146: if (g->_mp_alloc < gsize)
147: _mpz_realloc (g, gsize);
148: MPN_ZERO (g->_mp_d, g_zero_limbs);
149:
150: tp = g->_mp_d + g_zero_limbs;
151: cy_limb = mpn_lshift (tp, vp, vsize, g_zero_bits);
152: if (cy_limb != 0)
153: tp[vsize] = cy_limb;
154: }
155: else
156: {
157: if (g->_mp_alloc < gsize)
158: _mpz_realloc (g, gsize);
159: MPN_ZERO (g->_mp_d, g_zero_limbs);
160: MPN_COPY (g->_mp_d + g_zero_limbs, vp, vsize);
161: }
162:
163: g->_mp_size = gsize;
164: TMP_FREE (marker);
165: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to gcd.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / gmp / mpz