1 /* mpn_gcd_1 -- mpn and limb greatest common divisor.
3 Copyright 1994, 1996, 2000, 2001 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 Lesser General Public License as published by
9 the Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
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 Lesser General Public
15 License for more details.
17 You should have received a copy of the GNU Lesser General Public License
18 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */
25 #define GCD_1_METHOD 2
31 static const unsigned char zerotab[16] = {
32 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0
36 /* Does not work for U == 0 or V == 0. It would be tough to make it work for
37 V == 0 since gcd(x,0) = x, and U does not generally fit in an mp_limb_t.
39 The threshold for doing u%v when size==1 will vary by CPU according to
40 the speed of a division and the code generated for the main loop. Any
41 tuning for this is left to a CPU specific implementation. */
44 mpn_gcd_1 (mp_srcptr up, mp_size_t size, mp_limb_t vlimb)
47 unsigned long zero_bits, u_low_zero_bits;
51 ASSERT_MPN_NONZERO_P (up, size);
55 /* Need vlimb odd for modexact, want it odd to get common zeros. */
56 count_trailing_zeros (zero_bits, vlimb);
61 /* Must get common zeros before the mod reduction. If ulimb==0 then
62 vlimb already gives the common zeros. */
65 count_trailing_zeros (u_low_zero_bits, ulimb);
66 zero_bits = MIN (zero_bits, u_low_zero_bits);
69 ulimb = MPN_MOD_OR_MODEXACT_1_ODD (up, size, vlimb);
76 /* size==1, so up[0]!=0 */
77 count_trailing_zeros (u_low_zero_bits, ulimb);
78 ulimb >>= u_low_zero_bits;
79 zero_bits = MIN (zero_bits, u_low_zero_bits);
83 MP_LIMB_T_SWAP (ulimb, vlimb);
85 /* if u is much bigger than v, reduce using a division rather than
86 chipping away at it bit-by-bit */
87 if ((ulimb >> 16) > vlimb)
99 while (ulimb != vlimb)
114 while ((ulimb & 1) == 0);
116 else /* vlimb > ulimb. */
124 while ((vlimb & 1) == 0);
128 # if GCD_1_METHOD == 2
133 while (ulimb != vlimb)
136 mp_limb_t t = ulimb - vlimb;
137 mp_limb_t vgtu = LIMB_HIGHBIT_TO_MASK (t);
139 /* v <-- min (u, v) */
143 ulimb = (t ^ vgtu) - vgtu;
146 /* Number of trailing zeros is the same no matter if we look at
147 * t or ulimb, but using t gives more parallelism. */
150 while (UNLIKELY (c == 4))
157 c = zerotab[ulimb & 15];
166 count_trailing_zeros (c, t);
171 vlimb = (vlimb << 1) | 1;
173 # error Unknown GCD_1_METHOD
178 return vlimb << zero_bits;