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 "fmpq_mpoly.h"
13
14 /*
15 Given a polynomial tree with exponents stored in the keys and
16 coefficients stored in the data member,
17 the function _mpoly_rbnode_clear_evalall_tree_fmpz clears the tree
18 and stores the evaluation of the polynomial in "l".
19
20 poly = a0 x^0 + a1 x^1 + a2 x^2 + a3 x^3 + a4 x^4 + a5 x^5 + a6 x^6
21
22 tree = a3 x^3
23 a1 x^1 a5 x^5
24
25 a0 x^0 a2 x^2 a4 x^4 a6 x^6
26
27 l = a0*x^0 + x^1*(a1 + a2*x^1) + x^3*(a3 + a4*x^1 + x^2*(a5 + a6*x^1))
28
29 */
_mpoly_rbnode_clear_sp(mpoly_rbtree_t tree,mpoly_rbnode_t node,slong s,fmpq_t l,const fmpq_t x)30 static void _mpoly_rbnode_clear_sp(mpoly_rbtree_t tree, mpoly_rbnode_t node,
31 slong s, fmpq_t l, const fmpq_t x)
32 {
33 fmpq_t r, xp;
34 slong e = node->key;
35 FLINT_ASSERT(e >= s);
36
37 fmpq_init(r);
38 fmpq_zero(r);
39 if (node->right != tree->null)
40 _mpoly_rbnode_clear_sp(tree, node->right, e, r, x);
41
42 fmpq_zero(l);
43 if (node->left != tree->null)
44 _mpoly_rbnode_clear_sp(tree, node->left, s, l, x);
45
46 fmpq_init(xp);
47 fmpq_pow_si(xp, x, e - s);
48 fmpq_add(r, r, (fmpq*)(&node->data));
49 fmpq_addmul(l, xp, r);
50
51 fmpq_clear(r);
52 fmpq_clear(xp);
53 fmpq_clear((fmpq*)(&node->data));
54 flint_free(node);
55 }
56
57
58 /*
59 evaluate a f(xbar) at xbar = val,
60 */
_fmpz_mpoly_evaluate_all_tree_fmpq_sp(fmpq_t ev,const fmpz_mpoly_t poly,fmpq * const * vals,const fmpz_mpoly_ctx_t ctx)61 int _fmpz_mpoly_evaluate_all_tree_fmpq_sp(fmpq_t ev, const fmpz_mpoly_t poly,
62 fmpq * const * vals, const fmpz_mpoly_ctx_t ctx)
63 {
64 int success = 1;
65 int new;
66 slong i, j, k, N, bits, nvars = ctx->minfo->nvars;
67 slong main_exp, main_var, main_shift, main_off, shift, off;
68 ulong mask;
69 slong entries, k_len;
70 slong p_len;
71 fmpz * p_coeff;
72 ulong * p_exp;
73 slong * degrees;
74 slong * offs;
75 ulong * masks;
76 fmpq * powers;
77 fmpq_t t;
78 mpoly_rbtree_t tree;
79 mpoly_rbnode_struct * node;
80 TMP_INIT;
81
82 p_len = poly->length;
83 p_coeff = poly->coeffs;
84 p_exp = poly->exps;
85 bits = poly->bits;
86
87 FLINT_ASSERT(p_len > 0);
88
89 TMP_START;
90
91 degrees = (slong *) TMP_ALLOC(nvars*sizeof(slong));
92 fmpz_mpoly_degrees_si(degrees, poly, ctx);
93
94 /* pick main variable with highest degree */
95 main_var = 0;
96 for (i = 1; i < nvars; i++)
97 {
98 if (degrees[i] > degrees[main_var])
99 main_var = i;
100 }
101
102 /* compute how many masks are needed */
103 entries = 0;
104 for (i = 0; i < nvars; i++)
105 {
106 if (_fmpz_pow_ui_is_not_feasible(fmpq_height_bits(vals[i]), degrees[i]))
107 {
108 success = 0;
109 goto cleanup_degrees;
110 }
111
112 if (i == main_var)
113 continue;
114
115 entries += FLINT_BIT_COUNT(degrees[i]);
116 }
117 offs = (slong *) TMP_ALLOC(entries*sizeof(slong));
118 masks = (ulong *) TMP_ALLOC(entries*sizeof(slong));
119 powers = (fmpq *) TMP_ALLOC(entries*sizeof(fmpq));
120
121 N = mpoly_words_per_exp(bits, ctx->minfo);
122
123 /* store bit masks for each power of two of the non-main variables */
124 k = 0;
125 for (i = 0; i < nvars; i++)
126 {
127 flint_bitcnt_t varibits;
128
129 if (i == main_var)
130 continue;
131
132 mpoly_gen_offset_shift_sp(&off, &shift, i, bits, ctx->minfo);
133 varibits = FLINT_BIT_COUNT(degrees[i]);
134 for (j = 0; j < varibits; j++)
135 {
136 offs[k] = off;
137 masks[k] = UWORD(1) << (shift + j);
138 fmpq_init(powers + k);
139 if (j == 0)
140 fmpq_set(powers + k, vals[i]);
141 else
142 fmpq_mul(powers + k, powers + k - 1, powers + k - 1);
143 k++;
144 }
145 }
146 k_len = k;
147 FLINT_ASSERT(k_len == entries);
148
149 /* accumulate coefficients of the main variable */
150 mpoly_gen_offset_shift_sp(&main_off, &main_shift, main_var, bits, ctx->minfo);
151 mpoly_rbtree_init(tree);
152 fmpq_init(t);
153 mask = (-UWORD(1)) >> (FLINT_BITS - bits);
154 for (i = 0; i < p_len; i++)
155 {
156 main_exp = (p_exp[N*i + main_off] >> main_shift) & mask;
157 node = mpoly_rbtree_get(&new, tree, main_exp);
158 if (new)
159 {
160 fmpq_init((fmpq*)(&node->data));
161 fmpq_zero((fmpq*)(&node->data));
162 }
163
164 fmpz_set(fmpq_numref(t), p_coeff + i);
165 fmpz_one(fmpq_denref(t));
166 for (k = 0; k < k_len; k++)
167 {
168 if ((p_exp[N*i + offs[k]] & masks[k]) != WORD(0))
169 fmpq_mul(t, t, powers + k);
170 }
171 fmpq_add((fmpq*)(&node->data), (fmpq*)(&node->data), t);
172 }
173 fmpq_clear(t);
174
175 for (k = 0; k < k_len; k++)
176 fmpq_clear(powers + k);
177
178 /* use tree method to evaluate in the main variable */
179 _mpoly_rbnode_clear_sp(tree, tree->head->left, WORD(0), ev, vals[main_var]);
180
181 cleanup_degrees:
182
183 TMP_END;
184
185 return success;
186 }
187
188
189
190 /*
191 Given a polynomial tree with exponents stored in the keys and
192 coefficients stored in the data member,
193 the function _mpoly_rbnode_clear_evalall_tree_fmpz clears the tree
194 and stores the evaluation of the polynomial in "l".
195
196 poly = a0 x^0 + a1 x^1 + a2 x^2 + a3 x^3 + a4 x^4 + a5 x^5 + a6 x^6
197
198 tree = a3 x^3
199 a1 x^1 a5 x^5
200
201 a0 x^0 a2 x^2 a4 x^4 a6 x^6
202
203 l = a0*x^0 + x^1*(a1 + a2*x^1) + x^3*(a3 + a4*x^1 + x^2*(a5 + a6*x^1))
204
205 */
_mpoly_rbnode_clear_mp(mpoly_rbtree_t tree,mpoly_rbnode_t node,const fmpz_t s,fmpq_t l,const fmpq_t x)206 static int _mpoly_rbnode_clear_mp(mpoly_rbtree_t tree, mpoly_rbnode_t node,
207 const fmpz_t s, fmpq_t l, const fmpq_t x)
208 {
209 int success = 1;
210 fmpq_t r, xp;
211 FLINT_ASSERT(fmpz_cmp(&node->key, s) >= 0);
212
213 fmpq_init(r);
214 if (node->right != tree->null)
215 {
216 if (!_mpoly_rbnode_clear_mp(tree, node->right, &node->key, r, x))
217 success = 0;
218 }
219
220 fmpq_zero(l);
221 if (node->left != tree->null)
222 {
223 if (!_mpoly_rbnode_clear_mp(tree, node->left, s, l, x))
224 success = 0;
225 }
226
227 fmpq_init(xp);
228
229 fmpz_sub((fmpz*)(&node->key), (fmpz*)(&node->key), s);
230 if (!fmpq_pow_fmpz(xp, x, (fmpz*)(&node->key)))
231 success = 0;
232 fmpq_add(r, r, (fmpq*)(&node->data));
233 fmpq_addmul(l, xp, r);
234
235 fmpq_clear(r);
236 fmpq_clear(xp);
237 fmpq_clear((fmpq*)(&node->data));
238 fmpz_clear((fmpz*)(&node->key));
239 flint_free(node);
240
241 return success;
242 }
243
244
245 /*
246 evaluate a f(xbar) at xbar = val,
247 */
_fmpz_mpoly_evaluate_all_tree_fmpq_mp(fmpq_t ev,const fmpz_mpoly_t poly,fmpq * const * vals,const fmpz_mpoly_ctx_t ctx)248 int _fmpz_mpoly_evaluate_all_tree_fmpq_mp(fmpq_t ev, const fmpz_mpoly_t poly,
249 fmpq * const * vals, const fmpz_mpoly_ctx_t ctx)
250 {
251 int success = 1;
252 int new;
253 slong i, j, k, N, bits, nvars = ctx->minfo->nvars;
254 slong main_var, main_off, off;
255 fmpz_t main_exp;
256 slong entries, k_len;
257 slong p_len;
258 fmpz * p_coeff;
259 ulong * p_exp;
260 slong * degrees;
261 slong * offs;
262 ulong * masks;
263 fmpq * powers;
264 fmpq_t t;
265 fmpz_t s;
266 mpoly_rbtree_t tree;
267 mpoly_rbnode_struct * node;
268 TMP_INIT;
269
270 p_len = poly->length;
271 p_coeff = poly->coeffs;
272 p_exp = poly->exps;
273 bits = poly->bits;
274
275 FLINT_ASSERT(p_len > 0);
276
277 TMP_START;
278
279 degrees = _fmpz_vec_init(nvars);
280 mpoly_degrees_ffmpz(degrees, p_exp, p_len, bits, ctx->minfo);
281
282 /* pick main variable with highest degree */
283 main_var = 0;
284 for (i = 1; i < nvars; i++)
285 {
286 if (fmpz_cmp(degrees + i, degrees + main_var) > 0)
287 main_var = i;
288 }
289
290 /* compute how many masks are needed */
291 entries = 0;
292 for (i = 0; i < nvars; i++)
293 {
294 if (_fmpz_pow_fmpz_is_not_feasible(fmpq_height_bits(vals[i]), degrees + i))
295 {
296 success = 0;
297 goto cleanup_degrees;
298 }
299
300 if (i == main_var)
301 continue;
302
303 entries += fmpz_bits(degrees + i);
304 }
305 offs = (slong *) TMP_ALLOC(entries*sizeof(slong));
306 masks = (ulong *) TMP_ALLOC(entries*sizeof(slong));
307 powers = (fmpq *) TMP_ALLOC(entries*sizeof(fmpq));
308
309 N = mpoly_words_per_exp(bits, ctx->minfo);
310
311 /* store bit masks for each power of two of the non-main variables */
312 k = 0;
313 for (i = 0; i < nvars; i++)
314 {
315 flint_bitcnt_t varibits;
316
317 if (i == main_var)
318 continue;
319
320 off = mpoly_gen_offset_mp(i, bits, ctx->minfo);
321 varibits = fmpz_bits(degrees + i);
322 if (varibits >= FLINT_BITS && !fmpz_is_one(fmpq_denref(vals[i]))
323 && (fmpz_is_zero(fmpq_denref(vals[i]))
324 ))
325 {
326 flint_throw(FLINT_ERROR, "Exponent too large in fmpq_mpoly_evaluate_one_fmpq");
327 }
328 for (j = 0; j < varibits; j++)
329 {
330 offs[k] = off + (j / FLINT_BITS);
331 masks[k] = UWORD(1) << (j % FLINT_BITS);
332 fmpq_init(powers + k);
333 if (j == 0)
334 fmpq_set(powers + k, vals[i]);
335 else
336 fmpq_mul(powers + k, powers + k - 1, powers + k - 1);
337 k++;
338 }
339 }
340 k_len = k;
341 FLINT_ASSERT(k_len == entries);
342
343 /* accumulate coefficients of the main variable */
344 main_off = mpoly_gen_offset_mp(main_var, bits, ctx->minfo);
345 mpoly_rbtree_init(tree);
346 fmpq_init(t);
347 fmpz_init(main_exp);
348 for (i = 0; i < p_len; i++)
349 {
350 fmpz_set_ui_array(main_exp, p_exp + N*i + main_off, bits/FLINT_BITS);
351 node = mpoly_rbtree_get_fmpz(&new, tree, main_exp);
352 if (new)
353 {
354 fmpq_init((fmpq*)(&node->data));
355 fmpq_zero((fmpq*)(&node->data));
356 }
357
358 fmpz_set(fmpq_numref(t), p_coeff + i);
359 fmpz_one(fmpq_denref(t));
360 for (k = 0; k < k_len; k++)
361 {
362 if ((p_exp[N*i + offs[k]] & masks[k]) != WORD(0))
363 fmpq_mul(t, t, powers + k);
364 }
365 fmpq_add((fmpq*)(&node->data), (fmpq*)(&node->data), t);
366 }
367 fmpz_clear(main_exp);
368 fmpq_clear(t);
369
370 for (k = 0; k < k_len; k++)
371 fmpq_clear(powers + k);
372
373 /* use tree method to evaluate in the main variable */
374 fmpz_init(s);
375 if (!_mpoly_rbnode_clear_mp(tree, tree->head->left, s, ev, vals[main_var]))
376 success = 0;
377 fmpz_clear(s);
378
379 cleanup_degrees:
380
381 _fmpz_vec_clear(degrees, nvars);
382
383 TMP_END;
384
385 return success;
386 }
387
388
fmpq_mpoly_evaluate_all_fmpq(fmpq_t ev,const fmpq_mpoly_t A,fmpq * const * vals,const fmpq_mpoly_ctx_t ctx)389 int fmpq_mpoly_evaluate_all_fmpq(fmpq_t ev, const fmpq_mpoly_t A,
390 fmpq * const * vals, const fmpq_mpoly_ctx_t ctx)
391 {
392 int success;
393 fmpq_t t;
394
395 if (fmpq_mpoly_is_zero(A, ctx))
396 {
397 fmpq_zero(ev);
398 return 1;
399 }
400
401 fmpq_init(t);
402
403 if (A->zpoly->bits <= FLINT_BITS)
404 {
405 success = _fmpz_mpoly_evaluate_all_tree_fmpq_sp(t, A->zpoly, vals, ctx->zctx);
406 }
407 else
408 {
409 success = _fmpz_mpoly_evaluate_all_tree_fmpq_mp(t, A->zpoly, vals, ctx->zctx);
410 }
411
412 if (success)
413 fmpq_mul(ev, t, A->content);
414
415 fmpq_clear(t);
416 return success;
417 }
418