?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
isl_map_read_from_str: parse tuple entries with leading minus sign
[isl.git] / isl_polynomial.c
blobdf1898c48f8868e59e3b3f79b316c4b4040f007f
1 /*
2 * Copyright 2010 INRIA Saclay
3 *
4 * Use of this software is governed by the GNU LGPLv2.1 license
5 *
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
8 * 91893 Orsay, France
9 */
10
11 #include <stdlib.h>
12 #include <isl_seq.h>
13 #include <isl_union_map_private.h>
14 #include <isl_polynomial_private.h>
15 #include <isl_point_private.h>
16 #include <isl_dim_private.h>
17 #include <isl_map_private.h>
18
19 static unsigned pos(__isl_keep isl_dim *dim, enum isl_dim_type type)
20 {
21 switch (type) {
22 case isl_dim_param: return 0;
23 case isl_dim_in: return dim->nparam;
24 case isl_dim_out: return dim->nparam + dim->n_in;
25 default: return 0;
26 }
27 }
28
29 int isl_upoly_is_cst(__isl_keep struct isl_upoly *up)
30 {
31 if (!up)
32 return -1;
33
34 return up->var < 0;
35 }
36
37 __isl_keep struct isl_upoly_cst *isl_upoly_as_cst(__isl_keep struct isl_upoly *up)
38 {
39 if (!up)
40 return NULL;
41
42 isl_assert(up->ctx, up->var < 0, return NULL);
43
44 return (struct isl_upoly_cst *)up;
45 }
46
47 __isl_keep struct isl_upoly_rec *isl_upoly_as_rec(__isl_keep struct isl_upoly *up)
48 {
49 if (!up)
50 return NULL;
51
52 isl_assert(up->ctx, up->var >= 0, return NULL);
53
54 return (struct isl_upoly_rec *)up;
55 }
56
57 int isl_upoly_is_equal(__isl_keep struct isl_upoly *up1,
58 __isl_keep struct isl_upoly *up2)
59 {
60 int i;
61 struct isl_upoly_rec *rec1, *rec2;
62
63 if (!up1 || !up2)
64 return -1;
65 if (up1 == up2)
66 return 1;
67 if (up1->var != up2->var)
68 return 0;
69 if (isl_upoly_is_cst(up1)) {
70 struct isl_upoly_cst *cst1, *cst2;
71 cst1 = isl_upoly_as_cst(up1);
72 cst2 = isl_upoly_as_cst(up2);
73 if (!cst1 || !cst2)
74 return -1;
75 return isl_int_eq(cst1->n, cst2->n) &&
76 isl_int_eq(cst1->d, cst2->d);
77 }
78
79 rec1 = isl_upoly_as_rec(up1);
80 rec2 = isl_upoly_as_rec(up2);
81 if (!rec1 || !rec2)
82 return -1;
83
84 if (rec1->n != rec2->n)
85 return 0;
86
87 for (i = 0; i < rec1->n; ++i) {
88 int eq = isl_upoly_is_equal(rec1->p[i], rec2->p[i]);
89 if (eq < 0 || !eq)
90 return eq;
91 }
92
93 return 1;
94 }
95
96 int isl_upoly_is_zero(__isl_keep struct isl_upoly *up)
97 {
98 struct isl_upoly_cst *cst;
99
100 if (!up)
101 return -1;
102 if (!isl_upoly_is_cst(up))
103 return 0;
104
105 cst = isl_upoly_as_cst(up);
106 if (!cst)
107 return -1;
108
109 return isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d);
110 }
111
112 int isl_upoly_sgn(__isl_keep struct isl_upoly *up)
113 {
114 struct isl_upoly_cst *cst;
115
116 if (!up)
117 return 0;
118 if (!isl_upoly_is_cst(up))
119 return 0;
120
121 cst = isl_upoly_as_cst(up);
122 if (!cst)
123 return 0;
124
125 return isl_int_sgn(cst->n);
126 }
127
128 int isl_upoly_is_nan(__isl_keep struct isl_upoly *up)
129 {
130 struct isl_upoly_cst *cst;
131
132 if (!up)
133 return -1;
134 if (!isl_upoly_is_cst(up))
135 return 0;
136
137 cst = isl_upoly_as_cst(up);
138 if (!cst)
139 return -1;
140
141 return isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d);
142 }
143
144 int isl_upoly_is_infty(__isl_keep struct isl_upoly *up)
145 {
146 struct isl_upoly_cst *cst;
147
148 if (!up)
149 return -1;
150 if (!isl_upoly_is_cst(up))
151 return 0;
152
153 cst = isl_upoly_as_cst(up);
154 if (!cst)
155 return -1;
156
157 return isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d);
158 }
159
160 int isl_upoly_is_neginfty(__isl_keep struct isl_upoly *up)
161 {
162 struct isl_upoly_cst *cst;
163
164 if (!up)
165 return -1;
166 if (!isl_upoly_is_cst(up))
167 return 0;
168
169 cst = isl_upoly_as_cst(up);
170 if (!cst)
171 return -1;
172
173 return isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d);
174 }
175
176 int isl_upoly_is_one(__isl_keep struct isl_upoly *up)
177 {
178 struct isl_upoly_cst *cst;
179
180 if (!up)
181 return -1;
182 if (!isl_upoly_is_cst(up))
183 return 0;
184
185 cst = isl_upoly_as_cst(up);
186 if (!cst)
187 return -1;
188
189 return isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
190 }
191
192 int isl_upoly_is_negone(__isl_keep struct isl_upoly *up)
193 {
194 struct isl_upoly_cst *cst;
195
196 if (!up)
197 return -1;
198 if (!isl_upoly_is_cst(up))
199 return 0;
200
201 cst = isl_upoly_as_cst(up);
202 if (!cst)
203 return -1;
204
205 return isl_int_is_negone(cst->n) && isl_int_is_one(cst->d);
206 }
207
208 __isl_give struct isl_upoly_cst *isl_upoly_cst_alloc(struct isl_ctx *ctx)
209 {
210 struct isl_upoly_cst *cst;
211
212 cst = isl_alloc_type(ctx, struct isl_upoly_cst);
213 if (!cst)
214 return NULL;
215
216 cst->up.ref = 1;
217 cst->up.ctx = ctx;
218 isl_ctx_ref(ctx);
219 cst->up.var = -1;
220
221 isl_int_init(cst->n);
222 isl_int_init(cst->d);
223
224 return cst;
225 }
226
227 __isl_give struct isl_upoly *isl_upoly_zero(struct isl_ctx *ctx)
228 {
229 struct isl_upoly_cst *cst;
230
231 cst = isl_upoly_cst_alloc(ctx);
232 if (!cst)
233 return NULL;
234
235 isl_int_set_si(cst->n, 0);
236 isl_int_set_si(cst->d, 1);
237
238 return &cst->up;
239 }
240
241 __isl_give struct isl_upoly *isl_upoly_infty(struct isl_ctx *ctx)
242 {
243 struct isl_upoly_cst *cst;
244
245 cst = isl_upoly_cst_alloc(ctx);
246 if (!cst)
247 return NULL;
248
249 isl_int_set_si(cst->n, 1);
250 isl_int_set_si(cst->d, 0);
251
252 return &cst->up;
253 }
254
255 __isl_give struct isl_upoly *isl_upoly_neginfty(struct isl_ctx *ctx)
256 {
257 struct isl_upoly_cst *cst;
258
259 cst = isl_upoly_cst_alloc(ctx);
260 if (!cst)
261 return NULL;
262
263 isl_int_set_si(cst->n, -1);
264 isl_int_set_si(cst->d, 0);
265
266 return &cst->up;
267 }
268
269 __isl_give struct isl_upoly *isl_upoly_nan(struct isl_ctx *ctx)
270 {
271 struct isl_upoly_cst *cst;
272
273 cst = isl_upoly_cst_alloc(ctx);
274 if (!cst)
275 return NULL;
276
277 isl_int_set_si(cst->n, 0);
278 isl_int_set_si(cst->d, 0);
279
280 return &cst->up;
281 }
282
283 __isl_give struct isl_upoly *isl_upoly_rat_cst(struct isl_ctx *ctx,
284 isl_int n, isl_int d)
285 {
286 struct isl_upoly_cst *cst;
287
288 cst = isl_upoly_cst_alloc(ctx);
289 if (!cst)
290 return NULL;
291
292 isl_int_set(cst->n, n);
293 isl_int_set(cst->d, d);
294
295 return &cst->up;
296 }
297
298 __isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
299 int var, int size)
300 {
301 struct isl_upoly_rec *rec;
302
303 isl_assert(ctx, var >= 0, return NULL);
304 isl_assert(ctx, size >= 0, return NULL);
305 rec = isl_calloc(ctx, struct isl_upoly_rec,
306 sizeof(struct isl_upoly_rec) +
307 (size - 1) * sizeof(struct isl_upoly *));
308 if (!rec)
309 return NULL;
310
311 rec->up.ref = 1;
312 rec->up.ctx = ctx;
313 isl_ctx_ref(ctx);
314 rec->up.var = var;
315
316 rec->n = 0;
317 rec->size = size;
318
319 return rec;
320 }
321
322 __isl_give isl_qpolynomial *isl_qpolynomial_reset_dim(
323 __isl_take isl_qpolynomial *qp, __isl_take isl_dim *dim)
324 {
325 qp = isl_qpolynomial_cow(qp);
326 if (!qp || !dim)
327 goto error;
328
329 isl_dim_free(qp->dim);
330 qp->dim = dim;
331
332 return qp;
333 error:
334 isl_qpolynomial_free(qp);
335 isl_dim_free(dim);
336 return NULL;
337 }
338
339 isl_ctx *isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial *qp)
340 {
341 return qp ? qp->dim->ctx : NULL;
342 }
343
344 __isl_give isl_dim *isl_qpolynomial_get_dim(__isl_keep isl_qpolynomial *qp)
345 {
346 return qp ? isl_dim_copy(qp->dim) : NULL;
347 }
348
349 unsigned isl_qpolynomial_dim(__isl_keep isl_qpolynomial *qp,
350 enum isl_dim_type type)
351 {
352 return qp ? isl_dim_size(qp->dim, type) : 0;
353 }
354
355 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
356 {
357 return qp ? isl_upoly_is_zero(qp->upoly) : -1;
358 }
359
360 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
361 {
362 return qp ? isl_upoly_is_one(qp->upoly) : -1;
363 }
364
365 int isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial *qp)
366 {
367 return qp ? isl_upoly_is_nan(qp->upoly) : -1;
368 }
369
370 int isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial *qp)
371 {
372 return qp ? isl_upoly_is_infty(qp->upoly) : -1;
373 }
374
375 int isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial *qp)
376 {
377 return qp ? isl_upoly_is_neginfty(qp->upoly) : -1;
378 }
379
380 int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial *qp)
381 {
382 return qp ? isl_upoly_sgn(qp->upoly) : 0;
383 }
384
385 static void upoly_free_cst(__isl_take struct isl_upoly_cst *cst)
386 {
387 isl_int_clear(cst->n);
388 isl_int_clear(cst->d);
389 }
390
391 static void upoly_free_rec(__isl_take struct isl_upoly_rec *rec)
392 {
393 int i;
394
395 for (i = 0; i < rec->n; ++i)
396 isl_upoly_free(rec->p[i]);
397 }
398
399 __isl_give struct isl_upoly *isl_upoly_copy(__isl_keep struct isl_upoly *up)
400 {
401 if (!up)
402 return NULL;
403
404 up->ref++;
405 return up;
406 }
407
408 __isl_give struct isl_upoly *isl_upoly_dup_cst(__isl_keep struct isl_upoly *up)
409 {
410 struct isl_upoly_cst *cst;
411 struct isl_upoly_cst *dup;
412
413 cst = isl_upoly_as_cst(up);
414 if (!cst)
415 return NULL;
416
417 dup = isl_upoly_as_cst(isl_upoly_zero(up->ctx));
418 if (!dup)
419 return NULL;
420 isl_int_set(dup->n, cst->n);
421 isl_int_set(dup->d, cst->d);
422
423 return &dup->up;
424 }
425
426 __isl_give struct isl_upoly *isl_upoly_dup_rec(__isl_keep struct isl_upoly *up)
427 {
428 int i;
429 struct isl_upoly_rec *rec;
430 struct isl_upoly_rec *dup;
431
432 rec = isl_upoly_as_rec(up);
433 if (!rec)
434 return NULL;
435
436 dup = isl_upoly_alloc_rec(up->ctx, up->var, rec->n);
437 if (!dup)
438 return NULL;
439
440 for (i = 0; i < rec->n; ++i) {
441 dup->p[i] = isl_upoly_copy(rec->p[i]);
442 if (!dup->p[i])
443 goto error;
444 dup->n++;
445 }
446
447 return &dup->up;
448 error:
449 isl_upoly_free(&dup->up);
450 return NULL;
451 }
452
453 __isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
454 {
455 struct isl_upoly *dup;
456
457 if (!up)
458 return NULL;
459
460 if (isl_upoly_is_cst(up))
461 return isl_upoly_dup_cst(up);
462 else
463 return isl_upoly_dup_rec(up);
464 }
465
466 __isl_give struct isl_upoly *isl_upoly_cow(__isl_take struct isl_upoly *up)
467 {
468 if (!up)
469 return NULL;
470
471 if (up->ref == 1)
472 return up;
473 up->ref--;
474 return isl_upoly_dup(up);
475 }
476
477 void isl_upoly_free(__isl_take struct isl_upoly *up)
478 {
479 if (!up)
480 return;
481
482 if (--up->ref > 0)
483 return;
484
485 if (up->var < 0)
486 upoly_free_cst((struct isl_upoly_cst *)up);
487 else
488 upoly_free_rec((struct isl_upoly_rec *)up);
489
490 isl_ctx_deref(up->ctx);
491 free(up);
492 }
493
494 static void isl_upoly_cst_reduce(__isl_keep struct isl_upoly_cst *cst)
495 {
496 isl_int gcd;
497
498 isl_int_init(gcd);
499 isl_int_gcd(gcd, cst->n, cst->d);
500 if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
501 isl_int_divexact(cst->n, cst->n, gcd);
502 isl_int_divexact(cst->d, cst->d, gcd);
503 }
504 isl_int_clear(gcd);
505 }
506
507 __isl_give struct isl_upoly *isl_upoly_sum_cst(__isl_take struct isl_upoly *up1,
508 __isl_take struct isl_upoly *up2)
509 {
510 struct isl_upoly_cst *cst1;
511 struct isl_upoly_cst *cst2;
512
513 up1 = isl_upoly_cow(up1);
514 if (!up1 || !up2)
515 goto error;
516
517 cst1 = isl_upoly_as_cst(up1);
518 cst2 = isl_upoly_as_cst(up2);
519
520 if (isl_int_eq(cst1->d, cst2->d))
521 isl_int_add(cst1->n, cst1->n, cst2->n);
522 else {
523 isl_int_mul(cst1->n, cst1->n, cst2->d);
524 isl_int_addmul(cst1->n, cst2->n, cst1->d);
525 isl_int_mul(cst1->d, cst1->d, cst2->d);
526 }
527
528 isl_upoly_cst_reduce(cst1);
529
530 isl_upoly_free(up2);
531 return up1;
532 error:
533 isl_upoly_free(up1);
534 isl_upoly_free(up2);
535 return NULL;
536 }
537
538 static __isl_give struct isl_upoly *replace_by_zero(
539 __isl_take struct isl_upoly *up)
540 {
541 struct isl_ctx *ctx;
542
543 if (!up)
544 return NULL;
545 ctx = up->ctx;
546 isl_upoly_free(up);
547 return isl_upoly_zero(ctx);
548 }
549
550 static __isl_give struct isl_upoly *replace_by_constant_term(
551 __isl_take struct isl_upoly *up)
552 {
553 struct isl_upoly_rec *rec;
554 struct isl_upoly *cst;
555
556 if (!up)
557 return NULL;
558
559 rec = isl_upoly_as_rec(up);
560 if (!rec)
561 goto error;
562 cst = isl_upoly_copy(rec->p[0]);
563 isl_upoly_free(up);
564 return cst;
565 error:
566 isl_upoly_free(up);
567 return NULL;
568 }
569
570 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
571 __isl_take struct isl_upoly *up2)
572 {
573 int i;
574 struct isl_upoly_rec *rec1, *rec2;
575
576 if (!up1 || !up2)
577 goto error;
578
579 if (isl_upoly_is_nan(up1)) {
580 isl_upoly_free(up2);
581 return up1;
582 }
583
584 if (isl_upoly_is_nan(up2)) {
585 isl_upoly_free(up1);
586 return up2;
587 }
588
589 if (isl_upoly_is_zero(up1)) {
590 isl_upoly_free(up1);
591 return up2;
592 }
593
594 if (isl_upoly_is_zero(up2)) {
595 isl_upoly_free(up2);
596 return up1;
597 }
598
599 if (up1->var < up2->var)
600 return isl_upoly_sum(up2, up1);
601
602 if (up2->var < up1->var) {
603 struct isl_upoly_rec *rec;
604 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
605 isl_upoly_free(up1);
606 return up2;
607 }
608 up1 = isl_upoly_cow(up1);
609 rec = isl_upoly_as_rec(up1);
610 if (!rec)
611 goto error;
612 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
613 if (rec->n == 1)
614 up1 = replace_by_constant_term(up1);
615 return up1;
616 }
617
618 if (isl_upoly_is_cst(up1))
619 return isl_upoly_sum_cst(up1, up2);
620
621 rec1 = isl_upoly_as_rec(up1);
622 rec2 = isl_upoly_as_rec(up2);
623 if (!rec1 || !rec2)
624 goto error;
625
626 if (rec1->n < rec2->n)
627 return isl_upoly_sum(up2, up1);
628
629 up1 = isl_upoly_cow(up1);
630 rec1 = isl_upoly_as_rec(up1);
631 if (!rec1)
632 goto error;
633
634 for (i = rec2->n - 1; i >= 0; --i) {
635 rec1->p[i] = isl_upoly_sum(rec1->p[i],
636 isl_upoly_copy(rec2->p[i]));
637 if (!rec1->p[i])
638 goto error;
639 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
640 isl_upoly_free(rec1->p[i]);
641 rec1->n--;
642 }
643 }
644
645 if (rec1->n == 0)
646 up1 = replace_by_zero(up1);
647 else if (rec1->n == 1)
648 up1 = replace_by_constant_term(up1);
649
650 isl_upoly_free(up2);
651
652 return up1;
653 error:
654 isl_upoly_free(up1);
655 isl_upoly_free(up2);
656 return NULL;
657 }
658
659 __isl_give struct isl_upoly *isl_upoly_neg_cst(__isl_take struct isl_upoly *up)
660 {
661 struct isl_upoly_cst *cst;
662
663 if (isl_upoly_is_zero(up))
664 return up;
665
666 up = isl_upoly_cow(up);
667 if (!up)
668 return NULL;
669
670 cst = isl_upoly_as_cst(up);
671
672 isl_int_neg(cst->n, cst->n);
673
674 return up;
675 }
676
677 __isl_give struct isl_upoly *isl_upoly_neg(__isl_take struct isl_upoly *up)
678 {
679 int i;
680 struct isl_upoly_rec *rec;
681
682 if (!up)
683 return NULL;
684
685 if (isl_upoly_is_cst(up))
686 return isl_upoly_neg_cst(up);
687
688 up = isl_upoly_cow(up);
689 rec = isl_upoly_as_rec(up);
690 if (!rec)
691 goto error;
692
693 for (i = 0; i < rec->n; ++i) {
694 rec->p[i] = isl_upoly_neg(rec->p[i]);
695 if (!rec->p[i])
696 goto error;
697 }
698
699 return up;
700 error:
701 isl_upoly_free(up);
702 return NULL;
703 }
704
705 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
706 __isl_take struct isl_upoly *up2)
707 {
708 struct isl_upoly_cst *cst1;
709 struct isl_upoly_cst *cst2;
710
711 up1 = isl_upoly_cow(up1);
712 if (!up1 || !up2)
713 goto error;
714
715 cst1 = isl_upoly_as_cst(up1);
716 cst2 = isl_upoly_as_cst(up2);
717
718 isl_int_mul(cst1->n, cst1->n, cst2->n);
719 isl_int_mul(cst1->d, cst1->d, cst2->d);
720
721 isl_upoly_cst_reduce(cst1);
722
723 isl_upoly_free(up2);
724 return up1;
725 error:
726 isl_upoly_free(up1);
727 isl_upoly_free(up2);
728 return NULL;
729 }
730
731 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
732 __isl_take struct isl_upoly *up2)
733 {
734 struct isl_upoly_rec *rec1;
735 struct isl_upoly_rec *rec2;
736 struct isl_upoly_rec *res;
737 int i, j;
738 int size;
739
740 rec1 = isl_upoly_as_rec(up1);
741 rec2 = isl_upoly_as_rec(up2);
742 if (!rec1 || !rec2)
743 goto error;
744 size = rec1->n + rec2->n - 1;
745 res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
746 if (!res)
747 goto error;
748
749 for (i = 0; i < rec1->n; ++i) {
750 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
751 isl_upoly_copy(rec1->p[i]));
752 if (!res->p[i])
753 goto error;
754 res->n++;
755 }
756 for (; i < size; ++i) {
757 res->p[i] = isl_upoly_zero(up1->ctx);
758 if (!res->p[i])
759 goto error;
760 res->n++;
761 }
762 for (i = 0; i < rec1->n; ++i) {
763 for (j = 1; j < rec2->n; ++j) {
764 struct isl_upoly *up;
765 up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
766 isl_upoly_copy(rec1->p[i]));
767 res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
768 if (!res->p[i + j])
769 goto error;
770 }
771 }
772
773 isl_upoly_free(up1);
774 isl_upoly_free(up2);
775
776 return &res->up;
777 error:
778 isl_upoly_free(up1);
779 isl_upoly_free(up2);
780 isl_upoly_free(&res->up);
781 return NULL;
782 }
783
784 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
785 __isl_take struct isl_upoly *up2)
786 {
787 if (!up1 || !up2)
788 goto error;
789
790 if (isl_upoly_is_nan(up1)) {
791 isl_upoly_free(up2);
792 return up1;
793 }
794
795 if (isl_upoly_is_nan(up2)) {
796 isl_upoly_free(up1);
797 return up2;
798 }
799
800 if (isl_upoly_is_zero(up1)) {
801 isl_upoly_free(up2);
802 return up1;
803 }
804
805 if (isl_upoly_is_zero(up2)) {
806 isl_upoly_free(up1);
807 return up2;
808 }
809
810 if (isl_upoly_is_one(up1)) {
811 isl_upoly_free(up1);
812 return up2;
813 }
814
815 if (isl_upoly_is_one(up2)) {
816 isl_upoly_free(up2);
817 return up1;
818 }
819
820 if (up1->var < up2->var)
821 return isl_upoly_mul(up2, up1);
822
823 if (up2->var < up1->var) {
824 int i;
825 struct isl_upoly_rec *rec;
826 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
827 isl_ctx *ctx = up1->ctx;
828 isl_upoly_free(up1);
829 isl_upoly_free(up2);
830 return isl_upoly_nan(ctx);
831 }
832 up1 = isl_upoly_cow(up1);
833 rec = isl_upoly_as_rec(up1);
834 if (!rec)
835 goto error;
836
837 for (i = 0; i < rec->n; ++i) {
838 rec->p[i] = isl_upoly_mul(rec->p[i],
839 isl_upoly_copy(up2));
840 if (!rec->p[i])
841 goto error;
842 }
843 isl_upoly_free(up2);
844 return up1;
845 }
846
847 if (isl_upoly_is_cst(up1))
848 return isl_upoly_mul_cst(up1, up2);
849
850 return isl_upoly_mul_rec(up1, up2);
851 error:
852 isl_upoly_free(up1);
853 isl_upoly_free(up2);
854 return NULL;
855 }
856
857 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_dim *dim,
858 unsigned n_div, __isl_take struct isl_upoly *up)
859 {
860 struct isl_qpolynomial *qp = NULL;
861 unsigned total;
862
863 if (!dim || !up)
864 goto error;
865
866 total = isl_dim_total(dim);
867
868 qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
869 if (!qp)
870 goto error;
871
872 qp->ref = 1;
873 qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
874 if (!qp->div)
875 goto error;
876
877 qp->dim = dim;
878 qp->upoly = up;
879
880 return qp;
881 error:
882 isl_dim_free(dim);
883 isl_upoly_free(up);
884 isl_qpolynomial_free(qp);
885 return NULL;
886 }
887
888 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
889 {
890 if (!qp)
891 return NULL;
892
893 qp->ref++;
894 return qp;
895 }
896
897 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
898 {
899 struct isl_qpolynomial *dup;
900
901 if (!qp)
902 return NULL;
903
904 dup = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), qp->div->n_row,
905 isl_upoly_copy(qp->upoly));
906 if (!dup)
907 return NULL;
908 isl_mat_free(dup->div);
909 dup->div = isl_mat_copy(qp->div);
910 if (!dup->div)
911 goto error;
912
913 return dup;
914 error:
915 isl_qpolynomial_free(dup);
916 return NULL;
917 }
918
919 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
920 {
921 if (!qp)
922 return NULL;
923
924 if (qp->ref == 1)
925 return qp;
926 qp->ref--;
927 return isl_qpolynomial_dup(qp);
928 }
929
930 void isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
931 {
932 if (!qp)
933 return;
934
935 if (--qp->ref > 0)
936 return;
937
938 isl_dim_free(qp->dim);
939 isl_mat_free(qp->div);
940 isl_upoly_free(qp->upoly);
941
942 free(qp);
943 }
944
945 __isl_give struct isl_upoly *isl_upoly_pow(isl_ctx *ctx, int pos, int power)
946 {
947 int i;
948 struct isl_upoly *up;
949 struct isl_upoly_rec *rec;
950 struct isl_upoly_cst *cst;
951
952 rec = isl_upoly_alloc_rec(ctx, pos, 1 + power);
953 if (!rec)
954 return NULL;
955 for (i = 0; i < 1 + power; ++i) {
956 rec->p[i] = isl_upoly_zero(ctx);
957 if (!rec->p[i])
958 goto error;
959 rec->n++;
960 }
961 cst = isl_upoly_as_cst(rec->p[power]);
962 isl_int_set_si(cst->n, 1);
963
964 return &rec->up;
965 error:
966 isl_upoly_free(&rec->up);
967 return NULL;
968 }
969
970 /* r array maps original positions to new positions.
971 */
972 static __isl_give struct isl_upoly *reorder(__isl_take struct isl_upoly *up,
973 int *r)
974 {
975 int i;
976 struct isl_upoly_rec *rec;
977 struct isl_upoly *base;
978 struct isl_upoly *res;
979
980 if (isl_upoly_is_cst(up))
981 return up;
982
983 rec = isl_upoly_as_rec(up);
984 if (!rec)
985 goto error;
986
987 isl_assert(up->ctx, rec->n >= 1, goto error);
988
989 base = isl_upoly_pow(up->ctx, r[up->var], 1);
990 res = reorder(isl_upoly_copy(rec->p[rec->n - 1]), r);
991
992 for (i = rec->n - 2; i >= 0; --i) {
993 res = isl_upoly_mul(res, isl_upoly_copy(base));
994 res = isl_upoly_sum(res, reorder(isl_upoly_copy(rec->p[i]), r));
995 }
996
997 isl_upoly_free(base);
998 isl_upoly_free(up);
999
1000 return res;
1001 error:
1002 isl_upoly_free(up);
1003 return NULL;
1004 }
1005
1006 static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
1007 {
1008 int n_row, n_col;
1009 int equal;
1010
1011 isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
1012 div1->n_col >= div2->n_col, return -1);
1013
1014 if (div1->n_row == div2->n_row)
1015 return isl_mat_is_equal(div1, div2);
1016
1017 n_row = div1->n_row;
1018 n_col = div1->n_col;
1019 div1->n_row = div2->n_row;
1020 div1->n_col = div2->n_col;
1021
1022 equal = isl_mat_is_equal(div1, div2);
1023
1024 div1->n_row = n_row;
1025 div1->n_col = n_col;
1026
1027 return equal;
1028 }
1029
1030 static void expand_row(__isl_keep isl_mat *dst, int d,
1031 __isl_keep isl_mat *src, int s, int *exp)
1032 {
1033 int i;
1034 unsigned c = src->n_col - src->n_row;
1035
1036 isl_seq_cpy(dst->row[d], src->row[s], c);
1037 isl_seq_clr(dst->row[d] + c, dst->n_col - c);
1038
1039 for (i = 0; i < s; ++i)
1040 isl_int_set(dst->row[d][c + exp[i]], src->row[s][c + i]);
1041 }
1042
1043 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
1044 {
1045 int li, lj;
1046
1047 li = isl_seq_last_non_zero(div->row[i], div->n_col);
1048 lj = isl_seq_last_non_zero(div->row[j], div->n_col);
1049
1050 if (li != lj)
1051 return li - lj;
1052
1053 return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
1054 }
1055
1056 struct isl_div_sort_info {
1057 isl_mat *div;
1058 int row;
1059 };
1060
1061 static int div_sort_cmp(const void *p1, const void *p2)
1062 {
1063 const struct isl_div_sort_info *i1, *i2;
1064 i1 = (const struct isl_div_sort_info *) p1;
1065 i2 = (const struct isl_div_sort_info *) p2;
1066
1067 return cmp_row(i1->div, i1->row, i2->row);
1068 }
1069
1070 static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
1071 {
1072 int i;
1073 struct isl_div_sort_info *array = NULL;
1074 int *pos = NULL;
1075 int *reordering = NULL;
1076 unsigned div_pos;
1077
1078 if (!qp)
1079 return NULL;
1080 if (qp->div->n_row <= 1)
1081 return qp;
1082
1083 div_pos = isl_dim_total(qp->dim);
1084
1085 array = isl_alloc_array(qp->div->ctx, struct isl_div_sort_info,
1086 qp->div->n_row);
1087 pos = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1088 reordering = isl_alloc_array(qp->div->ctx, int, qp->div->n_col - 2);
1089 if (!array || !pos || !reordering)
1090 goto error;
1091
1092 for (i = 0; i < qp->div->n_row; ++i) {
1093 array[i].div = qp->div;
1094 array[i].row = i;
1095 pos[i] = i;
1096 }
1097
1098 qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info),
1099 div_sort_cmp);
1100
1101 for (i = 0; i < div_pos; ++i)
1102 reordering[i] = i;
1103
1104 for (i = 0; i < qp->div->n_row; ++i)
1105 reordering[div_pos + array[i].row] = div_pos + i;
1106
1107 for (i = 0; i < qp->div->n_row; ++i) {
1108 int t;
1109 if (pos[array[i].row] == i)
1110 continue;
1111 qp->div = isl_mat_cow(qp->div);
1112 qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]);
1113 t = pos[array[i].row];
1114 pos[array[i].row] = pos[i];
1115 pos[i] = t;
1116 }
1117
1118 qp->upoly = reorder(qp->upoly, reordering);
1119
1120 if (!qp->upoly || !qp->div)
1121 goto error;
1122
1123 free(pos);
1124 free(array);
1125 free(reordering);
1126
1127 return qp;
1128 error:
1129 free(pos);
1130 free(array);
1131 free(reordering);
1132 isl_qpolynomial_free(qp);
1133 return NULL;
1134 }
1135
1136 static __isl_give isl_mat *merge_divs(__isl_keep isl_mat *div1,
1137 __isl_keep isl_mat *div2, int *exp1, int *exp2)
1138 {
1139 int i, j, k;
1140 isl_mat *div = NULL;
1141 unsigned d = div1->n_col - div1->n_row;
1142
1143 div = isl_mat_alloc(div1->ctx, 1 + div1->n_row + div2->n_row,
1144 d + div1->n_row + div2->n_row);
1145 if (!div)
1146 return NULL;
1147
1148 for (i = 0, j = 0, k = 0; i < div1->n_row && j < div2->n_row; ++k) {
1149 int cmp;
1150
1151 expand_row(div, k, div1, i, exp1);
1152 expand_row(div, k + 1, div2, j, exp2);
1153
1154 cmp = cmp_row(div, k, k + 1);
1155 if (cmp == 0) {
1156 exp1[i++] = k;
1157 exp2[j++] = k;
1158 } else if (cmp < 0) {
1159 exp1[i++] = k;
1160 } else {
1161 exp2[j++] = k;
1162 isl_seq_cpy(div->row[k], div->row[k + 1], div->n_col);
1163 }
1164 }
1165 for (; i < div1->n_row; ++i, ++k) {
1166 expand_row(div, k, div1, i, exp1);
1167 exp1[i] = k;
1168 }
1169 for (; j < div2->n_row; ++j, ++k) {
1170 expand_row(div, k, div2, j, exp2);
1171 exp2[j] = k;
1172 }
1173
1174 div->n_row = k;
1175 div->n_col = d + k;
1176
1177 return div;
1178 }
1179
1180 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
1181 int *exp, int first)
1182 {
1183 int i;
1184 struct isl_upoly_rec *rec;
1185
1186 if (isl_upoly_is_cst(up))
1187 return up;
1188
1189 if (up->var < first)
1190 return up;
1191
1192 if (exp[up->var - first] == up->var - first)
1193 return up;
1194
1195 up = isl_upoly_cow(up);
1196 if (!up)
1197 goto error;
1198
1199 up->var = exp[up->var - first] + first;
1200
1201 rec = isl_upoly_as_rec(up);
1202 if (!rec)
1203 goto error;
1204
1205 for (i = 0; i < rec->n; ++i) {
1206 rec->p[i] = expand(rec->p[i], exp, first);
1207 if (!rec->p[i])
1208 goto error;
1209 }
1210
1211 return up;
1212 error:
1213 isl_upoly_free(up);
1214 return NULL;
1215 }
1216
1217 static __isl_give isl_qpolynomial *with_merged_divs(
1218 __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
1219 __isl_take isl_qpolynomial *qp2),
1220 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1221 {
1222 int *exp1 = NULL;
1223 int *exp2 = NULL;
1224 isl_mat *div = NULL;
1225
1226 qp1 = isl_qpolynomial_cow(qp1);
1227 qp2 = isl_qpolynomial_cow(qp2);
1228
1229 if (!qp1 || !qp2)
1230 goto error;
1231
1232 isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1233 qp1->div->n_col >= qp2->div->n_col, goto error);
1234
1235 exp1 = isl_alloc_array(qp1->div->ctx, int, qp1->div->n_row);
1236 exp2 = isl_alloc_array(qp2->div->ctx, int, qp2->div->n_row);
1237 if (!exp1 || !exp2)
1238 goto error;
1239
1240 div = merge_divs(qp1->div, qp2->div, exp1, exp2);
1241 if (!div)
1242 goto error;
1243
1244 isl_mat_free(qp1->div);
1245 qp1->div = isl_mat_copy(div);
1246 isl_mat_free(qp2->div);
1247 qp2->div = isl_mat_copy(div);
1248
1249 qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1250 qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1251
1252 if (!qp1->upoly || !qp2->upoly)
1253 goto error;
1254
1255 isl_mat_free(div);
1256 free(exp1);
1257 free(exp2);
1258
1259 return fn(qp1, qp2);
1260 error:
1261 isl_mat_free(div);
1262 free(exp1);
1263 free(exp2);
1264 isl_qpolynomial_free(qp1);
1265 isl_qpolynomial_free(qp2);
1266 return NULL;
1267 }
1268
1269 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1270 __isl_take isl_qpolynomial *qp2)
1271 {
1272 qp1 = isl_qpolynomial_cow(qp1);
1273
1274 if (!qp1 || !qp2)
1275 goto error;
1276
1277 if (qp1->div->n_row < qp2->div->n_row)
1278 return isl_qpolynomial_add(qp2, qp1);
1279
1280 isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1281 if (!compatible_divs(qp1->div, qp2->div))
1282 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1283
1284 qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1285 if (!qp1->upoly)
1286 goto error;
1287
1288 isl_qpolynomial_free(qp2);
1289
1290 return qp1;
1291 error:
1292 isl_qpolynomial_free(qp1);
1293 isl_qpolynomial_free(qp2);
1294 return NULL;
1295 }
1296
1297 __isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
1298 __isl_keep isl_set *dom,
1299 __isl_take isl_qpolynomial *qp1,
1300 __isl_take isl_qpolynomial *qp2)
1301 {
1302 return isl_qpolynomial_add(qp1, qp2);
1303 }
1304
1305 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1306 __isl_take isl_qpolynomial *qp2)
1307 {
1308 return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1309 }
1310
1311 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1312 {
1313 qp = isl_qpolynomial_cow(qp);
1314
1315 if (!qp)
1316 return NULL;
1317
1318 qp->upoly = isl_upoly_neg(qp->upoly);
1319 if (!qp->upoly)
1320 goto error;
1321
1322 return qp;
1323 error:
1324 isl_qpolynomial_free(qp);
1325 return NULL;
1326 }
1327
1328 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1329 __isl_take isl_qpolynomial *qp2)
1330 {
1331 qp1 = isl_qpolynomial_cow(qp1);
1332
1333 if (!qp1 || !qp2)
1334 goto error;
1335
1336 if (qp1->div->n_row < qp2->div->n_row)
1337 return isl_qpolynomial_mul(qp2, qp1);
1338
1339 isl_assert(qp1->dim->ctx, isl_dim_equal(qp1->dim, qp2->dim), goto error);
1340 if (!compatible_divs(qp1->div, qp2->div))
1341 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1342
1343 qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1344 if (!qp1->upoly)
1345 goto error;
1346
1347 isl_qpolynomial_free(qp2);
1348
1349 return qp1;
1350 error:
1351 isl_qpolynomial_free(qp1);
1352 isl_qpolynomial_free(qp2);
1353 return NULL;
1354 }
1355
1356 __isl_give isl_qpolynomial *isl_qpolynomial_zero(__isl_take isl_dim *dim)
1357 {
1358 return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1359 }
1360
1361 __isl_give isl_qpolynomial *isl_qpolynomial_infty(__isl_take isl_dim *dim)
1362 {
1363 return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
1364 }
1365
1366 __isl_give isl_qpolynomial *isl_qpolynomial_neginfty(__isl_take isl_dim *dim)
1367 {
1368 return isl_qpolynomial_alloc(dim, 0, isl_upoly_neginfty(dim->ctx));
1369 }
1370
1371 __isl_give isl_qpolynomial *isl_qpolynomial_nan(__isl_take isl_dim *dim)
1372 {
1373 return isl_qpolynomial_alloc(dim, 0, isl_upoly_nan(dim->ctx));
1374 }
1375
1376 __isl_give isl_qpolynomial *isl_qpolynomial_cst(__isl_take isl_dim *dim,
1377 isl_int v)
1378 {
1379 struct isl_qpolynomial *qp;
1380 struct isl_upoly_cst *cst;
1381
1382 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1383 if (!qp)
1384 return NULL;
1385
1386 cst = isl_upoly_as_cst(qp->upoly);
1387 isl_int_set(cst->n, v);
1388
1389 return qp;
1390 }
1391
1392 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1393 isl_int *n, isl_int *d)
1394 {
1395 struct isl_upoly_cst *cst;
1396
1397 if (!qp)
1398 return -1;
1399
1400 if (!isl_upoly_is_cst(qp->upoly))
1401 return 0;
1402
1403 cst = isl_upoly_as_cst(qp->upoly);
1404 if (!cst)
1405 return -1;
1406
1407 if (n)
1408 isl_int_set(*n, cst->n);
1409 if (d)
1410 isl_int_set(*d, cst->d);
1411
1412 return 1;
1413 }
1414
1415 int isl_upoly_is_affine(__isl_keep struct isl_upoly *up)
1416 {
1417 int is_cst;
1418 struct isl_upoly_rec *rec;
1419
1420 if (!up)
1421 return -1;
1422
1423 if (up->var < 0)
1424 return 1;
1425
1426 rec = isl_upoly_as_rec(up);
1427 if (!rec)
1428 return -1;
1429
1430 if (rec->n > 2)
1431 return 0;
1432
1433 isl_assert(up->ctx, rec->n > 1, return -1);
1434
1435 is_cst = isl_upoly_is_cst(rec->p[1]);
1436 if (is_cst < 0)
1437 return -1;
1438 if (!is_cst)
1439 return 0;
1440
1441 return isl_upoly_is_affine(rec->p[0]);
1442 }
1443
1444 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp)
1445 {
1446 if (!qp)
1447 return -1;
1448
1449 if (qp->div->n_row > 0)
1450 return 0;
1451
1452 return isl_upoly_is_affine(qp->upoly);
1453 }
1454
1455 static void update_coeff(__isl_keep isl_vec *aff,
1456 __isl_keep struct isl_upoly_cst *cst, int pos)
1457 {
1458 isl_int gcd;
1459 isl_int f;
1460
1461 if (isl_int_is_zero(cst->n))
1462 return;
1463
1464 isl_int_init(gcd);
1465 isl_int_init(f);
1466 isl_int_gcd(gcd, cst->d, aff->el[0]);
1467 isl_int_divexact(f, cst->d, gcd);
1468 isl_int_divexact(gcd, aff->el[0], gcd);
1469 isl_seq_scale(aff->el, aff->el, f, aff->size);
1470 isl_int_mul(aff->el[1 + pos], gcd, cst->n);
1471 isl_int_clear(gcd);
1472 isl_int_clear(f);
1473 }
1474
1475 int isl_upoly_update_affine(__isl_keep struct isl_upoly *up,
1476 __isl_keep isl_vec *aff)
1477 {
1478 struct isl_upoly_cst *cst;
1479 struct isl_upoly_rec *rec;
1480
1481 if (!up || !aff)
1482 return -1;
1483
1484 if (up->var < 0) {
1485 struct isl_upoly_cst *cst;
1486
1487 cst = isl_upoly_as_cst(up);
1488 if (!cst)
1489 return -1;
1490 update_coeff(aff, cst, 0);
1491 return 0;
1492 }
1493
1494 rec = isl_upoly_as_rec(up);
1495 if (!rec)
1496 return -1;
1497 isl_assert(up->ctx, rec->n == 2, return -1);
1498
1499 cst = isl_upoly_as_cst(rec->p[1]);
1500 if (!cst)
1501 return -1;
1502 update_coeff(aff, cst, 1 + up->var);
1503
1504 return isl_upoly_update_affine(rec->p[0], aff);
1505 }
1506
1507 __isl_give isl_vec *isl_qpolynomial_extract_affine(
1508 __isl_keep isl_qpolynomial *qp)
1509 {
1510 isl_vec *aff;
1511 unsigned d;
1512
1513 if (!qp)
1514 return NULL;
1515
1516 isl_assert(qp->div->ctx, qp->div->n_row == 0, return NULL);
1517 d = isl_dim_total(qp->dim);
1518 aff = isl_vec_alloc(qp->div->ctx, 2 + d);
1519 if (!aff)
1520 return NULL;
1521
1522 isl_seq_clr(aff->el + 1, 1 + d);
1523 isl_int_set_si(aff->el[0], 1);
1524
1525 if (isl_upoly_update_affine(qp->upoly, aff) < 0)
1526 goto error;
1527
1528 return aff;
1529 error:
1530 isl_vec_free(aff);
1531 return NULL;
1532 }
1533
1534 int isl_qpolynomial_is_equal(__isl_keep isl_qpolynomial *qp1,
1535 __isl_keep isl_qpolynomial *qp2)
1536 {
1537 if (!qp1 || !qp2)
1538 return -1;
1539
1540 return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
1541 }
1542
1543 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
1544 {
1545 int i;
1546 struct isl_upoly_rec *rec;
1547
1548 if (isl_upoly_is_cst(up)) {
1549 struct isl_upoly_cst *cst;
1550 cst = isl_upoly_as_cst(up);
1551 if (!cst)
1552 return;
1553 isl_int_lcm(*d, *d, cst->d);
1554 return;
1555 }
1556
1557 rec = isl_upoly_as_rec(up);
1558 if (!rec)
1559 return;
1560
1561 for (i = 0; i < rec->n; ++i)
1562 upoly_update_den(rec->p[i], d);
1563 }
1564
1565 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
1566 {
1567 isl_int_set_si(*d, 1);
1568 if (!qp)
1569 return;
1570 upoly_update_den(qp->upoly, d);
1571 }
1572
1573 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_dim *dim,
1574 int pos, int power)
1575 {
1576 struct isl_ctx *ctx;
1577
1578 if (!dim)
1579 return NULL;
1580
1581 ctx = dim->ctx;
1582
1583 return isl_qpolynomial_alloc(dim, 0, isl_upoly_pow(ctx, pos, power));
1584 }
1585
1586 __isl_give isl_qpolynomial *isl_qpolynomial_var(__isl_take isl_dim *dim,
1587 enum isl_dim_type type, unsigned pos)
1588 {
1589 if (!dim)
1590 return NULL;
1591
1592 isl_assert(dim->ctx, isl_dim_size(dim, isl_dim_in) == 0, goto error);
1593 isl_assert(dim->ctx, pos < isl_dim_size(dim, type), goto error);
1594
1595 if (type == isl_dim_set)
1596 pos += isl_dim_size(dim, isl_dim_param);
1597
1598 return isl_qpolynomial_pow(dim, pos, 1);
1599 error:
1600 isl_dim_free(dim);
1601 return NULL;
1602 }
1603
1604 __isl_give isl_qpolynomial *isl_qpolynomial_div_pow(__isl_take isl_div *div,
1605 int power)
1606 {
1607 struct isl_qpolynomial *qp = NULL;
1608 struct isl_upoly_rec *rec;
1609 struct isl_upoly_cst *cst;
1610 int i;
1611 int pos;
1612
1613 if (!div)
1614 return NULL;
1615 isl_assert(div->ctx, div->bmap->n_div == 1, goto error);
1616
1617 pos = isl_dim_total(div->bmap->dim);
1618 rec = isl_upoly_alloc_rec(div->ctx, pos, 1 + power);
1619 qp = isl_qpolynomial_alloc(isl_basic_map_get_dim(div->bmap), 1,
1620 &rec->up);
1621 if (!qp)
1622 goto error;
1623
1624 isl_seq_cpy(qp->div->row[0], div->line[0], qp->div->n_col - 1);
1625 isl_int_set_si(qp->div->row[0][qp->div->n_col - 1], 0);
1626
1627 for (i = 0; i < 1 + power; ++i) {
1628 rec->p[i] = isl_upoly_zero(div->ctx);
1629 if (!rec->p[i])
1630 goto error;
1631 rec->n++;
1632 }
1633 cst = isl_upoly_as_cst(rec->p[power]);
1634 isl_int_set_si(cst->n, 1);
1635
1636 isl_div_free(div);
1637
1638 return qp;
1639 error:
1640 isl_qpolynomial_free(qp);
1641 isl_div_free(div);
1642 return NULL;
1643 }
1644
1645 __isl_give isl_qpolynomial *isl_qpolynomial_div(__isl_take isl_div *div)
1646 {
1647 return isl_qpolynomial_div_pow(div, 1);
1648 }
1649
1650 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst(__isl_take isl_dim *dim,
1651 const isl_int n, const isl_int d)
1652 {
1653 struct isl_qpolynomial *qp;
1654 struct isl_upoly_cst *cst;
1655
1656 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1657 if (!qp)
1658 return NULL;
1659
1660 cst = isl_upoly_as_cst(qp->upoly);
1661 isl_int_set(cst->n, n);
1662 isl_int_set(cst->d, d);
1663
1664 return qp;
1665 }
1666
1667 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
1668 {
1669 struct isl_upoly_rec *rec;
1670 int i;
1671
1672 if (!up)
1673 return -1;
1674
1675 if (isl_upoly_is_cst(up))
1676 return 0;
1677
1678 if (up->var < d)
1679 active[up->var] = 1;
1680
1681 rec = isl_upoly_as_rec(up);
1682 for (i = 0; i < rec->n; ++i)
1683 if (up_set_active(rec->p[i], active, d) < 0)
1684 return -1;
1685
1686 return 0;
1687 }
1688
1689 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
1690 {
1691 int i, j;
1692 int d = isl_dim_total(qp->dim);
1693
1694 if (!qp || !active)
1695 return -1;
1696
1697 for (i = 0; i < d; ++i)
1698 for (j = 0; j < qp->div->n_row; ++j) {
1699 if (isl_int_is_zero(qp->div->row[j][2 + i]))
1700 continue;
1701 active[i] = 1;
1702 break;
1703 }
1704
1705 return up_set_active(qp->upoly, active, d);
1706 }
1707
1708 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
1709 enum isl_dim_type type, unsigned first, unsigned n)
1710 {
1711 int i;
1712 int *active = NULL;
1713 int involves = 0;
1714
1715 if (!qp)
1716 return -1;
1717 if (n == 0)
1718 return 0;
1719
1720 isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
1721 return -1);
1722 isl_assert(qp->dim->ctx, type == isl_dim_param ||
1723 type == isl_dim_set, return -1);
1724
1725 active = isl_calloc_array(set->ctx, int, isl_dim_total(qp->dim));
1726 if (set_active(qp, active) < 0)
1727 goto error;
1728
1729 if (type == isl_dim_set)
1730 first += isl_dim_size(qp->dim, isl_dim_param);
1731 for (i = 0; i < n; ++i)
1732 if (active[first + i]) {
1733 involves = 1;
1734 break;
1735 }
1736
1737 free(active);
1738
1739 return involves;
1740 error:
1741 free(active);
1742 return -1;
1743 }
1744
1745 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
1746 unsigned first, unsigned n)
1747 {
1748 int i;
1749 struct isl_upoly_rec *rec;
1750
1751 if (!up)
1752 return NULL;
1753 if (n == 0 || up->var < 0 || up->var < first)
1754 return up;
1755 if (up->var < first + n) {
1756 up = replace_by_constant_term(up);
1757 return isl_upoly_drop(up, first, n);
1758 }
1759 up = isl_upoly_cow(up);
1760 if (!up)
1761 return NULL;
1762 up->var -= n;
1763 rec = isl_upoly_as_rec(up);
1764 if (!rec)
1765 goto error;
1766
1767 for (i = 0; i < rec->n; ++i) {
1768 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
1769 if (!rec->p[i])
1770 goto error;
1771 }
1772
1773 return up;
1774 error:
1775 isl_upoly_free(up);
1776 return NULL;
1777 }
1778
1779 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
1780 __isl_take isl_qpolynomial *qp,
1781 enum isl_dim_type type, unsigned first, unsigned n)
1782 {
1783 if (!qp)
1784 return NULL;
1785 if (n == 0 && !isl_dim_get_tuple_name(qp->dim, type))
1786 return qp;
1787
1788 qp = isl_qpolynomial_cow(qp);
1789 if (!qp)
1790 return NULL;
1791
1792 isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
1793 goto error);
1794 isl_assert(qp->dim->ctx, type == isl_dim_param ||
1795 type == isl_dim_set, goto error);
1796
1797 qp->dim = isl_dim_drop(qp->dim, type, first, n);
1798 if (!qp->dim)
1799 goto error;
1800
1801 if (type == isl_dim_set)
1802 first += isl_dim_size(qp->dim, isl_dim_param);
1803
1804 qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
1805 if (!qp->div)
1806 goto error;
1807
1808 qp->upoly = isl_upoly_drop(qp->upoly, first, n);
1809 if (!qp->upoly)
1810 goto error;
1811
1812 return qp;
1813 error:
1814 isl_qpolynomial_free(qp);
1815 return NULL;
1816 }
1817
1818 __isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
1819 unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
1820 {
1821 int i;
1822 struct isl_upoly_rec *rec;
1823 struct isl_upoly *base, *res;
1824
1825 if (!up)
1826 return NULL;
1827
1828 if (isl_upoly_is_cst(up))
1829 return up;
1830
1831 if (up->var < first)
1832 return up;
1833
1834 rec = isl_upoly_as_rec(up);
1835 if (!rec)
1836 goto error;
1837
1838 isl_assert(up->ctx, rec->n >= 1, goto error);
1839
1840 if (up->var >= first + n)
1841 base = isl_upoly_pow(up->ctx, up->var, 1);
1842 else
1843 base = isl_upoly_copy(subs[up->var - first]);
1844
1845 res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
1846 for (i = rec->n - 2; i >= 0; --i) {
1847 struct isl_upoly *t;
1848 t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
1849 res = isl_upoly_mul(res, isl_upoly_copy(base));
1850 res = isl_upoly_sum(res, t);
1851 }
1852
1853 isl_upoly_free(base);
1854 isl_upoly_free(up);
1855
1856 return res;
1857 error:
1858 isl_upoly_free(up);
1859 return NULL;
1860 }
1861
1862 __isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
1863 isl_int denom, unsigned len)
1864 {
1865 int i;
1866 struct isl_upoly *up;
1867
1868 isl_assert(ctx, len >= 1, return NULL);
1869
1870 up = isl_upoly_rat_cst(ctx, f[0], denom);
1871 for (i = 0; i < len - 1; ++i) {
1872 struct isl_upoly *t;
1873 struct isl_upoly *c;
1874
1875 if (isl_int_is_zero(f[1 + i]))
1876 continue;
1877
1878 c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
1879 t = isl_upoly_pow(ctx, i, 1);
1880 t = isl_upoly_mul(c, t);
1881 up = isl_upoly_sum(up, t);
1882 }
1883
1884 return up;
1885 }
1886
1887 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
1888 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
1889 {
1890 int i, j, k;
1891 isl_int denom;
1892 unsigned total;
1893 struct isl_upoly *up;
1894
1895 if (!eq)
1896 goto error;
1897 if (eq->n_eq == 0) {
1898 isl_basic_set_free(eq);
1899 return qp;
1900 }
1901
1902 qp = isl_qpolynomial_cow(qp);
1903 if (!qp)
1904 goto error;
1905 qp->div = isl_mat_cow(qp->div);
1906 if (!qp->div)
1907 goto error;
1908
1909 total = 1 + isl_dim_total(eq->dim);
1910 isl_int_init(denom);
1911 for (i = 0; i < eq->n_eq; ++i) {
1912 j = isl_seq_last_non_zero(eq->eq[i], total);
1913 if (j < 0 || j == 0)
1914 continue;
1915
1916 for (k = 0; k < qp->div->n_row; ++k) {
1917 if (isl_int_is_zero(qp->div->row[k][1 + j]))
1918 continue;
1919 isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
1920 &qp->div->row[k][0]);
1921 isl_seq_normalize(qp->div->ctx,
1922 qp->div->row[k], 1 + total);
1923 }
1924
1925 if (isl_int_is_pos(eq->eq[i][j]))
1926 isl_seq_neg(eq->eq[i], eq->eq[i], total);
1927 isl_int_abs(denom, eq->eq[i][j]);
1928 isl_int_set_si(eq->eq[i][j], 0);
1929
1930 up = isl_upoly_from_affine(qp->dim->ctx,
1931 eq->eq[i], denom, total);
1932 qp->upoly = isl_upoly_subs(qp->upoly, j - 1, 1, &up);
1933 isl_upoly_free(up);
1934 }
1935 isl_int_clear(denom);
1936
1937 if (!qp->upoly)
1938 goto error;
1939
1940 isl_basic_set_free(eq);
1941
1942 qp = sort_divs(qp);
1943
1944 return qp;
1945 error:
1946 isl_basic_set_free(eq);
1947 isl_qpolynomial_free(qp);
1948 return NULL;
1949 }
1950
1951 #undef PW
1952 #define PW isl_pw_qpolynomial
1953 #undef EL
1954 #define EL isl_qpolynomial
1955 #undef IS_ZERO
1956 #define IS_ZERO is_zero
1957 #undef FIELD
1958 #define FIELD qp
1959
1960 #include <isl_pw_templ.c>
1961
1962 #undef UNION
1963 #define UNION isl_union_pw_qpolynomial
1964 #undef PART
1965 #define PART isl_pw_qpolynomial
1966 #undef PARTS
1967 #define PARTS pw_qpolynomial
1968
1969 #include <isl_union_templ.c>
1970
1971 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
1972 {
1973 if (!pwqp)
1974 return -1;
1975
1976 if (pwqp->n != -1)
1977 return 0;
1978
1979 if (!isl_set_fast_is_universe(pwqp->p[0].set))
1980 return 0;
1981
1982 return isl_qpolynomial_is_one(pwqp->p[0].qp);
1983 }
1984
1985 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
1986 __isl_take isl_pw_qpolynomial *pwqp1,
1987 __isl_take isl_pw_qpolynomial *pwqp2)
1988 {
1989 int i, j, n;
1990 struct isl_pw_qpolynomial *res;
1991 isl_set *set;
1992
1993 if (!pwqp1 || !pwqp2)
1994 goto error;
1995
1996 isl_assert(pwqp1->dim->ctx, isl_dim_equal(pwqp1->dim, pwqp2->dim),
1997 goto error);
1998
1999 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
2000 isl_pw_qpolynomial_free(pwqp2);
2001 return pwqp1;
2002 }
2003
2004 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
2005 isl_pw_qpolynomial_free(pwqp1);
2006 return pwqp2;
2007 }
2008
2009 if (isl_pw_qpolynomial_is_one(pwqp1)) {
2010 isl_pw_qpolynomial_free(pwqp1);
2011 return pwqp2;
2012 }
2013
2014 if (isl_pw_qpolynomial_is_one(pwqp2)) {
2015 isl_pw_qpolynomial_free(pwqp2);
2016 return pwqp1;
2017 }
2018
2019 n = pwqp1->n * pwqp2->n;
2020 res = isl_pw_qpolynomial_alloc_(isl_dim_copy(pwqp1->dim), n);
2021
2022 for (i = 0; i < pwqp1->n; ++i) {
2023 for (j = 0; j < pwqp2->n; ++j) {
2024 struct isl_set *common;
2025 struct isl_qpolynomial *prod;
2026 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
2027 isl_set_copy(pwqp2->p[j].set));
2028 if (isl_set_fast_is_empty(common)) {
2029 isl_set_free(common);
2030 continue;
2031 }
2032
2033 prod = isl_qpolynomial_mul(
2034 isl_qpolynomial_copy(pwqp1->p[i].qp),
2035 isl_qpolynomial_copy(pwqp2->p[j].qp));
2036
2037 res = isl_pw_qpolynomial_add_piece(res, common, prod);
2038 }
2039 }
2040
2041 isl_pw_qpolynomial_free(pwqp1);
2042 isl_pw_qpolynomial_free(pwqp2);
2043
2044 return res;
2045 error:
2046 isl_pw_qpolynomial_free(pwqp1);
2047 isl_pw_qpolynomial_free(pwqp2);
2048 return NULL;
2049 }
2050
2051 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_neg(
2052 __isl_take isl_pw_qpolynomial *pwqp)
2053 {
2054 int i, j, n;
2055 struct isl_pw_qpolynomial *res;
2056 isl_set *set;
2057
2058 if (!pwqp)
2059 return NULL;
2060
2061 if (isl_pw_qpolynomial_is_zero(pwqp))
2062 return pwqp;
2063
2064 pwqp = isl_pw_qpolynomial_cow(pwqp);
2065
2066 for (i = 0; i < pwqp->n; ++i) {
2067 pwqp->p[i].qp = isl_qpolynomial_neg(pwqp->p[i].qp);
2068 if (!pwqp->p[i].qp)
2069 goto error;
2070 }
2071
2072 return pwqp;
2073 error:
2074 isl_pw_qpolynomial_free(pwqp);
2075 return NULL;
2076 }
2077
2078 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_sub(
2079 __isl_take isl_pw_qpolynomial *pwqp1,
2080 __isl_take isl_pw_qpolynomial *pwqp2)
2081 {
2082 return isl_pw_qpolynomial_add(pwqp1, isl_pw_qpolynomial_neg(pwqp2));
2083 }
2084
2085 __isl_give struct isl_upoly *isl_upoly_eval(
2086 __isl_take struct isl_upoly *up, __isl_take isl_vec *vec)
2087 {
2088 int i;
2089 struct isl_upoly_rec *rec;
2090 struct isl_upoly *res;
2091 struct isl_upoly *base;
2092
2093 if (isl_upoly_is_cst(up)) {
2094 isl_vec_free(vec);
2095 return up;
2096 }
2097
2098 rec = isl_upoly_as_rec(up);
2099 if (!rec)
2100 goto error;
2101
2102 isl_assert(up->ctx, rec->n >= 1, goto error);
2103
2104 base = isl_upoly_rat_cst(up->ctx, vec->el[1 + up->var], vec->el[0]);
2105
2106 res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
2107 isl_vec_copy(vec));
2108
2109 for (i = rec->n - 2; i >= 0; --i) {
2110 res = isl_upoly_mul(res, isl_upoly_copy(base));
2111 res = isl_upoly_sum(res,
2112 isl_upoly_eval(isl_upoly_copy(rec->p[i]),
2113 isl_vec_copy(vec)));
2114 }
2115
2116 isl_upoly_free(base);
2117 isl_upoly_free(up);
2118 isl_vec_free(vec);
2119 return res;
2120 error:
2121 isl_upoly_free(up);
2122 isl_vec_free(vec);
2123 return NULL;
2124 }
2125
2126 __isl_give isl_qpolynomial *isl_qpolynomial_eval(
2127 __isl_take isl_qpolynomial *qp, __isl_take isl_point *pnt)
2128 {
2129 isl_vec *ext;
2130 struct isl_upoly *up;
2131 isl_dim *dim;
2132
2133 if (!qp || !pnt)
2134 goto error;
2135 isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, qp->dim), goto error);
2136
2137 if (qp->div->n_row == 0)
2138 ext = isl_vec_copy(pnt->vec);
2139 else {
2140 int i;
2141 unsigned dim = isl_dim_total(qp->dim);
2142 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
2143 if (!ext)
2144 goto error;
2145
2146 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
2147 for (i = 0; i < qp->div->n_row; ++i) {
2148 isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
2149 1 + dim + i, &ext->el[1+dim+i]);
2150 isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
2151 qp->div->row[i][0]);
2152 }
2153 }
2154
2155 up = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
2156 if (!up)
2157 goto error;
2158
2159 dim = isl_dim_copy(qp->dim);
2160 isl_qpolynomial_free(qp);
2161 isl_point_free(pnt);
2162
2163 return isl_qpolynomial_alloc(dim, 0, up);
2164 error:
2165 isl_qpolynomial_free(qp);
2166 isl_point_free(pnt);
2167 return NULL;
2168 }
2169
2170 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
2171 __isl_keep struct isl_upoly_cst *cst2)
2172 {
2173 int cmp;
2174 isl_int t;
2175 isl_int_init(t);
2176 isl_int_mul(t, cst1->n, cst2->d);
2177 isl_int_submul(t, cst2->n, cst1->d);
2178 cmp = isl_int_sgn(t);
2179 isl_int_clear(t);
2180 return cmp;
2181 }
2182
2183 int isl_qpolynomial_le_cst(__isl_keep isl_qpolynomial *qp1,
2184 __isl_keep isl_qpolynomial *qp2)
2185 {
2186 struct isl_upoly_cst *cst1, *cst2;
2187
2188 if (!qp1 || !qp2)
2189 return -1;
2190 isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), return -1);
2191 isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), return -1);
2192 if (isl_qpolynomial_is_nan(qp1))
2193 return -1;
2194 if (isl_qpolynomial_is_nan(qp2))
2195 return -1;
2196 cst1 = isl_upoly_as_cst(qp1->upoly);
2197 cst2 = isl_upoly_as_cst(qp2->upoly);
2198
2199 return isl_upoly_cmp(cst1, cst2) <= 0;
2200 }
2201
2202 __isl_give isl_qpolynomial *isl_qpolynomial_min_cst(
2203 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
2204 {
2205 struct isl_upoly_cst *cst1, *cst2;
2206 int cmp;
2207
2208 if (!qp1 || !qp2)
2209 goto error;
2210 isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
2211 isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
2212 cst1 = isl_upoly_as_cst(qp1->upoly);
2213 cst2 = isl_upoly_as_cst(qp2->upoly);
2214 cmp = isl_upoly_cmp(cst1, cst2);
2215
2216 if (cmp <= 0) {
2217 isl_qpolynomial_free(qp2);
2218 } else {
2219 isl_qpolynomial_free(qp1);
2220 qp1 = qp2;
2221 }
2222 return qp1;
2223 error:
2224 isl_qpolynomial_free(qp1);
2225 isl_qpolynomial_free(qp2);
2226 return NULL;
2227 }
2228
2229 __isl_give isl_qpolynomial *isl_qpolynomial_max_cst(
2230 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
2231 {
2232 struct isl_upoly_cst *cst1, *cst2;
2233 int cmp;
2234
2235 if (!qp1 || !qp2)
2236 goto error;
2237 isl_assert(qp1->dim->ctx, isl_upoly_is_cst(qp1->upoly), goto error);
2238 isl_assert(qp2->dim->ctx, isl_upoly_is_cst(qp2->upoly), goto error);
2239 cst1 = isl_upoly_as_cst(qp1->upoly);
2240 cst2 = isl_upoly_as_cst(qp2->upoly);
2241 cmp = isl_upoly_cmp(cst1, cst2);
2242
2243 if (cmp >= 0) {
2244 isl_qpolynomial_free(qp2);
2245 } else {
2246 isl_qpolynomial_free(qp1);
2247 qp1 = qp2;
2248 }
2249 return qp1;
2250 error:
2251 isl_qpolynomial_free(qp1);
2252 isl_qpolynomial_free(qp2);
2253 return NULL;
2254 }
2255
2256 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
2257 __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
2258 unsigned first, unsigned n)
2259 {
2260 unsigned total;
2261 unsigned g_pos;
2262 int *exp;
2263
2264 if (n == 0)
2265 return qp;
2266
2267 qp = isl_qpolynomial_cow(qp);
2268 if (!qp)
2269 return NULL;
2270
2271 isl_assert(qp->div->ctx, first <= isl_dim_size(qp->dim, type),
2272 goto error);
2273
2274 g_pos = pos(qp->dim, type) + first;
2275
2276 qp->div = isl_mat_insert_cols(qp->div, 2 + g_pos, n);
2277 if (!qp->div)
2278 goto error;
2279
2280 total = qp->div->n_col - 2;
2281 if (total > g_pos) {
2282 int i;
2283 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
2284 if (!exp)
2285 goto error;
2286 for (i = 0; i < total - g_pos; ++i)
2287 exp[i] = i + n;
2288 qp->upoly = expand(qp->upoly, exp, g_pos);
2289 free(exp);
2290 if (!qp->upoly)
2291 goto error;
2292 }
2293
2294 qp->dim = isl_dim_insert(qp->dim, type, first, n);
2295 if (!qp->dim)
2296 goto error;
2297
2298 return qp;
2299 error:
2300 isl_qpolynomial_free(qp);
2301 return NULL;
2302 }
2303
2304 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
2305 __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
2306 {
2307 unsigned pos;
2308
2309 pos = isl_qpolynomial_dim(qp, type);
2310
2311 return isl_qpolynomial_insert_dims(qp, type, pos, n);
2312 }
2313
2314 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
2315 __isl_take isl_pw_qpolynomial *pwqp,
2316 enum isl_dim_type type, unsigned n)
2317 {
2318 unsigned pos;
2319
2320 pos = isl_pw_qpolynomial_dim(pwqp, type);
2321
2322 return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
2323 }
2324
2325 static int *reordering_move(isl_ctx *ctx,
2326 unsigned len, unsigned dst, unsigned src, unsigned n)
2327 {
2328 int i;
2329 int *reordering;
2330
2331 reordering = isl_alloc_array(ctx, int, len);
2332 if (!reordering)
2333 return NULL;
2334
2335 if (dst <= src) {
2336 for (i = 0; i < dst; ++i)
2337 reordering[i] = i;
2338 for (i = 0; i < n; ++i)
2339 reordering[src + i] = dst + i;
2340 for (i = 0; i < src - dst; ++i)
2341 reordering[dst + i] = dst + n + i;
2342 for (i = 0; i < len - src - n; ++i)
2343 reordering[src + n + i] = src + n + i;
2344 } else {
2345 for (i = 0; i < src; ++i)
2346 reordering[i] = i;
2347 for (i = 0; i < n; ++i)
2348 reordering[src + i] = dst + i;
2349 for (i = 0; i < dst - src; ++i)
2350 reordering[src + n + i] = src + i;
2351 for (i = 0; i < len - dst - n; ++i)
2352 reordering[dst + n + i] = dst + n + i;
2353 }
2354
2355 return reordering;
2356 }
2357
2358 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
2359 __isl_take isl_qpolynomial *qp,
2360 enum isl_dim_type dst_type, unsigned dst_pos,
2361 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
2362 {
2363 unsigned g_dst_pos;
2364 unsigned g_src_pos;
2365 int *reordering;
2366
2367 qp = isl_qpolynomial_cow(qp);
2368 if (!qp)
2369 return NULL;
2370
2371 isl_assert(qp->dim->ctx, src_pos + n <= isl_dim_size(qp->dim, src_type),
2372 goto error);
2373
2374 g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
2375 g_src_pos = pos(qp->dim, src_type) + src_pos;
2376 if (dst_type > src_type)
2377 g_dst_pos -= n;
2378
2379 qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
2380 if (!qp->div)
2381 goto error;
2382 qp = sort_divs(qp);
2383 if (!qp)
2384 goto error;
2385
2386 reordering = reordering_move(qp->dim->ctx,
2387 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
2388 if (!reordering)
2389 goto error;
2390
2391 qp->upoly = reorder(qp->upoly, reordering);
2392 free(reordering);
2393 if (!qp->upoly)
2394 goto error;
2395
2396 qp->dim = isl_dim_move(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
2397 if (!qp->dim)
2398 goto error;
2399
2400 return qp;
2401 error:
2402 isl_qpolynomial_free(qp);
2403 return NULL;
2404 }
2405
2406 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_dim *dim,
2407 isl_int *f, isl_int denom)
2408 {
2409 struct isl_upoly *up;
2410
2411 if (!dim)
2412 return NULL;
2413
2414 up = isl_upoly_from_affine(dim->ctx, f, denom, 1 + isl_dim_total(dim));
2415
2416 return isl_qpolynomial_alloc(dim, 0, up);
2417 }
2418
2419 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
2420 __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
2421 {
2422 isl_int denom;
2423 isl_dim *dim;
2424 struct isl_upoly *up;
2425 isl_qpolynomial *qp;
2426 int sgn;
2427
2428 if (!c)
2429 return NULL;
2430
2431 isl_int_init(denom);
2432
2433 isl_constraint_get_coefficient(c, type, pos, &denom);
2434 isl_constraint_set_coefficient(c, type, pos, c->ctx->zero);
2435 sgn = isl_int_sgn(denom);
2436 isl_int_abs(denom, denom);
2437 up = isl_upoly_from_affine(c->ctx, c->line[0], denom,
2438 1 + isl_constraint_dim(c, isl_dim_all));
2439 if (sgn < 0)
2440 isl_int_neg(denom, denom);
2441 isl_constraint_set_coefficient(c, type, pos, denom);
2442
2443 dim = isl_dim_copy(c->bmap->dim);
2444
2445 isl_int_clear(denom);
2446 isl_constraint_free(c);
2447
2448 qp = isl_qpolynomial_alloc(dim, 0, up);
2449 if (sgn > 0)
2450 qp = isl_qpolynomial_neg(qp);
2451 return qp;
2452 }
2453
2454 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
2455 * in "qp" by subs[i].
2456 */
2457 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
2458 __isl_take isl_qpolynomial *qp,
2459 enum isl_dim_type type, unsigned first, unsigned n,
2460 __isl_keep isl_qpolynomial **subs)
2461 {
2462 int i;
2463 struct isl_upoly **ups;
2464
2465 if (n == 0)
2466 return qp;
2467
2468 qp = isl_qpolynomial_cow(qp);
2469 if (!qp)
2470 return NULL;
2471 for (i = 0; i < n; ++i)
2472 if (!subs[i])
2473 goto error;
2474
2475 isl_assert(qp->dim->ctx, first + n <= isl_dim_size(qp->dim, type),
2476 goto error);
2477
2478 for (i = 0; i < n; ++i)
2479 isl_assert(qp->dim->ctx, isl_dim_equal(qp->dim, subs[i]->dim),
2480 goto error);
2481
2482 isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
2483 for (i = 0; i < n; ++i)
2484 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
2485
2486 first += pos(qp->dim, type);
2487
2488 ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
2489 if (!ups)
2490 goto error;
2491 for (i = 0; i < n; ++i)
2492 ups[i] = subs[i]->upoly;
2493
2494 qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
2495
2496 free(ups);
2497
2498 if (!qp->upoly)
2499 goto error;
2500
2501 return qp;
2502 error:
2503 isl_qpolynomial_free(qp);
2504 return NULL;
2505 }
2506
2507 __isl_give isl_basic_set *add_div_constraints(__isl_take isl_basic_set *bset,
2508 __isl_take isl_mat *div)
2509 {
2510 int i;
2511 unsigned total;
2512
2513 if (!bset || !div)
2514 goto error;
2515
2516 bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
2517 if (!bset)
2518 goto error;
2519 total = isl_basic_set_total_dim(bset);
2520 for (i = 0; i < div->n_row; ++i)
2521 if (isl_basic_set_add_div_constraints_var(bset,
2522 total - div->n_row + i, div->row[i]) < 0)
2523 goto error;
2524
2525 isl_mat_free(div);
2526 return bset;
2527 error:
2528 isl_mat_free(div);
2529 isl_basic_set_free(bset);
2530 return NULL;
2531 }
2532
2533 /* Extend "bset" with extra set dimensions for each integer division
2534 * in "qp" and then call "fn" with the extended bset and the polynomial
2535 * that results from replacing each of the integer divisions by the
2536 * corresponding extra set dimension.
2537 */
2538 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
2539 __isl_keep isl_basic_set *bset,
2540 int (*fn)(__isl_take isl_basic_set *bset,
2541 __isl_take isl_qpolynomial *poly, void *user), void *user)
2542 {
2543 isl_dim *dim;
2544 isl_mat *div;
2545 isl_qpolynomial *poly;
2546
2547 if (!qp || !bset)
2548 goto error;
2549 if (qp->div->n_row == 0)
2550 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
2551 user);
2552
2553 div = isl_mat_copy(qp->div);
2554 dim = isl_dim_copy(qp->dim);
2555 dim = isl_dim_add(dim, isl_dim_set, qp->div->n_row);
2556 poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
2557 bset = isl_basic_set_copy(bset);
2558 bset = isl_basic_set_add(bset, isl_dim_set, qp->div->n_row);
2559 bset = add_div_constraints(bset, div);
2560
2561 return fn(bset, poly, user);
2562 error:
2563 return -1;
2564 }
2565
2566 /* Return total degree in variables first (inclusive) up to last (exclusive).
2567 */
2568 int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
2569 {
2570 int deg = -1;
2571 int i;
2572 struct isl_upoly_rec *rec;
2573
2574 if (!up)
2575 return -2;
2576 if (isl_upoly_is_zero(up))
2577 return -1;
2578 if (isl_upoly_is_cst(up) || up->var < first)
2579 return 0;
2580
2581 rec = isl_upoly_as_rec(up);
2582 if (!rec)
2583 return -2;
2584
2585 for (i = 0; i < rec->n; ++i) {
2586 int d;
2587
2588 if (isl_upoly_is_zero(rec->p[i]))
2589 continue;
2590 d = isl_upoly_degree(rec->p[i], first, last);
2591 if (up->var < last)
2592 d += i;
2593 if (d > deg)
2594 deg = d;
2595 }
2596
2597 return deg;
2598 }
2599
2600 /* Return total degree in set variables.
2601 */
2602 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
2603 {
2604 unsigned ovar;
2605 unsigned nvar;
2606
2607 if (!poly)
2608 return -2;
2609
2610 ovar = isl_dim_offset(poly->dim, isl_dim_set);
2611 nvar = isl_dim_size(poly->dim, isl_dim_set);
2612 return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
2613 }
2614
2615 __isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
2616 unsigned pos, int deg)
2617 {
2618 int i;
2619 struct isl_upoly_rec *rec;
2620
2621 if (!up)
2622 return NULL;
2623
2624 if (isl_upoly_is_cst(up) || up->var < pos) {
2625 if (deg == 0)
2626 return isl_upoly_copy(up);
2627 else
2628 return isl_upoly_zero(up->ctx);
2629 }
2630
2631 rec = isl_upoly_as_rec(up);
2632 if (!rec)
2633 return NULL;
2634
2635 if (up->var == pos) {
2636 if (deg < rec->n)
2637 return isl_upoly_copy(rec->p[deg]);
2638 else
2639 return isl_upoly_zero(up->ctx);
2640 }
2641
2642 up = isl_upoly_copy(up);
2643 up = isl_upoly_cow(up);
2644 rec = isl_upoly_as_rec(up);
2645 if (!rec)
2646 goto error;
2647
2648 for (i = 0; i < rec->n; ++i) {
2649 struct isl_upoly *t;
2650 t = isl_upoly_coeff(rec->p[i], pos, deg);
2651 if (!t)
2652 goto error;
2653 isl_upoly_free(rec->p[i]);
2654 rec->p[i] = t;
2655 }
2656
2657 return up;
2658 error:
2659 isl_upoly_free(up);
2660 return NULL;
2661 }
2662
2663 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
2664 */
2665 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
2666 __isl_keep isl_qpolynomial *qp,
2667 enum isl_dim_type type, unsigned t_pos, int deg)
2668 {
2669 unsigned g_pos;
2670 struct isl_upoly *up;
2671 isl_qpolynomial *c;
2672
2673 if (!qp)
2674 return NULL;
2675
2676 isl_assert(qp->div->ctx, t_pos < isl_dim_size(qp->dim, type),
2677 return NULL);
2678
2679 g_pos = pos(qp->dim, type) + t_pos;
2680 up = isl_upoly_coeff(qp->upoly, g_pos, deg);
2681
2682 c = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), qp->div->n_row, up);
2683 if (!c)
2684 return NULL;
2685 isl_mat_free(c->div);
2686 c->div = isl_mat_copy(qp->div);
2687 if (!c->div)
2688 goto error;
2689 return c;
2690 error:
2691 isl_qpolynomial_free(c);
2692 return NULL;
2693 }
2694
2695 /* Homogenize the polynomial in the variables first (inclusive) up to
2696 * last (exclusive) by inserting powers of variable first.
2697 * Variable first is assumed not to appear in the input.
2698 */
2699 __isl_give struct isl_upoly *isl_upoly_homogenize(
2700 __isl_take struct isl_upoly *up, int deg, int target,
2701 int first, int last)
2702 {
2703 int i;
2704 struct isl_upoly_rec *rec;
2705
2706 if (!up)
2707 return NULL;
2708 if (isl_upoly_is_zero(up))
2709 return up;
2710 if (deg == target)
2711 return up;
2712 if (isl_upoly_is_cst(up) || up->var < first) {
2713 struct isl_upoly *hom;
2714
2715 hom = isl_upoly_pow(up->ctx, first, target - deg);
2716 if (!hom)
2717 goto error;
2718 rec = isl_upoly_as_rec(hom);
2719 rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
2720
2721 return hom;
2722 }
2723
2724 up = isl_upoly_cow(up);
2725 rec = isl_upoly_as_rec(up);
2726 if (!rec)
2727 goto error;
2728
2729 for (i = 0; i < rec->n; ++i) {
2730 if (isl_upoly_is_zero(rec->p[i]))
2731 continue;
2732 rec->p[i] = isl_upoly_homogenize(rec->p[i],
2733 up->var < last ? deg + i : i, target,
2734 first, last);
2735 if (!rec->p[i])
2736 goto error;
2737 }
2738
2739 return up;
2740 error:
2741 isl_upoly_free(up);
2742 return NULL;
2743 }
2744
2745 /* Homogenize the polynomial in the set variables by introducing
2746 * powers of an extra set variable at position 0.
2747 */
2748 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
2749 __isl_take isl_qpolynomial *poly)
2750 {
2751 unsigned ovar;
2752 unsigned nvar;
2753 int deg = isl_qpolynomial_degree(poly);
2754
2755 if (deg < -1)
2756 goto error;
2757
2758 poly = isl_qpolynomial_insert_dims(poly, isl_dim_set, 0, 1);
2759 poly = isl_qpolynomial_cow(poly);
2760 if (!poly)
2761 goto error;
2762
2763 ovar = isl_dim_offset(poly->dim, isl_dim_set);
2764 nvar = isl_dim_size(poly->dim, isl_dim_set);
2765 poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
2766 ovar, ovar + nvar);
2767 if (!poly->upoly)
2768 goto error;
2769
2770 return poly;
2771 error:
2772 isl_qpolynomial_free(poly);
2773 return NULL;
2774 }
2775
2776 __isl_give isl_term *isl_term_alloc(__isl_take isl_dim *dim,
2777 __isl_take isl_mat *div)
2778 {
2779 isl_term *term;
2780 int n;
2781
2782 if (!dim || !div)
2783 goto error;
2784
2785 n = isl_dim_total(dim) + div->n_row;
2786
2787 term = isl_calloc(dim->ctx, struct isl_term,
2788 sizeof(struct isl_term) + (n - 1) * sizeof(int));
2789 if (!term)
2790 goto error;
2791
2792 term->ref = 1;
2793 term->dim = dim;
2794 term->div = div;
2795 isl_int_init(term->n);
2796 isl_int_init(term->d);
2797
2798 return term;
2799 error:
2800 isl_dim_free(dim);
2801 isl_mat_free(div);
2802 return NULL;
2803 }
2804
2805 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
2806 {
2807 if (!term)
2808 return NULL;
2809
2810 term->ref++;
2811 return term;
2812 }
2813
2814 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
2815 {
2816 int i;
2817 isl_term *dup;
2818 unsigned total;
2819
2820 if (term)
2821 return NULL;
2822
2823 total = isl_dim_total(term->dim) + term->div->n_row;
2824
2825 dup = isl_term_alloc(isl_dim_copy(term->dim), isl_mat_copy(term->div));
2826 if (!dup)
2827 return NULL;
2828
2829 isl_int_set(dup->n, term->n);
2830 isl_int_set(dup->d, term->d);
2831
2832 for (i = 0; i < total; ++i)
2833 dup->pow[i] = term->pow[i];
2834
2835 return dup;
2836 }
2837
2838 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
2839 {
2840 if (!term)
2841 return NULL;
2842
2843 if (term->ref == 1)
2844 return term;
2845 term->ref--;
2846 return isl_term_dup(term);
2847 }
2848
2849 void isl_term_free(__isl_take isl_term *term)
2850 {
2851 if (!term)
2852 return;
2853
2854 if (--term->ref > 0)
2855 return;
2856
2857 isl_dim_free(term->dim);
2858 isl_mat_free(term->div);
2859 isl_int_clear(term->n);
2860 isl_int_clear(term->d);
2861 free(term);
2862 }
2863
2864 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
2865 {
2866 if (!term)
2867 return 0;
2868
2869 switch (type) {
2870 case isl_dim_param:
2871 case isl_dim_in:
2872 case isl_dim_out: return isl_dim_size(term->dim, type);
2873 case isl_dim_div: return term->div->n_row;
2874 case isl_dim_all: return isl_dim_total(term->dim) + term->div->n_row;
2875 default: return 0;
2876 }
2877 }
2878
2879 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
2880 {
2881 return term ? term->dim->ctx : NULL;
2882 }
2883
2884 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
2885 {
2886 if (!term)
2887 return;
2888 isl_int_set(*n, term->n);
2889 }
2890
2891 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
2892 {
2893 if (!term)
2894 return;
2895 isl_int_set(*d, term->d);
2896 }
2897
2898 int isl_term_get_exp(__isl_keep isl_term *term,
2899 enum isl_dim_type type, unsigned pos)
2900 {
2901 if (!term)
2902 return -1;
2903
2904 isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
2905
2906 if (type >= isl_dim_set)
2907 pos += isl_dim_size(term->dim, isl_dim_param);
2908 if (type >= isl_dim_div)
2909 pos += isl_dim_size(term->dim, isl_dim_set);
2910
2911 return term->pow[pos];
2912 }
2913
2914 __isl_give isl_div *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
2915 {
2916 isl_basic_map *bmap;
2917 unsigned total;
2918 int k;
2919
2920 if (!term)
2921 return NULL;
2922
2923 isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
2924 return NULL);
2925
2926 total = term->div->n_col - term->div->n_row - 2;
2927 /* No nested divs for now */
2928 isl_assert(term->dim->ctx,
2929 isl_seq_first_non_zero(term->div->row[pos] + 2 + total,
2930 term->div->n_row) == -1,
2931 return NULL);
2932
2933 bmap = isl_basic_map_alloc_dim(isl_dim_copy(term->dim), 1, 0, 0);
2934 if ((k = isl_basic_map_alloc_div(bmap)) < 0)
2935 goto error;
2936
2937 isl_seq_cpy(bmap->div[k], term->div->row[pos], 2 + total);
2938
2939 return isl_basic_map_div(bmap, k);
2940 error:
2941 isl_basic_map_free(bmap);
2942 return NULL;
2943 }
2944
2945 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
2946 int (*fn)(__isl_take isl_term *term, void *user),
2947 __isl_take isl_term *term, void *user)
2948 {
2949 int i;
2950 struct isl_upoly_rec *rec;
2951
2952 if (!up || !term)
2953 goto error;
2954
2955 if (isl_upoly_is_zero(up))
2956 return term;
2957
2958 isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
2959 isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
2960 isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
2961
2962 if (isl_upoly_is_cst(up)) {
2963 struct isl_upoly_cst *cst;
2964 cst = isl_upoly_as_cst(up);
2965 if (!cst)
2966 goto error;
2967 term = isl_term_cow(term);
2968 if (!term)
2969 goto error;
2970 isl_int_set(term->n, cst->n);
2971 isl_int_set(term->d, cst->d);
2972 if (fn(isl_term_copy(term), user) < 0)
2973 goto error;
2974 return term;
2975 }
2976
2977 rec = isl_upoly_as_rec(up);
2978 if (!rec)
2979 goto error;
2980
2981 for (i = 0; i < rec->n; ++i) {
2982 term = isl_term_cow(term);
2983 if (!term)
2984 goto error;
2985 term->pow[up->var] = i;
2986 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
2987 if (!term)
2988 goto error;
2989 }
2990 term->pow[up->var] = 0;
2991
2992 return term;
2993 error:
2994 isl_term_free(term);
2995 return NULL;
2996 }
2997
2998 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
2999 int (*fn)(__isl_take isl_term *term, void *user), void *user)
3000 {
3001 isl_term *term;
3002
3003 if (!qp)
3004 return -1;
3005
3006 term = isl_term_alloc(isl_dim_copy(qp->dim), isl_mat_copy(qp->div));
3007 if (!term)
3008 return -1;
3009
3010 term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
3011
3012 isl_term_free(term);
3013
3014 return term ? 0 : -1;
3015 }
3016
3017 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
3018 {
3019 struct isl_upoly *up;
3020 isl_qpolynomial *qp;
3021 int i, n;
3022
3023 if (!term)
3024 return NULL;
3025
3026 n = isl_dim_total(term->dim) + term->div->n_row;
3027
3028 up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
3029 for (i = 0; i < n; ++i) {
3030 if (!term->pow[i])
3031 continue;
3032 up = isl_upoly_mul(up,
3033 isl_upoly_pow(term->dim->ctx, i, term->pow[i]));
3034 }
3035
3036 qp = isl_qpolynomial_alloc(isl_dim_copy(term->dim), term->div->n_row, up);
3037 if (!qp)
3038 goto error;
3039 isl_mat_free(qp->div);
3040 qp->div = isl_mat_copy(term->div);
3041 if (!qp->div)
3042 goto error;
3043
3044 isl_term_free(term);
3045 return qp;
3046 error:
3047 isl_qpolynomial_free(qp);
3048 isl_term_free(term);
3049 return NULL;
3050 }
3051
3052 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
3053 __isl_take isl_dim *dim)
3054 {
3055 int i;
3056 int extra;
3057 unsigned total;
3058
3059 if (!qp || !dim)
3060 goto error;
3061
3062 if (isl_dim_equal(qp->dim, dim)) {
3063 isl_dim_free(dim);
3064 return qp;
3065 }
3066
3067 qp = isl_qpolynomial_cow(qp);
3068 if (!qp)
3069 goto error;
3070
3071 extra = isl_dim_size(dim, isl_dim_set) -
3072 isl_dim_size(qp->dim, isl_dim_set);
3073 total = isl_dim_total(qp->dim);
3074 if (qp->div->n_row) {
3075 int *exp;
3076
3077 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
3078 if (!exp)
3079 goto error;
3080 for (i = 0; i < qp->div->n_row; ++i)
3081 exp[i] = extra + i;
3082 qp->upoly = expand(qp->upoly, exp, total);
3083 free(exp);
3084 if (!qp->upoly)
3085 goto error;
3086 }
3087 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
3088 if (!qp->div)
3089 goto error;
3090 for (i = 0; i < qp->div->n_row; ++i)
3091 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
3092
3093 isl_dim_free(qp->dim);
3094 qp->dim = dim;
3095
3096 return qp;
3097 error:
3098 isl_dim_free(dim);
3099 isl_qpolynomial_free(qp);
3100 return NULL;
3101 }
3102
3103 /* For each parameter or variable that does not appear in qp,
3104 * first eliminate the variable from all constraints and then set it to zero.
3105 */
3106 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
3107 __isl_keep isl_qpolynomial *qp)
3108 {
3109 int *active = NULL;
3110 int i;
3111 int d;
3112 unsigned nparam;
3113 unsigned nvar;
3114
3115 if (!set || !qp)
3116 goto error;
3117
3118 d = isl_dim_total(set->dim);
3119 active = isl_calloc_array(set->ctx, int, d);
3120 if (set_active(qp, active) < 0)
3121 goto error;
3122
3123 for (i = 0; i < d; ++i)
3124 if (!active[i])
3125 break;
3126
3127 if (i == d) {
3128 free(active);
3129 return set;
3130 }
3131
3132 nparam = isl_dim_size(set->dim, isl_dim_param);
3133 nvar = isl_dim_size(set->dim, isl_dim_set);
3134 for (i = 0; i < nparam; ++i) {
3135 if (active[i])
3136 continue;
3137 set = isl_set_eliminate(set, isl_dim_param, i, 1);
3138 set = isl_set_fix_si(set, isl_dim_param, i, 0);
3139 }
3140 for (i = 0; i < nvar; ++i) {
3141 if (active[nparam + i])
3142 continue;
3143 set = isl_set_eliminate(set, isl_dim_set, i, 1);
3144 set = isl_set_fix_si(set, isl_dim_set, i, 0);
3145 }
3146
3147 free(active);
3148
3149 return set;
3150 error:
3151 free(active);
3152 isl_set_free(set);
3153 return NULL;
3154 }
3155
3156 struct isl_opt_data {
3157 isl_qpolynomial *qp;
3158 int first;
3159 isl_qpolynomial *opt;
3160 int max;
3161 };
3162
3163 static int opt_fn(__isl_take isl_point *pnt, void *user)
3164 {
3165 struct isl_opt_data *data = (struct isl_opt_data *)user;
3166 isl_qpolynomial *val;
3167
3168 val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
3169 if (data->first) {
3170 data->first = 0;
3171 data->opt = val;
3172 } else if (data->max) {
3173 data->opt = isl_qpolynomial_max_cst(data->opt, val);
3174 } else {
3175 data->opt = isl_qpolynomial_min_cst(data->opt, val);
3176 }
3177
3178 return 0;
3179 }
3180
3181 __isl_give isl_qpolynomial *isl_qpolynomial_opt_on_domain(
3182 __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
3183 {
3184 struct isl_opt_data data = { NULL, 1, NULL, max };
3185
3186 if (!set || !qp)
3187 goto error;
3188
3189 if (isl_upoly_is_cst(qp->upoly)) {
3190 isl_set_free(set);
3191 return qp;
3192 }
3193
3194 set = fix_inactive(set, qp);
3195
3196 data.qp = qp;
3197 if (isl_set_foreach_point(set, opt_fn, &data) < 0)
3198 goto error;
3199
3200 if (data.first)
3201 data.opt = isl_qpolynomial_zero(isl_qpolynomial_get_dim(qp));
3202
3203 isl_set_free(set);
3204 isl_qpolynomial_free(qp);
3205 return data.opt;
3206 error:
3207 isl_set_free(set);
3208 isl_qpolynomial_free(qp);
3209 isl_qpolynomial_free(data.opt);
3210 return NULL;
3211 }
3212
3213 __isl_give isl_qpolynomial *isl_qpolynomial_morph(__isl_take isl_qpolynomial *qp,
3214 __isl_take isl_morph *morph)
3215 {
3216 int i;
3217 isl_ctx *ctx;
3218 struct isl_upoly *up;
3219 unsigned n_div;
3220 struct isl_upoly **subs;
3221 isl_mat *mat;
3222
3223 qp = isl_qpolynomial_cow(qp);
3224 if (!qp || !morph)
3225 goto error;
3226
3227 ctx = qp->dim->ctx;
3228 isl_assert(ctx, isl_dim_equal(qp->dim, morph->dom->dim), goto error);
3229
3230 subs = isl_calloc_array(ctx, struct isl_upoly *, morph->inv->n_row - 1);
3231 if (!subs)
3232 goto error;
3233
3234 for (i = 0; 1 + i < morph->inv->n_row; ++i)
3235 subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
3236 morph->inv->row[0][0], morph->inv->n_col);
3237
3238 qp->upoly = isl_upoly_subs(qp->upoly, 0, morph->inv->n_row - 1, subs);
3239
3240 for (i = 0; 1 + i < morph->inv->n_row; ++i)
3241 isl_upoly_free(subs[i]);
3242 free(subs);
3243
3244 mat = isl_mat_diagonal(isl_mat_identity(ctx, 1), isl_mat_copy(morph->inv));
3245 mat = isl_mat_diagonal(mat, isl_mat_identity(ctx, qp->div->n_row));
3246 qp->div = isl_mat_product(qp->div, mat);
3247 isl_dim_free(qp->dim);
3248 qp->dim = isl_dim_copy(morph->ran->dim);
3249
3250 if (!qp->upoly || !qp->div || !qp->dim)
3251 goto error;
3252
3253 isl_morph_free(morph);
3254
3255 return qp;
3256 error:
3257 isl_qpolynomial_free(qp);
3258 isl_morph_free(morph);
3259 return NULL;
3260 }
3261
3262 static int neg_entry(void **entry, void *user)
3263 {
3264 isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
3265
3266 *pwqp = isl_pw_qpolynomial_neg(*pwqp);
3267
3268 return *pwqp ? 0 : -1;
3269 }
3270
3271 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_neg(
3272 __isl_take isl_union_pw_qpolynomial *upwqp)
3273 {
3274 upwqp = isl_union_pw_qpolynomial_cow(upwqp);
3275 if (!upwqp)
3276 return NULL;
3277
3278 if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
3279 &neg_entry, NULL) < 0)
3280 goto error;
3281
3282 return upwqp;
3283 error:
3284 isl_union_pw_qpolynomial_free(upwqp);
3285 return NULL;
3286 }
3287
3288 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_sub(
3289 __isl_take isl_union_pw_qpolynomial *upwqp1,
3290 __isl_take isl_union_pw_qpolynomial *upwqp2)
3291 {
3292 return isl_union_pw_qpolynomial_add(upwqp1,
3293 isl_union_pw_qpolynomial_neg(upwqp2));
3294 }
3295
3296 static int mul_entry(void **entry, void *user)
3297 {
3298 struct isl_union_pw_qpolynomial_match_bin_data *data = user;
3299 uint32_t hash;
3300 struct isl_hash_table_entry *entry2;
3301 isl_pw_qpolynomial *pwpq = *entry;
3302 int empty;
3303
3304 hash = isl_dim_get_hash(pwpq->dim);
3305 entry2 = isl_hash_table_find(data->upwqp2->dim->ctx,
3306 &data->upwqp2->table,
3307 hash, &has_dim, pwpq->dim, 0);
3308 if (!entry2)
3309 return 0;
3310
3311 pwpq = isl_pw_qpolynomial_copy(pwpq);
3312 pwpq = isl_pw_qpolynomial_mul(pwpq,
3313 isl_pw_qpolynomial_copy(entry2->data));
3314
3315 empty = isl_pw_qpolynomial_is_zero(pwpq);
3316 if (empty < 0) {
3317 isl_pw_qpolynomial_free(pwpq);
3318 return -1;
3319 }
3320 if (empty) {
3321 isl_pw_qpolynomial_free(pwpq);
3322 return 0;
3323 }
3324
3325 data->res = isl_union_pw_qpolynomial_add_pw_qpolynomial(data->res, pwpq);
3326
3327 return 0;
3328 }
3329
3330 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
3331 __isl_take isl_union_pw_qpolynomial *upwqp1,
3332 __isl_take isl_union_pw_qpolynomial *upwqp2)
3333 {
3334 return match_bin_op(upwqp1, upwqp2, &mul_entry);
3335 }
Integer set library