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 /* exponents of B are not multiprecision */
_fmpq_mpoly_evaluate_one_fmpq_sp(fmpq_mpoly_t A,const fmpq_mpoly_t B,slong var,const fmpq_t val,const fmpq_mpoly_ctx_t ctx)15 int _fmpq_mpoly_evaluate_one_fmpq_sp(fmpq_mpoly_t A, const fmpq_mpoly_t B,
16 slong var, const fmpq_t val, const fmpq_mpoly_ctx_t ctx)
17 {
18 int success = 1;
19 int new;
20 slong i, j, N;
21 flint_bitcnt_t bits;
22 slong main_shift, main_off;
23 ulong main_exp, emin, emax;
24 ulong * cmpmask, * one;
25 slong Aalloc, Alen, Blen;
26 fmpz * Acoeff, * Bcoeff;
27 ulong * Aexp, * Bexp;
28 ulong * main_exps;
29 slong poweralloc;
30 fmpz * powers;
31 slong * inds;
32 ulong * exp_array, * exp;
33 slong next_loc;
34 slong heap_len;
35 mpoly_heap_s * heap;
36 mpoly_heap_t * x;
37 mpoly_heap_t * chain;
38 slong * store, * store_base;
39 mpoly_rbnode_struct ** stack;
40 slong stack_size;
41 mpoly_rbtree_t tree;
42 mpoly_rbnode_struct * node;
43 mpoly_rbnode_struct * root;
44 fmpq_t u;
45 TMP_INIT;
46
47 Blen = B->zpoly->length;
48 Bcoeff = B->zpoly->coeffs;
49 Bexp = B->zpoly->exps;
50 bits = B->zpoly->bits;
51
52 FLINT_ASSERT(A != B);
53 FLINT_ASSERT(bits <= FLINT_BITS);
54 FLINT_ASSERT(Blen > 0);
55
56 TMP_START;
57
58 fmpq_init(u);
59
60 N = mpoly_words_per_exp_sp(bits, ctx->zctx->minfo);
61 one = (ulong*) TMP_ALLOC(N*sizeof(ulong));
62 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
63 mpoly_gen_monomial_offset_shift_sp(one, &main_off, &main_shift,
64 var, bits, ctx->zctx->minfo);
65 mpoly_get_cmpmask(cmpmask, N, bits, ctx->zctx->minfo);
66
67 /* scan poly2 and put powers of var into tree */
68 /*
69 suppose x > y with lex order and poly2 is
70 x^3*y + x^2*y^3 + x*y + y + 1
71 and we are evaluating at y = 2
72
73 The data member of nodes in the tree of powers of y appearing in
74 the polynomial contain the term number of the first such appearance
75 data
76 y^0 4 (1 term)
77 y^1 0 (x^3*y term)
78 y^3 1 (x^2*y^3 term)
79
80 the inds array conains the term number of the next term with
81 the same power of y
82
83 inds[0] = 2 (x*y term)
84 inds[1] = -1
85 inds[2] = 3 (y term)
86 inds[3] = -1
87 inds[4] = -1
88 */
89 inds = (slong *) TMP_ALLOC(Blen*sizeof(slong));
90 mpoly_rbtree_init(tree);
91 emin = -UWORD(1);
92 emax = UWORD(0);
93
94 for (i = 0; i < Blen; i++)
95 {
96 ulong mask = (-UWORD(1)) >> (FLINT_BITS - bits);
97 main_exp = (Bexp[N*i + main_off] >> main_shift) & mask;
98 emin = FLINT_MIN(emin, main_exp);
99 emax = FLINT_MAX(emax, main_exp);
100 node = mpoly_rbtree_get(&new, tree, main_exp);
101 if (new)
102 {
103 node->data = (void *) i;
104 } else
105 {
106 inds[(slong) node->data2] = i;
107 }
108 node->data2 = (void *) i;
109 inds[i] = -WORD(1);
110 }
111 FLINT_ASSERT(emin <= emax);
112
113
114 success = success && !_fmpz_pow_ui_is_not_feasible(
115 fmpz_bits(fmpq_numref(val)), emin);
116 success = success && !_fmpz_pow_ui_is_not_feasible(
117 fmpz_bits(fmpq_denref(val)), emax);
118 if (success)
119 {
120 fmpz_pow_ui(fmpq_numref(u), fmpq_numref(val), emin);
121 fmpz_pow_ui(fmpq_denref(u), fmpq_denref(val), emax);
122 fmpq_mul(A->content, B->content, u);
123 }
124
125 /* manually traverse tree and add node data to heap */
126 heap_len = 1;
127 next_loc = tree->size + 4;
128 heap = (mpoly_heap_s *) TMP_ALLOC((tree->size + 1)*sizeof(mpoly_heap_s));
129 chain = (mpoly_heap_t *) TMP_ALLOC(tree->size*sizeof(mpoly_heap_t));
130 store = store_base = (slong *) TMP_ALLOC(2*tree->size*sizeof(slong));
131 exp_array = (ulong *) TMP_ALLOC(tree->size*N*sizeof(ulong));
132
133 i = 0;
134 stack_size = 0;
135 main_exps = (ulong *) TMP_ALLOC(tree->size*N*sizeof(ulong));
136 poweralloc = tree->size;
137 powers = _fmpz_vec_init(poweralloc);
138 stack = (mpoly_rbnode_struct **) TMP_ALLOC(tree->size
139 * sizeof(mpoly_rbnode_struct *));
140 root = tree->head->left;
141
142 looper:
143
144 while (root != tree->null)
145 {
146 stack[stack_size++] = root;
147 root = root->left;
148 }
149
150 if (stack_size == 0)
151 goto done;
152
153 root = stack[--stack_size];
154
155 FLINT_ASSERT(root->key >= emin);
156 FLINT_ASSERT(root->key <= emax);
157 mpoly_monomial_mul_ui(main_exps + N*i, one, N, root->key);
158
159 success = success && !_fmpz_pow_ui_is_not_feasible(
160 fmpz_bits(fmpq_numref(val)), root->key - emin);
161 success = success && !_fmpz_pow_ui_is_not_feasible(
162 fmpz_bits(fmpq_numref(val)), root->key - emin);
163 if (success)
164 {
165 fmpz_pow_ui(fmpq_numref(u), fmpq_numref(val), root->key - emin);
166 fmpz_pow_ui(fmpq_denref(u), fmpq_denref(val), emax - root->key);
167 fmpz_mul(powers + i, fmpq_numref(u), fmpq_denref(u));
168 }
169
170 x = chain + i;
171 x->i = i;
172 x->j = (slong) root->data;
173
174 x->next = NULL;
175 mpoly_monomial_sub(exp_array + N*i, Bexp + N*x->j, main_exps + N*i, N);
176 _mpoly_heap_insert(heap, exp_array + i*N, x,
177 &next_loc, &heap_len, N, cmpmask);
178
179 i++;
180 node = root->right;
181 flint_free(root);
182 root = node;
183
184 goto looper;
185
186 done:
187
188 FLINT_ASSERT(i == tree->size);
189
190 /* take from heap and put into A */
191 fmpz_mpoly_fit_bits(A->zpoly, bits, ctx->zctx);
192 A->zpoly->bits = bits;
193 Acoeff = A->zpoly->coeffs;
194 Aexp = A->zpoly->exps;
195 Aalloc = A->zpoly->alloc;
196 Alen = 0;
197 while (heap_len > 1)
198 {
199 exp = heap[1].exp;
200
201 _fmpz_mpoly_fit_length(&Acoeff, &Aexp, &Aalloc, Alen + 1, N);
202
203 mpoly_monomial_set(Aexp + Alen*N, exp, N);
204
205 fmpz_zero(Acoeff + Alen);
206 do {
207 x = _mpoly_heap_pop(heap, &heap_len, N, cmpmask);
208 do {
209 *store++ = x->i;
210 *store++ = x->j;
211 fmpz_addmul(Acoeff + Alen, powers + x->i, Bcoeff + x->j);
212 } while ((x = x->next) != NULL);
213 } while (heap_len > 1 && mpoly_monomial_equal(heap[1].exp, exp, N));
214
215 Alen += !fmpz_is_zero(Acoeff + Alen);
216
217 while (store > store_base)
218 {
219 j = *--store;
220 i = *--store;
221 if (inds[j] >= 0)
222 {
223 x = chain + i;
224 x->i = i;
225 x->j = inds[j];
226 x->next = NULL;
227 mpoly_monomial_sub(exp_array + N*i, Bexp + x->j*N,
228 main_exps + N*i, N);
229 _mpoly_heap_insert(heap, exp_array + i*N, x,
230 &next_loc, &heap_len, N, cmpmask);
231 }
232 }
233 }
234
235 A->zpoly->coeffs = Acoeff;
236 A->zpoly->exps = Aexp;
237 A->zpoly->alloc = Aalloc;
238 _fmpz_mpoly_set_length(A->zpoly, Alen, ctx->zctx);
239
240 fmpq_mpoly_reduce(A, ctx);
241
242 _fmpz_vec_clear(powers, poweralloc);
243
244 fmpq_clear(u);
245
246 TMP_END;
247
248 return success;
249 }
250
251 /* exponents of B are multiprecision */
_fmpq_mpoly_evaluate_one_fmpq_mp(fmpq_mpoly_t A,const fmpq_mpoly_t B,slong var,const fmpq_t val,const fmpq_mpoly_ctx_t ctx)252 int _fmpq_mpoly_evaluate_one_fmpq_mp(fmpq_mpoly_t A, const fmpq_mpoly_t B,
253 slong var, const fmpq_t val, const fmpq_mpoly_ctx_t ctx)
254 {
255 int success = 1;
256 int new;
257 slong i, j, N, bits;
258 slong main_off;
259 fmpz_t main_exp, emin, emax;
260 ulong * cmpmask, * main_one;
261 slong Aalloc, Alen, Blen;
262 fmpz * Acoeff, * Bcoeff;
263 ulong * Aexp, * Bexp;
264 ulong * main_exps;
265 slong poweralloc;
266 fmpz * powers;
267 slong * inds;
268 ulong * exp_array, * exp;
269 slong next_loc;
270 slong heap_len;
271 mpoly_heap_s * heap;
272 mpoly_heap_t * x;
273 mpoly_heap_t * chain;
274 slong * store, * store_base;
275 mpoly_rbnode_struct ** stack;
276 slong stack_size;
277 mpoly_rbtree_t tree;
278 mpoly_rbnode_struct * node;
279 mpoly_rbnode_struct * root;
280 fmpz_t t;
281 fmpq_t u;
282 TMP_INIT;
283
284 Blen = B->zpoly->length;
285 Bcoeff = B->zpoly->coeffs;
286 Bexp = B->zpoly->exps;
287 bits = B->zpoly->bits;
288
289 FLINT_ASSERT(A != B);
290 FLINT_ASSERT(bits > FLINT_BITS);
291 FLINT_ASSERT(Blen > 0);
292
293 TMP_START;
294
295 fmpz_init(t);
296 fmpq_init(u);
297 fmpz_init(emin);
298 fmpz_init(emax);
299
300 N = mpoly_words_per_exp(bits, ctx->zctx->minfo);
301 main_one = (ulong*) TMP_ALLOC(N*sizeof(ulong));
302 cmpmask = (ulong*) TMP_ALLOC(N*sizeof(ulong));
303 main_off = mpoly_gen_monomial_offset_mp(main_one, var, bits, ctx->zctx->minfo);
304 mpoly_get_cmpmask(cmpmask, N, bits, ctx->zctx->minfo);
305
306 /* scan poly2 and put powers of var into tree */
307 /*
308 suppose x > y with lex order and poly2 is
309 x^3*y + x^2*y^3 + x*y + y + 1
310 and we are evaluating at y = 2
311
312 The data member of nodes in the tree of powers of y appearing in
313 the polynomial contain the term number of the first such appearance
314 data
315 y^0 4 (1 term)
316 y^1 0 (x^3*y term)
317 y^3 1 (x^2*y^3 term)
318
319 the inds array conains the term number of the next term with
320 the same power of y
321
322 inds[0] = 2 (x*y term)
323 inds[1] = -1
324 inds[2] = 3 (y term)
325 inds[3] = -1
326 inds[4] = -1
327 */
328 inds = (slong *) TMP_ALLOC(Blen*sizeof(slong));
329 mpoly_rbtree_init(tree);
330 fmpz_init(main_exp);
331 for (i = 0; i < Blen; i++)
332 {
333 fmpz_set_ui_array(main_exp, Bexp + N*i + main_off, bits/FLINT_BITS);
334 node = mpoly_rbtree_get_fmpz(&new, tree, main_exp);
335 if (i == 0)
336 {
337 fmpz_set(emin, main_exp);
338 fmpz_set(emax, main_exp);
339 }
340 else
341 {
342 if (fmpz_cmp(emin, main_exp) > 0)
343 fmpz_set(emin, main_exp);
344 if (fmpz_cmp(emax, main_exp) < 0)
345 fmpz_set(emax, main_exp);
346 }
347 if (new)
348 {
349 node->data = (void *) i;
350 } else
351 {
352 inds[(slong) node->data2] = i;
353 }
354 node->data2 = (void *) i;
355 inds[i] = -WORD(1);
356 }
357 fmpz_clear(main_exp);
358 FLINT_ASSERT(fmpz_cmp(emax, emin) >= 0);
359
360 success = success && fmpz_pow_fmpz(fmpq_numref(u), fmpq_numref(val), emin);
361 success = success && fmpz_pow_fmpz(fmpq_denref(u), fmpq_denref(val), emax);
362 fmpq_mul(A->content, B->content, u);
363
364 /* manually traverse tree and add node data to heap */
365 heap_len = 1;
366 next_loc = tree->size + 4;
367 heap = (mpoly_heap_s *) TMP_ALLOC((tree->size + 1)*sizeof(mpoly_heap_s));
368 chain = (mpoly_heap_t *) TMP_ALLOC(tree->size*sizeof(mpoly_heap_t));
369 store = store_base = (slong *) TMP_ALLOC(2*tree->size*sizeof(slong));
370 exp_array = (ulong *) TMP_ALLOC(tree->size*N*sizeof(ulong));
371
372 i = 0;
373 stack_size = 0;
374 main_exps = (ulong *) TMP_ALLOC(tree->size*N*sizeof(ulong));
375 poweralloc = tree->size;
376 powers = _fmpz_vec_init(poweralloc);
377 stack = (mpoly_rbnode_struct **) TMP_ALLOC(tree->size
378 * sizeof(mpoly_rbnode_struct *));
379 root = tree->head->left;
380
381 looper:
382
383 while (root != tree->null)
384 {
385 stack[stack_size++] = root;
386 root = root->left;
387 }
388
389 if (stack_size == 0)
390 goto done;
391
392 root = stack[--stack_size];
393
394 mpoly_monomial_mul_fmpz(main_exps + N*i, main_one, N, &root->key);
395
396 fmpz_sub(t, &root->key, emin);
397 success = success && fmpz_pow_fmpz(fmpq_numref(u), fmpq_numref(val), t);
398 fmpz_sub(t, emax, &root->key);
399 success = success && fmpz_pow_fmpz(fmpq_denref(u), fmpq_denref(val), t);
400 fmpz_mul(powers + i, fmpq_numref(u), fmpq_denref(u));
401
402 x = chain + i;
403 x->i = i;
404 x->j = (slong) root->data;
405
406 x->next = NULL;
407 mpoly_monomial_sub_mp(exp_array + N*i, Bexp + N*x->j, main_exps + N*i, N);
408 _mpoly_heap_insert(heap, exp_array + i*N, x,
409 &next_loc, &heap_len, N, cmpmask);
410
411 i++;
412 node = root->right;
413 fmpz_clear(&root->key);
414 flint_free(root);
415 root = node;
416
417 goto looper;
418
419 done:
420
421 FLINT_ASSERT(i == tree->size);
422
423 /* take from heap and put into A */
424 fmpz_mpoly_fit_bits(A->zpoly, bits, ctx->zctx);
425 A->zpoly->bits = bits;
426 Acoeff = A->zpoly->coeffs;
427 Aexp = A->zpoly->exps;
428 Aalloc = A->zpoly->alloc;
429 Alen = 0;
430 while (heap_len > 1)
431 {
432 exp = heap[1].exp;
433
434 _fmpz_mpoly_fit_length(&Acoeff, &Aexp, &Aalloc, Alen + 1, N);
435
436 mpoly_monomial_set(Aexp + Alen*N, exp, N);
437
438 fmpz_zero(Acoeff + Alen);
439 do {
440 x = _mpoly_heap_pop(heap, &heap_len, N, cmpmask);
441 do {
442 *store++ = x->i;
443 *store++ = x->j;
444 fmpz_addmul(Acoeff + Alen, powers + x->i, Bcoeff + x->j);
445 } while ((x = x->next) != NULL);
446 } while (heap_len > 1 && mpoly_monomial_equal(heap[1].exp, exp, N));
447
448 Alen += !fmpz_is_zero(Acoeff + Alen);
449
450 while (store > store_base)
451 {
452 j = *--store;
453 i = *--store;
454 if (inds[j] >= 0)
455 {
456 x = chain + i;
457 x->i = i;
458 x->j = inds[j];
459 x->next = NULL;
460 mpoly_monomial_sub_mp(exp_array + N*i, Bexp + N*x->j,
461 main_exps + N*i, N);
462 _mpoly_heap_insert(heap, exp_array + N*i, x,
463 &next_loc, &heap_len, N, cmpmask);
464 }
465 }
466 }
467
468 A->zpoly->coeffs = Acoeff;
469 A->zpoly->exps = Aexp;
470 A->zpoly->alloc = Aalloc;
471 _fmpz_mpoly_set_length(A->zpoly, Alen, ctx->zctx);
472
473 _fmpz_vec_clear(powers, poweralloc);
474 fmpq_clear(u);
475 fmpz_clear(t);
476 fmpz_clear(emin);
477 fmpz_clear(emax);
478
479 fmpq_mpoly_reduce(A, ctx);
480
481 TMP_END;
482
483 return success;
484 }
485
fmpq_mpoly_evaluate_one_fmpq(fmpq_mpoly_t A,const fmpq_mpoly_t B,slong var,const fmpq_t val,const fmpq_mpoly_ctx_t ctx)486 int fmpq_mpoly_evaluate_one_fmpq(fmpq_mpoly_t A,
487 const fmpq_mpoly_t B, slong var, const fmpq_t val,
488 const fmpq_mpoly_ctx_t ctx)
489 {
490 if (B->zpoly->length == 0)
491 {
492 fmpq_mpoly_zero(A, ctx);
493 return 1;
494 }
495
496 if (A == B)
497 {
498 int success;
499 fmpq_mpoly_t T;
500 fmpq_mpoly_init(T, ctx);
501 success = fmpq_mpoly_evaluate_one_fmpq(T, B, var, val, ctx);
502 fmpq_mpoly_swap(A, T, ctx);
503 fmpq_mpoly_clear(T, ctx);
504 return success;
505 }
506
507 if (B->zpoly->bits <= FLINT_BITS)
508 {
509 return _fmpq_mpoly_evaluate_one_fmpq_sp(A, B, var, val, ctx);
510 }
511 else
512 {
513 return _fmpq_mpoly_evaluate_one_fmpq_mp(A, B, var, val, ctx);
514 }
515 }
516
517