1 /**
2 * Copyright 2015, SRI International.
3 *
4 * This file is part of LibPoly.
5 *
6 * LibPoly is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU Lesser General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * LibPoly is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public License
17 * along with LibPoly. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 #include <interval.h>
21 #include <upolynomial.h>
22 #include <variable_db.h>
23 #include <monomial.h>
24 #include <variable_list.h>
25
26 #include "polynomial/polynomial.h"
27 #include "polynomial/coefficient.h"
28 #include "polynomial/output.h"
29 #include "polynomial/gcd.h"
30
31 #include "upolynomial/upolynomial.h"
32
33 #include "number/integer.h"
34 #include "number/rational.h"
35 #include "number/value.h"
36 #include "interval/arithmetic.h"
37
38 #include "variable/variable_order.h"
39 #include "polynomial/polynomial_context.h"
40 #include "polynomial/polynomial_vector.h"
41
42 #include <assignment.h>
43
44 #include "utils/debug_trace.h"
45 #include "utils/statistics.h"
46
47 #include <assert.h>
48
49 static
50 void coefficient_resolve_algebraic(const lp_polynomial_context_t* ctx, const coefficient_t* A, const lp_assignment_t* m, coefficient_t* A_alg);
51
52 static
algebraic_interval_remember(const lp_variable_list_t * var_list,const lp_assignment_t * m)53 lp_dyadic_interval_t* algebraic_interval_remember(const lp_variable_list_t* var_list, const lp_assignment_t* m) {
54 // Remember the existing intervals so that we can recover later
55 lp_dyadic_interval_t* x_i_intervals = malloc(sizeof(lp_dyadic_interval_t)*var_list->list_size);
56 size_t i;
57 for (i = 0; i < var_list->list_size; ++ i) {
58 lp_variable_t x_i = var_list->list[i];
59 const lp_value_t* x_i_value = lp_assignment_get_value(m, x_i);
60 if (!lp_value_is_rational(x_i_value)) {
61 lp_dyadic_interval_construct_copy(x_i_intervals + i, &x_i_value->value.a.I);
62 } else {
63 lp_dyadic_interval_construct_zero(x_i_intervals + i);
64 }
65 }
66 return x_i_intervals;
67 }
68
69 static
algebraic_interval_restore(const lp_variable_list_t * var_list,lp_dyadic_interval_t * cache,const lp_assignment_t * m)70 void algebraic_interval_restore(const lp_variable_list_t* var_list, lp_dyadic_interval_t* cache, const lp_assignment_t* m) {
71 // Restore the intervals, and destroy the temps
72 size_t i;
73 for (i = 0; i < var_list->list_size; ++ i) {
74 lp_variable_t x_i = var_list->list[i];
75 const lp_value_t* x_i_value = lp_assignment_get_value(m, x_i);
76 if (x_i_value->type == LP_VALUE_ALGEBRAIC && !x_i_value->value.a.I.is_point) {
77 lp_algebraic_number_restore_interval_const(&x_i_value->value.a, cache + i);
78 }
79 lp_dyadic_interval_destruct(cache + i);
80 }
81 free(cache);
82 }
83
84
85 /**
86 * Make sure that the coefficient has the given capacity for the given variable.
87 * If the polynomial a constant or a polynomial in a smaller variable, we
88 * transform it into a polynomial and adjust the coefficent.
89 */
90 static void
91 coefficient_ensure_capacity(const lp_polynomial_context_t* ctx, coefficient_t* C, lp_variable_t x, size_t capacity);
92
93 /**
94 * Assumes:
95 * * Sub-coefficients have been normalized.
96 *
97 * Normalize the coefficient:
98 * * If the coefficient is a polynomial with only the constant coefficient, it
99 * should be upgraded
100 * * The highest coefficient should be non-zero
101 */
102 static void
103 coefficient_normalize(const lp_polynomial_context_t* ctx, coefficient_t* C);
104
105 /**
106 * Same as above, but uses the model.
107 */
108 static void
109 coefficient_normalize_m(const lp_polynomial_context_t* ctx, coefficient_t* C, const lp_assignment_t* m);
110
111 int
112 coefficient_is_normalized(const lp_polynomial_context_t* ctx, coefficient_t* C);
113
STAT_DECLARE(int,coefficient,construct)114 STAT_DECLARE(int, coefficient, construct)
115
116 void coefficient_construct(const lp_polynomial_context_t* ctx, coefficient_t* C) {
117 TRACE("coefficient::internal", "coefficient_construct()\n");
118 STAT_INCR(coefficient, construct)
119
120 C->type = COEFFICIENT_NUMERIC;
121 integer_construct_from_int(ctx->K, &C->value.num, 0);
122 }
123
STAT_DECLARE(int,coefficient,construct_from_int)124 STAT_DECLARE(int, coefficient, construct_from_int)
125
126 void coefficient_construct_from_int(const lp_polynomial_context_t* ctx, coefficient_t* C, long C_int) {
127 TRACE("coefficient::internal", "coefficient_construct_from_int()\n");
128 STAT_INCR(coefficient, construct_from_int)
129
130 C->type = COEFFICIENT_NUMERIC;
131 integer_construct_from_int(ctx->K, &C->value.num, C_int);
132 }
133
STAT_DECLARE(int,coefficient,construct_from_integer)134 STAT_DECLARE(int, coefficient, construct_from_integer)
135
136 void coefficient_construct_from_integer(const lp_polynomial_context_t* ctx, coefficient_t* C, const lp_integer_t* C_integer) {
137 TRACE("coefficient::internal", "coefficient_construct_from_integer()\n");
138 STAT_INCR(coefficient, construct_from_integer)
139
140 C->type = COEFFICIENT_NUMERIC;
141 integer_construct_copy(ctx->K, &C->value.num, C_integer);
142 }
143
STAT_DECLARE(int,coefficient,construct_rec)144 STAT_DECLARE(int, coefficient, construct_rec)
145
146 void coefficient_construct_rec(const lp_polynomial_context_t* ctx, coefficient_t* C, lp_variable_t x, size_t capacity) {
147 TRACE("coefficient::internal", "coefficient_construct_rec()\n");
148 STAT_INCR(coefficient, construct_rec)
149
150 C->type = COEFFICIENT_POLYNOMIAL;
151 C->value.rec.x = x;
152 C->value.rec.size = 0;
153 C->value.rec.capacity = 0;
154 C->value.rec.coefficients = 0;
155 coefficient_ensure_capacity(ctx, C, x, capacity);
156 }
157
STAT_DECLARE(int,coefficient,construct_simple_int)158 STAT_DECLARE(int, coefficient, construct_simple_int)
159
160 void coefficient_construct_simple_int(const lp_polynomial_context_t* ctx, coefficient_t* C, long a, lp_variable_t x, unsigned n) {
161 TRACE("coefficient::internal", "coefficient_construct_simple_int()\n");
162 STAT_INCR(coefficient, construct_simple_int)
163
164 if (n == 0) {
165 coefficient_construct_from_int(ctx, C, a);
166 } else {
167 // x^n
168 coefficient_construct_rec(ctx, C, x, n+1);
169 integer_assign_int(ctx->K, &COEFF(C, n)->value.num, a);
170 }
171 }
172
STAT_DECLARE(int,coefficient,construct_simple)173 STAT_DECLARE(int, coefficient, construct_simple)
174
175 void coefficient_construct_simple(const lp_polynomial_context_t* ctx, coefficient_t* C, const lp_integer_t* a, lp_variable_t x, unsigned n) {
176 TRACE("coefficient::internal", "coefficient_construct_simple()\n");
177 STAT_INCR(coefficient, construct_simple)
178
179 if (n == 0) {
180 coefficient_construct_from_integer(ctx, C, a);
181 } else {
182 // x^n
183 coefficient_construct_rec(ctx, C, x, n+1);
184 integer_assign(ctx->K, &COEFF(C, n)->value.num, a);
185 }
186 }
187
STAT_DECLARE(int,coefficient,construt_linear)188 STAT_DECLARE(int, coefficient, construt_linear)
189
190 void coefficient_construct_linear(const lp_polynomial_context_t* ctx, coefficient_t* C, const lp_integer_t* a, const lp_integer_t* b, lp_variable_t x) {
191 TRACE("coefficient::internal", "coefficient_construct_simple()\n");
192 STAT_INCR(coefficient, coefficient_construct_linear)
193
194 assert(integer_sgn(lp_Z, a) != 0);
195
196 // a*x + b
197 coefficient_construct_rec(ctx, C, x, 2);
198 integer_assign(ctx->K, &COEFF(C, 1)->value.num, a);
199 integer_assign(ctx->K, &COEFF(C, 0)->value.num, b);
200 }
201
STAT_DECLARE(int,coefficient,construct_copy)202 STAT_DECLARE(int, coefficient, construct_copy)
203
204 void coefficient_construct_copy(const lp_polynomial_context_t* ctx, coefficient_t* C, const coefficient_t* from) {
205 TRACE("coefficient::internal", "coefficient_construct_copy()\n");
206 STAT_INCR(coefficient, construct_copy)
207
208 size_t i;
209 switch(from->type) {
210 case COEFFICIENT_NUMERIC:
211 C->type = COEFFICIENT_NUMERIC;
212 integer_construct_copy(ctx->K, &C->value.num, &from->value.num);
213 break;
214 case COEFFICIENT_POLYNOMIAL:
215 C->type = COEFFICIENT_POLYNOMIAL;
216 C->value.rec.x = VAR(from);
217 C->value.rec.size = SIZE(from);
218 C->value.rec.capacity = SIZE(from);
219 C->value.rec.coefficients = malloc(SIZE(from) * sizeof(coefficient_t));
220 for (i = 0; i < SIZE(from); ++ i) {
221 coefficient_construct_copy(ctx, COEFF(C, i), COEFF(from, i));
222 }
223 break;
224 }
225 }
226
STAT_DECLARE(int,coefficient,construct_from_univariate)227 STAT_DECLARE(int, coefficient, construct_from_univariate)
228
229 void coefficient_construct_from_univariate(const lp_polynomial_context_t* ctx,
230 coefficient_t* C, const lp_upolynomial_t* C_u, lp_variable_t x) {
231
232 TRACE("coefficient::internal", "coefficient_construct_from_univariate()\n");
233 STAT_INCR(coefficient, construct_from_univariate)
234
235 // Get the coefficients
236 size_t C_u_deg = lp_upolynomial_degree(C_u);
237 lp_integer_t* coeff = malloc(sizeof(lp_integer_t)*(C_u_deg + 1));
238
239 size_t i;
240 for (i = 0; i <= C_u_deg; ++ i) {
241 integer_construct_from_int(ctx->K, coeff + i, 0);
242 }
243
244 lp_upolynomial_unpack(C_u, coeff);
245
246 // Construct the polynomial
247 coefficient_construct_rec(ctx, C, x, C_u_deg + 1);
248
249 // Move over the coefficients
250 for (i = 0; i <= C_u_deg; ++ i) {
251 integer_swap(&COEFF(C, i)->value.num, coeff + i);
252 integer_destruct(coeff + i);
253 }
254 free(coeff);
255
256 // Normalize (it might be constant)
257 coefficient_normalize(ctx, C);
258
259 assert(coefficient_is_normalized(ctx, C));
260 }
261
coefficient_destruct(coefficient_t * C)262 void coefficient_destruct(coefficient_t* C) {
263 TRACE("coefficient::internal", "coefficient_destruct()\n");
264
265 size_t i;
266 switch (C->type) {
267 case COEFFICIENT_NUMERIC:
268 integer_destruct(&C->value.num);
269 break;
270 case COEFFICIENT_POLYNOMIAL:
271 for (i = 0; i < CAPACITY(C); ++ i) {
272 coefficient_destruct(COEFF(C, i));
273 }
274 free(C->value.rec.coefficients);
275 break;
276 default:
277 assert(0);
278 }
279 }
280
STAT_DECLARE(int,coefficient,swap)281 STAT_DECLARE(int, coefficient, swap)
282
283 void coefficient_swap(coefficient_t* C1, coefficient_t* C2) {
284 TRACE("coefficient::internal", "coefficient_swap()\n");
285 STAT_INCR(coefficient, swap)
286 coefficient_t tmp = *C1;
287 *C1 = *C2;
288 *C2 = tmp;
289 }
290
STAT_DECLARE(int,coefficient,assign)291 STAT_DECLARE(int, coefficient, assign)
292
293 void coefficient_assign(const lp_polynomial_context_t* ctx, coefficient_t* C, const coefficient_t* from) {
294 TRACE("coefficient::internal", "coefficient_assign()\n");
295 STAT_INCR(coefficient, assign)
296
297 if (C != from) {
298 coefficient_t result;
299 switch(from->type) {
300 case COEFFICIENT_NUMERIC:
301 if (C->type == COEFFICIENT_POLYNOMIAL) {
302 coefficient_destruct(C);
303 coefficient_construct_copy(ctx, C, from);
304 } else {
305 integer_assign(ctx->K, &C->value.num, &from->value.num);
306 }
307 break;
308 case COEFFICIENT_POLYNOMIAL:
309 coefficient_construct_copy(ctx, &result, from);
310 coefficient_swap(&result, C);
311 coefficient_destruct(&result);
312 break;
313 }
314 }
315
316 assert(coefficient_is_normalized(ctx, C));
317 }
318
STAT_DECLARE(int,coefficient,assign_int)319 STAT_DECLARE(int, coefficient, assign_int)
320
321 void coefficient_assign_int(const lp_polynomial_context_t* ctx, coefficient_t* C, long x) {
322 TRACE("coefficient::internal", "coefficient_assign_int()\n");
323 STAT_INCR(coefficient, assign_int)
324
325 if (C->type == COEFFICIENT_POLYNOMIAL) {
326 coefficient_destruct(C);
327 coefficient_construct_from_int(ctx, C, x);
328 } else {
329 integer_assign_int(ctx->K, &C->value.num, x);
330 }
331
332 assert(coefficient_is_normalized(ctx, C));
333 }
334
STAT_DECLARE(int,coefficient,assign_integer)335 STAT_DECLARE(int, coefficient, assign_integer)
336
337 void coefficient_assign_integer(const lp_polynomial_context_t* ctx, coefficient_t* C, const lp_integer_t* x) {
338 TRACE("coefficient::internal", "coefficient_assign_int()\n");
339 STAT_INCR(coefficient, assign_integer)
340
341 if (C->type == COEFFICIENT_POLYNOMIAL) {
342 coefficient_destruct(C);
343 coefficient_construct_from_integer(ctx, C, x);
344 } else {
345 integer_assign(ctx->K, &C->value.num, x);
346 }
347
348 assert(coefficient_is_normalized(ctx, C));
349 }
350
coefficient_lc_safe(const lp_polynomial_context_t * ctx,const coefficient_t * C,lp_variable_t x)351 const coefficient_t* coefficient_lc_safe(const lp_polynomial_context_t* ctx, const coefficient_t* C, lp_variable_t x) {
352 __var_unused(ctx);
353 switch (C->type) {
354 case COEFFICIENT_NUMERIC:
355 return C;
356 break;
357 case COEFFICIENT_POLYNOMIAL:
358 if (VAR(C) == x) {
359 return COEFF(C, SIZE(C) - 1);
360 } else {
361 assert(lp_variable_order_cmp(ctx->var_order, x, VAR(C)) > 0);
362 return C;
363 }
364 default:
365 assert(0);
366 return 0;
367 }
368 }
369
coefficient_lc(const coefficient_t * C)370 const coefficient_t* coefficient_lc(const coefficient_t* C) {
371 switch (C->type) {
372 case COEFFICIENT_NUMERIC:
373 return C;
374 break;
375 case COEFFICIENT_POLYNOMIAL:
376 return COEFF(C, SIZE(C) - 1);
377 break;
378 }
379 assert(0);
380 return 0;
381 }
382
coefficient_lc_m(const lp_polynomial_context_t * ctx,const coefficient_t * C,const lp_assignment_t * M)383 const coefficient_t* coefficient_lc_m(const lp_polynomial_context_t* ctx, const coefficient_t* C, const lp_assignment_t* M) {
384 switch (C->type) {
385 case COEFFICIENT_NUMERIC:
386 return C;
387 break;
388 case COEFFICIENT_POLYNOMIAL: {
389 // Locate the first non-zero coefficient past the top one
390 int i = SIZE(C) - 1;
391 while (i > 0 && coefficient_sgn(ctx, COEFF(C, i), M) == 0) {
392 -- i;
393 }
394 return COEFF(C, i);
395 break;
396 }
397 }
398 assert(0);
399 return 0;
400 }
401
402
coefficient_reductum(const lp_polynomial_context_t * ctx,coefficient_t * R,const coefficient_t * C)403 void coefficient_reductum(const lp_polynomial_context_t* ctx, coefficient_t* R, const coefficient_t* C) {
404
405 assert(C->type == COEFFICIENT_POLYNOMIAL);
406
407 // Locate the first non-zero ceofficient past the top one
408 int i = SIZE(C) - 2;
409 while (i >= 0 && coefficient_is_zero(ctx, COEFF(C, i))) {
410 -- i;
411 }
412
413 if (i < 0) {
414 // All zero
415 coefficient_assign_int(ctx, R, 0);
416 return;
417 }
418
419 coefficient_t result;
420 coefficient_construct_rec(ctx, &result, VAR(C), i + 1);
421
422 // Copy the other coefficients
423 while (i >= 0) {
424 if (!coefficient_is_zero(ctx, COEFF(C, i))) {
425 coefficient_assign(ctx, COEFF(&result, i), COEFF(C, i));
426 }
427 -- i;
428 }
429
430 coefficient_normalize(ctx, &result);
431 coefficient_swap(R, &result);
432 coefficient_destruct(&result);
433 }
434
coefficient_reductum_m(const lp_polynomial_context_t * ctx,coefficient_t * R,const coefficient_t * C,const lp_assignment_t * m,lp_polynomial_vector_t * assumptions)435 void coefficient_reductum_m(const lp_polynomial_context_t* ctx, coefficient_t* R, const coefficient_t* C, const lp_assignment_t* m, lp_polynomial_vector_t* assumptions) {
436
437 assert(C->type == COEFFICIENT_POLYNOMIAL);
438
439 // Locate the first non-zero ceofficient (normal reductum is the next nonzero)
440 int i = SIZE(C) - 1;
441 while (i >= 0 && coefficient_sgn(ctx, COEFF(C, i), m) == 0) {
442 if (assumptions != 0 && !coefficient_is_constant(COEFF(C, i))) {
443 lp_polynomial_vector_push_back_coeff(assumptions, COEFF(C, i));
444 }
445 -- i;
446 }
447
448 if (i < 0) {
449 // All zero
450 coefficient_assign_int(ctx, R, 0);
451 return;
452 } else if (assumptions != 0 && !coefficient_is_constant(COEFF(C, i))) {
453 lp_polynomial_vector_push_back_coeff(assumptions, COEFF(C, i));
454 }
455
456 coefficient_t result;
457 coefficient_construct_rec(ctx, &result, VAR(C), i + 1);
458
459 // Copy the other coefficients
460 while (i >= 0) {
461 if (!coefficient_is_zero(ctx, COEFF(C, i))) {
462 coefficient_assign(ctx, COEFF(&result, i), COEFF(C, i));
463 }
464 -- i;
465 }
466
467 coefficient_normalize(ctx, &result);
468 coefficient_swap(R, &result);
469 coefficient_destruct(&result);
470 }
471
coefficient_is_constant(const coefficient_t * C)472 int coefficient_is_constant(const coefficient_t* C) {
473 return C->type == COEFFICIENT_NUMERIC;
474 }
475
coefficient_degree(const coefficient_t * C)476 size_t coefficient_degree(const coefficient_t* C) {
477 switch (C->type) {
478 case COEFFICIENT_NUMERIC:
479 return 0;
480 break;
481 case COEFFICIENT_POLYNOMIAL:
482 return SIZE(C) - 1;
483 break;
484 }
485 assert(0);
486 return 0;
487 }
488
coefficient_degree_m(const lp_polynomial_context_t * ctx,const coefficient_t * C,const lp_assignment_t * M)489 size_t coefficient_degree_m(const lp_polynomial_context_t* ctx, const coefficient_t* C, const lp_assignment_t* M) {
490
491 if (trace_is_enabled("coefficient::roots")) {
492 tracef("coefficient_degree_m("); coefficient_print(ctx, C, trace_out); tracef(")\n");
493 }
494
495 switch (C->type) {
496 case COEFFICIENT_NUMERIC:
497 return 0;
498 break;
499 case COEFFICIENT_POLYNOMIAL: {
500 // Locate the first non-zero coefficient past the top one
501 size_t i = SIZE(C) - 1;
502 while (i > 0 && coefficient_sgn(ctx, COEFF(C, i), M) == 0) {
503 -- i;
504 }
505 return i;
506 break;
507 }
508 }
509 assert(0);
510 return 0;
511
512 }
513
coefficient_degree_safe(const lp_polynomial_context_t * ctx,const coefficient_t * C,lp_variable_t x)514 size_t coefficient_degree_safe(const lp_polynomial_context_t* ctx, const coefficient_t* C, lp_variable_t x) {
515 __var_unused(ctx);
516 switch (C->type) {
517 case COEFFICIENT_NUMERIC:
518 return 0;
519 break;
520 case COEFFICIENT_POLYNOMIAL:
521 if (VAR(C) == x) {
522 return SIZE(C) - 1;
523 } else {
524 assert(lp_variable_order_cmp(ctx->var_order, x, VAR(C)) > 0);
525 return 0;
526 }
527 break;
528 }
529 assert(0);
530 return 0;
531 }
532
coefficient_top_variable(const coefficient_t * C)533 lp_variable_t coefficient_top_variable(const coefficient_t* C) {
534 if (C->type == COEFFICIENT_POLYNOMIAL) {
535 return VAR(C);
536 } else {
537 return lp_variable_null;
538 }
539 }
540
coefficient_get_coefficient(const coefficient_t * C,size_t d)541 const coefficient_t* coefficient_get_coefficient(const coefficient_t* C, size_t d) {
542
543 assert(d <= coefficient_degree(C));
544
545 switch(C->type) {
546 case COEFFICIENT_NUMERIC:
547 return C;
548 break;
549 case COEFFICIENT_POLYNOMIAL:
550 return COEFF(C, d);
551 break;
552 }
553
554 assert(0);
555 return 0;
556 }
557
558 static coefficient_t zero;
559 static int zero_initialized = 0;
560
get_zero()561 static const coefficient_t* get_zero() {
562 if (!zero_initialized) {
563 zero_initialized = 1;
564 zero.type = COEFFICIENT_NUMERIC;
565 integer_construct(&zero.value.num);
566 }
567 return &zero;
568 }
569
coefficient_get_coefficient_safe(const lp_polynomial_context_t * ctx,const coefficient_t * C,size_t d,lp_variable_t x)570 const coefficient_t* coefficient_get_coefficient_safe(const lp_polynomial_context_t* ctx, const coefficient_t* C, size_t d, lp_variable_t x) {
571 __var_unused(ctx);
572
573 if (d > coefficient_degree_safe(ctx, C, x)) {
574 return get_zero();
575 }
576
577 switch(C->type) {
578 case COEFFICIENT_NUMERIC:
579 return C;
580 break;
581 case COEFFICIENT_POLYNOMIAL:
582 if (VAR(C) == x) {
583 return COEFF(C, d);
584 } else {
585 assert(d == 0);
586 return C;
587 }
588 break;
589 }
590
591 assert(0);
592 return 0;
593 }
594
595
STAT_DECLARE(int,coefficient,is_zero)596 STAT_DECLARE(int, coefficient, is_zero)
597
598 int coefficient_is_zero(const lp_polynomial_context_t* ctx, const coefficient_t* C) {
599 STAT_INCR(coefficient, is_zero)
600 return C->type == COEFFICIENT_NUMERIC && integer_is_zero(ctx->K, &C->value.num);
601 }
602
STAT_DECLARE(int,coefficient,is_one)603 STAT_DECLARE(int, coefficient, is_one)
604
605 int coefficient_is_one(const lp_polynomial_context_t* ctx, const coefficient_t* C) {
606 STAT_INCR(coefficient, is_one)
607 return C->type == COEFFICIENT_NUMERIC && integer_cmp_int(ctx->K, &C->value.num, 1) == 0;
608 }
609
STAT_DECLARE(int,coefficient,is_minus_one)610 STAT_DECLARE(int, coefficient, is_minus_one)
611
612 int coefficient_is_minus_one(const lp_polynomial_context_t* ctx, const coefficient_t* C) {
613 STAT_INCR(coefficient, is_minus_one)
614 return C->type == COEFFICIENT_NUMERIC && integer_cmp_int(ctx->K, &C->value.num, -1) == 0;
615 }
616
coefficient_value_approx(const lp_polynomial_context_t * ctx,const coefficient_t * C,const lp_assignment_t * m,lp_rational_interval_t * value)617 void coefficient_value_approx(const lp_polynomial_context_t* ctx, const coefficient_t* C, const lp_assignment_t* m, lp_rational_interval_t* value) {
618
619 if (trace_is_enabled("coefficient")) {
620 tracef("coefficient_value_approx("); coefficient_print(ctx, C, trace_out); tracef(")\n");
621 }
622
623 if (C->type == COEFFICIENT_NUMERIC) {
624 lp_rational_interval_t result;
625 lp_rational_interval_construct_from_integer(&result, &C->value.num, 0, &C->value.num, 0);
626 lp_rational_interval_swap(value, &result);
627 lp_rational_interval_destruct(&result);
628 } else {
629
630 lp_rational_interval_t result, tmp1, tmp2, x_value;
631
632 lp_rational_interval_construct_zero(&result);
633 lp_rational_interval_construct_zero(&tmp1);
634 lp_rational_interval_construct_zero(&tmp2);
635 lp_rational_interval_construct_zero(&x_value);
636
637 if (trace_is_enabled("coefficient")) {
638 tracef("coefficient_value_approx(): x = %s\n", lp_variable_db_get_name(ctx->var_db, VAR(C)));
639 tracef("assignment = "); lp_assignment_print(m, trace_out); tracef("\n");
640 }
641
642 lp_assignment_get_value_approx(m, VAR(C), &x_value);
643
644 // Get the value of x
645 if (trace_is_enabled("coefficient")) {
646 tracef("coefficient_value_approx(): x_value = ");
647 lp_rational_interval_print(&x_value, trace_out);
648 tracef("\n");
649 }
650
651 // We compute using powers, just an attempt to compute better. For example
652 // if p = x^2 + x and x = [-1, 1] then
653 // a) x(x + 1) = [-1, 1]*[0, 2] = [-2, 2]
654 // b) x^2 + x = [0, 1] + [-1, 1] = [-1, 2]
655 // we choose to do it b) way
656
657 // Compute
658 size_t i;
659 for (i = 0; i < SIZE(C); ++ i) {
660 if (!coefficient_is_zero(ctx, COEFF(C, i))) {
661 coefficient_value_approx(ctx, COEFF(C, i), m, &tmp1);
662 rational_interval_pow(&tmp2, &x_value, i);
663 // tracef("tmp2 = "); lp_rational_interval_print(&tmp2, trace_out); tracef("\n");
664 // tracef("tmp1 = "); lp_rational_interval_print(&tmp1, trace_out); tracef("\n");
665 rational_interval_mul(&tmp2, &tmp2, &tmp1);
666 // tracef("tmp2 = "); lp_rational_interval_print(&tmp2, trace_out); tracef("\n");
667 // tracef("result = "); lp_rational_interval_print(&result, trace_out); tracef("\n");
668 rational_interval_add(&result, &result, &tmp2);
669 // tracef("result = "); lp_rational_interval_print(&result, trace_out); tracef("\n");
670 }
671 }
672
673 lp_rational_interval_swap(&result, value);
674 lp_rational_interval_destruct(&x_value);
675 lp_rational_interval_destruct(&tmp1);
676 lp_rational_interval_destruct(&tmp2);
677 lp_rational_interval_destruct(&result);
678 }
679
680 if (trace_is_enabled("coefficient")) {
681 tracef("coefficient_value_approx() => "); lp_rational_interval_print(value, trace_out); tracef("\n");
682 }
683
684 }
685
686
687 /**
688 * C is an univariate polynomial C(x), we compute a bound L = 1/2^k such that
689 * any root of C(x) that is not zero is outside of [L, -L].
690 *
691 * If C(x) = c_n x^n + ... + c_1 x + c_0 then the lower bound L on roots is an
692 * upper bound 1/L on roots of C(1/x) = c_n + ... + c0 x^n.
693 *
694 * We compute the upper bound using the Cauchy bound
695 *
696 * bound = 1 + max(|c_0|, ... |c_n|)/|c_0|
697 *
698 * We take L = 1/bound, and so we compute k = log(max/c0) = log(max) - log(c0)
699 */
coefficient_root_lower_bound(const coefficient_t * C)700 unsigned coefficient_root_lower_bound(const coefficient_t* C) {
701
702 assert(C->type == COEFFICIENT_POLYNOMIAL);
703 assert(coefficient_is_univariate(C));
704
705 size_t i = 0;
706
707 // Find the first non-zero coefficient
708 while (integer_is_zero(lp_Z, &COEFF(C, i)->value.num)) {
709 ++ i;
710 assert(i < SIZE(C));
711 }
712
713 // First one (modulo the initial zeroes)
714 unsigned log_c0 = integer_log2_abs(&COEFF(C, i)->value.num);
715
716 // Get thge max log
717 unsigned max_log = log_c0;
718 for (++ i; i < SIZE(C); ++ i) {
719 assert(COEFF(C, i)->type == COEFFICIENT_NUMERIC);
720 if (!integer_is_zero(lp_Z, &COEFF(C, i)->value.num)) {
721 unsigned current_log = integer_log2_abs(&COEFF(C, i)->value.num);
722 if (current_log > max_log) {
723 max_log = current_log;
724 }
725 }
726 }
727
728 // Return the bound:
729 // * max_log is upper bound approximation, that's good
730 // * log_c0 is upper bound approximation, so we add one
731 // * max_log >= log_c0, so we're safe
732 return max_log - log_c0 + 1;
733 }
734
coefficient_is_assigned(const lp_polynomial_context_t * ctx,const coefficient_t * C,const lp_assignment_t * m)735 int coefficient_is_assigned(const lp_polynomial_context_t* ctx, const coefficient_t* C, const lp_assignment_t* m) {
736 if (C->type == COEFFICIENT_POLYNOMIAL) {
737 if (lp_assignment_get_value(m, VAR(C))->type == LP_VALUE_NONE) {
738 return 0;
739 } else {
740 size_t i;
741 for (i = 0; i < SIZE(C); ++ i) {
742 if (!coefficient_is_assigned(ctx, COEFF(C, i), m)) {
743 // Not assigned
744 return 0;
745 }
746 }
747 }
748 }
749 return 1;
750 }
751
STAT_DECLARE(int,coefficient,sgn)752 STAT_DECLARE(int, coefficient, sgn)
753
754 int coefficient_sgn(const lp_polynomial_context_t* ctx, const coefficient_t* C, const lp_assignment_t* m) {
755
756 if (trace_is_enabled("coefficient::sgn")) {
757 tracef("coefficient_sgn("); coefficient_print(ctx, C, trace_out); tracef(")\n");
758 }
759 STAT_INCR(coefficient, sgn)
760
761 assert(ctx->K == lp_Z);
762
763 int sgn;
764
765 if (C->type == COEFFICIENT_NUMERIC) {
766 // For numeric coefficients we're done
767 sgn = integer_sgn(lp_Z, &C->value.num);
768 if (trace_is_enabled("coefficient::sgn")) {
769 tracef("coefficient_sgn(): constant => %d\n", sgn);
770 }
771 } else {
772 assert(C->type == COEFFICIENT_POLYNOMIAL);
773
774 if (trace_is_enabled("coefficient::sgn")) {
775 tracef("coefficient_sgn(): evaluating in rationals\n");
776 }
777
778 // Try to evaluate in the rationals
779 coefficient_t C_rat;
780 coefficient_construct(ctx, &C_rat);
781 lp_integer_t multiplier;
782 integer_construct(&multiplier);
783 coefficient_evaluate_rationals(ctx, C, m, &C_rat, &multiplier);
784
785 if (trace_is_enabled("coefficient::sgn")) {
786 tracef("coefficient_sgn(): C_rat = "); coefficient_print(ctx, &C_rat, trace_out); tracef("\n");
787 }
788
789 // If constant, we're done
790 if (C_rat.type == COEFFICIENT_NUMERIC) {
791 // val(C) = C_rat/multiplier with multiplier positive
792 sgn = integer_sgn(lp_Z, &C_rat.value.num);
793 if (trace_is_enabled("coefficient::sgn")) {
794 tracef("coefficient_sgn(): constant => %d\n", sgn);
795 }
796 } else {
797
798 // Approximate the value of C_rat
799 lp_rational_interval_t C_rat_approx;
800 lp_rational_interval_construct_zero(&C_rat_approx);
801
802 if (trace_is_enabled("coefficient::sgn")) {
803 tracef("coefficient_sgn(): approximating with intervals\n");
804 }
805
806 // Approximate the value by doing interval computation
807 coefficient_value_approx(ctx, &C_rat, m, &C_rat_approx);
808
809 if (trace_is_enabled("coefficient::sgn")) {
810 tracef("coefficient_sgn(): approx => "); lp_rational_interval_print(&C_rat_approx, trace_out); tracef("\n");
811 }
812
813 if (C_rat_approx.is_point || !lp_rational_interval_contains_zero(&C_rat_approx)) {
814 // Safe to give the sign based on the interval bound
815 sgn = lp_rational_interval_sgn(&C_rat_approx);
816 if (trace_is_enabled("coefficient::sgn")) {
817 tracef("coefficient_sgn(): interval is good => %d\n", sgn);
818 }
819 } else {
820
821 //
822 // At this point the value is most likely 0.
823 //
824
825 //
826 // We're still not sure, we need to evaluate the sign using resultants.
827 // We construct the polynomial
828 //
829 // A(z, x1, ..., xn) = z - C_rat(x1, ..., xn).
830 //
831 // If the value of C_rat is a, and values of x_i is a_i with defining polynomial p_i(x_i)
832 // then we know that
833 //
834 // B(z) = Resultant(A, p1, ..., pn)
835 //
836 // has a as at least one zero and value of C_rat corresponds to one of
837 // those zeros.
838 //
839 // We estimate the bound L on the smallest *non-zero* root of B. We
840 // then estimate C_rat until either
841 // * it doesn't contain 0, when we can return the sign
842 // * it doesn't belong to (-L, L) => doesn't contain 0 => return the sign
843 // * it belongs to (-L, L), and therefore the sign must be 0
844
845 // The temporary variable, we'll be using
846 lp_variable_t z = lp_polynomial_context_get_temp_variable(ctx);
847
848 // A = z - C_rat
849 coefficient_t A;
850 lp_integer_t A_lc;
851 integer_construct_from_int(lp_Z, &A_lc, 1);
852 coefficient_construct_simple(ctx, &A, &A_lc, z, 1);
853 coefficient_sub(ctx, &A, &A, &C_rat);
854 lp_integer_destruct(&A_lc);
855
856 // List of variables in C_rat;
857 lp_variable_list_t C_rat_vars;
858 lp_variable_list_construct(&C_rat_vars);
859 coefficient_get_variables(&C_rat, &C_rat_vars);
860
861 // Cache the assignment intervals
862 lp_dyadic_interval_t* rational_interval_cache = algebraic_interval_remember(&C_rat_vars, m);
863
864 // Compute B (we keep it in A) by resolving out all the variables except z
865 lp_variable_order_make_bot(ctx->var_order, z);
866 coefficient_order(ctx, &A);
867 coefficient_resolve_algebraic(ctx, &A, m, &A);
868 lp_variable_order_make_bot(ctx->var_order, lp_variable_null);
869
870 if (trace_is_enabled("coefficient::sgn")) {
871 tracef("coefficient_sgn(): A = ");
872 coefficient_print(ctx, &A, trace_out);
873 tracef("\n");
874 }
875
876 // Get the lower bound on the roots
877 unsigned k = coefficient_root_lower_bound(&A);
878 // Interval (-L, L)
879 lp_rational_interval_t L_interval;
880 lp_rational_t L, L_neg;
881 rational_construct_from_int(&L, 1, 1);
882 rational_div_2exp(&L, &L, k);
883 rational_construct(&L_neg);
884 rational_neg(&L_neg, &L);
885 lp_rational_interval_construct(&L_interval, &L_neg, 1, &L, 1);
886
887 // Refine until done, i.e. C_rat_approx = (l, u) not contains 0, or
888 //
889 for (;;) {
890
891 if (trace_is_enabled("coefficient::sgn")) {
892 tracef("coefficient_sgn(): C_rat_approx = "); lp_rational_interval_print(&C_rat_approx, trace_out); tracef("\n");
893 tracef("coefficient_sgn(): L_interval = "); lp_rational_interval_print(&L_interval, trace_out); tracef("\n");
894 }
895
896 // If a point, we're done
897 if (lp_rational_interval_is_point(&C_rat_approx)) {
898 break;
899 }
900
901 // If no zero in interval, we're done
902 if (!lp_rational_interval_contains_zero(&C_rat_approx)) {
903 break;
904 }
905
906 // If contained in L, we're also done
907 int contains_a = lp_rational_interval_contains_rational(&L_interval, &C_rat_approx.a);
908 int contains_b = lp_rational_interval_contains_rational(&L_interval, &C_rat_approx.b);
909 if (contains_a && contains_b) {
910 assert(lp_rational_interval_contains_zero(&C_rat_approx));
911 break;
912 }
913
914 // Refine the values
915 size_t i;
916 for (i = 0; i < C_rat_vars.list_size; ++ i) {
917 lp_variable_t x_i = C_rat_vars.list[i];
918 const lp_value_t* x_i_value = lp_assignment_get_value(m, x_i);
919 if (!lp_value_is_rational(x_i_value)) {
920 lp_algebraic_number_refine_const(&x_i_value->value.a);
921 }
922 }
923
924 // Approximate the value by doing interval computation
925 coefficient_value_approx(ctx, &C_rat, m, &C_rat_approx);
926 }
927
928 // Restore the cache (will free the cache)
929 algebraic_interval_restore(&C_rat_vars, rational_interval_cache, m);
930
931 // Release the temporary variable
932 lp_polynomial_context_release_temp_variable(ctx, z);
933
934 // Destruct temps
935 coefficient_destruct(&A);
936 lp_variable_list_destruct(&C_rat_vars);
937 lp_rational_interval_destruct(&L_interval);
938 lp_rational_destruct(&L);
939 lp_rational_destruct(&L_neg);
940
941 // Safe to give the sign based on the interval bound
942 // * interval_sgn returns 0 if 0 in interval, otherwise the sign
943 sgn = lp_rational_interval_sgn(&C_rat_approx);
944 if (trace_is_enabled("coefficient::sgn")) {
945 tracef("coefficient_sgn(): interval is good => %d\n", sgn);
946 }
947 }
948
949 // Destruct temps
950 lp_rational_interval_destruct(&C_rat_approx);
951 }
952
953 // Destruct temps
954 lp_integer_destruct(&multiplier);
955 coefficient_destruct(&C_rat);
956 }
957
958 if (trace_is_enabled("coefficient::sgn")) {
959 tracef("coefficient_sgn() => %d\n", sgn);
960 }
961
962 if (sgn < 0) { return -1; }
963 if (sgn > 0) { return 1; }
964
965 return 0;
966 }
967
STAT_DECLARE(int,coefficient,interval_value)968 STAT_DECLARE(int, coefficient, interval_value)
969
970 void coefficient_interval_value(const lp_polynomial_context_t* ctx, const coefficient_t* C, const lp_interval_assignment_t* m, lp_interval_t* out) {
971
972 if (trace_is_enabled("coefficient::interval")) {
973 tracef("coefficient_interval_value("); coefficient_print(ctx, C, trace_out); tracef(", "); lp_interval_assignment_print(m, trace_out); tracef(")\n");
974 }
975
976 if (C->type == COEFFICIENT_NUMERIC) {
977 lp_value_t value;
978 lp_value_construct(&value, LP_VALUE_INTEGER, &C->value.num);
979 lp_interval_t value_interval;
980 lp_interval_construct_point(&value_interval, &value);
981 lp_interval_swap(out, &value_interval);
982 lp_interval_destruct(&value_interval);
983 lp_value_destruct(&value);
984 } else {
985
986 lp_interval_t result, tmp1, tmp2;
987
988 lp_interval_construct_zero(&result);
989 lp_interval_construct_zero(&tmp1);
990 lp_interval_construct_zero(&tmp2);
991
992 if (trace_is_enabled("coefficient::interval")) {
993 tracef("coefficient_interval_value(): x = %s\n", lp_variable_db_get_name(ctx->var_db, VAR(C)));
994 tracef("assignment = "); lp_interval_assignment_print(m, trace_out); tracef("\n");
995 }
996
997 const lp_interval_t* x_value = lp_interval_assignment_get_interval(m, VAR(C));
998 assert(x_value);
999
1000 // Get the value of x
1001 if (trace_is_enabled("coefficient::interval")) {
1002 tracef("coefficient_interval_value(): x_value = ");
1003 lp_interval_print(x_value, trace_out);
1004 tracef("\n");
1005 }
1006
1007 // We compute using powers, just an attempt to compute better. For example
1008 // if p = x^2 + x and x = [-1, 1] then
1009 // a) x(x + 1) = [-1, 1]*[0, 2] = [-2, 2]
1010 // b) x^2 + x = [0, 1] + [-1, 1] = [-1, 2]
1011 // we choose to do it b) way
1012
1013 // Compute
1014 size_t i;
1015 for (i = 0; i < SIZE(C); ++ i) {
1016 if (!coefficient_is_zero(ctx, COEFF(C, i))) {
1017 // tracef("i = %zu\n", i);
1018 // tracef("x = "); lp_interval_print(x_value, trace_out); tracef("\n");
1019 /*
1020 * BD: this may have a side-effect on m (via lp_assignment_ensure_size)
1021 * which make x_value an invalid pointer.
1022 */
1023 coefficient_interval_value(ctx, COEFF(C, i), m, &tmp1);
1024 lp_interval_pow(&tmp2, x_value, i);
1025 // tracef("tmp2 = x^i = "); lp_interval_print(&tmp2, trace_out); tracef("\n");
1026 // tracef("tmp1 = "); lp_interval_print(&tmp1, trace_out); tracef("\n");
1027 lp_interval_mul(&tmp2, &tmp2, &tmp1);
1028 // tracef("tmp2 = tmp2 * tmp1 = "); lp_interval_print(&tmp2, trace_out); tracef("\n");
1029 // tracef("result = "); lp_interval_print(&result, trace_out); tracef("\n");
1030 lp_interval_add(&result, &result, &tmp2);
1031 // tracef("result = result + tmp2 = "); lp_interval_print(&result, trace_out); tracef("\n");
1032 }
1033 }
1034
1035 lp_interval_swap(&result, out);
1036 lp_interval_destruct(&tmp1);
1037 lp_interval_destruct(&tmp2);
1038 lp_interval_destruct(&result);
1039 }
1040
1041 if (trace_is_enabled("coefficient::interval")) {
1042 tracef("coefficient_value_approx() => "); lp_interval_print(out, trace_out); tracef("\n");
1043 }
1044 }
1045
STAT_DECLARE(int,coefficient,lc_sgn)1046 STAT_DECLARE(int, coefficient, lc_sgn)
1047
1048 int coefficient_lc_sgn(const lp_polynomial_context_t* ctx, const coefficient_t* C) {
1049 STAT_INCR(coefficient, lc_sgn)
1050
1051 while (C->type != COEFFICIENT_NUMERIC) {
1052 C = coefficient_lc(C);
1053 }
1054
1055 return integer_sgn(ctx->K, &C->value.num);
1056 }
1057
1058
STAT_DECLARE(int,coefficient,in_order)1059 STAT_DECLARE(int, coefficient, in_order)
1060
1061 int coefficient_in_order(const lp_polynomial_context_t* ctx, const coefficient_t* C) {
1062 TRACE("coefficient::internal", "coefficient_in_order()\n");
1063 STAT_INCR(coefficient, in_order)
1064
1065 size_t i;
1066 switch (C->type) {
1067 case COEFFICIENT_NUMERIC:
1068 return 1;
1069 break;
1070 case COEFFICIENT_POLYNOMIAL:
1071 // Check that the top is bigger than the top of coefficient and run
1072 // recursively
1073 for (i = 0; i < SIZE(C); ++ i) {
1074 const coefficient_t* C_i = COEFF(C, i);
1075 if (C_i->type == COEFFICIENT_POLYNOMIAL) {
1076 if (lp_variable_order_cmp(ctx->var_order, VAR(C), VAR(C_i)) <= 0) {
1077 // Top variable must be bigger than others
1078 return 0;
1079 } else if (!coefficient_in_order(ctx, C_i)) {
1080 return 0;
1081 }
1082 }
1083 }
1084
1085 break;
1086 }
1087
1088 return 1;
1089 }
1090
1091
coefficient_cmp_general(const lp_polynomial_context_t * ctx,const coefficient_t * C1,const coefficient_t * C2,int compare_values)1092 int coefficient_cmp_general(const lp_polynomial_context_t* ctx, const coefficient_t* C1, const coefficient_t* C2, int compare_values) {
1093 int cmp;
1094
1095 if (C1->type == COEFFICIENT_NUMERIC && C2->type == COEFFICIENT_NUMERIC) {
1096 if (compare_values) {
1097 cmp = integer_cmp(ctx->K, &C1->value.num, &C2->value.num);
1098 } else {
1099 cmp = 0;
1100 }
1101 } else if (C1->type == COEFFICIENT_NUMERIC) {
1102 // C1 is a constant, C1 is always smaller
1103 return -1;
1104 } else if (C2->type == COEFFICIENT_NUMERIC) {
1105 // C2 is a constant, C1 is always bigger
1106 return 1;
1107 } else {
1108 // Both are polynomials, compare the variable
1109 int var_cmp = lp_variable_order_cmp(ctx->var_order, VAR(C1), VAR(C2));
1110 if (var_cmp == 0) {
1111 if (compare_values) {
1112 // If the variables are the same, compare lexicographically
1113 int deg_cmp = ((int) SIZE(C1)) - ((int) SIZE(C2));
1114 if (deg_cmp == 0) {
1115 int i = C1->value.rec.size - 1;
1116 for (; i >= 0; -- i) {
1117 int coeff_cmp = coefficient_cmp_general(ctx, COEFF(C1, i), COEFF(C2, i), compare_values);
1118 if (coeff_cmp != 0) {
1119 cmp = coeff_cmp;
1120 break;
1121 }
1122 }
1123 if (i < 0) {
1124 // All coefficients equal
1125 cmp = 0;
1126 }
1127 } else {
1128 cmp = deg_cmp;
1129 }
1130 } else {
1131 return 0;
1132 }
1133 } else {
1134 // Variable comparison is enough
1135 cmp = var_cmp;
1136 }
1137 }
1138
1139 TRACE("coefficien::internal", "coefficient_cmp() => %d\n", cmp);
1140 return cmp;
1141 }
1142
STAT_DECLARE(int,coefficient,cmp)1143 STAT_DECLARE(int, coefficient, cmp)
1144
1145 int coefficient_cmp(const lp_polynomial_context_t* ctx, const coefficient_t* C1, const coefficient_t* C2) {
1146 TRACE("coefficient", "coefficient_cmp()\n");
1147 STAT_INCR(coefficient, cmp)
1148
1149 return coefficient_cmp_general(ctx, C1, C2, 1);
1150 }
1151
STAT_DECLARE(int,coefficient,cmp_type)1152 STAT_DECLARE(int, coefficient, cmp_type)
1153
1154 int coefficient_cmp_type(const lp_polynomial_context_t* ctx, const coefficient_t* C1, const coefficient_t* C2) {
1155 TRACE("coefficient::internal", "coefficient_cmp_type()\n");
1156 STAT_INCR(coefficient, cmp_type)
1157
1158 return coefficient_cmp_general(ctx, C1, C2, 0);
1159 }
1160
coefficient_traverse(const lp_polynomial_context_t * ctx,const coefficient_t * C,traverse_f f,lp_monomial_t * m,void * data)1161 void coefficient_traverse(const lp_polynomial_context_t* ctx, const coefficient_t* C, traverse_f f, lp_monomial_t* m, void* data) {
1162
1163 if (trace_is_enabled("coefficient::order")) {
1164 tracef("order = "); lp_variable_order_print(ctx->var_order, ctx->var_db, trace_out); tracef("\n");
1165 tracef("C = "); coefficient_print(ctx, C, trace_out); tracef("\n");
1166 tracef("m = "); monomial_print(ctx, m, trace_out); tracef("\n");
1167 }
1168
1169 size_t d;
1170 switch (C->type) {
1171 case COEFFICIENT_NUMERIC:
1172 integer_assign(ctx->K, &m->a, &C->value.num);
1173 (*f)(ctx, m, data);
1174 break;
1175 case COEFFICIENT_POLYNOMIAL:
1176 // The constant
1177 if (!coefficient_is_zero(ctx, COEFF(C, 0))) {
1178 coefficient_traverse(ctx, COEFF(C, 0), f, m, data);
1179 }
1180 // Power of x
1181 for (d = 1; d < SIZE(C); ++ d) {
1182 if (!coefficient_is_zero(ctx, COEFF(C, d))) {
1183 lp_monomial_push(m, VAR(C), d);
1184 coefficient_traverse(ctx, COEFF(C, d), f, m, data);
1185 lp_monomial_pop(m);
1186 }
1187 }
1188 break;
1189 }
1190 }
1191
1192 /**
1193 * Method called to add a monomial to C. The monomial should be ordered in the
1194 * same order as C, top variable at the m[0].
1195 */
coefficient_add_ordered_monomial(const lp_polynomial_context_t * ctx,lp_monomial_t * m,void * C_void)1196 void coefficient_add_ordered_monomial(const lp_polynomial_context_t* ctx, lp_monomial_t* m, void* C_void) {
1197
1198 coefficient_t* C = (coefficient_t*) C_void;
1199
1200 if (trace_is_enabled("coefficient::order")) {
1201 tracef("coefficient_add_monomial():\n");
1202 tracef("m = "); monomial_print(ctx, m, trace_out); tracef("\n");
1203 tracef("C = "); coefficient_print(ctx, C, trace_out); tracef("\n");
1204 }
1205
1206 if (m->n == 0) {
1207 // Just add a constant to C
1208 switch (C->type) {
1209 case COEFFICIENT_NUMERIC:
1210 integer_add(ctx->K, &C->value.num, &C->value.num, &m->a);
1211 break;
1212 case COEFFICIENT_POLYNOMIAL:
1213 coefficient_add_ordered_monomial(ctx, m, COEFF(C, 0));
1214 break;
1215 }
1216 } else {
1217 // Proper monomial (first variable is the top one)
1218 lp_variable_t x = m->p[0].x;
1219 unsigned d = m->p[0].d;
1220 // Compare the variables
1221 if (C->type == COEFFICIENT_NUMERIC || lp_variable_order_cmp(ctx->var_order, x, VAR(C)) >= 0) {
1222 coefficient_ensure_capacity(ctx, C, x, d+1);
1223 // Now, add the monomial to the right place
1224 m->p ++;
1225 m->n --;
1226 coefficient_add_ordered_monomial(ctx, m, COEFF(C, d));
1227 coefficient_normalize(ctx, C);
1228 m->p --;
1229 m->n ++;
1230 } else {
1231 coefficient_add_ordered_monomial(ctx, m, COEFF(C, 0));
1232 }
1233 }
1234
1235 assert(coefficient_is_normalized(ctx, C));
1236 }
1237
coefficient_order_and_add_monomial(const lp_polynomial_context_t * ctx,lp_monomial_t * m,void * C_void)1238 void coefficient_order_and_add_monomial(const lp_polynomial_context_t* ctx, lp_monomial_t* m, void* C_void) {
1239 lp_monomial_t m_ordered;
1240 lp_monomial_construct_copy(ctx, &m_ordered, m, /** sort */ 1);
1241 coefficient_add_ordered_monomial(ctx, &m_ordered, C_void);
1242 lp_monomial_destruct(&m_ordered);
1243 }
1244
coefficient_add_monomial(const lp_polynomial_context_t * ctx,coefficient_t * C,const lp_monomial_t * m)1245 void coefficient_add_monomial(const lp_polynomial_context_t* ctx, coefficient_t* C, const lp_monomial_t* m) {
1246 lp_monomial_t m_ordered;
1247 lp_monomial_construct_copy(ctx, &m_ordered, m, /** sort */ 1);
1248 coefficient_add_ordered_monomial(ctx, &m_ordered, C);
1249 lp_monomial_destruct(&m_ordered);
1250 }
1251
STAT_DECLARE(int,coefficient,order)1252 STAT_DECLARE(int, coefficient, order)
1253
1254 void coefficient_order(const lp_polynomial_context_t* ctx, coefficient_t* C) {
1255 TRACE("coefficient", "coefficient_order()\n");
1256 STAT_INCR(coefficient, order)
1257
1258 if (C->type == COEFFICIENT_NUMERIC) {
1259 // Numeric coefficients are always OK
1260 return;
1261 }
1262
1263 if (trace_is_enabled("coefficient::order")) {
1264 tracef("order = "); lp_variable_order_print(ctx->var_order, ctx->var_db, trace_out); tracef("\n");
1265 tracef("C = "); coefficient_print(ctx, C, trace_out); tracef("\n");
1266 }
1267
1268 // The coefficient we are building
1269 coefficient_t result;
1270 coefficient_construct(ctx, &result);
1271 // The monomials build in the original order
1272 lp_monomial_t m_tmp;
1273 lp_monomial_construct(ctx, &m_tmp);
1274 // For each monomial of C, add it to the result
1275 coefficient_traverse(ctx, C, coefficient_order_and_add_monomial, &m_tmp, &result);
1276 // Keep the result
1277 coefficient_swap(C, &result);
1278 // Destroy temps
1279 lp_monomial_destruct(&m_tmp);
1280 coefficient_destruct(&result);
1281
1282 assert(coefficient_is_normalized(ctx, C));
1283 }
1284
STAT_DECLARE(int,coefficient,add)1285 STAT_DECLARE(int, coefficient, add)
1286
1287 #define MAX(x, y) (x >= y ? x : y)
1288
1289 void coefficient_add(const lp_polynomial_context_t* ctx, coefficient_t* S, const coefficient_t* C1, const coefficient_t* C2) {
1290 TRACE("coefficient::arith", "coefficient_add()\n");
1291 STAT_INCR(coefficient, add)
1292
1293 if (trace_is_enabled("coefficient::arith")) {
1294 tracef("S = "); coefficient_print(ctx, S, trace_out); tracef("\n");
1295 tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
1296 tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
1297 }
1298
1299 coefficient_t result;
1300
1301 int type_cmp = coefficient_cmp_type(ctx, C1, C2);
1302
1303 if (type_cmp == 0) {
1304 if (C1->type == COEFFICIENT_NUMERIC) {
1305 assert(C2->type == COEFFICIENT_NUMERIC);
1306 // Add the integers
1307 integer_add(ctx->K, &S->value.num, &C1->value.num, &C2->value.num);
1308 } else {
1309 assert(C1->type == COEFFICIENT_POLYNOMIAL);
1310 assert(C2->type == COEFFICIENT_POLYNOMIAL);
1311 assert(VAR(C1) == VAR(C2));
1312 // Two polynomials over the same top variable
1313 size_t max_size = MAX(SIZE(C1), SIZE(C2));
1314 coefficient_construct_rec(ctx, &result, VAR(C1), max_size);
1315 size_t i;
1316 for (i = 0; i < max_size; ++ i) {
1317 if (i < SIZE(C1)) {
1318 if (i < SIZE(C2)) {
1319 // add C1 and C2
1320 coefficient_add(ctx, COEFF(&result, i), COEFF(C1, i), COEFF(C2, i));
1321 } else {
1322 // copy C1 coefficient
1323 coefficient_assign(ctx, COEFF(&result, i), COEFF(C1, i));
1324 }
1325 } else {
1326 // copy C2 coefficient
1327 coefficient_assign(ctx, COEFF(&result, i), COEFF(C2, i));
1328 }
1329 }
1330 coefficient_normalize(ctx, &result);
1331 coefficient_swap(&result, S);
1332 coefficient_destruct(&result);
1333 }
1334 } else if (type_cmp > 0) {
1335 // C1 > C2, add C2 into the constant of C1
1336 // We can't assign S to C1, since C2 might be S, so we use a temp
1337 coefficient_construct_copy(ctx, &result, C1);
1338 coefficient_add(ctx, COEFF(&result, 0), COEFF(C1, 0), C2);
1339 coefficient_swap(&result, S);
1340 coefficient_destruct(&result);
1341 // Since C1 is not a constant, no normalization needed, same size
1342 } else {
1343 // C1 < C2, add C1 into the constant of C2
1344 // We can't assign C2 to S1, since C1 might be S, so we use a temp
1345 coefficient_construct_copy(ctx, &result, C2);
1346 coefficient_add(ctx, COEFF(&result, 0), C1, COEFF(C2, 0));
1347 coefficient_swap(&result, S);
1348 coefficient_destruct(&result);
1349 // Since C2 is not a constant, no normalization needed, same size
1350 }
1351
1352 if (trace_is_enabled("coefficient::arith")) {
1353 tracef("add = "); coefficient_print(ctx, S, trace_out); tracef("\n");
1354 }
1355
1356 assert(coefficient_is_normalized(ctx, S));
1357 }
1358
STAT_DECLARE(int,coefficient,neg)1359 STAT_DECLARE(int, coefficient, neg)
1360
1361 void coefficient_neg(const lp_polynomial_context_t* ctx, coefficient_t* N, const coefficient_t* C) {
1362 TRACE("coefficient::arith", "coefficient_neg()\n");
1363 STAT_INCR(coefficient, neg)
1364
1365 size_t i;
1366 coefficient_t result;
1367
1368 switch (C->type) {
1369 case COEFFICIENT_NUMERIC:
1370 if (N->type == COEFFICIENT_POLYNOMIAL) {
1371 coefficient_destruct(N);
1372 coefficient_construct(ctx, N);
1373 }
1374 integer_neg(ctx->K, &N->value.num, &C->value.num);
1375 break;
1376 case COEFFICIENT_POLYNOMIAL:
1377 if (N != C) {
1378 coefficient_construct_rec(ctx, &result, VAR(C), SIZE(C));
1379 for (i = 0; i < SIZE(C); ++i) {
1380 if (!coefficient_is_zero(ctx, COEFF(C, i))) {
1381 coefficient_neg(ctx, COEFF(&result, i), COEFF(C, i));
1382 }
1383 }
1384 coefficient_normalize(ctx, &result);
1385 coefficient_swap(&result, N);
1386 coefficient_destruct(&result);
1387 } else {
1388 // In-place negation
1389 for (i = 0; i < SIZE(C); ++i) {
1390 if (!coefficient_is_zero(ctx, COEFF(C, i))) {
1391 coefficient_neg(ctx, COEFF(N, i), COEFF(C, i));
1392 }
1393 }
1394 }
1395 break;
1396 }
1397
1398 assert(coefficient_is_normalized(ctx, N));
1399 }
1400
STAT_DECLARE(int,coefficient,sub)1401 STAT_DECLARE(int, coefficient, sub)
1402
1403 void coefficient_sub(const lp_polynomial_context_t* ctx, coefficient_t* S, const coefficient_t* C1, const coefficient_t* C2) {
1404 TRACE("coefficient::arith", "coefficient_sub()\n");
1405 STAT_INCR(coefficient, sub)
1406
1407 if (trace_is_enabled("coefficient::arith")) {
1408 tracef("S = "); coefficient_print(ctx, S, trace_out); tracef("\n");
1409 tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
1410 tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
1411 }
1412
1413 coefficient_t result;
1414
1415 int type_cmp = coefficient_cmp_type(ctx, C1, C2);
1416
1417 if (type_cmp == 0) {
1418 if (C1->type == COEFFICIENT_NUMERIC) {
1419 assert(C2->type == COEFFICIENT_NUMERIC);
1420 // Subtract the integers
1421 integer_sub(ctx->K, &S->value.num, &C1->value.num, &C2->value.num);
1422 } else {
1423 assert(C1->type == COEFFICIENT_POLYNOMIAL);
1424 assert(C2->type == COEFFICIENT_POLYNOMIAL);
1425 // Two polynomials over the same top variable
1426 assert(VAR(C1) == VAR(C2));
1427 size_t max_size = MAX(SIZE(C1), SIZE(C2));
1428 coefficient_construct_rec(ctx, &result, VAR(C1), max_size);
1429 size_t i;
1430 for (i = 0; i < max_size; ++ i) {
1431 if (i < SIZE(C1)) {
1432 if (i < SIZE(C2)) {
1433 // subtract C1 and C2
1434 coefficient_sub(ctx, COEFF(&result, i), COEFF(C1, i), COEFF(C2, i));
1435 } else {
1436 // copy C1 coefficient
1437 coefficient_assign(ctx, COEFF(&result, i), COEFF(C1, i));
1438 }
1439 } else {
1440 // copy -C2 ceofficient
1441 coefficient_neg(ctx, COEFF(&result, i), COEFF(C2, i));
1442 }
1443 }
1444 coefficient_normalize(ctx, &result);
1445 coefficient_swap(&result, S);
1446 coefficient_destruct(&result);
1447 }
1448 } else if (type_cmp > 0) {
1449 // C1 > C2, subtract C2 into the constant of C1
1450 // Can't assign C1 to S, since C2 might be S
1451 coefficient_construct_copy(ctx, &result, C1);
1452 coefficient_sub(ctx, COEFF(&result, 0), COEFF(C1, 0), C2);
1453 coefficient_swap(&result, S);
1454 coefficient_destruct(&result);
1455 // Since C1 is not a constant, no normalization is needed
1456 } else {
1457 // C1 < C2, compute -(C2 - C1)
1458 coefficient_sub(ctx, S, C2, C1);
1459 coefficient_neg(ctx, S, S);
1460 // Since C2 is not a constant, no normalization is needed
1461 }
1462
1463 assert(coefficient_is_normalized(ctx, S));
1464 }
1465
1466
1467 void coefficient_add_mul(const lp_polynomial_context_t* ctx, coefficient_t* S, const coefficient_t* C1, const coefficient_t* C2);
1468
STAT_DECLARE(int,coefficient,mul)1469 STAT_DECLARE(int, coefficient, mul)
1470
1471 void coefficient_mul(const lp_polynomial_context_t* ctx, coefficient_t* P, const coefficient_t* C1, const coefficient_t* C2) {
1472 TRACE("coefficient::arith", "coefficient_mul()\n");
1473 STAT_INCR(coefficient, mul)
1474
1475 if (trace_is_enabled("coefficient::arith")) {
1476 tracef("P = "); coefficient_print(ctx, P, trace_out); tracef("\n");
1477 tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
1478 tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
1479 }
1480
1481 size_t i, j;
1482 coefficient_t result;
1483
1484 int type_cmp = coefficient_cmp_type(ctx, C1, C2);
1485
1486 if (type_cmp == 0) {
1487 if (C1->type == COEFFICIENT_NUMERIC) {
1488 assert(C2->type == COEFFICIENT_NUMERIC);
1489 // Multiply the integers
1490 integer_mul(ctx->K, &P->value.num, &C1->value.num, &C2->value.num);
1491 } else {
1492 assert(C1->type == COEFFICIENT_POLYNOMIAL);
1493 assert(C2->type == COEFFICIENT_POLYNOMIAL);
1494 // Two polynomials over the same top variable
1495 assert(VAR(C1) == VAR(C2));
1496 coefficient_construct_rec(ctx, &result, VAR(C1), SIZE(C1) + SIZE(C2) - 1);
1497 for (i = 0; i < SIZE(C1); ++ i) {
1498 if (!coefficient_is_zero(ctx, COEFF(C1, i))) {
1499 for (j = 0; j < SIZE(C2); ++ j) {
1500 if (!coefficient_is_zero(ctx, COEFF(C2, j))) {
1501 coefficient_add_mul(ctx, COEFF(&result, i + j), COEFF(C1, i), COEFF(C2, j));
1502 if (trace_is_enabled("coefficient::arith")) {
1503 tracef("result = "); coefficient_print(ctx, &result, trace_out); tracef("\n");
1504 }
1505 }
1506 }
1507 }
1508 }
1509 coefficient_normalize(ctx, &result);
1510 coefficient_swap(&result, P);
1511 coefficient_destruct(&result);
1512 }
1513 } else if (type_cmp > 0) {
1514 assert(C1->type == COEFFICIENT_POLYNOMIAL);
1515 // C1 > C2, multiply each coefficient of C1 with C2
1516 coefficient_construct_rec(ctx, &result, VAR(C1), SIZE(C1));
1517 for (i = 0; i < SIZE(C1); ++ i) {
1518 coefficient_mul(ctx, COEFF(&result, i), COEFF(C1, i), C2);
1519 }
1520 coefficient_normalize(ctx, &result);
1521 coefficient_swap(&result, P);
1522 coefficient_destruct(&result);
1523 } else {
1524 // C1 < C2, multiply each coefficient of C2 with C1
1525 coefficient_construct_rec(ctx, &result, VAR(C2), SIZE(C2));
1526 for (i = 0; i < SIZE(C2); ++ i) {
1527 if (!coefficient_is_zero(ctx, COEFF(C2, i))) {
1528 coefficient_mul(ctx, COEFF(&result, i), C1, COEFF(C2, i));
1529 }
1530 }
1531 coefficient_normalize(ctx, &result);
1532 coefficient_swap(&result, P);
1533 coefficient_destruct(&result);
1534 }
1535
1536 if (trace_is_enabled("coefficient::arith")) {
1537 tracef("mul = "); coefficient_print(ctx, P, trace_out); tracef("\n");
1538 }
1539
1540 assert(coefficient_is_normalized(ctx, P));
1541 }
1542
STAT_DECLARE(int,coefficient,mul_int)1543 STAT_DECLARE(int, coefficient, mul_int)
1544
1545 void coefficient_mul_int(const lp_polynomial_context_t* ctx, coefficient_t* P, const coefficient_t* C, long a) {
1546 TRACE("coefficient::arith", "coefficient_mul_int()\n");
1547 STAT_INCR(coefficient, mul_int)
1548
1549 if (trace_is_enabled("coefficient::arith")) {
1550 tracef("P = "); coefficient_print(ctx, P, trace_out); tracef("\n");
1551 tracef("C = "); coefficient_print(ctx, C, trace_out); tracef("\n");
1552 tracef("n = %ld\n", a);
1553 }
1554
1555 size_t i;
1556 coefficient_t result;
1557
1558 if (C->type == COEFFICIENT_NUMERIC) {
1559 if (P->type == COEFFICIENT_POLYNOMIAL) {
1560 coefficient_construct(ctx, &result);
1561 integer_mul_int(ctx->K, &result.value.num, &C->value.num, a);
1562 coefficient_swap(&result, P);
1563 coefficient_destruct(&result);
1564 } else {
1565 integer_mul_int(ctx->K, &P->value.num, &C->value.num, a);
1566 }
1567 } else {
1568 coefficient_construct_rec(ctx, &result, VAR(C), SIZE(C));
1569 for (i = 0; i < SIZE(C); ++ i) {
1570 if (!coefficient_is_zero(ctx, COEFF(C, i))) {
1571 coefficient_mul_int(ctx, COEFF(&result, i), COEFF(C, i), a);
1572 }
1573 }
1574 coefficient_normalize(ctx, &result);
1575 coefficient_swap(&result, P);
1576 coefficient_destruct(&result);
1577 }
1578
1579 if (trace_is_enabled("coefficient::arith")) {
1580 tracef("mul = "); coefficient_print(ctx, P, trace_out); tracef("\n");
1581 }
1582
1583 assert(coefficient_is_normalized(ctx, P));
1584 }
1585
STAT_DECLARE(int,coefficient,mul_integer)1586 STAT_DECLARE(int, coefficient, mul_integer)
1587
1588 void coefficient_mul_integer(const lp_polynomial_context_t* ctx, coefficient_t* P, const coefficient_t* C, const lp_integer_t* a) {
1589 TRACE("coefficient::arith", "coefficient_mul_int()\n");
1590 STAT_INCR(coefficient, mul_int)
1591
1592 if (trace_is_enabled("coefficient::arith")) {
1593 tracef("P = "); coefficient_print(ctx, P, trace_out); tracef("\n");
1594 tracef("C = "); coefficient_print(ctx, C, trace_out); tracef("\n");
1595 tracef("n = "); integer_print(a, trace_out);
1596 }
1597
1598 size_t i;
1599 coefficient_t result;
1600
1601 if (C->type == COEFFICIENT_NUMERIC) {
1602 if (P->type == COEFFICIENT_POLYNOMIAL) {
1603 coefficient_construct(ctx, &result);
1604 integer_mul(ctx->K, &result.value.num, &C->value.num, a);
1605 coefficient_swap(&result, P);
1606 coefficient_destruct(&result);
1607 } else {
1608 integer_mul(ctx->K, &P->value.num, &C->value.num, a);
1609 }
1610 } else {
1611 coefficient_construct_rec(ctx, &result, VAR(C), SIZE(C));
1612 for (i = 0; i < SIZE(C); ++ i) {
1613 if (!coefficient_is_zero(ctx, COEFF(C, i))) {
1614 coefficient_mul_integer(ctx, COEFF(&result, i), COEFF(C, i), a);
1615 }
1616 }
1617 coefficient_normalize(ctx, &result);
1618 coefficient_swap(&result, P);
1619 coefficient_destruct(&result);
1620 }
1621
1622 assert(coefficient_is_normalized(ctx, P));
1623 }
1624
1625
STAT_DECLARE(int,coefficient,shl)1626 STAT_DECLARE(int, coefficient, shl)
1627
1628 void coefficient_shl(const lp_polynomial_context_t* ctx, coefficient_t* S, const coefficient_t* C, lp_variable_t x, unsigned n) {
1629 TRACE("coefficient::arith", "coefficient_shl()\n");
1630 STAT_INCR(coefficient, shl)
1631
1632 if (trace_is_enabled("coefficient::arith")) {
1633 tracef("C = "); coefficient_print(ctx, C, trace_out); tracef("\n");
1634 tracef("x = %s\n", lp_variable_db_get_name(ctx->var_db, x));
1635 tracef("n = %u\n", n);
1636 }
1637
1638 coefficient_assign(ctx, S, C);
1639 if (!coefficient_is_zero(ctx, C) && n > 0) {
1640 int old_size = (S->type == COEFFICIENT_NUMERIC || VAR(S) != x) ? 1 : SIZE(S);
1641 coefficient_ensure_capacity(ctx, S, x, old_size + n);
1642 int i;
1643 for (i = old_size - 1; i >= 0; -- i) {
1644 if (!coefficient_is_zero(ctx, COEFF(S, i))) {
1645 coefficient_swap(COEFF(S, i + n), COEFF(S, i));
1646 }
1647 }
1648 }
1649
1650 if (trace_is_enabled("coefficient::arith")) {
1651 tracef("coefficient_shl() =>"); coefficient_print(ctx, S, trace_out); tracef("\n");
1652 }
1653
1654 assert(coefficient_is_normalized(ctx, S));
1655 }
1656
STAT_DECLARE(int,coefficient,shr)1657 STAT_DECLARE(int, coefficient, shr)
1658
1659 void coefficient_shr(const lp_polynomial_context_t* ctx, coefficient_t* S, const coefficient_t* C, lp_variable_t x, unsigned n) {
1660 TRACE("coefficient::arith", "coefficient_shr()\n");
1661 STAT_INCR(coefficient, shl)
1662
1663 if (trace_is_enabled("coefficient::arith")) {
1664 tracef("C = "); coefficient_print(ctx, C, trace_out); tracef("\n");
1665 tracef("x = %s\n", lp_variable_db_get_name(ctx->var_db, x));
1666 tracef("n = %u\n", n);
1667 }
1668
1669 if (n == 0) {
1670 coefficient_assign(ctx, S, C);
1671 return;
1672 }
1673
1674 if (C->type == COEFFICIENT_NUMERIC) {
1675 assert(coefficient_is_zero(ctx, C));
1676 coefficient_assign(ctx, S, C);
1677 return;
1678 }
1679
1680 assert(VAR(C) == x);
1681 assert(n + 1 <= SIZE(C));
1682
1683 if (n + 1 == SIZE(C)) {
1684 coefficient_t result;
1685 coefficient_construct_copy(ctx, &result, coefficient_lc(C));
1686 coefficient_swap(&result, S);
1687 coefficient_destruct(&result);
1688 } else {
1689 coefficient_t result;
1690 coefficient_construct_rec(ctx, &result, VAR(C), SIZE(C) - n);
1691 int i;
1692 for (i = 0; i < (int) SIZE(C) - (int) n; ++ i) {
1693 coefficient_assign(ctx, COEFF(&result, i), COEFF(C, i + n));
1694 }
1695 coefficient_swap(&result, S);
1696 coefficient_destruct(&result);
1697 }
1698
1699 if (trace_is_enabled("coefficient::arith")) {
1700 tracef("coefficient_shr() =>"); coefficient_print(ctx, S, trace_out); tracef("\n");
1701 }
1702
1703 assert(coefficient_is_normalized(ctx, S));
1704 }
1705
STAT_DECLARE(int,coefficient,pow)1706 STAT_DECLARE(int, coefficient, pow)
1707
1708 void coefficient_pow(const lp_polynomial_context_t* ctx, coefficient_t* P, const coefficient_t* C, unsigned n) {
1709 TRACE("coefficient", "coefficient_pow()\n");
1710 STAT_INCR(coefficient, pow)
1711
1712 if (trace_is_enabled("coefficient")) {
1713 tracef("P = "); coefficient_print(ctx, P, trace_out); tracef("\n");
1714 tracef("C = "); coefficient_print(ctx, C, trace_out); tracef("\n");
1715 }
1716
1717 if (n == 0) {
1718 coefficient_assign_int(ctx, P, 1);
1719 return;
1720 }
1721
1722 if (n == 1) {
1723 coefficient_assign(ctx, P, C);
1724 return;
1725 }
1726
1727 coefficient_t result, tmp;
1728
1729 switch(C->type) {
1730 case COEFFICIENT_NUMERIC:
1731 if (P->type == COEFFICIENT_POLYNOMIAL) {
1732 coefficient_construct(ctx, &result);
1733 integer_pow(ctx->K, &result.value.num, &C->value.num, n);
1734 coefficient_swap(P, &result);
1735 coefficient_destruct(&result);
1736 } else {
1737 integer_pow(ctx->K, &P->value.num, &C->value.num, n);
1738 }
1739 break;
1740 case COEFFICIENT_POLYNOMIAL:
1741 // Accumulator for C^n (start with 1)
1742 coefficient_construct_from_int(ctx, &result, 1);
1743 coefficient_ensure_capacity(ctx, &result, VAR(C), (SIZE(C)-1)*n + 1);
1744 // C^power of 2 (start with C)
1745 coefficient_construct_copy(ctx, &tmp, C);
1746 while (n) {
1747 if (n & 1) {
1748 coefficient_mul(ctx, &result, &result, &tmp);
1749 }
1750 coefficient_mul(ctx, &tmp, &tmp, &tmp);
1751 n >>= 1;
1752 }
1753 coefficient_normalize(ctx, &result);
1754 coefficient_swap(&result, P);
1755 coefficient_destruct(&tmp);
1756 coefficient_destruct(&result);
1757 break;
1758 }
1759
1760 assert(coefficient_is_normalized(ctx, P));
1761 }
1762
STAT_DECLARE(int,coefficient,add_mul)1763 STAT_DECLARE(int, coefficient, add_mul)
1764
1765 void coefficient_add_mul(const lp_polynomial_context_t* ctx, coefficient_t* S, const coefficient_t* C1, const coefficient_t* C2) {
1766 TRACE("coefficient::arith", "coefficient_add_mul()\n");
1767 STAT_INCR(coefficient, add_mul)
1768
1769 if (trace_is_enabled("coefficient::arith")) {
1770 tracef("S = "); coefficient_print(ctx, S, trace_out); tracef("\n");
1771 tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
1772 tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
1773 }
1774
1775 if (C1->type == COEFFICIENT_NUMERIC && C2->type == COEFFICIENT_NUMERIC && S->type == COEFFICIENT_NUMERIC) {
1776 integer_add_mul(ctx->K, &S->value.num, &C1->value.num, &C2->value.num);
1777 } else {
1778 coefficient_t mul;
1779 coefficient_construct(ctx, &mul);
1780 coefficient_mul(ctx, &mul, C1, C2);
1781 coefficient_add(ctx, S, S, &mul);
1782 coefficient_destruct(&mul);
1783 }
1784
1785 assert(coefficient_is_normalized(ctx, S));
1786 }
1787
STAT_DECLARE(int,coefficient,sub_mul)1788 STAT_DECLARE(int, coefficient, sub_mul)
1789
1790 void coefficient_sub_mul(const lp_polynomial_context_t* ctx, coefficient_t* S, const coefficient_t* C1, const coefficient_t* C2) {
1791 TRACE("coefficient::arith", "coefficient_sub_mul()\n");
1792 STAT_INCR(coefficient, sub_mul)
1793
1794 if (S->type == COEFFICIENT_NUMERIC && C1->type == COEFFICIENT_NUMERIC && C2->type == COEFFICIENT_NUMERIC) {
1795 integer_sub_mul(ctx->K, &S->value.num, &C1->value.num, &C2->value.num);
1796 } else {
1797 coefficient_t mul;
1798 coefficient_construct(ctx, &mul);
1799 coefficient_mul(ctx, &mul, C1, C2);
1800 coefficient_sub(ctx, S, S, &mul);
1801 coefficient_destruct(&mul);
1802 }
1803
1804 assert(coefficient_is_normalized(ctx, S));
1805 }
1806
STAT_DECLARE(int,coefficient,derivative)1807 STAT_DECLARE(int, coefficient, derivative)
1808
1809 void coefficient_derivative(const lp_polynomial_context_t* ctx, coefficient_t* C_d, const coefficient_t* C) {
1810 TRACE("coefficient", "coefficient_derivative()\n");
1811 STAT_INCR(coefficient, derivative)
1812
1813 size_t i;
1814 coefficient_t result;
1815
1816 switch(C->type) {
1817 case COEFFICIENT_NUMERIC:
1818 coefficient_construct(ctx, &result);
1819 break;
1820 case COEFFICIENT_POLYNOMIAL:
1821 coefficient_construct_rec(ctx, &result, VAR(C), SIZE(C));
1822 for (i = 1; i < SIZE(C); ++ i) {
1823 coefficient_mul_int(ctx, COEFF(&result, i-1), COEFF(C, i), i);
1824 }
1825 coefficient_normalize(ctx, &result);
1826 break;
1827 }
1828
1829 coefficient_swap(C_d, &result);
1830 coefficient_destruct(&result);
1831
1832 assert(coefficient_is_normalized(ctx, C_d));
1833 }
1834
coefficient_div_degrees(const lp_polynomial_context_t * ctx,coefficient_t * C,size_t p)1835 void coefficient_div_degrees(const lp_polynomial_context_t* ctx, coefficient_t* C, size_t p) {
1836 if (C->type == COEFFICIENT_POLYNOMIAL) {
1837 size_t i;
1838 for (i = 1; i < SIZE(C); ++ i) {
1839 if (!coefficient_is_zero(ctx, COEFF(C, i))) {
1840 assert(i % p == 0);
1841 assert(coefficient_is_zero(ctx, COEFF(C, i/ p)));
1842 coefficient_swap(COEFF(C, i), COEFF(C, i / p));
1843 }
1844 }
1845 coefficient_normalize(ctx, C);
1846 }
1847 }
1848
1849 //
1850 // Implementation of the division/reduction/gcd stuff
1851 //
1852
STAT_DECLARE(int,coefficient,reduce)1853 STAT_DECLARE(int, coefficient, reduce)
1854
1855 void coefficient_reduce(
1856 const lp_polynomial_context_t* ctx,
1857 const coefficient_t* A, const coefficient_t* B,
1858 coefficient_t* P, coefficient_t* Q, coefficient_t* R,
1859 remaindering_type_t type)
1860 {
1861 TRACE("coefficient", "coefficient_reduce()\n");
1862 STAT_INCR(coefficient, reduce)
1863
1864 if (trace_is_enabled("coefficient::reduce")) {
1865 tracef("A = "); coefficient_print(ctx, A, trace_out); tracef("\n");
1866 tracef("B = "); coefficient_print(ctx, B, trace_out); tracef("\n");
1867 }
1868
1869 assert(A->type == COEFFICIENT_POLYNOMIAL);
1870
1871 // Start with R = A, P = 1
1872 coefficient_t P_tmp, Q_tmp, R_tmp;
1873 if (P) {
1874 coefficient_construct_from_int(ctx, &P_tmp, 1);
1875 }
1876 if (Q) {
1877 coefficient_construct(ctx, &Q_tmp);
1878 }
1879 coefficient_construct_copy(ctx, &R_tmp, A);
1880
1881 // The main variable we are reducing
1882 lp_variable_t x = VAR(A);
1883
1884 // Keep track of the degrees of the reduct to account for dense/sparse
1885 int R_deg = coefficient_degree(&R_tmp);
1886 int R_deg_prev = R_deg;
1887 int B_deg = coefficient_degree_safe(ctx, B, x);
1888
1889 // Temporaries for computation
1890 coefficient_t lcm, r, b;
1891 coefficient_construct(ctx, &lcm);
1892 coefficient_construct_from_int(ctx, &r, 1);
1893 coefficient_construct(ctx, &b);
1894
1895 // Invariant:
1896 //
1897 // P*A = Q*B + R
1898 //
1899 // initially
1900 //
1901 // P = 1, Q = 0, R = A
1902 //
1903 do {
1904
1905 // Leading coefficient of B
1906 const coefficient_t* lc_B = coefficient_lc_safe(ctx, B, x);
1907
1908 // Account for the sparse operation
1909 switch (type) {
1910 case REMAINDERING_PSEUDO_DENSE:
1911 if (R_deg_prev - R_deg > 1) {
1912 // Multiply with the missed power of lc(B)
1913 int missed = R_deg < B_deg ?
1914 R_deg_prev - B_deg : R_deg_prev - R_deg - 1;
1915 assert(missed >= 0);
1916 if (missed > 0) {
1917 coefficient_t pow;
1918 coefficient_construct(ctx, &pow);
1919 coefficient_pow(ctx, &pow, lc_B, missed);
1920 if (P) {
1921 coefficient_mul(ctx, &P_tmp, &P_tmp, &pow);
1922 }
1923 if (Q) {
1924 coefficient_mul(ctx, &Q_tmp, &Q_tmp, &pow);
1925 }
1926 coefficient_mul(ctx, &R_tmp, &R_tmp, &pow);
1927 coefficient_destruct(&pow);
1928 }
1929 }
1930 break;
1931 default:
1932 break;
1933 }
1934
1935 // If we eliminated all of x we are done
1936 if (coefficient_is_zero(ctx, &R_tmp) || R_deg < B_deg) {
1937 break;
1938 }
1939
1940 // How much x do we need to supply
1941 int d = R_deg - B_deg;
1942
1943 // Eliminate the coefficient. We have
1944 //
1945 // R = r_i * x^i + red(R)
1946 // B = b_j * x^j + ref(B)
1947 //
1948 // and we compute
1949 //
1950 // lcm = lcm(r_i, b_j)
1951 // r = lcm/r_i
1952 // b = (lcm/b_j)*x^d
1953 //
1954 // R' = r*R - b*x^d*B
1955 //
1956 // thus eliminating the i-th power. As a consequence, the invariant is
1957 //
1958 // P*A = Q*B + R [*b]
1959 //
1960 // r*P*A = r*Q*B + rR = r*Q*B + R' + b*B
1961 // r*P*A = (r*Q + b)B + R'
1962 //
1963 // so we set (making sure that b is positive)
1964 //
1965 // P' = b*P
1966 // Q' = b*Q + c
1967 //
1968
1969 // Leading coefficient of R
1970 const coefficient_t* lc_R = coefficient_lc_safe(ctx, &R_tmp, x);
1971
1972 switch (type) {
1973 case REMAINDERING_EXACT_SPARSE:
1974 {
1975 // If we are exact, we assume that b_j divides r_i
1976 coefficient_div(ctx, &b, lc_R, lc_B);
1977 break;
1978 }
1979 case REMAINDERING_LCM_SPARSE:
1980 {
1981 // a, b, c
1982 coefficient_lcm(ctx, &lcm, lc_R, lc_B);
1983 coefficient_div(ctx, &r, &lcm, lc_R);
1984 coefficient_div(ctx, &b, &lcm, lc_B);
1985 if (coefficient_lc_sgn(ctx, &r) < 0) {
1986 coefficient_neg(ctx, &r, &r);
1987 coefficient_neg(ctx, &b, &b);
1988 }
1989 break;
1990 }
1991 case REMAINDERING_PSEUDO_DENSE:
1992 case REMAINDERING_PSEUDO_SPARSE:
1993 {
1994 coefficient_assign(ctx, &r, lc_B);
1995 coefficient_assign(ctx, &b, lc_R);
1996 break;
1997 }
1998 default:
1999 assert(0);
2000 }
2001
2002 // Add the power of x
2003 coefficient_shl(ctx, &b, &b, x, d);
2004
2005 assert(!coefficient_is_zero(ctx, &r));
2006 assert(!coefficient_is_zero(ctx, &b));
2007
2008 if (trace_is_enabled("coefficient::reduce")) {
2009 tracef("lcm = "); coefficient_print(ctx, &lcm, trace_out); tracef("\n");
2010 tracef("R = "); coefficient_print(ctx, &R_tmp, trace_out); tracef("\n");
2011 tracef("r = "); coefficient_print(ctx, &r, trace_out); tracef("\n");
2012 tracef("B = "); coefficient_print(ctx, B, trace_out); tracef("\n");
2013 tracef("b = "); coefficient_print(ctx, &b, trace_out); tracef("\n");
2014 }
2015
2016 // R'
2017 coefficient_mul(ctx, &R_tmp, &R_tmp, &r);
2018 coefficient_sub_mul(ctx, &R_tmp, B, &b);
2019
2020 if (trace_is_enabled("coefficient::reduce")) {
2021 tracef("R' = "); coefficient_print(ctx, &R_tmp, trace_out); tracef("\n");
2022 }
2023
2024 // Update the degrees of R
2025 R_deg_prev = R_deg;
2026 R_deg = coefficient_degree_safe(ctx, &R_tmp, x);
2027 assert(coefficient_is_zero(ctx, &R_tmp) || R_deg < R_deg_prev);
2028
2029
2030 // P' (if needed)
2031 if (P) {
2032 coefficient_mul(ctx, &P_tmp, &P_tmp, &r);
2033 }
2034
2035 // Q' (if needed)
2036 if (Q) {
2037 coefficient_mul(ctx, &Q_tmp, &Q_tmp, &r);
2038 coefficient_add(ctx, &Q_tmp, &Q_tmp, &b);
2039 }
2040
2041 if (trace_is_enabled("coefficient::reduce")) {
2042 if (P) {
2043 tracef("P = "); coefficient_print(ctx, &P_tmp, trace_out); tracef("\n");
2044 }
2045 if (Q) {
2046 tracef("Q = "); coefficient_print(ctx, &Q_tmp, trace_out); tracef("\n");
2047 }
2048 }
2049
2050 } while (1);
2051
2052 // Move the result out, and remove temps
2053 if (P) {
2054 coefficient_swap(P, &P_tmp);
2055 coefficient_destruct(&P_tmp);
2056 assert(coefficient_is_normalized(ctx, P));
2057 }
2058 if (Q) {
2059 coefficient_swap(Q, &Q_tmp);
2060 coefficient_destruct(&Q_tmp);
2061 assert(coefficient_is_normalized(ctx, Q));
2062 }
2063 if (R) {
2064 coefficient_swap(R, &R_tmp);
2065 assert(coefficient_is_normalized(ctx, R));
2066 }
2067 coefficient_destruct(&R_tmp);
2068
2069 coefficient_destruct(&lcm);
2070 coefficient_destruct(&r);
2071 coefficient_destruct(&b);
2072 }
2073
2074
coefficient_div_constant(const lp_polynomial_context_t * ctx,coefficient_t * C,const lp_integer_t * A)2075 void coefficient_div_constant(const lp_polynomial_context_t* ctx, coefficient_t* C, const lp_integer_t* A) {
2076
2077 size_t i ;
2078
2079 if (C->type == COEFFICIENT_NUMERIC) {
2080 integer_div_Z(&C->value.num, &C->value.num, A);
2081 } else {
2082 for (i = 0; i < SIZE(C); ++ i) {
2083 coefficient_div_constant(ctx, COEFF(C, i), A);
2084 }
2085 }
2086 }
2087
STAT_DECLARE(int,coefficient,div)2088 STAT_DECLARE(int, coefficient, div)
2089
2090 void coefficient_div(const lp_polynomial_context_t* ctx, coefficient_t* D, const coefficient_t* C1, const coefficient_t* C2) {
2091 TRACE("coefficient", "coefficient_div()\n");
2092 STAT_INCR(coefficient, div)
2093
2094 // 0/C2 = 0
2095 if (coefficient_is_zero(ctx, C1)) {
2096 coefficient_assign(ctx, D, C1);
2097 return;
2098 }
2099
2100 // C1/C1 = 1
2101 if (coefficient_cmp(ctx, C1, C2) == 0) {
2102 coefficient_assign_int(ctx, D, 1);
2103 return;
2104 }
2105
2106 // Special case for constants
2107 if (coefficient_is_constant(C2)) {
2108 coefficient_assign(ctx, D, C1);
2109 coefficient_div_constant(ctx, D, &C2->value.num);
2110 return;
2111 }
2112
2113 // If different variables
2114 if (VAR(C1) != VAR(C2)) {
2115 coefficient_t result;
2116 coefficient_construct_rec(ctx, &result, VAR(C1), SIZE(C1));
2117 size_t i;
2118 for (i = 0; i < SIZE(C1); ++ i) {
2119 coefficient_div(ctx, COEFF(&result, i), COEFF(C1, i), C2);
2120 }
2121 coefficient_swap(&result, D);
2122 coefficient_destruct(&result);
2123 return;
2124 }
2125
2126 // Both polynomials in the same variables, check if we can divide by x^k
2127 size_t i = 0;
2128 while (coefficient_is_zero(ctx, COEFF(C2, i)) && coefficient_is_zero(ctx, COEFF(C1, i))) {
2129 ++ i;
2130 }
2131 if (i > 0) {
2132 // i = first non-zero coefficient, shift by i
2133 coefficient_t C1_shifted, C2_shifted;
2134 lp_variable_t x = VAR(C2);
2135 coefficient_construct(ctx, &C1_shifted);
2136 coefficient_construct(ctx, &C2_shifted);
2137 coefficient_shr(ctx, &C1_shifted, C1, x, i);
2138 coefficient_shr(ctx, &C2_shifted, C2, x, i);
2139 coefficient_div(ctx, D, &C1_shifted, &C2_shifted);
2140 coefficient_destruct(&C1_shifted);
2141 coefficient_destruct(&C2_shifted);
2142 return;
2143 }
2144
2145 if (trace_is_enabled("coefficient") || trace_is_enabled("coefficient::div")) {
2146 tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
2147 tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
2148 }
2149
2150 assert(!coefficient_is_zero(ctx, C2));
2151
2152 if (trace_is_enabled("coefficient::check_division")) {
2153 coefficient_t R;
2154 coefficient_construct(ctx, &R);
2155 coefficient_reduce(ctx, C1, C2, 0, D, &R, REMAINDERING_EXACT_SPARSE);
2156 if (!coefficient_is_zero(ctx, &R)) {
2157 tracef("WRONG DIVISION!\n");
2158 tracef("P = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
2159 tracef("Q = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
2160 }
2161 coefficient_destruct(&R);
2162 } else {
2163 coefficient_reduce(ctx, C1, C2, 0, D, 0, REMAINDERING_EXACT_SPARSE);
2164 }
2165
2166 if (trace_is_enabled("coefficient")) {
2167 tracef("coefficient_div() => "); coefficient_print(ctx, D, trace_out); tracef("\n");
2168 }
2169
2170 assert(coefficient_is_normalized(ctx, D));
2171 }
2172
STAT_DECLARE(int,coefficient,rem)2173 STAT_DECLARE(int, coefficient, rem)
2174
2175 void coefficient_rem(const lp_polynomial_context_t* ctx, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2) {
2176 TRACE("coefficient", "coefficient_rem()\n");
2177 STAT_INCR(coefficient, rem)
2178
2179 if (trace_is_enabled("coefficient")) {
2180 tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
2181 tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
2182 }
2183
2184 assert(!coefficient_is_zero(ctx, C2));
2185
2186 int cmp_type = coefficient_cmp_type(ctx, C1, C2);
2187
2188 assert(cmp_type >= 0);
2189
2190 if (cmp_type == 0 && C1->type == COEFFICIENT_NUMERIC) {
2191 assert(C2->type == COEFFICIENT_NUMERIC);
2192 if (R->type == COEFFICIENT_POLYNOMIAL) {
2193 coefficient_destruct(R);
2194 coefficient_construct(ctx, R);
2195 }
2196 integer_rem_Z(&R->value.num, &C1->value.num, &C2->value.num);
2197 } else {
2198 coefficient_reduce(ctx, C1, C2, 0, 0, R, REMAINDERING_EXACT_SPARSE);
2199 }
2200
2201 if (trace_is_enabled("coefficient")) {
2202 tracef("coefficient_rem() => "); coefficient_print(ctx, R, trace_out); tracef("\n");
2203 }
2204
2205 assert(coefficient_is_normalized(ctx, R));
2206 }
2207
STAT_DECLARE(int,coefficient,sprem)2208 STAT_DECLARE(int, coefficient, sprem)
2209
2210 void coefficient_sprem(const lp_polynomial_context_t* ctx, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2) {
2211 TRACE("coefficient", "coefficient_sprem()\n");
2212 STAT_INCR(coefficient, sprem)
2213
2214 if (trace_is_enabled("coefficient")) {
2215 tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
2216 tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
2217 }
2218
2219 assert(!coefficient_is_zero(ctx, C2));
2220
2221 int cmp_type = coefficient_cmp_type(ctx, C1, C2);
2222
2223 assert(cmp_type >= 0);
2224
2225 if (cmp_type == 0 && C1->type == COEFFICIENT_NUMERIC) {
2226 assert(C2->type == COEFFICIENT_NUMERIC);
2227 if (R->type == COEFFICIENT_POLYNOMIAL) {
2228 coefficient_destruct(R);
2229 coefficient_construct(ctx, R);
2230 }
2231 integer_rem_Z(&R->value.num, &C1->value.num, &C2->value.num);
2232 } else {
2233 coefficient_reduce(ctx, C1, C2, 0, 0, R, REMAINDERING_PSEUDO_SPARSE);
2234 }
2235
2236 if (trace_is_enabled("coefficient")) {
2237 tracef("coefficient_sprem() => "); coefficient_print(ctx, R, trace_out); tracef("\n");
2238 }
2239
2240 assert(coefficient_is_normalized(ctx, R));
2241 }
2242
STAT_DECLARE(int,coefficient,prem)2243 STAT_DECLARE(int, coefficient, prem)
2244
2245 void coefficient_prem(const lp_polynomial_context_t* ctx, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2) {
2246 TRACE("coefficient", "coefficient_prem()\n");
2247 STAT_INCR(coefficient, rem)
2248
2249 if (trace_is_enabled("coefficient")) {
2250 tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
2251 tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
2252 }
2253
2254 assert(!coefficient_is_zero(ctx, C2));
2255
2256 int cmp_type = coefficient_cmp_type(ctx, C1, C2);
2257
2258 assert(cmp_type >= 0);
2259
2260 if (cmp_type == 0 && C1->type == COEFFICIENT_NUMERIC) {
2261 assert(C2->type == COEFFICIENT_NUMERIC);
2262 if (R->type == COEFFICIENT_POLYNOMIAL) {
2263 coefficient_destruct(R);
2264 coefficient_construct(ctx, R);
2265 }
2266 integer_rem_Z(&R->value.num, &C1->value.num, &C2->value.num);
2267 } else {
2268 coefficient_reduce(ctx, C1, C2, 0, 0, R, REMAINDERING_PSEUDO_DENSE);
2269 }
2270
2271 if (trace_is_enabled("coefficient")) {
2272 tracef("coefficient_prem() => "); coefficient_print(ctx, R, trace_out); tracef("\n");
2273 }
2274
2275 assert(coefficient_is_normalized(ctx, R));
2276 }
2277
2278
STAT_DECLARE(int,coefficient,divrem)2279 STAT_DECLARE(int, coefficient, divrem)
2280
2281 void coefficient_divrem(const lp_polynomial_context_t* ctx, coefficient_t* D, coefficient_t* R, const coefficient_t* C1, const coefficient_t* C2) {
2282 TRACE("coefficient", "coefficient_divrem()\n");
2283 STAT_INCR(coefficient, divrem)
2284
2285 if (trace_is_enabled("coefficient")) {
2286 tracef("C1 = "); coefficient_print(ctx, C1, trace_out); tracef("\n");
2287 tracef("C2 = "); coefficient_print(ctx, C2, trace_out); tracef("\n");
2288 }
2289
2290 assert(!coefficient_is_zero(ctx, C2));
2291
2292 int cmp_type = coefficient_cmp_type(ctx, C1, C2);
2293
2294 assert(cmp_type >= 0);
2295
2296 if (cmp_type == 0) {
2297 switch(C1->type) {
2298 case COEFFICIENT_NUMERIC:
2299 assert(C2->type == COEFFICIENT_NUMERIC);
2300 if (R->type == COEFFICIENT_POLYNOMIAL) {
2301 coefficient_destruct(R);
2302 coefficient_construct(ctx, R);
2303 }
2304 integer_rem_Z(&R->value.num, &C1->value.num, &C2->value.num);
2305 break;
2306 case COEFFICIENT_POLYNOMIAL:
2307 {
2308 coefficient_reduce(ctx, C1, C2, 0, D, R, REMAINDERING_EXACT_SPARSE);
2309 break;
2310 }
2311 default:
2312 assert(0);
2313 }
2314 } else {
2315 // Just use the regular methods
2316 coefficient_rem(ctx, R, COEFF(C1, 0), C2);
2317 coefficient_div(ctx, D, C1, C2);
2318 }
2319
2320 if (trace_is_enabled("coefficient")) {
2321 tracef("coefficient_divrem() => \n");
2322 tracef("D = "); coefficient_print(ctx, D, trace_out); tracef("\n");
2323 tracef("D = "); coefficient_print(ctx, R, trace_out); tracef("\n");
2324 }
2325
2326 assert(coefficient_is_normalized(ctx, R));
2327 }
2328
STAT_DECLARE(int,coefficient,divides)2329 STAT_DECLARE(int, coefficient, divides)
2330
2331 int coefficient_divides(const lp_polynomial_context_t* ctx, const coefficient_t* C1, const coefficient_t* C2) {
2332 TRACE("coefficient", "coefficient_divides()\n");
2333 STAT_INCR(coefficient, divides)
2334
2335 coefficient_t R;
2336 coefficient_construct(ctx, &R);
2337 coefficient_prem(ctx, &R, C2, C1);
2338 int divides = coefficient_is_zero(ctx, &R);
2339 coefficient_destruct(&R);
2340
2341 return divides;
2342 }
2343
coefficient_is_univariate(const coefficient_t * C)2344 int coefficient_is_univariate(const coefficient_t* C) {
2345 size_t i;
2346 if (C->type == COEFFICIENT_NUMERIC) {
2347 return 1;
2348 } else {
2349 for (i = 0; i < SIZE(C); ++ i) {
2350 if (COEFF(C, i)->type != COEFFICIENT_NUMERIC) {
2351 return 0;
2352 }
2353 }
2354 return 1;
2355 }
2356 }
2357
coefficient_is_linear(const coefficient_t * C)2358 int coefficient_is_linear(const coefficient_t* C) {
2359 if (C->type != COEFFICIENT_POLYNOMIAL) {
2360 return 0;
2361 }
2362 while (C->type == COEFFICIENT_POLYNOMIAL && coefficient_degree(C) == 1 && coefficient_lc(C)->type == COEFFICIENT_NUMERIC) {
2363 C = COEFF(C, 0);
2364 }
2365 return (C->type == COEFFICIENT_NUMERIC);
2366 }
2367
coefficient_get_constant(const coefficient_t * C)2368 const lp_integer_t* coefficient_get_constant(const coefficient_t* C) {
2369 while (C->type == COEFFICIENT_POLYNOMIAL) {
2370 C = COEFF(C, 0);
2371 }
2372 return &C->value.num;
2373 }
2374
coefficient_to_univariate(const lp_polynomial_context_t * ctx,const coefficient_t * C)2375 lp_upolynomial_t* coefficient_to_univariate(const lp_polynomial_context_t* ctx, const coefficient_t* C) {
2376 assert(C->type == COEFFICIENT_POLYNOMIAL);
2377
2378 lp_integer_t* coeff = malloc(sizeof(lp_integer_t)*SIZE(C));
2379
2380 size_t i;
2381 for (i = 0; i < SIZE(C); ++ i) {
2382 integer_construct_copy(ctx->K, coeff + i, coefficient_get_constant(COEFF(C, i)));
2383 }
2384
2385 lp_upolynomial_t* C_u = lp_upolynomial_construct(ctx->K, SIZE(C) - 1, coeff);
2386
2387 for (i = 0; i < SIZE(C); ++ i) {
2388 integer_destruct(coeff + i);
2389 }
2390 free(coeff);
2391
2392 return C_u;
2393 }
2394
STAT_DECLARE(int,coefficient,resultant)2395 STAT_DECLARE(int, coefficient, resultant)
2396
2397 void coefficient_resultant(const lp_polynomial_context_t* ctx, coefficient_t* res, const coefficient_t* A, const coefficient_t* B) {
2398
2399 if (trace_is_enabled("coefficient")) {
2400 tracef("coefficient_resultant("); coefficient_print(ctx, A, trace_out); tracef(", "); coefficient_print(ctx, B, trace_out); tracef(")\n");
2401 }
2402 STAT_INCR(coefficient, resultant)
2403
2404 if (trace_is_enabled("coefficient")) {
2405 tracef("A = "); coefficient_print(ctx, A, trace_out); tracef("\n");
2406 tracef("B = "); coefficient_print(ctx, B, trace_out); tracef("\n");
2407 }
2408
2409 assert(A->type == COEFFICIENT_POLYNOMIAL);
2410 assert(B->type == COEFFICIENT_POLYNOMIAL);
2411
2412 assert(VAR(B) == VAR(A));
2413
2414 size_t A_deg = coefficient_degree(A);
2415 size_t B_deg = coefficient_degree(B);
2416
2417 if (A_deg < B_deg) {
2418 coefficient_resultant(ctx, res, B, A);
2419 if ((A_deg % 2) && (B_deg % 2)) {
2420 coefficient_neg(ctx, res, res);
2421 }
2422 return;
2423 }
2424
2425 // Compute the PSC
2426 size_t psc_size = B_deg + 1;
2427 coefficient_t* psc = malloc(sizeof(coefficient_t)*psc_size);
2428 size_t i;
2429 for (i = 0; i < psc_size; ++ i) {
2430 coefficient_construct(ctx, psc + i);
2431 }
2432 coefficient_psc(ctx, psc, A, B);
2433
2434 // Resultant is PSC[0]
2435 coefficient_swap(res, psc);
2436
2437 // Remove temps
2438 for (i = 0; i < psc_size; ++ i) {
2439 coefficient_destruct(psc + i);
2440 }
2441 free(psc);
2442
2443 if (trace_is_enabled("coefficient")) {
2444 tracef("coefficient_resultant() => "); coefficient_print(ctx, res, trace_out); tracef("\n");
2445 }
2446
2447 }
2448
2449 ///
2450 /// Normalization that everyone is using
2451 ///
2452
2453 static void
coefficient_normalize(const lp_polynomial_context_t * ctx,coefficient_t * C)2454 coefficient_normalize(const lp_polynomial_context_t* ctx, coefficient_t* C) {
2455 if (C->type == COEFFICIENT_POLYNOMIAL) {
2456 assert(C->value.rec.size >= 1);
2457 size_t i = C->value.rec.size - 1;
2458 // Find the first non-zero coefficient
2459 while (i > 0 && coefficient_is_zero(ctx, COEFF(C, i))) {
2460 i --;
2461 }
2462 // If a constant, just upgrade it (it could be an actual constant or just
2463 // some other polynomial)
2464 if (i == 0) {
2465 coefficient_t result;
2466 coefficient_construct(ctx, &result);
2467 coefficient_swap(&result, COEFF(C, 0));
2468 coefficient_swap(&result, C);
2469 coefficient_destruct(&result);
2470 } else {
2471 // We are a proper polynomial, set the size
2472 C->value.rec.size = i + 1;
2473 }
2474 }
2475 }
2476
2477 static void
coefficient_normalize_m(const lp_polynomial_context_t * ctx,coefficient_t * C,const lp_assignment_t * m)2478 coefficient_normalize_m(const lp_polynomial_context_t* ctx, coefficient_t* C, const lp_assignment_t* m) {
2479 if (C->type == COEFFICIENT_POLYNOMIAL) {
2480 assert(C->value.rec.size >= 1);
2481 int i = C->value.rec.size - 1;
2482 // Find the first non-zero coefficient
2483 while (i >= 0 && coefficient_sgn(ctx, COEFF(C, i), m) == 0) {
2484 coefficient_assign_int(ctx, COEFF(C, i), 0);
2485 i --;
2486 }
2487 if (i < 0) {
2488 // Normalizes to 0
2489 coefficient_assign_int(ctx, C, 0);
2490 } else if (i == 0) {
2491 // If a constant, just upgrade it (it could be an actual constant or just
2492 // some other polynomial)
2493 coefficient_t result;
2494 coefficient_construct(ctx, &result);
2495 coefficient_swap(&result, COEFF(C, 0));
2496 coefficient_swap(&result, C);
2497 coefficient_destruct(&result);
2498 } else {
2499 // We are a proper polynomial, set the size
2500 C->value.rec.size = i + 1;
2501 }
2502 }
2503 }
2504
2505 int
coefficient_is_normalized(const lp_polynomial_context_t * ctx,coefficient_t * C)2506 coefficient_is_normalized(const lp_polynomial_context_t* ctx, coefficient_t* C) {
2507 switch (C->type) {
2508 case COEFFICIENT_NUMERIC:
2509 return 1;
2510 break;
2511 case COEFFICIENT_POLYNOMIAL:
2512 if (SIZE(C) <= 1) {
2513 return 0;
2514 }
2515 if (coefficient_is_zero(ctx, COEFF(C, SIZE(C) - 1))) {
2516 return 0;
2517 }
2518 return 1;
2519 break;
2520 default:
2521 assert(0);
2522 break;
2523 }
2524 return 0;
2525 }
2526
2527 static void
coefficient_ensure_capacity(const lp_polynomial_context_t * ctx,coefficient_t * C,lp_variable_t x,size_t capacity)2528 coefficient_ensure_capacity(const lp_polynomial_context_t* ctx, coefficient_t* C, lp_variable_t x, size_t capacity) {
2529 assert(capacity >= 1);
2530 coefficient_t tmp;
2531 switch(C->type) {
2532 case COEFFICIENT_NUMERIC:
2533 // Create a new recursive and put the current constant into constants place
2534 coefficient_construct_rec(ctx, &tmp, x, capacity);
2535 coefficient_swap(COEFF(&tmp, 0), C);
2536 coefficient_swap(C, &tmp);
2537 coefficient_destruct(&tmp);
2538 break;
2539 case COEFFICIENT_POLYNOMIAL:
2540 if (x != VAR(C)) {
2541 assert(lp_variable_order_cmp(ctx->var_order, x, VAR(C)) > 0);
2542 // Same as for constants above
2543 coefficient_construct_rec(ctx, &tmp, x, capacity);
2544 coefficient_swap(COEFF(&tmp, 0), C);
2545 coefficient_swap(C, &tmp);
2546 coefficient_destruct(&tmp);
2547 } else if (capacity > C->value.rec.capacity) {
2548 // Already recursive polynomial, resize up
2549 C->value.rec.coefficients = realloc(C->value.rec.coefficients, capacity * sizeof(coefficient_t));
2550 size_t i;
2551 for (i = C->value.rec.capacity; i < capacity; ++ i) {
2552 coefficient_construct(ctx, C->value.rec.coefficients + i);
2553 }
2554 C->value.rec.capacity = capacity;
2555 C->value.rec.size = capacity;
2556 }
2557 break;
2558 }
2559 }
2560
coefficient_evaluate_rationals(const lp_polynomial_context_t * ctx,const coefficient_t * C,const lp_assignment_t * M,coefficient_t * C_out,lp_integer_t * multiplier)2561 void coefficient_evaluate_rationals(const lp_polynomial_context_t* ctx, const coefficient_t* C, const lp_assignment_t* M, coefficient_t* C_out, lp_integer_t* multiplier) {
2562
2563 assert(multiplier);
2564 assert(ctx->K == lp_Z);
2565
2566 size_t i;
2567 lp_variable_t x;
2568
2569 // Start wit multiplier 1
2570 integer_assign_int(lp_Z, multiplier, 1);
2571
2572 if (C->type == COEFFICIENT_NUMERIC) {
2573 // Just a number, we're done
2574 coefficient_assign(ctx, C_out, C);
2575 } else {
2576 assert(C->type == COEFFICIENT_POLYNOMIAL);
2577
2578 // Temp for the result
2579 coefficient_t result;
2580
2581 // Get the variable and it's value, if any
2582 x = VAR(C);
2583 const lp_value_t* x_value = lp_assignment_get_value(M, x);
2584
2585 // The degree of the polynomial
2586 size_t size = SIZE(C);
2587
2588 // Check if the value is rational and we can substitute it
2589 if (!lp_value_is_rational(x_value))
2590 {
2591 // We can not substitute so
2592 //
2593 // C = a_n * x^n + ... + a_1 * x + a_0
2594 //
2595 // We substitute in all a_n obtaining a_k = b_n / m_k, m = lcm(m_1, ..., m_n)
2596 //
2597 // m * c = sum b_k * x^k * m / m_k
2598
2599 coefficient_construct_rec(ctx, &result, x, size);
2600
2601 // Compute the evaluation of the coefficients
2602 lp_integer_t* m = malloc(sizeof(coefficient_t)*size);
2603 for (i = 0; i < size; ++ i) {
2604 integer_construct(m + i);
2605 coefficient_evaluate_rationals(ctx, COEFF(C, i), M, COEFF(&result, i), m + i);
2606 }
2607
2608 // Compute the lcm of the m's
2609 lp_integer_assign(lp_Z, multiplier, m);
2610 for (i = 1; i < size; ++ i) {
2611 integer_lcm_Z(multiplier, multiplier, m + i);
2612 }
2613
2614 // Sum up
2615 lp_integer_t tmp;
2616 integer_construct(&tmp);
2617 for (i = 0; i < size; ++ i) {
2618 // m / m_k
2619 integer_div_exact(lp_Z, &tmp, multiplier, m + i);
2620 // b_i = b_i * R
2621 coefficient_mul_integer(ctx, COEFF(&result, i), COEFF(&result, i), &tmp);
2622 }
2623 integer_destruct(&tmp);
2624
2625 // Remove the temps
2626 for (i = 0; i < size; ++ i) {
2627 integer_destruct(m + i);
2628 }
2629 free(m);
2630
2631 } else {
2632
2633 coefficient_construct(ctx, &result);
2634
2635 // We have a value value = p/q
2636 lp_integer_t p, q;
2637 integer_construct(&p);
2638 integer_construct(&q);
2639 lp_value_get_num(x_value, &p);
2640 lp_value_get_den(x_value, &q);
2641
2642 // If we can substitute then
2643 //
2644 // C = a_n * (p/q)^n + ... + a_1 * (p/q) + a_0
2645 //
2646 // We substitute in all a_n obtaining a_k = b_n / m_k and get
2647 //
2648 // q^n * C = b_n * (p^n/m_n) + ... + b_1 * (p*q^n-1/m_1) + b_0 * q^n/m_0
2649 //
2650 // We get the m = lcm(m_1, ..., m_n) and get
2651 //
2652 // q^n * m * c = sum b_k * p^k * q^(n-k) * m / m_k
2653
2654
2655 // Compute the evaluation of the coefficients
2656 coefficient_t* b = malloc(sizeof(coefficient_t)*size);
2657 lp_integer_t* m = malloc(sizeof(lp_integer_t)*size);
2658 for (i = 0; i < size; ++ i) {
2659 coefficient_construct(ctx, b + i);
2660 integer_construct(m + i);
2661 coefficient_evaluate_rationals(ctx, COEFF(C, i), M, b + i, m + i);
2662 }
2663
2664 // Compute the lcm of the m's
2665 lp_integer_t m_lcm;
2666 lp_integer_construct_copy(lp_Z, &m_lcm, m);
2667 for (i = 1; i < size; ++ i) {
2668 integer_lcm_Z(&m_lcm, &m_lcm, m + i);
2669 }
2670
2671 // The powers
2672 lp_integer_t p_power, q_power;
2673 integer_construct_from_int(lp_Z, &p_power, 1);
2674 integer_construct(&q_power);
2675 integer_pow(lp_Z, &q_power, &q, size-1);
2676
2677 // Set the multiplier
2678 integer_mul(lp_Z, multiplier, &q_power, &m_lcm);
2679
2680 // Sum up
2681 lp_integer_t R;
2682 integer_construct(&R);
2683 for (i = 0; i < size; ++ i) {
2684 if (i) {
2685 // Update powers
2686 integer_mul(lp_Z, &p_power, &p_power, &p);
2687 integer_div_exact(lp_Z, &q_power, &q_power, &q);
2688 }
2689 // R = p^i * q^(n-i) * m / m_k
2690 integer_div_exact(lp_Z, &R, &m_lcm, m + i);
2691 integer_mul(lp_Z, &R, &R, &p_power);
2692 integer_mul(lp_Z, &R, &R, &q_power);
2693 // b_i = b_i * R
2694 coefficient_mul_integer(ctx, b + i, b + i, &R);
2695 // Add it
2696 coefficient_add(ctx, &result, &result, b + i);
2697 }
2698 integer_destruct(&R);
2699
2700 // Remove the temps
2701 for (i = 0; i < size; ++ i) {
2702 coefficient_destruct(b + i);
2703 integer_destruct(m + i);
2704 }
2705 free(b);
2706 free(m);
2707 integer_destruct(&m_lcm);
2708 integer_destruct(&p);
2709 integer_destruct(&q);
2710 integer_destruct(&p_power);
2711 integer_destruct(&q_power);
2712 }
2713
2714 // Finish up
2715 coefficient_normalize(ctx, &result);
2716 coefficient_swap(&result, C_out);
2717 coefficient_destruct(&result);
2718 }
2719
2720 assert(integer_sgn(lp_Z, multiplier) > 0);
2721
2722 }
2723
coefficient_get_variables(const coefficient_t * C,lp_variable_list_t * vars)2724 void coefficient_get_variables(const coefficient_t* C, lp_variable_list_t* vars) {
2725 if (C->type != COEFFICIENT_NUMERIC) {
2726 // Add the variable, if not already there
2727 lp_variable_t x = VAR(C);
2728 if (lp_variable_list_index(vars, x) < 0) {
2729 lp_variable_list_push(vars, x);
2730 }
2731 // Add children
2732 size_t i;
2733 for (i = 0; i < SIZE(C); ++ i) {
2734 coefficient_get_variables(COEFF(C, i), vars);
2735 }
2736 }
2737 }
2738
2739 /**
2740 * Isolate out the roots of a univariate polynomial.
2741 */
coefficient_roots_isolate_univariate(const lp_polynomial_context_t * ctx,const coefficient_t * A,lp_value_t * roots,size_t * roots_size)2742 void coefficient_roots_isolate_univariate(const lp_polynomial_context_t* ctx, const coefficient_t* A, lp_value_t* roots, size_t* roots_size) {
2743
2744 if (trace_is_enabled("coefficient::roots")) {
2745 tracef("coefficient_roots_isolate(): univariate, root finding\n");
2746 tracef("coefficient_roots_isolate(): A = "); coefficient_print(ctx, A, trace_out); tracef("\n");
2747 }
2748
2749 assert(coefficient_is_univariate(A));
2750
2751 // Special case for linear polynomials
2752 if (coefficient_degree(A) == 1) {
2753 // A_rat = ax + b => root = -b/a
2754 *roots_size = 1;
2755 lp_rational_t root;
2756 rational_construct_from_div(&root, &COEFF(A, 0)->value.num, &COEFF(A, 1)->value.num);
2757 rational_neg(&root, &root);
2758 lp_value_construct(roots, LP_VALUE_RATIONAL, &root);
2759 rational_destruct(&root);
2760 } else {
2761 // Get the univariate version
2762 lp_upolynomial_t* A_u = coefficient_to_univariate(ctx, A);
2763 assert(lp_upolynomial_degree(A_u) == coefficient_degree(A));
2764 // Make space for the algebraic numbers
2765 lp_algebraic_number_t* algebraic_roots = malloc(lp_upolynomial_degree(A_u) * sizeof(lp_algebraic_number_t));
2766
2767 if (trace_is_enabled("coefficient::roots")) {
2768 tracef("coefficient_roots_isolate(): A_u = "); lp_upolynomial_print(A_u, trace_out); tracef("\n");
2769 }
2770 // Isolate
2771 lp_upolynomial_roots_isolate(A_u, algebraic_roots, roots_size);
2772 assert(*roots_size <= coefficient_degree(A));
2773
2774 // Copy over the roots
2775 size_t i;
2776 for (i = 0; i < *roots_size; ++ i) {
2777 lp_value_construct(roots + i, LP_VALUE_ALGEBRAIC, algebraic_roots + i);
2778 lp_algebraic_number_destruct(algebraic_roots + i);
2779 }
2780 // Free the temps
2781 free(algebraic_roots);
2782 lp_upolynomial_delete(A_u);
2783 }
2784 }
2785
2786 /**
2787 * Resolve out the algebraic numbers assignment.
2788 */
2789 static
coefficient_resolve_algebraic(const lp_polynomial_context_t * ctx,const coefficient_t * A,const lp_assignment_t * m,coefficient_t * A_alg)2790 void coefficient_resolve_algebraic(const lp_polynomial_context_t* ctx, const coefficient_t* A, const lp_assignment_t* m, coefficient_t* A_alg) {
2791
2792 coefficient_assign(ctx, A_alg, A);
2793
2794 // Eliminate all variables
2795 for (;;) {
2796
2797 if (trace_is_enabled("coefficient")) {
2798 tracef("coefficient_resolve_algebraic(): A_alg = "); coefficient_print(ctx, A_alg, trace_out); tracef("\n");
2799 }
2800
2801 // If no more variables, it either has no zeroes, or it vanished
2802 if (A_alg->type == COEFFICIENT_NUMERIC) {
2803 break;
2804 }
2805
2806 // The variable to eliminate
2807 lp_variable_t y = VAR(A_alg);
2808
2809 if (trace_is_enabled("coefficient")) {
2810 tracef("coefficient_resolve_algebraic(): resolving out %s\n", lp_variable_db_get_name(ctx->var_db, y));
2811 }
2812
2813 // Get the algebraic number of the value of y as a polynomial in y
2814 const lp_value_t* y_value = lp_assignment_get_value(m, y);
2815 if (y_value->type == LP_VALUE_NONE) {
2816 assert(coefficient_is_univariate(A_alg));
2817 break;
2818 }
2819
2820 if (lp_value_is_rational(y_value)) {
2821 lp_integer_t multiplier;
2822 integer_construct(&multiplier);
2823 coefficient_evaluate_rationals(ctx, A_alg, m, A_alg, &multiplier);
2824 integer_destruct(&multiplier);
2825 continue;
2826 }
2827
2828 // Proper algebraic number
2829 const lp_upolynomial_t* y_upoly = y_value->value.a.f;
2830 assert(y_upoly != 0);
2831
2832 // Make it multivariate
2833 coefficient_t y_poly;
2834 coefficient_construct_from_univariate(ctx, &y_poly, y_upoly, y);
2835
2836 if (trace_is_enabled("coefficient")) {
2837 tracef("coefficient_resolve_algebraic(): y_poly = "); coefficient_print(ctx, &y_poly, trace_out); tracef("\n");
2838 }
2839
2840 // Resolve
2841 coefficient_resultant(ctx, A_alg, A_alg, &y_poly);
2842
2843 if (trace_is_enabled("coefficient")) {
2844 tracef("coefficient_resolve_algebraic(): resultant done = "); coefficient_print(ctx, A_alg, trace_out); tracef("\n");
2845 }
2846
2847 // Destroy temp
2848 coefficient_destruct(&y_poly);
2849 }
2850 }
2851
2852
coefficient_roots_isolate(const lp_polynomial_context_t * ctx,const coefficient_t * A,const lp_assignment_t * M,lp_value_t * roots,size_t * roots_size)2853 void coefficient_roots_isolate(const lp_polynomial_context_t* ctx, const coefficient_t* A, const lp_assignment_t* M, lp_value_t* roots, size_t* roots_size) {
2854
2855 if (trace_is_enabled("coefficient::roots")) {
2856 tracef("coefficient_roots_isolate("); coefficient_print(ctx, A, trace_out); tracef(")\n");
2857 }
2858
2859 assert(A->type == COEFFICIENT_POLYNOMIAL);
2860 assert(lp_assignment_get_value(M, coefficient_top_variable(A))->type == LP_VALUE_NONE);
2861
2862 // The variable
2863 lp_variable_t x = coefficient_top_variable(A);
2864 assert(x != lp_variable_null);
2865
2866 // Evaluate in the rationals
2867 coefficient_t A_rat;
2868 lp_integer_t multiplier;
2869 integer_construct(&multiplier);
2870 coefficient_construct(ctx, &A_rat);
2871 coefficient_evaluate_rationals(ctx, A, M, &A_rat, &multiplier);
2872
2873 if (trace_is_enabled("coefficient::roots")) {
2874 tracef("coefficient_roots_isolate(): rational evaluation: "); coefficient_print(ctx, &A_rat, trace_out); tracef("\n");
2875 }
2876
2877 // If this is a constant polynomial, no zeroes
2878 if (coefficient_degree_safe(ctx, &A_rat, x) == 0) {
2879 if (trace_is_enabled("coefficient::roots")) {
2880 tracef("coefficient_roots_isolate(): constant, no roots\n");
2881 }
2882 *roots_size = 0;
2883 } else {
2884
2885 // TODO: normalize A_rat by gcd of coefficients
2886
2887 // If univariate, just isolate the univariate roots
2888 if (coefficient_is_univariate(&A_rat)) {
2889 coefficient_roots_isolate_univariate(ctx, &A_rat, roots, roots_size);
2890 assert(*roots_size <= coefficient_degree(&A_rat));
2891 } else {
2892
2893 // Can not evaluate in Q, do the expensive stuff
2894 if (trace_is_enabled("coefficient::roots")) {
2895 tracef("coefficient_roots_isolate(): not univariate\n");
2896 }
2897
2898 // Not univariate, A_rat(x, y0, ..., y_n), with x the top variable. We
2899 // eliminate the y variables with resolvents of the algebraic numbers to
2900 // obtain a univariate polynomial that includes the zeroes of A_rat.
2901 // We reverse the order, and just eliminate one-by one, until we reach x
2902
2903
2904 // The variable can vanish => no roots
2905 if (A_rat.type != COEFFICIENT_POLYNOMIAL || VAR(&A_rat) != VAR(A)) {
2906 *roots_size = 0;
2907 } else {
2908
2909 // Top variable must be unassigned
2910 assert(lp_assignment_get_value(M, x)->type == LP_VALUE_NONE);
2911
2912 // Rerder to (y_n, ..., y0, x) and then restore the order
2913 lp_variable_order_make_bot(ctx->var_order, x);
2914 coefficient_order(ctx, &A_rat);
2915
2916 // The copy where we compute the resultants
2917 coefficient_t A_alg;
2918 coefficient_construct(ctx, &A_alg);
2919
2920 // Resolve algebraic assignments
2921 coefficient_resolve_algebraic(ctx, &A_rat, M, &A_alg);
2922
2923 // Restore the order
2924 lp_variable_order_make_bot(ctx->var_order, lp_variable_null);
2925 coefficient_order(ctx, &A_rat);
2926 coefficient_order(ctx, &A_alg);
2927
2928 if (trace_is_enabled("coefficient::roots")) {
2929 tracef("coefficient_roots_isolate(): A_alg = "); coefficient_print(ctx, &A_alg, trace_out); tracef("\n");
2930 }
2931
2932 if (coefficient_is_zero(ctx, &A_alg)) {
2933
2934 if (trace_is_enabled("coefficient::roots")) {
2935 tracef("coefficient_roots_isolate(): vanished :(\n");
2936 }
2937
2938 //
2939 // A_alg vanished, but this might be due to other roots of the polynomials
2940 // that define the assignment. We will first manually check if the polynomial
2941 // A_rat vanishes by trying to evaluate coefficient.
2942 //
2943 // If A_rat vanishes, we're done, there is no roots.
2944 //
2945 // If A_rat doesn't vanish, we have that (modulo the assignment)
2946 //
2947 // A_rat = c_k x^k + ... + c_0
2948 //
2949 // If c_k evaluates to a_k, we can extend the assignment with a new
2950 // variable y -> a_k and get the roots of the polynomial
2951 //
2952 // B = y x^k + ...
2953 //
2954 // This polynomial will not vanish by construction since 0 is not a zero
2955 // of the polynomial defining a_k (by properties of root isolation,
2956 // defining polynomials are don't have 0 as a root).
2957 //
2958
2959 // Reduce based on the model
2960 coefficient_normalize_m(ctx, &A_rat, M);
2961
2962 // If A_rat vanishes, or is constant, there is no zeroes
2963 if (coefficient_is_constant(&A_rat)) {
2964 if (trace_is_enabled("coefficient::roots")) {
2965 tracef("coefficient_roots_isolate(): truly vanished :)\n");
2966 }
2967 *roots_size = 0;
2968 } else {
2969
2970 // Get the value of the coefficient
2971 lp_value_t* lc_value = coefficient_evaluate(ctx, coefficient_lc_safe(ctx, &A_rat, x), M);
2972
2973 if (trace_is_enabled("coefficient::roots")) {
2974 tracef("lc_value = "); lp_value_print(lc_value, trace_out); tracef("\n");
2975 }
2976
2977 // The new variable
2978 lp_variable_t y = lp_polynomial_context_get_temp_variable(ctx);
2979
2980 // Set the value
2981 assert(lp_assignment_get_value(M, y)->type == LP_VALUE_NONE);
2982 lp_assignment_set_value((lp_assignment_t*) M, y, lc_value);
2983 lp_variable_order_push((lp_variable_order_t*) ctx->var_order, y);
2984
2985 // Make B = y*x^k + ...
2986 // x is unassigned, y is assigned, we push y to the order
2987 // x will be bigger so it is in the right order
2988 coefficient_t y_coeff;
2989 coefficient_construct_simple_int(ctx, &y_coeff, 1, y, 1);
2990 size_t A_rat_deg = coefficient_degree_safe(ctx, &A_rat, x);
2991 assert(A_rat_deg > 0);
2992 coefficient_swap(COEFF(&A_rat, A_rat_deg), &y_coeff);
2993 coefficient_destruct(&y_coeff);
2994
2995 if (trace_is_enabled("coefficient::roots")) {
2996 tracef("A_rat (with y) = "); coefficient_print(ctx, &A_rat, trace_out); tracef("\n");
2997 }
2998
2999 // Now, do it recursively (restore the order first)
3000 coefficient_roots_isolate(ctx, &A_rat, M, roots, roots_size);
3001 assert(*roots_size <= A_rat_deg);
3002
3003 if (trace_is_enabled("coefficient::roots")) {
3004 tracef("A_rat (with y) = "); coefficient_print(ctx, &A_rat, trace_out); tracef("\n");
3005 tracef("roots_size = %zu\n", *roots_size);
3006 }
3007
3008 // Undo local stuff
3009 lp_assignment_set_value((lp_assignment_t*) M, y, 0);
3010 assert(y == lp_variable_order_top(ctx->var_order));
3011 lp_variable_order_pop((lp_variable_order_t*) ctx->var_order);
3012 lp_polynomial_context_release_temp_variable(ctx, y);
3013 lp_value_delete(lc_value);
3014 }
3015 } else if (coefficient_is_constant(&A_alg)) {
3016 // If constant, we have no roots
3017 if (trace_is_enabled("coefficient::roots")) {
3018 tracef("coefficient_roots_isolate(): no roots, constant\n");
3019 }
3020 *roots_size = 0;
3021 } else {
3022
3023 // Isolate as univariate, then filter any extra roots obtained from the defining polynomials
3024 if (trace_is_enabled("coefficient::roots")) {
3025 tracef("coefficient_roots_isolate(): univariate reduct\n");
3026 coefficient_print(ctx, &A_alg, trace_out);
3027 tracef("\n");
3028 }
3029 assert(coefficient_degree_safe(ctx, &A_alg, x) >= 1);
3030
3031 // Get the univariate version
3032 lp_upolynomial_t* A_alg_u = coefficient_to_univariate(ctx, &A_alg);
3033 size_t A_alg_u_deg = lp_upolynomial_degree(A_alg_u);
3034 assert(coefficient_degree_safe(ctx, &A_alg, x) == A_alg_u_deg);
3035 // Make space for the algebraic numbers
3036 lp_algebraic_number_t* algebraic_roots = malloc(A_alg_u_deg*sizeof(lp_algebraic_number_t));
3037 lp_upolynomial_roots_isolate(A_alg_u, algebraic_roots, roots_size);
3038 assert(*roots_size <= A_alg_u_deg);
3039
3040 if (trace_is_enabled("coefficient::roots")) {
3041 tracef("coefficient_roots_isolate(): filtering roots\n");
3042 }
3043
3044 // Filter any bad roots
3045 size_t i, to_keep;
3046 for (i = 0, to_keep = 0; i < *roots_size; ++i) {
3047
3048 // Set the value of the variable
3049 assert(lp_assignment_get_value(M, x)->type == LP_VALUE_NONE);
3050 lp_value_t x_value;
3051 lp_value_construct(&x_value, LP_VALUE_ALGEBRAIC, algebraic_roots + i);
3052 lp_assignment_set_value((lp_assignment_t*) M, x, &x_value);
3053
3054 if (trace_is_enabled("coefficient::roots")) {
3055 tracef("coefficient_roots_isolate(): checking root: ");
3056 lp_value_print(&x_value, trace_out);
3057 tracef("\n");
3058 }
3059
3060 // If zero by border or full sign is 0 then keep it
3061 if (coefficient_sgn(ctx, &A_rat, M) == 0) {
3062 if (i != to_keep) {
3063 lp_algebraic_number_swap(algebraic_roots + to_keep, algebraic_roots + i);
3064 }
3065 to_keep++;
3066 }
3067 // Remove the value
3068 lp_assignment_set_value((lp_assignment_t*) M, x, 0);
3069 lp_value_destruct(&x_value);
3070 }
3071 // Destruct the bad roots
3072 for (i = to_keep; i < *roots_size; ++ i) {
3073 lp_algebraic_number_destruct(algebraic_roots + i);
3074 }
3075 *roots_size = to_keep;
3076 assert(*roots_size <= A_alg_u_deg);
3077
3078 // Copy over the roots
3079 for (i = 0; i < *roots_size; ++i) {
3080 lp_value_construct(roots + i, LP_VALUE_ALGEBRAIC, algebraic_roots + i);
3081 lp_algebraic_number_destruct(algebraic_roots + i);
3082 }
3083 // Free the temps
3084 free(algebraic_roots);
3085 lp_upolynomial_delete(A_alg_u);
3086
3087 }
3088
3089 // Remove temps
3090 coefficient_destruct(&A_alg);
3091 }
3092 }
3093 }
3094
3095 // Remove temps
3096 coefficient_destruct(&A_rat);
3097 integer_destruct(&multiplier);
3098 }
3099
coefficient_evaluate(const lp_polynomial_context_t * ctx,const coefficient_t * C,const lp_assignment_t * M)3100 lp_value_t* coefficient_evaluate(const lp_polynomial_context_t* ctx, const coefficient_t* C, const lp_assignment_t* M) {
3101
3102 if (trace_is_enabled("coefficient")) {
3103 tracef("coefficient_evaluate("); coefficient_print(ctx, C, trace_out); tracef(")\n");
3104 }
3105
3106 lp_value_t* result = 0;
3107
3108 if (trace_is_enabled("coefficient")) {
3109 tracef("coefficient_evaluate(): approximating\n");
3110 }
3111
3112 // Compute the initial value of the approximation
3113 lp_rational_interval_t C_approx;
3114 lp_rational_interval_construct_zero(&C_approx);
3115 coefficient_value_approx(ctx, C, M, &C_approx);
3116
3117 if (trace_is_enabled("coefficient")) {
3118 tracef("coefficient_evaluate(): approximation = >"); lp_rational_interval_print(&C_approx, trace_out); tracef("\n");
3119 }
3120
3121 if (lp_rational_interval_is_point(&C_approx)) {
3122 // If the approximation is a point we're done
3123 result = lp_value_new(LP_VALUE_RATIONAL, lp_rational_interval_get_point(&C_approx));
3124 if (trace_is_enabled("coefficient")) {
3125 tracef("coefficient_evaluate(): approximation is rational = >"); lp_value_print(result, trace_out); tracef("\n");
3126 }
3127 } else {
3128
3129 if (trace_is_enabled("coefficient")) {
3130 tracef("coefficient_evaluate(): approximation is not rational\n");
3131 }
3132
3133 //
3134 // To evaluate the polynomial C(x1, ..., xn) with values v1, ..., vn
3135 // we construct the polynomial
3136 //
3137 // A(y, x1, ..., xn) = y - C(x1, ..., xn)
3138 //
3139 // and resolve out the defining polynomials p1(x1), ..., pn(xn) to get B(y)
3140 //
3141 // B(y) doesn't vanish since y has coefficient 1.
3142 //
3143 // Value of C is one of the roots of B(y). We can isolate the roots and
3144 // pick the one that intersects with the approximate value of C.
3145
3146 // Get the variable
3147 lp_variable_t y = lp_polynomial_context_get_temp_variable(ctx);
3148 lp_variable_order_make_bot(ctx->var_order, y);
3149
3150 coefficient_t A;
3151 lp_integer_t one;
3152 integer_construct_from_int(lp_Z, &one, 1);
3153 coefficient_construct_simple(ctx, &A, &one, y, 1);
3154 coefficient_sub(ctx, &A, &A, C);
3155 assert(VAR(&A) != y); // y should be the bottom variable
3156
3157 if (trace_is_enabled("coefficient")) {
3158 tracef("coefficient_evaluate(): resolving algebraic numbers\n");
3159 tracef("coefficient_evaluate(): start = "); coefficient_print(ctx, &A, trace_out), tracef("\n");
3160 }
3161
3162 // Resolve the algebraic numbers
3163 coefficient_t B;
3164 coefficient_construct(ctx, &B);
3165 coefficient_resolve_algebraic(ctx, &A, M, &B);
3166 assert(coefficient_is_univariate(&B));
3167 assert(VAR(&B) == y);
3168
3169 if (trace_is_enabled("coefficient")) {
3170 tracef("coefficient_evaluate(): resolved = "); coefficient_print(ctx, &B, trace_out), tracef("\n");
3171 tracef("coefficient_evaluate(): univariate root isolation\n");
3172 }
3173
3174 size_t A_resolved_deg = coefficient_degree(&B);
3175 lp_value_t* roots = malloc(sizeof(lp_value_t)*A_resolved_deg);
3176 size_t roots_size = 0;
3177
3178 // Find the roots of A_resolved
3179 coefficient_roots_isolate_univariate(ctx, &B, roots, &roots_size);
3180
3181 if (trace_is_enabled("coefficient")) {
3182 tracef("coefficient_evaluate(): roots:\n");
3183 size_t i = 0;
3184 for (i = 0; i < roots_size; ++ i) {
3185 tracef("%zu: ", i);
3186 lp_value_print(roots + i, trace_out);
3187 tracef("\n");
3188 }
3189 }
3190
3191 // Cache the variables
3192 lp_variable_list_t C_vars;
3193 lp_variable_list_construct(&C_vars);
3194 coefficient_get_variables(C, &C_vars);
3195
3196 // Cache the intervals
3197 lp_dyadic_interval_t* rational_interval_cache = algebraic_interval_remember(&C_vars, M);
3198
3199 // Filter the roots until we get only one
3200 for (;;) {
3201
3202 assert(!lp_rational_interval_is_point(&C_approx));
3203
3204 // Filter
3205 size_t i, to_keep;
3206 for (i = 0, to_keep = 0; i < roots_size; ++ i) {
3207 if (lp_rational_interval_contains_value(&C_approx, roots + i)) {
3208 lp_value_assign(roots + to_keep, roots + i);
3209 to_keep ++;
3210 }
3211 }
3212 for (i = to_keep; i < roots_size; ++ i) {
3213 lp_value_destruct(roots + i);
3214 }
3215 roots_size = to_keep;
3216 assert(roots_size > 0);
3217
3218 if (trace_is_enabled("coefficient")) {
3219 tracef("coefficient_evaluate(): filtered roots:\n");
3220 size_t i = 0;
3221 for (i = 0; i < roots_size; ++ i) {
3222 tracef("%zu: ", i);
3223 lp_value_print(roots + i, trace_out);
3224 tracef("\n");
3225 }
3226 }
3227
3228 // Did we find it
3229 if (roots_size == 1) {
3230 break;
3231 }
3232
3233 // Nope, refine
3234 for (i = 0; i < C_vars.list_size; ++ i) {
3235 lp_variable_t x_i = C_vars.list[i];
3236 const lp_value_t* x_i_value = lp_assignment_get_value(M, x_i);
3237 if (!lp_value_is_rational(x_i_value)) {
3238 lp_algebraic_number_refine_const(&x_i_value->value.a);
3239 }
3240 }
3241
3242 // Approximate again
3243 coefficient_value_approx(ctx, C, M, &C_approx);
3244
3245 // If we reduced to a point, we're done
3246 if (lp_rational_interval_is_point(&C_approx)) {
3247 break;
3248 }
3249 }
3250
3251 if (lp_rational_interval_is_point(&C_approx)) {
3252 // Rational value approximation reduced to
3253 result = lp_value_new(LP_VALUE_RATIONAL, lp_rational_interval_get_point(&C_approx));
3254 } else {
3255 // We have the value in roots[0]
3256 result = lp_value_new_copy(roots);
3257 }
3258
3259 // Remove the roots
3260 assert(roots_size == 1);
3261 lp_value_destruct(roots);
3262 free(roots);
3263
3264 // Release the variable
3265 lp_polynomial_context_release_temp_variable(ctx, y);
3266 lp_variable_order_make_bot(ctx->var_order, lp_variable_null);
3267
3268 // Release the cache
3269 algebraic_interval_restore(&C_vars, rational_interval_cache, M);
3270
3271 // Remove temps
3272 integer_destruct(&one);
3273 coefficient_destruct(&A);
3274 coefficient_destruct(&B);
3275 lp_variable_list_destruct(&C_vars);
3276
3277 }
3278
3279 lp_rational_interval_destruct(&C_approx);
3280
3281 return result;
3282 }
3283
3284 static
hash_pair(size_t a,size_t b)3285 size_t hash_pair(size_t a, size_t b) {
3286 return a + 0x9e3779b9 + (b << 6) + (b >> 2);
3287 }
3288
coefficient_hash_traverse(const lp_polynomial_context_t * ctx,lp_monomial_t * p,void * hash_void)3289 void coefficient_hash_traverse(const lp_polynomial_context_t* ctx, lp_monomial_t* p, void* hash_void) {
3290 (void)(ctx);
3291 size_t* hash = (size_t*)(hash_void);
3292 *hash ^= integer_hash(&p->a);
3293 size_t i;
3294 for (i = 0; i < p->n; ++ i) {
3295 *hash ^= hash_pair(p->p[i].x, p->p[i].d);
3296 }
3297 }
3298
coefficient_hash(const lp_polynomial_context_t * ctx,const coefficient_t * A)3299 size_t coefficient_hash(const lp_polynomial_context_t* ctx, const coefficient_t* A) {
3300 size_t hash = 0;
3301
3302 // Traverse the polynomial and hash each of the monomials
3303 lp_monomial_t m_tmp;
3304 lp_monomial_construct(ctx, &m_tmp);
3305 coefficient_traverse(ctx, A, coefficient_hash_traverse, &m_tmp, &hash);
3306 lp_monomial_destruct(&m_tmp);
3307
3308 return hash;
3309 }
3310