1 /*
2 Copyright (C) 2020 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 <https://www.gnu.org/licenses/>.
10 */
11
12 #include "fq_nmod_mpoly.h"
13
_fq_nmod_mpoly_evaluate_one_fq_nmod_sp(fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,slong var,const fq_nmod_t val,const fq_nmod_mpoly_ctx_t ctx,n_poly_stack_t St)14 void _fq_nmod_mpoly_evaluate_one_fq_nmod_sp(
15 fq_nmod_mpoly_t A,
16 const fq_nmod_mpoly_t B,
17 slong var,
18 const fq_nmod_t val,
19 const fq_nmod_mpoly_ctx_t ctx,
20 n_poly_stack_t St)
21 {
22 slong d = fq_nmod_ctx_degree(ctx->fqctx);
23 slong i, N, off, shift;
24 ulong * cmpmask, * one;
25 slong Blen = B->length;
26 const mp_limb_t * Bcoeffs = B->coeffs;
27 const ulong * Bexps = B->exps;
28 flint_bitcnt_t bits = B->bits;
29 slong Alen;
30 mp_limb_t * Acoeffs;
31 ulong * Aexps;
32 ulong mask, k;
33 int need_sort = 0, cmp;
34 n_poly_struct * cache[3];
35 TMP_INIT;
36
37 FLINT_ASSERT(B->bits <= FLINT_BITS);
38
39 TMP_START;
40
41 n_poly_stack_fit_request(St, 3);
42 cache[0] = n_poly_stack_take_top(St);
43 cache[1] = n_poly_stack_take_top(St);
44 cache[2] = n_poly_stack_take_top(St);
45 n_fq_pow_cache_start_fq_nmod(val, cache[0], cache[1], cache[2], ctx->fqctx);
46
47 fq_nmod_mpoly_fit_length_reset_bits(A, Blen, bits, ctx);
48 Acoeffs = A->coeffs;
49 Aexps = A->exps;
50
51 mask = (-UWORD(1)) >> (FLINT_BITS - bits);
52 N = mpoly_words_per_exp(bits, ctx->minfo);
53 one = (ulong*) TMP_ALLOC(N*sizeof(ulong));
54 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
55 mpoly_gen_monomial_offset_shift_sp(one, &off, &shift, var, bits, ctx->minfo);
56 mpoly_get_cmpmask(cmpmask, N, bits, ctx->minfo);
57
58 Alen = 0;
59 for (i = 0; i < Blen; i++)
60 {
61 k = (Bexps[N*i + off] >> shift) & mask;
62 _n_fq_set(Acoeffs + d*Alen, Bcoeffs + d*i, d);
63 n_fq_pow_cache_mulpow_ui(Acoeffs + d*Alen, Bcoeffs + d*i, k,
64 cache[0], cache[1], cache[2], ctx->fqctx);
65 if (_n_fq_is_zero(Acoeffs + d*Alen, d))
66 continue;
67
68 mpoly_monomial_msub(Aexps + N*Alen, Bexps + N*i, k, one, N);
69 if (Alen < 1)
70 {
71 Alen = 1;
72 continue;
73 }
74 cmp = mpoly_monomial_cmp(Aexps + N*(Alen - 1), Aexps + N*Alen, N, cmpmask);
75 if (cmp != 0)
76 {
77 need_sort |= (cmp < 0);
78 Alen++;
79 continue;
80 }
81 _n_fq_add(Acoeffs + d*(Alen - 1), Acoeffs + d*(Alen - 1),
82 Acoeffs + d*Alen, d, fq_nmod_ctx_mod(ctx->fqctx));
83 Alen -= _n_fq_is_zero(Acoeffs + d*(Alen - 1), d);
84 }
85 A->length = Alen;
86
87 n_poly_stack_give_back(St, 3);
88
89 TMP_END;
90
91 if (need_sort)
92 {
93 fq_nmod_mpoly_sort_terms(A, ctx);
94 fq_nmod_mpoly_combine_like_terms(A, ctx);
95 }
96
97 FLINT_ASSERT(fq_nmod_mpoly_is_canonical(A, ctx));
98 }
99
100
101 /* exponents of B are multiprecision */
_fq_nmod_mpoly_evaluate_one_fq_nmod_mp(fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,slong var,const fq_nmod_t val,const fq_nmod_mpoly_ctx_t ctx,n_poly_stack_t St)102 static void _fq_nmod_mpoly_evaluate_one_fq_nmod_mp(
103 fq_nmod_mpoly_t A,
104 const fq_nmod_mpoly_t B,
105 slong var,
106 const fq_nmod_t val,
107 const fq_nmod_mpoly_ctx_t ctx,
108 n_poly_stack_t St)
109 {
110 slong d = fq_nmod_ctx_degree(ctx->fqctx);
111 slong i, N, off;
112 ulong * cmpmask, * one, * tmp;
113 slong Blen = B->length;
114 const mp_limb_t * Bcoeffs = B->coeffs;
115 const ulong * Bexps = B->exps;
116 flint_bitcnt_t bits = B->bits;
117 slong Alen;
118 mp_limb_t * Acoeffs;
119 ulong * Aexps;
120 fmpz_t k;
121 int need_sort = 0, cmp;
122 n_poly_struct * cache[3];
123 TMP_INIT;
124
125 FLINT_ASSERT((B->bits % FLINT_BITS) == 0);
126
127 TMP_START;
128
129 fmpz_init(k);
130
131 n_poly_stack_fit_request(St, 3);
132 cache[0] = n_poly_stack_take_top(St);
133 cache[1] = n_poly_stack_take_top(St);
134 cache[2] = n_poly_stack_take_top(St);
135 n_fq_pow_cache_start_fq_nmod(val, cache[0], cache[1], cache[2], ctx->fqctx);
136
137 fq_nmod_mpoly_fit_length_reset_bits(A, Blen, bits, ctx);
138 Acoeffs = A->coeffs;
139 Aexps = A->exps;
140
141 N = mpoly_words_per_exp(bits, ctx->minfo);
142 one = (ulong*) TMP_ALLOC(3*N*sizeof(ulong));
143 cmpmask = one + N;
144 tmp = cmpmask + N;
145 off = mpoly_gen_monomial_offset_mp(one, var, bits, ctx->minfo);
146 mpoly_get_cmpmask(cmpmask, N, bits, ctx->minfo);
147
148 Alen = 0;
149 for (i = 0; i < Blen; i++)
150 {
151 fmpz_set_ui_array(k, Bexps + N*i + off, bits/FLINT_BITS);
152 n_fq_pow_cache_mulpow_fmpz(Acoeffs + d*Alen, Bcoeffs + d*i, k,
153 cache[0], cache[1], cache[2], ctx->fqctx);
154 if (_n_fq_is_zero(Acoeffs + d*Alen, d))
155 continue;
156
157 mpoly_monomial_mul_fmpz(tmp, one, N, k);
158 mpoly_monomial_sub_mp(Aexps + N*Alen, Bexps + N*i, tmp, N);
159 if (Alen < 1)
160 {
161 Alen = 1;
162 continue;
163 }
164 cmp = mpoly_monomial_cmp(Aexps + N*(Alen - 1), Aexps + N*Alen, N, cmpmask);
165 if (cmp != 0)
166 {
167 need_sort |= (cmp < 0);
168 Alen++;
169 continue;
170 }
171 _n_fq_add(Acoeffs + d*(Alen - 1), Acoeffs + d*(Alen - 1),
172 Acoeffs + d*Alen, d, fq_nmod_ctx_mod(ctx->fqctx));
173 Alen -= _n_fq_is_zero(Acoeffs + d*(Alen - 1), d);
174 }
175 A->length = Alen;
176
177 n_poly_stack_give_back(St, 3);
178 fmpz_clear(k);
179
180 TMP_END;
181
182 if (need_sort)
183 {
184 fq_nmod_mpoly_sort_terms(A, ctx);
185 fq_nmod_mpoly_combine_like_terms(A, ctx);
186 }
187
188 FLINT_ASSERT(fq_nmod_mpoly_is_canonical(A, ctx));
189 }
190
fq_nmod_mpoly_evaluate_one_fq_nmod(fq_nmod_mpoly_t A,const fq_nmod_mpoly_t B,slong var,const fq_nmod_t val,const fq_nmod_mpoly_ctx_t ctx)191 void fq_nmod_mpoly_evaluate_one_fq_nmod(
192 fq_nmod_mpoly_t A,
193 const fq_nmod_mpoly_t B,
194 slong var,
195 const fq_nmod_t val,
196 const fq_nmod_mpoly_ctx_t ctx)
197 {
198 n_poly_stack_t St;
199
200 if (fq_nmod_mpoly_is_zero(B, ctx))
201 {
202 fq_nmod_mpoly_zero(A, ctx);
203 return;
204 }
205
206 n_poly_stack_init(St);
207
208 if (B->bits <= FLINT_BITS)
209 _fq_nmod_mpoly_evaluate_one_fq_nmod_sp(A, B, var, val, ctx, St);
210 else
211 _fq_nmod_mpoly_evaluate_one_fq_nmod_mp(A, B, var, val, ctx, St);
212
213 n_poly_stack_clear(St);
214 }
215