1 /* mpn_gcd_1 -- mpn and limb greatest common divisor. 2 3 Copyright 1994, 1996, 2000, 2001, 2009, 2012 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 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. 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 Lesser General Public 15 License for more details. 16 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/. */ 19 20 #include "gmp.h" 21 #include "gmp-impl.h" 22 #include "longlong.h" 23 24 #ifndef GCD_1_METHOD 25 #define GCD_1_METHOD 2 26 #endif 27 28 #define USE_ZEROTAB 0 29 30 #if USE_ZEROTAB 31 #define MAXSHIFT 4 32 #define MASK ((1 << MAXSHIFT) - 1) 33 static const unsigned char zerotab[1 << MAXSHIFT] = 34 { 35 #if MAXSHIFT > 4 36 5, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0, 37 #endif 38 4, 0, 1, 0, 2, 0, 1, 0, 3, 0, 1, 0, 2, 0, 1, 0 39 }; 40 #endif 41 42 /* Does not work for U == 0 or V == 0. It would be tough to make it work for 43 V == 0 since gcd(x,0) = x, and U does not generally fit in an mp_limb_t. 44 45 The threshold for doing u%v when size==1 will vary by CPU according to 46 the speed of a division and the code generated for the main loop. Any 47 tuning for this is left to a CPU specific implementation. */ 48 49 mp_limb_t 50 mpn_gcd_1 (mp_srcptr up, mp_size_t size, mp_limb_t vlimb) 51 { 52 mp_limb_t ulimb; 53 unsigned long zero_bits, u_low_zero_bits; 54 55 ASSERT (size >= 1); 56 ASSERT (vlimb != 0); 57 ASSERT_MPN_NONZERO_P (up, size); 58 59 ulimb = up[0]; 60 61 /* Need vlimb odd for modexact, want it odd to get common zeros. */ 62 count_trailing_zeros (zero_bits, vlimb); 63 vlimb >>= zero_bits; 64 65 if (size > 1) 66 { 67 /* Must get common zeros before the mod reduction. If ulimb==0 then 68 vlimb already gives the common zeros. */ 69 if (ulimb != 0) 70 { 71 count_trailing_zeros (u_low_zero_bits, ulimb); 72 zero_bits = MIN (zero_bits, u_low_zero_bits); 73 } 74 75 ulimb = MPN_MOD_OR_MODEXACT_1_ODD (up, size, vlimb); 76 if (ulimb == 0) 77 goto done; 78 79 goto strip_u_maybe; 80 } 81 82 /* size==1, so up[0]!=0 */ 83 count_trailing_zeros (u_low_zero_bits, ulimb); 84 ulimb >>= u_low_zero_bits; 85 zero_bits = MIN (zero_bits, u_low_zero_bits); 86 87 /* make u bigger */ 88 if (vlimb > ulimb) 89 MP_LIMB_T_SWAP (ulimb, vlimb); 90 91 /* if u is much bigger than v, reduce using a division rather than 92 chipping away at it bit-by-bit */ 93 if ((ulimb >> 16) > vlimb) 94 { 95 ulimb %= vlimb; 96 if (ulimb == 0) 97 goto done; 98 goto strip_u_maybe; 99 } 100 101 ASSERT (ulimb & 1); 102 ASSERT (vlimb & 1); 103 104 #if GCD_1_METHOD == 1 105 while (ulimb != vlimb) 106 { 107 ASSERT (ulimb & 1); 108 ASSERT (vlimb & 1); 109 110 if (ulimb > vlimb) 111 { 112 ulimb -= vlimb; 113 do 114 { 115 ulimb >>= 1; 116 ASSERT (ulimb != 0); 117 strip_u_maybe: 118 ; 119 } 120 while ((ulimb & 1) == 0); 121 } 122 else /* vlimb > ulimb. */ 123 { 124 vlimb -= ulimb; 125 do 126 { 127 vlimb >>= 1; 128 ASSERT (vlimb != 0); 129 } 130 while ((vlimb & 1) == 0); 131 } 132 } 133 #else 134 # if GCD_1_METHOD == 2 135 136 ulimb >>= 1; 137 vlimb >>= 1; 138 139 while (ulimb != vlimb) 140 { 141 int c; 142 mp_limb_t t; 143 mp_limb_t vgtu; 144 145 t = ulimb - vlimb; 146 vgtu = LIMB_HIGHBIT_TO_MASK (t); 147 148 /* v <-- min (u, v) */ 149 vlimb += (vgtu & t); 150 151 /* u <-- |u - v| */ 152 ulimb = (t ^ vgtu) - vgtu; 153 154 #if USE_ZEROTAB 155 /* Number of trailing zeros is the same no matter if we look at 156 * t or ulimb, but using t gives more parallelism. */ 157 c = zerotab[t & MASK]; 158 159 while (UNLIKELY (c == MAXSHIFT)) 160 { 161 ulimb >>= MAXSHIFT; 162 if (0) 163 strip_u_maybe: 164 vlimb >>= 1; 165 166 c = zerotab[ulimb & MASK]; 167 } 168 #else 169 if (0) 170 { 171 strip_u_maybe: 172 vlimb >>= 1; 173 t = ulimb; 174 } 175 count_trailing_zeros (c, t); 176 #endif 177 ulimb >>= (c + 1); 178 } 179 180 vlimb = (vlimb << 1) | 1; 181 # else 182 # error Unknown GCD_1_METHOD 183 # endif 184 #endif 185 186 done: 187 return vlimb << zero_bits; 188 } 189