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