?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
add isl_pw_qpolynomial_foreach_lifted_piece
[isl.git] / isl_polynomial.c
blob0b2e3c0c5ad9eb344c87d528234f7c53d8315b25
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 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
420 __isl_take struct isl_upoly *up2)
421 {
422 int i;
423 struct isl_upoly_rec *rec1, *rec2;
424
425 if (!up1 || !up2)
426 goto error;
427
428 if (isl_upoly_is_nan(up1)) {
429 isl_upoly_free(up2);
430 return up1;
431 }
432
433 if (isl_upoly_is_nan(up2)) {
434 isl_upoly_free(up1);
435 return up2;
436 }
437
438 if (isl_upoly_is_zero(up1)) {
439 isl_upoly_free(up1);
440 return up2;
441 }
442
443 if (isl_upoly_is_zero(up2)) {
444 isl_upoly_free(up2);
445 return up1;
446 }
447
448 if (up1->var < up2->var)
449 return isl_upoly_sum(up2, up1);
450
451 if (up2->var < up1->var) {
452 struct isl_upoly_rec *rec;
453 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
454 isl_upoly_free(up1);
455 return up2;
456 }
457 up1 = isl_upoly_cow(up1);
458 rec = isl_upoly_as_rec(up1);
459 if (!rec)
460 goto error;
461 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
462 if (rec->n == 1 && isl_upoly_is_zero(rec->p[0]))
463 up1 = replace_by_zero(up1);
464 return up1;
465 }
466
467 if (isl_upoly_is_cst(up1))
468 return isl_upoly_sum_cst(up1, up2);
469
470 rec1 = isl_upoly_as_rec(up1);
471 rec2 = isl_upoly_as_rec(up2);
472 if (!rec1 || !rec2)
473 goto error;
474
475 if (rec1->n < rec2->n)
476 return isl_upoly_sum(up2, up1);
477
478 up1 = isl_upoly_cow(up1);
479 rec1 = isl_upoly_as_rec(up1);
480 if (!rec1)
481 goto error;
482
483 for (i = rec2->n - 1; i >= 0; --i) {
484 rec1->p[i] = isl_upoly_sum(rec1->p[i],
485 isl_upoly_copy(rec2->p[i]));
486 if (!rec1->p[i])
487 goto error;
488 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
489 isl_upoly_free(rec1->p[i]);
490 rec1->n--;
491 }
492 }
493
494 if (rec1->n == 0)
495 up1 = replace_by_zero(up1);
496
497 isl_upoly_free(up2);
498
499 return up1;
500 error:
501 isl_upoly_free(up1);
502 isl_upoly_free(up2);
503 return NULL;
504 }
505
506 __isl_give struct isl_upoly *isl_upoly_neg_cst(__isl_take struct isl_upoly *up)
507 {
508 struct isl_upoly_cst *cst;
509
510 if (isl_upoly_is_zero(up))
511 return up;
512
513 up = isl_upoly_cow(up);
514 if (!up)
515 return NULL;
516
517 cst = isl_upoly_as_cst(up);
518
519 isl_int_neg(cst->n, cst->n);
520
521 return up;
522 }
523
524 __isl_give struct isl_upoly *isl_upoly_neg(__isl_take struct isl_upoly *up)
525 {
526 int i;
527 struct isl_upoly_rec *rec;
528
529 if (!up)
530 return NULL;
531
532 if (isl_upoly_is_cst(up))
533 return isl_upoly_neg_cst(up);
534
535 up = isl_upoly_cow(up);
536 rec = isl_upoly_as_rec(up);
537 if (!rec)
538 goto error;
539
540 for (i = 0; i < rec->n; ++i) {
541 rec->p[i] = isl_upoly_neg(rec->p[i]);
542 if (!rec->p[i])
543 goto error;
544 }
545
546 return up;
547 error:
548 isl_upoly_free(up);
549 return NULL;
550 }
551
552 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
553 __isl_take struct isl_upoly *up2)
554 {
555 struct isl_upoly_cst *cst1;
556 struct isl_upoly_cst *cst2;
557
558 up1 = isl_upoly_cow(up1);
559 if (!up1 || !up2)
560 goto error;
561
562 cst1 = isl_upoly_as_cst(up1);
563 cst2 = isl_upoly_as_cst(up2);
564
565 isl_int_mul(cst1->n, cst1->n, cst2->n);
566 isl_int_mul(cst1->d, cst1->d, cst2->d);
567
568 isl_upoly_cst_reduce(cst1);
569
570 isl_upoly_free(up2);
571 return up1;
572 error:
573 isl_upoly_free(up1);
574 isl_upoly_free(up2);
575 return NULL;
576 }
577
578 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
579 __isl_take struct isl_upoly *up2)
580 {
581 struct isl_upoly_rec *rec1;
582 struct isl_upoly_rec *rec2;
583 struct isl_upoly_rec *res;
584 int i, j;
585 int size;
586
587 rec1 = isl_upoly_as_rec(up1);
588 rec2 = isl_upoly_as_rec(up2);
589 if (!rec1 || !rec2)
590 goto error;
591 size = rec1->n + rec2->n - 1;
592 res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
593 if (!res)
594 goto error;
595
596 for (i = 0; i < rec1->n; ++i) {
597 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
598 isl_upoly_copy(rec1->p[i]));
599 if (!res->p[i])
600 goto error;
601 res->n++;
602 }
603 for (; i < size; ++i) {
604 res->p[i] = isl_upoly_zero(up1->ctx);
605 if (!res->p[i])
606 goto error;
607 res->n++;
608 }
609 for (i = 0; i < rec1->n; ++i) {
610 for (j = 1; j < rec2->n; ++j) {
611 struct isl_upoly *up;
612 up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
613 isl_upoly_copy(rec1->p[i]));
614 res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
615 if (!res->p[i + j])
616 goto error;
617 }
618 }
619
620 isl_upoly_free(up1);
621 isl_upoly_free(up2);
622
623 return &res->up;
624 error:
625 isl_upoly_free(up1);
626 isl_upoly_free(up2);
627 isl_upoly_free(&res->up);
628 return NULL;
629 }
630
631 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
632 __isl_take struct isl_upoly *up2)
633 {
634 if (!up1 || !up2)
635 goto error;
636
637 if (isl_upoly_is_nan(up1)) {
638 isl_upoly_free(up2);
639 return up1;
640 }
641
642 if (isl_upoly_is_nan(up2)) {
643 isl_upoly_free(up1);
644 return up2;
645 }
646
647 if (isl_upoly_is_zero(up1)) {
648 isl_upoly_free(up2);
649 return up1;
650 }
651
652 if (isl_upoly_is_zero(up2)) {
653 isl_upoly_free(up1);
654 return up2;
655 }
656
657 if (isl_upoly_is_one(up1)) {
658 isl_upoly_free(up1);
659 return up2;
660 }
661
662 if (isl_upoly_is_one(up2)) {
663 isl_upoly_free(up2);
664 return up1;
665 }
666
667 if (up1->var < up2->var)
668 return isl_upoly_mul(up2, up1);
669
670 if (up2->var < up1->var) {
671 int i;
672 struct isl_upoly_rec *rec;
673 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
674 isl_ctx *ctx = up1->ctx;
675 isl_upoly_free(up1);
676 isl_upoly_free(up2);
677 return isl_upoly_nan(ctx);
678 }
679 up1 = isl_upoly_cow(up1);
680 rec = isl_upoly_as_rec(up1);
681 if (!rec)
682 goto error;
683
684 for (i = 0; i < rec->n; ++i) {
685 rec->p[i] = isl_upoly_mul(rec->p[i],
686 isl_upoly_copy(up2));
687 if (!rec->p[i])
688 goto error;
689 }
690 isl_upoly_free(up2);
691 return up1;
692 }
693
694 if (isl_upoly_is_cst(up1))
695 return isl_upoly_mul_cst(up1, up2);
696
697 return isl_upoly_mul_rec(up1, up2);
698 error:
699 isl_upoly_free(up1);
700 isl_upoly_free(up2);
701 return NULL;
702 }
703
704 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_dim *dim,
705 unsigned n_div)
706 {
707 struct isl_qpolynomial *qp;
708 unsigned total;
709
710 if (!dim)
711 return NULL;
712
713 total = isl_dim_total(dim);
714
715 qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
716 if (!qp)
717 return NULL;
718
719 qp->ref = 1;
720 qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
721 if (!qp->div)
722 goto error;
723
724 qp->dim = dim;
725
726 return qp;
727 error:
728 isl_dim_free(dim);
729 isl_qpolynomial_free(qp);
730 return NULL;
731 }
732
733 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
734 {
735 if (!qp)
736 return NULL;
737
738 qp->ref++;
739 return qp;
740 }
741
742 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
743 {
744 struct isl_qpolynomial *dup;
745
746 if (!qp)
747 return NULL;
748
749 dup = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), qp->div->n_row);
750 if (!dup)
751 return NULL;
752 isl_mat_free(dup->div);
753 dup->div = isl_mat_copy(qp->div);
754 if (!dup->div)
755 goto error;
756 dup->upoly = isl_upoly_copy(qp->upoly);
757 if (!dup->upoly)
758 goto error;
759
760 return dup;
761 error:
762 isl_qpolynomial_free(dup);
763 return NULL;
764 }
765
766 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
767 {
768 if (!qp)
769 return NULL;
770
771 if (qp->ref == 1)
772 return qp;
773 qp->ref--;
774 return isl_qpolynomial_dup(qp);
775 }
776
777 void isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
778 {
779 if (!qp)
780 return;
781
782 if (--qp->ref > 0)
783 return;
784
785 isl_dim_free(qp->dim);
786 isl_mat_free(qp->div);
787 isl_upoly_free(qp->upoly);
788
789 free(qp);
790 }
791
792 static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
793 {
794 int n_row, n_col;
795 int equal;
796
797 isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
798 div1->n_col >= div2->n_col, return -1);
799
800 if (div1->n_row == div2->n_row)
801 return isl_mat_is_equal(div1, div2);
802
803 n_row = div1->n_row;
804 n_col = div1->n_col;
805 div1->n_row = div2->n_row;
806 div1->n_col = div2->n_col;
807
808 equal = isl_mat_is_equal(div1, div2);
809
810 div1->n_row = n_row;
811 div1->n_col = n_col;
812
813 return equal;
814 }
815
816 static void expand_row(__isl_keep isl_mat *dst, int d,
817 __isl_keep isl_mat *src, int s, int *exp)
818 {
819 int i;
820 unsigned c = src->n_col - src->n_row;
821
822 isl_seq_cpy(dst->row[d], src->row[s], c);
823 isl_seq_clr(dst->row[d] + c, dst->n_col - c);
824
825 for (i = 0; i < s; ++i)
826 isl_int_set(dst->row[d][c + exp[i]], src->row[s][c + i]);
827 }
828
829 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
830 {
831 int li, lj;
832
833 li = isl_seq_last_non_zero(div->row[i], div->n_col);
834 lj = isl_seq_last_non_zero(div->row[j], div->n_col);
835
836 if (li != lj)
837 return li - lj;
838
839 return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
840 }
841
842 static __isl_give isl_mat *merge_divs(__isl_keep isl_mat *div1,
843 __isl_keep isl_mat *div2, int *exp1, int *exp2)
844 {
845 int i, j, k;
846 isl_mat *div = NULL;
847 unsigned d = div1->n_col - div1->n_row;
848
849 div = isl_mat_alloc(div1->ctx, 1 + div1->n_row + div2->n_row,
850 d + div1->n_row + div2->n_row);
851 if (!div)
852 return NULL;
853
854 for (i = 0, j = 0, k = 0; i < div1->n_row && j < div2->n_row; ++k) {
855 int cmp;
856
857 expand_row(div, k, div1, i, exp1);
858 expand_row(div, k + 1, div2, j, exp2);
859
860 cmp = cmp_row(div, k, k + 1);
861 if (cmp == 0) {
862 exp1[i++] = k;
863 exp2[j++] = k;
864 } else if (cmp < 0) {
865 exp1[i++] = k;
866 } else {
867 exp2[j++] = k;
868 isl_seq_cpy(div->row[k], div->row[k + 1], div->n_col);
869 }
870 }
871 for (; i < div1->n_row; ++i, ++k) {
872 expand_row(div, k, div1, i, exp1);
873 exp1[i] = k;
874 }
875 for (; j < div2->n_row; ++j, ++k) {
876 expand_row(div, k, div2, j, exp2);
877 exp2[j] = k;
878 }
879
880 div->n_row = k;
881 div->n_col = d + k;
882
883 return div;
884 }
885
886 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
887 int *exp, int first)
888 {
889 int i;
890 struct isl_upoly_rec *rec;
891
892 if (isl_upoly_is_cst(up))
893 return up;
894
895 if (up->var < first)
896 return up;
897
898 if (exp[up->var - first] == up->var - first)
899 return up;
900
901 up = isl_upoly_cow(up);
902 if (!up)
903 goto error;
904
905 up->var = exp[up->var - first] + first;
906
907 rec = isl_upoly_as_rec(up);
908 if (!rec)
909 goto error;
910
911 for (i = 0; i < rec->n; ++i) {
912 rec->p[i] = expand(rec->p[i], exp, first);
913 if (!rec->p[i])
914 goto error;
915 }
916
917 return up;
918 error:
919 isl_upoly_free(up);
920 return NULL;
921 }
922
923 static __isl_give isl_qpolynomial *with_merged_divs(
924 __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
925 __isl_take isl_qpolynomial *qp2),
926 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
927 {
928 int *exp1 = NULL;
929 int *exp2 = NULL;
930 isl_mat *div = NULL;
931
932 qp1 = isl_qpolynomial_cow(qp1);
933 qp2 = isl_qpolynomial_cow(qp2);
934
935 if (!qp1 || !qp2)
936 goto error;
937
938 isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
939 qp1->div->n_col >= qp2->div->n_col, goto error);
940
941 exp1 = isl_alloc_array(qp1->div->ctx, int, qp1->div->n_row);
942 exp2 = isl_alloc_array(qp2->div->ctx, int, qp2->div->n_row);
943 if (!exp1 || !exp2)
944 goto error;
945
946 div = merge_divs(qp1->div, qp2->div, exp1, exp2);
947 if (!div)
948 goto error;
949
950 isl_mat_free(qp1->div);
951 qp1->div = isl_mat_copy(div);
952 isl_mat_free(qp2->div);
953 qp2->div = isl_mat_copy(div);
954
955 qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
956 qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
957
958 if (!qp1->upoly || !qp2->upoly)
959 goto error;
960
961 isl_mat_free(div);
962 free(exp1);
963 free(exp2);
964
965 return fn(qp1, qp2);
966 error:
967 isl_mat_free(div);
968 free(exp1);
969 free(exp2);
970 isl_qpolynomial_free(qp1);
971 isl_qpolynomial_free(qp2);
972 return NULL;
973 }
974
975 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
976 __isl_take isl_qpolynomial *qp2)
977 {
978 qp1 = isl_qpolynomial_cow(qp1);
979
980 if (!qp1 || !qp2)
981 goto error;
982
983 if (qp1->div->n_row < qp2->div->n_row)
984 return isl_qpolynomial_add(qp2, qp1);
985
986 isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
987 if (!compatible_divs(qp1->div, qp2->div))
988 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
989
990 qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
991 if (!qp1->upoly)
992 goto error;
993
994 isl_qpolynomial_free(qp2);
995
996 return qp1;
997 error:
998 isl_qpolynomial_free(qp1);
999 isl_qpolynomial_free(qp2);
1000 return NULL;
1001 }
1002
1003 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1004 __isl_take isl_qpolynomial *qp2)
1005 {
1006 return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1007 }
1008
1009 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1010 {
1011 qp = isl_qpolynomial_cow(qp);
1012
1013 if (!qp)
1014 return NULL;
1015
1016 qp->upoly = isl_upoly_neg(qp->upoly);
1017 if (!qp->upoly)
1018 goto error;
1019
1020 return qp;
1021 error:
1022 isl_qpolynomial_free(qp);
1023 return NULL;
1024 }
1025
1026 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1027 __isl_take isl_qpolynomial *qp2)
1028 {
1029 qp1 = isl_qpolynomial_cow(qp1);
1030
1031 if (!qp1 || !qp2)
1032 goto error;
1033
1034 if (qp1->div->n_row < qp2->div->n_row)
1035 return isl_qpolynomial_mul(qp2, qp1);
1036
1037 isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1038 if (!compatible_divs(qp1->div, qp2->div))
1039 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1040
1041 qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1042 if (!qp1->upoly)
1043 goto error;
1044
1045 isl_qpolynomial_free(qp2);
1046
1047 return qp1;
1048 error:
1049 isl_qpolynomial_free(qp1);
1050 isl_qpolynomial_free(qp2);
1051 return NULL;
1052 }
1053
1054 __isl_give isl_qpolynomial *isl_qpolynomial_zero(__isl_take isl_dim *dim)
1055 {
1056 struct isl_qpolynomial *qp;
1057 struct isl_upoly_cst *cst;
1058
1059 qp = isl_qpolynomial_alloc(dim, 0);
1060 if (!qp)
1061 return NULL;
1062
1063 qp->upoly = isl_upoly_zero(dim->ctx);
1064 if (!qp->upoly)
1065 goto error;
1066
1067 return qp;
1068 error:
1069 isl_qpolynomial_free(qp);
1070 return NULL;
1071 }
1072
1073 __isl_give isl_qpolynomial *isl_qpolynomial_infty(__isl_take isl_dim *dim)
1074 {
1075 struct isl_qpolynomial *qp;
1076 struct isl_upoly_cst *cst;
1077
1078 qp = isl_qpolynomial_alloc(dim, 0);
1079 if (!qp)
1080 return NULL;
1081
1082 qp->upoly = isl_upoly_infty(dim->ctx);
1083 if (!qp->upoly)
1084 goto error;
1085
1086 return qp;
1087 error:
1088 isl_qpolynomial_free(qp);
1089 return NULL;
1090 }
1091
1092 __isl_give isl_qpolynomial *isl_qpolynomial_nan(__isl_take isl_dim *dim)
1093 {
1094 struct isl_qpolynomial *qp;
1095 struct isl_upoly_cst *cst;
1096
1097 qp = isl_qpolynomial_alloc(dim, 0);
1098 if (!qp)
1099 return NULL;
1100
1101 qp->upoly = isl_upoly_nan(dim->ctx);
1102 if (!qp->upoly)
1103 goto error;
1104
1105 return qp;
1106 error:
1107 isl_qpolynomial_free(qp);
1108 return NULL;
1109 }
1110
1111 __isl_give isl_qpolynomial *isl_qpolynomial_cst(__isl_take isl_dim *dim,
1112 isl_int v)
1113 {
1114 struct isl_qpolynomial *qp;
1115 struct isl_upoly_cst *cst;
1116
1117 qp = isl_qpolynomial_alloc(dim, 0);
1118 if (!qp)
1119 return NULL;
1120
1121 qp->upoly = isl_upoly_zero(dim->ctx);
1122 if (!qp->upoly)
1123 goto error;
1124 cst = isl_upoly_as_cst(qp->upoly);
1125 isl_int_set(cst->n, v);
1126
1127 return qp;
1128 error:
1129 isl_qpolynomial_free(qp);
1130 return NULL;
1131 }
1132
1133 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1134 isl_int *n, isl_int *d)
1135 {
1136 struct isl_upoly_cst *cst;
1137
1138 if (!qp)
1139 return -1;
1140
1141 if (!isl_upoly_is_cst(qp->upoly))
1142 return 0;
1143
1144 cst = isl_upoly_as_cst(qp->upoly);
1145 if (!cst)
1146 return -1;
1147
1148 if (n)
1149 isl_int_set(*n, cst->n);
1150 if (d)
1151 isl_int_set(*d, cst->d);
1152
1153 return 1;
1154 }
1155
1156 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_dim *dim,
1157 int pos, int power)
1158 {
1159 int i;
1160 struct isl_qpolynomial *qp;
1161 struct isl_upoly_rec *rec;
1162 struct isl_upoly_cst *cst;
1163 struct isl_ctx *ctx;
1164
1165 if (!dim)
1166 return NULL;
1167
1168 ctx = dim->ctx;
1169
1170 qp = isl_qpolynomial_alloc(dim, 0);
1171 if (!qp)
1172 return NULL;
1173
1174 qp->upoly = &isl_upoly_alloc_rec(ctx, pos, 1 + power)->up;
1175 if (!qp->upoly)
1176 goto error;
1177 rec = isl_upoly_as_rec(qp->upoly);
1178 for (i = 0; i < 1 + power; ++i) {
1179 rec->p[i] = isl_upoly_zero(ctx);
1180 if (!rec->p[i])
1181 goto error;
1182 rec->n++;
1183 }
1184 cst = isl_upoly_as_cst(rec->p[power]);
1185 isl_int_set_si(cst->n, 1);
1186
1187 return qp;
1188 error:
1189 isl_qpolynomial_free(qp);
1190 return NULL;
1191 }
1192
1193 __isl_give isl_qpolynomial *isl_qpolynomial_var(__isl_take isl_dim *dim,
1194 enum isl_dim_type type, unsigned pos)
1195 {
1196 if (!dim)
1197 return NULL;
1198
1199 isl_assert(dim->ctx, isl_dim_size(dim, isl_dim_in) == 0, goto error);
1200 isl_assert(dim->ctx, pos < isl_dim_size(dim, type), goto error);
1201
1202 if (type == isl_dim_set)
1203 pos += isl_dim_size(dim, isl_dim_param);
1204
1205 return isl_qpolynomial_pow(dim, pos, 1);
1206 error:
1207 isl_dim_free(dim);
1208 return NULL;
1209 }
1210
1211 __isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
1212 int power)
1213 {
1214 struct isl_qpolynomial *qp = NULL;
1215 struct isl_upoly_rec *rec;
1216 struct isl_upoly_cst *cst;
1217 int i;
1218 int pos;
1219
1220 if (!div)
1221 return NULL;
1222 isl_assert(div->ctx, div->bmap->n_div == 1, goto error);
1223
1224 qp = isl_qpolynomial_alloc(isl_basic_map_get_dim(div->bmap), 1);
1225 if (!qp)
1226 goto error;
1227
1228 isl_seq_cpy(qp->div->row[0], div->line[0], qp->div->n_col - 1);
1229 isl_int_set_si(qp->div->row[0][qp->div->n_col - 1], 0);
1230
1231 pos = isl_dim_total(qp->dim);
1232 qp->upoly = &isl_upoly_alloc_rec(div->ctx, pos, 1 + power)->up;
1233 if (!qp->upoly)
1234 goto error;
1235 rec = isl_upoly_as_rec(qp->upoly);
1236 for (i = 0; i < 1 + power; ++i) {
1237 rec->p[i] = isl_upoly_zero(div->ctx);
1238 if (!rec->p[i])
1239 goto error;
1240 rec->n++;
1241 }
1242 cst = isl_upoly_as_cst(rec->p[power]);
1243 isl_int_set_si(cst->n, 1);
1244
1245 isl_div_free(div);
1246
1247 return qp;
1248 error:
1249 isl_qpolynomial_free(qp);
1250 isl_div_free(div);
1251 return NULL;
1252 }
1253
1254 __isl_give isl_qpolynomial *isl_qpolynomial_div(__isl_take isl_div *div)
1255 {
1256 return isl_qpolynomial_div_pow(div, 1);
1257 }
1258
1259 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst(__isl_take isl_dim *dim,
1260 const isl_int n, const isl_int d)
1261 {
1262 struct isl_qpolynomial *qp;
1263 struct isl_upoly_cst *cst;
1264
1265 qp = isl_qpolynomial_alloc(dim, 0);
1266 if (!qp)
1267 return NULL;
1268
1269 qp->upoly = isl_upoly_zero(dim->ctx);
1270 if (!qp->upoly)
1271 goto error;
1272 cst = isl_upoly_as_cst(qp->upoly);
1273 isl_int_set(cst->n, n);
1274 isl_int_set(cst->d, d);
1275
1276 return qp;
1277 error:
1278 isl_qpolynomial_free(qp);
1279 return NULL;
1280 }
1281
1282 static __isl_give isl_pw_qpolynomial *pw_qpolynomial_alloc(__isl_take isl_dim *dim,
1283 int n)
1284 {
1285 struct isl_pw_qpolynomial *pwqp;
1286
1287 if (!dim)
1288 return NULL;
1289 isl_assert(dim->ctx, n >= 0, return NULL);
1290 pwqp = isl_alloc(dim->ctx, struct isl_pw_qpolynomial,
1291 sizeof(struct isl_pw_qpolynomial) +
1292 (n - 1) * sizeof(struct isl_pw_qpolynomial_piece));
1293 if (!pwqp)
1294 goto error;
1295
1296 pwqp->ref = 1;
1297 pwqp->size = n;
1298 pwqp->n = 0;
1299 pwqp->dim = dim;
1300 return pwqp;
1301 error:
1302 isl_dim_free(dim);
1303 return NULL;
1304 }
1305
1306 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_piece(
1307 __isl_take isl_pw_qpolynomial *pwqp,
1308 __isl_take isl_set *set, __isl_take isl_qpolynomial *qp)
1309 {
1310 if (!pwqp || !set || !qp)
1311 goto error;
1312
1313 if (isl_set_fast_is_empty(set) || isl_qpolynomial_is_zero(qp)) {
1314 isl_set_free(set);
1315 isl_qpolynomial_free(qp);
1316 return pwqp;
1317 }
1318
1319 isl_assert(set->ctx, isl_dim_equal(pwqp->dim, qp->dim), goto error);
1320 isl_assert(set->ctx, pwqp->n < pwqp->size, goto error);
1321
1322 pwqp->p[pwqp->n].set = set;
1323 pwqp->p[pwqp->n].qp = qp;
1324 pwqp->n++;
1325
1326 return pwqp;
1327 error:
1328 isl_pw_qpolynomial_free(pwqp);
1329 isl_set_free(set);
1330 isl_qpolynomial_free(qp);
1331 return NULL;
1332 }
1333
1334 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_alloc(__isl_take isl_set *set,
1335 __isl_take isl_qpolynomial *qp)
1336 {
1337 struct isl_pw_qpolynomial *pwqp;
1338
1339 if (!set || !qp)
1340 goto error;
1341
1342 pwqp = pw_qpolynomial_alloc(isl_set_get_dim(set), 1);
1343
1344 return isl_pw_qpolynomial_add_piece(pwqp, set, qp);
1345 error:
1346 isl_set_free(set);
1347 isl_qpolynomial_free(qp);
1348 return NULL;
1349 }
1350
1351 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_zero(__isl_take isl_dim *dim)
1352 {
1353 return pw_qpolynomial_alloc(dim, 0);
1354 }
1355
1356 int isl_pw_qpolynomial_is_zero(__isl_keep isl_pw_qpolynomial *pwqp)
1357 {
1358 if (!pwqp)
1359 return -1;
1360
1361 return pwqp->n == 0;
1362 }
1363
1364 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
1365 {
1366 if (!pwqp)
1367 return -1;
1368
1369 if (pwqp->n != -1)
1370 return 0;
1371
1372 if (!isl_set_fast_is_universe(pwqp->p[0].set))
1373 return 0;
1374
1375 return isl_qpolynomial_is_one(pwqp->p[0].qp);
1376 }
1377
1378 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_copy(
1379 __isl_keep isl_pw_qpolynomial *pwqp)
1380 {
1381 if (!pwqp)
1382 return;
1383
1384 pwqp->ref++;
1385 return pwqp;
1386 }
1387
1388 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_dup(
1389 __isl_keep isl_pw_qpolynomial *pwqp)
1390 {
1391 int i;
1392 struct isl_pw_qpolynomial *dup;
1393
1394 if (!pwqp)
1395 return NULL;
1396
1397 dup = pw_qpolynomial_alloc(isl_dim_copy(pwqp->dim), pwqp->n);
1398 if (!dup)
1399 return NULL;
1400
1401 for (i = 0; i < pwqp->n; ++i)
1402 dup = isl_pw_qpolynomial_add_piece(dup,
1403 isl_set_copy(pwqp->p[i].set),
1404 isl_qpolynomial_copy(pwqp->p[i].qp));
1405
1406 return dup;
1407 error:
1408 isl_pw_qpolynomial_free(dup);
1409 return NULL;
1410 }
1411
1412 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_cow(
1413 __isl_take isl_pw_qpolynomial *pwqp)
1414 {
1415 if (!pwqp)
1416 return NULL;
1417
1418 if (pwqp->ref == 1)
1419 return pwqp;
1420 pwqp->ref--;
1421 return isl_pw_qpolynomial_dup(pwqp);
1422 }
1423
1424 void isl_pw_qpolynomial_free(__isl_take isl_pw_qpolynomial *pwqp)
1425 {
1426 int i;
1427
1428 if (!pwqp)
1429 return;
1430 if (--pwqp->ref > 0)
1431 return;
1432
1433 for (i = 0; i < pwqp->n; ++i) {
1434 isl_set_free(pwqp->p[i].set);
1435 isl_qpolynomial_free(pwqp->p[i].qp);
1436 }
1437 isl_dim_free(pwqp->dim);
1438 free(pwqp);
1439 }
1440
1441 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
1442 __isl_take isl_pw_qpolynomial *pwqp1,
1443 __isl_take isl_pw_qpolynomial *pwqp2)
1444 {
1445 int i, j, n;
1446 struct isl_pw_qpolynomial *res;
1447 isl_set *set;
1448
1449 if (!pwqp1 || !pwqp2)
1450 goto error;
1451
1452 isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1453 goto error);
1454
1455 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1456 isl_pw_qpolynomial_free(pwqp1);
1457 return pwqp2;
1458 }
1459
1460 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1461 isl_pw_qpolynomial_free(pwqp2);
1462 return pwqp1;
1463 }
1464
1465 n = (pwqp1->n + 1) * (pwqp2->n + 1);
1466 res = pw_qpolynomial_alloc(isl_dim_copy(pwqp1->dim), n);
1467
1468 for (i = 0; i < pwqp1->n; ++i) {
1469 set = isl_set_copy(pwqp1->p[i].set);
1470 for (j = 0; j < pwqp2->n; ++j) {
1471 struct isl_set *common;
1472 struct isl_qpolynomial *sum;
1473 set = isl_set_subtract(set,
1474 isl_set_copy(pwqp2->p[j].set));
1475 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
1476 isl_set_copy(pwqp2->p[j].set));
1477 if (isl_set_fast_is_empty(common)) {
1478 isl_set_free(common);
1479 continue;
1480 }
1481
1482 sum = isl_qpolynomial_add(
1483 isl_qpolynomial_copy(pwqp1->p[i].qp),
1484 isl_qpolynomial_copy(pwqp2->p[j].qp));
1485
1486 res = isl_pw_qpolynomial_add_piece(res, common, sum);
1487 }
1488 res = isl_pw_qpolynomial_add_piece(res, set,
1489 isl_qpolynomial_copy(pwqp1->p[i].qp));
1490 }
1491
1492 for (j = 0; j < pwqp2->n; ++j) {
1493 set = isl_set_copy(pwqp2->p[j].set);
1494 for (i = 0; i < pwqp1->n; ++i)
1495 set = isl_set_subtract(set,
1496 isl_set_copy(pwqp1->p[i].set));
1497 res = isl_pw_qpolynomial_add_piece(res, set,
1498 isl_qpolynomial_copy(pwqp2->p[j].qp));
1499 }
1500
1501 isl_pw_qpolynomial_free(pwqp1);
1502 isl_pw_qpolynomial_free(pwqp2);
1503
1504 return res;
1505 error:
1506 isl_pw_qpolynomial_free(pwqp1);
1507 isl_pw_qpolynomial_free(pwqp2);
1508 return NULL;
1509 }
1510
1511 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_disjoint(
1512 __isl_take isl_pw_qpolynomial *pwqp1,
1513 __isl_take isl_pw_qpolynomial *pwqp2)
1514 {
1515 int i;
1516 struct isl_pw_qpolynomial *res;
1517
1518 if (!pwqp1 || !pwqp2)
1519 goto error;
1520
1521 isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1522 goto error);
1523
1524 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1525 isl_pw_qpolynomial_free(pwqp1);
1526 return pwqp2;
1527 }
1528
1529 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1530 isl_pw_qpolynomial_free(pwqp2);
1531 return pwqp1;
1532 }
1533
1534 res = pw_qpolynomial_alloc(isl_dim_copy(pwqp1->dim), pwqp1->n + pwqp2->n);
1535
1536 for (i = 0; i < pwqp1->n; ++i)
1537 res = isl_pw_qpolynomial_add_piece(res,
1538 isl_set_copy(pwqp1->p[i].set),
1539 isl_qpolynomial_copy(pwqp1->p[i].qp));
1540
1541 for (i = 0; i < pwqp2->n; ++i)
1542 res = isl_pw_qpolynomial_add_piece(res,
1543 isl_set_copy(pwqp2->p[i].set),
1544 isl_qpolynomial_copy(pwqp2->p[i].qp));
1545
1546 isl_pw_qpolynomial_free(pwqp1);
1547 isl_pw_qpolynomial_free(pwqp2);
1548
1549 return res;
1550 error:
1551 isl_pw_qpolynomial_free(pwqp1);
1552 isl_pw_qpolynomial_free(pwqp2);
1553 return NULL;
1554 }
1555
1556 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
1557 __isl_take isl_pw_qpolynomial *pwqp1,
1558 __isl_take isl_pw_qpolynomial *pwqp2)
1559 {
1560 int i, j, n;
1561 struct isl_pw_qpolynomial *res;
1562 isl_set *set;
1563
1564 if (!pwqp1 || !pwqp2)
1565 goto error;
1566
1567 isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1568 goto error);
1569
1570 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
1571 isl_pw_qpolynomial_free(pwqp2);
1572 return pwqp1;
1573 }
1574
1575 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
1576 isl_pw_qpolynomial_free(pwqp1);
1577 return pwqp2;
1578 }
1579
1580 if (isl_pw_qpolynomial_is_one(pwqp1)) {
1581 isl_pw_qpolynomial_free(pwqp1);
1582 return pwqp2;
1583 }
1584
1585 if (isl_pw_qpolynomial_is_one(pwqp2)) {
1586 isl_pw_qpolynomial_free(pwqp2);
1587 return pwqp1;
1588 }
1589
1590 n = pwqp1->n * pwqp2->n;
1591 res = pw_qpolynomial_alloc(isl_dim_copy(pwqp1->dim), n);
1592
1593 for (i = 0; i < pwqp1->n; ++i) {
1594 for (j = 0; j < pwqp2->n; ++j) {
1595 struct isl_set *common;
1596 struct isl_qpolynomial *prod;
1597 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
1598 isl_set_copy(pwqp2->p[j].set));
1599 if (isl_set_fast_is_empty(common)) {
1600 isl_set_free(common);
1601 continue;
1602 }
1603
1604 prod = isl_qpolynomial_mul(
1605 isl_qpolynomial_copy(pwqp1->p[i].qp),
1606 isl_qpolynomial_copy(pwqp2->p[j].qp));
1607
1608 res = isl_pw_qpolynomial_add_piece(res, common, prod);
1609 }
1610 }
1611
1612 isl_pw_qpolynomial_free(pwqp1);
1613 isl_pw_qpolynomial_free(pwqp2);
1614
1615 return res;
1616 error:
1617 isl_pw_qpolynomial_free(pwqp1);
1618 isl_pw_qpolynomial_free(pwqp2);
1619 return NULL;
1620 }
1621
1622 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_neg(
1623 __isl_take isl_pw_qpolynomial *pwqp)
1624 {
1625 int i, j, n;
1626 struct isl_pw_qpolynomial *res;
1627 isl_set *set;
1628
1629 if (!pwqp)
1630 return NULL;
1631
1632 if (isl_pw_qpolynomial_is_zero(pwqp))
1633 return pwqp;
1634
1635 pwqp = isl_pw_qpolynomial_cow(pwqp);
1636
1637 for (i = 0; i < pwqp->n; ++i) {
1638 pwqp->p[i].qp = isl_qpolynomial_neg(pwqp->p[i].qp);
1639 if (!pwqp->p[i].qp)
1640 goto error;
1641 }
1642
1643 return pwqp;
1644 error:
1645 isl_pw_qpolynomial_free(pwqp);
1646 return NULL;
1647 }
1648
1649 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sub(
1650 __isl_take isl_pw_qpolynomial *pwqp1,
1651 __isl_take isl_pw_qpolynomial *pwqp2)
1652 {
1653 return isl_pw_qpolynomial_add(pwqp1, isl_pw_qpolynomial_neg(pwqp2));
1654 }
1655
1656 __isl_give struct isl_upoly *isl_upoly_eval(
1657 __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
1658 {
1659 int i;
1660 struct isl_upoly_rec *rec;
1661 struct isl_upoly *res;
1662 struct isl_upoly *base;
1663
1664 if (isl_upoly_is_cst(up)) {
1665 isl_vec_free(vec);
1666 return up;
1667 }
1668
1669 rec = isl_upoly_as_rec(up);
1670 if (!rec)
1671 goto error;
1672
1673 isl_assert(up->ctx, rec->n >= 1, goto error);
1674
1675 base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
1676
1677 res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
1678 isl_vec_copy(vec));
1679
1680 for (i = rec->n - 2; i >= 0; --i) {
1681 res = isl_upoly_mul(res, isl_upoly_copy(base));
1682 res = isl_upoly_sum(res,
1683 isl_upoly_eval(isl_upoly_copy(rec->p[i]),
1684 isl_vec_copy(vec)));
1685 }
1686
1687 isl_upoly_free(base);
1688 isl_upoly_free(up);
1689 isl_vec_free(vec);
1690 return res;
1691 error:
1692 isl_upoly_free(up);
1693 isl_vec_free(vec);
1694 return NULL;
1695 }
1696
1697 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
1698 __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
1699 {
1700 isl_vec *ext;
1701 struct isl_upoly *up;
1702 isl_dim *dim;
1703
1704 if (!qp || !pnt)
1705 goto error;
1706 isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, qp->dim), goto error);
1707
1708 if (qp->div->n_row == 0)
1709 ext = isl_vec_copy(pnt->vec);
1710 else {
1711 int i;
1712 unsigned dim = isl_dim_total(qp->dim);
1713 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
1714 if (!ext)
1715 goto error;
1716
1717 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
1718 for (i = 0; i < qp->div->n_row; ++i) {
1719 isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
1720 1 + dim + i, &ext->el[1+dim+i]);
1721 isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
1722 qp->div->row[i][0]);
1723 }
1724 }
1725
1726 up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
1727 if (!up)
1728 goto error;
1729
1730 dim = isl_dim_copy(qp->dim);
1731 isl_qpolynomial_free(qp);
1732 isl_point_free(pnt);
1733
1734 qp = isl_qpolynomial_alloc(dim, 0);
1735 if (!qp)
1736 isl_upoly_free(up);
1737 else
1738 qp->upoly = up;
1739
1740 return qp;
1741 error:
1742 isl_qpolynomial_free(qp);
1743 isl_point_free(pnt);
1744 return NULL;
1745 }
1746
1747 __isl_give isl_qpolynomial *isl_pw_qpolynomial_eval(
1748 __isl_take isl_pw_qpolynomial *pwqp, __isl_take isl_point *pnt)
1749 {
1750 int i;
1751 int found;
1752 isl_qpolynomial *qp;
1753
1754 if (!pwqp || !pnt)
1755 goto error;
1756 isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, pwqp->dim), goto error);
1757
1758 for (i = 0; i < pwqp->n; ++i) {
1759 found = isl_set_contains_point(pwqp->p[i].set, pnt);
1760 if (found < 0)
1761 goto error;
1762 if (found)
1763 break;
1764 }
1765 if (found)
1766 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(pwqp->p[i].qp),
1767 isl_point_copy(pnt));
1768 else
1769 qp = isl_qpolynomial_zero(isl_dim_copy(pwqp->dim));
1770 isl_pw_qpolynomial_free(pwqp);
1771 isl_point_free(pnt);
1772 return qp;
1773 error:
1774 isl_pw_qpolynomial_free(pwqp);
1775 isl_point_free(pnt);
1776 return NULL;
1777 }
1778
1779 __isl_give isl_qpolynomial *isl_qpolynomial_move(__isl_take isl_qpolynomial *qp,
1780 enum isl_dim_type dst_type, unsigned dst_pos,
1781 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1782 {
1783 if (!qp)
1784 return NULL;
1785
1786 qp->dim = isl_dim_move(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
1787 if (!qp->dim)
1788 goto error;
1789
1790 /* Do something to polynomials when needed; later */
1791
1792 return qp;
1793 error:
1794 isl_qpolynomial_free(qp);
1795 return NULL;
1796 }
1797
1798 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_move(
1799 __isl_take isl_pw_qpolynomial *pwqp,
1800 enum isl_dim_type dst_type, unsigned dst_pos,
1801 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
1802 {
1803 int i;
1804
1805 if (!pwqp)
1806 return NULL;
1807
1808 pwqp->dim = isl_dim_move(pwqp->dim,
1809 dst_type, dst_pos, src_type, src_pos, n);
1810 if (!pwqp->dim)
1811 goto error;
1812
1813 for (i = 0; i < pwqp->n; ++i) {
1814 pwqp->p[i].set = isl_set_move(pwqp->p[i].set, dst_type, dst_pos,
1815 src_type, src_pos, n);
1816 if (!pwqp->p[i].set)
1817 goto error;
1818 pwqp->p[i].qp = isl_qpolynomial_move(pwqp->p[i].qp,
1819 dst_type, dst_pos, src_type, src_pos, n);
1820 if (!pwqp->p[i].qp)
1821 goto error;
1822 }
1823
1824 return pwqp;
1825 error:
1826 isl_pw_qpolynomial_free(pwqp);
1827 return NULL;
1828 }
1829
1830 __isl_give isl_dim *isl_pw_qpolynomial_get_dim(
1831 __isl_keep isl_pw_qpolynomial *pwqp)
1832 {
1833 if (!pwqp)
1834 return NULL;
1835
1836 return isl_dim_copy(pwqp->dim);
1837 }
1838
1839 unsigned isl_pw_qpolynomial_dim(__isl_keep isl_pw_qpolynomial *pwqp,
1840 enum isl_dim_type type)
1841 {
1842 return pwqp ? isl_dim_size(pwqp->dim, type) : 0;
1843 }
1844
1845 __isl_give isl_term *isl_term_alloc(__isl_take isl_dim *dim,
1846 __isl_take isl_mat *div)
1847 {
1848 isl_term *term;
1849 int n;
1850
1851 if (!dim || !div)
1852 goto error;
1853
1854 n = isl_dim_total(dim) + div->n_row;
1855
1856 term = isl_calloc(dim->ctx, struct isl_term,
1857 sizeof(struct isl_term) + (n - 1) * sizeof(int));
1858 if (!term)
1859 goto error;
1860
1861 term->ref = 1;
1862 term->dim = dim;
1863 term->div = div;
1864 isl_int_init(term->n);
1865 isl_int_init(term->d);
1866
1867 return term;
1868 error:
1869 isl_dim_free(dim);
1870 isl_mat_free(div);
1871 return NULL;
1872 }
1873
1874 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
1875 {
1876 if (!term)
1877 return NULL;
1878
1879 term->ref++;
1880 return term;
1881 }
1882
1883 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
1884 {
1885 int i;
1886 isl_term *dup;
1887 unsigned total;
1888
1889 if (term)
1890 return NULL;
1891
1892 total = isl_dim_total(term->dim) + term->div->n_row;
1893
1894 dup = isl_term_alloc(isl_dim_copy(term->dim), isl_mat_copy(term->div));
1895 if (!dup)
1896 return NULL;
1897
1898 isl_int_set(dup->n, term->n);
1899 isl_int_set(dup->d, term->d);
1900
1901 for (i = 0; i < total; ++i)
1902 dup->pow[i] = term->pow[i];
1903
1904 return dup;
1905 }
1906
1907 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
1908 {
1909 if (!term)
1910 return NULL;
1911
1912 if (term->ref == 1)
1913 return term;
1914 term->ref--;
1915 return isl_term_dup(term);
1916 }
1917
1918 void isl_term_free(__isl_take isl_term *term)
1919 {
1920 if (!term)
1921 return;
1922
1923 if (--term->ref > 0)
1924 return;
1925
1926 isl_dim_free(term->dim);
1927 isl_mat_free(term->div);
1928 isl_int_clear(term->n);
1929 isl_int_clear(term->d);
1930 free(term);
1931 }
1932
1933 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
1934 {
1935 if (!term)
1936 return 0;
1937
1938 switch (type) {
1939 case isl_dim_param:
1940 case isl_dim_in:
1941 case isl_dim_out: return isl_dim_size(term->dim, type);
1942 case isl_dim_div: return term->div->n_row;
1943 case isl_dim_all: return isl_dim_total(term->dim) + term->div->n_row;
1944 }
1945 }
1946
1947 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
1948 {
1949 return term ? term->dim->ctx : NULL;
1950 }
1951
1952 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
1953 {
1954 if (!term)
1955 return;
1956 isl_int_set(*n, term->n);
1957 }
1958
1959 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
1960 {
1961 if (!term)
1962 return;
1963 isl_int_set(*d, term->d);
1964 }
1965
1966 int isl_term_get_exp(__isl_keep isl_term *term,
1967 enum isl_dim_type type, unsigned pos)
1968 {
1969 if (!term)
1970 return -1;
1971
1972 isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
1973
1974 if (type >= isl_dim_set)
1975 pos += isl_dim_size(term->dim, isl_dim_param);
1976 if (type >= isl_dim_div)
1977 pos += isl_dim_size(term->dim, isl_dim_set);
1978
1979 return term->pow[pos];
1980 }
1981
1982 __isl_give isl_div *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
1983 {
1984 isl_basic_map *bmap;
1985 unsigned total;
1986 int k;
1987
1988 if (!term)
1989 return NULL;
1990
1991 isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
1992 return NULL);
1993
1994 total = term->div->n_col - term->div->n_row - 2;
1995 /* No nested divs for now */
1996 isl_assert(term->dim->ctx,
1997 isl_seq_first_non_zero(term->div->row[pos] + 2 + total,
1998 term->div->n_row) == -1,
1999 return NULL);
2000
2001 bmap = isl_basic_map_alloc_dim(isl_dim_copy(term->dim), 1, 0, 0);
2002 if ((k = isl_basic_map_alloc_div(bmap)) < 0)
2003 goto error;
2004
2005 isl_seq_cpy(bmap->div[k], term->div->row[pos], 2 + total);
2006
2007 return isl_basic_map_div(bmap, k);
2008 error:
2009 isl_basic_map_free(bmap);
2010 return NULL;
2011 }
2012
2013 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
2014 int (*fn)(__isl_take isl_term *term, void *user),
2015 __isl_take isl_term *term, void *user)
2016 {
2017 int i;
2018 struct isl_upoly_rec *rec;
2019
2020 if (!up || !term)
2021 goto error;
2022
2023 if (isl_upoly_is_zero(up))
2024 return term;
2025
2026 isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
2027 isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
2028 isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
2029
2030 if (isl_upoly_is_cst(up)) {
2031 struct isl_upoly_cst *cst;
2032 cst = isl_upoly_as_cst(up);
2033 if (!cst)
2034 goto error;
2035 term = isl_term_cow(term);
2036 if (!term)
2037 goto error;
2038 isl_int_set(term->n, cst->n);
2039 isl_int_set(term->d, cst->d);
2040 if (fn(isl_term_copy(term), user) < 0)
2041 goto error;
2042 return term;
2043 }
2044
2045 rec = isl_upoly_as_rec(up);
2046 if (!rec)
2047 goto error;
2048
2049 for (i = 0; i < rec->n; ++i) {
2050 term = isl_term_cow(term);
2051 if (!term)
2052 goto error;
2053 term->pow[up->var] = i;
2054 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
2055 if (!term)
2056 goto error;
2057 }
2058 term->pow[up->var] = 0;
2059
2060 return term;
2061 error:
2062 isl_term_free(term);
2063 return NULL;
2064 }
2065
2066 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
2067 int (*fn)(__isl_take isl_term *term, void *user), void *user)
2068 {
2069 isl_term *term;
2070
2071 if (!qp)
2072 return -1;
2073
2074 term = isl_term_alloc(isl_dim_copy(qp->dim), isl_mat_copy(qp->div));
2075 if (!term)
2076 return -1;
2077
2078 term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
2079
2080 isl_term_free(term);
2081
2082 return term ? 0 : -1;
2083 }
2084
2085 int isl_pw_qpolynomial_foreach_piece(__isl_keep isl_pw_qpolynomial *pwqp,
2086 int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2087 void *user), void *user)
2088 {
2089 int i;
2090
2091 if (!pwqp)
2092 return -1;
2093
2094 for (i = 0; i < pwqp->n; ++i)
2095 if (fn(isl_set_copy(pwqp->p[i].set),
2096 isl_qpolynomial_copy(pwqp->p[i].qp), user) < 0)
2097 return -1;
2098
2099 return 0;
2100 }
2101
2102 static int any_divs(__isl_keep isl_set *set)
2103 {
2104 int i;
2105
2106 if (!set)
2107 return -1;
2108
2109 for (i = 0; i < set->n; ++i)
2110 if (set->p[i]->n_div > 0)
2111 return 1;
2112
2113 return 0;
2114 }
2115
2116 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
2117 __isl_take isl_dim *dim)
2118 {
2119 if (!qp || !dim)
2120 goto error;
2121
2122 if (isl_dim_equal(qp->dim, dim)) {
2123 isl_dim_free(dim);
2124 return qp;
2125 }
2126
2127 qp = isl_qpolynomial_cow(qp);
2128 if (!qp)
2129 goto error;
2130
2131 if (qp->div->n_row) {
2132 int i;
2133 int extra;
2134 unsigned total;
2135 int *exp;
2136
2137 extra = isl_dim_size(dim, isl_dim_set) -
2138 isl_dim_size(qp->dim, isl_dim_set);
2139 total = isl_dim_total(qp->dim);
2140 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
2141 if (!exp)
2142 goto error;
2143 for (i = 0; i < qp->div->n_row; ++i)
2144 exp[i] = extra + i;
2145 qp->upoly = expand(qp->upoly, exp, total);
2146 free(exp);
2147 if (!qp->upoly)
2148 goto error;
2149 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
2150 if (!qp->div)
2151 goto error;
2152 for (i = 0; i < qp->div->n_row; ++i)
2153 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
2154 }
2155
2156 isl_dim_free(qp->dim);
2157 qp->dim = dim;
2158
2159 return qp;
2160 error:
2161 isl_dim_free(dim);
2162 isl_qpolynomial_free(qp);
2163 return NULL;
2164 }
2165
2166 static int foreach_lifted_subset(__isl_take isl_set *set,
2167 __isl_take isl_qpolynomial *qp,
2168 int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2169 void *user), void *user)
2170 {
2171 int i;
2172
2173 if (!set || !qp)
2174 goto error;
2175
2176 for (i = 0; i < set->n; ++i) {
2177 isl_set *lift;
2178 isl_qpolynomial *copy;
2179
2180 lift = isl_set_from_basic_set(isl_basic_set_copy(set->p[i]));
2181 lift = isl_set_lift(lift);
2182
2183 copy = isl_qpolynomial_copy(qp);
2184 copy = isl_qpolynomial_lift(copy, isl_set_get_dim(lift));
2185
2186 if (fn(lift, copy, user) < 0)
2187 goto error;
2188 }
2189
2190 isl_set_free(set);
2191 isl_qpolynomial_free(qp);
2192
2193 return 0;
2194 error:
2195 isl_set_free(set);
2196 isl_qpolynomial_free(qp);
2197 return -1;
2198 }
2199
2200 int isl_pw_qpolynomial_foreach_lifted_piece(__isl_keep isl_pw_qpolynomial *pwqp,
2201 int (*fn)(__isl_take isl_set *set, __isl_take isl_qpolynomial *qp,
2202 void *user), void *user)
2203 {
2204 int i;
2205
2206 if (!pwqp)
2207 return -1;
2208
2209 for (i = 0; i < pwqp->n; ++i) {
2210 isl_set *set;
2211 isl_qpolynomial *qp;
2212
2213 set = isl_set_copy(pwqp->p[i].set);
2214 qp = isl_qpolynomial_copy(pwqp->p[i].qp);
2215 if (!any_divs(set)) {
2216 if (fn(set, qp, user) < 0)
2217 return -1;
2218 continue;
2219 }
2220 if (foreach_lifted_subset(set, qp, fn, user) < 0)
2221 return -1;
2222 }
2223
2224 return 0;
2225 }
Integer set library