1 /* mpn_divisible_p -- mpn by mpn divisibility test 2 3 THE FUNCTIONS IN THIS FILE ARE FOR INTERNAL USE ONLY. THEY'RE ALMOST 4 CERTAIN TO BE SUBJECT TO INCOMPATIBLE CHANGES OR DISAPPEAR COMPLETELY IN 5 FUTURE GNU MP RELEASES. 6 7 Copyright 2001, 2002, 2005, 2009 Free Software Foundation, Inc. 8 9 This file is part of the GNU MP Library. 10 11 The GNU MP Library is free software; you can redistribute it and/or modify 12 it under the terms of the GNU Lesser General Public License as published by 13 the Free Software Foundation; either version 3 of the License, or (at your 14 option) any later version. 15 16 The GNU MP Library is distributed in the hope that it will be useful, but 17 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 18 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 19 License for more details. 20 21 You should have received a copy of the GNU Lesser General Public License 22 along with the GNU MP Library. If not, see http://www.gnu.org/licenses/. */ 23 24 #include "gmp.h" 25 #include "gmp-impl.h" 26 #include "longlong.h" 27 28 29 /* Determine whether {ap,an} is divisible by {dp,dn}. Must have both 30 operands normalized, meaning high limbs non-zero, except that an==0 is 31 allowed. 32 33 There usually won't be many low zero bits on d, but the checks for this 34 are fast and might pick up a few operand combinations, in particular they 35 might reduce d to fit the single-limb mod_1/modexact_1 code. 36 37 Future: 38 39 Getting the remainder limb by limb would make an early exit possible on 40 finding a non-zero. This would probably have to be bdivmod style so 41 there's no addback, but it would need a multi-precision inverse and so 42 might be slower than the plain method (on small sizes at least). 43 44 When d must be normalized (shifted to high bit set), it's possible to 45 just append a low zero limb to "a" rather than bit-shifting as 46 mpn_tdiv_qr does internally, so long as it's already been checked that a 47 has at least as many trailing zeros bits as d. Or equivalently, pass 48 qxn==1 to mpn_tdiv_qr, if/when it accepts that. */ 49 50 int 51 mpn_divisible_p (mp_srcptr ap, mp_size_t an, 52 mp_srcptr dp, mp_size_t dn) 53 { 54 mp_limb_t alow, dlow, dmask; 55 mp_ptr qp, rp, tp; 56 mp_size_t i; 57 mp_limb_t di; 58 unsigned twos; 59 TMP_DECL; 60 61 ASSERT (an >= 0); 62 ASSERT (an == 0 || ap[an-1] != 0); 63 ASSERT (dn >= 1); 64 ASSERT (dp[dn-1] != 0); 65 ASSERT_MPN (ap, an); 66 ASSERT_MPN (dp, dn); 67 68 /* When a<d only a==0 is divisible. 69 Notice this test covers all cases of an==0. */ 70 if (an < dn) 71 return (an == 0); 72 73 /* Strip low zero limbs from d, requiring a==0 on those. */ 74 for (;;) 75 { 76 alow = *ap; 77 dlow = *dp; 78 79 if (dlow != 0) 80 break; 81 82 if (alow != 0) 83 return 0; /* a has fewer low zero limbs than d, so not divisible */ 84 85 /* a!=0 and d!=0 so won't get to n==0 */ 86 an--; ASSERT (an >= 1); 87 dn--; ASSERT (dn >= 1); 88 ap++; 89 dp++; 90 } 91 92 /* a must have at least as many low zero bits as d */ 93 dmask = LOW_ZEROS_MASK (dlow); 94 if ((alow & dmask) != 0) 95 return 0; 96 97 if (dn == 1) 98 { 99 if (ABOVE_THRESHOLD (an, BMOD_1_TO_MOD_1_THRESHOLD)) 100 return mpn_mod_1 (ap, an, dlow) == 0; 101 102 count_trailing_zeros (twos, dlow); 103 dlow >>= twos; 104 return mpn_modexact_1_odd (ap, an, dlow) == 0; 105 } 106 107 if (dn == 2) 108 { 109 mp_limb_t dsecond = dp[1]; 110 if (dsecond <= dmask) 111 { 112 count_trailing_zeros (twos, dlow); 113 dlow = (dlow >> twos) | (dsecond << (GMP_NUMB_BITS-twos)); 114 ASSERT_LIMB (dlow); 115 return MPN_MOD_OR_MODEXACT_1_ODD (ap, an, dlow) == 0; 116 } 117 } 118 119 /* Should we compute Q = A * D^(-1) mod B^k, 120 R = A - Q * D mod B^k 121 here, for some small values of k? Then check if R = 0 (mod B^k). */ 122 123 /* We could also compute A' = A mod T and D' = D mod P, for some 124 P = 3 * 5 * 7 * 11 ..., and then check if any prime factor from P 125 dividing D' also divides A'. */ 126 127 TMP_MARK; 128 129 rp = TMP_ALLOC_LIMBS (an + 1); 130 qp = TMP_ALLOC_LIMBS (an - dn + 1); /* FIXME: Could we avoid this */ 131 132 count_trailing_zeros (twos, dp[0]); 133 134 if (twos != 0) 135 { 136 tp = TMP_ALLOC_LIMBS (dn); 137 ASSERT_NOCARRY (mpn_rshift (tp, dp, dn, twos)); 138 dp = tp; 139 140 ASSERT_NOCARRY (mpn_rshift (rp, ap, an, twos)); 141 } 142 else 143 { 144 MPN_COPY (rp, ap, an); 145 } 146 if (rp[an - 1] >= dp[dn - 1]) 147 { 148 rp[an] = 0; 149 an++; 150 } 151 else if (an == dn) 152 { 153 TMP_FREE; 154 return 0; 155 } 156 157 ASSERT (an > dn); /* requirement of functions below */ 158 159 if (BELOW_THRESHOLD (dn, DC_BDIV_QR_THRESHOLD) || 160 BELOW_THRESHOLD (an - dn, DC_BDIV_QR_THRESHOLD)) 161 { 162 binvert_limb (di, dp[0]); 163 mpn_sbpi1_bdiv_qr (qp, rp, an, dp, dn, -di); 164 rp += an - dn; 165 } 166 else if (BELOW_THRESHOLD (dn, MU_BDIV_QR_THRESHOLD)) 167 { 168 binvert_limb (di, dp[0]); 169 mpn_dcpi1_bdiv_qr (qp, rp, an, dp, dn, -di); 170 rp += an - dn; 171 } 172 else 173 { 174 tp = TMP_ALLOC_LIMBS (mpn_mu_bdiv_qr_itch (an, dn)); 175 mpn_mu_bdiv_qr (qp, rp, rp, an, dp, dn, tp); 176 } 177 178 /* test for {rp,dn} zero or non-zero */ 179 i = 0; 180 do 181 { 182 if (rp[i] != 0) 183 { 184 TMP_FREE; 185 return 0; 186 } 187 } 188 while (++i < dn); 189 190 TMP_FREE; 191 return 1; 192 } 193