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 
_nmod_mpoly_mul_dense(nmod_mpoly_t P,const nmod_mpoly_t A,fmpz * maxAfields,const nmod_mpoly_t B,fmpz * maxBfields,const nmod_mpoly_ctx_t ctx)15 int _nmod_mpoly_mul_dense(nmod_mpoly_t P,
16                                  const nmod_mpoly_t A, fmpz * maxAfields,
17                                  const nmod_mpoly_t B, fmpz * maxBfields,
18                                                     const nmod_mpoly_ctx_t ctx)
19 {
20     int success = 1;
21     slong i;
22     slong nvars = ctx->minfo->nvars;
23     nmod_mpolyd_ctx_t dctx;
24     nmod_mpolyd_t Ad, Bd, Pd;
25     nmod_poly_t Au, Bu, Pu;
26     slong * Abounds, * Bbounds, * Pbounds;
27     TMP_INIT;
28 
29     FLINT_ASSERT(A->length != 0);
30     FLINT_ASSERT(B->length != 0);
31 
32     FLINT_ASSERT(A->bits <= FLINT_BITS);
33     FLINT_ASSERT(B->bits <= FLINT_BITS);
34 
35     TMP_START;
36 
37     /* set the ordering of variables for the dense representation */
38     nmod_mpolyd_ctx_init(dctx, nvars);
39 
40     /*
41         for each variable v except for the outermost variable,
42         we need to pack to degree deg_v(A) + deg_v(B)
43     */
44     Abounds = (slong *) TMP_ALLOC(ctx->minfo->nvars*sizeof(slong));
45     Bbounds = (slong *) TMP_ALLOC(ctx->minfo->nvars*sizeof(slong));
46     Pbounds = (slong *) TMP_ALLOC(ctx->minfo->nvars*sizeof(slong));
47     mpoly_get_monomial_ui_unpacked_ffmpz((ulong *)Abounds, maxAfields, ctx->minfo);
48     mpoly_get_monomial_ui_unpacked_ffmpz((ulong *)Bbounds, maxBfields, ctx->minfo);
49     for (i = 0; i < ctx->minfo->nvars; i++)
50     {
51         Abounds[i] = Abounds[i] + 1;
52         Bbounds[i] = Bbounds[i] + 1;
53         Pbounds[i] = Abounds[i] + Bbounds[i] - 1;
54         if ((Abounds[i] | Bbounds[i] | Pbounds[i]) < WORD(0))
55         {
56             goto failed_stage1;
57         }
58         if (i != dctx->perm[0])
59         {
60             /* variable of index i is not the outermost */
61             Abounds[i] = Pbounds[i];
62             Bbounds[i] = Pbounds[i];
63         }
64     }
65 
66     nmod_mpolyd_init(Ad, nvars);
67     nmod_mpolyd_init(Bd, nvars);
68     nmod_mpolyd_init(Pd, nvars);
69 
70     success = 1;
71     success = success && nmod_mpolyd_set_degbounds_perm(Ad, dctx, Abounds);
72     success = success && nmod_mpolyd_set_degbounds_perm(Bd, dctx, Bbounds);
73     success = success && nmod_mpolyd_set_degbounds_perm(Pd, dctx, Pbounds);
74     if (!success)
75     {
76         goto failed_stage2;
77     }
78 
79     nmod_mpoly_convert_to_nmod_mpolyd_degbound(Ad, dctx, A, ctx);
80     nmod_mpoly_convert_to_nmod_mpolyd_degbound(Bd, dctx, B, ctx);
81 
82     /* let Au and Bu borrow Ad and Bd */
83     Au->alloc  = Ad->coeff_alloc;
84     Au->coeffs = Ad->coeffs;
85     Au->length = nmod_mpolyd_length(Ad);
86     Au->mod.n    = ctx->ffinfo->mod.n;
87     Au->mod.ninv = ctx->ffinfo->mod.ninv;
88     Au->mod.norm = ctx->ffinfo->mod.norm;
89 
90     Bu->alloc  = Bd->coeff_alloc;
91     Bu->coeffs = Bd->coeffs;
92     Bu->length = nmod_mpolyd_length(Bd);
93     Bu->mod.n    = ctx->ffinfo->mod.n;
94     Bu->mod.ninv = ctx->ffinfo->mod.ninv;
95     Bu->mod.norm = ctx->ffinfo->mod.norm;
96 
97     /* manually move P to Pu */
98     Pu->alloc  = Pd->coeff_alloc;
99     Pu->coeffs = Pd->coeffs;
100     Pu->length = 0;
101     Pu->mod.n    = ctx->ffinfo->mod.n;
102     Pu->mod.ninv = ctx->ffinfo->mod.ninv;
103     Pu->mod.norm = ctx->ffinfo->mod.norm;
104 
105     nmod_poly_mul(Pu, Au, Bu);
106 
107     /* manually move Pu to P */
108     Pd->coeff_alloc = Pu->alloc;
109     Pd->coeffs      = Pu->coeffs;
110     for (i = Pu->length; i < Pd->coeff_alloc; i++)
111         Pd->coeffs[i] = UWORD(0);
112 
113     nmod_mpolyd_clear(Bd);
114     nmod_mpolyd_clear(Ad);
115     nmod_mpoly_convert_from_nmod_mpolyd(P, ctx, Pd, dctx);
116     nmod_mpolyd_clear(Pd);
117 
118     nmod_mpolyd_ctx_clear(dctx);
119 
120 done:
121     TMP_END;
122     return success;
123 
124 failed_stage2:
125     nmod_mpolyd_clear(Ad);
126     nmod_mpolyd_clear(Bd);
127     nmod_mpolyd_clear(Pd);
128 
129 failed_stage1:
130     nmod_mpolyd_ctx_clear(dctx);
131     success = 0;
132     goto done;
133 }
134 
135 
nmod_mpoly_mul_dense(nmod_mpoly_t A,const nmod_mpoly_t B,const nmod_mpoly_t C,const nmod_mpoly_ctx_t ctx)136 int nmod_mpoly_mul_dense(nmod_mpoly_t A, const nmod_mpoly_t B,
137                               const nmod_mpoly_t C, const nmod_mpoly_ctx_t ctx)
138 {
139     slong i;
140     int success;
141     fmpz * maxBfields, * maxCfields;
142     TMP_INIT;
143 
144     if (B->length == 0 || C->length == 0)
145     {
146         nmod_mpoly_zero(A, ctx);
147         return 1;
148     }
149 
150     if (B->bits > FLINT_BITS || C->bits > FLINT_BITS)
151     {
152         return 0;
153     }
154 
155     TMP_START;
156 
157     maxBfields = (fmpz *) TMP_ALLOC(ctx->minfo->nfields*sizeof(fmpz));
158     maxCfields = (fmpz *) TMP_ALLOC(ctx->minfo->nfields*sizeof(fmpz));
159     for (i = 0; i < ctx->minfo->nfields; i++)
160     {
161         fmpz_init(maxBfields + i);
162         fmpz_init(maxCfields + i);
163     }
164     mpoly_max_fields_fmpz(maxBfields, B->exps, B->length, B->bits, ctx->minfo);
165     mpoly_max_fields_fmpz(maxCfields, C->exps, C->length, C->bits, ctx->minfo);
166 
167     success = _nmod_mpoly_mul_dense(A, B, maxBfields, C, maxCfields, ctx);
168 
169     for (i = 0; i < ctx->minfo->nfields; i++)
170     {
171         fmpz_clear(maxBfields + i);
172         fmpz_clear(maxCfields + i);
173     }
174 
175     TMP_END;
176     return success;
177 }
178