?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
isl_schedule_constraints: split proximity constraints into coincidence/proximity
[isl.git] / isl_polynomial.c
blob1058e8fb41c73fff560d0195938b68868529f3cd
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 int n_div1, n_div2;
1374
1375 qp1 = isl_qpolynomial_cow(qp1);
1376 qp2 = isl_qpolynomial_cow(qp2);
1377
1378 if (!qp1 || !qp2)
1379 goto error;
1380
1381 isl_assert(qp1->div->ctx, qp1->div->n_row >= qp2->div->n_row &&
1382 qp1->div->n_col >= qp2->div->n_col, goto error);
1383
1384 n_div1 = qp1->div->n_row;
1385 n_div2 = qp2->div->n_row;
1386 exp1 = isl_alloc_array(qp1->div->ctx, int, n_div1);
1387 exp2 = isl_alloc_array(qp2->div->ctx, int, n_div2);
1388 if ((n_div1 && !exp1) || (n_div2 && !exp2))
1389 goto error;
1390
1391 div = isl_merge_divs(qp1->div, qp2->div, exp1, exp2);
1392 if (!div)
1393 goto error;
1394
1395 isl_mat_free(qp1->div);
1396 qp1->div = isl_mat_copy(div);
1397 isl_mat_free(qp2->div);
1398 qp2->div = isl_mat_copy(div);
1399
1400 qp1->upoly = expand(qp1->upoly, exp1, div->n_col - div->n_row - 2);
1401 qp2->upoly = expand(qp2->upoly, exp2, div->n_col - div->n_row - 2);
1402
1403 if (!qp1->upoly || !qp2->upoly)
1404 goto error;
1405
1406 isl_mat_free(div);
1407 free(exp1);
1408 free(exp2);
1409
1410 return fn(qp1, qp2);
1411 error:
1412 isl_mat_free(div);
1413 free(exp1);
1414 free(exp2);
1415 isl_qpolynomial_free(qp1);
1416 isl_qpolynomial_free(qp2);
1417 return NULL;
1418 }
1419
1420 __isl_give isl_qpolynomial *isl_qpolynomial_add(__isl_take isl_qpolynomial *qp1,
1421 __isl_take isl_qpolynomial *qp2)
1422 {
1423 qp1 = isl_qpolynomial_cow(qp1);
1424
1425 if (!qp1 || !qp2)
1426 goto error;
1427
1428 if (qp1->div->n_row < qp2->div->n_row)
1429 return isl_qpolynomial_add(qp2, qp1);
1430
1431 isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1432 if (!compatible_divs(qp1->div, qp2->div))
1433 return with_merged_divs(isl_qpolynomial_add, qp1, qp2);
1434
1435 qp1->upoly = isl_upoly_sum(qp1->upoly, isl_upoly_copy(qp2->upoly));
1436 if (!qp1->upoly)
1437 goto error;
1438
1439 isl_qpolynomial_free(qp2);
1440
1441 return qp1;
1442 error:
1443 isl_qpolynomial_free(qp1);
1444 isl_qpolynomial_free(qp2);
1445 return NULL;
1446 }
1447
1448 __isl_give isl_qpolynomial *isl_qpolynomial_add_on_domain(
1449 __isl_keep isl_set *dom,
1450 __isl_take isl_qpolynomial *qp1,
1451 __isl_take isl_qpolynomial *qp2)
1452 {
1453 qp1 = isl_qpolynomial_add(qp1, qp2);
1454 qp1 = isl_qpolynomial_gist(qp1, isl_set_copy(dom));
1455 return qp1;
1456 }
1457
1458 __isl_give isl_qpolynomial *isl_qpolynomial_sub(__isl_take isl_qpolynomial *qp1,
1459 __isl_take isl_qpolynomial *qp2)
1460 {
1461 return isl_qpolynomial_add(qp1, isl_qpolynomial_neg(qp2));
1462 }
1463
1464 __isl_give isl_qpolynomial *isl_qpolynomial_add_isl_int(
1465 __isl_take isl_qpolynomial *qp, isl_int v)
1466 {
1467 if (isl_int_is_zero(v))
1468 return qp;
1469
1470 qp = isl_qpolynomial_cow(qp);
1471 if (!qp)
1472 return NULL;
1473
1474 qp->upoly = isl_upoly_add_isl_int(qp->upoly, v);
1475 if (!qp->upoly)
1476 goto error;
1477
1478 return qp;
1479 error:
1480 isl_qpolynomial_free(qp);
1481 return NULL;
1482
1483 }
1484
1485 __isl_give isl_qpolynomial *isl_qpolynomial_neg(__isl_take isl_qpolynomial *qp)
1486 {
1487 if (!qp)
1488 return NULL;
1489
1490 return isl_qpolynomial_mul_isl_int(qp, qp->dim->ctx->negone);
1491 }
1492
1493 __isl_give isl_qpolynomial *isl_qpolynomial_mul_isl_int(
1494 __isl_take isl_qpolynomial *qp, isl_int v)
1495 {
1496 if (isl_int_is_one(v))
1497 return qp;
1498
1499 if (qp && isl_int_is_zero(v)) {
1500 isl_qpolynomial *zero;
1501 zero = isl_qpolynomial_zero_on_domain(isl_space_copy(qp->dim));
1502 isl_qpolynomial_free(qp);
1503 return zero;
1504 }
1505
1506 qp = isl_qpolynomial_cow(qp);
1507 if (!qp)
1508 return NULL;
1509
1510 qp->upoly = isl_upoly_mul_isl_int(qp->upoly, v);
1511 if (!qp->upoly)
1512 goto error;
1513
1514 return qp;
1515 error:
1516 isl_qpolynomial_free(qp);
1517 return NULL;
1518 }
1519
1520 __isl_give isl_qpolynomial *isl_qpolynomial_scale(
1521 __isl_take isl_qpolynomial *qp, isl_int v)
1522 {
1523 return isl_qpolynomial_mul_isl_int(qp, v);
1524 }
1525
1526 /* Multiply "qp" by "v".
1527 */
1528 __isl_give isl_qpolynomial *isl_qpolynomial_scale_val(
1529 __isl_take isl_qpolynomial *qp, __isl_take isl_val *v)
1530 {
1531 if (!qp || !v)
1532 goto error;
1533
1534 if (!isl_val_is_rat(v))
1535 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
1536 "expecting rational factor", goto error);
1537
1538 if (isl_val_is_one(v)) {
1539 isl_val_free(v);
1540 return qp;
1541 }
1542
1543 if (isl_val_is_zero(v)) {
1544 isl_space *space;
1545
1546 space = isl_qpolynomial_get_domain_space(qp);
1547 isl_qpolynomial_free(qp);
1548 isl_val_free(v);
1549 return isl_qpolynomial_zero_on_domain(space);
1550 }
1551
1552 qp = isl_qpolynomial_cow(qp);
1553 if (!qp)
1554 goto error;
1555
1556 qp->upoly = isl_upoly_scale_val(qp->upoly, v);
1557 if (!qp->upoly)
1558 qp = isl_qpolynomial_free(qp);
1559
1560 isl_val_free(v);
1561 return qp;
1562 error:
1563 isl_val_free(v);
1564 isl_qpolynomial_free(qp);
1565 return NULL;
1566 }
1567
1568 __isl_give isl_qpolynomial *isl_qpolynomial_mul(__isl_take isl_qpolynomial *qp1,
1569 __isl_take isl_qpolynomial *qp2)
1570 {
1571 qp1 = isl_qpolynomial_cow(qp1);
1572
1573 if (!qp1 || !qp2)
1574 goto error;
1575
1576 if (qp1->div->n_row < qp2->div->n_row)
1577 return isl_qpolynomial_mul(qp2, qp1);
1578
1579 isl_assert(qp1->dim->ctx, isl_space_is_equal(qp1->dim, qp2->dim), goto error);
1580 if (!compatible_divs(qp1->div, qp2->div))
1581 return with_merged_divs(isl_qpolynomial_mul, qp1, qp2);
1582
1583 qp1->upoly = isl_upoly_mul(qp1->upoly, isl_upoly_copy(qp2->upoly));
1584 if (!qp1->upoly)
1585 goto error;
1586
1587 isl_qpolynomial_free(qp2);
1588
1589 return qp1;
1590 error:
1591 isl_qpolynomial_free(qp1);
1592 isl_qpolynomial_free(qp2);
1593 return NULL;
1594 }
1595
1596 __isl_give isl_qpolynomial *isl_qpolynomial_pow(__isl_take isl_qpolynomial *qp,
1597 unsigned power)
1598 {
1599 qp = isl_qpolynomial_cow(qp);
1600
1601 if (!qp)
1602 return NULL;
1603
1604 qp->upoly = isl_upoly_pow(qp->upoly, power);
1605 if (!qp->upoly)
1606 goto error;
1607
1608 return qp;
1609 error:
1610 isl_qpolynomial_free(qp);
1611 return NULL;
1612 }
1613
1614 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_pow(
1615 __isl_take isl_pw_qpolynomial *pwqp, unsigned power)
1616 {
1617 int i;
1618
1619 if (power == 1)
1620 return pwqp;
1621
1622 pwqp = isl_pw_qpolynomial_cow(pwqp);
1623 if (!pwqp)
1624 return NULL;
1625
1626 for (i = 0; i < pwqp->n; ++i) {
1627 pwqp->p[i].qp = isl_qpolynomial_pow(pwqp->p[i].qp, power);
1628 if (!pwqp->p[i].qp)
1629 return isl_pw_qpolynomial_free(pwqp);
1630 }
1631
1632 return pwqp;
1633 }
1634
1635 __isl_give isl_qpolynomial *isl_qpolynomial_zero_on_domain(
1636 __isl_take isl_space *dim)
1637 {
1638 if (!dim)
1639 return NULL;
1640 return isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1641 }
1642
1643 __isl_give isl_qpolynomial *isl_qpolynomial_one_on_domain(
1644 __isl_take isl_space *dim)
1645 {
1646 if (!dim)
1647 return NULL;
1648 return isl_qpolynomial_alloc(dim, 0, isl_upoly_one(dim->ctx));
1649 }
1650
1651 __isl_give isl_qpolynomial *isl_qpolynomial_infty_on_domain(
1652 __isl_take isl_space *dim)
1653 {
1654 if (!dim)
1655 return NULL;
1656 return isl_qpolynomial_alloc(dim, 0, isl_upoly_infty(dim->ctx));
1657 }
1658
1659 __isl_give isl_qpolynomial *isl_qpolynomial_neginfty_on_domain(
1660 __isl_take isl_space *dim)
1661 {
1662 if (!dim)
1663 return NULL;
1664 return isl_qpolynomial_alloc(dim, 0, isl_upoly_neginfty(dim->ctx));
1665 }
1666
1667 __isl_give isl_qpolynomial *isl_qpolynomial_nan_on_domain(
1668 __isl_take isl_space *dim)
1669 {
1670 if (!dim)
1671 return NULL;
1672 return isl_qpolynomial_alloc(dim, 0, isl_upoly_nan(dim->ctx));
1673 }
1674
1675 __isl_give isl_qpolynomial *isl_qpolynomial_cst_on_domain(
1676 __isl_take isl_space *dim,
1677 isl_int v)
1678 {
1679 struct isl_qpolynomial *qp;
1680 struct isl_upoly_cst *cst;
1681
1682 if (!dim)
1683 return NULL;
1684
1685 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
1686 if (!qp)
1687 return NULL;
1688
1689 cst = isl_upoly_as_cst(qp->upoly);
1690 isl_int_set(cst->n, v);
1691
1692 return qp;
1693 }
1694
1695 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial *qp,
1696 isl_int *n, isl_int *d)
1697 {
1698 struct isl_upoly_cst *cst;
1699
1700 if (!qp)
1701 return -1;
1702
1703 if (!isl_upoly_is_cst(qp->upoly))
1704 return 0;
1705
1706 cst = isl_upoly_as_cst(qp->upoly);
1707 if (!cst)
1708 return -1;
1709
1710 if (n)
1711 isl_int_set(*n, cst->n);
1712 if (d)
1713 isl_int_set(*d, cst->d);
1714
1715 return 1;
1716 }
1717
1718 /* Return the constant term of "up".
1719 */
1720 static __isl_give isl_val *isl_upoly_get_constant_val(
1721 __isl_keep struct isl_upoly *up)
1722 {
1723 struct isl_upoly_cst *cst;
1724
1725 if (!up)
1726 return NULL;
1727
1728 while (!isl_upoly_is_cst(up)) {
1729 struct isl_upoly_rec *rec;
1730
1731 rec = isl_upoly_as_rec(up);
1732 if (!rec)
1733 return NULL;
1734 up = rec->p[0];
1735 }
1736
1737 cst = isl_upoly_as_cst(up);
1738 if (!cst)
1739 return NULL;
1740 return isl_val_rat_from_isl_int(cst->up.ctx, cst->n, cst->d);
1741 }
1742
1743 /* Return the constant term of "qp".
1744 */
1745 __isl_give isl_val *isl_qpolynomial_get_constant_val(
1746 __isl_keep isl_qpolynomial *qp)
1747 {
1748 if (!qp)
1749 return NULL;
1750
1751 return isl_upoly_get_constant_val(qp->upoly);
1752 }
1753
1754 int isl_upoly_is_affine(__isl_keep struct isl_upoly *up)
1755 {
1756 int is_cst;
1757 struct isl_upoly_rec *rec;
1758
1759 if (!up)
1760 return -1;
1761
1762 if (up->var < 0)
1763 return 1;
1764
1765 rec = isl_upoly_as_rec(up);
1766 if (!rec)
1767 return -1;
1768
1769 if (rec->n > 2)
1770 return 0;
1771
1772 isl_assert(up->ctx, rec->n > 1, return -1);
1773
1774 is_cst = isl_upoly_is_cst(rec->p[1]);
1775 if (is_cst < 0)
1776 return -1;
1777 if (!is_cst)
1778 return 0;
1779
1780 return isl_upoly_is_affine(rec->p[0]);
1781 }
1782
1783 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial *qp)
1784 {
1785 if (!qp)
1786 return -1;
1787
1788 if (qp->div->n_row > 0)
1789 return 0;
1790
1791 return isl_upoly_is_affine(qp->upoly);
1792 }
1793
1794 static void update_coeff(__isl_keep isl_vec *aff,
1795 __isl_keep struct isl_upoly_cst *cst, int pos)
1796 {
1797 isl_int gcd;
1798 isl_int f;
1799
1800 if (isl_int_is_zero(cst->n))
1801 return;
1802
1803 isl_int_init(gcd);
1804 isl_int_init(f);
1805 isl_int_gcd(gcd, cst->d, aff->el[0]);
1806 isl_int_divexact(f, cst->d, gcd);
1807 isl_int_divexact(gcd, aff->el[0], gcd);
1808 isl_seq_scale(aff->el, aff->el, f, aff->size);
1809 isl_int_mul(aff->el[1 + pos], gcd, cst->n);
1810 isl_int_clear(gcd);
1811 isl_int_clear(f);
1812 }
1813
1814 int isl_upoly_update_affine(__isl_keep struct isl_upoly *up,
1815 __isl_keep isl_vec *aff)
1816 {
1817 struct isl_upoly_cst *cst;
1818 struct isl_upoly_rec *rec;
1819
1820 if (!up || !aff)
1821 return -1;
1822
1823 if (up->var < 0) {
1824 struct isl_upoly_cst *cst;
1825
1826 cst = isl_upoly_as_cst(up);
1827 if (!cst)
1828 return -1;
1829 update_coeff(aff, cst, 0);
1830 return 0;
1831 }
1832
1833 rec = isl_upoly_as_rec(up);
1834 if (!rec)
1835 return -1;
1836 isl_assert(up->ctx, rec->n == 2, return -1);
1837
1838 cst = isl_upoly_as_cst(rec->p[1]);
1839 if (!cst)
1840 return -1;
1841 update_coeff(aff, cst, 1 + up->var);
1842
1843 return isl_upoly_update_affine(rec->p[0], aff);
1844 }
1845
1846 __isl_give isl_vec *isl_qpolynomial_extract_affine(
1847 __isl_keep isl_qpolynomial *qp)
1848 {
1849 isl_vec *aff;
1850 unsigned d;
1851
1852 if (!qp)
1853 return NULL;
1854
1855 d = isl_space_dim(qp->dim, isl_dim_all);
1856 aff = isl_vec_alloc(qp->div->ctx, 2 + d + qp->div->n_row);
1857 if (!aff)
1858 return NULL;
1859
1860 isl_seq_clr(aff->el + 1, 1 + d + qp->div->n_row);
1861 isl_int_set_si(aff->el[0], 1);
1862
1863 if (isl_upoly_update_affine(qp->upoly, aff) < 0)
1864 goto error;
1865
1866 return aff;
1867 error:
1868 isl_vec_free(aff);
1869 return NULL;
1870 }
1871
1872 int isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial *qp1,
1873 __isl_keep isl_qpolynomial *qp2)
1874 {
1875 int equal;
1876
1877 if (!qp1 || !qp2)
1878 return -1;
1879
1880 equal = isl_space_is_equal(qp1->dim, qp2->dim);
1881 if (equal < 0 || !equal)
1882 return equal;
1883
1884 equal = isl_mat_is_equal(qp1->div, qp2->div);
1885 if (equal < 0 || !equal)
1886 return equal;
1887
1888 return isl_upoly_is_equal(qp1->upoly, qp2->upoly);
1889 }
1890
1891 static void upoly_update_den(__isl_keep struct isl_upoly *up, isl_int *d)
1892 {
1893 int i;
1894 struct isl_upoly_rec *rec;
1895
1896 if (isl_upoly_is_cst(up)) {
1897 struct isl_upoly_cst *cst;
1898 cst = isl_upoly_as_cst(up);
1899 if (!cst)
1900 return;
1901 isl_int_lcm(*d, *d, cst->d);
1902 return;
1903 }
1904
1905 rec = isl_upoly_as_rec(up);
1906 if (!rec)
1907 return;
1908
1909 for (i = 0; i < rec->n; ++i)
1910 upoly_update_den(rec->p[i], d);
1911 }
1912
1913 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial *qp, isl_int *d)
1914 {
1915 isl_int_set_si(*d, 1);
1916 if (!qp)
1917 return;
1918 upoly_update_den(qp->upoly, d);
1919 }
1920
1921 __isl_give isl_qpolynomial *isl_qpolynomial_var_pow_on_domain(
1922 __isl_take isl_space *dim, int pos, int power)
1923 {
1924 struct isl_ctx *ctx;
1925
1926 if (!dim)
1927 return NULL;
1928
1929 ctx = dim->ctx;
1930
1931 return isl_qpolynomial_alloc(dim, 0, isl_upoly_var_pow(ctx, pos, power));
1932 }
1933
1934 __isl_give isl_qpolynomial *isl_qpolynomial_var_on_domain(__isl_take isl_space *dim,
1935 enum isl_dim_type type, unsigned pos)
1936 {
1937 if (!dim)
1938 return NULL;
1939
1940 isl_assert(dim->ctx, isl_space_dim(dim, isl_dim_in) == 0, goto error);
1941 isl_assert(dim->ctx, pos < isl_space_dim(dim, type), goto error);
1942
1943 if (type == isl_dim_set)
1944 pos += isl_space_dim(dim, isl_dim_param);
1945
1946 return isl_qpolynomial_var_pow_on_domain(dim, pos, 1);
1947 error:
1948 isl_space_free(dim);
1949 return NULL;
1950 }
1951
1952 __isl_give struct isl_upoly *isl_upoly_subs(__isl_take struct isl_upoly *up,
1953 unsigned first, unsigned n, __isl_keep struct isl_upoly **subs)
1954 {
1955 int i;
1956 struct isl_upoly_rec *rec;
1957 struct isl_upoly *base, *res;
1958
1959 if (!up)
1960 return NULL;
1961
1962 if (isl_upoly_is_cst(up))
1963 return up;
1964
1965 if (up->var < first)
1966 return up;
1967
1968 rec = isl_upoly_as_rec(up);
1969 if (!rec)
1970 goto error;
1971
1972 isl_assert(up->ctx, rec->n >= 1, goto error);
1973
1974 if (up->var >= first + n)
1975 base = isl_upoly_var_pow(up->ctx, up->var, 1);
1976 else
1977 base = isl_upoly_copy(subs[up->var - first]);
1978
1979 res = isl_upoly_subs(isl_upoly_copy(rec->p[rec->n - 1]), first, n, subs);
1980 for (i = rec->n - 2; i >= 0; --i) {
1981 struct isl_upoly *t;
1982 t = isl_upoly_subs(isl_upoly_copy(rec->p[i]), first, n, subs);
1983 res = isl_upoly_mul(res, isl_upoly_copy(base));
1984 res = isl_upoly_sum(res, t);
1985 }
1986
1987 isl_upoly_free(base);
1988 isl_upoly_free(up);
1989
1990 return res;
1991 error:
1992 isl_upoly_free(up);
1993 return NULL;
1994 }
1995
1996 __isl_give struct isl_upoly *isl_upoly_from_affine(isl_ctx *ctx, isl_int *f,
1997 isl_int denom, unsigned len)
1998 {
1999 int i;
2000 struct isl_upoly *up;
2001
2002 isl_assert(ctx, len >= 1, return NULL);
2003
2004 up = isl_upoly_rat_cst(ctx, f[0], denom);
2005 for (i = 0; i < len - 1; ++i) {
2006 struct isl_upoly *t;
2007 struct isl_upoly *c;
2008
2009 if (isl_int_is_zero(f[1 + i]))
2010 continue;
2011
2012 c = isl_upoly_rat_cst(ctx, f[1 + i], denom);
2013 t = isl_upoly_var_pow(ctx, i, 1);
2014 t = isl_upoly_mul(c, t);
2015 up = isl_upoly_sum(up, t);
2016 }
2017
2018 return up;
2019 }
2020
2021 /* Remove common factor of non-constant terms and denominator.
2022 */
2023 static void normalize_div(__isl_keep isl_qpolynomial *qp, int div)
2024 {
2025 isl_ctx *ctx = qp->div->ctx;
2026 unsigned total = qp->div->n_col - 2;
2027
2028 isl_seq_gcd(qp->div->row[div] + 2, total, &ctx->normalize_gcd);
2029 isl_int_gcd(ctx->normalize_gcd,
2030 ctx->normalize_gcd, qp->div->row[div][0]);
2031 if (isl_int_is_one(ctx->normalize_gcd))
2032 return;
2033
2034 isl_seq_scale_down(qp->div->row[div] + 2, qp->div->row[div] + 2,
2035 ctx->normalize_gcd, total);
2036 isl_int_divexact(qp->div->row[div][0], qp->div->row[div][0],
2037 ctx->normalize_gcd);
2038 isl_int_fdiv_q(qp->div->row[div][1], qp->div->row[div][1],
2039 ctx->normalize_gcd);
2040 }
2041
2042 /* Replace the integer division identified by "div" by the polynomial "s".
2043 * The integer division is assumed not to appear in the definition
2044 * of any other integer divisions.
2045 */
2046 static __isl_give isl_qpolynomial *substitute_div(
2047 __isl_take isl_qpolynomial *qp,
2048 int div, __isl_take struct isl_upoly *s)
2049 {
2050 int i;
2051 int total;
2052 int *reordering;
2053
2054 if (!qp || !s)
2055 goto error;
2056
2057 qp = isl_qpolynomial_cow(qp);
2058 if (!qp)
2059 goto error;
2060
2061 total = isl_space_dim(qp->dim, isl_dim_all);
2062 qp->upoly = isl_upoly_subs(qp->upoly, total + div, 1, &s);
2063 if (!qp->upoly)
2064 goto error;
2065
2066 reordering = isl_alloc_array(qp->dim->ctx, int, total + qp->div->n_row);
2067 if (!reordering)
2068 goto error;
2069 for (i = 0; i < total + div; ++i)
2070 reordering[i] = i;
2071 for (i = total + div + 1; i < total + qp->div->n_row; ++i)
2072 reordering[i] = i - 1;
2073 qp->div = isl_mat_drop_rows(qp->div, div, 1);
2074 qp->div = isl_mat_drop_cols(qp->div, 2 + total + div, 1);
2075 qp->upoly = reorder(qp->upoly, reordering);
2076 free(reordering);
2077
2078 if (!qp->upoly || !qp->div)
2079 goto error;
2080
2081 isl_upoly_free(s);
2082 return qp;
2083 error:
2084 isl_qpolynomial_free(qp);
2085 isl_upoly_free(s);
2086 return NULL;
2087 }
2088
2089 /* Replace all integer divisions [e/d] that turn out to not actually be integer
2090 * divisions because d is equal to 1 by their definition, i.e., e.
2091 */
2092 static __isl_give isl_qpolynomial *substitute_non_divs(
2093 __isl_take isl_qpolynomial *qp)
2094 {
2095 int i, j;
2096 int total;
2097 struct isl_upoly *s;
2098
2099 if (!qp)
2100 return NULL;
2101
2102 total = isl_space_dim(qp->dim, isl_dim_all);
2103 for (i = 0; qp && i < qp->div->n_row; ++i) {
2104 if (!isl_int_is_one(qp->div->row[i][0]))
2105 continue;
2106 for (j = i + 1; j < qp->div->n_row; ++j) {
2107 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
2108 continue;
2109 isl_seq_combine(qp->div->row[j] + 1,
2110 qp->div->ctx->one, qp->div->row[j] + 1,
2111 qp->div->row[j][2 + total + i],
2112 qp->div->row[i] + 1, 1 + total + i);
2113 isl_int_set_si(qp->div->row[j][2 + total + i], 0);
2114 normalize_div(qp, j);
2115 }
2116 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
2117 qp->div->row[i][0], qp->div->n_col - 1);
2118 qp = substitute_div(qp, i, s);
2119 --i;
2120 }
2121
2122 return qp;
2123 }
2124
2125 /* Reduce the coefficients of div "div" to lie in the interval [0, d-1],
2126 * with d the denominator. When replacing the coefficient e of x by
2127 * d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
2128 * inside the division, so we need to add floor(e/d) * x outside.
2129 * That is, we replace q by q' + floor(e/d) * x and we therefore need
2130 * to adjust the coefficient of x in each later div that depends on the
2131 * current div "div" and also in the affine expression "aff"
2132 * (if it too depends on "div").
2133 */
2134 static void reduce_div(__isl_keep isl_qpolynomial *qp, int div,
2135 __isl_keep isl_vec *aff)
2136 {
2137 int i, j;
2138 isl_int v;
2139 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2140
2141 isl_int_init(v);
2142 for (i = 0; i < 1 + total + div; ++i) {
2143 if (isl_int_is_nonneg(qp->div->row[div][1 + i]) &&
2144 isl_int_lt(qp->div->row[div][1 + i], qp->div->row[div][0]))
2145 continue;
2146 isl_int_fdiv_q(v, qp->div->row[div][1 + i], qp->div->row[div][0]);
2147 isl_int_fdiv_r(qp->div->row[div][1 + i],
2148 qp->div->row[div][1 + i], qp->div->row[div][0]);
2149 if (!isl_int_is_zero(aff->el[1 + total + div]))
2150 isl_int_addmul(aff->el[i], v, aff->el[1 + total + div]);
2151 for (j = div + 1; j < qp->div->n_row; ++j) {
2152 if (isl_int_is_zero(qp->div->row[j][2 + total + div]))
2153 continue;
2154 isl_int_addmul(qp->div->row[j][1 + i],
2155 v, qp->div->row[j][2 + total + div]);
2156 }
2157 }
2158 isl_int_clear(v);
2159 }
2160
2161 /* Check if the last non-zero coefficient is bigger that half of the
2162 * denominator. If so, we will invert the div to further reduce the number
2163 * of distinct divs that may appear.
2164 * If the last non-zero coefficient is exactly half the denominator,
2165 * then we continue looking for earlier coefficients that are bigger
2166 * than half the denominator.
2167 */
2168 static int needs_invert(__isl_keep isl_mat *div, int row)
2169 {
2170 int i;
2171 int cmp;
2172
2173 for (i = div->n_col - 1; i >= 1; --i) {
2174 if (isl_int_is_zero(div->row[row][i]))
2175 continue;
2176 isl_int_mul_ui(div->row[row][i], div->row[row][i], 2);
2177 cmp = isl_int_cmp(div->row[row][i], div->row[row][0]);
2178 isl_int_divexact_ui(div->row[row][i], div->row[row][i], 2);
2179 if (cmp)
2180 return cmp > 0;
2181 if (i == 1)
2182 return 1;
2183 }
2184
2185 return 0;
2186 }
2187
2188 /* Replace div "div" q = [e/d] by -[(-e+(d-1))/d].
2189 * We only invert the coefficients of e (and the coefficient of q in
2190 * later divs and in "aff"). After calling this function, the
2191 * coefficients of e should be reduced again.
2192 */
2193 static void invert_div(__isl_keep isl_qpolynomial *qp, int div,
2194 __isl_keep isl_vec *aff)
2195 {
2196 unsigned total = qp->div->n_col - qp->div->n_row - 2;
2197
2198 isl_seq_neg(qp->div->row[div] + 1,
2199 qp->div->row[div] + 1, qp->div->n_col - 1);
2200 isl_int_sub_ui(qp->div->row[div][1], qp->div->row[div][1], 1);
2201 isl_int_add(qp->div->row[div][1],
2202 qp->div->row[div][1], qp->div->row[div][0]);
2203 if (!isl_int_is_zero(aff->el[1 + total + div]))
2204 isl_int_neg(aff->el[1 + total + div], aff->el[1 + total + div]);
2205 isl_mat_col_mul(qp->div, 2 + total + div,
2206 qp->div->ctx->negone, 2 + total + div);
2207 }
2208
2209 /* Assuming "qp" is a monomial, reduce all its divs to have coefficients
2210 * in the interval [0, d-1], with d the denominator and such that the
2211 * last non-zero coefficient that is not equal to d/2 is smaller than d/2.
2212 *
2213 * After the reduction, some divs may have become redundant or identical,
2214 * so we call substitute_non_divs and sort_divs. If these functions
2215 * eliminate divs or merge two or more divs into one, the coefficients
2216 * of the enclosing divs may have to be reduced again, so we call
2217 * ourselves recursively if the number of divs decreases.
2218 */
2219 static __isl_give isl_qpolynomial *reduce_divs(__isl_take isl_qpolynomial *qp)
2220 {
2221 int i;
2222 isl_vec *aff = NULL;
2223 struct isl_upoly *s;
2224 unsigned n_div;
2225
2226 if (!qp)
2227 return NULL;
2228
2229 aff = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
2230 aff = isl_vec_clr(aff);
2231 if (!aff)
2232 goto error;
2233
2234 isl_int_set_si(aff->el[1 + qp->upoly->var], 1);
2235
2236 for (i = 0; i < qp->div->n_row; ++i) {
2237 normalize_div(qp, i);
2238 reduce_div(qp, i, aff);
2239 if (needs_invert(qp->div, i)) {
2240 invert_div(qp, i, aff);
2241 reduce_div(qp, i, aff);
2242 }
2243 }
2244
2245 s = isl_upoly_from_affine(qp->div->ctx, aff->el,
2246 qp->div->ctx->one, aff->size);
2247 qp->upoly = isl_upoly_subs(qp->upoly, qp->upoly->var, 1, &s);
2248 isl_upoly_free(s);
2249 if (!qp->upoly)
2250 goto error;
2251
2252 isl_vec_free(aff);
2253
2254 n_div = qp->div->n_row;
2255 qp = substitute_non_divs(qp);
2256 qp = sort_divs(qp);
2257 if (qp && qp->div->n_row < n_div)
2258 return reduce_divs(qp);
2259
2260 return qp;
2261 error:
2262 isl_qpolynomial_free(qp);
2263 isl_vec_free(aff);
2264 return NULL;
2265 }
2266
2267 __isl_give isl_qpolynomial *isl_qpolynomial_rat_cst_on_domain(
2268 __isl_take isl_space *dim, const isl_int n, const isl_int d)
2269 {
2270 struct isl_qpolynomial *qp;
2271 struct isl_upoly_cst *cst;
2272
2273 if (!dim)
2274 return NULL;
2275
2276 qp = isl_qpolynomial_alloc(dim, 0, isl_upoly_zero(dim->ctx));
2277 if (!qp)
2278 return NULL;
2279
2280 cst = isl_upoly_as_cst(qp->upoly);
2281 isl_int_set(cst->n, n);
2282 isl_int_set(cst->d, d);
2283
2284 return qp;
2285 }
2286
2287 /* Return an isl_qpolynomial that is equal to "val" on domain space "domain".
2288 */
2289 __isl_give isl_qpolynomial *isl_qpolynomial_val_on_domain(
2290 __isl_take isl_space *domain, __isl_take isl_val *val)
2291 {
2292 isl_qpolynomial *qp;
2293 struct isl_upoly_cst *cst;
2294
2295 if (!domain || !val)
2296 goto error;
2297
2298 qp = isl_qpolynomial_alloc(isl_space_copy(domain), 0,
2299 isl_upoly_zero(domain->ctx));
2300 if (!qp)
2301 goto error;
2302
2303 cst = isl_upoly_as_cst(qp->upoly);
2304 isl_int_set(cst->n, val->n);
2305 isl_int_set(cst->d, val->d);
2306
2307 isl_space_free(domain);
2308 isl_val_free(val);
2309 return qp;
2310 error:
2311 isl_space_free(domain);
2312 isl_val_free(val);
2313 return NULL;
2314 }
2315
2316 static int up_set_active(__isl_keep struct isl_upoly *up, int *active, int d)
2317 {
2318 struct isl_upoly_rec *rec;
2319 int i;
2320
2321 if (!up)
2322 return -1;
2323
2324 if (isl_upoly_is_cst(up))
2325 return 0;
2326
2327 if (up->var < d)
2328 active[up->var] = 1;
2329
2330 rec = isl_upoly_as_rec(up);
2331 for (i = 0; i < rec->n; ++i)
2332 if (up_set_active(rec->p[i], active, d) < 0)
2333 return -1;
2334
2335 return 0;
2336 }
2337
2338 static int set_active(__isl_keep isl_qpolynomial *qp, int *active)
2339 {
2340 int i, j;
2341 int d = isl_space_dim(qp->dim, isl_dim_all);
2342
2343 if (!qp || !active)
2344 return -1;
2345
2346 for (i = 0; i < d; ++i)
2347 for (j = 0; j < qp->div->n_row; ++j) {
2348 if (isl_int_is_zero(qp->div->row[j][2 + i]))
2349 continue;
2350 active[i] = 1;
2351 break;
2352 }
2353
2354 return up_set_active(qp->upoly, active, d);
2355 }
2356
2357 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial *qp,
2358 enum isl_dim_type type, unsigned first, unsigned n)
2359 {
2360 int i;
2361 int *active = NULL;
2362 int involves = 0;
2363
2364 if (!qp)
2365 return -1;
2366 if (n == 0)
2367 return 0;
2368
2369 isl_assert(qp->dim->ctx,
2370 first + n <= isl_qpolynomial_dim(qp, type), return -1);
2371 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2372 type == isl_dim_in, return -1);
2373
2374 active = isl_calloc_array(qp->dim->ctx, int,
2375 isl_space_dim(qp->dim, isl_dim_all));
2376 if (set_active(qp, active) < 0)
2377 goto error;
2378
2379 if (type == isl_dim_in)
2380 first += isl_space_dim(qp->dim, isl_dim_param);
2381 for (i = 0; i < n; ++i)
2382 if (active[first + i]) {
2383 involves = 1;
2384 break;
2385 }
2386
2387 free(active);
2388
2389 return involves;
2390 error:
2391 free(active);
2392 return -1;
2393 }
2394
2395 /* Remove divs that do not appear in the quasi-polynomial, nor in any
2396 * of the divs that do appear in the quasi-polynomial.
2397 */
2398 static __isl_give isl_qpolynomial *remove_redundant_divs(
2399 __isl_take isl_qpolynomial *qp)
2400 {
2401 int i, j;
2402 int d;
2403 int len;
2404 int skip;
2405 int *active = NULL;
2406 int *reordering = NULL;
2407 int redundant = 0;
2408 int n_div;
2409 isl_ctx *ctx;
2410
2411 if (!qp)
2412 return NULL;
2413 if (qp->div->n_row == 0)
2414 return qp;
2415
2416 d = isl_space_dim(qp->dim, isl_dim_all);
2417 len = qp->div->n_col - 2;
2418 ctx = isl_qpolynomial_get_ctx(qp);
2419 active = isl_calloc_array(ctx, int, len);
2420 if (!active)
2421 goto error;
2422
2423 if (up_set_active(qp->upoly, active, len) < 0)
2424 goto error;
2425
2426 for (i = qp->div->n_row - 1; i >= 0; --i) {
2427 if (!active[d + i]) {
2428 redundant = 1;
2429 continue;
2430 }
2431 for (j = 0; j < i; ++j) {
2432 if (isl_int_is_zero(qp->div->row[i][2 + d + j]))
2433 continue;
2434 active[d + j] = 1;
2435 break;
2436 }
2437 }
2438
2439 if (!redundant) {
2440 free(active);
2441 return qp;
2442 }
2443
2444 reordering = isl_alloc_array(qp->div->ctx, int, len);
2445 if (!reordering)
2446 goto error;
2447
2448 for (i = 0; i < d; ++i)
2449 reordering[i] = i;
2450
2451 skip = 0;
2452 n_div = qp->div->n_row;
2453 for (i = 0; i < n_div; ++i) {
2454 if (!active[d + i]) {
2455 qp->div = isl_mat_drop_rows(qp->div, i - skip, 1);
2456 qp->div = isl_mat_drop_cols(qp->div,
2457 2 + d + i - skip, 1);
2458 skip++;
2459 }
2460 reordering[d + i] = d + i - skip;
2461 }
2462
2463 qp->upoly = reorder(qp->upoly, reordering);
2464
2465 if (!qp->upoly || !qp->div)
2466 goto error;
2467
2468 free(active);
2469 free(reordering);
2470
2471 return qp;
2472 error:
2473 free(active);
2474 free(reordering);
2475 isl_qpolynomial_free(qp);
2476 return NULL;
2477 }
2478
2479 __isl_give struct isl_upoly *isl_upoly_drop(__isl_take struct isl_upoly *up,
2480 unsigned first, unsigned n)
2481 {
2482 int i;
2483 struct isl_upoly_rec *rec;
2484
2485 if (!up)
2486 return NULL;
2487 if (n == 0 || up->var < 0 || up->var < first)
2488 return up;
2489 if (up->var < first + n) {
2490 up = replace_by_constant_term(up);
2491 return isl_upoly_drop(up, first, n);
2492 }
2493 up = isl_upoly_cow(up);
2494 if (!up)
2495 return NULL;
2496 up->var -= n;
2497 rec = isl_upoly_as_rec(up);
2498 if (!rec)
2499 goto error;
2500
2501 for (i = 0; i < rec->n; ++i) {
2502 rec->p[i] = isl_upoly_drop(rec->p[i], first, n);
2503 if (!rec->p[i])
2504 goto error;
2505 }
2506
2507 return up;
2508 error:
2509 isl_upoly_free(up);
2510 return NULL;
2511 }
2512
2513 __isl_give isl_qpolynomial *isl_qpolynomial_set_dim_name(
2514 __isl_take isl_qpolynomial *qp,
2515 enum isl_dim_type type, unsigned pos, const char *s)
2516 {
2517 qp = isl_qpolynomial_cow(qp);
2518 if (!qp)
2519 return NULL;
2520 qp->dim = isl_space_set_dim_name(qp->dim, type, pos, s);
2521 if (!qp->dim)
2522 goto error;
2523 return qp;
2524 error:
2525 isl_qpolynomial_free(qp);
2526 return NULL;
2527 }
2528
2529 __isl_give isl_qpolynomial *isl_qpolynomial_drop_dims(
2530 __isl_take isl_qpolynomial *qp,
2531 enum isl_dim_type type, unsigned first, unsigned n)
2532 {
2533 if (!qp)
2534 return NULL;
2535 if (type == isl_dim_out)
2536 isl_die(qp->dim->ctx, isl_error_invalid,
2537 "cannot drop output/set dimension",
2538 goto error);
2539 if (type == isl_dim_in)
2540 type = isl_dim_set;
2541 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2542 return qp;
2543
2544 qp = isl_qpolynomial_cow(qp);
2545 if (!qp)
2546 return NULL;
2547
2548 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
2549 goto error);
2550 isl_assert(qp->dim->ctx, type == isl_dim_param ||
2551 type == isl_dim_set, goto error);
2552
2553 qp->dim = isl_space_drop_dims(qp->dim, type, first, n);
2554 if (!qp->dim)
2555 goto error;
2556
2557 if (type == isl_dim_set)
2558 first += isl_space_dim(qp->dim, isl_dim_param);
2559
2560 qp->div = isl_mat_drop_cols(qp->div, 2 + first, n);
2561 if (!qp->div)
2562 goto error;
2563
2564 qp->upoly = isl_upoly_drop(qp->upoly, first, n);
2565 if (!qp->upoly)
2566 goto error;
2567
2568 return qp;
2569 error:
2570 isl_qpolynomial_free(qp);
2571 return NULL;
2572 }
2573
2574 /* Project the domain of the quasi-polynomial onto its parameter space.
2575 * The quasi-polynomial may not involve any of the domain dimensions.
2576 */
2577 __isl_give isl_qpolynomial *isl_qpolynomial_project_domain_on_params(
2578 __isl_take isl_qpolynomial *qp)
2579 {
2580 isl_space *space;
2581 unsigned n;
2582 int involves;
2583
2584 n = isl_qpolynomial_dim(qp, isl_dim_in);
2585 involves = isl_qpolynomial_involves_dims(qp, isl_dim_in, 0, n);
2586 if (involves < 0)
2587 return isl_qpolynomial_free(qp);
2588 if (involves)
2589 isl_die(isl_qpolynomial_get_ctx(qp), isl_error_invalid,
2590 "polynomial involves some of the domain dimensions",
2591 return isl_qpolynomial_free(qp));
2592 qp = isl_qpolynomial_drop_dims(qp, isl_dim_in, 0, n);
2593 space = isl_qpolynomial_get_domain_space(qp);
2594 space = isl_space_params(space);
2595 qp = isl_qpolynomial_reset_domain_space(qp, space);
2596 return qp;
2597 }
2598
2599 static __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities_lifted(
2600 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2601 {
2602 int i, j, k;
2603 isl_int denom;
2604 unsigned total;
2605 unsigned n_div;
2606 struct isl_upoly *up;
2607
2608 if (!eq)
2609 goto error;
2610 if (eq->n_eq == 0) {
2611 isl_basic_set_free(eq);
2612 return qp;
2613 }
2614
2615 qp = isl_qpolynomial_cow(qp);
2616 if (!qp)
2617 goto error;
2618 qp->div = isl_mat_cow(qp->div);
2619 if (!qp->div)
2620 goto error;
2621
2622 total = 1 + isl_space_dim(eq->dim, isl_dim_all);
2623 n_div = eq->n_div;
2624 isl_int_init(denom);
2625 for (i = 0; i < eq->n_eq; ++i) {
2626 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
2627 if (j < 0 || j == 0 || j >= total)
2628 continue;
2629
2630 for (k = 0; k < qp->div->n_row; ++k) {
2631 if (isl_int_is_zero(qp->div->row[k][1 + j]))
2632 continue;
2633 isl_seq_elim(qp->div->row[k] + 1, eq->eq[i], j, total,
2634 &qp->div->row[k][0]);
2635 normalize_div(qp, k);
2636 }
2637
2638 if (isl_int_is_pos(eq->eq[i][j]))
2639 isl_seq_neg(eq->eq[i], eq->eq[i], total);
2640 isl_int_abs(denom, eq->eq[i][j]);
2641 isl_int_set_si(eq->eq[i][j], 0);
2642
2643 up = isl_upoly_from_affine(qp->dim->ctx,
2644 eq->eq[i], denom, total);
2645 qp->upoly = isl_upoly_subs(qp->upoly, j - 1, 1, &up);
2646 isl_upoly_free(up);
2647 }
2648 isl_int_clear(denom);
2649
2650 if (!qp->upoly)
2651 goto error;
2652
2653 isl_basic_set_free(eq);
2654
2655 qp = substitute_non_divs(qp);
2656 qp = sort_divs(qp);
2657
2658 return qp;
2659 error:
2660 isl_basic_set_free(eq);
2661 isl_qpolynomial_free(qp);
2662 return NULL;
2663 }
2664
2665 /* Exploit the equalities in "eq" to simplify the quasi-polynomial.
2666 */
2667 __isl_give isl_qpolynomial *isl_qpolynomial_substitute_equalities(
2668 __isl_take isl_qpolynomial *qp, __isl_take isl_basic_set *eq)
2669 {
2670 if (!qp || !eq)
2671 goto error;
2672 if (qp->div->n_row > 0)
2673 eq = isl_basic_set_add_dims(eq, isl_dim_set, qp->div->n_row);
2674 return isl_qpolynomial_substitute_equalities_lifted(qp, eq);
2675 error:
2676 isl_basic_set_free(eq);
2677 isl_qpolynomial_free(qp);
2678 return NULL;
2679 }
2680
2681 static __isl_give isl_basic_set *add_div_constraints(
2682 __isl_take isl_basic_set *bset, __isl_take isl_mat *div)
2683 {
2684 int i;
2685 unsigned total;
2686
2687 if (!bset || !div)
2688 goto error;
2689
2690 bset = isl_basic_set_extend_constraints(bset, 0, 2 * div->n_row);
2691 if (!bset)
2692 goto error;
2693 total = isl_basic_set_total_dim(bset);
2694 for (i = 0; i < div->n_row; ++i)
2695 if (isl_basic_set_add_div_constraints_var(bset,
2696 total - div->n_row + i, div->row[i]) < 0)
2697 goto error;
2698
2699 isl_mat_free(div);
2700 return bset;
2701 error:
2702 isl_mat_free(div);
2703 isl_basic_set_free(bset);
2704 return NULL;
2705 }
2706
2707 /* Look for equalities among the variables shared by context and qp
2708 * and the integer divisions of qp, if any.
2709 * The equalities are then used to eliminate variables and/or integer
2710 * divisions from qp.
2711 */
2712 __isl_give isl_qpolynomial *isl_qpolynomial_gist(
2713 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2714 {
2715 isl_basic_set *aff;
2716
2717 if (!qp)
2718 goto error;
2719 if (qp->div->n_row > 0) {
2720 isl_basic_set *bset;
2721 context = isl_set_add_dims(context, isl_dim_set,
2722 qp->div->n_row);
2723 bset = isl_basic_set_universe(isl_set_get_space(context));
2724 bset = add_div_constraints(bset, isl_mat_copy(qp->div));
2725 context = isl_set_intersect(context,
2726 isl_set_from_basic_set(bset));
2727 }
2728
2729 aff = isl_set_affine_hull(context);
2730 return isl_qpolynomial_substitute_equalities_lifted(qp, aff);
2731 error:
2732 isl_qpolynomial_free(qp);
2733 isl_set_free(context);
2734 return NULL;
2735 }
2736
2737 __isl_give isl_qpolynomial *isl_qpolynomial_gist_params(
2738 __isl_take isl_qpolynomial *qp, __isl_take isl_set *context)
2739 {
2740 isl_space *space = isl_qpolynomial_get_domain_space(qp);
2741 isl_set *dom_context = isl_set_universe(space);
2742 dom_context = isl_set_intersect_params(dom_context, context);
2743 return isl_qpolynomial_gist(qp, dom_context);
2744 }
2745
2746 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_qpolynomial(
2747 __isl_take isl_qpolynomial *qp)
2748 {
2749 isl_set *dom;
2750
2751 if (!qp)
2752 return NULL;
2753 if (isl_qpolynomial_is_zero(qp)) {
2754 isl_space *dim = isl_qpolynomial_get_space(qp);
2755 isl_qpolynomial_free(qp);
2756 return isl_pw_qpolynomial_zero(dim);
2757 }
2758
2759 dom = isl_set_universe(isl_qpolynomial_get_domain_space(qp));
2760 return isl_pw_qpolynomial_alloc(dom, qp);
2761 }
2762
2763 #undef PW
2764 #define PW isl_pw_qpolynomial
2765 #undef EL
2766 #define EL isl_qpolynomial
2767 #undef EL_IS_ZERO
2768 #define EL_IS_ZERO is_zero
2769 #undef ZERO
2770 #define ZERO zero
2771 #undef IS_ZERO
2772 #define IS_ZERO is_zero
2773 #undef FIELD
2774 #define FIELD qp
2775 #undef DEFAULT_IS_ZERO
2776 #define DEFAULT_IS_ZERO 1
2777
2778 #define NO_PULLBACK
2779
2780 #include <isl_pw_templ.c>
2781
2782 #undef UNION
2783 #define UNION isl_union_pw_qpolynomial
2784 #undef PART
2785 #define PART isl_pw_qpolynomial
2786 #undef PARTS
2787 #define PARTS pw_qpolynomial
2788 #define ALIGN_DOMAIN
2789
2790 #include <isl_union_templ.c>
2791
2792 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
2793 {
2794 if (!pwqp)
2795 return -1;
2796
2797 if (pwqp->n != -1)
2798 return 0;
2799
2800 if (!isl_set_plain_is_universe(pwqp->p[0].set))
2801 return 0;
2802
2803 return isl_qpolynomial_is_one(pwqp->p[0].qp);
2804 }
2805
2806 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
2807 __isl_take isl_pw_qpolynomial *pwqp1,
2808 __isl_take isl_pw_qpolynomial *pwqp2)
2809 {
2810 return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2);
2811 }
2812
2813 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
2814 __isl_take isl_pw_qpolynomial *pwqp1,
2815 __isl_take isl_pw_qpolynomial *pwqp2)
2816 {
2817 int i, j, n;
2818 struct isl_pw_qpolynomial *res;
2819
2820 if (!pwqp1 || !pwqp2)
2821 goto error;
2822
2823 isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim),
2824 goto error);
2825
2826 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
2827 isl_pw_qpolynomial_free(pwqp2);
2828 return pwqp1;
2829 }
2830
2831 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
2832 isl_pw_qpolynomial_free(pwqp1);
2833 return pwqp2;
2834 }
2835
2836 if (isl_pw_qpolynomial_is_one(pwqp1)) {
2837 isl_pw_qpolynomial_free(pwqp1);
2838 return pwqp2;
2839 }
2840
2841 if (isl_pw_qpolynomial_is_one(pwqp2)) {
2842 isl_pw_qpolynomial_free(pwqp2);
2843 return pwqp1;
2844 }
2845
2846 n = pwqp1->n * pwqp2->n;
2847 res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n);
2848
2849 for (i = 0; i < pwqp1->n; ++i) {
2850 for (j = 0; j < pwqp2->n; ++j) {
2851 struct isl_set *common;
2852 struct isl_qpolynomial *prod;
2853 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
2854 isl_set_copy(pwqp2->p[j].set));
2855 if (isl_set_plain_is_empty(common)) {
2856 isl_set_free(common);
2857 continue;
2858 }
2859
2860 prod = isl_qpolynomial_mul(
2861 isl_qpolynomial_copy(pwqp1->p[i].qp),
2862 isl_qpolynomial_copy(pwqp2->p[j].qp));
2863
2864 res = isl_pw_qpolynomial_add_piece(res, common, prod);
2865 }
2866 }
2867
2868 isl_pw_qpolynomial_free(pwqp1);
2869 isl_pw_qpolynomial_free(pwqp2);
2870
2871 return res;
2872 error:
2873 isl_pw_qpolynomial_free(pwqp1);
2874 isl_pw_qpolynomial_free(pwqp2);
2875 return NULL;
2876 }
2877
2878 __isl_give isl_val *isl_upoly_eval(__isl_take struct isl_upoly *up,
2879 __isl_take isl_vec *vec)
2880 {
2881 int i;
2882 struct isl_upoly_rec *rec;
2883 isl_val *res;
2884 isl_val *base;
2885
2886 if (isl_upoly_is_cst(up)) {
2887 isl_vec_free(vec);
2888 res = isl_upoly_get_constant_val(up);
2889 isl_upoly_free(up);
2890 return res;
2891 }
2892
2893 rec = isl_upoly_as_rec(up);
2894 if (!rec)
2895 goto error;
2896
2897 isl_assert(up->ctx, rec->n >= 1, goto error);
2898
2899 base = isl_val_rat_from_isl_int(up->ctx,
2900 vec->el[1 + up->var], vec->el[0]);
2901
2902 res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
2903 isl_vec_copy(vec));
2904
2905 for (i = rec->n - 2; i >= 0; --i) {
2906 res = isl_val_mul(res, isl_val_copy(base));
2907 res = isl_val_add(res,
2908 isl_upoly_eval(isl_upoly_copy(rec->p[i]),
2909 isl_vec_copy(vec)));
2910 }
2911
2912 isl_val_free(base);
2913 isl_upoly_free(up);
2914 isl_vec_free(vec);
2915 return res;
2916 error:
2917 isl_upoly_free(up);
2918 isl_vec_free(vec);
2919 return NULL;
2920 }
2921
2922 __isl_give isl_val *isl_qpolynomial_eval(__isl_take isl_qpolynomial *qp,
2923 __isl_take isl_point *pnt)
2924 {
2925 isl_vec *ext;
2926 isl_val *v;
2927
2928 if (!qp || !pnt)
2929 goto error;
2930 isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error);
2931
2932 if (qp->div->n_row == 0)
2933 ext = isl_vec_copy(pnt->vec);
2934 else {
2935 int i;
2936 unsigned dim = isl_space_dim(qp->dim, isl_dim_all);
2937 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
2938 if (!ext)
2939 goto error;
2940
2941 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
2942 for (i = 0; i < qp->div->n_row; ++i) {
2943 isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
2944 1 + dim + i, &ext->el[1+dim+i]);
2945 isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
2946 qp->div->row[i][0]);
2947 }
2948 }
2949
2950 v = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
2951
2952 isl_qpolynomial_free(qp);
2953 isl_point_free(pnt);
2954
2955 return v;
2956 error:
2957 isl_qpolynomial_free(qp);
2958 isl_point_free(pnt);
2959 return NULL;
2960 }
2961
2962 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
2963 __isl_keep struct isl_upoly_cst *cst2)
2964 {
2965 int cmp;
2966 isl_int t;
2967 isl_int_init(t);
2968 isl_int_mul(t, cst1->n, cst2->d);
2969 isl_int_submul(t, cst2->n, cst1->d);
2970 cmp = isl_int_sgn(t);
2971 isl_int_clear(t);
2972 return cmp;
2973 }
2974
2975 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
2976 __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
2977 unsigned first, unsigned n)
2978 {
2979 unsigned total;
2980 unsigned g_pos;
2981 int *exp;
2982
2983 if (!qp)
2984 return NULL;
2985 if (type == isl_dim_out)
2986 isl_die(qp->div->ctx, isl_error_invalid,
2987 "cannot insert output/set dimensions",
2988 goto error);
2989 if (type == isl_dim_in)
2990 type = isl_dim_set;
2991 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
2992 return qp;
2993
2994 qp = isl_qpolynomial_cow(qp);
2995 if (!qp)
2996 return NULL;
2997
2998 isl_assert(qp->div->ctx, first <= isl_space_dim(qp->dim, type),
2999 goto error);
3000
3001 g_pos = pos(qp->dim, type) + first;
3002
3003 qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
3004 if (!qp->div)
3005 goto error;
3006
3007 total = qp->div->n_col - 2;
3008 if (total > g_pos) {
3009 int i;
3010 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
3011 if (!exp)
3012 goto error;
3013 for (i = 0; i < total - g_pos; ++i)
3014 exp[i] = i + n;
3015 qp->upoly = expand(qp->upoly, exp, g_pos);
3016 free(exp);
3017 if (!qp->upoly)
3018 goto error;
3019 }
3020
3021 qp->dim = isl_space_insert_dims(qp->dim, type, first, n);
3022 if (!qp->dim)
3023 goto error;
3024
3025 return qp;
3026 error:
3027 isl_qpolynomial_free(qp);
3028 return NULL;
3029 }
3030
3031 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
3032 __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
3033 {
3034 unsigned pos;
3035
3036 pos = isl_qpolynomial_dim(qp, type);
3037
3038 return isl_qpolynomial_insert_dims(qp, type, pos, n);
3039 }
3040
3041 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
3042 __isl_take isl_pw_qpolynomial *pwqp,
3043 enum isl_dim_type type, unsigned n)
3044 {
3045 unsigned pos;
3046
3047 pos = isl_pw_qpolynomial_dim(pwqp, type);
3048
3049 return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
3050 }
3051
3052 static int *reordering_move(isl_ctx *ctx,
3053 unsigned len, unsigned dst, unsigned src, unsigned n)
3054 {
3055 int i;
3056 int *reordering;
3057
3058 reordering = isl_alloc_array(ctx, int, len);
3059 if (!reordering)
3060 return NULL;
3061
3062 if (dst <= src) {
3063 for (i = 0; i < dst; ++i)
3064 reordering[i] = i;
3065 for (i = 0; i < n; ++i)
3066 reordering[src + i] = dst + i;
3067 for (i = 0; i < src - dst; ++i)
3068 reordering[dst + i] = dst + n + i;
3069 for (i = 0; i < len - src - n; ++i)
3070 reordering[src + n + i] = src + n + i;
3071 } else {
3072 for (i = 0; i < src; ++i)
3073 reordering[i] = i;
3074 for (i = 0; i < n; ++i)
3075 reordering[src + i] = dst + i;
3076 for (i = 0; i < dst - src; ++i)
3077 reordering[src + n + i] = src + i;
3078 for (i = 0; i < len - dst - n; ++i)
3079 reordering[dst + n + i] = dst + n + i;
3080 }
3081
3082 return reordering;
3083 }
3084
3085 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
3086 __isl_take isl_qpolynomial *qp,
3087 enum isl_dim_type dst_type, unsigned dst_pos,
3088 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
3089 {
3090 unsigned g_dst_pos;
3091 unsigned g_src_pos;
3092 int *reordering;
3093
3094 if (n == 0)
3095 return qp;
3096
3097 qp = isl_qpolynomial_cow(qp);
3098 if (!qp)
3099 return NULL;
3100
3101 if (dst_type == isl_dim_out || src_type == isl_dim_out)
3102 isl_die(qp->dim->ctx, isl_error_invalid,
3103 "cannot move output/set dimension",
3104 goto error);
3105 if (dst_type == isl_dim_in)
3106 dst_type = isl_dim_set;
3107 if (src_type == isl_dim_in)
3108 src_type = isl_dim_set;
3109
3110 isl_assert(qp->dim->ctx, src_pos + n <= isl_space_dim(qp->dim, src_type),
3111 goto error);
3112
3113 g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
3114 g_src_pos = pos(qp->dim, src_type) + src_pos;
3115 if (dst_type > src_type)
3116 g_dst_pos -= n;
3117
3118 qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
3119 if (!qp->div)
3120 goto error;
3121 qp = sort_divs(qp);
3122 if (!qp)
3123 goto error;
3124
3125 reordering = reordering_move(qp->dim->ctx,
3126 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
3127 if (!reordering)
3128 goto error;
3129
3130 qp->upoly = reorder(qp->upoly, reordering);
3131 free(reordering);
3132 if (!qp->upoly)
3133 goto error;
3134
3135 qp->dim = isl_space_move_dims(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
3136 if (!qp->dim)
3137 goto error;
3138
3139 return qp;
3140 error:
3141 isl_qpolynomial_free(qp);
3142 return NULL;
3143 }
3144
3145 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_space *dim,
3146 isl_int *f, isl_int denom)
3147 {
3148 struct isl_upoly *up;
3149
3150 dim = isl_space_domain(dim);
3151 if (!dim)
3152 return NULL;
3153
3154 up = isl_upoly_from_affine(dim->ctx, f, denom,
3155 1 + isl_space_dim(dim, isl_dim_all));
3156
3157 return isl_qpolynomial_alloc(dim, 0, up);
3158 }
3159
3160 __isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
3161 {
3162 isl_ctx *ctx;
3163 struct isl_upoly *up;
3164 isl_qpolynomial *qp;
3165
3166 if (!aff)
3167 return NULL;
3168
3169 ctx = isl_aff_get_ctx(aff);
3170 up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
3171 aff->v->size - 1);
3172
3173 qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff),
3174 aff->ls->div->n_row, up);
3175 if (!qp)
3176 goto error;
3177
3178 isl_mat_free(qp->div);
3179 qp->div = isl_mat_copy(aff->ls->div);
3180 qp->div = isl_mat_cow(qp->div);
3181 if (!qp->div)
3182 goto error;
3183
3184 isl_aff_free(aff);
3185 qp = reduce_divs(qp);
3186 qp = remove_redundant_divs(qp);
3187 return qp;
3188 error:
3189 isl_aff_free(aff);
3190 return NULL;
3191 }
3192
3193 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff(
3194 __isl_take isl_pw_aff *pwaff)
3195 {
3196 int i;
3197 isl_pw_qpolynomial *pwqp;
3198
3199 if (!pwaff)
3200 return NULL;
3201
3202 pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff),
3203 pwaff->n);
3204
3205 for (i = 0; i < pwaff->n; ++i) {
3206 isl_set *dom;
3207 isl_qpolynomial *qp;
3208
3209 dom = isl_set_copy(pwaff->p[i].set);
3210 qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff));
3211 pwqp = isl_pw_qpolynomial_add_piece(pwqp, dom, qp);
3212 }
3213
3214 isl_pw_aff_free(pwaff);
3215 return pwqp;
3216 }
3217
3218 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
3219 __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
3220 {
3221 isl_aff *aff;
3222
3223 aff = isl_constraint_get_bound(c, type, pos);
3224 isl_constraint_free(c);
3225 return isl_qpolynomial_from_aff(aff);
3226 }
3227
3228 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
3229 * in "qp" by subs[i].
3230 */
3231 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
3232 __isl_take isl_qpolynomial *qp,
3233 enum isl_dim_type type, unsigned first, unsigned n,
3234 __isl_keep isl_qpolynomial **subs)
3235 {
3236 int i;
3237 struct isl_upoly **ups;
3238
3239 if (n == 0)
3240 return qp;
3241
3242 qp = isl_qpolynomial_cow(qp);
3243 if (!qp)
3244 return NULL;
3245
3246 if (type == isl_dim_out)
3247 isl_die(qp->dim->ctx, isl_error_invalid,
3248 "cannot substitute output/set dimension",
3249 goto error);
3250 if (type == isl_dim_in)
3251 type = isl_dim_set;
3252
3253 for (i = 0; i < n; ++i)
3254 if (!subs[i])
3255 goto error;
3256
3257 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
3258 goto error);
3259
3260 for (i = 0; i < n; ++i)
3261 isl_assert(qp->dim->ctx, isl_space_is_equal(qp->dim, subs[i]->dim),
3262 goto error);
3263
3264 isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
3265 for (i = 0; i < n; ++i)
3266 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
3267
3268 first += pos(qp->dim, type);
3269
3270 ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
3271 if (!ups)
3272 goto error;
3273 for (i = 0; i < n; ++i)
3274 ups[i] = subs[i]->upoly;
3275
3276 qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
3277
3278 free(ups);
3279
3280 if (!qp->upoly)
3281 goto error;
3282
3283 return qp;
3284 error:
3285 isl_qpolynomial_free(qp);
3286 return NULL;
3287 }
3288
3289 /* Extend "bset" with extra set dimensions for each integer division
3290 * in "qp" and then call "fn" with the extended bset and the polynomial
3291 * that results from replacing each of the integer divisions by the
3292 * corresponding extra set dimension.
3293 */
3294 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
3295 __isl_keep isl_basic_set *bset,
3296 int (*fn)(__isl_take isl_basic_set *bset,
3297 __isl_take isl_qpolynomial *poly, void *user), void *user)
3298 {
3299 isl_space *dim;
3300 isl_mat *div;
3301 isl_qpolynomial *poly;
3302
3303 if (!qp || !bset)
3304 goto error;
3305 if (qp->div->n_row == 0)
3306 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
3307 user);
3308
3309 div = isl_mat_copy(qp->div);
3310 dim = isl_space_copy(qp->dim);
3311 dim = isl_space_add_dims(dim, isl_dim_set, qp->div->n_row);
3312 poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
3313 bset = isl_basic_set_copy(bset);
3314 bset = isl_basic_set_add_dims(bset, isl_dim_set, qp->div->n_row);
3315 bset = add_div_constraints(bset, div);
3316
3317 return fn(bset, poly, user);
3318 error:
3319 return -1;
3320 }
3321
3322 /* Return total degree in variables first (inclusive) up to last (exclusive).
3323 */
3324 int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
3325 {
3326 int deg = -1;
3327 int i;
3328 struct isl_upoly_rec *rec;
3329
3330 if (!up)
3331 return -2;
3332 if (isl_upoly_is_zero(up))
3333 return -1;
3334 if (isl_upoly_is_cst(up) || up->var < first)
3335 return 0;
3336
3337 rec = isl_upoly_as_rec(up);
3338 if (!rec)
3339 return -2;
3340
3341 for (i = 0; i < rec->n; ++i) {
3342 int d;
3343
3344 if (isl_upoly_is_zero(rec->p[i]))
3345 continue;
3346 d = isl_upoly_degree(rec->p[i], first, last);
3347 if (up->var < last)
3348 d += i;
3349 if (d > deg)
3350 deg = d;
3351 }
3352
3353 return deg;
3354 }
3355
3356 /* Return total degree in set variables.
3357 */
3358 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
3359 {
3360 unsigned ovar;
3361 unsigned nvar;
3362
3363 if (!poly)
3364 return -2;
3365
3366 ovar = isl_space_offset(poly->dim, isl_dim_set);
3367 nvar = isl_space_dim(poly->dim, isl_dim_set);
3368 return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
3369 }
3370
3371 __isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
3372 unsigned pos, int deg)
3373 {
3374 int i;
3375 struct isl_upoly_rec *rec;
3376
3377 if (!up)
3378 return NULL;
3379
3380 if (isl_upoly_is_cst(up) || up->var < pos) {
3381 if (deg == 0)
3382 return isl_upoly_copy(up);
3383 else
3384 return isl_upoly_zero(up->ctx);
3385 }
3386
3387 rec = isl_upoly_as_rec(up);
3388 if (!rec)
3389 return NULL;
3390
3391 if (up->var == pos) {
3392 if (deg < rec->n)
3393 return isl_upoly_copy(rec->p[deg]);
3394 else
3395 return isl_upoly_zero(up->ctx);
3396 }
3397
3398 up = isl_upoly_copy(up);
3399 up = isl_upoly_cow(up);
3400 rec = isl_upoly_as_rec(up);
3401 if (!rec)
3402 goto error;
3403
3404 for (i = 0; i < rec->n; ++i) {
3405 struct isl_upoly *t;
3406 t = isl_upoly_coeff(rec->p[i], pos, deg);
3407 if (!t)
3408 goto error;
3409 isl_upoly_free(rec->p[i]);
3410 rec->p[i] = t;
3411 }
3412
3413 return up;
3414 error:
3415 isl_upoly_free(up);
3416 return NULL;
3417 }
3418
3419 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
3420 */
3421 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
3422 __isl_keep isl_qpolynomial *qp,
3423 enum isl_dim_type type, unsigned t_pos, int deg)
3424 {
3425 unsigned g_pos;
3426 struct isl_upoly *up;
3427 isl_qpolynomial *c;
3428
3429 if (!qp)
3430 return NULL;
3431
3432 if (type == isl_dim_out)
3433 isl_die(qp->div->ctx, isl_error_invalid,
3434 "output/set dimension does not have a coefficient",
3435 return NULL);
3436 if (type == isl_dim_in)
3437 type = isl_dim_set;
3438
3439 isl_assert(qp->div->ctx, t_pos < isl_space_dim(qp->dim, type),
3440 return NULL);
3441
3442 g_pos = pos(qp->dim, type) + t_pos;
3443 up = isl_upoly_coeff(qp->upoly, g_pos, deg);
3444
3445 c = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row, up);
3446 if (!c)
3447 return NULL;
3448 isl_mat_free(c->div);
3449 c->div = isl_mat_copy(qp->div);
3450 if (!c->div)
3451 goto error;
3452 return c;
3453 error:
3454 isl_qpolynomial_free(c);
3455 return NULL;
3456 }
3457
3458 /* Homogenize the polynomial in the variables first (inclusive) up to
3459 * last (exclusive) by inserting powers of variable first.
3460 * Variable first is assumed not to appear in the input.
3461 */
3462 __isl_give struct isl_upoly *isl_upoly_homogenize(
3463 __isl_take struct isl_upoly *up, int deg, int target,
3464 int first, int last)
3465 {
3466 int i;
3467 struct isl_upoly_rec *rec;
3468
3469 if (!up)
3470 return NULL;
3471 if (isl_upoly_is_zero(up))
3472 return up;
3473 if (deg == target)
3474 return up;
3475 if (isl_upoly_is_cst(up) || up->var < first) {
3476 struct isl_upoly *hom;
3477
3478 hom = isl_upoly_var_pow(up->ctx, first, target - deg);
3479 if (!hom)
3480 goto error;
3481 rec = isl_upoly_as_rec(hom);
3482 rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
3483
3484 return hom;
3485 }
3486
3487 up = isl_upoly_cow(up);
3488 rec = isl_upoly_as_rec(up);
3489 if (!rec)
3490 goto error;
3491
3492 for (i = 0; i < rec->n; ++i) {
3493 if (isl_upoly_is_zero(rec->p[i]))
3494 continue;
3495 rec->p[i] = isl_upoly_homogenize(rec->p[i],
3496 up->var < last ? deg + i : i, target,
3497 first, last);
3498 if (!rec->p[i])
3499 goto error;
3500 }
3501
3502 return up;
3503 error:
3504 isl_upoly_free(up);
3505 return NULL;
3506 }
3507
3508 /* Homogenize the polynomial in the set variables by introducing
3509 * powers of an extra set variable at position 0.
3510 */
3511 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
3512 __isl_take isl_qpolynomial *poly)
3513 {
3514 unsigned ovar;
3515 unsigned nvar;
3516 int deg = isl_qpolynomial_degree(poly);
3517
3518 if (deg < -1)
3519 goto error;
3520
3521 poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1);
3522 poly = isl_qpolynomial_cow(poly);
3523 if (!poly)
3524 goto error;
3525
3526 ovar = isl_space_offset(poly->dim, isl_dim_set);
3527 nvar = isl_space_dim(poly->dim, isl_dim_set);
3528 poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
3529 ovar, ovar + nvar);
3530 if (!poly->upoly)
3531 goto error;
3532
3533 return poly;
3534 error:
3535 isl_qpolynomial_free(poly);
3536 return NULL;
3537 }
3538
3539 __isl_give isl_term *isl_term_alloc(__isl_take isl_space *dim,
3540 __isl_take isl_mat *div)
3541 {
3542 isl_term *term;
3543 int n;
3544
3545 if (!dim || !div)
3546 goto error;
3547
3548 n = isl_space_dim(dim, isl_dim_all) + div->n_row;
3549
3550 term = isl_calloc(dim->ctx, struct isl_term,
3551 sizeof(struct isl_term) + (n - 1) * sizeof(int));
3552 if (!term)
3553 goto error;
3554
3555 term->ref = 1;
3556 term->dim = dim;
3557 term->div = div;
3558 isl_int_init(term->n);
3559 isl_int_init(term->d);
3560
3561 return term;
3562 error:
3563 isl_space_free(dim);
3564 isl_mat_free(div);
3565 return NULL;
3566 }
3567
3568 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
3569 {
3570 if (!term)
3571 return NULL;
3572
3573 term->ref++;
3574 return term;
3575 }
3576
3577 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
3578 {
3579 int i;
3580 isl_term *dup;
3581 unsigned total;
3582
3583 if (!term)
3584 return NULL;
3585
3586 total = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3587
3588 dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div));
3589 if (!dup)
3590 return NULL;
3591
3592 isl_int_set(dup->n, term->n);
3593 isl_int_set(dup->d, term->d);
3594
3595 for (i = 0; i < total; ++i)
3596 dup->pow[i] = term->pow[i];
3597
3598 return dup;
3599 }
3600
3601 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
3602 {
3603 if (!term)
3604 return NULL;
3605
3606 if (term->ref == 1)
3607 return term;
3608 term->ref--;
3609 return isl_term_dup(term);
3610 }
3611
3612 void isl_term_free(__isl_take isl_term *term)
3613 {
3614 if (!term)
3615 return;
3616
3617 if (--term->ref > 0)
3618 return;
3619
3620 isl_space_free(term->dim);
3621 isl_mat_free(term->div);
3622 isl_int_clear(term->n);
3623 isl_int_clear(term->d);
3624 free(term);
3625 }
3626
3627 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
3628 {
3629 if (!term)
3630 return 0;
3631
3632 switch (type) {
3633 case isl_dim_param:
3634 case isl_dim_in:
3635 case isl_dim_out: return isl_space_dim(term->dim, type);
3636 case isl_dim_div: return term->div->n_row;
3637 case isl_dim_all: return isl_space_dim(term->dim, isl_dim_all) +
3638 term->div->n_row;
3639 default: return 0;
3640 }
3641 }
3642
3643 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
3644 {
3645 return term ? term->dim->ctx : NULL;
3646 }
3647
3648 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
3649 {
3650 if (!term)
3651 return;
3652 isl_int_set(*n, term->n);
3653 }
3654
3655 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
3656 {
3657 if (!term)
3658 return;
3659 isl_int_set(*d, term->d);
3660 }
3661
3662 /* Return the coefficient of the term "term".
3663 */
3664 __isl_give isl_val *isl_term_get_coefficient_val(__isl_keep isl_term *term)
3665 {
3666 if (!term)
3667 return NULL;
3668
3669 return isl_val_rat_from_isl_int(isl_term_get_ctx(term),
3670 term->n, term->d);
3671 }
3672
3673 int isl_term_get_exp(__isl_keep isl_term *term,
3674 enum isl_dim_type type, unsigned pos)
3675 {
3676 if (!term)
3677 return -1;
3678
3679 isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
3680
3681 if (type >= isl_dim_set)
3682 pos += isl_space_dim(term->dim, isl_dim_param);
3683 if (type >= isl_dim_div)
3684 pos += isl_space_dim(term->dim, isl_dim_set);
3685
3686 return term->pow[pos];
3687 }
3688
3689 __isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
3690 {
3691 isl_local_space *ls;
3692 isl_aff *aff;
3693
3694 if (!term)
3695 return NULL;
3696
3697 isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
3698 return NULL);
3699
3700 ls = isl_local_space_alloc_div(isl_space_copy(term->dim),
3701 isl_mat_copy(term->div));
3702 aff = isl_aff_alloc(ls);
3703 if (!aff)
3704 return NULL;
3705
3706 isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size);
3707
3708 aff = isl_aff_normalize(aff);
3709
3710 return aff;
3711 }
3712
3713 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
3714 int (*fn)(__isl_take isl_term *term, void *user),
3715 __isl_take isl_term *term, void *user)
3716 {
3717 int i;
3718 struct isl_upoly_rec *rec;
3719
3720 if (!up || !term)
3721 goto error;
3722
3723 if (isl_upoly_is_zero(up))
3724 return term;
3725
3726 isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
3727 isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
3728 isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
3729
3730 if (isl_upoly_is_cst(up)) {
3731 struct isl_upoly_cst *cst;
3732 cst = isl_upoly_as_cst(up);
3733 if (!cst)
3734 goto error;
3735 term = isl_term_cow(term);
3736 if (!term)
3737 goto error;
3738 isl_int_set(term->n, cst->n);
3739 isl_int_set(term->d, cst->d);
3740 if (fn(isl_term_copy(term), user) < 0)
3741 goto error;
3742 return term;
3743 }
3744
3745 rec = isl_upoly_as_rec(up);
3746 if (!rec)
3747 goto error;
3748
3749 for (i = 0; i < rec->n; ++i) {
3750 term = isl_term_cow(term);
3751 if (!term)
3752 goto error;
3753 term->pow[up->var] = i;
3754 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
3755 if (!term)
3756 goto error;
3757 }
3758 term->pow[up->var] = 0;
3759
3760 return term;
3761 error:
3762 isl_term_free(term);
3763 return NULL;
3764 }
3765
3766 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
3767 int (*fn)(__isl_take isl_term *term, void *user), void *user)
3768 {
3769 isl_term *term;
3770
3771 if (!qp)
3772 return -1;
3773
3774 term = isl_term_alloc(isl_space_copy(qp->dim), isl_mat_copy(qp->div));
3775 if (!term)
3776 return -1;
3777
3778 term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
3779
3780 isl_term_free(term);
3781
3782 return term ? 0 : -1;
3783 }
3784
3785 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
3786 {
3787 struct isl_upoly *up;
3788 isl_qpolynomial *qp;
3789 int i, n;
3790
3791 if (!term)
3792 return NULL;
3793
3794 n = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3795
3796 up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
3797 for (i = 0; i < n; ++i) {
3798 if (!term->pow[i])
3799 continue;
3800 up = isl_upoly_mul(up,
3801 isl_upoly_var_pow(term->dim->ctx, i, term->pow[i]));
3802 }
3803
3804 qp = isl_qpolynomial_alloc(isl_space_copy(term->dim), term->div->n_row, up);
3805 if (!qp)
3806 goto error;
3807 isl_mat_free(qp->div);
3808 qp->div = isl_mat_copy(term->div);
3809 if (!qp->div)
3810 goto error;
3811
3812 isl_term_free(term);
3813 return qp;
3814 error:
3815 isl_qpolynomial_free(qp);
3816 isl_term_free(term);
3817 return NULL;
3818 }
3819
3820 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
3821 __isl_take isl_space *dim)
3822 {
3823 int i;
3824 int extra;
3825 unsigned total;
3826
3827 if (!qp || !dim)
3828 goto error;
3829
3830 if (isl_space_is_equal(qp->dim, dim)) {
3831 isl_space_free(dim);
3832 return qp;
3833 }
3834
3835 qp = isl_qpolynomial_cow(qp);
3836 if (!qp)
3837 goto error;
3838
3839 extra = isl_space_dim(dim, isl_dim_set) -
3840 isl_space_dim(qp->dim, isl_dim_set);
3841 total = isl_space_dim(qp->dim, isl_dim_all);
3842 if (qp->div->n_row) {
3843 int *exp;
3844
3845 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
3846 if (!exp)
3847 goto error;
3848 for (i = 0; i < qp->div->n_row; ++i)
3849 exp[i] = extra + i;
3850 qp->upoly = expand(qp->upoly, exp, total);
3851 free(exp);
3852 if (!qp->upoly)
3853 goto error;
3854 }
3855 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
3856 if (!qp->div)
3857 goto error;
3858 for (i = 0; i < qp->div->n_row; ++i)
3859 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
3860
3861 isl_space_free(qp->dim);
3862 qp->dim = dim;
3863
3864 return qp;
3865 error:
3866 isl_space_free(dim);
3867 isl_qpolynomial_free(qp);
3868 return NULL;
3869 }
3870
3871 /* For each parameter or variable that does not appear in qp,
3872 * first eliminate the variable from all constraints and then set it to zero.
3873 */
3874 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
3875 __isl_keep isl_qpolynomial *qp)
3876 {
3877 int *active = NULL;
3878 int i;
3879 int d;
3880 unsigned nparam;
3881 unsigned nvar;
3882
3883 if (!set || !qp)
3884 goto error;
3885
3886 d = isl_space_dim(set->dim, isl_dim_all);
3887 active = isl_calloc_array(set->ctx, int, d);
3888 if (set_active(qp, active) < 0)
3889 goto error;
3890
3891 for (i = 0; i < d; ++i)
3892 if (!active[i])
3893 break;
3894
3895 if (i == d) {
3896 free(active);
3897 return set;
3898 }
3899
3900 nparam = isl_space_dim(set->dim, isl_dim_param);
3901 nvar = isl_space_dim(set->dim, isl_dim_set);
3902 for (i = 0; i < nparam; ++i) {
3903 if (active[i])
3904 continue;
3905 set = isl_set_eliminate(set, isl_dim_param, i, 1);
3906 set = isl_set_fix_si(set, isl_dim_param, i, 0);
3907 }
3908 for (i = 0; i < nvar; ++i) {
3909 if (active[nparam + i])
3910 continue;
3911 set = isl_set_eliminate(set, isl_dim_set, i, 1);
3912 set = isl_set_fix_si(set, isl_dim_set, i, 0);
3913 }
3914
3915 free(active);
3916
3917 return set;
3918 error:
3919 free(active);
3920 isl_set_free(set);
3921 return NULL;
3922 }
3923
3924 struct isl_opt_data {
3925 isl_qpolynomial *qp;
3926 int first;
3927 isl_val *opt;
3928 int max;
3929 };
3930
3931 static int opt_fn(__isl_take isl_point *pnt, void *user)
3932 {
3933 struct isl_opt_data *data = (struct isl_opt_data *)user;
3934 isl_val *val;
3935
3936 val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
3937 if (data->first) {
3938 data->first = 0;
3939 data->opt = val;
3940 } else if (data->max) {
3941 data->opt = isl_val_max(data->opt, val);
3942 } else {
3943 data->opt = isl_val_min(data->opt, val);
3944 }
3945
3946 return 0;
3947 }
3948
3949 __isl_give isl_val *isl_qpolynomial_opt_on_domain(
3950 __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
3951 {
3952 struct isl_opt_data data = { NULL, 1, NULL, max };
3953
3954 if (!set || !qp)
3955 goto error;
3956
3957 if (isl_upoly_is_cst(qp->upoly)) {
3958 isl_set_free(set);
3959 data.opt = isl_qpolynomial_get_constant_val(qp);
3960 isl_qpolynomial_free(qp);
3961 return data.opt;
3962 }
3963
3964 set = fix_inactive(set, qp);
3965
3966 data.qp = qp;
3967 if (isl_set_foreach_point(set, opt_fn, &data) < 0)
3968 goto error;
3969
3970 if (data.first)
3971 data.opt = isl_val_zero(isl_set_get_ctx(set));
3972
3973 isl_set_free(set);
3974 isl_qpolynomial_free(qp);
3975 return data.opt;
3976 error:
3977 isl_set_free(set);
3978 isl_qpolynomial_free(qp);
3979 isl_val_free(data.opt);
3980 return NULL;
3981 }
3982
3983 __isl_give isl_qpolynomial *isl_qpolynomial_morph_domain(
3984 __isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph)
3985 {
3986 int i;
3987 int n_sub;
3988 isl_ctx *ctx;
3989 struct isl_upoly **subs;
3990 isl_mat *mat, *diag;
3991
3992 qp = isl_qpolynomial_cow(qp);
3993 if (!qp || !morph)
3994 goto error;
3995
3996 ctx = qp->dim->ctx;
3997 isl_assert(ctx, isl_space_is_equal(qp->dim, morph->dom->dim), goto error);
3998
3999 n_sub = morph->inv->n_row - 1;
4000 if (morph->inv->n_row != morph->inv->n_col)
4001 n_sub += qp->div->n_row;
4002 subs = isl_calloc_array(ctx, struct isl_upoly *, n_sub);
4003 if (n_sub && !subs)
4004 goto error;
4005
4006 for (i = 0; 1 + i < morph->inv->n_row; ++i)
4007 subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
4008 morph->inv->row[0][0], morph->inv->n_col);
4009 if (morph->inv->n_row != morph->inv->n_col)
4010 for (i = 0; i < qp->div->n_row; ++i)
4011 subs[morph->inv->n_row - 1 + i] =
4012 isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
4013
4014 qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
4015
4016 for (i = 0; i < n_sub; ++i)
4017 isl_upoly_free(subs[i]);
4018 free(subs);
4019
4020 diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
4021 mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
4022 diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
4023 mat = isl_mat_diagonal(mat, diag);
4024 qp->div = isl_mat_product(qp->div, mat);
4025 isl_space_free(qp->dim);
4026 qp->dim = isl_space_copy(morph->ran->dim);
4027
4028 if (!qp->upoly || !qp->div || !qp->dim)
4029 goto error;
4030
4031 isl_morph_free(morph);
4032
4033 return qp;
4034 error:
4035 isl_qpolynomial_free(qp);
4036 isl_morph_free(morph);
4037 return NULL;
4038 }
4039
4040 static int neg_entry(void **entry, void *user)
4041 {
4042 isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
4043
4044 *pwqp = isl_pw_qpolynomial_neg(*pwqp);
4045
4046 return *pwqp ? 0 : -1;
4047 }
4048
4049 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_neg(
4050 __isl_take isl_union_pw_qpolynomial *upwqp)
4051 {
4052 upwqp = isl_union_pw_qpolynomial_cow(upwqp);
4053 if (!upwqp)
4054 return NULL;
4055
4056 if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
4057 &neg_entry, NULL) < 0)
4058 goto error;
4059
4060 return upwqp;
4061 error:
4062 isl_union_pw_qpolynomial_free(upwqp);
4063 return NULL;
4064 }
4065
4066 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
4067 __isl_take isl_union_pw_qpolynomial *upwqp1,
4068 __isl_take isl_union_pw_qpolynomial *upwqp2)
4069 {
4070 return match_bin_op(upwqp1, upwqp2, &isl_pw_qpolynomial_mul);
4071 }
4072
4073 /* Reorder the columns of the given div definitions according to the
4074 * given reordering.
4075 */
4076 static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
4077 __isl_take isl_reordering *r)
4078 {
4079 int i, j;
4080 isl_mat *mat;
4081 int extra;
4082
4083 if (!div || !r)
4084 goto error;
4085
4086 extra = isl_space_dim(r->dim, isl_dim_all) + div->n_row - r->len;
4087 mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
4088 if (!mat)
4089 goto error;
4090
4091 for (i = 0; i < div->n_row; ++i) {
4092 isl_seq_cpy(mat->row[i], div->row[i], 2);
4093 isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
4094 for (j = 0; j < r->len; ++j)
4095 isl_int_set(mat->row[i][2 + r->pos[j]],
4096 div->row[i][2 + j]);
4097 }
4098
4099 isl_reordering_free(r);
4100 isl_mat_free(div);
4101 return mat;
4102 error:
4103 isl_reordering_free(r);
4104 isl_mat_free(div);
4105 return NULL;
4106 }
4107
4108 /* Reorder the dimension of "qp" according to the given reordering.
4109 */
4110 __isl_give isl_qpolynomial *isl_qpolynomial_realign_domain(
4111 __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
4112 {
4113 qp = isl_qpolynomial_cow(qp);
4114 if (!qp)
4115 goto error;
4116
4117 r = isl_reordering_extend(r, qp->div->n_row);
4118 if (!r)
4119 goto error;
4120
4121 qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
4122 if (!qp->div)
4123 goto error;
4124
4125 qp->upoly = reorder(qp->upoly, r->pos);
4126 if (!qp->upoly)
4127 goto error;
4128
4129 qp = isl_qpolynomial_reset_domain_space(qp, isl_space_copy(r->dim));
4130
4131 isl_reordering_free(r);
4132 return qp;
4133 error:
4134 isl_qpolynomial_free(qp);
4135 isl_reordering_free(r);
4136 return NULL;
4137 }
4138
4139 __isl_give isl_qpolynomial *isl_qpolynomial_align_params(
4140 __isl_take isl_qpolynomial *qp, __isl_take isl_space *model)
4141 {
4142 if (!qp || !model)
4143 goto error;
4144
4145 if (!isl_space_match(qp->dim, isl_dim_param, model, isl_dim_param)) {
4146 isl_reordering *exp;
4147
4148 model = isl_space_drop_dims(model, isl_dim_in,
4149 0, isl_space_dim(model, isl_dim_in));
4150 model = isl_space_drop_dims(model, isl_dim_out,
4151 0, isl_space_dim(model, isl_dim_out));
4152 exp = isl_parameter_alignment_reordering(qp->dim, model);
4153 exp = isl_reordering_extend_space(exp,
4154 isl_qpolynomial_get_domain_space(qp));
4155 qp = isl_qpolynomial_realign_domain(qp, exp);
4156 }
4157
4158 isl_space_free(model);
4159 return qp;
4160 error:
4161 isl_space_free(model);
4162 isl_qpolynomial_free(qp);
4163 return NULL;
4164 }
4165
4166 struct isl_split_periods_data {
4167 int max_periods;
4168 isl_pw_qpolynomial *res;
4169 };
4170
4171 /* Create a slice where the integer division "div" has the fixed value "v".
4172 * In particular, if "div" refers to floor(f/m), then create a slice
4173 *
4174 * m v <= f <= m v + (m - 1)
4175 *
4176 * or
4177 *
4178 * f - m v >= 0
4179 * -f + m v + (m - 1) >= 0
4180 */
4181 static __isl_give isl_set *set_div_slice(__isl_take isl_space *dim,
4182 __isl_keep isl_qpolynomial *qp, int div, isl_int v)
4183 {
4184 int total;
4185 isl_basic_set *bset = NULL;
4186 int k;
4187
4188 if (!dim || !qp)
4189 goto error;
4190
4191 total = isl_space_dim(dim, isl_dim_all);
4192 bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0, 0, 2);
4193
4194 k = isl_basic_set_alloc_inequality(bset);
4195 if (k < 0)
4196 goto error;
4197 isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4198 isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
4199
4200 k = isl_basic_set_alloc_inequality(bset);
4201 if (k < 0)
4202 goto error;
4203 isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4204 isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
4205 isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
4206 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4207
4208 isl_space_free(dim);
4209 return isl_set_from_basic_set(bset);
4210 error:
4211 isl_basic_set_free(bset);
4212 isl_space_free(dim);
4213 return NULL;
4214 }
4215
4216 static int split_periods(__isl_take isl_set *set,
4217 __isl_take isl_qpolynomial *qp, void *user);
4218
4219 /* Create a slice of the domain "set" such that integer division "div"
4220 * has the fixed value "v" and add the results to data->res,
4221 * replacing the integer division by "v" in "qp".
4222 */
4223 static int set_div(__isl_take isl_set *set,
4224 __isl_take isl_qpolynomial *qp, int div, isl_int v,
4225 struct isl_split_periods_data *data)
4226 {
4227 int i;
4228 int total;
4229 isl_set *slice;
4230 struct isl_upoly *cst;
4231
4232 slice = set_div_slice(isl_set_get_space(set), qp, div, v);
4233 set = isl_set_intersect(set, slice);
4234
4235 if (!qp)
4236 goto error;
4237
4238 total = isl_space_dim(qp->dim, isl_dim_all);
4239
4240 for (i = div + 1; i < qp->div->n_row; ++i) {
4241 if (isl_int_is_zero(qp->div->row[i][2 + total + div]))
4242 continue;
4243 isl_int_addmul(qp->div->row[i][1],
4244 qp->div->row[i][2 + total + div], v);
4245 isl_int_set_si(qp->div->row[i][2 + total + div], 0);
4246 }
4247
4248 cst = isl_upoly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
4249 qp = substitute_div(qp, div, cst);
4250
4251 return split_periods(set, qp, data);
4252 error:
4253 isl_set_free(set);
4254 isl_qpolynomial_free(qp);
4255 return -1;
4256 }
4257
4258 /* Split the domain "set" such that integer division "div"
4259 * has a fixed value (ranging from "min" to "max") on each slice
4260 * and add the results to data->res.
4261 */
4262 static int split_div(__isl_take isl_set *set,
4263 __isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
4264 struct isl_split_periods_data *data)
4265 {
4266 for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
4267 isl_set *set_i = isl_set_copy(set);
4268 isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
4269
4270 if (set_div(set_i, qp_i, div, min, data) < 0)
4271 goto error;
4272 }
4273 isl_set_free(set);
4274 isl_qpolynomial_free(qp);
4275 return 0;
4276 error:
4277 isl_set_free(set);
4278 isl_qpolynomial_free(qp);
4279 return -1;
4280 }
4281
4282 /* If "qp" refers to any integer division
4283 * that can only attain "max_periods" distinct values on "set"
4284 * then split the domain along those distinct values.
4285 * Add the results (or the original if no splitting occurs)
4286 * to data->res.
4287 */
4288 static int split_periods(__isl_take isl_set *set,
4289 __isl_take isl_qpolynomial *qp, void *user)
4290 {
4291 int i;
4292 isl_pw_qpolynomial *pwqp;
4293 struct isl_split_periods_data *data;
4294 isl_int min, max;
4295 int total;
4296 int r = 0;
4297
4298 data = (struct isl_split_periods_data *)user;
4299
4300 if (!set || !qp)
4301 goto error;
4302
4303 if (qp->div->n_row == 0) {
4304 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4305 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4306 return 0;
4307 }
4308
4309 isl_int_init(min);
4310 isl_int_init(max);
4311 total = isl_space_dim(qp->dim, isl_dim_all);
4312 for (i = 0; i < qp->div->n_row; ++i) {
4313 enum isl_lp_result lp_res;
4314
4315 if (isl_seq_first_non_zero(qp->div->row[i] + 2 + total,
4316 qp->div->n_row) != -1)
4317 continue;
4318
4319 lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
4320 set->ctx->one, &min, 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(min, min, qp->div->row[i][0]);
4326
4327 lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
4328 set->ctx->one, &max, NULL, NULL);
4329 if (lp_res == isl_lp_error)
4330 goto error2;
4331 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4332 continue;
4333 isl_int_fdiv_q(max, max, qp->div->row[i][0]);
4334
4335 isl_int_sub(max, max, min);
4336 if (isl_int_cmp_si(max, data->max_periods) < 0) {
4337 isl_int_add(max, max, min);
4338 break;
4339 }
4340 }
4341
4342 if (i < qp->div->n_row) {
4343 r = split_div(set, qp, i, min, max, data);
4344 } else {
4345 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4346 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4347 }
4348
4349 isl_int_clear(max);
4350 isl_int_clear(min);
4351
4352 return r;
4353 error2:
4354 isl_int_clear(max);
4355 isl_int_clear(min);
4356 error:
4357 isl_set_free(set);
4358 isl_qpolynomial_free(qp);
4359 return -1;
4360 }
4361
4362 /* If any quasi-polynomial in pwqp refers to any integer division
4363 * that can only attain "max_periods" distinct values on its domain
4364 * then split the domain along those distinct values.
4365 */
4366 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
4367 __isl_take isl_pw_qpolynomial *pwqp, int max_periods)
4368 {
4369 struct isl_split_periods_data data;
4370
4371 data.max_periods = max_periods;
4372 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4373
4374 if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
4375 goto error;
4376
4377 isl_pw_qpolynomial_free(pwqp);
4378
4379 return data.res;
4380 error:
4381 isl_pw_qpolynomial_free(data.res);
4382 isl_pw_qpolynomial_free(pwqp);
4383 return NULL;
4384 }
4385
4386 /* Construct a piecewise quasipolynomial that is constant on the given
4387 * domain. In particular, it is
4388 * 0 if cst == 0
4389 * 1 if cst == 1
4390 * infinity if cst == -1
4391 */
4392 static __isl_give isl_pw_qpolynomial *constant_on_domain(
4393 __isl_take isl_basic_set *bset, int cst)
4394 {
4395 isl_space *dim;
4396 isl_qpolynomial *qp;
4397
4398 if (!bset)
4399 return NULL;
4400
4401 bset = isl_basic_set_params(bset);
4402 dim = isl_basic_set_get_space(bset);
4403 if (cst < 0)
4404 qp = isl_qpolynomial_infty_on_domain(dim);
4405 else if (cst == 0)
4406 qp = isl_qpolynomial_zero_on_domain(dim);
4407 else
4408 qp = isl_qpolynomial_one_on_domain(dim);
4409 return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
4410 }
4411
4412 /* Factor bset, call fn on each of the factors and return the product.
4413 *
4414 * If no factors can be found, simply call fn on the input.
4415 * Otherwise, construct the factors based on the factorizer,
4416 * call fn on each factor and compute the product.
4417 */
4418 static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
4419 __isl_take isl_basic_set *bset,
4420 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4421 {
4422 int i, n;
4423 isl_space *dim;
4424 isl_set *set;
4425 isl_factorizer *f;
4426 isl_qpolynomial *qp;
4427 isl_pw_qpolynomial *pwqp;
4428 unsigned nparam;
4429 unsigned nvar;
4430
4431 f = isl_basic_set_factorizer(bset);
4432 if (!f)
4433 goto error;
4434 if (f->n_group == 0) {
4435 isl_factorizer_free(f);
4436 return fn(bset);
4437 }
4438
4439 nparam = isl_basic_set_dim(bset, isl_dim_param);
4440 nvar = isl_basic_set_dim(bset, isl_dim_set);
4441
4442 dim = isl_basic_set_get_space(bset);
4443 dim = isl_space_domain(dim);
4444 set = isl_set_universe(isl_space_copy(dim));
4445 qp = isl_qpolynomial_one_on_domain(dim);
4446 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4447
4448 bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
4449
4450 for (i = 0, n = 0; i < f->n_group; ++i) {
4451 isl_basic_set *bset_i;
4452 isl_pw_qpolynomial *pwqp_i;
4453
4454 bset_i = isl_basic_set_copy(bset);
4455 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4456 nparam + n + f->len[i], nvar - n - f->len[i]);
4457 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4458 nparam, n);
4459 bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
4460 n + f->len[i], nvar - n - f->len[i]);
4461 bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
4462
4463 pwqp_i = fn(bset_i);
4464 pwqp = isl_pw_qpolynomial_mul(pwqp, pwqp_i);
4465
4466 n += f->len[i];
4467 }
4468
4469 isl_basic_set_free(bset);
4470 isl_factorizer_free(f);
4471
4472 return pwqp;
4473 error:
4474 isl_basic_set_free(bset);
4475 return NULL;
4476 }
4477
4478 /* Factor bset, call fn on each of the factors and return the product.
4479 * The function is assumed to evaluate to zero on empty domains,
4480 * to one on zero-dimensional domains and to infinity on unbounded domains
4481 * and will not be called explicitly on zero-dimensional or unbounded domains.
4482 *
4483 * We first check for some special cases and remove all equalities.
4484 * Then we hand over control to compressed_multiplicative_call.
4485 */
4486 __isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
4487 __isl_take isl_basic_set *bset,
4488 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4489 {
4490 int bounded;
4491 isl_morph *morph;
4492 isl_pw_qpolynomial *pwqp;
4493
4494 if (!bset)
4495 return NULL;
4496
4497 if (isl_basic_set_plain_is_empty(bset))
4498 return constant_on_domain(bset, 0);
4499
4500 if (isl_basic_set_dim(bset, isl_dim_set) == 0)
4501 return constant_on_domain(bset, 1);
4502
4503 bounded = isl_basic_set_is_bounded(bset);
4504 if (bounded < 0)
4505 goto error;
4506 if (!bounded)
4507 return constant_on_domain(bset, -1);
4508
4509 if (bset->n_eq == 0)
4510 return compressed_multiplicative_call(bset, fn);
4511
4512 morph = isl_basic_set_full_compression(bset);
4513 bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
4514
4515 pwqp = compressed_multiplicative_call(bset, fn);
4516
4517 morph = isl_morph_dom_params(morph);
4518 morph = isl_morph_ran_params(morph);
4519 morph = isl_morph_inverse(morph);
4520
4521 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph);
4522
4523 return pwqp;
4524 error:
4525 isl_basic_set_free(bset);
4526 return NULL;
4527 }
4528
4529 /* Drop all floors in "qp", turning each integer division [a/m] into
4530 * a rational division a/m. If "down" is set, then the integer division
4531 * is replaced by (a-(m-1))/m instead.
4532 */
4533 static __isl_give isl_qpolynomial *qp_drop_floors(
4534 __isl_take isl_qpolynomial *qp, int down)
4535 {
4536 int i;
4537 struct isl_upoly *s;
4538
4539 if (!qp)
4540 return NULL;
4541 if (qp->div->n_row == 0)
4542 return qp;
4543
4544 qp = isl_qpolynomial_cow(qp);
4545 if (!qp)
4546 return NULL;
4547
4548 for (i = qp->div->n_row - 1; i >= 0; --i) {
4549 if (down) {
4550 isl_int_sub(qp->div->row[i][1],
4551 qp->div->row[i][1], qp->div->row[i][0]);
4552 isl_int_add_ui(qp->div->row[i][1],
4553 qp->div->row[i][1], 1);
4554 }
4555 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
4556 qp->div->row[i][0], qp->div->n_col - 1);
4557 qp = substitute_div(qp, i, s);
4558 if (!qp)
4559 return NULL;
4560 }
4561
4562 return qp;
4563 }
4564
4565 /* Drop all floors in "pwqp", turning each integer division [a/m] into
4566 * a rational division a/m.
4567 */
4568 static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
4569 __isl_take isl_pw_qpolynomial *pwqp)
4570 {
4571 int i;
4572
4573 if (!pwqp)
4574 return NULL;
4575
4576 if (isl_pw_qpolynomial_is_zero(pwqp))
4577 return pwqp;
4578
4579 pwqp = isl_pw_qpolynomial_cow(pwqp);
4580 if (!pwqp)
4581 return NULL;
4582
4583 for (i = 0; i < pwqp->n; ++i) {
4584 pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
4585 if (!pwqp->p[i].qp)
4586 goto error;
4587 }
4588
4589 return pwqp;
4590 error:
4591 isl_pw_qpolynomial_free(pwqp);
4592 return NULL;
4593 }
4594
4595 /* Adjust all the integer divisions in "qp" such that they are at least
4596 * one over the given orthant (identified by "signs"). This ensures
4597 * that they will still be non-negative even after subtracting (m-1)/m.
4598 *
4599 * In particular, f is replaced by f' + v, changing f = [a/m]
4600 * to f' = [(a - m v)/m].
4601 * If the constant term k in a is smaller than m,
4602 * the constant term of v is set to floor(k/m) - 1.
4603 * For any other term, if the coefficient c and the variable x have
4604 * the same sign, then no changes are needed.
4605 * Otherwise, if the variable is positive (and c is negative),
4606 * then the coefficient of x in v is set to floor(c/m).
4607 * If the variable is negative (and c is positive),
4608 * then the coefficient of x in v is set to ceil(c/m).
4609 */
4610 static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
4611 int *signs)
4612 {
4613 int i, j;
4614 int total;
4615 isl_vec *v = NULL;
4616 struct isl_upoly *s;
4617
4618 qp = isl_qpolynomial_cow(qp);
4619 if (!qp)
4620 return NULL;
4621 qp->div = isl_mat_cow(qp->div);
4622 if (!qp->div)
4623 goto error;
4624
4625 total = isl_space_dim(qp->dim, isl_dim_all);
4626 v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
4627
4628 for (i = 0; i < qp->div->n_row; ++i) {
4629 isl_int *row = qp->div->row[i];
4630 v = isl_vec_clr(v);
4631 if (!v)
4632 goto error;
4633 if (isl_int_lt(row[1], row[0])) {
4634 isl_int_fdiv_q(v->el[0], row[1], row[0]);
4635 isl_int_sub_ui(v->el[0], v->el[0], 1);
4636 isl_int_submul(row[1], row[0], v->el[0]);
4637 }
4638 for (j = 0; j < total; ++j) {
4639 if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
4640 continue;
4641 if (signs[j] < 0)
4642 isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
4643 else
4644 isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
4645 isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
4646 }
4647 for (j = 0; j < i; ++j) {
4648 if (isl_int_sgn(row[2 + total + j]) >= 0)
4649 continue;
4650 isl_int_fdiv_q(v->el[1 + total + j],
4651 row[2 + total + j], row[0]);
4652 isl_int_submul(row[2 + total + j],
4653 row[0], v->el[1 + total + j]);
4654 }
4655 for (j = i + 1; j < qp->div->n_row; ++j) {
4656 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
4657 continue;
4658 isl_seq_combine(qp->div->row[j] + 1,
4659 qp->div->ctx->one, qp->div->row[j] + 1,
4660 qp->div->row[j][2 + total + i], v->el, v->size);
4661 }
4662 isl_int_set_si(v->el[1 + total + i], 1);
4663 s = isl_upoly_from_affine(qp->dim->ctx, v->el,
4664 qp->div->ctx->one, v->size);
4665 qp->upoly = isl_upoly_subs(qp->upoly, total + i, 1, &s);
4666 isl_upoly_free(s);
4667 if (!qp->upoly)
4668 goto error;
4669 }
4670
4671 isl_vec_free(v);
4672 return qp;
4673 error:
4674 isl_vec_free(v);
4675 isl_qpolynomial_free(qp);
4676 return NULL;
4677 }
4678
4679 struct isl_to_poly_data {
4680 int sign;
4681 isl_pw_qpolynomial *res;
4682 isl_qpolynomial *qp;
4683 };
4684
4685 /* Appoximate data->qp by a polynomial on the orthant identified by "signs".
4686 * We first make all integer divisions positive and then split the
4687 * quasipolynomials into terms with sign data->sign (the direction
4688 * of the requested approximation) and terms with the opposite sign.
4689 * In the first set of terms, each integer division [a/m] is
4690 * overapproximated by a/m, while in the second it is underapproximated
4691 * by (a-(m-1))/m.
4692 */
4693 static int to_polynomial_on_orthant(__isl_take isl_set *orthant, int *signs,
4694 void *user)
4695 {
4696 struct isl_to_poly_data *data = user;
4697 isl_pw_qpolynomial *t;
4698 isl_qpolynomial *qp, *up, *down;
4699
4700 qp = isl_qpolynomial_copy(data->qp);
4701 qp = make_divs_pos(qp, signs);
4702
4703 up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
4704 up = qp_drop_floors(up, 0);
4705 down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
4706 down = qp_drop_floors(down, 1);
4707
4708 isl_qpolynomial_free(qp);
4709 qp = isl_qpolynomial_add(up, down);
4710
4711 t = isl_pw_qpolynomial_alloc(orthant, qp);
4712 data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
4713
4714 return 0;
4715 }
4716
4717 /* Approximate each quasipolynomial by a polynomial. If "sign" is positive,
4718 * the polynomial will be an overapproximation. If "sign" is negative,
4719 * it will be an underapproximation. If "sign" is zero, the approximation
4720 * will lie somewhere in between.
4721 *
4722 * In particular, is sign == 0, we simply drop the floors, turning
4723 * the integer divisions into rational divisions.
4724 * Otherwise, we split the domains into orthants, make all integer divisions
4725 * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
4726 * depending on the requested sign and the sign of the term in which
4727 * the integer division appears.
4728 */
4729 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
4730 __isl_take isl_pw_qpolynomial *pwqp, int sign)
4731 {
4732 int i;
4733 struct isl_to_poly_data data;
4734
4735 if (sign == 0)
4736 return pwqp_drop_floors(pwqp);
4737
4738 if (!pwqp)
4739 return NULL;
4740
4741 data.sign = sign;
4742 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4743
4744 for (i = 0; i < pwqp->n; ++i) {
4745 if (pwqp->p[i].qp->div->n_row == 0) {
4746 isl_pw_qpolynomial *t;
4747 t = isl_pw_qpolynomial_alloc(
4748 isl_set_copy(pwqp->p[i].set),
4749 isl_qpolynomial_copy(pwqp->p[i].qp));
4750 data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
4751 continue;
4752 }
4753 data.qp = pwqp->p[i].qp;
4754 if (isl_set_foreach_orthant(pwqp->p[i].set,
4755 &to_polynomial_on_orthant, &data) < 0)
4756 goto error;
4757 }
4758
4759 isl_pw_qpolynomial_free(pwqp);
4760
4761 return data.res;
4762 error:
4763 isl_pw_qpolynomial_free(pwqp);
4764 isl_pw_qpolynomial_free(data.res);
4765 return NULL;
4766 }
4767
4768 static int poly_entry(void **entry, void *user)
4769 {
4770 int *sign = user;
4771 isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
4772
4773 *pwqp = isl_pw_qpolynomial_to_polynomial(*pwqp, *sign);
4774
4775 return *pwqp ? 0 : -1;
4776 }
4777
4778 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
4779 __isl_take isl_union_pw_qpolynomial *upwqp, int sign)
4780 {
4781 upwqp = isl_union_pw_qpolynomial_cow(upwqp);
4782 if (!upwqp)
4783 return NULL;
4784
4785 if (isl_hash_table_foreach(upwqp->dim->ctx, &upwqp->table,
4786 &poly_entry, &sign) < 0)
4787 goto error;
4788
4789 return upwqp;
4790 error:
4791 isl_union_pw_qpolynomial_free(upwqp);
4792 return NULL;
4793 }
4794
4795 __isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
4796 __isl_take isl_qpolynomial *qp)
4797 {
4798 int i, k;
4799 isl_space *dim;
4800 isl_vec *aff = NULL;
4801 isl_basic_map *bmap = NULL;
4802 unsigned pos;
4803 unsigned n_div;
4804
4805 if (!qp)
4806 return NULL;
4807 if (!isl_upoly_is_affine(qp->upoly))
4808 isl_die(qp->dim->ctx, isl_error_invalid,
4809 "input quasi-polynomial not affine", goto error);
4810 aff = isl_qpolynomial_extract_affine(qp);
4811 if (!aff)
4812 goto error;
4813 dim = isl_qpolynomial_get_space(qp);
4814 pos = 1 + isl_space_offset(dim, isl_dim_out);
4815 n_div = qp->div->n_row;
4816 bmap = isl_basic_map_alloc_space(dim, n_div, 1, 2 * n_div);
4817
4818 for (i = 0; i < n_div; ++i) {
4819 k = isl_basic_map_alloc_div(bmap);
4820 if (k < 0)
4821 goto error;
4822 isl_seq_cpy(bmap->div[k], qp->div->row[i], qp->div->n_col);
4823 isl_int_set_si(bmap->div[k][qp->div->n_col], 0);
4824 if (isl_basic_map_add_div_constraints(bmap, k) < 0)
4825 goto error;
4826 }
4827 k = isl_basic_map_alloc_equality(bmap);
4828 if (k < 0)
4829 goto error;
4830 isl_int_neg(bmap->eq[k][pos], aff->el[0]);
4831 isl_seq_cpy(bmap->eq[k], aff->el + 1, pos);
4832 isl_seq_cpy(bmap->eq[k] + pos + 1, aff->el + 1 + pos, n_div);
4833
4834 isl_vec_free(aff);
4835 isl_qpolynomial_free(qp);
4836 bmap = isl_basic_map_finalize(bmap);
4837 return bmap;
4838 error:
4839 isl_vec_free(aff);
4840 isl_qpolynomial_free(qp);
4841 isl_basic_map_free(bmap);
4842 return NULL;
4843 }
Integer set library