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