?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
isl_qpolynomial_add: replace 0-degree polynomials by their constant terms
[isl.git] / isl_polynomial.c
blobe0f5adbe1d3bb1aacd5f0c6e113edd7c50ea61d2
1 #include <isl_seq.h>
2 #include <isl_polynomial_private.h>
3 #include <isl_point_private.h>
4 #include <isl_dim_private.h>
5 #include <isl_map_private.h>
6
7 int isl_upoly_is_cst(__isl_keep struct isl_upoly *up)
8 {
9 if (!up)
10 return -1;
11
12 return up->var < 0;
13 }
14
15 __isl_keep struct isl_upoly_cst *isl_upoly_as_cst(__isl_keep struct isl_upoly *up)
16 {
17 if (!up)
18 return NULL;
19
20 isl_assert(up->ctx, up->var < 0, return NULL);
21
22 return (struct isl_upoly_cst *)up;
23 }
24
25 __isl_keep struct isl_upoly_rec *isl_upoly_as_rec(__isl_keep struct isl_upoly *up)
26 {
27 if (!up)
28 return NULL;
29
30 isl_assert(up->ctx, up->var >= 0, return NULL);
31
32 return (struct isl_upoly_rec *)up;
33 }
34
35 int isl_upoly_is_zero(__isl_keep struct isl_upoly *up)
36 {
37 struct isl_upoly_cst *cst;
38
39 if (!up)
40 return -1;
41 if (!isl_upoly_is_cst(up))
42 return 0;
43
44 cst = isl_upoly_as_cst(up);
45 if (!cst)
46 return -1;
47
48 return isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d);
49 }
50
51 int isl_upoly_is_nan(__isl_keep struct isl_upoly *up)
52 {
53 struct isl_upoly_cst *cst;
54
55 if (!up)
56 return -1;
57 if (!isl_upoly_is_cst(up))
58 return 0;
59
60 cst = isl_upoly_as_cst(up);
61 if (!cst)
62 return -1;
63
64 return isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d);
65 }
66
67 int isl_upoly_is_infty(__isl_keep struct isl_upoly *up)
68 {
69 struct isl_upoly_cst *cst;
70
71 if (!up)
72 return -1;
73 if (!isl_upoly_is_cst(up))
74 return 0;
75
76 cst = isl_upoly_as_cst(up);
77 if (!cst)
78 return -1;
79
80 return isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d);
81 }
82
83 int isl_upoly_is_neginfty(__isl_keep struct isl_upoly *up)
84 {
85 struct isl_upoly_cst *cst;
86
87 if (!up)
88 return -1;
89 if (!isl_upoly_is_cst(up))
90 return 0;
91
92 cst = isl_upoly_as_cst(up);
93 if (!cst)
94 return -1;
95
96 return isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d);
97 }
98
99 int isl_upoly_is_one(__isl_keep struct isl_upoly *up)
100 {
101 struct isl_upoly_cst *cst;
102
103 if (!up)
104 return -1;
105 if (!isl_upoly_is_cst(up))
106 return 0;
107
108 cst = isl_upoly_as_cst(up);
109 if (!cst)
110 return -1;
111
112 return isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
113 }
114
115 int isl_upoly_is_negone(__isl_keep struct isl_upoly *up)
116 {
117 struct isl_upoly_cst *cst;
118
119 if (!up)
120 return -1;
121 if (!isl_upoly_is_cst(up))
122 return 0;
123
124 cst = isl_upoly_as_cst(up);
125 if (!cst)
126 return -1;
127
128 return isl_int_is_negone(cst->n) && isl_int_is_one(cst->d);
129 }
130
131 __isl_give struct isl_upoly_cst *isl_upoly_cst_alloc(struct isl_ctx *ctx)
132 {
133 struct isl_upoly_cst *cst;
134
135 cst = isl_alloc_type(ctx, struct isl_upoly_cst);
136 if (!cst)
137 return NULL;
138
139 cst->up.ref = 1;
140 cst->up.ctx = ctx;
141 isl_ctx_ref(ctx);
142 cst->up.var = -1;
143
144 isl_int_init(cst->n);
145 isl_int_init(cst->d);
146
147 return cst;
148 }
149
150 __isl_give struct isl_upoly *isl_upoly_zero(struct isl_ctx *ctx)
151 {
152 struct isl_upoly_cst *cst;
153
154 cst = isl_upoly_cst_alloc(ctx);
155 if (!cst)
156 return NULL;
157
158 isl_int_set_si(cst->n, 0);
159 isl_int_set_si(cst->d, 1);
160
161 return &cst->up;
162 }
163
164 __isl_give struct isl_upoly *isl_upoly_infty(struct isl_ctx *ctx)
165 {
166 struct isl_upoly_cst *cst;
167
168 cst = isl_upoly_cst_alloc(ctx);
169 if (!cst)
170 return NULL;
171
172 isl_int_set_si(cst->n, 1);
173 isl_int_set_si(cst->d, 0);
174
175 return &cst->up;
176 }
177
178 __isl_give struct isl_upoly *isl_upoly_nan(struct isl_ctx *ctx)
179 {
180 struct isl_upoly_cst *cst;
181
182 cst = isl_upoly_cst_alloc(ctx);
183 if (!cst)
184 return NULL;
185
186 isl_int_set_si(cst->n, 0);
187 isl_int_set_si(cst->d, 0);
188
189 return &cst->up;
190 }
191
192 __isl_give struct isl_upoly *isl_upoly_rat_cst(struct isl_ctx *ctx,
193 isl_int n, isl_int d)
194 {
195 struct isl_upoly_cst *cst;
196
197 cst = isl_upoly_cst_alloc(ctx);
198 if (!cst)
199 return NULL;
200
201 isl_int_set(cst->n, n);
202 isl_int_set(cst->d, d);
203
204 return &cst->up;
205 }
206
207 __isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
208 int var, int size)
209 {
210 struct isl_upoly_rec *rec;
211
212 isl_assert(ctx, var >= 0, return NULL);
213 isl_assert(ctx, size >= 0, return NULL);
214 rec = isl_calloc(dim->ctx, struct isl_upoly_rec,
215 sizeof(struct isl_upoly_rec) +
216 (size - 1) * sizeof(struct isl_upoly *));
217 if (!rec)
218 return NULL;
219
220 rec->up.ref = 1;
221 rec->up.ctx = ctx;
222 isl_ctx_ref(ctx);
223 rec->up.var = var;
224
225 rec->n = 0;
226 rec->size = size;
227
228 return rec;
229 error:
230 isl_upoly_free(&rec->up);
231 return NULL;
232 }
233
234 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
235 {
236 struct isl_upoly_cst *cst;
237
238 if (!qp)
239 return -1;
240
241 return isl_upoly_is_zero(qp->upoly);
242 }
243
244 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
245 {
246 struct isl_upoly_cst *cst;
247
248 if (!qp)
249 return -1;
250
251 return isl_upoly_is_one(qp->upoly);
252 }
253
254 static void upoly_free_cst(__isl_take struct isl_upoly_cst *cst)
255 {
256 isl_int_clear(cst->n);
257 isl_int_clear(cst->d);
258 }
259
260 static void upoly_free_rec(__isl_take struct isl_upoly_rec *rec)
261 {
262 int i;
263
264 for (i = 0; i < rec->n; ++i)
265 isl_upoly_free(rec->p[i]);
266 }
267
268 __isl_give struct isl_upoly *isl_upoly_copy(__isl_keep struct isl_upoly *up)
269 {
270 if (!up)
271 return NULL;
272
273 up->ref++;
274 return up;
275 }
276
277 __isl_give struct isl_upoly *isl_upoly_dup_cst(__isl_keep struct isl_upoly *up)
278 {
279 struct isl_upoly_cst *cst;
280 struct isl_upoly_cst *dup;
281
282 cst = isl_upoly_as_cst(up);
283 if (!cst)
284 return NULL;
285
286 dup = isl_upoly_as_cst(isl_upoly_zero(up->ctx));
287 if (!dup)
288 return NULL;
289 isl_int_set(dup->n, cst->n);
290 isl_int_set(dup->d, cst->d);
291
292 return &dup->up;
293 }
294
295 __isl_give struct isl_upoly *isl_upoly_dup_rec(__isl_keep struct isl_upoly *up)
296 {
297 int i;
298 struct isl_upoly_rec *rec;
299 struct isl_upoly_rec *dup;
300
301 rec = isl_upoly_as_rec(up);
302 if (!rec)
303 return NULL;
304
305 dup = isl_upoly_alloc_rec(up->ctx, up->var, rec->n);
306 if (!dup)
307 return NULL;
308
309 for (i = 0; i < rec->n; ++i) {
310 dup->p[i] = isl_upoly_copy(rec->p[i]);
311 if (!dup->p[i])
312 goto error;
313 dup->n++;
314 }
315
316 return &dup->up;
317 error:
318 isl_upoly_free(&dup->up);
319 return NULL;
320 }
321
322 __isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
323 {
324 struct isl_upoly *dup;
325
326 if (!up)
327 return NULL;
328
329 if (isl_upoly_is_cst(up))
330 return isl_upoly_dup_cst(up);
331 else
332 return isl_upoly_dup_rec(up);
333 }
334
335 __isl_give struct isl_upoly *isl_upoly_cow(__isl_take struct isl_upoly *up)
336 {
337 if (!up)
338 return NULL;
339
340 if (up->ref == 1)
341 return up;
342 up->ref--;
343 return isl_upoly_dup(up);
344 }
345
346 void isl_upoly_free(__isl_take struct isl_upoly *up)
347 {
348 if (!up)
349 return;
350
351 if (--up->ref > 0)
352 return;
353
354 if (up->var < 0)
355 upoly_free_cst((struct isl_upoly_cst *)up);
356 else
357 upoly_free_rec((struct isl_upoly_rec *)up);
358
359 isl_ctx_deref(up->ctx);
360 free(up);
361 }
362
363 static void isl_upoly_cst_reduce(__isl_keep struct isl_upoly_cst *cst)
364 {
365 isl_int gcd;
366
367 isl_int_init(gcd);
368 isl_int_gcd(gcd, cst->n, cst->d);
369 if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
370 isl_int_divexact(cst->n, cst->n, gcd);
371 isl_int_divexact(cst->d, cst->d, gcd);
372 }
373 isl_int_clear(gcd);
374 }
375
376 __isl_give struct isl_upoly *isl_upoly_sum_cst(__isl_take struct isl_upoly *up1,
377 __isl_take struct isl_upoly *up2)
378 {
379 struct isl_upoly_cst *cst1;
380 struct isl_upoly_cst *cst2;
381
382 up1 = isl_upoly_cow(up1);
383 if (!up1 || !up2)
384 goto error;
385
386 cst1 = isl_upoly_as_cst(up1);
387 cst2 = isl_upoly_as_cst(up2);
388
389 if (isl_int_eq(cst1->d, cst2->d))
390 isl_int_add(cst1->n, cst1->n, cst2->n);
391 else {
392 isl_int_mul(cst1->n, cst1->n, cst2->d);
393 isl_int_addmul(cst1->n, cst2->n, cst1->d);
394 isl_int_mul(cst1->d, cst1->d, cst2->d);
395 }
396
397 isl_upoly_cst_reduce(cst1);
398
399 isl_upoly_free(up2);
400 return up1;
401 error:
402 isl_upoly_free(up1);
403 isl_upoly_free(up2);
404 return NULL;
405 }
406
407 static __isl_give struct isl_upoly *replace_by_zero(
408 __isl_take struct isl_upoly *up)
409 {
410 struct isl_ctx *ctx;
411
412 if (!up)
413 return NULL;
414 ctx = up->ctx;
415 isl_upoly_free(up);
416 return isl_upoly_zero(ctx);
417 }
418
419 static __isl_give struct isl_upoly *replace_by_constant_term(
420 __isl_take struct isl_upoly *up)
421 {
422 struct isl_upoly_rec *rec;
423 struct isl_upoly *cst;
424
425 if (!up)
426 return NULL;
427
428 rec = isl_upoly_as_rec(up);
429 if (!rec)
430 goto error;
431 isl_assert(up->ctx, rec->n == 1, goto error);
432 cst = isl_upoly_copy(rec->p[0]);
433 isl_upoly_free(up);
434 return cst;
435 error:
436 isl_upoly_free(up);
437 return NULL;
438 }
439
440 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
441 __isl_take struct isl_upoly *up2)
442 {
443 int i;
444 struct isl_upoly_rec *rec1, *rec2;
445
446 if (!up1 || !up2)
447 goto error;
448
449 if (isl_upoly_is_nan(up1)) {
450 isl_upoly_free(up2);
451 return up1;
452 }
453
454 if (isl_upoly_is_nan(up2)) {
455 isl_upoly_free(up1);
456 return up2;
457 }
458
459 if (isl_upoly_is_zero(up1)) {
460 isl_upoly_free(up1);
461 return up2;
462 }
463
464 if (isl_upoly_is_zero(up2)) {
465 isl_upoly_free(up2);
466 return up1;
467 }
468
469 if (up1->var < up2->var)
470 return isl_upoly_sum(up2, up1);
471
472 if (up2->var < up1->var) {
473 struct isl_upoly_rec *rec;
474 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
475 isl_upoly_free(up1);
476 return up2;
477 }
478 up1 = isl_upoly_cow(up1);
479 rec = isl_upoly_as_rec(up1);
480 if (!rec)
481 goto error;
482 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
483 if (rec->n == 1)
484 up1 = replace_by_constant_term(up1);
485 return up1;
486 }
487
488 if (isl_upoly_is_cst(up1))
489 return isl_upoly_sum_cst(up1, up2);
490
491 rec1 = isl_upoly_as_rec(up1);
492 rec2 = isl_upoly_as_rec(up2);
493 if (!rec1 || !rec2)
494 goto error;
495
496 if (rec1->n < rec2->n)
497 return isl_upoly_sum(up2, up1);
498
499 up1 = isl_upoly_cow(up1);
500 rec1 = isl_upoly_as_rec(up1);
501 if (!rec1)
502 goto error;
503
504 for (i = rec2->n - 1; i >= 0; --i) {
505 rec1->p[i] = isl_upoly_sum(rec1->p[i],
506 isl_upoly_copy(rec2->p[i]));
507 if (!rec1->p[i])
508 goto error;
509 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
510 isl_upoly_free(rec1->p[i]);
511 rec1->n--;
512 }
513 }
514
515 if (rec1->n == 0)
516 up1 = replace_by_zero(up1);
517 else if (rec1->n == 1)
518 up1 = replace_by_constant_term(up1);
519
520 isl_upoly_free(up2);
521
522 return up1;
523 error:
524 isl_upoly_free(up1);
525 isl_upoly_free(up2);
526 return NULL;
527 }
528
529 __isl_give struct isl_upoly *isl_upoly_neg_cst(__isl_take struct isl_upoly *up)
530 {
531 struct isl_upoly_cst *cst;
532
533 if (isl_upoly_is_zero(up))
534 return up;
535
536 up = isl_upoly_cow(up);
537 if (!up)
538 return NULL;
539
540 cst = isl_upoly_as_cst(up);
541
542 isl_int_neg(cst->n, cst->n);
543
544 return up;
545 }
546
547 __isl_give struct isl_upoly *isl_upoly_neg(__isl_take struct isl_upoly *up)
548 {
549 int i;
550 struct isl_upoly_rec *rec;
551
552 if (!up)
553 return NULL;
554
555 if (isl_upoly_is_cst(up))
556 return isl_upoly_neg_cst(up);
557
558 up = isl_upoly_cow(up);
559 rec = isl_upoly_as_rec(up);
560 if (!rec)
561 goto error;
562
563 for (i = 0; i < rec->n; ++i) {
564 rec->p[i] = isl_upoly_neg(rec->p[i]);
565 if (!rec->p[i])
566 goto error;
567 }
568
569 return up;
570 error:
571 isl_upoly_free(up);
572 return NULL;
573 }
574
575 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
576 __isl_take struct isl_upoly *up2)
577 {
578 struct isl_upoly_cst *cst1;
579 struct isl_upoly_cst *cst2;
580
581 up1 = isl_upoly_cow(up1);
582 if (!up1 || !up2)
583 goto error;
584
585 cst1 = isl_upoly_as_cst(up1);
586 cst2 = isl_upoly_as_cst(up2);
587
588 isl_int_mul(cst1->n, cst1->n, cst2->n);
589 isl_int_mul(cst1->d, cst1->d, cst2->d);
590
591 isl_upoly_cst_reduce(cst1);
592
593 isl_upoly_free(up2);
594 return up1;
595 error:
596 isl_upoly_free(up1);
597 isl_upoly_free(up2);
598 return NULL;
599 }
600
601 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
602 __isl_take struct isl_upoly *up2)
603 {
604 struct isl_upoly_rec *rec1;
605 struct isl_upoly_rec *rec2;
606 struct isl_upoly_rec *res;
607 int i, j;
608 int size;
609
610 rec1 = isl_upoly_as_rec(up1);
611 rec2 = isl_upoly_as_rec(up2);
612 if (!rec1 || !rec2)
613 goto error;
614 size = rec1->n + rec2->n - 1;
615 res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
616 if (!res)
617 goto error;
618
619 for (i = 0; i < rec1->n; ++i) {
620 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
621 isl_upoly_copy(rec1->p[i]));
622 if (!res->p[i])
623 goto error;
624 res->n++;
625 }
626 for (; i < size; ++i) {
627 res->p[i] = isl_upoly_zero(up1->ctx);
628 if (!res->p[i])
629 goto error;
630 res->n++;
631 }
632 for (i = 0; i < rec1->n; ++i) {
633 for (j = 1; j < rec2->n; ++j) {
634 struct isl_upoly *up;
635 up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
636 isl_upoly_copy(rec1->p[i]));
637 res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
638 if (!res->p[i + j])
639 goto error;
640 }
641 }
642
643 isl_upoly_free(up1);
644 isl_upoly_free(up2);
645
646 return &res->up;
647 error:
648 isl_upoly_free(up1);
649 isl_upoly_free(up2);
650 isl_upoly_free(&res->up);
651 return NULL;
652 }
653
654 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
655 __isl_take struct isl_upoly *up2)
656 {
657 if (!up1 || !up2)
658 goto error;
659
660 if (isl_upoly_is_nan(up1)) {
661 isl_upoly_free(up2);
662 return up1;
663 }
664
665 if (isl_upoly_is_nan(up2)) {
666 isl_upoly_free(up1);
667 return up2;
668 }
669
670 if (isl_upoly_is_zero(up1)) {
671 isl_upoly_free(up2);
672 return up1;
673 }
674
675 if (isl_upoly_is_zero(up2)) {
676 isl_upoly_free(up1);
677 return up2;
678 }
679
680 if (isl_upoly_is_one(up1)) {
681 isl_upoly_free(up1);
682 return up2;
683 }
684
685 if (isl_upoly_is_one(up2)) {
686 isl_upoly_free(up2);
687 return up1;
688 }
689
690 if (up1->var < up2->var)
691 return isl_upoly_mul(up2, up1);
692
693 if (up2->var < up1->var) {
694 int i;
695 struct isl_upoly_rec *rec;
696 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
697 isl_ctx *ctx = up1->ctx;
698 isl_upoly_free(up1);
699 isl_upoly_free(up2);
700 return isl_upoly_nan(ctx);
701 }
702 up1 = isl_upoly_cow(up1);
703 rec = isl_upoly_as_rec(up1);
704 if (!rec)
705 goto error;
706
707 for (i = 0; i < rec->n; ++i) {
708 rec->p[i] = isl_upoly_mul(rec->p[i],
709 isl_upoly_copy(up2));
710 if (!rec->p[i])
711 goto error;
712 }
713 isl_upoly_free(up2);
714 return up1;
715 }
716
717 if (isl_upoly_is_cst(up1))
718 return isl_upoly_mul_cst(up1, up2);
719
720 return isl_upoly_mul_rec(up1, up2);
721 error:
722 isl_upoly_free(up1);
723 isl_upoly_free(up2);
724 return NULL;
725 }
726
727 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_dim *dim,
728 unsigned n_div)
729 {
730 struct isl_qpolynomial *qp;
731 unsigned total;
732
733 if (!dim)
734 return NULL;
735
736 total = isl_dim_total(dim);
737
738 qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
739 if (!qp)
740 return NULL;
741
742 qp->ref = 1;
743 qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
744 if (!qp->div)
745 goto error;
746
747 qp->dim = dim;
748
749 return qp;
750 error:
751 isl_dim_free(dim);
752 isl_qpolynomial_free(qp);
753 return NULL;
754 }
755
756 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
757 {
758 if (!qp)
759 return NULL;
760
761 qp->ref++;
762 return qp;
763 }
764
765 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
766 {
767 struct isl_qpolynomial *dup;
768
769 if (!qp)
770 return NULL;
771
772 dup = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), qp->div->n_row);
773 if (!dup)
774 return NULL;
775 isl_mat_free(dup->div);
776 dup->div = isl_mat_copy(qp->div);
777 if (!dup->div)
778 goto error;
779 dup->upoly = isl_upoly_copy(qp->upoly);
780 if (!dup->upoly)
781 goto error;
782
783 return dup;
784 error:
785 isl_qpolynomial_free(dup);
786 return NULL;
787 }
788
789 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
790 {
791 if (!qp)
792 return NULL;
793
794 if (qp->ref == 1)
795 return qp;
796 qp->ref--;
797 return isl_qpolynomial_dup(qp);
798 }
799
800 void isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
801 {
802 if (!qp)
803 return;
804
805 if (--qp->ref > 0)
806 return;
807
808 isl_dim_free(qp->dim);
809 isl_mat_free(qp->div);
810 isl_upoly_free(qp->upoly);
811
812 free(qp);
813 }
814
815 static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
816 {
817 int n_row, n_col;
818 int equal;
819
820 isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
821 div1->n_col >= div2->n_col, return -1);
822
823 if (div1->n_row == div2->n_row)
824 return isl_mat_is_equal(div1, div2);
825
826 n_row = div1->n_row;
827 n_col = div1->n_col;
828 div1->n_row = div2->n_row;
829 div1->n_col = div2->n_col;
830
831 equal = isl_mat_is_equal(div1, div2);
832
833 div1->n_row = n_row;
834 div1->n_col = n_col;
835
836 return equal;
837 }
838
839 static void expand_row(__isl_keep isl_mat *dst, int d,
840 __isl_keep isl_mat *src, int s, int *exp)
841 {
842 int i;
843 unsigned c = src->n_col - src->n_row;
844
845 isl_seq_cpy(dst->row[d], src->row[s], c);
846 isl_seq_clr(dst->row[d] + c, dst->n_col - c);
847
848 for (i = 0; i < s; ++i)
849 isl_int_set(dst->row[d][c + exp[i]], src->row[s][c + i]);
850 }
851
852 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
853 {
854 int li, lj;
855
856 li = isl_seq_last_non_zero(div->row[i], div->n_col);
857 lj = isl_seq_last_non_zero(div->row[j], div->n_col);
858
859 if (li != lj)
860 return li - lj;
861
862 return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
863 }
864
865 static __isl_give isl_mat *merge_divs(__isl_keep isl_mat *div1,
866 __isl_keep isl_mat *div2, int *exp1, int *exp2)
867 {
868 int i, j, k;
869 isl_mat *div = NULL;
870 unsigned d = div1->n_col - div1->n_row;
871
872 div = isl_mat_alloc(div1->ctx, 1 + div1->n_row + div2->n_row,
873 d + div1->n_row + div2->n_row);
874 if (!div)
875 return NULL;
876
877 for (i = 0, j = 0, k = 0; i < div1->n_row && j < div2->n_row; ++k) {
878 int cmp;
879
880 expand_row(div, k, div1, i, exp1);
881 expand_row(div, k + 1, div2, j, exp2);
882
883 cmp = cmp_row(div, k, k + 1);
884 if (cmp == 0) {
885 exp1[i++] = k;
886 exp2[j++] = k;
887 } else if (cmp < 0) {
888 exp1[i++] = k;
889 } else {
890 exp2[j++] = k;
891 isl_seq_cpy(div->row[k], div->row[k + 1], div->n_col);
892 }
893 }
894 for (; i < div1->n_row; ++i, ++k) {
895 expand_row(div, k, div1, i, exp1);
896 exp1[i] = k;
897 }
898 for (; j < div2->n_row; ++j, ++k) {
899 expand_row(div, k, div2, j, exp2);
900 exp2[j] = k;
901 }
902
903 div->n_row = k;
904 div->n_col = d + k;
905
906 return div;
907 }
908
909 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
910 int *exp, int first)
911 {
912 int i;
913 struct isl_upoly_rec *rec;
914
915 if (isl_upoly_is_cst(up))
916 return up;
917
918 if (up->var < first)
919 return up;
920
921 if (exp[up->var - first] == up->var - first)
922 return up;
923
924 up = isl_upoly_cow(up);
925 if (!up)
926 goto error;
927
928 up->var = exp[up->var - first] + first;
929
930 rec = isl_upoly_as_rec(up);
931 if (!rec)
932 goto error;
933
934 for (i = 0; i < rec->n; ++i) {
935 rec->p[i] = expand(rec->p[i], exp, first);
936 if (!rec->p[i])
937 goto error;
938 }
939
940 return up;
941 error:
942 isl_upoly_free(up);
943 return NULL;
944 }
945
946 static __isl_give isl_qpolynomial *with_merged_divs(
947 __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
948 __isl_take isl_qpolynomial *qp2),
949 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
950 {
951 int *exp1 = NULL;
952 int *exp2 = NULL;
953 isl_mat *div = NULL;
954
955 qp1 = isl_qpolynomial_cow(qp1);
956 qp2 = isl_qpolynomial_cow(qp2);
957
958 if (!qp1 || !qp2)
959 goto error;
960
961 isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
962 qp1->div->n_col >= qp2->div->n_col, goto error);
963
964 exp1 = isl_alloc_array(qp1->div->ctx, int, qp1->div->n_row);
965 exp2 = isl_alloc_array(qp2->div->ctx, int, qp2->div->n_row);
966 if (!exp1 || !exp2)
967 goto error;
968
969 div = merge_divs(qp1->div, qp2->div, exp1, exp2);
970 if (!div)
971 goto error;
972
973 isl_mat_free(qp1->div);
974 qp1->div = isl_mat_copy(div);
975 isl_mat_free(qp2->div);
976 qp2->div = isl_mat_copy(div);
977
978 qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
979 qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
980
981 if (!qp1->upoly || !qp2->upoly)
982 goto error;
983
984 isl_mat_free(div);
985 free(exp1);
986 free(exp2);
987
988 return fn(qp1, qp2);
989 error:
990 isl_mat_free(div);
991 free(exp1);
992 free(exp2);
993 isl_qpolynomial_free(qp1);
994 isl_qpolynomial_free(qp2);
995 return NULL;
996 }
997
998 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
999 __isl_take isl_qpolynomial *qp2)
1000 {
1001 qp1 = isl_qpolynomial_cow(qp1);
1002
1003 if (!qp1 || !qp2)
1004 goto error;
1005
1006 if (qp1->div->n_row < qp2->div->n_row)
1007 return isl_qpolynomial_add(qp2, qp1);
1008
1009 isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1010 if (!compatible_divs(qp1->div, qp2->div))
1011 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1012
1013 qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1014 if (!qp1->upoly)
1015 goto error;
1016
1017 isl_qpolynomial_free(qp2);
1018
1019 return qp1;
1020 error:
1021 isl_qpolynomial_free(qp1);
1022 isl_qpolynomial_free(qp2);
1023 return NULL;
1024 }
1025
1026 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1027 __isl_take isl_qpolynomial *qp2)
1028 {
1029 return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1030 }
1031
1032 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1033 {
1034 qp = isl_qpolynomial_cow(qp);
1035
1036 if (!qp)
1037 return NULL;
1038
1039 qp->upoly = isl_upoly_neg(qp->upoly);
1040 if (!qp->upoly)
1041 goto error;
1042
1043 return qp;
1044 error:
1045 isl_qpolynomial_free(qp);
1046 return NULL;
1047 }
1048
1049 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1050 __isl_take isl_qpolynomial *qp2)
1051 {
1052 qp1 = isl_qpolynomial_cow(qp1);
1053
1054 if (!qp1 || !qp2)
1055 goto error;
1056
1057 if (qp1->div->n_row < qp2->div->n_row)
1058 return isl_qpolynomial_mul(qp2, qp1);
1059
1060 isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1061 if (!compatible_divs(qp1->div, qp2->div))
1062 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1063
1064 qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1065 if (!qp1->upoly)
1066 goto error;
1067
1068 isl_qpolynomial_free(qp2);
1069
1070 return qp1;
1071 error:
1072 isl_qpolynomial_free(qp1);
1073 isl_qpolynomial_free(qp2);
1074 return NULL;
1075 }
1076
1077 __isl_give isl_qpolynomial *isl_qpolynomial_zero(__isl_take isl_dim *dim)
1078 {
1079 struct isl_qpolynomial *qp;
1080 struct isl_upoly_cst *cst;
1081
1082 qp = isl_qpolynomial_alloc(dim, 0);
1083 if (!qp)
1084 return NULL;
1085
1086 qp->upoly = isl_upoly_zero(dim->ctx);
1087 if (!qp->upoly)
1088 goto error;
1089
1090 return qp;
1091 error:
1092 isl_qpolynomial_free(qp);
1093 return NULL;
1094 }
1095
1096 __isl_give isl_qpolynomial *isl_qpolynomial_infty(__isl_take isl_dim *dim)
1097 {
1098 struct isl_qpolynomial *qp;
1099 struct isl_upoly_cst *cst;
1100
1101 qp = isl_qpolynomial_alloc(dim, 0);
1102 if (!qp)
1103 return NULL;
1104
1105 qp->upoly = isl_upoly_infty(dim->ctx);
1106 if (!qp->upoly)
1107 goto error;
1108
1109 return qp;
1110 error:
1111 isl_qpolynomial_free(qp);
1112 return NULL;
1113 }
1114
1115 __isl_give isl_qpolynomial *isl_qpolynomial_nan(__isl_take isl_dim *dim)
1116 {
1117 struct isl_qpolynomial *qp;
1118 struct isl_upoly_cst *cst;
1119
1120 qp = isl_qpolynomial_alloc(dim, 0);
1121 if (!qp)
1122 return NULL;
1123
1124 qp->upoly = isl_upoly_nan(dim->ctx);
1125 if (!qp->upoly)
1126 goto error;
1127
1128 return qp;
1129 error:
1130 isl_qpolynomial_free(qp);
1131 return NULL;
1132 }
1133
1134 __isl_give isl_qpolynomial *isl_qpolynomial_cst(__isl_take isl_dim *dim,
1135 isl_int v)
1136 {
1137 struct isl_qpolynomial *qp;
1138 struct isl_upoly_cst *cst;
1139
1140 qp = isl_qpolynomial_alloc(dim, 0);
1141 if (!qp)
1142 return NULL;
1143
1144 qp->upoly = isl_upoly_zero(dim->ctx);
1145 if (!qp->upoly)
1146 goto error;
1147 cst = isl_upoly_as_cst(qp->upoly);
1148 isl_int_set(cst->n, v);
1149
1150 return qp;
1151 error:
1152 isl_qpolynomial_free(qp);
1153 return NULL;
1154 }
1155
1156 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1157 isl_int *n, isl_int *d)
1158 {
1159 struct isl_upoly_cst *cst;
1160
1161 if (!qp)
1162 return -1;
1163
1164 if (!isl_upoly_is_cst(qp->upoly))
1165 return 0;
1166
1167 cst = isl_upoly_as_cst(qp->upoly);
1168 if (!cst)
1169 return -1;
1170
1171 if (n)
1172 isl_int_set(*n, cst->n);
1173 if (d)
1174 isl_int_set(*d, cst->d);
1175
1176 return 1;
1177 }
1178
1179 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_dim *dim,
1180 int pos, int power)
1181 {
1182 int i;
1183 struct isl_qpolynomial *qp;
1184 struct isl_upoly_rec *rec;
1185 struct isl_upoly_cst *cst;
1186 struct isl_ctx *ctx;
1187
1188 if (!dim)
1189 return NULL;
1190
1191 ctx = dim->ctx;
1192
1193 qp = isl_qpolynomial_alloc(dim, 0);
1194 if (!qp)
1195 return NULL;
1196
1197 qp->upoly = &isl_upoly_alloc_rec(ctx, pos, 1 + power)->up;
1198 if (!qp->upoly)
1199 goto error;
1200 rec = isl_upoly_as_rec(qp->upoly);
1201 for (i = 0; i < 1 + power; ++i) {
1202 rec->p[i] = isl_upoly_zero(ctx);
1203 if (!rec->p[i])
1204 goto error;
1205 rec->n++;
1206 }
1207 cst = isl_upoly_as_cst(rec->p[power]);
1208 isl_int_set_si(cst->n, 1);
1209
1210 return qp;
1211 error:
1212 isl_qpolynomial_free(qp);
1213 return NULL;
1214 }
1215
1216 __isl_give isl_qpolynomial *isl_qpolynomial_var(__isl_take isl_dim *dim,
1217 enum isl_dim_type type, unsigned pos)
1218 {
1219 if (!dim)
1220 return NULL;
1221
1222 isl_assert(dim->ctx, isl_dim_size(dim, isl_dim_in) == 0, goto error);
1223 isl_assert(dim->ctx, pos < isl_dim_size(dim, type), goto error);
1224
1225 if (type == isl_dim_set)
1226 pos += isl_dim_size(dim, isl_dim_param);
1227
1228 return isl_qpolynomial_pow(dim, pos, 1);
1229 error:
1230 isl_dim_free(dim);
1231 return NULL;
1232 }
1233
1234 __isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
1235 int power)
1236 {
1237 struct isl_qpolynomial *qp = NULL;
1238 struct isl_upoly_rec *rec;
1239 struct isl_upoly_cst *cst;
1240 int i;
1241 int pos;
1242
1243 if (!div)
1244 return NULL;
1245 isl_assert(div->ctx, div->bmap->n_div == 1, goto error);
1246
1247 qp = isl_qpolynomial_alloc(isl_basic_map_get_dim(div->bmap), 1);
1248 if (!qp)
1249 goto error;
1250
1251 isl_seq_cpy(qp->div->row[0], div->line[0], qp->div->n_col - 1);
1252 isl_int_set_si(qp->div->row[0][qp->div->n_col - 1], 0);
1253
1254 pos = isl_dim_total(qp->dim);
1255 qp->upoly = &isl_upoly_alloc_rec(div->ctx, pos, 1 + power)->up;
1256 if (!qp->upoly)
1257 goto error;
1258 rec = isl_upoly_as_rec(qp->upoly);
1259 for (i = 0; i < 1 + power; ++i) {
1260 rec->p[i] = isl_upoly_zero(div->ctx);
1261 if (!rec->p[i])
1262 goto error;
1263 rec->n++;
1264 }
1265 cst = isl_upoly_as_cst(rec->p[power]);
1266 isl_int_set_si(cst->n, 1);
1267
1268 isl_div_free(div);
1269
1270 return qp;
1271 error:
1272 isl_qpolynomial_free(qp);
1273 isl_div_free(div);
1274 return NULL;
1275 }
1276
1277 __isl_give isl_qpolynomial *isl_qpolynomial_div(__isl_take isl_div *div)
1278 {
1279 return isl_qpolynomial_div_pow(div, 1);
1280 }
1281
1282 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst(__isl_take isl_dim *dim,
1283 const isl_int n, const isl_int d)
1284 {
1285 struct isl_qpolynomial *qp;
1286 struct isl_upoly_cst *cst;
1287
1288 qp = isl_qpolynomial_alloc(dim, 0);
1289 if (!qp)
1290 return NULL;
1291
1292 qp->upoly = isl_upoly_zero(dim->ctx);
1293 if (!qp->upoly)
1294 goto error;
1295 cst = isl_upoly_as_cst(qp->upoly);
1296 isl_int_set(cst->n, n);
1297 isl_int_set(cst->d, d);
1298
1299 return qp;
1300 error:
1301 isl_qpolynomial_free(qp);
1302 return NULL;
1303 }
1304
1305 #undef PW
1306 #define PW isl_pw_qpolynomial
1307 #undef EL
1308 #define EL isl_qpolynomial
1309 #undef IS_ZERO
1310 #define IS_ZERO is_zero
1311 #undef FIELD
1312 #define FIELD qp
1313 #undef ADD
1314 #define ADD add
1315
1316 #include <isl_pw_templ.c>
1317
1318 #undef PW
1319 #define PW isl_pw_qpolynomial_fold
1320 #undef EL
1321 #define EL isl_qpolynomial_fold
1322 #undef IS_ZERO
1323 #define IS_ZERO is_empty
1324 #undef FIELD
1325 #define FIELD fold
1326 #undef ADD
1327 #define ADD fold
1328
1329 #include <isl_pw_templ.c>
1330
1331 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
1332 {
1333 if (!pwqp)
1334 return -1;
1335
1336 if (pwqp->n != -1)
1337 return 0;
1338
1339 if (!isl_set_fast_is_universe(pwqp->p[0].set))
1340 return 0;
1341
1342 return isl_qpolynomial_is_one(pwqp->p[0].qp);
1343 }
1344
1345 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
1346 __isl_take isl_pw_qpolynomial *pwqp1,
1347 __isl_take isl_pw_qpolynomial *pwqp2)
1348 {
1349 int i, j, n;
1350 struct isl_pw_qpolynomial *res;
1351 isl_set *set;
1352
1353 if (!pwqp1 || !pwqp2)
1354 goto error;
1355
1356 isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1357 goto error);
1358
1359 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1360 isl_pw_qpolynomial_free(pwqp2);
1361 return pwqp1;
1362 }
1363
1364 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1365 isl_pw_qpolynomial_free(pwqp1);
1366 return pwqp2;
1367 }
1368
1369 if (isl_pw_qpolynomial_is_one(pwqp1)) {
1370 isl_pw_qpolynomial_free(pwqp1);
1371 return pwqp2;
1372 }
1373
1374 if (isl_pw_qpolynomial_is_one(pwqp2)) {
1375 isl_pw_qpolynomial_free(pwqp2);
1376 return pwqp1;
1377 }
1378
1379 n = pwqp1->n * pwqp2->n;
1380 res = isl_pw_qpolynomial_alloc_(isl_dim_copy(pwqp1->dim), n);
1381
1382 for (i = 0; i < pwqp1->n; ++i) {
1383 for (j = 0; j < pwqp2->n; ++j) {
1384 struct isl_set *common;
1385 struct isl_qpolynomial *prod;
1386 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
1387 isl_set_copy(pwqp2->p[j].set));
1388 if (isl_set_fast_is_empty(common)) {
1389 isl_set_free(common);
1390 continue;
1391 }
1392
1393 prod = isl_qpolynomial_mul(
1394 isl_qpolynomial_copy(pwqp1->p[i].qp),
1395 isl_qpolynomial_copy(pwqp2->p[j].qp));
1396
1397 res = isl_pw_qpolynomial_add_piece(res, common, prod);
1398 }
1399 }
1400
1401 isl_pw_qpolynomial_free(pwqp1);
1402 isl_pw_qpolynomial_free(pwqp2);
1403
1404 return res;
1405 error:
1406 isl_pw_qpolynomial_free(pwqp1);
1407 isl_pw_qpolynomial_free(pwqp2);
1408 return NULL;
1409 }
1410
1411 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_neg(
1412 __isl_take isl_pw_qpolynomial *pwqp)
1413 {
1414 int i, j, n;
1415 struct isl_pw_qpolynomial *res;
1416 isl_set *set;
1417
1418 if (!pwqp)
1419 return NULL;
1420
1421 if (isl_pw_qpolynomial_is_zero(pwqp))
1422 return pwqp;
1423
1424 pwqp = isl_pw_qpolynomial_cow(pwqp);
1425
1426 for (i = 0; i < pwqp->n; ++i) {
1427 pwqp->p[i].qp = isl_qpolynomial_neg(pwqp->p[i].qp);
1428 if (!pwqp->p[i].qp)
1429 goto error;
1430 }
1431
1432 return pwqp;
1433 error:
1434 isl_pw_qpolynomial_free(pwqp);
1435 return NULL;
1436 }
1437
1438 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sub(
1439 __isl_take isl_pw_qpolynomial *pwqp1,
1440 __isl_take isl_pw_qpolynomial *pwqp2)
1441 {
1442 return isl_pw_qpolynomial_add(pwqp1, isl_pw_qpolynomial_neg(pwqp2));
1443 }
1444
1445 __isl_give struct isl_upoly *isl_upoly_eval(
1446 __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
1447 {
1448 int i;
1449 struct isl_upoly_rec *rec;
1450 struct isl_upoly *res;
1451 struct isl_upoly *base;
1452
1453 if (isl_upoly_is_cst(up)) {
1454 isl_vec_free(vec);
1455 return up;
1456 }
1457
1458 rec = isl_upoly_as_rec(up);
1459 if (!rec)
1460 goto error;
1461
1462 isl_assert(up->ctx, rec->n >= 1, goto error);
1463
1464 base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
1465
1466 res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
1467 isl_vec_copy(vec));
1468
1469 for (i = rec->n - 2; i >= 0; --i) {
1470 res = isl_upoly_mul(res, isl_upoly_copy(base));
1471 res = isl_upoly_sum(res,
1472 isl_upoly_eval(isl_upoly_copy(rec->p[i]),
1473 isl_vec_copy(vec)));
1474 }
1475
1476 isl_upoly_free(base);
1477 isl_upoly_free(up);
1478 isl_vec_free(vec);
1479 return res;
1480 error:
1481 isl_upoly_free(up);
1482 isl_vec_free(vec);
1483 return NULL;
1484 }
1485
1486 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
1487 __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
1488 {
1489 isl_vec *ext;
1490 struct isl_upoly *up;
1491 isl_dim *dim;
1492
1493 if (!qp || !pnt)
1494 goto error;
1495 isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, qp->dim), goto error);
1496
1497 if (qp->div->n_row == 0)
1498 ext = isl_vec_copy(pnt->vec);
1499 else {
1500 int i;
1501 unsigned dim = isl_dim_total(qp->dim);
1502 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
1503 if (!ext)
1504 goto error;
1505
1506 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
1507 for (i = 0; i < qp->div->n_row; ++i) {
1508 isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
1509 1 + dim + i, &ext->el[1+dim+i]);
1510 isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
1511 qp->div->row[i][0]);
1512 }
1513 }
1514
1515 up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
1516 if (!up)
1517 goto error;
1518
1519 dim = isl_dim_copy(qp->dim);
1520 isl_qpolynomial_free(qp);
1521 isl_point_free(pnt);
1522
1523 qp = isl_qpolynomial_alloc(dim, 0);
1524 if (!qp)
1525 isl_upoly_free(up);
1526 else
1527 qp->upoly = up;
1528
1529 return qp;
1530 error:
1531 isl_qpolynomial_free(qp);
1532 isl_point_free(pnt);
1533 return NULL;
1534 }
1535
1536 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
1537 __isl_keep struct isl_upoly_cst *cst2)
1538 {
1539 int cmp;
1540 isl_int t;
1541 isl_int_init(t);
1542 isl_int_mul(t, cst1->n, cst2->d);
1543 isl_int_submul(t, cst2->n, cst1->d);
1544 cmp = isl_int_sgn(t);
1545 isl_int_clear(t);
1546 return cmp;
1547 }
1548
1549 static __isl_give isl_qpolynomial *qpolynomial_min(
1550 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1551 {
1552 struct isl_upoly_cst *cst1, *cst2;
1553 int cmp;
1554
1555 if (!qp1 || !qp2)
1556 goto error;
1557 isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
1558 isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
1559 cst1 = isl_upoly_as_cst(qp1->upoly);
1560 cst2 = isl_upoly_as_cst(qp2->upoly);
1561 cmp = isl_upoly_cmp(cst1, cst2);
1562
1563 if (cmp <= 0) {
1564 isl_qpolynomial_free(qp2);
1565 } else {
1566 isl_qpolynomial_free(qp1);
1567 qp1 = qp2;
1568 }
1569 return qp1;
1570 error:
1571 isl_qpolynomial_free(qp1);
1572 isl_qpolynomial_free(qp2);
1573 return NULL;
1574 }
1575
1576 static __isl_give isl_qpolynomial *qpolynomial_max(
1577 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1578 {
1579 struct isl_upoly_cst *cst1, *cst2;
1580 int cmp;
1581
1582 if (!qp1 || !qp2)
1583 goto error;
1584 isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
1585 isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
1586 cst1 = isl_upoly_as_cst(qp1->upoly);
1587 cst2 = isl_upoly_as_cst(qp2->upoly);
1588 cmp = isl_upoly_cmp(cst1, cst2);
1589
1590 if (cmp >= 0) {
1591 isl_qpolynomial_free(qp2);
1592 } else {
1593 isl_qpolynomial_free(qp1);
1594 qp1 = qp2;
1595 }
1596 return qp1;
1597 error:
1598 isl_qpolynomial_free(qp1);
1599 isl_qpolynomial_free(qp2);
1600 return NULL;
1601 }
1602
1603 __isl_give isl_qpolynomial *isl_qpolynomial_fold_eval(
1604 __isl_take isl_qpolynomial_fold *fold, __isl_take isl_point *pnt)
1605 {
1606 isl_qpolynomial *qp;
1607
1608 if (!fold || !pnt)
1609 goto error;
1610 isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, fold->dim), goto error);
1611 isl_assert(pnt->dim->ctx,
1612 fold->type == isl_fold_max || fold->type == isl_fold_min,
1613 goto error);
1614
1615 if (fold->n == 0)
1616 qp = isl_qpolynomial_zero(isl_dim_copy(fold->dim));
1617 else {
1618 int i;
1619 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(fold->qp[0]),
1620 isl_point_copy(pnt));
1621 for (i = 1; i < fold->n; ++i) {
1622 isl_qpolynomial *qp_i;
1623 qp_i = isl_qpolynomial_eval(
1624 isl_qpolynomial_copy(fold->qp[i]),
1625 isl_point_copy(pnt));
1626 if (fold->type == isl_fold_max)
1627 qp = qpolynomial_max(qp, qp_i);
1628 else
1629 qp = qpolynomial_min(qp, qp_i);
1630 }
1631 }
1632 isl_qpolynomial_fold_free(fold);
1633 isl_point_free(pnt);
1634
1635 return qp;
1636 error:
1637 isl_qpolynomial_fold_free(fold);
1638 isl_point_free(pnt);
1639 return NULL;
1640 }
1641
1642 __isl_give isl_qpolynomial *isl_qpolynomial_move(__isl_take isl_qpolynomial *qp,
1643 enum isl_dim_type dst_type, unsigned dst_pos,
1644 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1645 {
1646 if (!qp)
1647 return NULL;
1648
1649 qp->dim = isl_dim_move(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
1650 if (!qp->dim)
1651 goto error;
1652
1653 /* Do something to polynomials when needed; later */
1654
1655 return qp;
1656 error:
1657 isl_qpolynomial_free(qp);
1658 return NULL;
1659 }
1660
1661 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_move(
1662 __isl_take isl_pw_qpolynomial *pwqp,
1663 enum isl_dim_type dst_type, unsigned dst_pos,
1664 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1665 {
1666 int i;
1667
1668 if (!pwqp)
1669 return NULL;
1670
1671 pwqp->dim = isl_dim_move(pwqp->dim,
1672 dst_type, dst_pos, src_type, src_pos, n);
1673 if (!pwqp->dim)
1674 goto error;
1675
1676 for (i = 0; i < pwqp->n; ++i) {
1677 pwqp->p[i].set = isl_set_move(pwqp->p[i].set, dst_type, dst_pos,
1678 src_type, src_pos, n);
1679 if (!pwqp->p[i].set)
1680 goto error;
1681 pwqp->p[i].qp = isl_qpolynomial_move(pwqp->p[i].qp,
1682 dst_type, dst_pos, src_type, src_pos, n);
1683 if (!pwqp->p[i].qp)
1684 goto error;
1685 }
1686
1687 return pwqp;
1688 error:
1689 isl_pw_qpolynomial_free(pwqp);
1690 return NULL;
1691 }
1692
1693 __isl_give isl_dim *isl_pw_qpolynomial_get_dim(
1694 __isl_keep isl_pw_qpolynomial *pwqp)
1695 {
1696 if (!pwqp)
1697 return NULL;
1698
1699 return isl_dim_copy(pwqp->dim);
1700 }
1701
1702 unsigned isl_pw_qpolynomial_dim(__isl_keep isl_pw_qpolynomial *pwqp,
1703 enum isl_dim_type type)
1704 {
1705 return pwqp ? isl_dim_size(pwqp->dim, type) : 0;
1706 }
1707
1708 __isl_give isl_term *isl_term_alloc(__isl_take isl_dim *dim,
1709 __isl_take isl_mat *div)
1710 {
1711 isl_term *term;
1712 int n;
1713
1714 if (!dim || !div)
1715 goto error;
1716
1717 n = isl_dim_total(dim) + div->n_row;
1718
1719 term = isl_calloc(dim->ctx, struct isl_term,
1720 sizeof(struct isl_term) + (n - 1) * sizeof(int));
1721 if (!term)
1722 goto error;
1723
1724 term->ref = 1;
1725 term->dim = dim;
1726 term->div = div;
1727 isl_int_init(term->n);
1728 isl_int_init(term->d);
1729
1730 return term;
1731 error:
1732 isl_dim_free(dim);
1733 isl_mat_free(div);
1734 return NULL;
1735 }
1736
1737 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
1738 {
1739 if (!term)
1740 return NULL;
1741
1742 term->ref++;
1743 return term;
1744 }
1745
1746 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
1747 {
1748 int i;
1749 isl_term *dup;
1750 unsigned total;
1751
1752 if (term)
1753 return NULL;
1754
1755 total = isl_dim_total(term->dim) + term->div->n_row;
1756
1757 dup = isl_term_alloc(isl_dim_copy(term->dim), isl_mat_copy(term->div));
1758 if (!dup)
1759 return NULL;
1760
1761 isl_int_set(dup->n, term->n);
1762 isl_int_set(dup->d, term->d);
1763
1764 for (i = 0; i < total; ++i)
1765 dup->pow[i] = term->pow[i];
1766
1767 return dup;
1768 }
1769
1770 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
1771 {
1772 if (!term)
1773 return NULL;
1774
1775 if (term->ref == 1)
1776 return term;
1777 term->ref--;
1778 return isl_term_dup(term);
1779 }
1780
1781 void isl_term_free(__isl_take isl_term *term)
1782 {
1783 if (!term)
1784 return;
1785
1786 if (--term->ref > 0)
1787 return;
1788
1789 isl_dim_free(term->dim);
1790 isl_mat_free(term->div);
1791 isl_int_clear(term->n);
1792 isl_int_clear(term->d);
1793 free(term);
1794 }
1795
1796 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
1797 {
1798 if (!term)
1799 return 0;
1800
1801 switch (type) {
1802 case isl_dim_param:
1803 case isl_dim_in:
1804 case isl_dim_out: return isl_dim_size(term->dim, type);
1805 case isl_dim_div: return term->div->n_row;
1806 case isl_dim_all: return isl_dim_total(term->dim) + term->div->n_row;
1807 }
1808 }
1809
1810 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
1811 {
1812 return term ? term->dim->ctx : NULL;
1813 }
1814
1815 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
1816 {
1817 if (!term)
1818 return;
1819 isl_int_set(*n, term->n);
1820 }
1821
1822 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
1823 {
1824 if (!term)
1825 return;
1826 isl_int_set(*d, term->d);
1827 }
1828
1829 int isl_term_get_exp(__isl_keep isl_term *term,
1830 enum isl_dim_type type, unsigned pos)
1831 {
1832 if (!term)
1833 return -1;
1834
1835 isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
1836
1837 if (type >= isl_dim_set)
1838 pos += isl_dim_size(term->dim, isl_dim_param);
1839 if (type >= isl_dim_div)
1840 pos += isl_dim_size(term->dim, isl_dim_set);
1841
1842 return term->pow[pos];
1843 }
1844
1845 __isl_give isl_div *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
1846 {
1847 isl_basic_map *bmap;
1848 unsigned total;
1849 int k;
1850
1851 if (!term)
1852 return NULL;
1853
1854 isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
1855 return NULL);
1856
1857 total = term->div->n_col - term->div->n_row - 2;
1858 /* No nested divs for now */
1859 isl_assert(term->dim->ctx,
1860 isl_seq_first_non_zero(term->div->row[pos] + 2 + total,
1861 term->div->n_row) == -1,
1862 return NULL);
1863
1864 bmap = isl_basic_map_alloc_dim(isl_dim_copy(term->dim), 1, 0, 0);
1865 if ((k = isl_basic_map_alloc_div(bmap)) < 0)
1866 goto error;
1867
1868 isl_seq_cpy(bmap->div[k], term->div->row[pos], 2 + total);
1869
1870 return isl_basic_map_div(bmap, k);
1871 error:
1872 isl_basic_map_free(bmap);
1873 return NULL;
1874 }
1875
1876 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
1877 int (*fn)(__isl_take isl_term *term, void *user),
1878 __isl_take isl_term *term, void *user)
1879 {
1880 int i;
1881 struct isl_upoly_rec *rec;
1882
1883 if (!up || !term)
1884 goto error;
1885
1886 if (isl_upoly_is_zero(up))
1887 return term;
1888
1889 isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
1890 isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
1891 isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
1892
1893 if (isl_upoly_is_cst(up)) {
1894 struct isl_upoly_cst *cst;
1895 cst = isl_upoly_as_cst(up);
1896 if (!cst)
1897 goto error;
1898 term = isl_term_cow(term);
1899 if (!term)
1900 goto error;
1901 isl_int_set(term->n, cst->n);
1902 isl_int_set(term->d, cst->d);
1903 if (fn(isl_term_copy(term), user) < 0)
1904 goto error;
1905 return term;
1906 }
1907
1908 rec = isl_upoly_as_rec(up);
1909 if (!rec)
1910 goto error;
1911
1912 for (i = 0; i < rec->n; ++i) {
1913 term = isl_term_cow(term);
1914 if (!term)
1915 goto error;
1916 term->pow[up->var] = i;
1917 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
1918 if (!term)
1919 goto error;
1920 }
1921 term->pow[up->var] = 0;
1922
1923 return term;
1924 error:
1925 isl_term_free(term);
1926 return NULL;
1927 }
1928
1929 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
1930 int (*fn)(__isl_take isl_term *term, void *user), void *user)
1931 {
1932 isl_term *term;
1933
1934 if (!qp)
1935 return -1;
1936
1937 term = isl_term_alloc(isl_dim_copy(qp->dim), isl_mat_copy(qp->div));
1938 if (!term)
1939 return -1;
1940
1941 term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
1942
1943 isl_term_free(term);
1944
1945 return term ? 0 : -1;
1946 }
1947
1948 int isl_pw_qpolynomial_foreach_piece(__isl_keep isl_pw_qpolynomial *pwqp,
1949 int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
1950 void *user), void *user)
1951 {
1952 int i;
1953
1954 if (!pwqp)
1955 return -1;
1956
1957 for (i = 0; i < pwqp->n; ++i)
1958 if (fn(isl_set_copy(pwqp->p[i].set),
1959 isl_qpolynomial_copy(pwqp->p[i].qp), user) < 0)
1960 return -1;
1961
1962 return 0;
1963 }
1964
1965 static int any_divs(__isl_keep isl_set *set)
1966 {
1967 int i;
1968
1969 if (!set)
1970 return -1;
1971
1972 for (i = 0; i < set->n; ++i)
1973 if (set->p[i]->n_div > 0)
1974 return 1;
1975
1976 return 0;
1977 }
1978
1979 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
1980 __isl_take isl_dim *dim)
1981 {
1982 if (!qp || !dim)
1983 goto error;
1984
1985 if (isl_dim_equal(qp->dim, dim)) {
1986 isl_dim_free(dim);
1987 return qp;
1988 }
1989
1990 qp = isl_qpolynomial_cow(qp);
1991 if (!qp)
1992 goto error;
1993
1994 if (qp->div->n_row) {
1995 int i;
1996 int extra;
1997 unsigned total;
1998 int *exp;
1999
2000 extra = isl_dim_size(dim, isl_dim_set) -
2001 isl_dim_size(qp->dim, isl_dim_set);
2002 total = isl_dim_total(qp->dim);
2003 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
2004 if (!exp)
2005 goto error;
2006 for (i = 0; i < qp->div->n_row; ++i)
2007 exp[i] = extra + i;
2008 qp->upoly = expand(qp->upoly, exp, total);
2009 free(exp);
2010 if (!qp->upoly)
2011 goto error;
2012 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
2013 if (!qp->div)
2014 goto error;
2015 for (i = 0; i < qp->div->n_row; ++i)
2016 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
2017 }
2018
2019 isl_dim_free(qp->dim);
2020 qp->dim = dim;
2021
2022 return qp;
2023 error:
2024 isl_dim_free(dim);
2025 isl_qpolynomial_free(qp);
2026 return NULL;
2027 }
2028
2029 static int foreach_lifted_subset(__isl_take isl_set *set,
2030 __isl_take isl_qpolynomial *qp,
2031 int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2032 void *user), void *user)
2033 {
2034 int i;
2035
2036 if (!set || !qp)
2037 goto error;
2038
2039 for (i = 0; i < set->n; ++i) {
2040 isl_set *lift;
2041 isl_qpolynomial *copy;
2042
2043 lift = isl_set_from_basic_set(isl_basic_set_copy(set->p[i]));
2044 lift = isl_set_lift(lift);
2045
2046 copy = isl_qpolynomial_copy(qp);
2047 copy = isl_qpolynomial_lift(copy, isl_set_get_dim(lift));
2048
2049 if (fn(lift, copy, user) < 0)
2050 goto error;
2051 }
2052
2053 isl_set_free(set);
2054 isl_qpolynomial_free(qp);
2055
2056 return 0;
2057 error:
2058 isl_set_free(set);
2059 isl_qpolynomial_free(qp);
2060 return -1;
2061 }
2062
2063 int isl_pw_qpolynomial_foreach_lifted_piece(__isl_keep isl_pw_qpolynomial *pwqp,
2064 int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2065 void *user), void *user)
2066 {
2067 int i;
2068
2069 if (!pwqp)
2070 return -1;
2071
2072 for (i = 0; i < pwqp->n; ++i) {
2073 isl_set *set;
2074 isl_qpolynomial *qp;
2075
2076 set = isl_set_copy(pwqp->p[i].set);
2077 qp = isl_qpolynomial_copy(pwqp->p[i].qp);
2078 if (!any_divs(set)) {
2079 if (fn(set, qp, user) < 0)
2080 return -1;
2081 continue;
2082 }
2083 if (foreach_lifted_subset(set, qp, fn, user) < 0)
2084 return -1;
2085 }
2086
2087 return 0;
2088 }
2089
2090 static __isl_give isl_qpolynomial_fold *qpolynomial_fold_alloc(
2091 enum isl_fold type, __isl_take isl_dim *dim, int n)
2092 {
2093 isl_qpolynomial_fold *fold;
2094
2095 if (!dim)
2096 goto error;
2097
2098 isl_assert(dim->ctx, n >= 0, goto error);
2099 fold = isl_alloc(dim->ctx, struct isl_qpolynomial_fold,
2100 sizeof(struct isl_qpolynomial_fold) +
2101 (n - 1) * sizeof(struct isl_qpolynomial *));
2102 if (!fold)
2103 goto error;
2104
2105 fold->ref = 1;
2106 fold->size = n;
2107 fold->n = 0;
2108 fold->type = type;
2109 fold->dim = dim;
2110
2111 return fold;
2112 error:
2113 isl_dim_free(dim);
2114 return NULL;
2115 }
2116
2117 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_empty(enum isl_fold type,
2118 __isl_take isl_dim *dim)
2119 {
2120 return qpolynomial_fold_alloc(type, dim, 0);
2121 }
2122
2123 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_alloc(
2124 enum isl_fold type, __isl_take isl_qpolynomial *qp)
2125 {
2126 isl_qpolynomial_fold *fold;
2127
2128 if (!qp)
2129 return NULL;
2130
2131 fold = qpolynomial_fold_alloc(type, isl_dim_copy(qp->dim), 1);
2132 if (!fold)
2133 goto error;
2134
2135 fold->qp[0] = qp;
2136 fold->n++;
2137
2138 return fold;
2139 error:
2140 isl_qpolynomial_fold_free(fold);
2141 isl_qpolynomial_free(qp);
2142 return NULL;
2143 }
2144
2145 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(
2146 __isl_keep isl_qpolynomial_fold *fold)
2147 {
2148 if (!fold)
2149 return NULL;
2150
2151 fold->ref++;
2152 return fold;
2153 }
2154
2155 void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold *fold)
2156 {
2157 int i;
2158
2159 if (!fold)
2160 return;
2161 if (--fold->ref > 0)
2162 return;
2163
2164 for (i = 0; i < fold->n; ++i)
2165 isl_qpolynomial_free(fold->qp[i]);
2166 isl_dim_free(fold->dim);
2167 free(fold);
2168 }
2169
2170 int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold *fold)
2171 {
2172 if (!fold)
2173 return -1;
2174
2175 return fold->n == 0;
2176 }
2177
2178 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold(
2179 __isl_take isl_qpolynomial_fold *fold1,
2180 __isl_take isl_qpolynomial_fold *fold2)
2181 {
2182 int i;
2183 struct isl_qpolynomial_fold *res = NULL;
2184
2185 if (!fold1 || !fold2)
2186 goto error;
2187
2188 isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
2189 isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
2190 goto error);
2191
2192 if (isl_qpolynomial_fold_is_empty(fold1)) {
2193 isl_qpolynomial_fold_free(fold1);
2194 return fold2;
2195 }
2196
2197 if (isl_qpolynomial_fold_is_empty(fold2)) {
2198 isl_qpolynomial_fold_free(fold2);
2199 return fold1;
2200 }
2201
2202 res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
2203 fold1->n + fold2->n);
2204 if (!res)
2205 goto error;
2206
2207 for (i = 0; i < fold1->n; ++i) {
2208 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
2209 if (!res->qp[res->n])
2210 goto error;
2211 res->n++;
2212 }
2213
2214 for (i = 0; i < fold2->n; ++i) {
2215 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
2216 if (!res->qp[res->n])
2217 goto error;
2218 res->n++;
2219 }
2220
2221 isl_qpolynomial_fold_free(fold1);
2222 isl_qpolynomial_fold_free(fold2);
2223
2224 return res;
2225 error:
2226 isl_qpolynomial_fold_free(res);
2227 isl_qpolynomial_fold_free(fold1);
2228 isl_qpolynomial_fold_free(fold2);
2229 return NULL;
2230 }
2231
2232 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_from_pw_qpolynomial(
2233 enum isl_fold type, __isl_take isl_pw_qpolynomial *pwqp)
2234 {
2235 int i;
2236 isl_pw_qpolynomial_fold *pwf;
2237
2238 if (!pwqp)
2239 return NULL;
2240
2241 pwf = isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pwqp->dim), pwqp->n);
2242
2243 for (i = 0; i < pwqp->n; ++i)
2244 pwf = isl_pw_qpolynomial_fold_add_piece(pwf,
2245 isl_set_copy(pwqp->p[i].set),
2246 isl_qpolynomial_fold_alloc(type,
2247 isl_qpolynomial_copy(pwqp->p[i].qp)));
2248
2249 isl_pw_qpolynomial_free(pwqp);
2250
2251 return pwf;
2252 }
Integer set library