1 /*
2 * This file is part of the Yices SMT Solver.
3 * Copyright (C) 2017 SRI International.
4 *
5 * Yices is free software: you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation, either version 3 of the License, or
8 * (at your option) any later version.
9 *
10 * Yices is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with Yices. If not, see <http://www.gnu.org/licenses/>.
17 */
18
19 /*
20 * Products of variables: x_1^d_1 ... x_n^d_n
21 *
22 * These are used to represent polynomials
23 * Variables are just integer indices
24 */
25
26 #include <assert.h>
27 #include <stdbool.h>
28
29 #include "terms/power_products.h"
30 #include "utils/hash_functions.h"
31 #include "utils/memalloc.h"
32 #include "utils/prng.h"
33
34 #include "yices_limits.h"
35
36
37 /*
38 * Initialization and deletion of buffers
39 */
init_pp_buffer(pp_buffer_t * b,uint32_t n)40 void init_pp_buffer(pp_buffer_t *b, uint32_t n) {
41 if (n >= PPROD_MAX_LENGTH) {
42 out_of_memory();
43 }
44 b->size = n;
45 b->len = 0;
46 b->prod = (varexp_t *) safe_malloc(n * sizeof(varexp_t));
47 }
48
delete_pp_buffer(pp_buffer_t * b)49 void delete_pp_buffer(pp_buffer_t *b) {
50 safe_free(b->prod);
51 b->prod = NULL;
52 }
53
54
55 /*
56 * Auxiliary functions
57 */
58 // make buffer 50% larger
pp_buffer_extend(pp_buffer_t * b)59 static void pp_buffer_extend(pp_buffer_t *b) {
60 uint32_t n;
61
62 n = b->size + 1;
63 n += n >> 1;
64 if (n >= PPROD_MAX_LENGTH) {
65 out_of_memory();
66 }
67 b->prod = (varexp_t *) safe_realloc(b->prod, n * sizeof(varexp_t));
68 b->size = n;
69 }
70
71 // make it large enough for at least n elements
pp_buffer_resize(pp_buffer_t * b,uint32_t n)72 static void pp_buffer_resize(pp_buffer_t *b, uint32_t n) {
73 uint32_t new_size;
74
75 if (b->size < n) {
76 new_size = b->size + 1;
77 new_size += new_size >> 1;
78 if (new_size < n) new_size = n;
79
80 if (new_size >= PPROD_MAX_LENGTH) {
81 out_of_memory();
82 }
83 b->prod = (varexp_t *) safe_realloc(b->prod, new_size * sizeof(varexp_t));
84 b->size = new_size;
85 }
86 }
87
88 // add pair v^d to the buffer
pp_buffer_pushback(pp_buffer_t * b,int32_t v,uint32_t d)89 static void pp_buffer_pushback(pp_buffer_t *b, int32_t v, uint32_t d) {
90 uint32_t i;
91
92 i = b->len;
93 if (i == b->size) {
94 pp_buffer_extend(b);
95 }
96 assert(i < b->size);
97 b->prod[i].var = v;
98 b->prod[i].exp = d;
99 b->len = i + 1;
100 }
101
102 // add pairs v[0]^d[0] ... v[n-1]^d[n-1]
pp_buffer_pusharray(pp_buffer_t * b,uint32_t n,int32_t * v,uint32_t * d)103 static void pp_buffer_pusharray(pp_buffer_t *b, uint32_t n, int32_t *v, uint32_t *d) {
104 uint32_t i, l;
105
106 l = b->len;
107 pp_buffer_resize(b, l + n);
108 for (i=0; i<n; i++) {
109 b->prod[l+i].var = v[i];
110 b->prod[l+i].exp = d[i];
111 }
112 b->len = l + n;
113 }
114
115 // add pairs v[0]^1 ... v[n-1]^1
pp_buffer_pushvars(pp_buffer_t * b,uint32_t n,int32_t * v)116 static void pp_buffer_pushvars(pp_buffer_t *b, uint32_t n, int32_t *v) {
117 uint32_t i, l;
118
119 l = b->len;
120 pp_buffer_resize(b, l + n);
121 for (i=0; i<n; i++) {
122 b->prod[l+i].var = v[i];
123 b->prod[l+i].exp = 1;
124 }
125 b->len = l + n;
126 }
127
128 // variant: pairs are in varexp array a
pp_buffer_push_varexp_array(pp_buffer_t * b,uint32_t n,varexp_t * a)129 static void pp_buffer_push_varexp_array(pp_buffer_t *b, uint32_t n, varexp_t *a) {
130 uint32_t i, l;
131
132 l = b->len;
133 pp_buffer_resize(b, l + n);
134 for (i=0; i<n; i++) {
135 b->prod[l + i] = a[i];
136 }
137 b->len = l + n;
138 }
139
140
141
142 /*
143 * NORMALIZATION
144 */
145
146 /*
147 * Sort in increasing order of variables:
148 * - a = array of pairs <variable, exponent>
149 * - n = size of the array.
150 */
151 static void isort_varexp_array(varexp_t *a, uint32_t n);
152 static void qsort_varexp_array(varexp_t *a, uint32_t n);
153
sort_varexp_array(varexp_t * a,uint32_t n)154 static inline void sort_varexp_array(varexp_t *a, uint32_t n) {
155 if (n <= 10) {
156 isort_varexp_array(a, n);
157 } else {
158 qsort_varexp_array(a, n);
159 }
160 }
161
isort_varexp_array(varexp_t * a,uint32_t n)162 static void isort_varexp_array(varexp_t *a, uint32_t n) {
163 uint32_t i, j, d, e;
164 int32_t v, w;
165
166 for (i=1; i<n; i++) {
167 v = a[i].var;
168 d = a[i].exp;
169 j = 0;
170 while (a[j].var < v) j ++;
171 while (j < i) {
172 w = a[j].var; a[j].var = v; v = w;
173 e = a[j].exp; a[j].exp = d; d = e;
174 j ++;
175 }
176 a[j].var = v;
177 a[j].exp = d;
178 }
179 }
180
qsort_varexp_array(varexp_t * a,uint32_t n)181 static void qsort_varexp_array(varexp_t *a, uint32_t n) {
182 uint32_t i, j;
183 int32_t pivot;
184 varexp_t aux;
185 uint32_t seed = PRNG_DEFAULT_SEED;
186
187 // random pivot
188 i = random_uint(&seed, n);
189 aux = a[i];
190 pivot = a[i].var;
191
192 // swap with a[0];
193 a[i] = a[0];
194 a[0] = aux;
195
196 i = 0;
197 j = n;
198
199 // The test i <= j in the second loop is required for termination
200 // if all elements are smaller than the pivot.
201 do { j--; } while (a[j].var > pivot);
202 do { i++; } while (i <= j && a[i].var < pivot);
203
204 while (i < j) {
205 aux = a[i]; a[i] = a[j]; a[j] = aux;
206
207 do { j--; } while (a[j].var > pivot);
208 do { i++; } while (a[i].var < pivot);
209 }
210
211 // swap pivot = a[0] and a[j]
212 aux = a[0]; a[0] = a[j]; a[j] = aux;
213
214 // sort a[0 ... j-1] and a[j+1 ... n-1]
215 sort_varexp_array(a, j);
216 j ++;
217 sort_varexp_array(a + j, n - j);
218 }
219
220
221
222 /*
223 * Sort then merge pairs <var, d> with the same var.
224 * Also remove pairs <var, d> with d = 0.
225 * return the number of elements remaining after normalization
226 */
normalize_varexp_array(varexp_t * a,uint32_t n)227 static uint32_t normalize_varexp_array(varexp_t *a, uint32_t n) {
228 uint32_t i, j, d;
229 int32_t v;
230
231 if (n == 0) return n;
232
233 sort_varexp_array(a, n);
234
235 j = 0;
236 d = a[0].exp;
237 v = a[0].var;
238 for (i=1; i<n; i++) {
239 if (a[i].var == v) {
240 d += a[i].exp;
241 } else {
242 if (d != 0) {
243 a[j].var = v;
244 a[j].exp = d;
245 j ++;
246 }
247 v = a[i].var;
248 d = a[i].exp;
249 }
250 }
251
252 if (d != 0) {
253 a[j].var = v;
254 a[j].exp = d;
255 j ++;
256 }
257
258 return j;
259 }
260
261
pp_buffer_normalize(pp_buffer_t * b)262 void pp_buffer_normalize(pp_buffer_t *b) {
263 b->len = normalize_varexp_array(b->prod, b->len);
264 }
265
266
267 /*
268 * API: exported operations
269 */
pp_buffer_reset(pp_buffer_t * b)270 void pp_buffer_reset(pp_buffer_t *b) {
271 b->len = 0;
272 }
273
pp_buffer_set_var(pp_buffer_t * b,int32_t v)274 void pp_buffer_set_var(pp_buffer_t *b, int32_t v) {
275 b->len = 0;
276 pp_buffer_pushback(b, v, 1);
277 }
278
pp_buffer_set_varexp(pp_buffer_t * b,int32_t v,uint32_t d)279 void pp_buffer_set_varexp(pp_buffer_t *b, int32_t v, uint32_t d) {
280 b->len = 0;
281 pp_buffer_pushback(b, v, d);
282 }
283
pp_buffer_set_vars(pp_buffer_t * b,uint32_t n,int32_t * v)284 void pp_buffer_set_vars(pp_buffer_t *b, uint32_t n, int32_t *v) {
285 b->len = 0;
286 pp_buffer_pushvars(b, n, v);
287 }
288
pp_buffer_set_varexps(pp_buffer_t * b,uint32_t n,int32_t * v,uint32_t * d)289 void pp_buffer_set_varexps(pp_buffer_t *b, uint32_t n, int32_t *v, uint32_t *d) {
290 b->len = 0;
291 pp_buffer_pusharray(b, n, v, d);
292 }
293
pp_buffer_set_pprod(pp_buffer_t * b,pprod_t * p)294 void pp_buffer_set_pprod(pp_buffer_t *b, pprod_t *p) {
295 assert(p != end_pp);
296
297 b->len = 0;
298 if (pp_is_var(p)) {
299 pp_buffer_pushback(b, var_of_pp(p), 1);
300 } else if (!pp_is_empty(p)) {
301 pp_buffer_push_varexp_array(b, p->len, p->prod);
302 }
303 }
304
pp_buffer_copy(pp_buffer_t * b,pp_buffer_t * src)305 void pp_buffer_copy(pp_buffer_t *b, pp_buffer_t *src) {
306 b->len = 0;
307 pp_buffer_push_varexp_array(b, src->len, src->prod);
308 }
309
310
pp_buffer_mul_var(pp_buffer_t * b,int32_t v)311 void pp_buffer_mul_var(pp_buffer_t *b, int32_t v) {
312 pp_buffer_pushback(b, v, 1);
313 pp_buffer_normalize(b);
314 }
315
pp_buffer_mul_varexp(pp_buffer_t * b,int32_t v,uint32_t d)316 void pp_buffer_mul_varexp(pp_buffer_t *b, int32_t v, uint32_t d) {
317 pp_buffer_pushback(b, v, d);
318 pp_buffer_normalize(b);
319 }
320
pp_buffer_mul_vars(pp_buffer_t * b,uint32_t n,int32_t * v)321 void pp_buffer_mul_vars(pp_buffer_t *b, uint32_t n, int32_t *v) {
322 pp_buffer_pushvars(b, n, v);
323 pp_buffer_normalize(b);
324 }
325
pp_buffer_mul_varexps(pp_buffer_t * b,uint32_t n,int32_t * v,uint32_t * d)326 void pp_buffer_mul_varexps(pp_buffer_t *b, uint32_t n, int32_t *v, uint32_t *d) {
327 pp_buffer_pusharray(b, n, v, d);
328 pp_buffer_normalize(b);
329 }
330
pp_buffer_mul_pprod(pp_buffer_t * b,pprod_t * p)331 void pp_buffer_mul_pprod(pp_buffer_t *b, pprod_t *p) {
332 assert(p != end_pp);
333
334 if (pp_is_var(p)) {
335 pp_buffer_pushback(b, var_of_pp(p), 1);
336 } else if (! pp_is_empty(p)) {
337 pp_buffer_push_varexp_array(b, p->len, p->prod);
338 }
339 pp_buffer_normalize(b);
340 }
341
pp_buffer_mul_buffer(pp_buffer_t * b,pp_buffer_t * src)342 void pp_buffer_mul_buffer(pp_buffer_t *b, pp_buffer_t *src) {
343 pp_buffer_push_varexp_array(b, src->len, src->prod);
344 pp_buffer_normalize(b);
345 }
346
347
pp_buffer_push_var(pp_buffer_t * b,int32_t v)348 void pp_buffer_push_var(pp_buffer_t *b, int32_t v) {
349 pp_buffer_pushback(b, v, 1);
350 }
351
pp_buffer_push_varexp(pp_buffer_t * b,int32_t v,uint32_t d)352 void pp_buffer_push_varexp(pp_buffer_t *b, int32_t v, uint32_t d) {
353 pp_buffer_pushback(b, v, d);
354 }
355
356
357 /*
358 * Exponentiation
359 */
pp_buffer_exponentiate(pp_buffer_t * b,uint32_t d)360 void pp_buffer_exponentiate(pp_buffer_t *b, uint32_t d) {
361 uint32_t i, n;
362
363 n = b->len;
364 for (i=0; i<n; i++) {
365 b->prod[i].exp *= d;
366 }
367 pp_buffer_normalize(b);
368 }
369
370
371 /*
372 * Degree
373 */
varexp_array_degree(varexp_t * a,int32_t n)374 static uint32_t varexp_array_degree(varexp_t *a, int32_t n) {
375 uint32_t d, i;
376
377 d = 0;
378 for (i=0; i<n; i++) {
379 d += a[i].exp;
380 }
381 return d;
382 }
383
pp_buffer_degree(pp_buffer_t * b)384 uint32_t pp_buffer_degree(pp_buffer_t *b) {
385 return varexp_array_degree(b->prod, b->len);
386 }
387
pprod_degree(pprod_t * p)388 uint32_t pprod_degree(pprod_t *p) {
389 uint32_t d;
390
391 assert(p != end_pp);
392
393 d = 0;
394 if (pp_is_var(p)) {
395 d = 1;
396 } else if (! pp_is_empty(p)) {
397 d = p->degree;
398 }
399 return d;
400 }
401
402
403 /*
404 * Check that degree is less than YICES_MAX_DEGREE
405 */
below_max_degree(varexp_t * a,uint32_t n)406 static bool below_max_degree(varexp_t *a, uint32_t n) {
407 uint32_t d, i, e;
408
409 d = 0;
410 for (i=0; i<n; i++) {
411 e = a[i].exp;
412 if (e >= YICES_MAX_DEGREE) return false;
413 d += e;
414 if (d >= YICES_MAX_DEGREE) return false;
415 }
416
417 return true;
418 }
419
pp_buffer_below_max_degree(pp_buffer_t * b)420 bool pp_buffer_below_max_degree(pp_buffer_t *b) {
421 return below_max_degree(b->prod, b->len);
422 }
423
pprod_below_max_degree(pprod_t * p)424 bool pprod_below_max_degree(pprod_t *p) {
425 return pprod_degree(p) < YICES_MAX_DEGREE;
426 }
427
428
429
430 /*
431 * Degree of a variable x in a product p
432 * - we use a sequential search since n is usually small
433 */
varexp_array_var_degree(varexp_t * a,uint32_t n,int32_t x)434 static uint32_t varexp_array_var_degree(varexp_t *a, uint32_t n, int32_t x) {
435 uint32_t d, i;
436
437 d = 0;
438 for (i=0; i<n; i++) {
439 if (a[i].var == x) {
440 d = a[i].exp;
441 break;
442 }
443 }
444 return d;
445 }
446
pp_buffer_var_degree(pp_buffer_t * b,int32_t x)447 uint32_t pp_buffer_var_degree(pp_buffer_t *b, int32_t x) {
448 return varexp_array_var_degree(b->prod, b->len, x);
449 }
450
pprod_var_degree(pprod_t * p,int32_t x)451 uint32_t pprod_var_degree(pprod_t *p, int32_t x) {
452 uint32_t d;
453
454 assert(p != end_pp);
455
456 d = 0;
457 if (pp_is_var(p)) {
458 d = (var_of_pp(p) == x);
459 } else if (! pp_is_empty(p)) {
460 d =varexp_array_var_degree(p->prod, p->len, x);
461 }
462 return d;
463 }
464
465
466
467 /*
468 * Build a pprod object by making a copy of a
469 * - a must be normalized
470 * - n = length of a
471 * - n must be less than PPROD_MAX_LENGTH
472 */
make_pprod(varexp_t * a,uint32_t n)473 pprod_t *make_pprod(varexp_t *a, uint32_t n) {
474 pprod_t *p;
475 uint32_t i;
476
477 assert(0 <= n && n < PPROD_MAX_LENGTH);
478
479 if (n == 0) {
480 p = empty_pp;
481 } else if (n == 1 && a[0].exp == 1) {
482 p = var_pp(a[0].var);
483 } else {
484 p = (pprod_t *) safe_malloc(sizeof(pprod_t) + n * sizeof(varexp_t));
485 p->len = n;
486 p->degree = varexp_array_degree(a, n);
487 for (i=0; i<n; i++) {
488 p->prod[i] = a[i];
489 }
490 }
491
492 return p;
493 }
494
495
496 /*
497 * Construct a pprod object from buffer b
498 */
pp_buffer_getprod(pp_buffer_t * b)499 pprod_t *pp_buffer_getprod(pp_buffer_t *b) {
500 return make_pprod(b->prod, b->len);
501 }
502
503
504
505 /*
506 * Free the object constructed by make_pprod
507 */
free_pprod(pprod_t * p)508 void free_pprod(pprod_t *p) {
509 if (! pp_is_var(p) && ! pp_is_empty(p)) {
510 safe_free(p);
511 }
512 }
513
514
515 /*
516 * Comparison: two arrays of the same length n
517 * Both arrays must be normalized.
518 */
varexp_array_equal(varexp_t * a,varexp_t * b,uint32_t n)519 bool varexp_array_equal(varexp_t *a, varexp_t *b, uint32_t n) {
520 uint32_t i;
521
522 for (i=0; i<n; i++) {
523 if (a[i].var != b[i].var || a[i].exp != b[i].exp) {
524 return false;
525 }
526 }
527 return true;
528 }
529
530
531 /*
532 * Check whether b1 and b2 are equal
533 * - both must be normalized
534 */
pp_buffer_equal(pp_buffer_t * b1,pp_buffer_t * b2)535 bool pp_buffer_equal(pp_buffer_t *b1, pp_buffer_t *b2) {
536 return b1->len == b2->len && varexp_array_equal(b1->prod, b2->prod, b1->len);
537 }
538
539
540 /*
541 * Minimum of two exponent
542 */
min_exp(uint32_t e1,uint32_t e2)543 static inline uint32_t min_exp(uint32_t e1, uint32_t e2) {
544 return (e1 < e2) ? e1 : e2;
545 }
546
547 /*
548 * Compute the common factors of b and b1. Store the result in b.
549 * - b and b1 must be normalized
550 * - the result is normalized (on exit)
551 */
pp_buffer_gcd(pp_buffer_t * b,pp_buffer_t * b1)552 void pp_buffer_gcd(pp_buffer_t *b, pp_buffer_t *b1) {
553 uint32_t i, j, n, d;
554 int32_t x;
555
556 n = b->len;
557 j = 0;
558 for (i=0; i<n; i++) {
559 assert(b->prod[i].exp > 0);
560 x = b->prod[i].var;
561 d = pp_buffer_var_degree(b1, x);
562 if (d == 0) continue; // remove x
563 b->prod[j].var = x;
564 b->prod[j].exp = min_exp(d, b->prod[i].exp);
565 j ++;
566 }
567 b->len = j;
568 }
569
570
571 /*
572 * Divide b by b1. Store the result in b
573 * - b and b1 must be normalized.
574 *
575 * The actual operation is more general:
576 * We replace x^e in b by x^max(0, e - d) where d = exponent of x in b1.
577 * So x disappears if e <= d, and x has exponent (e - d) in the result otherwise.
578 */
pp_buffer_divide(pp_buffer_t * b,pp_buffer_t * b1)579 void pp_buffer_divide(pp_buffer_t *b, pp_buffer_t *b1) {
580 uint32_t i, j, n, d;
581 int32_t x;
582
583 n = b->len;
584 j = 0;
585 for (i=0; i<n; i++) {
586 assert(b->prod[i].exp > 0);
587 x = b->prod[i].var;
588 d = pp_buffer_var_degree(b1, x);
589 if (d >= b->prod[i].exp) continue;
590 b->prod[j].var = x;
591 b->prod[j].exp = b->prod[i].exp - d;
592 j ++;
593 }
594 b->len = j;
595 }
596
597
598 /*
599 * Divide p by variable x
600 */
pp_buffer_divide_by_var(pp_buffer_t * b,int32_t x)601 void pp_buffer_divide_by_var(pp_buffer_t *b, int32_t x) {
602 uint32_t i, n;
603
604 n = b->len;
605 for (i=0; i<n; i++) {
606 assert(b->prod[i].exp > 0);
607 if (b->prod[i].var == x) goto found;
608 }
609
610 return; // x does not occur in b.
611
612 found:
613 assert(i < n && b->prod[i].var == x && b->prod[i].exp > 0);
614 b->prod[i].exp --;
615 if (b->prod[i].exp == 0) {
616 i ++;
617 while (i < n) {
618 b->prod[i-1] = b->prod[i];
619 i ++;
620 }
621 b->len = n-1;
622 }
623 }
624
625
626 /*
627 * The ordering between power products must be compatible with the
628 * product and with the ordering on variables. That is, we want
629 * 1) a < b => a * c < b * c for any c
630 * 2) x < y as variables => x^1 y^0 < x^0 y^1
631 *
632 * The lexical ordering defined as follows works:
633 * Let a = x_1^d_1 .... x_n^d_n and b = x_1^c_1 ... x_n^c_n
634 * then a < b if for some i we have
635 * d_1 = c_1 and ... and d_{i-1} = c_{i-1} and d_i > c_i
636 *
637 * Input:
638 * - a and b must be normalized, na = length of a, nb = length of b
639 * Output:
640 * - a negative number if a < b
641 * - zero if a == b
642 * - a positive number if a > b
643 */
varexp_array_lexcmp(varexp_t * a,uint32_t na,varexp_t * b,uint32_t nb)644 static int32_t varexp_array_lexcmp(varexp_t *a, uint32_t na, varexp_t *b, uint32_t nb) {
645 uint32_t i, n;
646 int32_t x;
647
648 n = (na < nb) ? na : nb;
649
650 i = 0;
651 while (i<n && a[i].var == b[i].var && a[i].exp == b[i].exp) {
652 i++;
653 }
654
655 if (i < n) {
656 /*
657 * let x_b = b[i].var and x_a = a[i].var
658 * d_b = b[i].exp and d_a = a[i].exp
659 *
660 * cases:
661 * x_a == x_b and d_a < d_b --> b smaller than a
662 * x_a == x_b and d_a > d_b --> a smaller than b
663 * x_a < x_b --> a smaller than b
664 * x_a > x_b --> b smaller than a
665 */
666 x = a[i].var - b[i].var;
667 if (x == 0) {
668 // the conversion to int32_t is safe if the total
669 // degree of a and b is less than YICES_MAX_DEGREE
670 assert(((int32_t) b[i].exp >= 0) && ((int32_t) a[i].exp >= 0));
671 return (int32_t) (b[i].exp - a[i].exp);
672 } else {
673 return x;
674 }
675 } else {
676 /*
677 * if na < nb, a is a prefix of b so a < b
678 * if nb < na, b is a prefix of a so b < a
679 */
680 return na - nb;
681 }
682 }
683
684
685 /*
686 * Lexicographic comparison: p1 and p2 must be normalized
687 */
pprod_lex_cmp(pprod_t * p1,pprod_t * p2)688 int32_t pprod_lex_cmp(pprod_t *p1, pprod_t *p2) {
689 varexp_t aux1, aux2;
690 varexp_t *a, *b;
691 uint32_t na, nb;
692
693 assert(p1 != end_pp && p2 != end_pp);
694
695 na = 0;
696 a = NULL;
697 if (pp_is_var(p1)) {
698 aux1.var = var_of_pp(p1);
699 aux1.exp = 1;
700 a = &aux1;
701 na = 1;
702 } else if (! pp_is_empty(p1)) {
703 a = p1->prod;
704 na = p1->len;
705 }
706
707 nb = 0;
708 b = NULL;
709 if (pp_is_var(p2)) {
710 aux2.var = var_of_pp(p2);
711 aux2.exp = 1;
712 b = &aux2;
713 nb = 1;
714 } else if (! pp_is_empty(p2)) {
715 b = p2->prod;
716 nb = p2->len;
717 }
718
719 return varexp_array_lexcmp(a, na, b, nb);
720 }
721
722
723
724 /*
725 * Check whether p1 precedes p2 (strictly) in the degree + lex ordering
726 */
pprod_precedes(pprod_t * p1,pprod_t * p2)727 bool pprod_precedes(pprod_t *p1, pprod_t *p2) {
728 uint32_t d1, d2;
729
730 if (p1 == p2) return false;
731 if (p1 == end_pp) return false;
732 if (p2 == end_pp) return true;
733
734 d1 = pprod_degree(p1);
735 d2 = pprod_degree(p2);
736 return (d1 < d2) || (d1 == d2 && pprod_lex_cmp(p1, p2) < 0);
737 }
738
739
740 /*
741 * Test whether a divides b
742 */
varexp_array_divides(varexp_t * a,uint32_t na,varexp_t * b,uint32_t nb)743 static bool varexp_array_divides(varexp_t *a, uint32_t na, varexp_t *b, uint32_t nb) {
744 uint32_t i, j;
745 int32_t v;
746
747 if (na > nb) return false;
748
749 j = 0;
750 for (i=0; i<na; i++) {
751 v = a[i].var;
752 while (j < nb && b[j].var < v) {
753 j++;
754 }
755 if (j == nb || b[j].var > v || a[i].exp > b[j].exp) {
756 return false;
757 }
758 }
759
760 return true;
761 }
762
763
764 /*
765 * Test whether p1 divides p2
766 */
pprod_divides(pprod_t * p1,pprod_t * p2)767 bool pprod_divides(pprod_t *p1, pprod_t *p2) {
768 varexp_t aux1, aux2;
769 varexp_t *a, *b;
770 uint32_t na, nb;
771
772 assert(p1 != end_pp && p2 != end_pp);
773
774 na = 0;
775 a = NULL;
776 if (pp_is_var(p1)) {
777 aux1.var = var_of_pp(p1);
778 aux1.exp = 1;
779 a = &aux1;
780 na = 1;
781 } else if (! pp_is_empty(p1)) {
782 a = p1->prod;
783 na = p1->len;
784 }
785
786 nb = 0;
787 b = NULL;
788 if (pp_is_var(p2)) {
789 aux2.var = var_of_pp(p2);
790 aux2.exp = 1;
791 b = &aux2;
792 nb = 1;
793 } else if (! pp_is_empty(p2)) {
794 b = p2->prod;
795 nb = p2->len;
796 }
797
798 return varexp_array_divides(a, na, b, nb);
799 }
800
801
802
803
804 /*
805 * Check whether p1 divides p2 and if so store (p2/p1) in b
806 */
pprod_divisor(pp_buffer_t * b,pprod_t * p1,pprod_t * p2)807 bool pprod_divisor(pp_buffer_t *b, pprod_t *p1, pprod_t *p2) {
808 uint32_t i, j, n;
809 int32_t v;
810 varexp_t *a1, *a2;
811 varexp_t aux;
812
813 assert(p1 != end_pp && p2 != end_pp);
814
815 if (pp_is_empty(p1)) {
816 n = 0;
817 a1 = NULL;
818 } else if (pp_is_var(p1)) {
819 aux.var = var_of_pp(p1);
820 aux.exp = 1;
821 n = 1;
822 a1 = &aux;
823 } else {
824 n = p1->len;
825 a1 = p1->prod;
826 }
827
828 pp_buffer_set_pprod(b, p2);
829 if (n > b->len) return false;
830
831 pp_buffer_pushback(b, INT32_MAX, UINT32_MAX); // add an end marker
832
833 a2 = b->prod;
834
835 j =0;
836 for (i=0; i<n; i++) {
837 v = a1[i].var;
838 while (a2[j].var < v) j++;
839 if (a2[j].var > v || a2[j].exp < a1[i].exp) return false;
840 a2[j].exp -= a1[i].exp;
841 j ++;
842 }
843
844 // normalization: remove zero exponents and end marker
845 b->len --;
846 i = 0;
847 for (j=0; j<b->len; j++) {
848 if (a2[j].exp > 0) {
849 a2[i] = a2[j];
850 i ++;
851 }
852 }
853
854 b->len = i;
855
856 return true;
857 }
858
859
860 /*
861 * Check whether p1 and p2 are equal
862 */
pprod_equal(pprod_t * p1,pprod_t * p2)863 bool pprod_equal(pprod_t *p1, pprod_t *p2) {
864 assert(p1 != end_pp && p2 != end_pp);
865
866 if (p1 == p2) return true;
867
868 if (pp_is_var(p1) | pp_is_var(p2)) {
869 // p1 and p2 are distinct variables
870 // or only one of them is a variable
871 return false;
872 }
873
874 if (pp_is_empty(p1) | pp_is_empty(p2)) {
875 // one empty and the other is not
876 return false;
877 }
878
879 // two non-empty, non variable power products
880 return p1->len == p2->len && varexp_array_equal(p1->prod, p2->prod, p1->len);
881 }
882
883