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_mpoly.h"
13 #include "fmpz_mpoly.h"
14 
nmod_mpoly_pow_ui(nmod_mpoly_t A,const nmod_mpoly_t B,ulong k,const nmod_mpoly_ctx_t ctx)15 int nmod_mpoly_pow_ui(nmod_mpoly_t A, const nmod_mpoly_t B,
16                                            ulong k, const nmod_mpoly_ctx_t ctx)
17 {
18     int success = 1;
19     slong i, exp_bits, N;
20     fmpz * maxBfields;
21     ulong * cmpmask;
22     ulong * Bexp;
23     int freeBexp;
24     fmpz_mpoly_t T;
25     slong Tlen;
26     fmpz * Bcoeffs_fmpz;
27     TMP_INIT;
28 
29     if (k == 0)
30     {
31         nmod_mpoly_set_ui(A, ctx->ffinfo->mod.n > 1, ctx);
32         return 1;
33     }
34 
35     if (B->length == 0)
36     {
37         nmod_mpoly_zero(A, ctx);
38         return 1;
39     }
40 
41     if (k == 1)
42     {
43         nmod_mpoly_set(A, B, ctx);
44         return 1;
45     }
46 
47     if (k == 2)
48     {
49         nmod_mpoly_mul_johnson(A, B, B, ctx);
50         return 1;
51     }
52 
53     if (A == B)
54     {
55         nmod_mpoly_t T;
56         nmod_mpoly_init(T, ctx);
57         success = nmod_mpoly_pow_ui(T, B, k, ctx);
58         nmod_mpoly_swap(A, T, ctx);
59         nmod_mpoly_clear(T, ctx);
60         return success;
61     }
62 
63     TMP_START;
64 
65     maxBfields = (fmpz *) TMP_ALLOC(ctx->minfo->nfields*sizeof(fmpz));
66     for (i = 0; i < ctx->minfo->nfields; i++)
67         fmpz_init(maxBfields + i);
68 
69     mpoly_max_fields_fmpz(maxBfields, B->exps, B->length, B->bits, ctx->minfo);
70 
71     _fmpz_vec_scalar_mul_ui(maxBfields, maxBfields, ctx->minfo->nfields, k);
72 
73     exp_bits = _fmpz_vec_max_bits(maxBfields, ctx->minfo->nfields);
74     exp_bits = FLINT_MAX(MPOLY_MIN_BITS, exp_bits + 1);
75     exp_bits = FLINT_MAX(exp_bits, B->bits);
76     exp_bits = mpoly_fix_bits(exp_bits, ctx->minfo);
77 
78     for (i = 0; i < ctx->minfo->nfields; i++)
79         fmpz_clear(maxBfields + i);
80 
81     N = mpoly_words_per_exp(exp_bits, ctx->minfo);
82 
83     freeBexp = 0;
84     Bexp = B->exps;
85     if (exp_bits > B->bits)
86     {
87         freeBexp = 1;
88         Bexp = (ulong *) flint_malloc(N*B->length*sizeof(ulong));
89         mpoly_repack_monomials(Bexp, exp_bits, B->exps, B->bits,
90                                                         B->length, ctx->minfo);
91     }
92 
93     if (B->length == 1)
94     {
95         /* powering a monomial */
96         nmod_mpoly_fit_length(A, WORD(1), ctx);
97         nmod_mpoly_fit_bits(A, exp_bits, ctx);
98         A->bits = exp_bits;
99 
100         if (exp_bits <= FLINT_BITS)
101             mpoly_monomial_mul_ui(A->exps, Bexp, N, k);
102         else
103             mpoly_monomial_mul_ui_mp(A->exps, Bexp, N, k);
104 
105         A->coeffs[0] = nmod_pow_ui(B->coeffs[0], k, ctx->ffinfo->mod);
106         _nmod_mpoly_set_length(A, WORD(A->coeffs[0] != 0), ctx);
107     }
108     else
109     {
110         cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
111         mpoly_get_cmpmask(cmpmask, N, exp_bits, ctx->minfo);
112 
113         T->alloc = k*(B->length - 1) + 1;
114         T->coeffs = (fmpz *) flint_calloc(T->alloc, sizeof(fmpz));
115         T->exps   = (ulong *) flint_malloc(T->alloc*N*sizeof(ulong));
116         T->length = 0;
117         T->bits = exp_bits;
118 
119         Bcoeffs_fmpz = _fmpz_vec_init(B->length);
120         _fmpz_vec_set_nmod_vec(Bcoeffs_fmpz, B->coeffs, B->length, ctx->ffinfo->mod);
121 
122         if (ctx->ffinfo->mod.n > 99999 || !n_is_prime(ctx->ffinfo->mod.n))
123         {
124             slong Alen, Tlen;
125             ulong limit = (ulong)(WORD_MAX)/(ulong)(2*sizeof(fmpz));
126 
127             if (B->length > 1 && k > limit/(ulong)(B->length - 1))
128             {
129                 success = 0;
130             }
131             else
132             {
133                 Tlen = _fmpz_mpoly_pow_fps(&T->coeffs, &T->exps, &T->alloc,
134                            Bcoeffs_fmpz, Bexp, B->length, k, exp_bits, N, cmpmask);
135 
136                 nmod_mpoly_fit_length(A, Tlen, ctx);
137                 nmod_mpoly_fit_bits(A, exp_bits, ctx);
138                 A->bits = exp_bits;
139 
140                 Alen = 0;
141                 for (i = 0; i < Tlen; i++)
142                 {
143                     A->coeffs[Alen] = fmpz_fdiv_ui(T->coeffs + i, ctx->ffinfo->mod.n);
144                     mpoly_monomial_set(A->exps + N*Alen, T->exps + N*i, N);
145                     Alen += (A->coeffs[Alen] != UWORD(0));
146                 }
147 
148                 _nmod_mpoly_set_length(A, Alen, ctx);
149             }
150         }
151         else
152         {
153             ulong ne;
154             slong Slen, Ulen;
155             nmod_mpoly_t S, R, U;
156 
157             nmod_mpoly_init3(S, B->length, exp_bits, ctx);
158             nmod_mpoly_init3(U, B->length, exp_bits, ctx);
159             nmod_mpoly_init3(R, B->length, exp_bits, ctx);
160             nmod_mpoly_one(R, ctx);
161 
162             ne = UWORD(1);
163 
164             while (k > 0)
165             {
166                 ulong kmodn;
167                 NMOD_RED(kmodn, k, ctx->ffinfo->mod);
168 
169                 if (kmodn > 0)
170                 {
171                     /* will accomplish R *= B^(n^e*(k%n)) */
172 
173                     if (kmodn <= 2)
174                     {
175                         nmod_mpoly_fit_length(S, B->length, ctx);
176 
177                         /* S = B^(n^e) */
178                         for (i = 0; i < B->length; i++)
179                         {
180                             S->coeffs[i] = B->coeffs[i];
181                             mpoly_monomial_mul_ui_mp(S->exps + N*i,
182                                                         Bexp + N*i, N, ne);
183                         }
184                         _nmod_mpoly_set_length(S, B->length, ctx);
185 
186                         if (kmodn == 2)
187                         {
188                             /* R *= S */
189                             Ulen = _nmod_mpoly_mul_johnson(
190                                             &U->coeffs, &U->exps, &U->alloc,
191                                              R->coeffs, R->exps, R->length,
192                                              S->coeffs, S->exps, S->length,
193                                            exp_bits,  N, cmpmask, ctx->ffinfo);
194                             _nmod_mpoly_set_length(U, Ulen, ctx);
195                             nmod_mpoly_swap(R, U, ctx);
196                         }
197                     }
198                     else
199                     {
200                         /* S = B^(n^e*(k%n)) */
201                         Tlen = _fmpz_mpoly_pow_fps(
202                                               &T->coeffs, &T->exps, &T->alloc,
203                                                 Bcoeffs_fmpz, Bexp, B->length,
204                                                   kmodn, exp_bits, N, cmpmask);
205                         nmod_mpoly_fit_length(S, Tlen, ctx);
206                         Slen = 0;
207                         for (i = 0; i < Tlen; i++)
208                         {
209                             S->coeffs[Slen] = fmpz_fdiv_ui(T->coeffs + i,
210                                                            ctx->ffinfo->mod.n);
211                             mpoly_monomial_mul_ui_mp(S->exps + N*Slen,
212                                                      T->exps + N*i, N, ne);
213                             Slen += (S->coeffs[Slen] != UWORD(0));
214                         }
215                         _nmod_mpoly_set_length(S, Slen, ctx);
216                     }
217 
218                     /* R = R * S */
219                     if (nmod_mpoly_is_one(R, ctx))
220                     {
221                         nmod_mpoly_swap(R, S, ctx);
222                     }
223                     else
224                     {
225                         Ulen = _nmod_mpoly_mul_johnson(
226                                             &U->coeffs, &U->exps, &U->alloc,
227                                              R->coeffs, R->exps, R->length,
228                                              S->coeffs, S->exps, S->length,
229                                            exp_bits,  N, cmpmask, ctx->ffinfo);
230                         _nmod_mpoly_set_length(U, Ulen, ctx);
231                         nmod_mpoly_swap(R, U, ctx);
232                     }
233                 }
234 
235                 k = k/ctx->ffinfo->mod.n;
236                 ne = ne * ctx->ffinfo->mod.n;
237             }
238 
239             nmod_mpoly_swap(A, R, ctx);
240 
241             nmod_mpoly_clear(S, ctx);
242             nmod_mpoly_clear(U, ctx);
243             nmod_mpoly_clear(R, ctx);
244         }
245 
246         _fmpz_vec_clear(Bcoeffs_fmpz, B->length);
247 
248         for (i = 0; i < T->alloc; i++)
249             fmpz_clear(T->coeffs + i);
250 
251         flint_free(T->coeffs);
252         flint_free(T->exps);
253     }
254 
255     if (freeBexp)
256     {
257         flint_free(Bexp);
258     }
259 
260     TMP_END;
261 
262     return success;
263 }
264