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
mpn_divisible_p(mp_srcptr ap,mp_size_t an,mp_srcptr dp,mp_size_t dn)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