Annotation of OpenXM_contrib/gmp/mpn/generic/gcd.c, Revision 1.1.1.1
1.1 maekawa 1: /* mpn/gcd.c: mpn_gcd for gcd of two odd integers.
2:
3: Copyright (C) 1991, 1993, 1994, 1995, 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: /* Integer greatest common divisor of two unsigned integers, using
23: the accelerated algorithm (see reference below).
24:
25: mp_size_t mpn_gcd (vp, vsize, up, usize).
26:
27: Preconditions [U = (up, usize) and V = (vp, vsize)]:
28:
29: 1. V is odd.
30: 2. numbits(U) >= numbits(V).
31:
32: Both U and V are destroyed by the operation. The result is left at vp,
33: and its size is returned.
34:
35: Ken Weber (kweber@mat.ufrgs.br, kweber@mcs.kent.edu)
36:
37: Funding for this work has been partially provided by Conselho Nacional
38: de Desenvolvimento Cienti'fico e Tecnolo'gico (CNPq) do Brazil, Grant
39: 301314194-2, and was done while I was a visiting reseacher in the Instituto
40: de Matema'tica at Universidade Federal do Rio Grande do Sul (UFRGS).
41:
42: Refer to
43: K. Weber, The accelerated integer GCD algorithm, ACM Transactions on
44: Mathematical Software, v. 21 (March), 1995, pp. 111-122. */
45:
46: #include "gmp.h"
47: #include "gmp-impl.h"
48: #include "longlong.h"
49:
50: /* If MIN (usize, vsize) > ACCEL_THRESHOLD, then the accelerated algorithm is
51: used, otherwise the binary algorithm is used. This may be adjusted for
52: different architectures. */
53: #ifndef ACCEL_THRESHOLD
54: #define ACCEL_THRESHOLD 4
55: #endif
56:
57: /* When U and V differ in size by more than BMOD_THRESHOLD, the accelerated
58: algorithm reduces using the bmod operation. Otherwise, the k-ary reduction
59: is used. 0 <= BMOD_THRESHOLD < BITS_PER_MP_LIMB. */
60: enum
61: {
62: BMOD_THRESHOLD = BITS_PER_MP_LIMB/2
63: };
64:
65: #define SIGN_BIT (~(~(mp_limb_t)0 >> 1))
66:
67:
68: #define SWAP_LIMB(UL, VL) do{mp_limb_t __l=(UL);(UL)=(VL);(VL)=__l;}while(0)
69: #define SWAP_PTR(UP, VP) do{mp_ptr __p=(UP);(UP)=(VP);(VP)=__p;}while(0)
70: #define SWAP_SZ(US, VS) do{mp_size_t __s=(US);(US)=(VS);(VS)=__s;}while(0)
71: #define SWAP_MPN(UP, US, VP, VS) do{SWAP_PTR(UP,VP);SWAP_SZ(US,VS);}while(0)
72:
73: /* Use binary algorithm to compute V <-- GCD (V, U) for usize, vsize == 2.
74: Both U and V must be odd. */
75: static __gmp_inline mp_size_t
76: #if __STDC__
77: gcd_2 (mp_ptr vp, mp_srcptr up)
78: #else
79: gcd_2 (vp, up)
80: mp_ptr vp;
81: mp_srcptr up;
82: #endif
83: {
84: mp_limb_t u0, u1, v0, v1;
85: mp_size_t vsize;
86:
87: u0 = up[0], u1 = up[1], v0 = vp[0], v1 = vp[1];
88:
89: while (u1 != v1 && u0 != v0)
90: {
91: unsigned long int r;
92: if (u1 > v1)
93: {
94: u1 -= v1 + (u0 < v0), u0 -= v0;
95: count_trailing_zeros (r, u0);
96: u0 = u1 << (BITS_PER_MP_LIMB - r) | u0 >> r;
97: u1 >>= r;
98: }
99: else /* u1 < v1. */
100: {
101: v1 -= u1 + (v0 < u0), v0 -= u0;
102: count_trailing_zeros (r, v0);
103: v0 = v1 << (BITS_PER_MP_LIMB - r) | v0 >> r;
104: v1 >>= r;
105: }
106: }
107:
108: vp[0] = v0, vp[1] = v1, vsize = 1 + (v1 != 0);
109:
110: /* If U == V == GCD, done. Otherwise, compute GCD (V, |U - V|). */
111: if (u1 == v1 && u0 == v0)
112: return vsize;
113:
114: v0 = (u0 == v0) ? (u1 > v1) ? u1-v1 : v1-u1 : (u0 > v0) ? u0-v0 : v0-u0;
115: vp[0] = mpn_gcd_1 (vp, vsize, v0);
116:
117: return 1;
118: }
119:
120: /* The function find_a finds 0 < N < 2^BITS_PER_MP_LIMB such that there exists
121: 0 < |D| < 2^BITS_PER_MP_LIMB, and N == D * C mod 2^(2*BITS_PER_MP_LIMB).
122: In the reference article, D was computed along with N, but it is better to
123: compute D separately as D <-- N / C mod 2^(BITS_PER_MP_LIMB + 1), treating
124: the result as a twos' complement signed integer.
125:
126: Initialize N1 to C mod 2^(2*BITS_PER_MP_LIMB). According to the reference
127: article, N2 should be initialized to 2^(2*BITS_PER_MP_LIMB), but we use
128: 2^(2*BITS_PER_MP_LIMB) - N1 to start the calculations within double
129: precision. If N2 > N1 initially, the first iteration of the while loop
130: will swap them. In all other situations, N1 >= N2 is maintained. */
131:
132: static __gmp_inline mp_limb_t
133: #if __STDC__
134: find_a (mp_srcptr cp)
135: #else
136: find_a (cp)
137: mp_srcptr cp;
138: #endif
139: {
140: unsigned long int leading_zero_bits = 0;
141:
142: mp_limb_t n1_l = cp[0]; /* N1 == n1_h * 2^BITS_PER_MP_LIMB + n1_l. */
143: mp_limb_t n1_h = cp[1];
144:
145: mp_limb_t n2_l = -n1_l; /* N2 == n2_h * 2^BITS_PER_MP_LIMB + n2_l. */
146: mp_limb_t n2_h = ~n1_h;
147:
148: /* Main loop. */
149: while (n2_h) /* While N2 >= 2^BITS_PER_MP_LIMB. */
150: {
151: /* N1 <-- N1 % N2. */
152: if ((SIGN_BIT >> leading_zero_bits & n2_h) == 0)
153: {
154: unsigned long int i;
155: count_leading_zeros (i, n2_h);
156: i -= leading_zero_bits, leading_zero_bits += i;
157: n2_h = n2_h<<i | n2_l>>(BITS_PER_MP_LIMB - i), n2_l <<= i;
158: do
159: {
160: if (n1_h > n2_h || (n1_h == n2_h && n1_l >= n2_l))
161: n1_h -= n2_h + (n1_l < n2_l), n1_l -= n2_l;
162: n2_l = n2_l>>1 | n2_h<<(BITS_PER_MP_LIMB - 1), n2_h >>= 1;
163: i -= 1;
164: }
165: while (i);
166: }
167: if (n1_h > n2_h || (n1_h == n2_h && n1_l >= n2_l))
168: n1_h -= n2_h + (n1_l < n2_l), n1_l -= n2_l;
169:
170: SWAP_LIMB (n1_h, n2_h);
171: SWAP_LIMB (n1_l, n2_l);
172: }
173:
174: return n2_l;
175: }
176:
177: mp_size_t
178: #if __STDC__
179: mpn_gcd (mp_ptr gp, mp_ptr vp, mp_size_t vsize, mp_ptr up, mp_size_t usize)
180: #else
181: mpn_gcd (gp, vp, vsize, up, usize)
182: mp_ptr gp;
183: mp_ptr vp;
184: mp_size_t vsize;
185: mp_ptr up;
186: mp_size_t usize;
187: #endif
188: {
189: mp_ptr orig_vp = vp;
190: mp_size_t orig_vsize = vsize;
191: int binary_gcd_ctr; /* Number of times binary gcd will execute. */
192: TMP_DECL (marker);
193:
194: TMP_MARK (marker);
195:
196: /* Use accelerated algorithm if vsize is over ACCEL_THRESHOLD.
197: Two EXTRA limbs for U and V are required for kary reduction. */
198: if (vsize > ACCEL_THRESHOLD)
199: {
200: unsigned long int vbitsize, d;
201: mp_ptr orig_up = up;
202: mp_size_t orig_usize = usize;
203: mp_ptr anchor_up = (mp_ptr) TMP_ALLOC ((usize + 2) * BYTES_PER_MP_LIMB);
204:
205: MPN_COPY (anchor_up, orig_up, usize);
206: up = anchor_up;
207:
208: count_leading_zeros (d, up[usize-1]);
209: d = usize * BITS_PER_MP_LIMB - d;
210: count_leading_zeros (vbitsize, vp[vsize-1]);
211: vbitsize = vsize * BITS_PER_MP_LIMB - vbitsize;
212: d = d - vbitsize + 1;
213:
214: /* Use bmod reduction to quickly discover whether V divides U. */
215: up[usize++] = 0; /* Insert leading zero. */
216: mpn_bdivmod (up, up, usize, vp, vsize, d);
217:
218: /* Now skip U/V mod 2^d and any low zero limbs. */
219: d /= BITS_PER_MP_LIMB, up += d, usize -= d;
220: while (usize != 0 && up[0] == 0)
221: up++, usize--;
222:
223: if (usize == 0) /* GCD == ORIG_V. */
224: goto done;
225:
226: vp = (mp_ptr) TMP_ALLOC ((vsize + 2) * BYTES_PER_MP_LIMB);
227: MPN_COPY (vp, orig_vp, vsize);
228:
229: do /* Main loop. */
230: {
231: if (up[usize-1] & SIGN_BIT) /* U < 0; take twos' compl. */
232: {
233: mp_size_t i;
234: anchor_up[0] = -up[0];
235: for (i = 1; i < usize; i++)
236: anchor_up[i] = ~up[i];
237: up = anchor_up;
238: }
239:
240: MPN_NORMALIZE_NOT_ZERO (up, usize);
241:
242: if ((up[0] & 1) == 0) /* Result even; remove twos. */
243: {
244: unsigned long int r;
245: count_trailing_zeros (r, up[0]);
246: mpn_rshift (anchor_up, up, usize, r);
247: usize -= (anchor_up[usize-1] == 0);
248: }
249: else if (anchor_up != up)
250: MPN_COPY (anchor_up, up, usize);
251:
252: SWAP_MPN (anchor_up, usize, vp, vsize);
253: up = anchor_up;
254:
255: if (vsize <= 2) /* Kary can't handle < 2 limbs and */
256: break; /* isn't efficient for == 2 limbs. */
257:
258: d = vbitsize;
259: count_leading_zeros (vbitsize, vp[vsize-1]);
260: vbitsize = vsize * BITS_PER_MP_LIMB - vbitsize;
261: d = d - vbitsize + 1;
262:
263: if (d > BMOD_THRESHOLD) /* Bmod reduction. */
264: {
265: up[usize++] = 0;
266: mpn_bdivmod (up, up, usize, vp, vsize, d);
267: d /= BITS_PER_MP_LIMB, up += d, usize -= d;
268: }
269: else /* Kary reduction. */
270: {
271: mp_limb_t bp[2], cp[2];
272:
273: /* C <-- V/U mod 2^(2*BITS_PER_MP_LIMB). */
274: cp[0] = vp[0], cp[1] = vp[1];
275: mpn_bdivmod (cp, cp, 2, up, 2, 2*BITS_PER_MP_LIMB);
276:
277: /* U <-- find_a (C) * U. */
278: up[usize] = mpn_mul_1 (up, up, usize, find_a (cp));
279: usize++;
280:
281: /* B <-- A/C == U/V mod 2^(BITS_PER_MP_LIMB + 1).
282: bp[0] <-- U/V mod 2^BITS_PER_MP_LIMB and
283: bp[1] <-- ( (U - bp[0] * V)/2^BITS_PER_MP_LIMB ) / V mod 2 */
284: bp[0] = up[0], bp[1] = up[1];
285: mpn_bdivmod (bp, bp, 2, vp, 2, BITS_PER_MP_LIMB);
286: bp[1] &= 1; /* Since V is odd, division is unnecessary. */
287:
288: up[usize++] = 0;
289: if (bp[1]) /* B < 0: U <-- U + (-B) * V. */
290: {
291: mp_limb_t c = mpn_addmul_1 (up, vp, vsize, -bp[0]);
292: mpn_add_1 (up + vsize, up + vsize, usize - vsize, c);
293: }
294: else /* B >= 0: U <-- U - B * V. */
295: {
296: mp_limb_t b = mpn_submul_1 (up, vp, vsize, bp[0]);
297: mpn_sub_1 (up + vsize, up + vsize, usize - vsize, b);
298: }
299:
300: up += 2, usize -= 2; /* At least two low limbs are zero. */
301: }
302:
303: /* Must remove low zero limbs before complementing. */
304: while (usize != 0 && up[0] == 0)
305: up++, usize--;
306: }
307: while (usize);
308:
309: /* Compute GCD (ORIG_V, GCD (ORIG_U, V)). Binary will execute twice. */
310: up = orig_up, usize = orig_usize;
311: binary_gcd_ctr = 2;
312: }
313: else
314: binary_gcd_ctr = 1;
315:
316: /* Finish up with the binary algorithm. Executes once or twice. */
317: for ( ; binary_gcd_ctr--; up = orig_vp, usize = orig_vsize)
318: {
319: if (usize > 2) /* First make U close to V in size. */
320: {
321: unsigned long int vbitsize, d;
322: count_leading_zeros (d, up[usize-1]);
323: d = usize * BITS_PER_MP_LIMB - d;
324: count_leading_zeros (vbitsize, vp[vsize-1]);
325: vbitsize = vsize * BITS_PER_MP_LIMB - vbitsize;
326: d = d - vbitsize - 1;
327: if (d != -(unsigned long int)1 && d > 2)
328: {
329: mpn_bdivmod (up, up, usize, vp, vsize, d); /* Result > 0. */
330: d /= (unsigned long int)BITS_PER_MP_LIMB, up += d, usize -= d;
331: }
332: }
333:
334: /* Start binary GCD. */
335: do
336: {
337: mp_size_t zeros;
338:
339: /* Make sure U is odd. */
340: MPN_NORMALIZE (up, usize);
341: while (up[0] == 0)
342: up += 1, usize -= 1;
343: if ((up[0] & 1) == 0)
344: {
345: unsigned long int r;
346: count_trailing_zeros (r, up[0]);
347: mpn_rshift (up, up, usize, r);
348: usize -= (up[usize-1] == 0);
349: }
350:
351: /* Keep usize >= vsize. */
352: if (usize < vsize)
353: SWAP_MPN (up, usize, vp, vsize);
354:
355: if (usize <= 2) /* Double precision. */
356: {
357: if (vsize == 1)
358: vp[0] = mpn_gcd_1 (up, usize, vp[0]);
359: else
360: vsize = gcd_2 (vp, up);
361: break; /* Binary GCD done. */
362: }
363:
364: /* Count number of low zero limbs of U - V. */
365: for (zeros = 0; up[zeros] == vp[zeros] && ++zeros != vsize; )
366: continue;
367:
368: /* If U < V, swap U and V; in any case, subtract V from U. */
369: if (zeros == vsize) /* Subtract done. */
370: up += zeros, usize -= zeros;
371: else if (usize == vsize)
372: {
373: mp_size_t size = vsize;
374: do
375: size--;
376: while (up[size] == vp[size]);
377: if (up[size] < vp[size]) /* usize == vsize. */
378: SWAP_PTR (up, vp);
379: up += zeros, usize = size + 1 - zeros;
380: mpn_sub_n (up, up, vp + zeros, usize);
381: }
382: else
383: {
384: mp_size_t size = vsize - zeros;
385: up += zeros, usize -= zeros;
386: if (mpn_sub_n (up, up, vp + zeros, size))
387: {
388: while (up[size] == 0) /* Propagate borrow. */
389: up[size++] = -(mp_limb_t)1;
390: up[size] -= 1;
391: }
392: }
393: }
394: while (usize); /* End binary GCD. */
395: }
396:
397: done:
398: if (vp != gp)
399: MPN_COPY (gp, vp, vsize);
400: TMP_FREE (marker);
401: return vsize;
402: }
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 / mpn / generic