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