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