1 /*
2 Copyright (C) 2018 Daniel Schultz
3
4 This file is part of FLINT.
5
6 FLINT is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "nmod_poly.h"
13 #include "nmod_mpoly.h"
14
15 /*
16 Convert B to A if the degrees of A are <= expected_deg
17 If not, return 0 and set A to 0.
18 */
nmod_mpoly_convert_from_nmod_mpolyd_degbound(nmod_mpoly_t A,const nmod_mpoly_ctx_t ctx,const nmod_mpolyd_t B,const nmod_mpolyd_ctx_t dctx,slong * expect_deg)19 int nmod_mpoly_convert_from_nmod_mpolyd_degbound(
20 nmod_mpoly_t A, const nmod_mpoly_ctx_t ctx,
21 const nmod_mpolyd_t B, const nmod_mpolyd_ctx_t dctx,
22 slong * expect_deg)
23 {
24 int ret;
25 slong off, j, k, N;
26 slong bits, nvars = ctx->minfo->nvars;
27 slong Alen;
28 slong * rexpect_deg, * perm = dctx->perm;
29 slong perm_nontrivial;
30 ulong topmask, outrange;
31 ulong * exps, * pcurexp, * pexps, * rangemask;
32 TMP_INIT;
33
34 FLINT_ASSERT(nvars == B->nvars);
35 FLINT_ASSERT(nvars <= FLINT_BITS);
36
37 TMP_START;
38 exps = (ulong *) TMP_ALLOC(nvars*sizeof(ulong));
39 rexpect_deg = (slong *) TMP_ALLOC(nvars*sizeof(slong));
40
41 /* find bits needed for the result */
42 off = 1;
43 perm_nontrivial = 0;
44 for (j = 0; j < nvars; j++)
45 {
46 rexpect_deg[j] = expect_deg[perm[j]];
47 FLINT_ASSERT(rexpect_deg[j] >= 0);
48 off *= B->deg_bounds[j];
49 exps[perm[j]] = rexpect_deg[j];
50 perm_nontrivial |= j ^ perm[j];
51 }
52
53 FLINT_ASSERT(off <= B->coeff_alloc);
54
55 bits = mpoly_exp_bits_required_ui(exps, ctx->minfo);
56 bits = mpoly_fix_bits(bits, ctx->minfo);
57 N = mpoly_words_per_exp(bits, ctx->minfo);
58
59 /* we are going to push back terms manually */
60 Alen = 0;
61 nmod_mpoly_zero(A, ctx);
62 nmod_mpoly_fit_bits(A, bits, ctx);
63 A->bits = bits;
64
65 /* find exponent vector for all variables */
66 pexps = (ulong *) TMP_ALLOC(N*nvars*sizeof(ulong));
67 for (k = 0; k < nvars; k++)
68 {
69 for (j = 0; j < nvars; j++)
70 exps[perm[j]] = (j == k);
71 mpoly_set_monomial_ui(pexps + k*N, exps, bits, ctx->minfo);
72 }
73
74 /* get most significant exponent in exps and its vector in ptempexp */
75 off--;
76 pcurexp = (ulong *) TMP_ALLOC(N*sizeof(ulong));
77 rangemask = (ulong *) TMP_ALLOC(nvars*sizeof(ulong));
78 outrange = 0;
79 mpoly_monomial_zero(pcurexp, N);
80 k = off;
81 for (j = nvars - 1; j >= 0; j--)
82 {
83 exps[j] = k % B->deg_bounds[j];
84 rangemask[j] = UWORD(1) << j;
85 outrange ^= (outrange ^ FLINT_SIGN_EXT(rexpect_deg[j]
86 - exps[j])) & rangemask[j];
87 k = k / B->deg_bounds[j];
88 mpoly_monomial_madd_inplace_mp(pcurexp, exps[j], pexps + N*j, N);
89 }
90
91 /* scan down through the exponents */
92 topmask = 0;
93 for (; off >= 0; off--)
94 {
95 if (B->coeffs[off] != UWORD(0))
96 {
97 if (outrange)
98 goto failed_out_range;
99
100 _nmod_mpoly_fit_length(&A->coeffs, &A->exps, &A->alloc, Alen + 1, N);
101 A->coeffs[Alen] = B->coeffs[off];
102 mpoly_monomial_set(A->exps + N*Alen, pcurexp, N);
103 topmask |= (A->exps + N*Alen)[N - 1];
104 Alen++;
105 }
106
107 j = nvars - 1;
108 do {
109 --exps[j];
110 outrange ^= (outrange ^ FLINT_SIGN_EXT(rexpect_deg[j]
111 - exps[j])) & rangemask[j];
112 if ((slong)(exps[j]) < WORD(0))
113 {
114 FLINT_ASSERT(off == 0 || j > 0);
115 FLINT_ASSERT(exps[j] == -UWORD(1));
116 exps[j] = B->deg_bounds[j] - 1;
117 outrange ^= (outrange ^ FLINT_SIGN_EXT(rexpect_deg[j]
118 - exps[j])) & rangemask[j];
119 mpoly_monomial_madd_inplace_mp(pcurexp, exps[j], pexps + N*j, N);
120 } else
121 {
122 mpoly_monomial_sub_mp(pcurexp, pcurexp, pexps + N*j, N);
123 break;
124 }
125 } while (--j >= 0);
126 }
127 _nmod_mpoly_set_length(A, Alen, ctx);
128
129 /* sort the exponents if needed */
130 if (ctx->minfo->ord != ORD_LEX || perm_nontrivial != WORD(0))
131 {
132 slong msb;
133 mpoly_get_cmpmask(pcurexp, N, bits, ctx->minfo);
134 if (topmask != UWORD(0))
135 {
136 count_leading_zeros(msb, topmask);
137 msb = (FLINT_BITS - 1)^msb;
138 } else
139 {
140 msb = -WORD(1);
141 }
142 if (N == 1) {
143 if (msb >= WORD(0))
144 {
145 _nmod_mpoly_radix_sort1(A, 0, A->length,
146 msb, pcurexp[0], topmask);
147 }
148 } else {
149 _nmod_mpoly_radix_sort(A, 0, A->length,
150 (N - 1)*FLINT_BITS + msb, N, pcurexp);
151 }
152 }
153
154 ret = 1;
155
156 done:
157 TMP_END;
158 return ret;
159
160 failed_out_range:
161 ret = 0;
162 _nmod_mpoly_set_length(A, WORD(0), ctx);
163 goto done;
164 }
165
166 /*
167 return -1 : function failed
168 0 : B does not divide A
169 1 : B divides A and quotient is in Q
170 */
nmod_mpoly_divides_dense(nmod_mpoly_t Q,const nmod_mpoly_t A,const nmod_mpoly_t B,const nmod_mpoly_ctx_t ctx)171 int nmod_mpoly_divides_dense(nmod_mpoly_t Q,
172 const nmod_mpoly_t A, const nmod_mpoly_t B,
173 const nmod_mpoly_ctx_t ctx)
174 {
175 int ret, success;
176 slong i;
177 slong nvars = ctx->minfo->nvars;
178 nmod_mpolyd_ctx_t dctx;
179 nmod_mpolyd_t Ad, Bd, Qd;
180 nmod_poly_t Au, Bu, Qu, Ru;
181 slong * Abounds, * Bbounds, * Qbounds, * Edegs;
182 TMP_INIT;
183
184 if (B->length == 0)
185 {
186 if (A->length == 0 || nmod_mpoly_ctx_modulus(ctx) == 1)
187 {
188 nmod_mpoly_set(Q, A, ctx);
189 return 1;
190 } else
191 flint_throw(FLINT_DIVZERO, "Divide by zero in nmod_mpoly_divides_dense");
192 }
193
194 if (A->length == 0)
195 {
196 nmod_mpoly_zero(Q, ctx);
197 return 1;
198 }
199
200 if (A->bits > FLINT_BITS || B->bits > FLINT_BITS || ctx->minfo->nvars > FLINT_BITS)
201 {
202 return -1;
203 }
204
205 TMP_START;
206
207 /* set the ordering of variables for the dense representation */
208 nmod_mpolyd_ctx_init(dctx, nvars);
209
210 /*
211 for each variable v we need to pack to degree deg_v(A)
212 except for the outermost variable
213 */
214 Abounds = (slong *) TMP_ALLOC(ctx->minfo->nvars*sizeof(slong));
215 Bbounds = (slong *) TMP_ALLOC(ctx->minfo->nvars*sizeof(slong));
216 Qbounds = (slong *) TMP_ALLOC(ctx->minfo->nvars*sizeof(slong));
217 Edegs = (slong *) TMP_ALLOC(ctx->minfo->nvars*sizeof(slong));
218
219 mpoly_degrees_si(Abounds, A->exps, A->length, A->bits, ctx->minfo);
220 mpoly_degrees_si(Bbounds, B->exps, B->length, B->bits, ctx->minfo);
221 for (i = 0; i < ctx->minfo->nvars; i++)
222 {
223 /* if divides, expected degrees */
224 Edegs[i] = Abounds[i] - Bbounds[i];
225
226 if (Abounds[i] < Bbounds[i])
227 {
228 ret = 0;
229 nmod_mpoly_zero(Q, ctx);
230 goto cleanup_stage1;
231 }
232
233 if (i != dctx->perm[0])
234 {
235 /* variable of index i is not the outermost */
236 Qbounds[i] = Abounds[i] + 1;
237 Bbounds[i] = Abounds[i] + 1;
238 } else
239 {
240 /* variable of index i is the outermost */
241 Qbounds[i] = Abounds[i] - Bbounds[i] + 1;
242 Bbounds[i] = Bbounds[i] + 1;
243 }
244 Abounds[i] = Abounds[i] + 1;
245
246 if (Abounds[i] < WORD(0))
247 {
248 ret = -1;
249 nmod_mpoly_zero(Q, ctx);
250 goto cleanup_stage1;
251 }
252 }
253
254 nmod_mpolyd_init(Ad, nvars);
255 nmod_mpolyd_init(Bd, nvars);
256 nmod_mpolyd_init(Qd, nvars);
257
258 success = 1;
259 success = success && nmod_mpolyd_set_degbounds_perm(Ad, dctx, Abounds);
260 success = success && nmod_mpolyd_set_degbounds_perm(Bd, dctx, Bbounds);
261 success = success && nmod_mpolyd_set_degbounds_perm(Qd, dctx, Qbounds);
262 if (!success)
263 {
264 ret = -1;
265 goto cleanup_stage2;
266 }
267
268 nmod_mpoly_convert_to_nmod_mpolyd_degbound(Ad, dctx, A, ctx);
269 nmod_mpoly_convert_to_nmod_mpolyd_degbound(Bd, dctx, B, ctx);
270
271 /* let Au and Bu borrow Ad and Bd */
272 Au->alloc = Ad->coeff_alloc;
273 Au->coeffs = Ad->coeffs;
274 Au->length = nmod_mpolyd_length(Ad);
275 Au->mod.n = ctx->ffinfo->mod.n;
276 Au->mod.ninv = ctx->ffinfo->mod.ninv;
277 Au->mod.norm = ctx->ffinfo->mod.norm;
278
279 Bu->alloc = Bd->coeff_alloc;
280 Bu->coeffs = Bd->coeffs;
281 Bu->length = nmod_mpolyd_length(Bd);
282 Bu->mod.n = ctx->ffinfo->mod.n;
283 Bu->mod.ninv = ctx->ffinfo->mod.ninv;
284 Bu->mod.norm = ctx->ffinfo->mod.norm;
285
286 /* manually move Qd to Qu */
287 Qu->alloc = Qd->coeff_alloc;
288 Qu->coeffs = Qd->coeffs;
289 Qu->length = 0;
290 Qu->mod.n = ctx->ffinfo->mod.n;
291 Qu->mod.ninv = ctx->ffinfo->mod.ninv;
292 Qu->mod.norm = ctx->ffinfo->mod.norm;
293
294 nmod_poly_init(Ru, ctx->ffinfo->mod.n);
295 nmod_poly_divrem(Qu, Ru, Au, Bu);
296 if (!nmod_poly_is_zero(Ru))
297 {
298 ret = 0;
299 goto cleanup_stage3;
300 }
301 nmod_poly_clear(Ru);
302
303 /* manually move Qu to Qd */
304 Qd->coeff_alloc = Qu->alloc;
305 Qd->coeffs = Qu->coeffs;
306 for (i = Qu->length; i < Qd->coeff_alloc; i++)
307 Qd->coeffs[i] = UWORD(0);
308
309 nmod_mpolyd_clear(Bd);
310 nmod_mpolyd_clear(Ad);
311 ret = nmod_mpoly_convert_from_nmod_mpolyd_degbound(Q, ctx, Qd, dctx, Edegs);
312 nmod_mpolyd_clear(Qd);
313 nmod_mpolyd_ctx_clear(dctx);
314
315 done:
316 TMP_END;
317 return ret;
318
319 cleanup_stage3:
320 nmod_poly_clear(Ru);
321
322 cleanup_stage2:
323 nmod_mpolyd_clear(Ad);
324 nmod_mpolyd_clear(Bd);
325 nmod_mpolyd_clear(Qd);
326
327 cleanup_stage1:
328 nmod_mpolyd_ctx_clear(dctx);
329 goto done;
330 }
331