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, 2014, 2017, 2018 Free Software
8 Foundation, Inc.
9
10 This file is part of the GNU MP Library.
11
12 The GNU MP Library is free software; you can redistribute it and/or modify
13 it under the terms of either:
14
15 * the GNU Lesser General Public License as published by the Free
16 Software Foundation; either version 3 of the License, or (at your
17 option) any later version.
18
19 or
20
21 * the GNU General Public License as published by the Free Software
22 Foundation; either version 2 of the License, or (at your option) any
23 later version.
24
25 or both in parallel, as here.
26
27 The GNU MP Library is distributed in the hope that it will be useful, but
28 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
29 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
30 for more details.
31
32 You should have received copies of the GNU General Public License and the
33 GNU Lesser General Public License along with the GNU MP Library. If not,
34 see https://www.gnu.org/licenses/. */
35
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
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_limb_t di;
66 unsigned twos;
67 int c;
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 count_trailing_zeros (twos, dlow);
117 if (dn == 2)
118 {
119 mp_limb_t dsecond = dp[1];
120 if (dsecond <= dmask)
121 {
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 TMP_ALLOC_LIMBS_2 (rp, an + 1,
139 qp, an - dn + 1); /* FIXME: Could we avoid this? */
140
141 if (twos != 0)
142 {
143 tp = TMP_ALLOC_LIMBS (dn);
144 ASSERT_NOCARRY (mpn_rshift (tp, dp, dn, twos));
145 dp = tp;
146
147 ASSERT_NOCARRY (mpn_rshift (rp, ap, an, twos));
148 }
149 else
150 {
151 MPN_COPY (rp, ap, an);
152 }
153 if (rp[an - 1] >= dp[dn - 1])
154 {
155 rp[an] = 0;
156 an++;
157 }
158 else if (an == dn)
159 {
160 TMP_FREE;
161 return 0;
162 }
163
164 ASSERT (an > dn); /* requirement of functions below */
165
166 if (BELOW_THRESHOLD (dn, DC_BDIV_QR_THRESHOLD) ||
167 BELOW_THRESHOLD (an - dn, DC_BDIV_QR_THRESHOLD))
168 {
169 binvert_limb (di, dp[0]);
170 mpn_sbpi1_bdiv_qr (qp, rp, an, dp, dn, -di);
171 rp += an - dn;
172 }
173 else if (BELOW_THRESHOLD (dn, MU_BDIV_QR_THRESHOLD))
174 {
175 binvert_limb (di, dp[0]);
176 mpn_dcpi1_bdiv_qr (qp, rp, an, dp, dn, -di);
177 rp += an - dn;
178 }
179 else
180 {
181 tp = TMP_ALLOC_LIMBS (mpn_mu_bdiv_qr_itch (an, dn));
182 mpn_mu_bdiv_qr (qp, rp, rp, an, dp, dn, tp);
183 }
184
185 /* In general, bdiv may return either R = 0 or R = D when D divides
186 A. But R = 0 can happen only when A = 0, which we already have
187 excluded. Furthermore, R == D (mod B^{dn}) implies no carry, so
188 we don't need to check the carry returned from bdiv. */
189
190 MPN_CMP (c, rp, dp, dn);
191
192 TMP_FREE;
193 return c == 0;
194 }