?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
isl_*_eval: return an isl_val instead of an isl_qpolynomial
[isl.git] / isl_polynomial.c
bloba69b654a39ad05c0275252a93a9d7085ada5ebe9
1 /*
2 * Copyright 2010 INRIA Saclay
3 *
4 * Use of this software is governed by the MIT 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 #define ISL_DIM_H
13 #include <isl_ctx_private.h>
14 #include <isl_map_private.h>
15 #include <isl_factorization.h>
16 #include <isl_lp_private.h>
17 #include <isl_seq.h>
18 #include <isl_union_map_private.h>
19 #include <isl_constraint_private.h>
20 #include <isl_polynomial_private.h>
21 #include <isl_point_private.h>
22 #include <isl_space_private.h>
23 #include <isl_mat_private.h>
24 #include <isl_vec_private.h>
25 #include <isl_range.h>
26 #include <isl_local_space_private.h>
27 #include <isl_aff_private.h>
28 #include <isl_val_private.h>
29 #include <isl_config.h>
30 #include <isl/deprecated/polynomial_int.h>
31
32 static unsigned pos(__isl_keep isl_space *dim, enum isl_dim_type type)
33 {
34 switch (type) {
35 case isl_dim_param: return 0;
36 case isl_dim_in: return dim->nparam;
37 case isl_dim_out: return dim->nparam + dim->n_in;
38 default: return 0;
39 }
40 }
41
42 int isl_upoly_is_cst(__isl_keep struct isl_upoly *up)
43 {
44 if (!up)
45 return -1;
46
47 return up->var < 0;
48 }
49
50 __isl_keep struct isl_upoly_cst *isl_upoly_as_cst(__isl_keep struct isl_upoly *up)
51 {
52 if (!up)
53 return NULL;
54
55 isl_assert(up->ctx, up->var < 0, return NULL);
56
57 return (struct isl_upoly_cst *)up;
58 }
59
60 __isl_keep struct isl_upoly_rec *isl_upoly_as_rec(__isl_keep struct isl_upoly *up)
61 {
62 if (!up)
63 return NULL;
64
65 isl_assert(up->ctx, up->var >= 0, return NULL);
66
67 return (struct isl_upoly_rec *)up;
68 }
69
70 int isl_upoly_is_equal(__isl_keep struct isl_upoly *up1,
71 __isl_keep struct isl_upoly *up2)
72 {
73 int i;
74 struct isl_upoly_rec *rec1, *rec2;
75
76 if (!up1 || !up2)
77 return -1;
78 if (up1 == up2)
79 return 1;
80 if (up1->var != up2->var)
81 return 0;
82 if (isl_upoly_is_cst(up1)) {
83 struct isl_upoly_cst *cst1, *cst2;
84 cst1 = isl_upoly_as_cst(up1);
85 cst2 = isl_upoly_as_cst(up2);
86 if (!cst1 || !cst2)
87 return -1;
88 return isl_int_eq(cst1->n, cst2->n) &&
89 isl_int_eq(cst1->d, cst2->d);
90 }
91
92 rec1 = isl_upoly_as_rec(up1);
93 rec2 = isl_upoly_as_rec(up2);
94 if (!rec1 || !rec2)
95 return -1;
96
97 if (rec1->n != rec2->n)
98 return 0;
99
100 for (i = 0; i < rec1->n; ++i) {
101 int eq = isl_upoly_is_equal(rec1->p[i], rec2->p[i]);
102 if (eq < 0 || !eq)
103 return eq;
104 }
105
106 return 1;
107 }
108
109 int isl_upoly_is_zero(__isl_keep struct isl_upoly *up)
110 {
111 struct isl_upoly_cst *cst;
112
113 if (!up)
114 return -1;
115 if (!isl_upoly_is_cst(up))
116 return 0;
117
118 cst = isl_upoly_as_cst(up);
119 if (!cst)
120 return -1;
121
122 return isl_int_is_zero(cst->n) && isl_int_is_pos(cst->d);
123 }
124
125 int isl_upoly_sgn(__isl_keep struct isl_upoly *up)
126 {
127 struct isl_upoly_cst *cst;
128
129 if (!up)
130 return 0;
131 if (!isl_upoly_is_cst(up))
132 return 0;
133
134 cst = isl_upoly_as_cst(up);
135 if (!cst)
136 return 0;
137
138 return isl_int_sgn(cst->n);
139 }
140
141 int isl_upoly_is_nan(__isl_keep struct isl_upoly *up)
142 {
143 struct isl_upoly_cst *cst;
144
145 if (!up)
146 return -1;
147 if (!isl_upoly_is_cst(up))
148 return 0;
149
150 cst = isl_upoly_as_cst(up);
151 if (!cst)
152 return -1;
153
154 return isl_int_is_zero(cst->n) && isl_int_is_zero(cst->d);
155 }
156
157 int isl_upoly_is_infty(__isl_keep struct isl_upoly *up)
158 {
159 struct isl_upoly_cst *cst;
160
161 if (!up)
162 return -1;
163 if (!isl_upoly_is_cst(up))
164 return 0;
165
166 cst = isl_upoly_as_cst(up);
167 if (!cst)
168 return -1;
169
170 return isl_int_is_pos(cst->n) && isl_int_is_zero(cst->d);
171 }
172
173 int isl_upoly_is_neginfty(__isl_keep struct isl_upoly *up)
174 {
175 struct isl_upoly_cst *cst;
176
177 if (!up)
178 return -1;
179 if (!isl_upoly_is_cst(up))
180 return 0;
181
182 cst = isl_upoly_as_cst(up);
183 if (!cst)
184 return -1;
185
186 return isl_int_is_neg(cst->n) && isl_int_is_zero(cst->d);
187 }
188
189 int isl_upoly_is_one(__isl_keep struct isl_upoly *up)
190 {
191 struct isl_upoly_cst *cst;
192
193 if (!up)
194 return -1;
195 if (!isl_upoly_is_cst(up))
196 return 0;
197
198 cst = isl_upoly_as_cst(up);
199 if (!cst)
200 return -1;
201
202 return isl_int_eq(cst->n, cst->d) && isl_int_is_pos(cst->d);
203 }
204
205 int isl_upoly_is_negone(__isl_keep struct isl_upoly *up)
206 {
207 struct isl_upoly_cst *cst;
208
209 if (!up)
210 return -1;
211 if (!isl_upoly_is_cst(up))
212 return 0;
213
214 cst = isl_upoly_as_cst(up);
215 if (!cst)
216 return -1;
217
218 return isl_int_is_negone(cst->n) && isl_int_is_one(cst->d);
219 }
220
221 __isl_give struct isl_upoly_cst *isl_upoly_cst_alloc(struct isl_ctx *ctx)
222 {
223 struct isl_upoly_cst *cst;
224
225 cst = isl_alloc_type(ctx, struct isl_upoly_cst);
226 if (!cst)
227 return NULL;
228
229 cst->up.ref = 1;
230 cst->up.ctx = ctx;
231 isl_ctx_ref(ctx);
232 cst->up.var = -1;
233
234 isl_int_init(cst->n);
235 isl_int_init(cst->d);
236
237 return cst;
238 }
239
240 __isl_give struct isl_upoly *isl_upoly_zero(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, 0);
249 isl_int_set_si(cst->d, 1);
250
251 return &cst->up;
252 }
253
254 __isl_give struct isl_upoly *isl_upoly_one(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, 1);
264
265 return &cst->up;
266 }
267
268 __isl_give struct isl_upoly *isl_upoly_infty(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, 1);
277 isl_int_set_si(cst->d, 0);
278
279 return &cst->up;
280 }
281
282 __isl_give struct isl_upoly *isl_upoly_neginfty(struct isl_ctx *ctx)
283 {
284 struct isl_upoly_cst *cst;
285
286 cst = isl_upoly_cst_alloc(ctx);
287 if (!cst)
288 return NULL;
289
290 isl_int_set_si(cst->n, -1);
291 isl_int_set_si(cst->d, 0);
292
293 return &cst->up;
294 }
295
296 __isl_give struct isl_upoly *isl_upoly_nan(struct isl_ctx *ctx)
297 {
298 struct isl_upoly_cst *cst;
299
300 cst = isl_upoly_cst_alloc(ctx);
301 if (!cst)
302 return NULL;
303
304 isl_int_set_si(cst->n, 0);
305 isl_int_set_si(cst->d, 0);
306
307 return &cst->up;
308 }
309
310 __isl_give struct isl_upoly *isl_upoly_rat_cst(struct isl_ctx *ctx,
311 isl_int n, isl_int d)
312 {
313 struct isl_upoly_cst *cst;
314
315 cst = isl_upoly_cst_alloc(ctx);
316 if (!cst)
317 return NULL;
318
319 isl_int_set(cst->n, n);
320 isl_int_set(cst->d, d);
321
322 return &cst->up;
323 }
324
325 __isl_give struct isl_upoly_rec *isl_upoly_alloc_rec(struct isl_ctx *ctx,
326 int var, int size)
327 {
328 struct isl_upoly_rec *rec;
329
330 isl_assert(ctx, var >= 0, return NULL);
331 isl_assert(ctx, size >= 0, return NULL);
332 rec = isl_calloc(ctx, struct isl_upoly_rec,
333 sizeof(struct isl_upoly_rec) +
334 size * sizeof(struct isl_upoly *));
335 if (!rec)
336 return NULL;
337
338 rec->up.ref = 1;
339 rec->up.ctx = ctx;
340 isl_ctx_ref(ctx);
341 rec->up.var = var;
342
343 rec->n = 0;
344 rec->size = size;
345
346 return rec;
347 }
348
349 __isl_give isl_qpolynomial *isl_qpolynomial_reset_domain_space(
350 __isl_take isl_qpolynomial *qp, __isl_take isl_space *dim)
351 {
352 qp = isl_qpolynomial_cow(qp);
353 if (!qp || !dim)
354 goto error;
355
356 isl_space_free(qp->dim);
357 qp->dim = dim;
358
359 return qp;
360 error:
361 isl_qpolynomial_free(qp);
362 isl_space_free(dim);
363 return NULL;
364 }
365
366 /* Reset the space of "qp". This function is called from isl_pw_templ.c
367 * and doesn't know if the space of an element object is represented
368 * directly or through its domain. It therefore passes along both.
369 */
370 __isl_give isl_qpolynomial *isl_qpolynomial_reset_space_and_domain(
371 __isl_take isl_qpolynomial *qp, __isl_take isl_space *space,
372 __isl_take isl_space *domain)
373 {
374 isl_space_free(space);
375 return isl_qpolynomial_reset_domain_space(qp, domain);
376 }
377
378 isl_ctx *isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial *qp)
379 {
380 return qp ? qp->dim->ctx : NULL;
381 }
382
383 __isl_give isl_space *isl_qpolynomial_get_domain_space(
384 __isl_keep isl_qpolynomial *qp)
385 {
386 return qp ? isl_space_copy(qp->dim) : NULL;
387 }
388
389 __isl_give isl_space *isl_qpolynomial_get_space(__isl_keep isl_qpolynomial *qp)
390 {
391 isl_space *space;
392 if (!qp)
393 return NULL;
394 space = isl_space_copy(qp->dim);
395 space = isl_space_from_domain(space);
396 space = isl_space_add_dims(space, isl_dim_out, 1);
397 return space;
398 }
399
400 /* Externally, an isl_qpolynomial has a map space, but internally, the
401 * ls field corresponds to the domain of that space.
402 */
403 unsigned isl_qpolynomial_dim(__isl_keep isl_qpolynomial *qp,
404 enum isl_dim_type type)
405 {
406 if (!qp)
407 return 0;
408 if (type == isl_dim_out)
409 return 1;
410 if (type == isl_dim_in)
411 type = isl_dim_set;
412 return isl_space_dim(qp->dim, type);
413 }
414
415 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial *qp)
416 {
417 return qp ? isl_upoly_is_zero(qp->upoly) : -1;
418 }
419
420 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial *qp)
421 {
422 return qp ? isl_upoly_is_one(qp->upoly) : -1;
423 }
424
425 int isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial *qp)
426 {
427 return qp ? isl_upoly_is_nan(qp->upoly) : -1;
428 }
429
430 int isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial *qp)
431 {
432 return qp ? isl_upoly_is_infty(qp->upoly) : -1;
433 }
434
435 int isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial *qp)
436 {
437 return qp ? isl_upoly_is_neginfty(qp->upoly) : -1;
438 }
439
440 int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial *qp)
441 {
442 return qp ? isl_upoly_sgn(qp->upoly) : 0;
443 }
444
445 static void upoly_free_cst(__isl_take struct isl_upoly_cst *cst)
446 {
447 isl_int_clear(cst->n);
448 isl_int_clear(cst->d);
449 }
450
451 static void upoly_free_rec(__isl_take struct isl_upoly_rec *rec)
452 {
453 int i;
454
455 for (i = 0; i < rec->n; ++i)
456 isl_upoly_free(rec->p[i]);
457 }
458
459 __isl_give struct isl_upoly *isl_upoly_copy(__isl_keep struct isl_upoly *up)
460 {
461 if (!up)
462 return NULL;
463
464 up->ref++;
465 return up;
466 }
467
468 __isl_give struct isl_upoly *isl_upoly_dup_cst(__isl_keep struct isl_upoly *up)
469 {
470 struct isl_upoly_cst *cst;
471 struct isl_upoly_cst *dup;
472
473 cst = isl_upoly_as_cst(up);
474 if (!cst)
475 return NULL;
476
477 dup = isl_upoly_as_cst(isl_upoly_zero(up->ctx));
478 if (!dup)
479 return NULL;
480 isl_int_set(dup->n, cst->n);
481 isl_int_set(dup->d, cst->d);
482
483 return &dup->up;
484 }
485
486 __isl_give struct isl_upoly *isl_upoly_dup_rec(__isl_keep struct isl_upoly *up)
487 {
488 int i;
489 struct isl_upoly_rec *rec;
490 struct isl_upoly_rec *dup;
491
492 rec = isl_upoly_as_rec(up);
493 if (!rec)
494 return NULL;
495
496 dup = isl_upoly_alloc_rec(up->ctx, up->var, rec->n);
497 if (!dup)
498 return NULL;
499
500 for (i = 0; i < rec->n; ++i) {
501 dup->p[i] = isl_upoly_copy(rec->p[i]);
502 if (!dup->p[i])
503 goto error;
504 dup->n++;
505 }
506
507 return &dup->up;
508 error:
509 isl_upoly_free(&dup->up);
510 return NULL;
511 }
512
513 __isl_give struct isl_upoly *isl_upoly_dup(__isl_keep struct isl_upoly *up)
514 {
515 if (!up)
516 return NULL;
517
518 if (isl_upoly_is_cst(up))
519 return isl_upoly_dup_cst(up);
520 else
521 return isl_upoly_dup_rec(up);
522 }
523
524 __isl_give struct isl_upoly *isl_upoly_cow(__isl_take struct isl_upoly *up)
525 {
526 if (!up)
527 return NULL;
528
529 if (up->ref == 1)
530 return up;
531 up->ref--;
532 return isl_upoly_dup(up);
533 }
534
535 void isl_upoly_free(__isl_take struct isl_upoly *up)
536 {
537 if (!up)
538 return;
539
540 if (--up->ref > 0)
541 return;
542
543 if (up->var < 0)
544 upoly_free_cst((struct isl_upoly_cst *)up);
545 else
546 upoly_free_rec((struct isl_upoly_rec *)up);
547
548 isl_ctx_deref(up->ctx);
549 free(up);
550 }
551
552 static void isl_upoly_cst_reduce(__isl_keep struct isl_upoly_cst *cst)
553 {
554 isl_int gcd;
555
556 isl_int_init(gcd);
557 isl_int_gcd(gcd, cst->n, cst->d);
558 if (!isl_int_is_zero(gcd) && !isl_int_is_one(gcd)) {
559 isl_int_divexact(cst->n, cst->n, gcd);
560 isl_int_divexact(cst->d, cst->d, gcd);
561 }
562 isl_int_clear(gcd);
563 }
564
565 __isl_give struct isl_upoly *isl_upoly_sum_cst(__isl_take struct isl_upoly *up1,
566 __isl_take struct isl_upoly *up2)
567 {
568 struct isl_upoly_cst *cst1;
569 struct isl_upoly_cst *cst2;
570
571 up1 = isl_upoly_cow(up1);
572 if (!up1 || !up2)
573 goto error;
574
575 cst1 = isl_upoly_as_cst(up1);
576 cst2 = isl_upoly_as_cst(up2);
577
578 if (isl_int_eq(cst1->d, cst2->d))
579 isl_int_add(cst1->n, cst1->n, cst2->n);
580 else {
581 isl_int_mul(cst1->n, cst1->n, cst2->d);
582 isl_int_addmul(cst1->n, cst2->n, cst1->d);
583 isl_int_mul(cst1->d, cst1->d, cst2->d);
584 }
585
586 isl_upoly_cst_reduce(cst1);
587
588 isl_upoly_free(up2);
589 return up1;
590 error:
591 isl_upoly_free(up1);
592 isl_upoly_free(up2);
593 return NULL;
594 }
595
596 static __isl_give struct isl_upoly *replace_by_zero(
597 __isl_take struct isl_upoly *up)
598 {
599 struct isl_ctx *ctx;
600
601 if (!up)
602 return NULL;
603 ctx = up->ctx;
604 isl_upoly_free(up);
605 return isl_upoly_zero(ctx);
606 }
607
608 static __isl_give struct isl_upoly *replace_by_constant_term(
609 __isl_take struct isl_upoly *up)
610 {
611 struct isl_upoly_rec *rec;
612 struct isl_upoly *cst;
613
614 if (!up)
615 return NULL;
616
617 rec = isl_upoly_as_rec(up);
618 if (!rec)
619 goto error;
620 cst = isl_upoly_copy(rec->p[0]);
621 isl_upoly_free(up);
622 return cst;
623 error:
624 isl_upoly_free(up);
625 return NULL;
626 }
627
628 __isl_give struct isl_upoly *isl_upoly_sum(__isl_take struct isl_upoly *up1,
629 __isl_take struct isl_upoly *up2)
630 {
631 int i;
632 struct isl_upoly_rec *rec1, *rec2;
633
634 if (!up1 || !up2)
635 goto error;
636
637 if (isl_upoly_is_nan(up1)) {
638 isl_upoly_free(up2);
639 return up1;
640 }
641
642 if (isl_upoly_is_nan(up2)) {
643 isl_upoly_free(up1);
644 return up2;
645 }
646
647 if (isl_upoly_is_zero(up1)) {
648 isl_upoly_free(up1);
649 return up2;
650 }
651
652 if (isl_upoly_is_zero(up2)) {
653 isl_upoly_free(up2);
654 return up1;
655 }
656
657 if (up1->var < up2->var)
658 return isl_upoly_sum(up2, up1);
659
660 if (up2->var < up1->var) {
661 struct isl_upoly_rec *rec;
662 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
663 isl_upoly_free(up1);
664 return up2;
665 }
666 up1 = isl_upoly_cow(up1);
667 rec = isl_upoly_as_rec(up1);
668 if (!rec)
669 goto error;
670 rec->p[0] = isl_upoly_sum(rec->p[0], up2);
671 if (rec->n == 1)
672 up1 = replace_by_constant_term(up1);
673 return up1;
674 }
675
676 if (isl_upoly_is_cst(up1))
677 return isl_upoly_sum_cst(up1, up2);
678
679 rec1 = isl_upoly_as_rec(up1);
680 rec2 = isl_upoly_as_rec(up2);
681 if (!rec1 || !rec2)
682 goto error;
683
684 if (rec1->n < rec2->n)
685 return isl_upoly_sum(up2, up1);
686
687 up1 = isl_upoly_cow(up1);
688 rec1 = isl_upoly_as_rec(up1);
689 if (!rec1)
690 goto error;
691
692 for (i = rec2->n - 1; i >= 0; --i) {
693 rec1->p[i] = isl_upoly_sum(rec1->p[i],
694 isl_upoly_copy(rec2->p[i]));
695 if (!rec1->p[i])
696 goto error;
697 if (i == rec1->n - 1 && isl_upoly_is_zero(rec1->p[i])) {
698 isl_upoly_free(rec1->p[i]);
699 rec1->n--;
700 }
701 }
702
703 if (rec1->n == 0)
704 up1 = replace_by_zero(up1);
705 else if (rec1->n == 1)
706 up1 = replace_by_constant_term(up1);
707
708 isl_upoly_free(up2);
709
710 return up1;
711 error:
712 isl_upoly_free(up1);
713 isl_upoly_free(up2);
714 return NULL;
715 }
716
717 __isl_give struct isl_upoly *isl_upoly_cst_add_isl_int(
718 __isl_take struct isl_upoly *up, isl_int v)
719 {
720 struct isl_upoly_cst *cst;
721
722 up = isl_upoly_cow(up);
723 if (!up)
724 return NULL;
725
726 cst = isl_upoly_as_cst(up);
727
728 isl_int_addmul(cst->n, cst->d, v);
729
730 return up;
731 }
732
733 __isl_give struct isl_upoly *isl_upoly_add_isl_int(
734 __isl_take struct isl_upoly *up, isl_int v)
735 {
736 struct isl_upoly_rec *rec;
737
738 if (!up)
739 return NULL;
740
741 if (isl_upoly_is_cst(up))
742 return isl_upoly_cst_add_isl_int(up, v);
743
744 up = isl_upoly_cow(up);
745 rec = isl_upoly_as_rec(up);
746 if (!rec)
747 goto error;
748
749 rec->p[0] = isl_upoly_add_isl_int(rec->p[0], v);
750 if (!rec->p[0])
751 goto error;
752
753 return up;
754 error:
755 isl_upoly_free(up);
756 return NULL;
757 }
758
759 __isl_give struct isl_upoly *isl_upoly_cst_mul_isl_int(
760 __isl_take struct isl_upoly *up, isl_int v)
761 {
762 struct isl_upoly_cst *cst;
763
764 if (isl_upoly_is_zero(up))
765 return up;
766
767 up = isl_upoly_cow(up);
768 if (!up)
769 return NULL;
770
771 cst = isl_upoly_as_cst(up);
772
773 isl_int_mul(cst->n, cst->n, v);
774
775 return up;
776 }
777
778 __isl_give struct isl_upoly *isl_upoly_mul_isl_int(
779 __isl_take struct isl_upoly *up, isl_int v)
780 {
781 int i;
782 struct isl_upoly_rec *rec;
783
784 if (!up)
785 return NULL;
786
787 if (isl_upoly_is_cst(up))
788 return isl_upoly_cst_mul_isl_int(up, v);
789
790 up = isl_upoly_cow(up);
791 rec = isl_upoly_as_rec(up);
792 if (!rec)
793 goto error;
794
795 for (i = 0; i < rec->n; ++i) {
796 rec->p[i] = isl_upoly_mul_isl_int(rec->p[i], v);
797 if (!rec->p[i])
798 goto error;
799 }
800
801 return up;
802 error:
803 isl_upoly_free(up);
804 return NULL;
805 }
806
807 /* Multiply the constant polynomial "up" by "v".
808 */
809 static __isl_give struct isl_upoly *isl_upoly_cst_scale_val(
810 __isl_take struct isl_upoly *up, __isl_keep isl_val *v)
811 {
812 struct isl_upoly_cst *cst;
813
814 if (isl_upoly_is_zero(up))
815 return up;
816
817 up = isl_upoly_cow(up);
818 if (!up)
819 return NULL;
820
821 cst = isl_upoly_as_cst(up);
822
823 isl_int_mul(cst->n, cst->n, v->n);
824 isl_int_mul(cst->d, cst->d, v->d);
825 isl_upoly_cst_reduce(cst);
826
827 return up;
828 }
829
830 /* Multiply the polynomial "up" by "v".
831 */
832 static __isl_give struct isl_upoly *isl_upoly_scale_val(
833 __isl_take struct isl_upoly *up, __isl_keep isl_val *v)
834 {
835 int i;
836 struct isl_upoly_rec *rec;
837
838 if (!up)
839 return NULL;
840
841 if (isl_upoly_is_cst(up))
842 return isl_upoly_cst_scale_val(up, v);
843
844 up = isl_upoly_cow(up);
845 rec = isl_upoly_as_rec(up);
846 if (!rec)
847 goto error;
848
849 for (i = 0; i < rec->n; ++i) {
850 rec->p[i] = isl_upoly_scale_val(rec->p[i], v);
851 if (!rec->p[i])
852 goto error;
853 }
854
855 return up;
856 error:
857 isl_upoly_free(up);
858 return NULL;
859 }
860
861 __isl_give struct isl_upoly *isl_upoly_mul_cst(__isl_take struct isl_upoly *up1,
862 __isl_take struct isl_upoly *up2)
863 {
864 struct isl_upoly_cst *cst1;
865 struct isl_upoly_cst *cst2;
866
867 up1 = isl_upoly_cow(up1);
868 if (!up1 || !up2)
869 goto error;
870
871 cst1 = isl_upoly_as_cst(up1);
872 cst2 = isl_upoly_as_cst(up2);
873
874 isl_int_mul(cst1->n, cst1->n, cst2->n);
875 isl_int_mul(cst1->d, cst1->d, cst2->d);
876
877 isl_upoly_cst_reduce(cst1);
878
879 isl_upoly_free(up2);
880 return up1;
881 error:
882 isl_upoly_free(up1);
883 isl_upoly_free(up2);
884 return NULL;
885 }
886
887 __isl_give struct isl_upoly *isl_upoly_mul_rec(__isl_take struct isl_upoly *up1,
888 __isl_take struct isl_upoly *up2)
889 {
890 struct isl_upoly_rec *rec1;
891 struct isl_upoly_rec *rec2;
892 struct isl_upoly_rec *res = NULL;
893 int i, j;
894 int size;
895
896 rec1 = isl_upoly_as_rec(up1);
897 rec2 = isl_upoly_as_rec(up2);
898 if (!rec1 || !rec2)
899 goto error;
900 size = rec1->n + rec2->n - 1;
901 res = isl_upoly_alloc_rec(up1->ctx, up1->var, size);
902 if (!res)
903 goto error;
904
905 for (i = 0; i < rec1->n; ++i) {
906 res->p[i] = isl_upoly_mul(isl_upoly_copy(rec2->p[0]),
907 isl_upoly_copy(rec1->p[i]));
908 if (!res->p[i])
909 goto error;
910 res->n++;
911 }
912 for (; i < size; ++i) {
913 res->p[i] = isl_upoly_zero(up1->ctx);
914 if (!res->p[i])
915 goto error;
916 res->n++;
917 }
918 for (i = 0; i < rec1->n; ++i) {
919 for (j = 1; j < rec2->n; ++j) {
920 struct isl_upoly *up;
921 up = isl_upoly_mul(isl_upoly_copy(rec2->p[j]),
922 isl_upoly_copy(rec1->p[i]));
923 res->p[i + j] = isl_upoly_sum(res->p[i + j], up);
924 if (!res->p[i + j])
925 goto error;
926 }
927 }
928
929 isl_upoly_free(up1);
930 isl_upoly_free(up2);
931
932 return &res->up;
933 error:
934 isl_upoly_free(up1);
935 isl_upoly_free(up2);
936 isl_upoly_free(&res->up);
937 return NULL;
938 }
939
940 __isl_give struct isl_upoly *isl_upoly_mul(__isl_take struct isl_upoly *up1,
941 __isl_take struct isl_upoly *up2)
942 {
943 if (!up1 || !up2)
944 goto error;
945
946 if (isl_upoly_is_nan(up1)) {
947 isl_upoly_free(up2);
948 return up1;
949 }
950
951 if (isl_upoly_is_nan(up2)) {
952 isl_upoly_free(up1);
953 return up2;
954 }
955
956 if (isl_upoly_is_zero(up1)) {
957 isl_upoly_free(up2);
958 return up1;
959 }
960
961 if (isl_upoly_is_zero(up2)) {
962 isl_upoly_free(up1);
963 return up2;
964 }
965
966 if (isl_upoly_is_one(up1)) {
967 isl_upoly_free(up1);
968 return up2;
969 }
970
971 if (isl_upoly_is_one(up2)) {
972 isl_upoly_free(up2);
973 return up1;
974 }
975
976 if (up1->var < up2->var)
977 return isl_upoly_mul(up2, up1);
978
979 if (up2->var < up1->var) {
980 int i;
981 struct isl_upoly_rec *rec;
982 if (isl_upoly_is_infty(up2) || isl_upoly_is_neginfty(up2)) {
983 isl_ctx *ctx = up1->ctx;
984 isl_upoly_free(up1);
985 isl_upoly_free(up2);
986 return isl_upoly_nan(ctx);
987 }
988 up1 = isl_upoly_cow(up1);
989 rec = isl_upoly_as_rec(up1);
990 if (!rec)
991 goto error;
992
993 for (i = 0; i < rec->n; ++i) {
994 rec->p[i] = isl_upoly_mul(rec->p[i],
995 isl_upoly_copy(up2));
996 if (!rec->p[i])
997 goto error;
998 }
999 isl_upoly_free(up2);
1000 return up1;
1001 }
1002
1003 if (isl_upoly_is_cst(up1))
1004 return isl_upoly_mul_cst(up1, up2);
1005
1006 return isl_upoly_mul_rec(up1, up2);
1007 error:
1008 isl_upoly_free(up1);
1009 isl_upoly_free(up2);
1010 return NULL;
1011 }
1012
1013 __isl_give struct isl_upoly *isl_upoly_pow(__isl_take struct isl_upoly *up,
1014 unsigned power)
1015 {
1016 struct isl_upoly *res;
1017
1018 if (!up)
1019 return NULL;
1020 if (power == 1)
1021 return up;
1022
1023 if (power % 2)
1024 res = isl_upoly_copy(up);
1025 else
1026 res = isl_upoly_one(up->ctx);
1027
1028 while (power >>= 1) {
1029 up = isl_upoly_mul(up, isl_upoly_copy(up));
1030 if (power % 2)
1031 res = isl_upoly_mul(res, isl_upoly_copy(up));
1032 }
1033
1034 isl_upoly_free(up);
1035 return res;
1036 }
1037
1038 __isl_give isl_qpolynomial *isl_qpolynomial_alloc(__isl_take isl_space *dim,
1039 unsigned n_div, __isl_take struct isl_upoly *up)
1040 {
1041 struct isl_qpolynomial *qp = NULL;
1042 unsigned total;
1043
1044 if (!dim || !up)
1045 goto error;
1046
1047 if (!isl_space_is_set(dim))
1048 isl_die(isl_space_get_ctx(dim), isl_error_invalid,
1049 "domain of polynomial should be a set", goto error);
1050
1051 total = isl_space_dim(dim, isl_dim_all);
1052
1053 qp = isl_calloc_type(dim->ctx, struct isl_qpolynomial);
1054 if (!qp)
1055 goto error;
1056
1057 qp->ref = 1;
1058 qp->div = isl_mat_alloc(dim->ctx, n_div, 1 + 1 + total + n_div);
1059 if (!qp->div)
1060 goto error;
1061
1062 qp->dim = dim;
1063 qp->upoly = up;
1064
1065 return qp;
1066 error:
1067 isl_space_free(dim);
1068 isl_upoly_free(up);
1069 isl_qpolynomial_free(qp);
1070 return NULL;
1071 }
1072
1073 __isl_give isl_qpolynomial *isl_qpolynomial_copy(__isl_keep isl_qpolynomial *qp)
1074 {
1075 if (!qp)
1076 return NULL;
1077
1078 qp->ref++;
1079 return qp;
1080 }
1081
1082 __isl_give isl_qpolynomial *isl_qpolynomial_dup(__isl_keep isl_qpolynomial *qp)
1083 {
1084 struct isl_qpolynomial *dup;
1085
1086 if (!qp)
1087 return NULL;
1088
1089 dup = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row,
1090 isl_upoly_copy(qp->upoly));
1091 if (!dup)
1092 return NULL;
1093 isl_mat_free(dup->div);
1094 dup->div = isl_mat_copy(qp->div);
1095 if (!dup->div)
1096 goto error;
1097
1098 return dup;
1099 error:
1100 isl_qpolynomial_free(dup);
1101 return NULL;
1102 }
1103
1104 __isl_give isl_qpolynomial *isl_qpolynomial_cow(__isl_take isl_qpolynomial *qp)
1105 {
1106 if (!qp)
1107 return NULL;
1108
1109 if (qp->ref == 1)
1110 return qp;
1111 qp->ref--;
1112 return isl_qpolynomial_dup(qp);
1113 }
1114
1115 void *isl_qpolynomial_free(__isl_take isl_qpolynomial *qp)
1116 {
1117 if (!qp)
1118 return NULL;
1119
1120 if (--qp->ref > 0)
1121 return NULL;
1122
1123 isl_space_free(qp->dim);
1124 isl_mat_free(qp->div);
1125 isl_upoly_free(qp->upoly);
1126
1127 free(qp);
1128 return NULL;
1129 }
1130
1131 __isl_give struct isl_upoly *isl_upoly_var_pow(isl_ctx *ctx, int pos, int power)
1132 {
1133 int i;
1134 struct isl_upoly_rec *rec;
1135 struct isl_upoly_cst *cst;
1136
1137 rec = isl_upoly_alloc_rec(ctx, pos, 1 + power);
1138 if (!rec)
1139 return NULL;
1140 for (i = 0; i < 1 + power; ++i) {
1141 rec->p[i] = isl_upoly_zero(ctx);
1142 if (!rec->p[i])
1143 goto error;
1144 rec->n++;
1145 }
1146 cst = isl_upoly_as_cst(rec->p[power]);
1147 isl_int_set_si(cst->n, 1);
1148
1149 return &rec->up;
1150 error:
1151 isl_upoly_free(&rec->up);
1152 return NULL;
1153 }
1154
1155 /* r array maps original positions to new positions.
1156 */
1157 static __isl_give struct isl_upoly *reorder(__isl_take struct isl_upoly *up,
1158 int *r)
1159 {
1160 int i;
1161 struct isl_upoly_rec *rec;
1162 struct isl_upoly *base;
1163 struct isl_upoly *res;
1164
1165 if (isl_upoly_is_cst(up))
1166 return up;
1167
1168 rec = isl_upoly_as_rec(up);
1169 if (!rec)
1170 goto error;
1171
1172 isl_assert(up->ctx, rec->n >= 1, goto error);
1173
1174 base = isl_upoly_var_pow(up->ctx, r[up->var], 1);
1175 res = reorder(isl_upoly_copy(rec->p[rec->n - 1]), r);
1176
1177 for (i = rec->n - 2; i >= 0; --i) {
1178 res = isl_upoly_mul(res, isl_upoly_copy(base));
1179 res = isl_upoly_sum(res, reorder(isl_upoly_copy(rec->p[i]), r));
1180 }
1181
1182 isl_upoly_free(base);
1183 isl_upoly_free(up);
1184
1185 return res;
1186 error:
1187 isl_upoly_free(up);
1188 return NULL;
1189 }
1190
1191 static int compatible_divs(__isl_keep isl_mat *div1, __isl_keep isl_mat *div2)
1192 {
1193 int n_row, n_col;
1194 int equal;
1195
1196 isl_assert(div1->ctx, div1->n_row >= div2->n_row &&
1197 div1->n_col >= div2->n_col, return -1);
1198
1199 if (div1->n_row == div2->n_row)
1200 return isl_mat_is_equal(div1, div2);
1201
1202 n_row = div1->n_row;
1203 n_col = div1->n_col;
1204 div1->n_row = div2->n_row;
1205 div1->n_col = div2->n_col;
1206
1207 equal = isl_mat_is_equal(div1, div2);
1208
1209 div1->n_row = n_row;
1210 div1->n_col = n_col;
1211
1212 return equal;
1213 }
1214
1215 static int cmp_row(__isl_keep isl_mat *div, int i, int j)
1216 {
1217 int li, lj;
1218
1219 li = isl_seq_last_non_zero(div->row[i], div->n_col);
1220 lj = isl_seq_last_non_zero(div->row[j], div->n_col);
1221
1222 if (li != lj)
1223 return li - lj;
1224
1225 return isl_seq_cmp(div->row[i], div->row[j], div->n_col);
1226 }
1227
1228 struct isl_div_sort_info {
1229 isl_mat *div;
1230 int row;
1231 };
1232
1233 static int div_sort_cmp(const void *p1, const void *p2)
1234 {
1235 const struct isl_div_sort_info *i1, *i2;
1236 i1 = (const struct isl_div_sort_info *) p1;
1237 i2 = (const struct isl_div_sort_info *) p2;
1238
1239 return cmp_row(i1->div, i1->row, i2->row);
1240 }
1241
1242 /* Sort divs and remove duplicates.
1243 */
1244 static __isl_give isl_qpolynomial *sort_divs(__isl_take isl_qpolynomial *qp)
1245 {
1246 int i;
1247 int skip;
1248 int len;
1249 struct isl_div_sort_info *array = NULL;
1250 int *pos = NULL, *at = NULL;
1251 int *reordering = NULL;
1252 unsigned div_pos;
1253
1254 if (!qp)
1255 return NULL;
1256 if (qp->div->n_row <= 1)
1257 return qp;
1258
1259 div_pos = isl_space_dim(qp->dim, isl_dim_all);
1260
1261 array = isl_alloc_array(qp->div->ctx, struct isl_div_sort_info,
1262 qp->div->n_row);
1263 pos = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1264 at = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
1265 len = qp->div->n_col - 2;
1266 reordering = isl_alloc_array(qp->div->ctx, int, len);
1267 if (!array || !pos || !at || !reordering)
1268 goto error;
1269
1270 for (i = 0; i < qp->div->n_row; ++i) {
1271 array[i].div = qp->div;
1272 array[i].row = i;
1273 pos[i] = i;
1274 at[i] = i;
1275 }
1276
1277 qsort(array, qp->div->n_row, sizeof(struct isl_div_sort_info),
1278 div_sort_cmp);
1279
1280 for (i = 0; i < div_pos; ++i)
1281 reordering[i] = i;
1282
1283 for (i = 0; i < qp->div->n_row; ++i) {
1284 if (pos[array[i].row] == i)
1285 continue;
1286 qp->div = isl_mat_swap_rows(qp->div, i, pos[array[i].row]);
1287 pos[at[i]] = pos[array[i].row];
1288 at[pos[array[i].row]] = at[i];
1289 at[i] = array[i].row;
1290 pos[array[i].row] = i;
1291 }
1292
1293 skip = 0;
1294 for (i = 0; i < len - div_pos; ++i) {
1295 if (i > 0 &&
1296 isl_seq_eq(qp->div->row[i - skip - 1],
1297 qp->div->row[i - skip], qp->div->n_col)) {
1298 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
1299 isl_mat_col_add(qp->div, 2 + div_pos + i - skip - 1,
1300 2 + div_pos + i - skip);
1301 qp->div = isl_mat_drop_cols(qp->div,
1302 2 + div_pos + i - skip, 1);
1303 skip++;
1304 }
1305 reordering[div_pos + array[i].row] = div_pos + i - skip;
1306 }
1307
1308 qp->upoly = reorder(qp->upoly, reordering);
1309
1310 if (!qp->upoly || !qp->div)
1311 goto error;
1312
1313 free(at);
1314 free(pos);
1315 free(array);
1316 free(reordering);
1317
1318 return qp;
1319 error:
1320 free(at);
1321 free(pos);
1322 free(array);
1323 free(reordering);
1324 isl_qpolynomial_free(qp);
1325 return NULL;
1326 }
1327
1328 static __isl_give struct isl_upoly *expand(__isl_take struct isl_upoly *up,
1329 int *exp, int first)
1330 {
1331 int i;
1332 struct isl_upoly_rec *rec;
1333
1334 if (isl_upoly_is_cst(up))
1335 return up;
1336
1337 if (up->var < first)
1338 return up;
1339
1340 if (exp[up->var - first] == up->var - first)
1341 return up;
1342
1343 up = isl_upoly_cow(up);
1344 if (!up)
1345 goto error;
1346
1347 up->var = exp[up->var - first] + first;
1348
1349 rec = isl_upoly_as_rec(up);
1350 if (!rec)
1351 goto error;
1352
1353 for (i = 0; i < rec->n; ++i) {
1354 rec->p[i] = expand(rec->p[i], exp, first);
1355 if (!rec->p[i])
1356 goto error;
1357 }
1358
1359 return up;
1360 error:
1361 isl_upoly_free(up);
1362 return NULL;
1363 }
1364
1365 static __isl_give isl_qpolynomial *with_merged_divs(
1366 __isl_give isl_qpolynomial *(*fn)(__isl_take isl_qpolynomial *qp1,
1367 __isl_take isl_qpolynomial *qp2),
1368 __isl_take isl_qpolynomial *qp1, __isl_take isl_qpolynomial *qp2)
1369 {
1370 int *exp1 = NULL;
1371 int *exp2 = NULL;
1372 isl_mat *div = NULL;
1373
1374 qp1 = isl_qpolynomial_cow(qp1);
1375 qp2 = isl_qpolynomial_cow(qp2);
1376
1377 if (!qp1 || !qp2)
1378 goto error;
1379
1380 isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1381 qp1->div->n_col >= qp2->div->n_col, goto error);
1382
1383 exp1 = isl_alloc_array(qp1->div->ctx, int, qp1->div->n_row);
1384 exp2 = isl_alloc_array(qp2->div->ctx, int, qp2->div->n_row);
1385 if (!exp1 || !exp2)
1386 goto error;
1387
1388 div = isl_merge_divs(qp1->div, qp2->div, exp1, exp2);
1389 if (!div)
1390 goto error;
1391
1392 isl_mat_free(qp1->div);
1393 qp1->div = isl_mat_copy(div);
1394 isl_mat_free(qp2->div);
1395 qp2->div = isl_mat_copy(div);
1396
1397 qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1398 qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1399
1400 if (!qp1->upoly || !qp2->upoly)
1401 goto error;
1402
1403 isl_mat_free(div);
1404 free(exp1);
1405 free(exp2);
1406
1407 return fn(qp1, qp2);
1408 error:
1409 isl_mat_free(div);
1410 free(exp1);
1411 free(exp2);
1412 isl_qpolynomial_free(qp1);
1413 isl_qpolynomial_free(qp2);
1414 return NULL;
1415 }
1416
1417 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1418 __isl_take isl_qpolynomial *qp2)
1419 {
1420 qp1 = isl_qpolynomial_cow(qp1);
1421
1422 if (!qp1 || !qp2)
1423 goto error;
1424
1425 if (qp1->div->n_row < qp2->div->n_row)
1426 return isl_qpolynomial_add(qp2, qp1);
1427
1428 isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1429 if (!compatible_divs(qp1->div, qp2->div))
1430 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1431
1432 qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1433 if (!qp1->upoly)
1434 goto error;
1435
1436 isl_qpolynomial_free(qp2);
1437
1438 return qp1;
1439 error:
1440 isl_qpolynomial_free(qp1);
1441 isl_qpolynomial_free(qp2);
1442 return NULL;
1443 }
1444
1445 __isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
1446 __isl_keep isl_set *dom,
1447 __isl_take isl_qpolynomial *qp1,
1448 __isl_take isl_qpolynomial *qp2)
1449 {
1450 qp1 = isl_qpolynomial_add(qp1, qp2);
1451 qp1 = isl_qpolynomial_gist(qp1, isl_set_copy(dom));
1452 return qp1;
1453 }
1454
1455 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1456 __isl_take isl_qpolynomial *qp2)
1457 {
1458 return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1459 }
1460
1461 __isl_give isl_qpolynomial *isl_qpolynomial_add_isl_int(
1462 __isl_take isl_qpolynomial *qp, isl_int v)
1463 {
1464 if (isl_int_is_zero(v))
1465 return qp;
1466
1467 qp = isl_qpolynomial_cow(qp);
1468 if (!qp)
1469 return NULL;
1470
1471 qp->upoly = isl_upoly_add_isl_int(qp->upoly, v);
1472 if (!qp->upoly)
1473 goto error;
1474
1475 return qp;
1476 error:
1477 isl_qpolynomial_free(qp);
1478 return NULL;
1479
1480 }
1481
1482 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1483 {
1484 if (!qp)
1485 return NULL;
1486
1487 return isl_qpolynomial_mul_isl_int(qp, qp->dim->ctx->negone);
1488 }
1489
1490 __isl_give isl_qpolynomial *isl_qpolynomial_mul_isl_int(
1491 __isl_take isl_qpolynomial *qp, isl_int v)
1492 {
1493 if (isl_int_is_one(v))
1494 return qp;
1495
1496 if (qp && isl_int_is_zero(v)) {
1497 isl_qpolynomial *zero;
1498 zero = isl_qpolynomial_zero_on_domain(isl_space_copy(qp->dim));
1499 isl_qpolynomial_free(qp);
1500 return zero;
1501 }
1502
1503 qp = isl_qpolynomial_cow(qp);
1504 if (!qp)
1505 return NULL;
1506
1507 qp->upoly = isl_upoly_mul_isl_int(qp->upoly, v);
1508 if (!qp->upoly)
1509 goto error;
1510
1511 return qp;
1512 error:
1513 isl_qpolynomial_free(qp);
1514 return NULL;
1515 }
1516
1517 __isl_give isl_qpolynomial *isl_qpolynomial_scale(
1518 __isl_take isl_qpolynomial *qp, isl_int v)
1519 {
1520 return isl_qpolynomial_mul_isl_int(qp, v);
1521 }
1522
1523 /* Multiply "qp" by "v".
1524 */
1525 __isl_give isl_qpolynomial *isl_qpolynomial_scale_val(
1526 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v)
1527 {
1528 if (!qp || !v)
1529 goto error;
1530
1531 if (!isl_val_is_rat(v))
1532 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
1533 "expecting rational factor", goto error);
1534
1535 if (isl_val_is_one(v)) {
1536 isl_val_free(v);
1537 return qp;
1538 }
1539
1540 if (isl_val_is_zero(v)) {
1541 isl_space *space;
1542
1543 space = isl_qpolynomial_get_domain_space(qp);
1544 isl_qpolynomial_free(qp);
1545 isl_val_free(v);
1546 return isl_qpolynomial_zero_on_domain(space);
1547 }
1548
1549 qp = isl_qpolynomial_cow(qp);
1550 if (!qp)
1551 goto error;
1552
1553 qp->upoly = isl_upoly_scale_val(qp->upoly, v);
1554 if (!qp->upoly)
1555 qp = isl_qpolynomial_free(qp);
1556
1557 isl_val_free(v);
1558 return qp;
1559 error:
1560 isl_val_free(v);
1561 isl_qpolynomial_free(qp);
1562 return NULL;
1563 }
1564
1565 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1566 __isl_take isl_qpolynomial *qp2)
1567 {
1568 qp1 = isl_qpolynomial_cow(qp1);
1569
1570 if (!qp1 || !qp2)
1571 goto error;
1572
1573 if (qp1->div->n_row < qp2->div->n_row)
1574 return isl_qpolynomial_mul(qp2, qp1);
1575
1576 isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1577 if (!compatible_divs(qp1->div, qp2->div))
1578 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1579
1580 qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1581 if (!qp1->upoly)
1582 goto error;
1583
1584 isl_qpolynomial_free(qp2);
1585
1586 return qp1;
1587 error:
1588 isl_qpolynomial_free(qp1);
1589 isl_qpolynomial_free(qp2);
1590 return NULL;
1591 }
1592
1593 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_qpolynomial *qp,
1594 unsigned power)
1595 {
1596 qp = isl_qpolynomial_cow(qp);
1597
1598 if (!qp)
1599 return NULL;
1600
1601 qp->upoly = isl_upoly_pow(qp->upoly, power);
1602 if (!qp->upoly)
1603 goto error;
1604
1605 return qp;
1606 error:
1607 isl_qpolynomial_free(qp);
1608 return NULL;
1609 }
1610
1611 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_pow(
1612 __isl_take isl_pw_qpolynomial *pwqp, unsigned power)
1613 {
1614 int i;
1615
1616 if (power == 1)
1617 return pwqp;
1618
1619 pwqp = isl_pw_qpolynomial_cow(pwqp);
1620 if (!pwqp)
1621 return NULL;
1622
1623 for (i = 0; i < pwqp->n; ++i) {
1624 pwqp->p[i].qp = isl_qpolynomial_pow(pwqp->p[i].qp, power);
1625 if (!pwqp->p[i].qp)
1626 return isl_pw_qpolynomial_free(pwqp);
1627 }
1628
1629 return pwqp;
1630 }
1631
1632 __isl_give isl_qpolynomial *isl_qpolynomial_zero_on_domain(
1633 __isl_take isl_space *dim)
1634 {
1635 if (!dim)
1636 return NULL;
1637 return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1638 }
1639
1640 __isl_give isl_qpolynomial *isl_qpolynomial_one_on_domain(
1641 __isl_take isl_space *dim)
1642 {
1643 if (!dim)
1644 return NULL;
1645 return isl_qpolynomial_alloc(dim, 0, isl_upoly_one(dim->ctx));
1646 }
1647
1648 __isl_give isl_qpolynomial *isl_qpolynomial_infty_on_domain(
1649 __isl_take isl_space *dim)
1650 {
1651 if (!dim)
1652 return NULL;
1653 return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
1654 }
1655
1656 __isl_give isl_qpolynomial *isl_qpolynomial_neginfty_on_domain(
1657 __isl_take isl_space *dim)
1658 {
1659 if (!dim)
1660 return NULL;
1661 return isl_qpolynomial_alloc(dim, 0, isl_upoly_neginfty(dim->ctx));
1662 }
1663
1664 __isl_give isl_qpolynomial *isl_qpolynomial_nan_on_domain(
1665 __isl_take isl_space *dim)
1666 {
1667 if (!dim)
1668 return NULL;
1669 return isl_qpolynomial_alloc(dim, 0, isl_upoly_nan(dim->ctx));
1670 }
1671
1672 __isl_give isl_qpolynomial *isl_qpolynomial_cst_on_domain(
1673 __isl_take isl_space *dim,
1674 isl_int v)
1675 {
1676 struct isl_qpolynomial *qp;
1677 struct isl_upoly_cst *cst;
1678
1679 if (!dim)
1680 return NULL;
1681
1682 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1683 if (!qp)
1684 return NULL;
1685
1686 cst = isl_upoly_as_cst(qp->upoly);
1687 isl_int_set(cst->n, v);
1688
1689 return qp;
1690 }
1691
1692 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1693 isl_int *n, isl_int *d)
1694 {
1695 struct isl_upoly_cst *cst;
1696
1697 if (!qp)
1698 return -1;
1699
1700 if (!isl_upoly_is_cst(qp->upoly))
1701 return 0;
1702
1703 cst = isl_upoly_as_cst(qp->upoly);
1704 if (!cst)
1705 return -1;
1706
1707 if (n)
1708 isl_int_set(*n, cst->n);
1709 if (d)
1710 isl_int_set(*d, cst->d);
1711
1712 return 1;
1713 }
1714
1715 /* Return the constant term of "up".
1716 */
1717 static __isl_give isl_val *isl_upoly_get_constant_val(
1718 __isl_keep struct isl_upoly *up)
1719 {
1720 struct isl_upoly_cst *cst;
1721
1722 if (!up)
1723 return NULL;
1724
1725 while (!isl_upoly_is_cst(up)) {
1726 struct isl_upoly_rec *rec;
1727
1728 rec = isl_upoly_as_rec(up);
1729 if (!rec)
1730 return NULL;
1731 up = rec->p[0];
1732 }
1733
1734 cst = isl_upoly_as_cst(up);
1735 if (!cst)
1736 return NULL;
1737 return isl_val_rat_from_isl_int(cst->up.ctx, cst->n, cst->d);
1738 }
1739
1740 /* Return the constant term of "qp".
1741 */
1742 __isl_give isl_val *isl_qpolynomial_get_constant_val(
1743 __isl_keep isl_qpolynomial *qp)
1744 {
1745 if (!qp)
1746 return NULL;
1747
1748 return isl_upoly_get_constant_val(qp->upoly);
1749 }
1750
1751 int isl_upoly_is_affine(__isl_keep struct isl_upoly *up)
1752 {
1753 int is_cst;
1754 struct isl_upoly_rec *rec;
1755
1756 if (!up)
1757 return -1;
1758
1759 if (up->var < 0)
1760 return 1;
1761
1762 rec = isl_upoly_as_rec(up);
1763 if (!rec)
1764 return -1;
1765
1766 if (rec->n > 2)
1767 return 0;
1768
1769 isl_assert(up->ctx, rec->n > 1, return -1);
1770
1771 is_cst = isl_upoly_is_cst(rec->p[1]);
1772 if (is_cst < 0)
1773 return -1;
1774 if (!is_cst)
1775 return 0;
1776
1777 return isl_upoly_is_affine(rec->p[0]);
1778 }
1779
1780 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp)
1781 {
1782 if (!qp)
1783 return -1;
1784
1785 if (qp->div->n_row > 0)
1786 return 0;
1787
1788 return isl_upoly_is_affine(qp->upoly);
1789 }
1790
1791 static void update_coeff(__isl_keep isl_vec *aff,
1792 __isl_keep struct isl_upoly_cst *cst, int pos)
1793 {
1794 isl_int gcd;
1795 isl_int f;
1796
1797 if (isl_int_is_zero(cst->n))
1798 return;
1799
1800 isl_int_init(gcd);
1801 isl_int_init(f);
1802 isl_int_gcd(gcd, cst->d, aff->el[0]);
1803 isl_int_divexact(f, cst->d, gcd);
1804 isl_int_divexact(gcd, aff->el[0], gcd);
1805 isl_seq_scale(aff->el, aff->el, f, aff->size);
1806 isl_int_mul(aff->el[1 + pos], gcd, cst->n);
1807 isl_int_clear(gcd);
1808 isl_int_clear(f);
1809 }
1810
1811 int isl_upoly_update_affine(__isl_keep struct isl_upoly *up,
1812 __isl_keep isl_vec *aff)
1813 {
1814 struct isl_upoly_cst *cst;
1815 struct isl_upoly_rec *rec;
1816
1817 if (!up || !aff)
1818 return -1;
1819
1820 if (up->var < 0) {
1821 struct isl_upoly_cst *cst;
1822
1823 cst = isl_upoly_as_cst(up);
1824 if (!cst)
1825 return -1;
1826 update_coeff(aff, cst, 0);
1827 return 0;
1828 }
1829
1830 rec = isl_upoly_as_rec(up);
1831 if (!rec)
1832 return -1;
1833 isl_assert(up->ctx, rec->n == 2, return -1);
1834
1835 cst = isl_upoly_as_cst(rec->p[1]);
1836 if (!cst)
1837 return -1;
1838 update_coeff(aff, cst, 1 + up->var);
1839
1840 return isl_upoly_update_affine(rec->p[0], aff);
1841 }
1842
1843 __isl_give isl_vec *isl_qpolynomial_extract_affine(
1844 __isl_keep isl_qpolynomial *qp)
1845 {
1846 isl_vec *aff;
1847 unsigned d;
1848
1849 if (!qp)
1850 return NULL;
1851
1852 d = isl_space_dim(qp->dim, isl_dim_all);
1853 aff = isl_vec_alloc(qp->div->ctx, 2 + d + qp->div->n_row);
1854 if (!aff)
1855 return NULL;
1856
1857 isl_seq_clr(aff->el + 1, 1 + d + qp->div->n_row);
1858 isl_int_set_si(aff->el[0], 1);
1859
1860 if (isl_upoly_update_affine(qp->upoly, aff) < 0)
1861 goto error;
1862
1863 return aff;
1864 error:
1865 isl_vec_free(aff);
1866 return NULL;
1867 }
1868
1869 int isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial *qp1,
1870 __isl_keep isl_qpolynomial *qp2)
1871 {
1872 int equal;
1873
1874 if (!qp1 || !qp2)
1875 return -1;
1876
1877 equal = isl_space_is_equal(qp1->dim, qp2->dim);
1878 if (equal < 0 || !equal)
1879 return equal;
1880
1881 equal = isl_mat_is_equal(qp1->div, qp2->div);
1882 if (equal < 0 || !equal)
1883 return equal;
1884
1885 return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
1886 }
1887
1888 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
1889 {
1890 int i;
1891 struct isl_upoly_rec *rec;
1892
1893 if (isl_upoly_is_cst(up)) {
1894 struct isl_upoly_cst *cst;
1895 cst = isl_upoly_as_cst(up);
1896 if (!cst)
1897 return;
1898 isl_int_lcm(*d, *d, cst->d);
1899 return;
1900 }
1901
1902 rec = isl_upoly_as_rec(up);
1903 if (!rec)
1904 return;
1905
1906 for (i = 0; i < rec->n; ++i)
1907 upoly_update_den(rec->p[i], d);
1908 }
1909
1910 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
1911 {
1912 isl_int_set_si(*d, 1);
1913 if (!qp)
1914 return;
1915 upoly_update_den(qp->upoly, d);
1916 }
1917
1918 __isl_give isl_qpolynomial *isl_qpolynomial_var_pow_on_domain(
1919 __isl_take isl_space *dim, int pos, int power)
1920 {
1921 struct isl_ctx *ctx;
1922
1923 if (!dim)
1924 return NULL;
1925
1926 ctx = dim->ctx;
1927
1928 return isl_qpolynomial_alloc(dim, 0, isl_upoly_var_pow(ctx, pos, power));
1929 }
1930
1931 __isl_give isl_qpolynomial *isl_qpolynomial_var_on_domain(__isl_take isl_space *dim,
1932 enum isl_dim_type type, unsigned pos)
1933 {
1934 if (!dim)
1935 return NULL;
1936
1937 isl_assert(dim->ctx, isl_space_dim(dim, isl_dim_in) == 0, goto error);
1938 isl_assert(dim->ctx, pos < isl_space_dim(dim, type), goto error);
1939
1940 if (type == isl_dim_set)
1941 pos += isl_space_dim(dim, isl_dim_param);
1942
1943 return isl_qpolynomial_var_pow_on_domain(dim, pos, 1);
1944 error:
1945 isl_space_free(dim);
1946 return NULL;
1947 }
1948
1949 __isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
1950 unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
1951 {
1952 int i;
1953 struct isl_upoly_rec *rec;
1954 struct isl_upoly *base, *res;
1955
1956 if (!up)
1957 return NULL;
1958
1959 if (isl_upoly_is_cst(up))
1960 return up;
1961
1962 if (up->var < first)
1963 return up;
1964
1965 rec = isl_upoly_as_rec(up);
1966 if (!rec)
1967 goto error;
1968
1969 isl_assert(up->ctx, rec->n >= 1, goto error);
1970
1971 if (up->var >= first + n)
1972 base = isl_upoly_var_pow(up->ctx, up->var, 1);
1973 else
1974 base = isl_upoly_copy(subs[up->var - first]);
1975
1976 res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
1977 for (i = rec->n - 2; i >= 0; --i) {
1978 struct isl_upoly *t;
1979 t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
1980 res = isl_upoly_mul(res, isl_upoly_copy(base));
1981 res = isl_upoly_sum(res, t);
1982 }
1983
1984 isl_upoly_free(base);
1985 isl_upoly_free(up);
1986
1987 return res;
1988 error:
1989 isl_upoly_free(up);
1990 return NULL;
1991 }
1992
1993 __isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
1994 isl_int denom, unsigned len)
1995 {
1996 int i;
1997 struct isl_upoly *up;
1998
1999 isl_assert(ctx, len >= 1, return NULL);
2000
2001 up = isl_upoly_rat_cst(ctx, f[0], denom);
2002 for (i = 0; i < len - 1; ++i) {
2003 struct isl_upoly *t;
2004 struct isl_upoly *c;
2005
2006 if (isl_int_is_zero(f[1 + i]))
2007 continue;
2008
2009 c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
2010 t = isl_upoly_var_pow(ctx, i, 1);
2011 t = isl_upoly_mul(c, t);
2012 up = isl_upoly_sum(up, t);
2013 }
2014
2015 return up;
2016 }
2017
2018 /* Remove common factor of non-constant terms and denominator.
2019 */
2020 static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
2021 {
2022 isl_ctx *ctx = qp->div->ctx;
2023 unsigned total = qp->div->n_col - 2;
2024
2025 isl_seq_gcd(qp->div->row[div] + 2, total, &ctx->normalize_gcd);
2026 isl_int_gcd(ctx->normalize_gcd,
2027 ctx->normalize_gcd, qp->div->row[div][0]);
2028 if (isl_int_is_one(ctx->normalize_gcd))
2029 return;
2030
2031 isl_seq_scale_down(qp->div->row[div] + 2, qp->div->row[div] + 2,
2032 ctx->normalize_gcd, total);
2033 isl_int_divexact(qp->div->row[div][0], qp->div->row[div][0],
2034 ctx->normalize_gcd);
2035 isl_int_fdiv_q(qp->div->row[div][1], qp->div->row[div][1],
2036 ctx->normalize_gcd);
2037 }
2038
2039 /* Replace the integer division identified by "div" by the polynomial "s".
2040 * The integer division is assumed not to appear in the definition
2041 * of any other integer divisions.
2042 */
2043 static __isl_give isl_qpolynomial *substitute_div(
2044 __isl_take isl_qpolynomial *qp,
2045 int div, __isl_take struct isl_upoly *s)
2046 {
2047 int i;
2048 int total;
2049 int *reordering;
2050
2051 if (!qp || !s)
2052 goto error;
2053
2054 qp = isl_qpolynomial_cow(qp);
2055 if (!qp)
2056 goto error;
2057
2058 total = isl_space_dim(qp->dim, isl_dim_all);
2059 qp->upoly = isl_upoly_subs(qp->upoly, total + div, 1, &s);
2060 if (!qp->upoly)
2061 goto error;
2062
2063 reordering = isl_alloc_array(qp->dim->ctx, int, total + qp->div->n_row);
2064 if (!reordering)
2065 goto error;
2066 for (i = 0; i < total + div; ++i)
2067 reordering[i] = i;
2068 for (i = total + div + 1; i < total + qp->div->n_row; ++i)
2069 reordering[i] = i - 1;
2070 qp->div = isl_mat_drop_rows(qp->div, div, 1);
2071 qp->div = isl_mat_drop_cols(qp->div, 2 + total + div, 1);
2072 qp->upoly = reorder(qp->upoly, reordering);
2073 free(reordering);
2074
2075 if (!qp->upoly || !qp->div)
2076 goto error;
2077
2078 isl_upoly_free(s);
2079 return qp;
2080 error:
2081 isl_qpolynomial_free(qp);
2082 isl_upoly_free(s);
2083 return NULL;
2084 }
2085
2086 /* Replace all integer divisions [e/d] that turn out to not actually be integer
2087 * divisions because d is equal to 1 by their definition, i.e., e.
2088 */
2089 static __isl_give isl_qpolynomial *substitute_non_divs(
2090 __isl_take isl_qpolynomial *qp)
2091 {
2092 int i, j;
2093 int total;
2094 struct isl_upoly *s;
2095
2096 if (!qp)
2097 return NULL;
2098
2099 total = isl_space_dim(qp->dim, isl_dim_all);
2100 for (i = 0; qp && i < qp->div->n_row; ++i) {
2101 if (!isl_int_is_one(qp->div->row[i][0]))
2102 continue;
2103 for (j = i + 1; j < qp->div->n_row; ++j) {
2104 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
2105 continue;
2106 isl_seq_combine(qp->div->row[j] + 1,
2107 qp->div->ctx->one, qp->div->row[j] + 1,
2108 qp->div->row[j][2 + total + i],
2109 qp->div->row[i] + 1, 1 + total + i);
2110 isl_int_set_si(qp->div->row[j][2 + total + i], 0);
2111 normalize_div(qp, j);
2112 }
2113 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
2114 qp->div->row[i][0], qp->div->n_col - 1);
2115 qp = substitute_div(qp, i, s);
2116 --i;
2117 }
2118
2119 return qp;
2120 }
2121
2122 /* Reduce the coefficients of div "div" to lie in the interval [0, d-1],
2123 * with d the denominator. When replacing the coefficient e of x by
2124 * d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
2125 * inside the division, so we need to add floor(e/d) * x outside.
2126 * That is, we replace q by q' + floor(e/d) * x and we therefore need
2127 * to adjust the coefficient of x in each later div that depends on the
2128 * current div "div" and also in the affine expression "aff"
2129 * (if it too depends on "div").
2130 */
2131 static void reduce_div(__isl_keep isl_qpolynomial *qp, int div,
2132 __isl_keep isl_vec *aff)
2133 {
2134 int i, j;
2135 isl_int v;
2136 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2137
2138 isl_int_init(v);
2139 for (i = 0; i < 1 + total + div; ++i) {
2140 if (isl_int_is_nonneg(qp->div->row[div][1 + i]) &&
2141 isl_int_lt(qp->div->row[div][1 + i], qp->div->row[div][0]))
2142 continue;
2143 isl_int_fdiv_q(v, qp->div->row[div][1 + i], qp->div->row[div][0]);
2144 isl_int_fdiv_r(qp->div->row[div][1 + i],
2145 qp->div->row[div][1 + i], qp->div->row[div][0]);
2146 if (!isl_int_is_zero(aff->el[1 + total + div]))
2147 isl_int_addmul(aff->el[i], v, aff->el[1 + total + div]);
2148 for (j = div + 1; j < qp->div->n_row; ++j) {
2149 if (isl_int_is_zero(qp->div->row[j][2 + total + div]))
2150 continue;
2151 isl_int_addmul(qp->div->row[j][1 + i],
2152 v, qp->div->row[j][2 + total + div]);
2153 }
2154 }
2155 isl_int_clear(v);
2156 }
2157
2158 /* Check if the last non-zero coefficient is bigger that half of the
2159 * denominator. If so, we will invert the div to further reduce the number
2160 * of distinct divs that may appear.
2161 * If the last non-zero coefficient is exactly half the denominator,
2162 * then we continue looking for earlier coefficients that are bigger
2163 * than half the denominator.
2164 */
2165 static int needs_invert(__isl_keep isl_mat *div, int row)
2166 {
2167 int i;
2168 int cmp;
2169
2170 for (i = div->n_col - 1; i >= 1; --i) {
2171 if (isl_int_is_zero(div->row[row][i]))
2172 continue;
2173 isl_int_mul_ui(div->row[row][i], div->row[row][i], 2);
2174 cmp = isl_int_cmp(div->row[row][i], div->row[row][0]);
2175 isl_int_divexact_ui(div->row[row][i], div->row[row][i], 2);
2176 if (cmp)
2177 return cmp > 0;
2178 if (i == 1)
2179 return 1;
2180 }
2181
2182 return 0;
2183 }
2184
2185 /* Replace div "div" q = [e/d] by -[(-e+(d-1))/d].
2186 * We only invert the coefficients of e (and the coefficient of q in
2187 * later divs and in "aff"). After calling this function, the
2188 * coefficients of e should be reduced again.
2189 */
2190 static void invert_div(__isl_keep isl_qpolynomial *qp, int div,
2191 __isl_keep isl_vec *aff)
2192 {
2193 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2194
2195 isl_seq_neg(qp->div->row[div] + 1,
2196 qp->div->row[div] + 1, qp->div->n_col - 1);
2197 isl_int_sub_ui(qp->div->row[div][1], qp->div->row[div][1], 1);
2198 isl_int_add(qp->div->row[div][1],
2199 qp->div->row[div][1], qp->div->row[div][0]);
2200 if (!isl_int_is_zero(aff->el[1 + total + div]))
2201 isl_int_neg(aff->el[1 + total + div], aff->el[1 + total + div]);
2202 isl_mat_col_mul(qp->div, 2 + total + div,
2203 qp->div->ctx->negone, 2 + total + div);
2204 }
2205
2206 /* Assuming "qp" is a monomial, reduce all its divs to have coefficients
2207 * in the interval [0, d-1], with d the denominator and such that the
2208 * last non-zero coefficient that is not equal to d/2 is smaller than d/2.
2209 *
2210 * After the reduction, some divs may have become redundant or identical,
2211 * so we call substitute_non_divs and sort_divs. If these functions
2212 * eliminate divs or merge two or more divs into one, the coefficients
2213 * of the enclosing divs may have to be reduced again, so we call
2214 * ourselves recursively if the number of divs decreases.
2215 */
2216 static __isl_give isl_qpolynomial *reduce_divs(__isl_take isl_qpolynomial *qp)
2217 {
2218 int i;
2219 isl_vec *aff = NULL;
2220 struct isl_upoly *s;
2221 unsigned n_div;
2222
2223 if (!qp)
2224 return NULL;
2225
2226 aff = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
2227 aff = isl_vec_clr(aff);
2228 if (!aff)
2229 goto error;
2230
2231 isl_int_set_si(aff->el[1 + qp->upoly->var], 1);
2232
2233 for (i = 0; i < qp->div->n_row; ++i) {
2234 normalize_div(qp, i);
2235 reduce_div(qp, i, aff);
2236 if (needs_invert(qp->div, i)) {
2237 invert_div(qp, i, aff);
2238 reduce_div(qp, i, aff);
2239 }
2240 }
2241
2242 s = isl_upoly_from_affine(qp->div->ctx, aff->el,
2243 qp->div->ctx->one, aff->size);
2244 qp->upoly = isl_upoly_subs(qp->upoly, qp->upoly->var, 1, &s);
2245 isl_upoly_free(s);
2246 if (!qp->upoly)
2247 goto error;
2248
2249 isl_vec_free(aff);
2250
2251 n_div = qp->div->n_row;
2252 qp = substitute_non_divs(qp);
2253 qp = sort_divs(qp);
2254 if (qp && qp->div->n_row < n_div)
2255 return reduce_divs(qp);
2256
2257 return qp;
2258 error:
2259 isl_qpolynomial_free(qp);
2260 isl_vec_free(aff);
2261 return NULL;
2262 }
2263
2264 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst_on_domain(
2265 __isl_take isl_space *dim, const isl_int n, const isl_int d)
2266 {
2267 struct isl_qpolynomial *qp;
2268 struct isl_upoly_cst *cst;
2269
2270 if (!dim)
2271 return NULL;
2272
2273 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
2274 if (!qp)
2275 return NULL;
2276
2277 cst = isl_upoly_as_cst(qp->upoly);
2278 isl_int_set(cst->n, n);
2279 isl_int_set(cst->d, d);
2280
2281 return qp;
2282 }
2283
2284 /* Return an isl_qpolynomial that is equal to "val" on domain space "domain".
2285 */
2286 __isl_give isl_qpolynomial *isl_qpolynomial_val_on_domain(
2287 __isl_take isl_space *domain, __isl_take isl_val *val)
2288 {
2289 isl_qpolynomial *qp;
2290 struct isl_upoly_cst *cst;
2291
2292 if (!domain || !val)
2293 goto error;
2294
2295 qp = isl_qpolynomial_alloc(domain, 0, isl_upoly_zero(domain->ctx));
2296 if (!qp)
2297 goto error;
2298
2299 cst = isl_upoly_as_cst(qp->upoly);
2300 isl_int_set(cst->n, val->n);
2301 isl_int_set(cst->d, val->d);
2302
2303 isl_val_free(val);
2304 return qp;
2305 error:
2306 isl_space_free(domain);
2307 isl_val_free(val);
2308 return NULL;
2309 }
2310
2311 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
2312 {
2313 struct isl_upoly_rec *rec;
2314 int i;
2315
2316 if (!up)
2317 return -1;
2318
2319 if (isl_upoly_is_cst(up))
2320 return 0;
2321
2322 if (up->var < d)
2323 active[up->var] = 1;
2324
2325 rec = isl_upoly_as_rec(up);
2326 for (i = 0; i < rec->n; ++i)
2327 if (up_set_active(rec->p[i], active, d) < 0)
2328 return -1;
2329
2330 return 0;
2331 }
2332
2333 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
2334 {
2335 int i, j;
2336 int d = isl_space_dim(qp->dim, isl_dim_all);
2337
2338 if (!qp || !active)
2339 return -1;
2340
2341 for (i = 0; i < d; ++i)
2342 for (j = 0; j < qp->div->n_row; ++j) {
2343 if (isl_int_is_zero(qp->div->row[j][2 + i]))
2344 continue;
2345 active[i] = 1;
2346 break;
2347 }
2348
2349 return up_set_active(qp->upoly, active, d);
2350 }
2351
2352 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
2353 enum isl_dim_type type, unsigned first, unsigned n)
2354 {
2355 int i;
2356 int *active = NULL;
2357 int involves = 0;
2358
2359 if (!qp)
2360 return -1;
2361 if (n == 0)
2362 return 0;
2363
2364 isl_assert(qp->dim->ctx,
2365 first + n <= isl_qpolynomial_dim(qp, type), return -1);
2366 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2367 type == isl_dim_in, return -1);
2368
2369 active = isl_calloc_array(qp->dim->ctx, int,
2370 isl_space_dim(qp->dim, isl_dim_all));
2371 if (set_active(qp, active) < 0)
2372 goto error;
2373
2374 if (type == isl_dim_in)
2375 first += isl_space_dim(qp->dim, isl_dim_param);
2376 for (i = 0; i < n; ++i)
2377 if (active[first + i]) {
2378 involves = 1;
2379 break;
2380 }
2381
2382 free(active);
2383
2384 return involves;
2385 error:
2386 free(active);
2387 return -1;
2388 }
2389
2390 /* Remove divs that do not appear in the quasi-polynomial, nor in any
2391 * of the divs that do appear in the quasi-polynomial.
2392 */
2393 static __isl_give isl_qpolynomial *remove_redundant_divs(
2394 __isl_take isl_qpolynomial *qp)
2395 {
2396 int i, j;
2397 int d;
2398 int len;
2399 int skip;
2400 int *active = NULL;
2401 int *reordering = NULL;
2402 int redundant = 0;
2403 int n_div;
2404 isl_ctx *ctx;
2405
2406 if (!qp)
2407 return NULL;
2408 if (qp->div->n_row == 0)
2409 return qp;
2410
2411 d = isl_space_dim(qp->dim, isl_dim_all);
2412 len = qp->div->n_col - 2;
2413 ctx = isl_qpolynomial_get_ctx(qp);
2414 active = isl_calloc_array(ctx, int, len);
2415 if (!active)
2416 goto error;
2417
2418 if (up_set_active(qp->upoly, active, len) < 0)
2419 goto error;
2420
2421 for (i = qp->div->n_row - 1; i >= 0; --i) {
2422 if (!active[d + i]) {
2423 redundant = 1;
2424 continue;
2425 }
2426 for (j = 0; j < i; ++j) {
2427 if (isl_int_is_zero(qp->div->row[i][2 + d + j]))
2428 continue;
2429 active[d + j] = 1;
2430 break;
2431 }
2432 }
2433
2434 if (!redundant) {
2435 free(active);
2436 return qp;
2437 }
2438
2439 reordering = isl_alloc_array(qp->div->ctx, int, len);
2440 if (!reordering)
2441 goto error;
2442
2443 for (i = 0; i < d; ++i)
2444 reordering[i] = i;
2445
2446 skip = 0;
2447 n_div = qp->div->n_row;
2448 for (i = 0; i < n_div; ++i) {
2449 if (!active[d + i]) {
2450 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
2451 qp->div = isl_mat_drop_cols(qp->div,
2452 2 + d + i - skip, 1);
2453 skip++;
2454 }
2455 reordering[d + i] = d + i - skip;
2456 }
2457
2458 qp->upoly = reorder(qp->upoly, reordering);
2459
2460 if (!qp->upoly || !qp->div)
2461 goto error;
2462
2463 free(active);
2464 free(reordering);
2465
2466 return qp;
2467 error:
2468 free(active);
2469 free(reordering);
2470 isl_qpolynomial_free(qp);
2471 return NULL;
2472 }
2473
2474 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
2475 unsigned first, unsigned n)
2476 {
2477 int i;
2478 struct isl_upoly_rec *rec;
2479
2480 if (!up)
2481 return NULL;
2482 if (n == 0 || up->var < 0 || up->var < first)
2483 return up;
2484 if (up->var < first + n) {
2485 up = replace_by_constant_term(up);
2486 return isl_upoly_drop(up, first, n);
2487 }
2488 up = isl_upoly_cow(up);
2489 if (!up)
2490 return NULL;
2491 up->var -= n;
2492 rec = isl_upoly_as_rec(up);
2493 if (!rec)
2494 goto error;
2495
2496 for (i = 0; i < rec->n; ++i) {
2497 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
2498 if (!rec->p[i])
2499 goto error;
2500 }
2501
2502 return up;
2503 error:
2504 isl_upoly_free(up);
2505 return NULL;
2506 }
2507
2508 __isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
2509 __isl_take isl_qpolynomial *qp,
2510 enum isl_dim_type type, unsigned pos, const char *s)
2511 {
2512 qp = isl_qpolynomial_cow(qp);
2513 if (!qp)
2514 return NULL;
2515 qp->dim = isl_space_set_dim_name(qp->dim, type, pos, s);
2516 if (!qp->dim)
2517 goto error;
2518 return qp;
2519 error:
2520 isl_qpolynomial_free(qp);
2521 return NULL;
2522 }
2523
2524 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
2525 __isl_take isl_qpolynomial *qp,
2526 enum isl_dim_type type, unsigned first, unsigned n)
2527 {
2528 if (!qp)
2529 return NULL;
2530 if (type == isl_dim_out)
2531 isl_die(qp->dim->ctx, isl_error_invalid,
2532 "cannot drop output/set dimension",
2533 goto error);
2534 if (type == isl_dim_in)
2535 type = isl_dim_set;
2536 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2537 return qp;
2538
2539 qp = isl_qpolynomial_cow(qp);
2540 if (!qp)
2541 return NULL;
2542
2543 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
2544 goto error);
2545 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2546 type == isl_dim_set, goto error);
2547
2548 qp->dim = isl_space_drop_dims(qp->dim, type, first, n);
2549 if (!qp->dim)
2550 goto error;
2551
2552 if (type == isl_dim_set)
2553 first += isl_space_dim(qp->dim, isl_dim_param);
2554
2555 qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
2556 if (!qp->div)
2557 goto error;
2558
2559 qp->upoly = isl_upoly_drop(qp->upoly, first, n);
2560 if (!qp->upoly)
2561 goto error;
2562
2563 return qp;
2564 error:
2565 isl_qpolynomial_free(qp);
2566 return NULL;
2567 }
2568
2569 /* Project the domain of the quasi-polynomial onto its parameter space.
2570 * The quasi-polynomial may not involve any of the domain dimensions.
2571 */
2572 __isl_give isl_qpolynomial *isl_qpolynomial_project_domain_on_params(
2573 __isl_take isl_qpolynomial *qp)
2574 {
2575 isl_space *space;
2576 unsigned n;
2577 int involves;
2578
2579 n = isl_qpolynomial_dim(qp, isl_dim_in);
2580 involves = isl_qpolynomial_involves_dims(qp, isl_dim_in, 0, n);
2581 if (involves < 0)
2582 return isl_qpolynomial_free(qp);
2583 if (involves)
2584 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
2585 "polynomial involves some of the domain dimensions",
2586 return isl_qpolynomial_free(qp));
2587 qp = isl_qpolynomial_drop_dims(qp, isl_dim_in, 0, n);
2588 space = isl_qpolynomial_get_domain_space(qp);
2589 space = isl_space_params(space);
2590 qp = isl_qpolynomial_reset_domain_space(qp, space);
2591 return qp;
2592 }
2593
2594 static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted(
2595 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2596 {
2597 int i, j, k;
2598 isl_int denom;
2599 unsigned total;
2600 unsigned n_div;
2601 struct isl_upoly *up;
2602
2603 if (!eq)
2604 goto error;
2605 if (eq->n_eq == 0) {
2606 isl_basic_set_free(eq);
2607 return qp;
2608 }
2609
2610 qp = isl_qpolynomial_cow(qp);
2611 if (!qp)
2612 goto error;
2613 qp->div = isl_mat_cow(qp->div);
2614 if (!qp->div)
2615 goto error;
2616
2617 total = 1 + isl_space_dim(eq->dim, isl_dim_all);
2618 n_div = eq->n_div;
2619 isl_int_init(denom);
2620 for (i = 0; i < eq->n_eq; ++i) {
2621 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
2622 if (j < 0 || j == 0 || j >= total)
2623 continue;
2624
2625 for (k = 0; k < qp->div->n_row; ++k) {
2626 if (isl_int_is_zero(qp->div->row[k][1 + j]))
2627 continue;
2628 isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
2629 &qp->div->row[k][0]);
2630 normalize_div(qp, k);
2631 }
2632
2633 if (isl_int_is_pos(eq->eq[i][j]))
2634 isl_seq_neg(eq->eq[i], eq->eq[i], total);
2635 isl_int_abs(denom, eq->eq[i][j]);
2636 isl_int_set_si(eq->eq[i][j], 0);
2637
2638 up = isl_upoly_from_affine(qp->dim->ctx,
2639 eq->eq[i], denom, total);
2640 qp->upoly = isl_upoly_subs(qp->upoly, j - 1, 1, &up);
2641 isl_upoly_free(up);
2642 }
2643 isl_int_clear(denom);
2644
2645 if (!qp->upoly)
2646 goto error;
2647
2648 isl_basic_set_free(eq);
2649
2650 qp = substitute_non_divs(qp);
2651 qp = sort_divs(qp);
2652
2653 return qp;
2654 error:
2655 isl_basic_set_free(eq);
2656 isl_qpolynomial_free(qp);
2657 return NULL;
2658 }
2659
2660 /* Exploit the equalities in "eq" to simplify the quasi-polynomial.
2661 */
2662 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
2663 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2664 {
2665 if (!qp || !eq)
2666 goto error;
2667 if (qp->div->n_row > 0)
2668 eq = isl_basic_set_add_dims(eq, isl_dim_set, qp->div->n_row);
2669 return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
2670 error:
2671 isl_basic_set_free(eq);
2672 isl_qpolynomial_free(qp);
2673 return NULL;
2674 }
2675
2676 static __isl_give isl_basic_set *add_div_constraints(
2677 __isl_take isl_basic_set *bset, __isl_take isl_mat *div)
2678 {
2679 int i;
2680 unsigned total;
2681
2682 if (!bset || !div)
2683 goto error;
2684
2685 bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
2686 if (!bset)
2687 goto error;
2688 total = isl_basic_set_total_dim(bset);
2689 for (i = 0; i < div->n_row; ++i)
2690 if (isl_basic_set_add_div_constraints_var(bset,
2691 total - div->n_row + i, div->row[i]) < 0)
2692 goto error;
2693
2694 isl_mat_free(div);
2695 return bset;
2696 error:
2697 isl_mat_free(div);
2698 isl_basic_set_free(bset);
2699 return NULL;
2700 }
2701
2702 /* Look for equalities among the variables shared by context and qp
2703 * and the integer divisions of qp, if any.
2704 * The equalities are then used to eliminate variables and/or integer
2705 * divisions from qp.
2706 */
2707 __isl_give isl_qpolynomial *isl_qpolynomial_gist(
2708 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2709 {
2710 isl_basic_set *aff;
2711
2712 if (!qp)
2713 goto error;
2714 if (qp->div->n_row > 0) {
2715 isl_basic_set *bset;
2716 context = isl_set_add_dims(context, isl_dim_set,
2717 qp->div->n_row);
2718 bset = isl_basic_set_universe(isl_set_get_space(context));
2719 bset = add_div_constraints(bset, isl_mat_copy(qp->div));
2720 context = isl_set_intersect(context,
2721 isl_set_from_basic_set(bset));
2722 }
2723
2724 aff = isl_set_affine_hull(context);
2725 return isl_qpolynomial_substitute_equalities_lifted(qp, aff);
2726 error:
2727 isl_qpolynomial_free(qp);
2728 isl_set_free(context);
2729 return NULL;
2730 }
2731
2732 __isl_give isl_qpolynomial *isl_qpolynomial_gist_params(
2733 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2734 {
2735 isl_space *space = isl_qpolynomial_get_domain_space(qp);
2736 isl_set *dom_context = isl_set_universe(space);
2737 dom_context = isl_set_intersect_params(dom_context, context);
2738 return isl_qpolynomial_gist(qp, dom_context);
2739 }
2740
2741 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_qpolynomial(
2742 __isl_take isl_qpolynomial *qp)
2743 {
2744 isl_set *dom;
2745
2746 if (!qp)
2747 return NULL;
2748 if (isl_qpolynomial_is_zero(qp)) {
2749 isl_space *dim = isl_qpolynomial_get_space(qp);
2750 isl_qpolynomial_free(qp);
2751 return isl_pw_qpolynomial_zero(dim);
2752 }
2753
2754 dom = isl_set_universe(isl_qpolynomial_get_domain_space(qp));
2755 return isl_pw_qpolynomial_alloc(dom, qp);
2756 }
2757
2758 #undef PW
2759 #define PW isl_pw_qpolynomial
2760 #undef EL
2761 #define EL isl_qpolynomial
2762 #undef EL_IS_ZERO
2763 #define EL_IS_ZERO is_zero
2764 #undef ZERO
2765 #define ZERO zero
2766 #undef IS_ZERO
2767 #define IS_ZERO is_zero
2768 #undef FIELD
2769 #define FIELD qp
2770 #undef DEFAULT_IS_ZERO
2771 #define DEFAULT_IS_ZERO 1
2772
2773 #define NO_PULLBACK
2774
2775 #include <isl_pw_templ.c>
2776
2777 #undef UNION
2778 #define UNION isl_union_pw_qpolynomial
2779 #undef PART
2780 #define PART isl_pw_qpolynomial
2781 #undef PARTS
2782 #define PARTS pw_qpolynomial
2783 #define ALIGN_DOMAIN
2784
2785 #include <isl_union_templ.c>
2786
2787 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
2788 {
2789 if (!pwqp)
2790 return -1;
2791
2792 if (pwqp->n != -1)
2793 return 0;
2794
2795 if (!isl_set_plain_is_universe(pwqp->p[0].set))
2796 return 0;
2797
2798 return isl_qpolynomial_is_one(pwqp->p[0].qp);
2799 }
2800
2801 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
2802 __isl_take isl_pw_qpolynomial *pwqp1,
2803 __isl_take isl_pw_qpolynomial *pwqp2)
2804 {
2805 return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2);
2806 }
2807
2808 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
2809 __isl_take isl_pw_qpolynomial *pwqp1,
2810 __isl_take isl_pw_qpolynomial *pwqp2)
2811 {
2812 int i, j, n;
2813 struct isl_pw_qpolynomial *res;
2814
2815 if (!pwqp1 || !pwqp2)
2816 goto error;
2817
2818 isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim),
2819 goto error);
2820
2821 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
2822 isl_pw_qpolynomial_free(pwqp2);
2823 return pwqp1;
2824 }
2825
2826 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
2827 isl_pw_qpolynomial_free(pwqp1);
2828 return pwqp2;
2829 }
2830
2831 if (isl_pw_qpolynomial_is_one(pwqp1)) {
2832 isl_pw_qpolynomial_free(pwqp1);
2833 return pwqp2;
2834 }
2835
2836 if (isl_pw_qpolynomial_is_one(pwqp2)) {
2837 isl_pw_qpolynomial_free(pwqp2);
2838 return pwqp1;
2839 }
2840
2841 n = pwqp1->n * pwqp2->n;
2842 res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n);
2843
2844 for (i = 0; i < pwqp1->n; ++i) {
2845 for (j = 0; j < pwqp2->n; ++j) {
2846 struct isl_set *common;
2847 struct isl_qpolynomial *prod;
2848 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
2849 isl_set_copy(pwqp2->p[j].set));
2850 if (isl_set_plain_is_empty(common)) {
2851 isl_set_free(common);
2852 continue;
2853 }
2854
2855 prod = isl_qpolynomial_mul(
2856 isl_qpolynomial_copy(pwqp1->p[i].qp),
2857 isl_qpolynomial_copy(pwqp2->p[j].qp));
2858
2859 res = isl_pw_qpolynomial_add_piece(res, common, prod);
2860 }
2861 }
2862
2863 isl_pw_qpolynomial_free(pwqp1);
2864 isl_pw_qpolynomial_free(pwqp2);
2865
2866 return res;
2867 error:
2868 isl_pw_qpolynomial_free(pwqp1);
2869 isl_pw_qpolynomial_free(pwqp2);
2870 return NULL;
2871 }
2872
2873 __isl_give isl_val *isl_upoly_eval(__isl_take struct isl_upoly *up,
2874 __isl_take isl_vec *vec)
2875 {
2876 int i;
2877 struct isl_upoly_rec *rec;
2878 isl_val *res;
2879 isl_val *base;
2880
2881 if (isl_upoly_is_cst(up)) {
2882 isl_vec_free(vec);
2883 res = isl_upoly_get_constant_val(up);
2884 isl_upoly_free(up);
2885 return res;
2886 }
2887
2888 rec = isl_upoly_as_rec(up);
2889 if (!rec)
2890 goto error;
2891
2892 isl_assert(up->ctx, rec->n >= 1, goto error);
2893
2894 base = isl_val_rat_from_isl_int(up->ctx,
2895 vec->el[1 + up->var], vec->el[0]);
2896
2897 res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
2898 isl_vec_copy(vec));
2899
2900 for (i = rec->n - 2; i >= 0; --i) {
2901 res = isl_val_mul(res, isl_val_copy(base));
2902 res = isl_val_add(res,
2903 isl_upoly_eval(isl_upoly_copy(rec->p[i]),
2904 isl_vec_copy(vec)));
2905 }
2906
2907 isl_val_free(base);
2908 isl_upoly_free(up);
2909 isl_vec_free(vec);
2910 return res;
2911 error:
2912 isl_upoly_free(up);
2913 isl_vec_free(vec);
2914 return NULL;
2915 }
2916
2917 __isl_give isl_val *isl_qpolynomial_eval(__isl_take isl_qpolynomial *qp,
2918 __isl_take isl_point *pnt)
2919 {
2920 isl_vec *ext;
2921 isl_val *v;
2922
2923 if (!qp || !pnt)
2924 goto error;
2925 isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error);
2926
2927 if (qp->div->n_row == 0)
2928 ext = isl_vec_copy(pnt->vec);
2929 else {
2930 int i;
2931 unsigned dim = isl_space_dim(qp->dim, isl_dim_all);
2932 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
2933 if (!ext)
2934 goto error;
2935
2936 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
2937 for (i = 0; i < qp->div->n_row; ++i) {
2938 isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
2939 1 + dim + i, &ext->el[1+dim+i]);
2940 isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
2941 qp->div->row[i][0]);
2942 }
2943 }
2944
2945 v = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
2946
2947 isl_qpolynomial_free(qp);
2948 isl_point_free(pnt);
2949
2950 return v;
2951 error:
2952 isl_qpolynomial_free(qp);
2953 isl_point_free(pnt);
2954 return NULL;
2955 }
2956
2957 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
2958 __isl_keep struct isl_upoly_cst *cst2)
2959 {
2960 int cmp;
2961 isl_int t;
2962 isl_int_init(t);
2963 isl_int_mul(t, cst1->n, cst2->d);
2964 isl_int_submul(t, cst2->n, cst1->d);
2965 cmp = isl_int_sgn(t);
2966 isl_int_clear(t);
2967 return cmp;
2968 }
2969
2970 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
2971 __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
2972 unsigned first, unsigned n)
2973 {
2974 unsigned total;
2975 unsigned g_pos;
2976 int *exp;
2977
2978 if (!qp)
2979 return NULL;
2980 if (type == isl_dim_out)
2981 isl_die(qp->div->ctx, isl_error_invalid,
2982 "cannot insert output/set dimensions",
2983 goto error);
2984 if (type == isl_dim_in)
2985 type = isl_dim_set;
2986 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2987 return qp;
2988
2989 qp = isl_qpolynomial_cow(qp);
2990 if (!qp)
2991 return NULL;
2992
2993 isl_assert(qp->div->ctx, first <= isl_space_dim(qp->dim, type),
2994 goto error);
2995
2996 g_pos = pos(qp->dim, type) + first;
2997
2998 qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
2999 if (!qp->div)
3000 goto error;
3001
3002 total = qp->div->n_col - 2;
3003 if (total > g_pos) {
3004 int i;
3005 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
3006 if (!exp)
3007 goto error;
3008 for (i = 0; i < total - g_pos; ++i)
3009 exp[i] = i + n;
3010 qp->upoly = expand(qp->upoly, exp, g_pos);
3011 free(exp);
3012 if (!qp->upoly)
3013 goto error;
3014 }
3015
3016 qp->dim = isl_space_insert_dims(qp->dim, type, first, n);
3017 if (!qp->dim)
3018 goto error;
3019
3020 return qp;
3021 error:
3022 isl_qpolynomial_free(qp);
3023 return NULL;
3024 }
3025
3026 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
3027 __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
3028 {
3029 unsigned pos;
3030
3031 pos = isl_qpolynomial_dim(qp, type);
3032
3033 return isl_qpolynomial_insert_dims(qp, type, pos, n);
3034 }
3035
3036 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
3037 __isl_take isl_pw_qpolynomial *pwqp,
3038 enum isl_dim_type type, unsigned n)
3039 {
3040 unsigned pos;
3041
3042 pos = isl_pw_qpolynomial_dim(pwqp, type);
3043
3044 return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
3045 }
3046
3047 static int *reordering_move(isl_ctx *ctx,
3048 unsigned len, unsigned dst, unsigned src, unsigned n)
3049 {
3050 int i;
3051 int *reordering;
3052
3053 reordering = isl_alloc_array(ctx, int, len);
3054 if (!reordering)
3055 return NULL;
3056
3057 if (dst <= src) {
3058 for (i = 0; i < dst; ++i)
3059 reordering[i] = i;
3060 for (i = 0; i < n; ++i)
3061 reordering[src + i] = dst + i;
3062 for (i = 0; i < src - dst; ++i)
3063 reordering[dst + i] = dst + n + i;
3064 for (i = 0; i < len - src - n; ++i)
3065 reordering[src + n + i] = src + n + i;
3066 } else {
3067 for (i = 0; i < src; ++i)
3068 reordering[i] = i;
3069 for (i = 0; i < n; ++i)
3070 reordering[src + i] = dst + i;
3071 for (i = 0; i < dst - src; ++i)
3072 reordering[src + n + i] = src + i;
3073 for (i = 0; i < len - dst - n; ++i)
3074 reordering[dst + n + i] = dst + n + i;
3075 }
3076
3077 return reordering;
3078 }
3079
3080 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
3081 __isl_take isl_qpolynomial *qp,
3082 enum isl_dim_type dst_type, unsigned dst_pos,
3083 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
3084 {
3085 unsigned g_dst_pos;
3086 unsigned g_src_pos;
3087 int *reordering;
3088
3089 qp = isl_qpolynomial_cow(qp);
3090 if (!qp)
3091 return NULL;
3092
3093 if (dst_type == isl_dim_out || src_type == isl_dim_out)
3094 isl_die(qp->dim->ctx, isl_error_invalid,
3095 "cannot move output/set dimension",
3096 goto error);
3097 if (dst_type == isl_dim_in)
3098 dst_type = isl_dim_set;
3099 if (src_type == isl_dim_in)
3100 src_type = isl_dim_set;
3101
3102 isl_assert(qp->dim->ctx, src_pos + n <= isl_space_dim(qp->dim, src_type),
3103 goto error);
3104
3105 g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
3106 g_src_pos = pos(qp->dim, src_type) + src_pos;
3107 if (dst_type > src_type)
3108 g_dst_pos -= n;
3109
3110 qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
3111 if (!qp->div)
3112 goto error;
3113 qp = sort_divs(qp);
3114 if (!qp)
3115 goto error;
3116
3117 reordering = reordering_move(qp->dim->ctx,
3118 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
3119 if (!reordering)
3120 goto error;
3121
3122 qp->upoly = reorder(qp->upoly, reordering);
3123 free(reordering);
3124 if (!qp->upoly)
3125 goto error;
3126
3127 qp->dim = isl_space_move_dims(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
3128 if (!qp->dim)
3129 goto error;
3130
3131 return qp;
3132 error:
3133 isl_qpolynomial_free(qp);
3134 return NULL;
3135 }
3136
3137 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_space *dim,
3138 isl_int *f, isl_int denom)
3139 {
3140 struct isl_upoly *up;
3141
3142 dim = isl_space_domain(dim);
3143 if (!dim)
3144 return NULL;
3145
3146 up = isl_upoly_from_affine(dim->ctx, f, denom,
3147 1 + isl_space_dim(dim, isl_dim_all));
3148
3149 return isl_qpolynomial_alloc(dim, 0, up);
3150 }
3151
3152 __isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
3153 {
3154 isl_ctx *ctx;
3155 struct isl_upoly *up;
3156 isl_qpolynomial *qp;
3157
3158 if (!aff)
3159 return NULL;
3160
3161 ctx = isl_aff_get_ctx(aff);
3162 up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
3163 aff->v->size - 1);
3164
3165 qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff),
3166 aff->ls->div->n_row, up);
3167 if (!qp)
3168 goto error;
3169
3170 isl_mat_free(qp->div);
3171 qp->div = isl_mat_copy(aff->ls->div);
3172 qp->div = isl_mat_cow(qp->div);
3173 if (!qp->div)
3174 goto error;
3175
3176 isl_aff_free(aff);
3177 qp = reduce_divs(qp);
3178 qp = remove_redundant_divs(qp);
3179 return qp;
3180 error:
3181 isl_aff_free(aff);
3182 return NULL;
3183 }
3184
3185 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff(
3186 __isl_take isl_pw_aff *pwaff)
3187 {
3188 int i;
3189 isl_pw_qpolynomial *pwqp;
3190
3191 if (!pwaff)
3192 return NULL;
3193
3194 pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff),
3195 pwaff->n);
3196
3197 for (i = 0; i < pwaff->n; ++i) {
3198 isl_set *dom;
3199 isl_qpolynomial *qp;
3200
3201 dom = isl_set_copy(pwaff->p[i].set);
3202 qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff));
3203 pwqp = isl_pw_qpolynomial_add_piece(pwqp, dom, qp);
3204 }
3205
3206 isl_pw_aff_free(pwaff);
3207 return pwqp;
3208 }
3209
3210 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
3211 __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
3212 {
3213 isl_aff *aff;
3214
3215 aff = isl_constraint_get_bound(c, type, pos);
3216 isl_constraint_free(c);
3217 return isl_qpolynomial_from_aff(aff);
3218 }
3219
3220 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
3221 * in "qp" by subs[i].
3222 */
3223 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
3224 __isl_take isl_qpolynomial *qp,
3225 enum isl_dim_type type, unsigned first, unsigned n,
3226 __isl_keep isl_qpolynomial **subs)
3227 {
3228 int i;
3229 struct isl_upoly **ups;
3230
3231 if (n == 0)
3232 return qp;
3233
3234 qp = isl_qpolynomial_cow(qp);
3235 if (!qp)
3236 return NULL;
3237
3238 if (type == isl_dim_out)
3239 isl_die(qp->dim->ctx, isl_error_invalid,
3240 "cannot substitute output/set dimension",
3241 goto error);
3242 if (type == isl_dim_in)
3243 type = isl_dim_set;
3244
3245 for (i = 0; i < n; ++i)
3246 if (!subs[i])
3247 goto error;
3248
3249 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
3250 goto error);
3251
3252 for (i = 0; i < n; ++i)
3253 isl_assert(qp->dim->ctx, isl_space_is_equal(qp->dim, subs[i]->dim),
3254 goto error);
3255
3256 isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
3257 for (i = 0; i < n; ++i)
3258 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
3259
3260 first += pos(qp->dim, type);
3261
3262 ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
3263 if (!ups)
3264 goto error;
3265 for (i = 0; i < n; ++i)
3266 ups[i] = subs[i]->upoly;
3267
3268 qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
3269
3270 free(ups);
3271
3272 if (!qp->upoly)
3273 goto error;
3274
3275 return qp;
3276 error:
3277 isl_qpolynomial_free(qp);
3278 return NULL;
3279 }
3280
3281 /* Extend "bset" with extra set dimensions for each integer division
3282 * in "qp" and then call "fn" with the extended bset and the polynomial
3283 * that results from replacing each of the integer divisions by the
3284 * corresponding extra set dimension.
3285 */
3286 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
3287 __isl_keep isl_basic_set *bset,
3288 int (*fn)(__isl_take isl_basic_set *bset,
3289 __isl_take isl_qpolynomial *poly, void *user), void *user)
3290 {
3291 isl_space *dim;
3292 isl_mat *div;
3293 isl_qpolynomial *poly;
3294
3295 if (!qp || !bset)
3296 goto error;
3297 if (qp->div->n_row == 0)
3298 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
3299 user);
3300
3301 div = isl_mat_copy(qp->div);
3302 dim = isl_space_copy(qp->dim);
3303 dim = isl_space_add_dims(dim, isl_dim_set, qp->div->n_row);
3304 poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
3305 bset = isl_basic_set_copy(bset);
3306 bset = isl_basic_set_add_dims(bset, isl_dim_set, qp->div->n_row);
3307 bset = add_div_constraints(bset, div);
3308
3309 return fn(bset, poly, user);
3310 error:
3311 return -1;
3312 }
3313
3314 /* Return total degree in variables first (inclusive) up to last (exclusive).
3315 */
3316 int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
3317 {
3318 int deg = -1;
3319 int i;
3320 struct isl_upoly_rec *rec;
3321
3322 if (!up)
3323 return -2;
3324 if (isl_upoly_is_zero(up))
3325 return -1;
3326 if (isl_upoly_is_cst(up) || up->var < first)
3327 return 0;
3328
3329 rec = isl_upoly_as_rec(up);
3330 if (!rec)
3331 return -2;
3332
3333 for (i = 0; i < rec->n; ++i) {
3334 int d;
3335
3336 if (isl_upoly_is_zero(rec->p[i]))
3337 continue;
3338 d = isl_upoly_degree(rec->p[i], first, last);
3339 if (up->var < last)
3340 d += i;
3341 if (d > deg)
3342 deg = d;
3343 }
3344
3345 return deg;
3346 }
3347
3348 /* Return total degree in set variables.
3349 */
3350 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
3351 {
3352 unsigned ovar;
3353 unsigned nvar;
3354
3355 if (!poly)
3356 return -2;
3357
3358 ovar = isl_space_offset(poly->dim, isl_dim_set);
3359 nvar = isl_space_dim(poly->dim, isl_dim_set);
3360 return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
3361 }
3362
3363 __isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
3364 unsigned pos, int deg)
3365 {
3366 int i;
3367 struct isl_upoly_rec *rec;
3368
3369 if (!up)
3370 return NULL;
3371
3372 if (isl_upoly_is_cst(up) || up->var < pos) {
3373 if (deg == 0)
3374 return isl_upoly_copy(up);
3375 else
3376 return isl_upoly_zero(up->ctx);
3377 }
3378
3379 rec = isl_upoly_as_rec(up);
3380 if (!rec)
3381 return NULL;
3382
3383 if (up->var == pos) {
3384 if (deg < rec->n)
3385 return isl_upoly_copy(rec->p[deg]);
3386 else
3387 return isl_upoly_zero(up->ctx);
3388 }
3389
3390 up = isl_upoly_copy(up);
3391 up = isl_upoly_cow(up);
3392 rec = isl_upoly_as_rec(up);
3393 if (!rec)
3394 goto error;
3395
3396 for (i = 0; i < rec->n; ++i) {
3397 struct isl_upoly *t;
3398 t = isl_upoly_coeff(rec->p[i], pos, deg);
3399 if (!t)
3400 goto error;
3401 isl_upoly_free(rec->p[i]);
3402 rec->p[i] = t;
3403 }
3404
3405 return up;
3406 error:
3407 isl_upoly_free(up);
3408 return NULL;
3409 }
3410
3411 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
3412 */
3413 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
3414 __isl_keep isl_qpolynomial *qp,
3415 enum isl_dim_type type, unsigned t_pos, int deg)
3416 {
3417 unsigned g_pos;
3418 struct isl_upoly *up;
3419 isl_qpolynomial *c;
3420
3421 if (!qp)
3422 return NULL;
3423
3424 if (type == isl_dim_out)
3425 isl_die(qp->div->ctx, isl_error_invalid,
3426 "output/set dimension does not have a coefficient",
3427 return NULL);
3428 if (type == isl_dim_in)
3429 type = isl_dim_set;
3430
3431 isl_assert(qp->div->ctx, t_pos < isl_space_dim(qp->dim, type),
3432 return NULL);
3433
3434 g_pos = pos(qp->dim, type) + t_pos;
3435 up = isl_upoly_coeff(qp->upoly, g_pos, deg);
3436
3437 c = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row, up);
3438 if (!c)
3439 return NULL;
3440 isl_mat_free(c->div);
3441 c->div = isl_mat_copy(qp->div);
3442 if (!c->div)
3443 goto error;
3444 return c;
3445 error:
3446 isl_qpolynomial_free(c);
3447 return NULL;
3448 }
3449
3450 /* Homogenize the polynomial in the variables first (inclusive) up to
3451 * last (exclusive) by inserting powers of variable first.
3452 * Variable first is assumed not to appear in the input.
3453 */
3454 __isl_give struct isl_upoly *isl_upoly_homogenize(
3455 __isl_take struct isl_upoly *up, int deg, int target,
3456 int first, int last)
3457 {
3458 int i;
3459 struct isl_upoly_rec *rec;
3460
3461 if (!up)
3462 return NULL;
3463 if (isl_upoly_is_zero(up))
3464 return up;
3465 if (deg == target)
3466 return up;
3467 if (isl_upoly_is_cst(up) || up->var < first) {
3468 struct isl_upoly *hom;
3469
3470 hom = isl_upoly_var_pow(up->ctx, first, target - deg);
3471 if (!hom)
3472 goto error;
3473 rec = isl_upoly_as_rec(hom);
3474 rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
3475
3476 return hom;
3477 }
3478
3479 up = isl_upoly_cow(up);
3480 rec = isl_upoly_as_rec(up);
3481 if (!rec)
3482 goto error;
3483
3484 for (i = 0; i < rec->n; ++i) {
3485 if (isl_upoly_is_zero(rec->p[i]))
3486 continue;
3487 rec->p[i] = isl_upoly_homogenize(rec->p[i],
3488 up->var < last ? deg + i : i, target,
3489 first, last);
3490 if (!rec->p[i])
3491 goto error;
3492 }
3493
3494 return up;
3495 error:
3496 isl_upoly_free(up);
3497 return NULL;
3498 }
3499
3500 /* Homogenize the polynomial in the set variables by introducing
3501 * powers of an extra set variable at position 0.
3502 */
3503 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
3504 __isl_take isl_qpolynomial *poly)
3505 {
3506 unsigned ovar;
3507 unsigned nvar;
3508 int deg = isl_qpolynomial_degree(poly);
3509
3510 if (deg < -1)
3511 goto error;
3512
3513 poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1);
3514 poly = isl_qpolynomial_cow(poly);
3515 if (!poly)
3516 goto error;
3517
3518 ovar = isl_space_offset(poly->dim, isl_dim_set);
3519 nvar = isl_space_dim(poly->dim, isl_dim_set);
3520 poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
3521 ovar, ovar + nvar);
3522 if (!poly->upoly)
3523 goto error;
3524
3525 return poly;
3526 error:
3527 isl_qpolynomial_free(poly);
3528 return NULL;
3529 }
3530
3531 __isl_give isl_term *isl_term_alloc(__isl_take isl_space *dim,
3532 __isl_take isl_mat *div)
3533 {
3534 isl_term *term;
3535 int n;
3536
3537 if (!dim || !div)
3538 goto error;
3539
3540 n = isl_space_dim(dim, isl_dim_all) + div->n_row;
3541
3542 term = isl_calloc(dim->ctx, struct isl_term,
3543 sizeof(struct isl_term) + (n - 1) * sizeof(int));
3544 if (!term)
3545 goto error;
3546
3547 term->ref = 1;
3548 term->dim = dim;
3549 term->div = div;
3550 isl_int_init(term->n);
3551 isl_int_init(term->d);
3552
3553 return term;
3554 error:
3555 isl_space_free(dim);
3556 isl_mat_free(div);
3557 return NULL;
3558 }
3559
3560 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
3561 {
3562 if (!term)
3563 return NULL;
3564
3565 term->ref++;
3566 return term;
3567 }
3568
3569 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
3570 {
3571 int i;
3572 isl_term *dup;
3573 unsigned total;
3574
3575 if (!term)
3576 return NULL;
3577
3578 total = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3579
3580 dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div));
3581 if (!dup)
3582 return NULL;
3583
3584 isl_int_set(dup->n, term->n);
3585 isl_int_set(dup->d, term->d);
3586
3587 for (i = 0; i < total; ++i)
3588 dup->pow[i] = term->pow[i];
3589
3590 return dup;
3591 }
3592
3593 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
3594 {
3595 if (!term)
3596 return NULL;
3597
3598 if (term->ref == 1)
3599 return term;
3600 term->ref--;
3601 return isl_term_dup(term);
3602 }
3603
3604 void isl_term_free(__isl_take isl_term *term)
3605 {
3606 if (!term)
3607 return;
3608
3609 if (--term->ref > 0)
3610 return;
3611
3612 isl_space_free(term->dim);
3613 isl_mat_free(term->div);
3614 isl_int_clear(term->n);
3615 isl_int_clear(term->d);
3616 free(term);
3617 }
3618
3619 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
3620 {
3621 if (!term)
3622 return 0;
3623
3624 switch (type) {
3625 case isl_dim_param:
3626 case isl_dim_in:
3627 case isl_dim_out: return isl_space_dim(term->dim, type);
3628 case isl_dim_div: return term->div->n_row;
3629 case isl_dim_all: return isl_space_dim(term->dim, isl_dim_all) +
3630 term->div->n_row;
3631 default: return 0;
3632 }
3633 }
3634
3635 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
3636 {
3637 return term ? term->dim->ctx : NULL;
3638 }
3639
3640 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
3641 {
3642 if (!term)
3643 return;
3644 isl_int_set(*n, term->n);
3645 }
3646
3647 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
3648 {
3649 if (!term)
3650 return;
3651 isl_int_set(*d, term->d);
3652 }
3653
3654 /* Return the coefficient of the term "term".
3655 */
3656 __isl_give isl_val *isl_term_get_coefficient_val(__isl_keep isl_term *term)
3657 {
3658 if (!term)
3659 return NULL;
3660
3661 return isl_val_rat_from_isl_int(isl_term_get_ctx(term),
3662 term->n, term->d);
3663 }
3664
3665 int isl_term_get_exp(__isl_keep isl_term *term,
3666 enum isl_dim_type type, unsigned pos)
3667 {
3668 if (!term)
3669 return -1;
3670
3671 isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
3672
3673 if (type >= isl_dim_set)
3674 pos += isl_space_dim(term->dim, isl_dim_param);
3675 if (type >= isl_dim_div)
3676 pos += isl_space_dim(term->dim, isl_dim_set);
3677
3678 return term->pow[pos];
3679 }
3680
3681 __isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
3682 {
3683 isl_local_space *ls;
3684 isl_aff *aff;
3685
3686 if (!term)
3687 return NULL;
3688
3689 isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
3690 return NULL);
3691
3692 ls = isl_local_space_alloc_div(isl_space_copy(term->dim),
3693 isl_mat_copy(term->div));
3694 aff = isl_aff_alloc(ls);
3695 if (!aff)
3696 return NULL;
3697
3698 isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size);
3699
3700 aff = isl_aff_normalize(aff);
3701
3702 return aff;
3703 }
3704
3705 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
3706 int (*fn)(__isl_take isl_term *term, void *user),
3707 __isl_take isl_term *term, void *user)
3708 {
3709 int i;
3710 struct isl_upoly_rec *rec;
3711
3712 if (!up || !term)
3713 goto error;
3714
3715 if (isl_upoly_is_zero(up))
3716 return term;
3717
3718 isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
3719 isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
3720 isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
3721
3722 if (isl_upoly_is_cst(up)) {
3723 struct isl_upoly_cst *cst;
3724 cst = isl_upoly_as_cst(up);
3725 if (!cst)
3726 goto error;
3727 term = isl_term_cow(term);
3728 if (!term)
3729 goto error;
3730 isl_int_set(term->n, cst->n);
3731 isl_int_set(term->d, cst->d);
3732 if (fn(isl_term_copy(term), user) < 0)
3733 goto error;
3734 return term;
3735 }
3736
3737 rec = isl_upoly_as_rec(up);
3738 if (!rec)
3739 goto error;
3740
3741 for (i = 0; i < rec->n; ++i) {
3742 term = isl_term_cow(term);
3743 if (!term)
3744 goto error;
3745 term->pow[up->var] = i;
3746 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
3747 if (!term)
3748 goto error;
3749 }
3750 term->pow[up->var] = 0;
3751
3752 return term;
3753 error:
3754 isl_term_free(term);
3755 return NULL;
3756 }
3757
3758 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
3759 int (*fn)(__isl_take isl_term *term, void *user), void *user)
3760 {
3761 isl_term *term;
3762
3763 if (!qp)
3764 return -1;
3765
3766 term = isl_term_alloc(isl_space_copy(qp->dim), isl_mat_copy(qp->div));
3767 if (!term)
3768 return -1;
3769
3770 term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
3771
3772 isl_term_free(term);
3773
3774 return term ? 0 : -1;
3775 }
3776
3777 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
3778 {
3779 struct isl_upoly *up;
3780 isl_qpolynomial *qp;
3781 int i, n;
3782
3783 if (!term)
3784 return NULL;
3785
3786 n = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3787
3788 up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
3789 for (i = 0; i < n; ++i) {
3790 if (!term->pow[i])
3791 continue;
3792 up = isl_upoly_mul(up,
3793 isl_upoly_var_pow(term->dim->ctx, i, term->pow[i]));
3794 }
3795
3796 qp = isl_qpolynomial_alloc(isl_space_copy(term->dim), term->div->n_row, up);
3797 if (!qp)
3798 goto error;
3799 isl_mat_free(qp->div);
3800 qp->div = isl_mat_copy(term->div);
3801 if (!qp->div)
3802 goto error;
3803
3804 isl_term_free(term);
3805 return qp;
3806 error:
3807 isl_qpolynomial_free(qp);
3808 isl_term_free(term);
3809 return NULL;
3810 }
3811
3812 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
3813 __isl_take isl_space *dim)
3814 {
3815 int i;
3816 int extra;
3817 unsigned total;
3818
3819 if (!qp || !dim)
3820 goto error;
3821
3822 if (isl_space_is_equal(qp->dim, dim)) {
3823 isl_space_free(dim);
3824 return qp;
3825 }
3826
3827 qp = isl_qpolynomial_cow(qp);
3828 if (!qp)
3829 goto error;
3830
3831 extra = isl_space_dim(dim, isl_dim_set) -
3832 isl_space_dim(qp->dim, isl_dim_set);
3833 total = isl_space_dim(qp->dim, isl_dim_all);
3834 if (qp->div->n_row) {
3835 int *exp;
3836
3837 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
3838 if (!exp)
3839 goto error;
3840 for (i = 0; i < qp->div->n_row; ++i)
3841 exp[i] = extra + i;
3842 qp->upoly = expand(qp->upoly, exp, total);
3843 free(exp);
3844 if (!qp->upoly)
3845 goto error;
3846 }
3847 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
3848 if (!qp->div)
3849 goto error;
3850 for (i = 0; i < qp->div->n_row; ++i)
3851 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
3852
3853 isl_space_free(qp->dim);
3854 qp->dim = dim;
3855
3856 return qp;
3857 error:
3858 isl_space_free(dim);
3859 isl_qpolynomial_free(qp);
3860 return NULL;
3861 }
3862
3863 /* For each parameter or variable that does not appear in qp,
3864 * first eliminate the variable from all constraints and then set it to zero.
3865 */
3866 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
3867 __isl_keep isl_qpolynomial *qp)
3868 {
3869 int *active = NULL;
3870 int i;
3871 int d;
3872 unsigned nparam;
3873 unsigned nvar;
3874
3875 if (!set || !qp)
3876 goto error;
3877
3878 d = isl_space_dim(set->dim, isl_dim_all);
3879 active = isl_calloc_array(set->ctx, int, d);
3880 if (set_active(qp, active) < 0)
3881 goto error;
3882
3883 for (i = 0; i < d; ++i)
3884 if (!active[i])
3885 break;
3886
3887 if (i == d) {
3888 free(active);
3889 return set;
3890 }
3891
3892 nparam = isl_space_dim(set->dim, isl_dim_param);
3893 nvar = isl_space_dim(set->dim, isl_dim_set);
3894 for (i = 0; i < nparam; ++i) {
3895 if (active[i])
3896 continue;
3897 set = isl_set_eliminate(set, isl_dim_param, i, 1);
3898 set = isl_set_fix_si(set, isl_dim_param, i, 0);
3899 }
3900 for (i = 0; i < nvar; ++i) {
3901 if (active[nparam + i])
3902 continue;
3903 set = isl_set_eliminate(set, isl_dim_set, i, 1);
3904 set = isl_set_fix_si(set, isl_dim_set, i, 0);
3905 }
3906
3907 free(active);
3908
3909 return set;
3910 error:
3911 free(active);
3912 isl_set_free(set);
3913 return NULL;
3914 }
3915
3916 struct isl_opt_data {
3917 isl_qpolynomial *qp;
3918 int first;
3919 isl_val *opt;
3920 int max;
3921 };
3922
3923 static int opt_fn(__isl_take isl_point *pnt, void *user)
3924 {
3925 struct isl_opt_data *data = (struct isl_opt_data *)user;
3926 isl_val *val;
3927
3928 val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
3929 if (data->first) {
3930 data->first = 0;
3931 data->opt = val;
3932 } else if (data->max) {
3933 data->opt = isl_val_max(data->opt, val);
3934 } else {
3935 data->opt = isl_val_min(data->opt, val);
3936 }
3937
3938 return 0;
3939 }
3940
3941 __isl_give isl_val *isl_qpolynomial_opt_on_domain(
3942 __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
3943 {
3944 struct isl_opt_data data = { NULL, 1, NULL, max };
3945
3946 if (!set || !qp)
3947 goto error;
3948
3949 if (isl_upoly_is_cst(qp->upoly)) {
3950 isl_set_free(set);
3951 data.opt = isl_qpolynomial_get_constant_val(qp);
3952 isl_qpolynomial_free(qp);
3953 return data.opt;
3954 }
3955
3956 set = fix_inactive(set, qp);
3957
3958 data.qp = qp;
3959 if (isl_set_foreach_point(set, opt_fn, &data) < 0)
3960 goto error;
3961
3962 if (data.first)
3963 data.opt = isl_val_zero(isl_set_get_ctx(set));
3964
3965 isl_set_free(set);
3966 isl_qpolynomial_free(qp);
3967 return data.opt;
3968 error:
3969 isl_set_free(set);
3970 isl_qpolynomial_free(qp);
3971 isl_val_free(data.opt);
3972 return NULL;
3973 }
3974
3975 __isl_give isl_qpolynomial *isl_qpolynomial_morph_domain(
3976 __isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph)
3977 {
3978 int i;
3979 int n_sub;
3980 isl_ctx *ctx;
3981 struct isl_upoly **subs;
3982 isl_mat *mat, *diag;
3983
3984 qp = isl_qpolynomial_cow(qp);
3985 if (!qp || !morph)
3986 goto error;
3987
3988 ctx = qp->dim->ctx;
3989 isl_assert(ctx, isl_space_is_equal(qp->dim, morph->dom->dim), goto error);
3990
3991 n_sub = morph->inv->n_row - 1;
3992 if (morph->inv->n_row != morph->inv->n_col)
3993 n_sub += qp->div->n_row;
3994 subs = isl_calloc_array(ctx, struct isl_upoly *, n_sub);
3995 if (!subs)
3996 goto error;
3997
3998 for (i = 0; 1 + i < morph->inv->n_row; ++i)
3999 subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
4000 morph->inv->row[0][0], morph->inv->n_col);
4001 if (morph->inv->n_row != morph->inv->n_col)
4002 for (i = 0; i < qp->div->n_row; ++i)
4003 subs[morph->inv->n_row - 1 + i] =
4004 isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
4005
4006 qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
4007
4008 for (i = 0; i < n_sub; ++i)
4009 isl_upoly_free(subs[i]);
4010 free(subs);
4011
4012 diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
4013 mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
4014 diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
4015 mat = isl_mat_diagonal(mat, diag);
4016 qp->div = isl_mat_product(qp->div, mat);
4017 isl_space_free(qp->dim);
4018 qp->dim = isl_space_copy(morph->ran->dim);
4019
4020 if (!qp->upoly || !qp->div || !qp->dim)
4021 goto error;
4022
4023 isl_morph_free(morph);
4024
4025 return qp;
4026 error:
4027 isl_qpolynomial_free(qp);
4028 isl_morph_free(morph);
4029 return NULL;
4030 }
4031
4032 static int neg_entry(void **entry, void *user)
4033 {
4034 isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
4035
4036 *pwqp = isl_pw_qpolynomial_neg(*pwqp);
4037
4038 return *pwqp ? 0 : -1;
4039 }
4040
4041 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_neg(
4042 __isl_take isl_union_pw_qpolynomial *upwqp)
4043 {
4044 upwqp = isl_union_pw_qpolynomial_cow(upwqp);
4045 if (!upwqp)
4046 return NULL;
4047
4048 if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
4049 &neg_entry, NULL) < 0)
4050 goto error;
4051
4052 return upwqp;
4053 error:
4054 isl_union_pw_qpolynomial_free(upwqp);
4055 return NULL;
4056 }
4057
4058 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
4059 __isl_take isl_union_pw_qpolynomial *upwqp1,
4060 __isl_take isl_union_pw_qpolynomial *upwqp2)
4061 {
4062 return match_bin_op(upwqp1, upwqp2, &isl_pw_qpolynomial_mul);
4063 }
4064
4065 /* Reorder the columns of the given div definitions according to the
4066 * given reordering.
4067 */
4068 static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
4069 __isl_take isl_reordering *r)
4070 {
4071 int i, j;
4072 isl_mat *mat;
4073 int extra;
4074
4075 if (!div || !r)
4076 goto error;
4077
4078 extra = isl_space_dim(r->dim, isl_dim_all) + div->n_row - r->len;
4079 mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
4080 if (!mat)
4081 goto error;
4082
4083 for (i = 0; i < div->n_row; ++i) {
4084 isl_seq_cpy(mat->row[i], div->row[i], 2);
4085 isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
4086 for (j = 0; j < r->len; ++j)
4087 isl_int_set(mat->row[i][2 + r->pos[j]],
4088 div->row[i][2 + j]);
4089 }
4090
4091 isl_reordering_free(r);
4092 isl_mat_free(div);
4093 return mat;
4094 error:
4095 isl_reordering_free(r);
4096 isl_mat_free(div);
4097 return NULL;
4098 }
4099
4100 /* Reorder the dimension of "qp" according to the given reordering.
4101 */
4102 __isl_give isl_qpolynomial *isl_qpolynomial_realign_domain(
4103 __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
4104 {
4105 qp = isl_qpolynomial_cow(qp);
4106 if (!qp)
4107 goto error;
4108
4109 r = isl_reordering_extend(r, qp->div->n_row);
4110 if (!r)
4111 goto error;
4112
4113 qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
4114 if (!qp->div)
4115 goto error;
4116
4117 qp->upoly = reorder(qp->upoly, r->pos);
4118 if (!qp->upoly)
4119 goto error;
4120
4121 qp = isl_qpolynomial_reset_domain_space(qp, isl_space_copy(r->dim));
4122
4123 isl_reordering_free(r);
4124 return qp;
4125 error:
4126 isl_qpolynomial_free(qp);
4127 isl_reordering_free(r);
4128 return NULL;
4129 }
4130
4131 __isl_give isl_qpolynomial *isl_qpolynomial_align_params(
4132 __isl_take isl_qpolynomial *qp, __isl_take isl_space *model)
4133 {
4134 if (!qp || !model)
4135 goto error;
4136
4137 if (!isl_space_match(qp->dim, isl_dim_param, model, isl_dim_param)) {
4138 isl_reordering *exp;
4139
4140 model = isl_space_drop_dims(model, isl_dim_in,
4141 0, isl_space_dim(model, isl_dim_in));
4142 model = isl_space_drop_dims(model, isl_dim_out,
4143 0, isl_space_dim(model, isl_dim_out));
4144 exp = isl_parameter_alignment_reordering(qp->dim, model);
4145 exp = isl_reordering_extend_space(exp,
4146 isl_qpolynomial_get_domain_space(qp));
4147 qp = isl_qpolynomial_realign_domain(qp, exp);
4148 }
4149
4150 isl_space_free(model);
4151 return qp;
4152 error:
4153 isl_space_free(model);
4154 isl_qpolynomial_free(qp);
4155 return NULL;
4156 }
4157
4158 struct isl_split_periods_data {
4159 int max_periods;
4160 isl_pw_qpolynomial *res;
4161 };
4162
4163 /* Create a slice where the integer division "div" has the fixed value "v".
4164 * In particular, if "div" refers to floor(f/m), then create a slice
4165 *
4166 * m v <= f <= m v + (m - 1)
4167 *
4168 * or
4169 *
4170 * f - m v >= 0
4171 * -f + m v + (m - 1) >= 0
4172 */
4173 static __isl_give isl_set *set_div_slice(__isl_take isl_space *dim,
4174 __isl_keep isl_qpolynomial *qp, int div, isl_int v)
4175 {
4176 int total;
4177 isl_basic_set *bset = NULL;
4178 int k;
4179
4180 if (!dim || !qp)
4181 goto error;
4182
4183 total = isl_space_dim(dim, isl_dim_all);
4184 bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0, 0, 2);
4185
4186 k = isl_basic_set_alloc_inequality(bset);
4187 if (k < 0)
4188 goto error;
4189 isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4190 isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
4191
4192 k = isl_basic_set_alloc_inequality(bset);
4193 if (k < 0)
4194 goto error;
4195 isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4196 isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
4197 isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
4198 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4199
4200 isl_space_free(dim);
4201 return isl_set_from_basic_set(bset);
4202 error:
4203 isl_basic_set_free(bset);
4204 isl_space_free(dim);
4205 return NULL;
4206 }
4207
4208 static int split_periods(__isl_take isl_set *set,
4209 __isl_take isl_qpolynomial *qp, void *user);
4210
4211 /* Create a slice of the domain "set" such that integer division "div"
4212 * has the fixed value "v" and add the results to data->res,
4213 * replacing the integer division by "v" in "qp".
4214 */
4215 static int set_div(__isl_take isl_set *set,
4216 __isl_take isl_qpolynomial *qp, int div, isl_int v,
4217 struct isl_split_periods_data *data)
4218 {
4219 int i;
4220 int total;
4221 isl_set *slice;
4222 struct isl_upoly *cst;
4223
4224 slice = set_div_slice(isl_set_get_space(set), qp, div, v);
4225 set = isl_set_intersect(set, slice);
4226
4227 if (!qp)
4228 goto error;
4229
4230 total = isl_space_dim(qp->dim, isl_dim_all);
4231
4232 for (i = div + 1; i < qp->div->n_row; ++i) {
4233 if (isl_int_is_zero(qp->div->row[i][2 + total + div]))
4234 continue;
4235 isl_int_addmul(qp->div->row[i][1],
4236 qp->div->row[i][2 + total + div], v);
4237 isl_int_set_si(qp->div->row[i][2 + total + div], 0);
4238 }
4239
4240 cst = isl_upoly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
4241 qp = substitute_div(qp, div, cst);
4242
4243 return split_periods(set, qp, data);
4244 error:
4245 isl_set_free(set);
4246 isl_qpolynomial_free(qp);
4247 return -1;
4248 }
4249
4250 /* Split the domain "set" such that integer division "div"
4251 * has a fixed value (ranging from "min" to "max") on each slice
4252 * and add the results to data->res.
4253 */
4254 static int split_div(__isl_take isl_set *set,
4255 __isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
4256 struct isl_split_periods_data *data)
4257 {
4258 for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
4259 isl_set *set_i = isl_set_copy(set);
4260 isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
4261
4262 if (set_div(set_i, qp_i, div, min, data) < 0)
4263 goto error;
4264 }
4265 isl_set_free(set);
4266 isl_qpolynomial_free(qp);
4267 return 0;
4268 error:
4269 isl_set_free(set);
4270 isl_qpolynomial_free(qp);
4271 return -1;
4272 }
4273
4274 /* If "qp" refers to any integer division
4275 * that can only attain "max_periods" distinct values on "set"
4276 * then split the domain along those distinct values.
4277 * Add the results (or the original if no splitting occurs)
4278 * to data->res.
4279 */
4280 static int split_periods(__isl_take isl_set *set,
4281 __isl_take isl_qpolynomial *qp, void *user)
4282 {
4283 int i;
4284 isl_pw_qpolynomial *pwqp;
4285 struct isl_split_periods_data *data;
4286 isl_int min, max;
4287 int total;
4288 int r = 0;
4289
4290 data = (struct isl_split_periods_data *)user;
4291
4292 if (!set || !qp)
4293 goto error;
4294
4295 if (qp->div->n_row == 0) {
4296 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4297 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4298 return 0;
4299 }
4300
4301 isl_int_init(min);
4302 isl_int_init(max);
4303 total = isl_space_dim(qp->dim, isl_dim_all);
4304 for (i = 0; i < qp->div->n_row; ++i) {
4305 enum isl_lp_result lp_res;
4306
4307 if (isl_seq_first_non_zero(qp->div->row[i] + 2 + total,
4308 qp->div->n_row) != -1)
4309 continue;
4310
4311 lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
4312 set->ctx->one, &min, NULL, NULL);
4313 if (lp_res == isl_lp_error)
4314 goto error2;
4315 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4316 continue;
4317 isl_int_fdiv_q(min, min, qp->div->row[i][0]);
4318
4319 lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
4320 set->ctx->one, &max, NULL, NULL);
4321 if (lp_res == isl_lp_error)
4322 goto error2;
4323 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4324 continue;
4325 isl_int_fdiv_q(max, max, qp->div->row[i][0]);
4326
4327 isl_int_sub(max, max, min);
4328 if (isl_int_cmp_si(max, data->max_periods) < 0) {
4329 isl_int_add(max, max, min);
4330 break;
4331 }
4332 }
4333
4334 if (i < qp->div->n_row) {
4335 r = split_div(set, qp, i, min, max, data);
4336 } else {
4337 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4338 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4339 }
4340
4341 isl_int_clear(max);
4342 isl_int_clear(min);
4343
4344 return r;
4345 error2:
4346 isl_int_clear(max);
4347 isl_int_clear(min);
4348 error:
4349 isl_set_free(set);
4350 isl_qpolynomial_free(qp);
4351 return -1;
4352 }
4353
4354 /* If any quasi-polynomial in pwqp refers to any integer division
4355 * that can only attain "max_periods" distinct values on its domain
4356 * then split the domain along those distinct values.
4357 */
4358 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
4359 __isl_take isl_pw_qpolynomial *pwqp, int max_periods)
4360 {
4361 struct isl_split_periods_data data;
4362
4363 data.max_periods = max_periods;
4364 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4365
4366 if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
4367 goto error;
4368
4369 isl_pw_qpolynomial_free(pwqp);
4370
4371 return data.res;
4372 error:
4373 isl_pw_qpolynomial_free(data.res);
4374 isl_pw_qpolynomial_free(pwqp);
4375 return NULL;
4376 }
4377
4378 /* Construct a piecewise quasipolynomial that is constant on the given
4379 * domain. In particular, it is
4380 * 0 if cst == 0
4381 * 1 if cst == 1
4382 * infinity if cst == -1
4383 */
4384 static __isl_give isl_pw_qpolynomial *constant_on_domain(
4385 __isl_take isl_basic_set *bset, int cst)
4386 {
4387 isl_space *dim;
4388 isl_qpolynomial *qp;
4389
4390 if (!bset)
4391 return NULL;
4392
4393 bset = isl_basic_set_params(bset);
4394 dim = isl_basic_set_get_space(bset);
4395 if (cst < 0)
4396 qp = isl_qpolynomial_infty_on_domain(dim);
4397 else if (cst == 0)
4398 qp = isl_qpolynomial_zero_on_domain(dim);
4399 else
4400 qp = isl_qpolynomial_one_on_domain(dim);
4401 return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
4402 }
4403
4404 /* Factor bset, call fn on each of the factors and return the product.
4405 *
4406 * If no factors can be found, simply call fn on the input.
4407 * Otherwise, construct the factors based on the factorizer,
4408 * call fn on each factor and compute the product.
4409 */
4410 static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
4411 __isl_take isl_basic_set *bset,
4412 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4413 {
4414 int i, n;
4415 isl_space *dim;
4416 isl_set *set;
4417 isl_factorizer *f;
4418 isl_qpolynomial *qp;
4419 isl_pw_qpolynomial *pwqp;
4420 unsigned nparam;
4421 unsigned nvar;
4422
4423 f = isl_basic_set_factorizer(bset);
4424 if (!f)
4425 goto error;
4426 if (f->n_group == 0) {
4427 isl_factorizer_free(f);
4428 return fn(bset);
4429 }
4430
4431 nparam = isl_basic_set_dim(bset, isl_dim_param);
4432 nvar = isl_basic_set_dim(bset, isl_dim_set);
4433
4434 dim = isl_basic_set_get_space(bset);
4435 dim = isl_space_domain(dim);
4436 set = isl_set_universe(isl_space_copy(dim));
4437 qp = isl_qpolynomial_one_on_domain(dim);
4438 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4439
4440 bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
4441
4442 for (i = 0, n = 0; i < f->n_group; ++i) {
4443 isl_basic_set *bset_i;
4444 isl_pw_qpolynomial *pwqp_i;
4445
4446 bset_i = isl_basic_set_copy(bset);
4447 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4448 nparam + n + f->len[i], nvar - n - f->len[i]);
4449 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4450 nparam, n);
4451 bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
4452 n + f->len[i], nvar - n - f->len[i]);
4453 bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
4454
4455 pwqp_i = fn(bset_i);
4456 pwqp = isl_pw_qpolynomial_mul(pwqp, pwqp_i);
4457
4458 n += f->len[i];
4459 }
4460
4461 isl_basic_set_free(bset);
4462 isl_factorizer_free(f);
4463
4464 return pwqp;
4465 error:
4466 isl_basic_set_free(bset);
4467 return NULL;
4468 }
4469
4470 /* Factor bset, call fn on each of the factors and return the product.
4471 * The function is assumed to evaluate to zero on empty domains,
4472 * to one on zero-dimensional domains and to infinity on unbounded domains
4473 * and will not be called explicitly on zero-dimensional or unbounded domains.
4474 *
4475 * We first check for some special cases and remove all equalities.
4476 * Then we hand over control to compressed_multiplicative_call.
4477 */
4478 __isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
4479 __isl_take isl_basic_set *bset,
4480 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4481 {
4482 int bounded;
4483 isl_morph *morph;
4484 isl_pw_qpolynomial *pwqp;
4485
4486 if (!bset)
4487 return NULL;
4488
4489 if (isl_basic_set_plain_is_empty(bset))
4490 return constant_on_domain(bset, 0);
4491
4492 if (isl_basic_set_dim(bset, isl_dim_set) == 0)
4493 return constant_on_domain(bset, 1);
4494
4495 bounded = isl_basic_set_is_bounded(bset);
4496 if (bounded < 0)
4497 goto error;
4498 if (!bounded)
4499 return constant_on_domain(bset, -1);
4500
4501 if (bset->n_eq == 0)
4502 return compressed_multiplicative_call(bset, fn);
4503
4504 morph = isl_basic_set_full_compression(bset);
4505 bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
4506
4507 pwqp = compressed_multiplicative_call(bset, fn);
4508
4509 morph = isl_morph_dom_params(morph);
4510 morph = isl_morph_ran_params(morph);
4511 morph = isl_morph_inverse(morph);
4512
4513 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph);
4514
4515 return pwqp;
4516 error:
4517 isl_basic_set_free(bset);
4518 return NULL;
4519 }
4520
4521 /* Drop all floors in "qp", turning each integer division [a/m] into
4522 * a rational division a/m. If "down" is set, then the integer division
4523 * is replaced by (a-(m-1))/m instead.
4524 */
4525 static __isl_give isl_qpolynomial *qp_drop_floors(
4526 __isl_take isl_qpolynomial *qp, int down)
4527 {
4528 int i;
4529 struct isl_upoly *s;
4530
4531 if (!qp)
4532 return NULL;
4533 if (qp->div->n_row == 0)
4534 return qp;
4535
4536 qp = isl_qpolynomial_cow(qp);
4537 if (!qp)
4538 return NULL;
4539
4540 for (i = qp->div->n_row - 1; i >= 0; --i) {
4541 if (down) {
4542 isl_int_sub(qp->div->row[i][1],
4543 qp->div->row[i][1], qp->div->row[i][0]);
4544 isl_int_add_ui(qp->div->row[i][1],
4545 qp->div->row[i][1], 1);
4546 }
4547 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
4548 qp->div->row[i][0], qp->div->n_col - 1);
4549 qp = substitute_div(qp, i, s);
4550 if (!qp)
4551 return NULL;
4552 }
4553
4554 return qp;
4555 }
4556
4557 /* Drop all floors in "pwqp", turning each integer division [a/m] into
4558 * a rational division a/m.
4559 */
4560 static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
4561 __isl_take isl_pw_qpolynomial *pwqp)
4562 {
4563 int i;
4564
4565 if (!pwqp)
4566 return NULL;
4567
4568 if (isl_pw_qpolynomial_is_zero(pwqp))
4569 return pwqp;
4570
4571 pwqp = isl_pw_qpolynomial_cow(pwqp);
4572 if (!pwqp)
4573 return NULL;
4574
4575 for (i = 0; i < pwqp->n; ++i) {
4576 pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
4577 if (!pwqp->p[i].qp)
4578 goto error;
4579 }
4580
4581 return pwqp;
4582 error:
4583 isl_pw_qpolynomial_free(pwqp);
4584 return NULL;
4585 }
4586
4587 /* Adjust all the integer divisions in "qp" such that they are at least
4588 * one over the given orthant (identified by "signs"). This ensures
4589 * that they will still be non-negative even after subtracting (m-1)/m.
4590 *
4591 * In particular, f is replaced by f' + v, changing f = [a/m]
4592 * to f' = [(a - m v)/m].
4593 * If the constant term k in a is smaller than m,
4594 * the constant term of v is set to floor(k/m) - 1.
4595 * For any other term, if the coefficient c and the variable x have
4596 * the same sign, then no changes are needed.
4597 * Otherwise, if the variable is positive (and c is negative),
4598 * then the coefficient of x in v is set to floor(c/m).
4599 * If the variable is negative (and c is positive),
4600 * then the coefficient of x in v is set to ceil(c/m).
4601 */
4602 static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
4603 int *signs)
4604 {
4605 int i, j;
4606 int total;
4607 isl_vec *v = NULL;
4608 struct isl_upoly *s;
4609
4610 qp = isl_qpolynomial_cow(qp);
4611 if (!qp)
4612 return NULL;
4613 qp->div = isl_mat_cow(qp->div);
4614 if (!qp->div)
4615 goto error;
4616
4617 total = isl_space_dim(qp->dim, isl_dim_all);
4618 v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
4619
4620 for (i = 0; i < qp->div->n_row; ++i) {
4621 isl_int *row = qp->div->row[i];
4622 v = isl_vec_clr(v);
4623 if (!v)
4624 goto error;
4625 if (isl_int_lt(row[1], row[0])) {
4626 isl_int_fdiv_q(v->el[0], row[1], row[0]);
4627 isl_int_sub_ui(v->el[0], v->el[0], 1);
4628 isl_int_submul(row[1], row[0], v->el[0]);
4629 }
4630 for (j = 0; j < total; ++j) {
4631 if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
4632 continue;
4633 if (signs[j] < 0)
4634 isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
4635 else
4636 isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
4637 isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
4638 }
4639 for (j = 0; j < i; ++j) {
4640 if (isl_int_sgn(row[2 + total + j]) >= 0)
4641 continue;
4642 isl_int_fdiv_q(v->el[1 + total + j],
4643 row[2 + total + j], row[0]);
4644 isl_int_submul(row[2 + total + j],
4645 row[0], v->el[1 + total + j]);
4646 }
4647 for (j = i + 1; j < qp->div->n_row; ++j) {
4648 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
4649 continue;
4650 isl_seq_combine(qp->div->row[j] + 1,
4651 qp->div->ctx->one, qp->div->row[j] + 1,
4652 qp->div->row[j][2 + total + i], v->el, v->size);
4653 }
4654 isl_int_set_si(v->el[1 + total + i], 1);
4655 s = isl_upoly_from_affine(qp->dim->ctx, v->el,
4656 qp->div->ctx->one, v->size);
4657 qp->upoly = isl_upoly_subs(qp->upoly, total + i, 1, &s);
4658 isl_upoly_free(s);
4659 if (!qp->upoly)
4660 goto error;
4661 }
4662
4663 isl_vec_free(v);
4664 return qp;
4665 error:
4666 isl_vec_free(v);
4667 isl_qpolynomial_free(qp);
4668 return NULL;
4669 }
4670
4671 struct isl_to_poly_data {
4672 int sign;
4673 isl_pw_qpolynomial *res;
4674 isl_qpolynomial *qp;
4675 };
4676
4677 /* Appoximate data->qp by a polynomial on the orthant identified by "signs".
4678 * We first make all integer divisions positive and then split the
4679 * quasipolynomials into terms with sign data->sign (the direction
4680 * of the requested approximation) and terms with the opposite sign.
4681 * In the first set of terms, each integer division [a/m] is
4682 * overapproximated by a/m, while in the second it is underapproximated
4683 * by (a-(m-1))/m.
4684 */
4685 static int to_polynomial_on_orthant(__isl_take isl_set *orthant, int *signs,
4686 void *user)
4687 {
4688 struct isl_to_poly_data *data = user;
4689 isl_pw_qpolynomial *t;
4690 isl_qpolynomial *qp, *up, *down;
4691
4692 qp = isl_qpolynomial_copy(data->qp);
4693 qp = make_divs_pos(qp, signs);
4694
4695 up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
4696 up = qp_drop_floors(up, 0);
4697 down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
4698 down = qp_drop_floors(down, 1);
4699
4700 isl_qpolynomial_free(qp);
4701 qp = isl_qpolynomial_add(up, down);
4702
4703 t = isl_pw_qpolynomial_alloc(orthant, qp);
4704 data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
4705
4706 return 0;
4707 }
4708
4709 /* Approximate each quasipolynomial by a polynomial. If "sign" is positive,
4710 * the polynomial will be an overapproximation. If "sign" is negative,
4711 * it will be an underapproximation. If "sign" is zero, the approximation
4712 * will lie somewhere in between.
4713 *
4714 * In particular, is sign == 0, we simply drop the floors, turning
4715 * the integer divisions into rational divisions.
4716 * Otherwise, we split the domains into orthants, make all integer divisions
4717 * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
4718 * depending on the requested sign and the sign of the term in which
4719 * the integer division appears.
4720 */
4721 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
4722 __isl_take isl_pw_qpolynomial *pwqp, int sign)
4723 {
4724 int i;
4725 struct isl_to_poly_data data;
4726
4727 if (sign == 0)
4728 return pwqp_drop_floors(pwqp);
4729
4730 if (!pwqp)
4731 return NULL;
4732
4733 data.sign = sign;
4734 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4735
4736 for (i = 0; i < pwqp->n; ++i) {
4737 if (pwqp->p[i].qp->div->n_row == 0) {
4738 isl_pw_qpolynomial *t;
4739 t = isl_pw_qpolynomial_alloc(
4740 isl_set_copy(pwqp->p[i].set),
4741 isl_qpolynomial_copy(pwqp->p[i].qp));
4742 data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
4743 continue;
4744 }
4745 data.qp = pwqp->p[i].qp;
4746 if (isl_set_foreach_orthant(pwqp->p[i].set,
4747 &to_polynomial_on_orthant, &data) < 0)
4748 goto error;
4749 }
4750
4751 isl_pw_qpolynomial_free(pwqp);
4752
4753 return data.res;
4754 error:
4755 isl_pw_qpolynomial_free(pwqp);
4756 isl_pw_qpolynomial_free(data.res);
4757 return NULL;
4758 }
4759
4760 static int poly_entry(void **entry, void *user)
4761 {
4762 int *sign = user;
4763 isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
4764
4765 *pwqp = isl_pw_qpolynomial_to_polynomial(*pwqp, *sign);
4766
4767 return *pwqp ? 0 : -1;
4768 }
4769
4770 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
4771 __isl_take isl_union_pw_qpolynomial *upwqp, int sign)
4772 {
4773 upwqp = isl_union_pw_qpolynomial_cow(upwqp);
4774 if (!upwqp)
4775 return NULL;
4776
4777 if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
4778 &poly_entry, &sign) < 0)
4779 goto error;
4780
4781 return upwqp;
4782 error:
4783 isl_union_pw_qpolynomial_free(upwqp);
4784 return NULL;
4785 }
4786
4787 __isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
4788 __isl_take isl_qpolynomial *qp)
4789 {
4790 int i, k;
4791 isl_space *dim;
4792 isl_vec *aff = NULL;
4793 isl_basic_map *bmap = NULL;
4794 unsigned pos;
4795 unsigned n_div;
4796
4797 if (!qp)
4798 return NULL;
4799 if (!isl_upoly_is_affine(qp->upoly))
4800 isl_die(qp->dim->ctx, isl_error_invalid,
4801 "input quasi-polynomial not affine", goto error);
4802 aff = isl_qpolynomial_extract_affine(qp);
4803 if (!aff)
4804 goto error;
4805 dim = isl_qpolynomial_get_space(qp);
4806 pos = 1 + isl_space_offset(dim, isl_dim_out);
4807 n_div = qp->div->n_row;
4808 bmap = isl_basic_map_alloc_space(dim, n_div, 1, 2 * n_div);
4809
4810 for (i = 0; i < n_div; ++i) {
4811 k = isl_basic_map_alloc_div(bmap);
4812 if (k < 0)
4813 goto error;
4814 isl_seq_cpy(bmap->div[k], qp->div->row[i], qp->div->n_col);
4815 isl_int_set_si(bmap->div[k][qp->div->n_col], 0);
4816 if (isl_basic_map_add_div_constraints(bmap, k) < 0)
4817 goto error;
4818 }
4819 k = isl_basic_map_alloc_equality(bmap);
4820 if (k < 0)
4821 goto error;
4822 isl_int_neg(bmap->eq[k][pos], aff->el[0]);
4823 isl_seq_cpy(bmap->eq[k], aff->el + 1, pos);
4824 isl_seq_cpy(bmap->eq[k] + pos + 1, aff->el + 1 + pos, n_div);
4825
4826 isl_vec_free(aff);
4827 isl_qpolynomial_free(qp);
4828 bmap = isl_basic_map_finalize(bmap);
4829 return bmap;
4830 error:
4831 isl_vec_free(aff);
4832 isl_qpolynomial_free(qp);
4833 isl_basic_map_free(bmap);
4834 return NULL;
4835 }
Integer set library