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