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