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