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