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