?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
isl_ast_build_expr.c: isl_ast_expr_add_term: allow stealing from constant term
[isl.git] / isl_polynomial.c
blob6b73bb397d098ada7d8948d5250624260d5c495c
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 #define ALIGN_DOMAIN
2812
2813 #include <isl_union_templ.c>
2814
2815 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial *pwqp)
2816 {
2817 if (!pwqp)
2818 return -1;
2819
2820 if (pwqp->n != -1)
2821 return 0;
2822
2823 if (!isl_set_plain_is_universe(pwqp->p[0].set))
2824 return 0;
2825
2826 return isl_qpolynomial_is_one(pwqp->p[0].qp);
2827 }
2828
2829 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add(
2830 __isl_take isl_pw_qpolynomial *pwqp1,
2831 __isl_take isl_pw_qpolynomial *pwqp2)
2832 {
2833 return isl_pw_qpolynomial_union_add_(pwqp1, pwqp2);
2834 }
2835
2836 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_mul(
2837 __isl_take isl_pw_qpolynomial *pwqp1,
2838 __isl_take isl_pw_qpolynomial *pwqp2)
2839 {
2840 int i, j, n;
2841 struct isl_pw_qpolynomial *res;
2842
2843 if (!pwqp1 || !pwqp2)
2844 goto error;
2845
2846 isl_assert(pwqp1->dim->ctx, isl_space_is_equal(pwqp1->dim, pwqp2->dim),
2847 goto error);
2848
2849 if (isl_pw_qpolynomial_is_zero(pwqp1)) {
2850 isl_pw_qpolynomial_free(pwqp2);
2851 return pwqp1;
2852 }
2853
2854 if (isl_pw_qpolynomial_is_zero(pwqp2)) {
2855 isl_pw_qpolynomial_free(pwqp1);
2856 return pwqp2;
2857 }
2858
2859 if (isl_pw_qpolynomial_is_one(pwqp1)) {
2860 isl_pw_qpolynomial_free(pwqp1);
2861 return pwqp2;
2862 }
2863
2864 if (isl_pw_qpolynomial_is_one(pwqp2)) {
2865 isl_pw_qpolynomial_free(pwqp2);
2866 return pwqp1;
2867 }
2868
2869 n = pwqp1->n * pwqp2->n;
2870 res = isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1->dim), n);
2871
2872 for (i = 0; i < pwqp1->n; ++i) {
2873 for (j = 0; j < pwqp2->n; ++j) {
2874 struct isl_set *common;
2875 struct isl_qpolynomial *prod;
2876 common = isl_set_intersect(isl_set_copy(pwqp1->p[i].set),
2877 isl_set_copy(pwqp2->p[j].set));
2878 if (isl_set_plain_is_empty(common)) {
2879 isl_set_free(common);
2880 continue;
2881 }
2882
2883 prod = isl_qpolynomial_mul(
2884 isl_qpolynomial_copy(pwqp1->p[i].qp),
2885 isl_qpolynomial_copy(pwqp2->p[j].qp));
2886
2887 res = isl_pw_qpolynomial_add_piece(res, common, prod);
2888 }
2889 }
2890
2891 isl_pw_qpolynomial_free(pwqp1);
2892 isl_pw_qpolynomial_free(pwqp2);
2893
2894 return res;
2895 error:
2896 isl_pw_qpolynomial_free(pwqp1);
2897 isl_pw_qpolynomial_free(pwqp2);
2898 return NULL;
2899 }
2900
2901 __isl_give isl_val *isl_upoly_eval(__isl_take struct isl_upoly *up,
2902 __isl_take isl_vec *vec)
2903 {
2904 int i;
2905 struct isl_upoly_rec *rec;
2906 isl_val *res;
2907 isl_val *base;
2908
2909 if (isl_upoly_is_cst(up)) {
2910 isl_vec_free(vec);
2911 res = isl_upoly_get_constant_val(up);
2912 isl_upoly_free(up);
2913 return res;
2914 }
2915
2916 rec = isl_upoly_as_rec(up);
2917 if (!rec)
2918 goto error;
2919
2920 isl_assert(up->ctx, rec->n >= 1, goto error);
2921
2922 base = isl_val_rat_from_isl_int(up->ctx,
2923 vec->el[1 + up->var], vec->el[0]);
2924
2925 res = isl_upoly_eval(isl_upoly_copy(rec->p[rec->n - 1]),
2926 isl_vec_copy(vec));
2927
2928 for (i = rec->n - 2; i >= 0; --i) {
2929 res = isl_val_mul(res, isl_val_copy(base));
2930 res = isl_val_add(res,
2931 isl_upoly_eval(isl_upoly_copy(rec->p[i]),
2932 isl_vec_copy(vec)));
2933 }
2934
2935 isl_val_free(base);
2936 isl_upoly_free(up);
2937 isl_vec_free(vec);
2938 return res;
2939 error:
2940 isl_upoly_free(up);
2941 isl_vec_free(vec);
2942 return NULL;
2943 }
2944
2945 __isl_give isl_val *isl_qpolynomial_eval(__isl_take isl_qpolynomial *qp,
2946 __isl_take isl_point *pnt)
2947 {
2948 isl_vec *ext;
2949 isl_val *v;
2950
2951 if (!qp || !pnt)
2952 goto error;
2953 isl_assert(pnt->dim->ctx, isl_space_is_equal(pnt->dim, qp->dim), goto error);
2954
2955 if (qp->div->n_row == 0)
2956 ext = isl_vec_copy(pnt->vec);
2957 else {
2958 int i;
2959 unsigned dim = isl_space_dim(qp->dim, isl_dim_all);
2960 ext = isl_vec_alloc(qp->dim->ctx, 1 + dim + qp->div->n_row);
2961 if (!ext)
2962 goto error;
2963
2964 isl_seq_cpy(ext->el, pnt->vec->el, pnt->vec->size);
2965 for (i = 0; i < qp->div->n_row; ++i) {
2966 isl_seq_inner_product(qp->div->row[i] + 1, ext->el,
2967 1 + dim + i, &ext->el[1+dim+i]);
2968 isl_int_fdiv_q(ext->el[1+dim+i], ext->el[1+dim+i],
2969 qp->div->row[i][0]);
2970 }
2971 }
2972
2973 v = isl_upoly_eval(isl_upoly_copy(qp->upoly), ext);
2974
2975 isl_qpolynomial_free(qp);
2976 isl_point_free(pnt);
2977
2978 return v;
2979 error:
2980 isl_qpolynomial_free(qp);
2981 isl_point_free(pnt);
2982 return NULL;
2983 }
2984
2985 int isl_upoly_cmp(__isl_keep struct isl_upoly_cst *cst1,
2986 __isl_keep struct isl_upoly_cst *cst2)
2987 {
2988 int cmp;
2989 isl_int t;
2990 isl_int_init(t);
2991 isl_int_mul(t, cst1->n, cst2->d);
2992 isl_int_submul(t, cst2->n, cst1->d);
2993 cmp = isl_int_sgn(t);
2994 isl_int_clear(t);
2995 return cmp;
2996 }
2997
2998 __isl_give isl_qpolynomial *isl_qpolynomial_insert_dims(
2999 __isl_take isl_qpolynomial *qp, enum isl_dim_type type,
3000 unsigned first, unsigned n)
3001 {
3002 unsigned total;
3003 unsigned g_pos;
3004 int *exp;
3005
3006 if (!qp)
3007 return NULL;
3008 if (type == isl_dim_out)
3009 isl_die(qp->div->ctx, isl_error_invalid,
3010 "cannot insert output/set dimensions",
3011 goto error);
3012 if (type == isl_dim_in)
3013 type = isl_dim_set;
3014 if (n == 0 && !isl_space_is_named_or_nested(qp->dim, type))
3015 return qp;
3016
3017 qp = isl_qpolynomial_cow(qp);
3018 if (!qp)
3019 return NULL;
3020
3021 isl_assert(qp->div->ctx, first <= isl_space_dim(qp->dim, type),
3022 goto error);
3023
3024 g_pos = pos(qp->dim, type) + first;
3025
3026 qp->div = isl_mat_insert_zero_cols(qp->div, 2 + g_pos, n);
3027 if (!qp->div)
3028 goto error;
3029
3030 total = qp->div->n_col - 2;
3031 if (total > g_pos) {
3032 int i;
3033 exp = isl_alloc_array(qp->div->ctx, int, total - g_pos);
3034 if (!exp)
3035 goto error;
3036 for (i = 0; i < total - g_pos; ++i)
3037 exp[i] = i + n;
3038 qp->upoly = expand(qp->upoly, exp, g_pos);
3039 free(exp);
3040 if (!qp->upoly)
3041 goto error;
3042 }
3043
3044 qp->dim = isl_space_insert_dims(qp->dim, type, first, n);
3045 if (!qp->dim)
3046 goto error;
3047
3048 return qp;
3049 error:
3050 isl_qpolynomial_free(qp);
3051 return NULL;
3052 }
3053
3054 __isl_give isl_qpolynomial *isl_qpolynomial_add_dims(
3055 __isl_take isl_qpolynomial *qp, enum isl_dim_type type, unsigned n)
3056 {
3057 unsigned pos;
3058
3059 pos = isl_qpolynomial_dim(qp, type);
3060
3061 return isl_qpolynomial_insert_dims(qp, type, pos, n);
3062 }
3063
3064 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_add_dims(
3065 __isl_take isl_pw_qpolynomial *pwqp,
3066 enum isl_dim_type type, unsigned n)
3067 {
3068 unsigned pos;
3069
3070 pos = isl_pw_qpolynomial_dim(pwqp, type);
3071
3072 return isl_pw_qpolynomial_insert_dims(pwqp, type, pos, n);
3073 }
3074
3075 static int *reordering_move(isl_ctx *ctx,
3076 unsigned len, unsigned dst, unsigned src, unsigned n)
3077 {
3078 int i;
3079 int *reordering;
3080
3081 reordering = isl_alloc_array(ctx, int, len);
3082 if (!reordering)
3083 return NULL;
3084
3085 if (dst <= src) {
3086 for (i = 0; i < dst; ++i)
3087 reordering[i] = i;
3088 for (i = 0; i < n; ++i)
3089 reordering[src + i] = dst + i;
3090 for (i = 0; i < src - dst; ++i)
3091 reordering[dst + i] = dst + n + i;
3092 for (i = 0; i < len - src - n; ++i)
3093 reordering[src + n + i] = src + n + i;
3094 } else {
3095 for (i = 0; i < src; ++i)
3096 reordering[i] = i;
3097 for (i = 0; i < n; ++i)
3098 reordering[src + i] = dst + i;
3099 for (i = 0; i < dst - src; ++i)
3100 reordering[src + n + i] = src + i;
3101 for (i = 0; i < len - dst - n; ++i)
3102 reordering[dst + n + i] = dst + n + i;
3103 }
3104
3105 return reordering;
3106 }
3107
3108 __isl_give isl_qpolynomial *isl_qpolynomial_move_dims(
3109 __isl_take isl_qpolynomial *qp,
3110 enum isl_dim_type dst_type, unsigned dst_pos,
3111 enum isl_dim_type src_type, unsigned src_pos, unsigned n)
3112 {
3113 unsigned g_dst_pos;
3114 unsigned g_src_pos;
3115 int *reordering;
3116
3117 if (n == 0)
3118 return qp;
3119
3120 qp = isl_qpolynomial_cow(qp);
3121 if (!qp)
3122 return NULL;
3123
3124 if (dst_type == isl_dim_out || src_type == isl_dim_out)
3125 isl_die(qp->dim->ctx, isl_error_invalid,
3126 "cannot move output/set dimension",
3127 goto error);
3128 if (dst_type == isl_dim_in)
3129 dst_type = isl_dim_set;
3130 if (src_type == isl_dim_in)
3131 src_type = isl_dim_set;
3132
3133 isl_assert(qp->dim->ctx, src_pos + n <= isl_space_dim(qp->dim, src_type),
3134 goto error);
3135
3136 g_dst_pos = pos(qp->dim, dst_type) + dst_pos;
3137 g_src_pos = pos(qp->dim, src_type) + src_pos;
3138 if (dst_type > src_type)
3139 g_dst_pos -= n;
3140
3141 qp->div = isl_mat_move_cols(qp->div, 2 + g_dst_pos, 2 + g_src_pos, n);
3142 if (!qp->div)
3143 goto error;
3144 qp = sort_divs(qp);
3145 if (!qp)
3146 goto error;
3147
3148 reordering = reordering_move(qp->dim->ctx,
3149 qp->div->n_col - 2, g_dst_pos, g_src_pos, n);
3150 if (!reordering)
3151 goto error;
3152
3153 qp->upoly = reorder(qp->upoly, reordering);
3154 free(reordering);
3155 if (!qp->upoly)
3156 goto error;
3157
3158 qp->dim = isl_space_move_dims(qp->dim, dst_type, dst_pos, src_type, src_pos, n);
3159 if (!qp->dim)
3160 goto error;
3161
3162 return qp;
3163 error:
3164 isl_qpolynomial_free(qp);
3165 return NULL;
3166 }
3167
3168 __isl_give isl_qpolynomial *isl_qpolynomial_from_affine(__isl_take isl_space *dim,
3169 isl_int *f, isl_int denom)
3170 {
3171 struct isl_upoly *up;
3172
3173 dim = isl_space_domain(dim);
3174 if (!dim)
3175 return NULL;
3176
3177 up = isl_upoly_from_affine(dim->ctx, f, denom,
3178 1 + isl_space_dim(dim, isl_dim_all));
3179
3180 return isl_qpolynomial_alloc(dim, 0, up);
3181 }
3182
3183 __isl_give isl_qpolynomial *isl_qpolynomial_from_aff(__isl_take isl_aff *aff)
3184 {
3185 isl_ctx *ctx;
3186 struct isl_upoly *up;
3187 isl_qpolynomial *qp;
3188
3189 if (!aff)
3190 return NULL;
3191
3192 ctx = isl_aff_get_ctx(aff);
3193 up = isl_upoly_from_affine(ctx, aff->v->el + 1, aff->v->el[0],
3194 aff->v->size - 1);
3195
3196 qp = isl_qpolynomial_alloc(isl_aff_get_domain_space(aff),
3197 aff->ls->div->n_row, up);
3198 if (!qp)
3199 goto error;
3200
3201 isl_mat_free(qp->div);
3202 qp->div = isl_mat_copy(aff->ls->div);
3203 qp->div = isl_mat_cow(qp->div);
3204 if (!qp->div)
3205 goto error;
3206
3207 isl_aff_free(aff);
3208 qp = reduce_divs(qp);
3209 qp = remove_redundant_divs(qp);
3210 return qp;
3211 error:
3212 isl_aff_free(aff);
3213 return isl_qpolynomial_free(qp);
3214 }
3215
3216 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_from_pw_aff(
3217 __isl_take isl_pw_aff *pwaff)
3218 {
3219 int i;
3220 isl_pw_qpolynomial *pwqp;
3221
3222 if (!pwaff)
3223 return NULL;
3224
3225 pwqp = isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff),
3226 pwaff->n);
3227
3228 for (i = 0; i < pwaff->n; ++i) {
3229 isl_set *dom;
3230 isl_qpolynomial *qp;
3231
3232 dom = isl_set_copy(pwaff->p[i].set);
3233 qp = isl_qpolynomial_from_aff(isl_aff_copy(pwaff->p[i].aff));
3234 pwqp = isl_pw_qpolynomial_add_piece(pwqp, dom, qp);
3235 }
3236
3237 isl_pw_aff_free(pwaff);
3238 return pwqp;
3239 }
3240
3241 __isl_give isl_qpolynomial *isl_qpolynomial_from_constraint(
3242 __isl_take isl_constraint *c, enum isl_dim_type type, unsigned pos)
3243 {
3244 isl_aff *aff;
3245
3246 aff = isl_constraint_get_bound(c, type, pos);
3247 isl_constraint_free(c);
3248 return isl_qpolynomial_from_aff(aff);
3249 }
3250
3251 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
3252 * in "qp" by subs[i].
3253 */
3254 __isl_give isl_qpolynomial *isl_qpolynomial_substitute(
3255 __isl_take isl_qpolynomial *qp,
3256 enum isl_dim_type type, unsigned first, unsigned n,
3257 __isl_keep isl_qpolynomial **subs)
3258 {
3259 int i;
3260 struct isl_upoly **ups;
3261
3262 if (n == 0)
3263 return qp;
3264
3265 qp = isl_qpolynomial_cow(qp);
3266 if (!qp)
3267 return NULL;
3268
3269 if (type == isl_dim_out)
3270 isl_die(qp->dim->ctx, isl_error_invalid,
3271 "cannot substitute output/set dimension",
3272 goto error);
3273 if (type == isl_dim_in)
3274 type = isl_dim_set;
3275
3276 for (i = 0; i < n; ++i)
3277 if (!subs[i])
3278 goto error;
3279
3280 isl_assert(qp->dim->ctx, first + n <= isl_space_dim(qp->dim, type),
3281 goto error);
3282
3283 for (i = 0; i < n; ++i)
3284 isl_assert(qp->dim->ctx, isl_space_is_equal(qp->dim, subs[i]->dim),
3285 goto error);
3286
3287 isl_assert(qp->dim->ctx, qp->div->n_row == 0, goto error);
3288 for (i = 0; i < n; ++i)
3289 isl_assert(qp->dim->ctx, subs[i]->div->n_row == 0, goto error);
3290
3291 first += pos(qp->dim, type);
3292
3293 ups = isl_alloc_array(qp->dim->ctx, struct isl_upoly *, n);
3294 if (!ups)
3295 goto error;
3296 for (i = 0; i < n; ++i)
3297 ups[i] = subs[i]->upoly;
3298
3299 qp->upoly = isl_upoly_subs(qp->upoly, first, n, ups);
3300
3301 free(ups);
3302
3303 if (!qp->upoly)
3304 goto error;
3305
3306 return qp;
3307 error:
3308 isl_qpolynomial_free(qp);
3309 return NULL;
3310 }
3311
3312 /* Extend "bset" with extra set dimensions for each integer division
3313 * in "qp" and then call "fn" with the extended bset and the polynomial
3314 * that results from replacing each of the integer divisions by the
3315 * corresponding extra set dimension.
3316 */
3317 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial *qp,
3318 __isl_keep isl_basic_set *bset,
3319 int (*fn)(__isl_take isl_basic_set *bset,
3320 __isl_take isl_qpolynomial *poly, void *user), void *user)
3321 {
3322 isl_space *dim;
3323 isl_mat *div;
3324 isl_qpolynomial *poly;
3325
3326 if (!qp || !bset)
3327 goto error;
3328 if (qp->div->n_row == 0)
3329 return fn(isl_basic_set_copy(bset), isl_qpolynomial_copy(qp),
3330 user);
3331
3332 div = isl_mat_copy(qp->div);
3333 dim = isl_space_copy(qp->dim);
3334 dim = isl_space_add_dims(dim, isl_dim_set, qp->div->n_row);
3335 poly = isl_qpolynomial_alloc(dim, 0, isl_upoly_copy(qp->upoly));
3336 bset = isl_basic_set_copy(bset);
3337 bset = isl_basic_set_add_dims(bset, isl_dim_set, qp->div->n_row);
3338 bset = add_div_constraints(bset, div);
3339
3340 return fn(bset, poly, user);
3341 error:
3342 return -1;
3343 }
3344
3345 /* Return total degree in variables first (inclusive) up to last (exclusive).
3346 */
3347 int isl_upoly_degree(__isl_keep struct isl_upoly *up, int first, int last)
3348 {
3349 int deg = -1;
3350 int i;
3351 struct isl_upoly_rec *rec;
3352
3353 if (!up)
3354 return -2;
3355 if (isl_upoly_is_zero(up))
3356 return -1;
3357 if (isl_upoly_is_cst(up) || up->var < first)
3358 return 0;
3359
3360 rec = isl_upoly_as_rec(up);
3361 if (!rec)
3362 return -2;
3363
3364 for (i = 0; i < rec->n; ++i) {
3365 int d;
3366
3367 if (isl_upoly_is_zero(rec->p[i]))
3368 continue;
3369 d = isl_upoly_degree(rec->p[i], first, last);
3370 if (up->var < last)
3371 d += i;
3372 if (d > deg)
3373 deg = d;
3374 }
3375
3376 return deg;
3377 }
3378
3379 /* Return total degree in set variables.
3380 */
3381 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial *poly)
3382 {
3383 unsigned ovar;
3384 unsigned nvar;
3385
3386 if (!poly)
3387 return -2;
3388
3389 ovar = isl_space_offset(poly->dim, isl_dim_set);
3390 nvar = isl_space_dim(poly->dim, isl_dim_set);
3391 return isl_upoly_degree(poly->upoly, ovar, ovar + nvar);
3392 }
3393
3394 __isl_give struct isl_upoly *isl_upoly_coeff(__isl_keep struct isl_upoly *up,
3395 unsigned pos, int deg)
3396 {
3397 int i;
3398 struct isl_upoly_rec *rec;
3399
3400 if (!up)
3401 return NULL;
3402
3403 if (isl_upoly_is_cst(up) || up->var < pos) {
3404 if (deg == 0)
3405 return isl_upoly_copy(up);
3406 else
3407 return isl_upoly_zero(up->ctx);
3408 }
3409
3410 rec = isl_upoly_as_rec(up);
3411 if (!rec)
3412 return NULL;
3413
3414 if (up->var == pos) {
3415 if (deg < rec->n)
3416 return isl_upoly_copy(rec->p[deg]);
3417 else
3418 return isl_upoly_zero(up->ctx);
3419 }
3420
3421 up = isl_upoly_copy(up);
3422 up = isl_upoly_cow(up);
3423 rec = isl_upoly_as_rec(up);
3424 if (!rec)
3425 goto error;
3426
3427 for (i = 0; i < rec->n; ++i) {
3428 struct isl_upoly *t;
3429 t = isl_upoly_coeff(rec->p[i], pos, deg);
3430 if (!t)
3431 goto error;
3432 isl_upoly_free(rec->p[i]);
3433 rec->p[i] = t;
3434 }
3435
3436 return up;
3437 error:
3438 isl_upoly_free(up);
3439 return NULL;
3440 }
3441
3442 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
3443 */
3444 __isl_give isl_qpolynomial *isl_qpolynomial_coeff(
3445 __isl_keep isl_qpolynomial *qp,
3446 enum isl_dim_type type, unsigned t_pos, int deg)
3447 {
3448 unsigned g_pos;
3449 struct isl_upoly *up;
3450 isl_qpolynomial *c;
3451
3452 if (!qp)
3453 return NULL;
3454
3455 if (type == isl_dim_out)
3456 isl_die(qp->div->ctx, isl_error_invalid,
3457 "output/set dimension does not have a coefficient",
3458 return NULL);
3459 if (type == isl_dim_in)
3460 type = isl_dim_set;
3461
3462 isl_assert(qp->div->ctx, t_pos < isl_space_dim(qp->dim, type),
3463 return NULL);
3464
3465 g_pos = pos(qp->dim, type) + t_pos;
3466 up = isl_upoly_coeff(qp->upoly, g_pos, deg);
3467
3468 c = isl_qpolynomial_alloc(isl_space_copy(qp->dim), qp->div->n_row, up);
3469 if (!c)
3470 return NULL;
3471 isl_mat_free(c->div);
3472 c->div = isl_mat_copy(qp->div);
3473 if (!c->div)
3474 goto error;
3475 return c;
3476 error:
3477 isl_qpolynomial_free(c);
3478 return NULL;
3479 }
3480
3481 /* Homogenize the polynomial in the variables first (inclusive) up to
3482 * last (exclusive) by inserting powers of variable first.
3483 * Variable first is assumed not to appear in the input.
3484 */
3485 __isl_give struct isl_upoly *isl_upoly_homogenize(
3486 __isl_take struct isl_upoly *up, int deg, int target,
3487 int first, int last)
3488 {
3489 int i;
3490 struct isl_upoly_rec *rec;
3491
3492 if (!up)
3493 return NULL;
3494 if (isl_upoly_is_zero(up))
3495 return up;
3496 if (deg == target)
3497 return up;
3498 if (isl_upoly_is_cst(up) || up->var < first) {
3499 struct isl_upoly *hom;
3500
3501 hom = isl_upoly_var_pow(up->ctx, first, target - deg);
3502 if (!hom)
3503 goto error;
3504 rec = isl_upoly_as_rec(hom);
3505 rec->p[target - deg] = isl_upoly_mul(rec->p[target - deg], up);
3506
3507 return hom;
3508 }
3509
3510 up = isl_upoly_cow(up);
3511 rec = isl_upoly_as_rec(up);
3512 if (!rec)
3513 goto error;
3514
3515 for (i = 0; i < rec->n; ++i) {
3516 if (isl_upoly_is_zero(rec->p[i]))
3517 continue;
3518 rec->p[i] = isl_upoly_homogenize(rec->p[i],
3519 up->var < last ? deg + i : i, target,
3520 first, last);
3521 if (!rec->p[i])
3522 goto error;
3523 }
3524
3525 return up;
3526 error:
3527 isl_upoly_free(up);
3528 return NULL;
3529 }
3530
3531 /* Homogenize the polynomial in the set variables by introducing
3532 * powers of an extra set variable at position 0.
3533 */
3534 __isl_give isl_qpolynomial *isl_qpolynomial_homogenize(
3535 __isl_take isl_qpolynomial *poly)
3536 {
3537 unsigned ovar;
3538 unsigned nvar;
3539 int deg = isl_qpolynomial_degree(poly);
3540
3541 if (deg < -1)
3542 goto error;
3543
3544 poly = isl_qpolynomial_insert_dims(poly, isl_dim_in, 0, 1);
3545 poly = isl_qpolynomial_cow(poly);
3546 if (!poly)
3547 goto error;
3548
3549 ovar = isl_space_offset(poly->dim, isl_dim_set);
3550 nvar = isl_space_dim(poly->dim, isl_dim_set);
3551 poly->upoly = isl_upoly_homogenize(poly->upoly, 0, deg,
3552 ovar, ovar + nvar);
3553 if (!poly->upoly)
3554 goto error;
3555
3556 return poly;
3557 error:
3558 isl_qpolynomial_free(poly);
3559 return NULL;
3560 }
3561
3562 __isl_give isl_term *isl_term_alloc(__isl_take isl_space *dim,
3563 __isl_take isl_mat *div)
3564 {
3565 isl_term *term;
3566 int n;
3567
3568 if (!dim || !div)
3569 goto error;
3570
3571 n = isl_space_dim(dim, isl_dim_all) + div->n_row;
3572
3573 term = isl_calloc(dim->ctx, struct isl_term,
3574 sizeof(struct isl_term) + (n - 1) * sizeof(int));
3575 if (!term)
3576 goto error;
3577
3578 term->ref = 1;
3579 term->dim = dim;
3580 term->div = div;
3581 isl_int_init(term->n);
3582 isl_int_init(term->d);
3583
3584 return term;
3585 error:
3586 isl_space_free(dim);
3587 isl_mat_free(div);
3588 return NULL;
3589 }
3590
3591 __isl_give isl_term *isl_term_copy(__isl_keep isl_term *term)
3592 {
3593 if (!term)
3594 return NULL;
3595
3596 term->ref++;
3597 return term;
3598 }
3599
3600 __isl_give isl_term *isl_term_dup(__isl_keep isl_term *term)
3601 {
3602 int i;
3603 isl_term *dup;
3604 unsigned total;
3605
3606 if (!term)
3607 return NULL;
3608
3609 total = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3610
3611 dup = isl_term_alloc(isl_space_copy(term->dim), isl_mat_copy(term->div));
3612 if (!dup)
3613 return NULL;
3614
3615 isl_int_set(dup->n, term->n);
3616 isl_int_set(dup->d, term->d);
3617
3618 for (i = 0; i < total; ++i)
3619 dup->pow[i] = term->pow[i];
3620
3621 return dup;
3622 }
3623
3624 __isl_give isl_term *isl_term_cow(__isl_take isl_term *term)
3625 {
3626 if (!term)
3627 return NULL;
3628
3629 if (term->ref == 1)
3630 return term;
3631 term->ref--;
3632 return isl_term_dup(term);
3633 }
3634
3635 void isl_term_free(__isl_take isl_term *term)
3636 {
3637 if (!term)
3638 return;
3639
3640 if (--term->ref > 0)
3641 return;
3642
3643 isl_space_free(term->dim);
3644 isl_mat_free(term->div);
3645 isl_int_clear(term->n);
3646 isl_int_clear(term->d);
3647 free(term);
3648 }
3649
3650 unsigned isl_term_dim(__isl_keep isl_term *term, enum isl_dim_type type)
3651 {
3652 if (!term)
3653 return 0;
3654
3655 switch (type) {
3656 case isl_dim_param:
3657 case isl_dim_in:
3658 case isl_dim_out: return isl_space_dim(term->dim, type);
3659 case isl_dim_div: return term->div->n_row;
3660 case isl_dim_all: return isl_space_dim(term->dim, isl_dim_all) +
3661 term->div->n_row;
3662 default: return 0;
3663 }
3664 }
3665
3666 isl_ctx *isl_term_get_ctx(__isl_keep isl_term *term)
3667 {
3668 return term ? term->dim->ctx : NULL;
3669 }
3670
3671 void isl_term_get_num(__isl_keep isl_term *term, isl_int *n)
3672 {
3673 if (!term)
3674 return;
3675 isl_int_set(*n, term->n);
3676 }
3677
3678 void isl_term_get_den(__isl_keep isl_term *term, isl_int *d)
3679 {
3680 if (!term)
3681 return;
3682 isl_int_set(*d, term->d);
3683 }
3684
3685 /* Return the coefficient of the term "term".
3686 */
3687 __isl_give isl_val *isl_term_get_coefficient_val(__isl_keep isl_term *term)
3688 {
3689 if (!term)
3690 return NULL;
3691
3692 return isl_val_rat_from_isl_int(isl_term_get_ctx(term),
3693 term->n, term->d);
3694 }
3695
3696 int isl_term_get_exp(__isl_keep isl_term *term,
3697 enum isl_dim_type type, unsigned pos)
3698 {
3699 if (!term)
3700 return -1;
3701
3702 isl_assert(term->dim->ctx, pos < isl_term_dim(term, type), return -1);
3703
3704 if (type >= isl_dim_set)
3705 pos += isl_space_dim(term->dim, isl_dim_param);
3706 if (type >= isl_dim_div)
3707 pos += isl_space_dim(term->dim, isl_dim_set);
3708
3709 return term->pow[pos];
3710 }
3711
3712 __isl_give isl_aff *isl_term_get_div(__isl_keep isl_term *term, unsigned pos)
3713 {
3714 isl_local_space *ls;
3715 isl_aff *aff;
3716
3717 if (!term)
3718 return NULL;
3719
3720 isl_assert(term->dim->ctx, pos < isl_term_dim(term, isl_dim_div),
3721 return NULL);
3722
3723 ls = isl_local_space_alloc_div(isl_space_copy(term->dim),
3724 isl_mat_copy(term->div));
3725 aff = isl_aff_alloc(ls);
3726 if (!aff)
3727 return NULL;
3728
3729 isl_seq_cpy(aff->v->el, term->div->row[pos], aff->v->size);
3730
3731 aff = isl_aff_normalize(aff);
3732
3733 return aff;
3734 }
3735
3736 __isl_give isl_term *isl_upoly_foreach_term(__isl_keep struct isl_upoly *up,
3737 int (*fn)(__isl_take isl_term *term, void *user),
3738 __isl_take isl_term *term, void *user)
3739 {
3740 int i;
3741 struct isl_upoly_rec *rec;
3742
3743 if (!up || !term)
3744 goto error;
3745
3746 if (isl_upoly_is_zero(up))
3747 return term;
3748
3749 isl_assert(up->ctx, !isl_upoly_is_nan(up), goto error);
3750 isl_assert(up->ctx, !isl_upoly_is_infty(up), goto error);
3751 isl_assert(up->ctx, !isl_upoly_is_neginfty(up), goto error);
3752
3753 if (isl_upoly_is_cst(up)) {
3754 struct isl_upoly_cst *cst;
3755 cst = isl_upoly_as_cst(up);
3756 if (!cst)
3757 goto error;
3758 term = isl_term_cow(term);
3759 if (!term)
3760 goto error;
3761 isl_int_set(term->n, cst->n);
3762 isl_int_set(term->d, cst->d);
3763 if (fn(isl_term_copy(term), user) < 0)
3764 goto error;
3765 return term;
3766 }
3767
3768 rec = isl_upoly_as_rec(up);
3769 if (!rec)
3770 goto error;
3771
3772 for (i = 0; i < rec->n; ++i) {
3773 term = isl_term_cow(term);
3774 if (!term)
3775 goto error;
3776 term->pow[up->var] = i;
3777 term = isl_upoly_foreach_term(rec->p[i], fn, term, user);
3778 if (!term)
3779 goto error;
3780 }
3781 term->pow[up->var] = 0;
3782
3783 return term;
3784 error:
3785 isl_term_free(term);
3786 return NULL;
3787 }
3788
3789 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial *qp,
3790 int (*fn)(__isl_take isl_term *term, void *user), void *user)
3791 {
3792 isl_term *term;
3793
3794 if (!qp)
3795 return -1;
3796
3797 term = isl_term_alloc(isl_space_copy(qp->dim), isl_mat_copy(qp->div));
3798 if (!term)
3799 return -1;
3800
3801 term = isl_upoly_foreach_term(qp->upoly, fn, term, user);
3802
3803 isl_term_free(term);
3804
3805 return term ? 0 : -1;
3806 }
3807
3808 __isl_give isl_qpolynomial *isl_qpolynomial_from_term(__isl_take isl_term *term)
3809 {
3810 struct isl_upoly *up;
3811 isl_qpolynomial *qp;
3812 int i, n;
3813
3814 if (!term)
3815 return NULL;
3816
3817 n = isl_space_dim(term->dim, isl_dim_all) + term->div->n_row;
3818
3819 up = isl_upoly_rat_cst(term->dim->ctx, term->n, term->d);
3820 for (i = 0; i < n; ++i) {
3821 if (!term->pow[i])
3822 continue;
3823 up = isl_upoly_mul(up,
3824 isl_upoly_var_pow(term->dim->ctx, i, term->pow[i]));
3825 }
3826
3827 qp = isl_qpolynomial_alloc(isl_space_copy(term->dim), term->div->n_row, up);
3828 if (!qp)
3829 goto error;
3830 isl_mat_free(qp->div);
3831 qp->div = isl_mat_copy(term->div);
3832 if (!qp->div)
3833 goto error;
3834
3835 isl_term_free(term);
3836 return qp;
3837 error:
3838 isl_qpolynomial_free(qp);
3839 isl_term_free(term);
3840 return NULL;
3841 }
3842
3843 __isl_give isl_qpolynomial *isl_qpolynomial_lift(__isl_take isl_qpolynomial *qp,
3844 __isl_take isl_space *dim)
3845 {
3846 int i;
3847 int extra;
3848 unsigned total;
3849
3850 if (!qp || !dim)
3851 goto error;
3852
3853 if (isl_space_is_equal(qp->dim, dim)) {
3854 isl_space_free(dim);
3855 return qp;
3856 }
3857
3858 qp = isl_qpolynomial_cow(qp);
3859 if (!qp)
3860 goto error;
3861
3862 extra = isl_space_dim(dim, isl_dim_set) -
3863 isl_space_dim(qp->dim, isl_dim_set);
3864 total = isl_space_dim(qp->dim, isl_dim_all);
3865 if (qp->div->n_row) {
3866 int *exp;
3867
3868 exp = isl_alloc_array(qp->div->ctx, int, qp->div->n_row);
3869 if (!exp)
3870 goto error;
3871 for (i = 0; i < qp->div->n_row; ++i)
3872 exp[i] = extra + i;
3873 qp->upoly = expand(qp->upoly, exp, total);
3874 free(exp);
3875 if (!qp->upoly)
3876 goto error;
3877 }
3878 qp->div = isl_mat_insert_cols(qp->div, 2 + total, extra);
3879 if (!qp->div)
3880 goto error;
3881 for (i = 0; i < qp->div->n_row; ++i)
3882 isl_seq_clr(qp->div->row[i] + 2 + total, extra);
3883
3884 isl_space_free(qp->dim);
3885 qp->dim = dim;
3886
3887 return qp;
3888 error:
3889 isl_space_free(dim);
3890 isl_qpolynomial_free(qp);
3891 return NULL;
3892 }
3893
3894 /* For each parameter or variable that does not appear in qp,
3895 * first eliminate the variable from all constraints and then set it to zero.
3896 */
3897 static __isl_give isl_set *fix_inactive(__isl_take isl_set *set,
3898 __isl_keep isl_qpolynomial *qp)
3899 {
3900 int *active = NULL;
3901 int i;
3902 int d;
3903 unsigned nparam;
3904 unsigned nvar;
3905
3906 if (!set || !qp)
3907 goto error;
3908
3909 d = isl_space_dim(set->dim, isl_dim_all);
3910 active = isl_calloc_array(set->ctx, int, d);
3911 if (set_active(qp, active) < 0)
3912 goto error;
3913
3914 for (i = 0; i < d; ++i)
3915 if (!active[i])
3916 break;
3917
3918 if (i == d) {
3919 free(active);
3920 return set;
3921 }
3922
3923 nparam = isl_space_dim(set->dim, isl_dim_param);
3924 nvar = isl_space_dim(set->dim, isl_dim_set);
3925 for (i = 0; i < nparam; ++i) {
3926 if (active[i])
3927 continue;
3928 set = isl_set_eliminate(set, isl_dim_param, i, 1);
3929 set = isl_set_fix_si(set, isl_dim_param, i, 0);
3930 }
3931 for (i = 0; i < nvar; ++i) {
3932 if (active[nparam + i])
3933 continue;
3934 set = isl_set_eliminate(set, isl_dim_set, i, 1);
3935 set = isl_set_fix_si(set, isl_dim_set, i, 0);
3936 }
3937
3938 free(active);
3939
3940 return set;
3941 error:
3942 free(active);
3943 isl_set_free(set);
3944 return NULL;
3945 }
3946
3947 struct isl_opt_data {
3948 isl_qpolynomial *qp;
3949 int first;
3950 isl_val *opt;
3951 int max;
3952 };
3953
3954 static int opt_fn(__isl_take isl_point *pnt, void *user)
3955 {
3956 struct isl_opt_data *data = (struct isl_opt_data *)user;
3957 isl_val *val;
3958
3959 val = isl_qpolynomial_eval(isl_qpolynomial_copy(data->qp), pnt);
3960 if (data->first) {
3961 data->first = 0;
3962 data->opt = val;
3963 } else if (data->max) {
3964 data->opt = isl_val_max(data->opt, val);
3965 } else {
3966 data->opt = isl_val_min(data->opt, val);
3967 }
3968
3969 return 0;
3970 }
3971
3972 __isl_give isl_val *isl_qpolynomial_opt_on_domain(
3973 __isl_take isl_qpolynomial *qp, __isl_take isl_set *set, int max)
3974 {
3975 struct isl_opt_data data = { NULL, 1, NULL, max };
3976
3977 if (!set || !qp)
3978 goto error;
3979
3980 if (isl_upoly_is_cst(qp->upoly)) {
3981 isl_set_free(set);
3982 data.opt = isl_qpolynomial_get_constant_val(qp);
3983 isl_qpolynomial_free(qp);
3984 return data.opt;
3985 }
3986
3987 set = fix_inactive(set, qp);
3988
3989 data.qp = qp;
3990 if (isl_set_foreach_point(set, opt_fn, &data) < 0)
3991 goto error;
3992
3993 if (data.first)
3994 data.opt = isl_val_zero(isl_set_get_ctx(set));
3995
3996 isl_set_free(set);
3997 isl_qpolynomial_free(qp);
3998 return data.opt;
3999 error:
4000 isl_set_free(set);
4001 isl_qpolynomial_free(qp);
4002 isl_val_free(data.opt);
4003 return NULL;
4004 }
4005
4006 __isl_give isl_qpolynomial *isl_qpolynomial_morph_domain(
4007 __isl_take isl_qpolynomial *qp, __isl_take isl_morph *morph)
4008 {
4009 int i;
4010 int n_sub;
4011 isl_ctx *ctx;
4012 struct isl_upoly **subs;
4013 isl_mat *mat, *diag;
4014
4015 qp = isl_qpolynomial_cow(qp);
4016 if (!qp || !morph)
4017 goto error;
4018
4019 ctx = qp->dim->ctx;
4020 isl_assert(ctx, isl_space_is_equal(qp->dim, morph->dom->dim), goto error);
4021
4022 n_sub = morph->inv->n_row - 1;
4023 if (morph->inv->n_row != morph->inv->n_col)
4024 n_sub += qp->div->n_row;
4025 subs = isl_calloc_array(ctx, struct isl_upoly *, n_sub);
4026 if (n_sub && !subs)
4027 goto error;
4028
4029 for (i = 0; 1 + i < morph->inv->n_row; ++i)
4030 subs[i] = isl_upoly_from_affine(ctx, morph->inv->row[1 + i],
4031 morph->inv->row[0][0], morph->inv->n_col);
4032 if (morph->inv->n_row != morph->inv->n_col)
4033 for (i = 0; i < qp->div->n_row; ++i)
4034 subs[morph->inv->n_row - 1 + i] =
4035 isl_upoly_var_pow(ctx, morph->inv->n_col - 1 + i, 1);
4036
4037 qp->upoly = isl_upoly_subs(qp->upoly, 0, n_sub, subs);
4038
4039 for (i = 0; i < n_sub; ++i)
4040 isl_upoly_free(subs[i]);
4041 free(subs);
4042
4043 diag = isl_mat_diag(ctx, 1, morph->inv->row[0][0]);
4044 mat = isl_mat_diagonal(diag, isl_mat_copy(morph->inv));
4045 diag = isl_mat_diag(ctx, qp->div->n_row, morph->inv->row[0][0]);
4046 mat = isl_mat_diagonal(mat, diag);
4047 qp->div = isl_mat_product(qp->div, mat);
4048 isl_space_free(qp->dim);
4049 qp->dim = isl_space_copy(morph->ran->dim);
4050
4051 if (!qp->upoly || !qp->div || !qp->dim)
4052 goto error;
4053
4054 isl_morph_free(morph);
4055
4056 return qp;
4057 error:
4058 isl_qpolynomial_free(qp);
4059 isl_morph_free(morph);
4060 return NULL;
4061 }
4062
4063 static int neg_entry(void **entry, void *user)
4064 {
4065 isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
4066
4067 *pwqp = isl_pw_qpolynomial_neg(*pwqp);
4068
4069 return *pwqp ? 0 : -1;
4070 }
4071
4072 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_neg(
4073 __isl_take isl_union_pw_qpolynomial *upwqp)
4074 {
4075 upwqp = isl_union_pw_qpolynomial_cow(upwqp);
4076 if (!upwqp)
4077 return NULL;
4078
4079 if (isl_hash_table_foreach(upwqp->space->ctx, &upwqp->table,
4080 &neg_entry, NULL) < 0)
4081 goto error;
4082
4083 return upwqp;
4084 error:
4085 isl_union_pw_qpolynomial_free(upwqp);
4086 return NULL;
4087 }
4088
4089 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_mul(
4090 __isl_take isl_union_pw_qpolynomial *upwqp1,
4091 __isl_take isl_union_pw_qpolynomial *upwqp2)
4092 {
4093 return match_bin_op(upwqp1, upwqp2, &isl_pw_qpolynomial_mul);
4094 }
4095
4096 /* Reorder the columns of the given div definitions according to the
4097 * given reordering.
4098 */
4099 static __isl_give isl_mat *reorder_divs(__isl_take isl_mat *div,
4100 __isl_take isl_reordering *r)
4101 {
4102 int i, j;
4103 isl_mat *mat;
4104 int extra;
4105
4106 if (!div || !r)
4107 goto error;
4108
4109 extra = isl_space_dim(r->dim, isl_dim_all) + div->n_row - r->len;
4110 mat = isl_mat_alloc(div->ctx, div->n_row, div->n_col + extra);
4111 if (!mat)
4112 goto error;
4113
4114 for (i = 0; i < div->n_row; ++i) {
4115 isl_seq_cpy(mat->row[i], div->row[i], 2);
4116 isl_seq_clr(mat->row[i] + 2, mat->n_col - 2);
4117 for (j = 0; j < r->len; ++j)
4118 isl_int_set(mat->row[i][2 + r->pos[j]],
4119 div->row[i][2 + j]);
4120 }
4121
4122 isl_reordering_free(r);
4123 isl_mat_free(div);
4124 return mat;
4125 error:
4126 isl_reordering_free(r);
4127 isl_mat_free(div);
4128 return NULL;
4129 }
4130
4131 /* Reorder the dimension of "qp" according to the given reordering.
4132 */
4133 __isl_give isl_qpolynomial *isl_qpolynomial_realign_domain(
4134 __isl_take isl_qpolynomial *qp, __isl_take isl_reordering *r)
4135 {
4136 qp = isl_qpolynomial_cow(qp);
4137 if (!qp)
4138 goto error;
4139
4140 r = isl_reordering_extend(r, qp->div->n_row);
4141 if (!r)
4142 goto error;
4143
4144 qp->div = reorder_divs(qp->div, isl_reordering_copy(r));
4145 if (!qp->div)
4146 goto error;
4147
4148 qp->upoly = reorder(qp->upoly, r->pos);
4149 if (!qp->upoly)
4150 goto error;
4151
4152 qp = isl_qpolynomial_reset_domain_space(qp, isl_space_copy(r->dim));
4153
4154 isl_reordering_free(r);
4155 return qp;
4156 error:
4157 isl_qpolynomial_free(qp);
4158 isl_reordering_free(r);
4159 return NULL;
4160 }
4161
4162 __isl_give isl_qpolynomial *isl_qpolynomial_align_params(
4163 __isl_take isl_qpolynomial *qp, __isl_take isl_space *model)
4164 {
4165 if (!qp || !model)
4166 goto error;
4167
4168 if (!isl_space_match(qp->dim, isl_dim_param, model, isl_dim_param)) {
4169 isl_reordering *exp;
4170
4171 model = isl_space_drop_dims(model, isl_dim_in,
4172 0, isl_space_dim(model, isl_dim_in));
4173 model = isl_space_drop_dims(model, isl_dim_out,
4174 0, isl_space_dim(model, isl_dim_out));
4175 exp = isl_parameter_alignment_reordering(qp->dim, model);
4176 exp = isl_reordering_extend_space(exp,
4177 isl_qpolynomial_get_domain_space(qp));
4178 qp = isl_qpolynomial_realign_domain(qp, exp);
4179 }
4180
4181 isl_space_free(model);
4182 return qp;
4183 error:
4184 isl_space_free(model);
4185 isl_qpolynomial_free(qp);
4186 return NULL;
4187 }
4188
4189 struct isl_split_periods_data {
4190 int max_periods;
4191 isl_pw_qpolynomial *res;
4192 };
4193
4194 /* Create a slice where the integer division "div" has the fixed value "v".
4195 * In particular, if "div" refers to floor(f/m), then create a slice
4196 *
4197 * m v <= f <= m v + (m - 1)
4198 *
4199 * or
4200 *
4201 * f - m v >= 0
4202 * -f + m v + (m - 1) >= 0
4203 */
4204 static __isl_give isl_set *set_div_slice(__isl_take isl_space *dim,
4205 __isl_keep isl_qpolynomial *qp, int div, isl_int v)
4206 {
4207 int total;
4208 isl_basic_set *bset = NULL;
4209 int k;
4210
4211 if (!dim || !qp)
4212 goto error;
4213
4214 total = isl_space_dim(dim, isl_dim_all);
4215 bset = isl_basic_set_alloc_space(isl_space_copy(dim), 0, 0, 2);
4216
4217 k = isl_basic_set_alloc_inequality(bset);
4218 if (k < 0)
4219 goto error;
4220 isl_seq_cpy(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4221 isl_int_submul(bset->ineq[k][0], v, qp->div->row[div][0]);
4222
4223 k = isl_basic_set_alloc_inequality(bset);
4224 if (k < 0)
4225 goto error;
4226 isl_seq_neg(bset->ineq[k], qp->div->row[div] + 1, 1 + total);
4227 isl_int_addmul(bset->ineq[k][0], v, qp->div->row[div][0]);
4228 isl_int_add(bset->ineq[k][0], bset->ineq[k][0], qp->div->row[div][0]);
4229 isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
4230
4231 isl_space_free(dim);
4232 return isl_set_from_basic_set(bset);
4233 error:
4234 isl_basic_set_free(bset);
4235 isl_space_free(dim);
4236 return NULL;
4237 }
4238
4239 static int split_periods(__isl_take isl_set *set,
4240 __isl_take isl_qpolynomial *qp, void *user);
4241
4242 /* Create a slice of the domain "set" such that integer division "div"
4243 * has the fixed value "v" and add the results to data->res,
4244 * replacing the integer division by "v" in "qp".
4245 */
4246 static int set_div(__isl_take isl_set *set,
4247 __isl_take isl_qpolynomial *qp, int div, isl_int v,
4248 struct isl_split_periods_data *data)
4249 {
4250 int i;
4251 int total;
4252 isl_set *slice;
4253 struct isl_upoly *cst;
4254
4255 slice = set_div_slice(isl_set_get_space(set), qp, div, v);
4256 set = isl_set_intersect(set, slice);
4257
4258 if (!qp)
4259 goto error;
4260
4261 total = isl_space_dim(qp->dim, isl_dim_all);
4262
4263 for (i = div + 1; i < qp->div->n_row; ++i) {
4264 if (isl_int_is_zero(qp->div->row[i][2 + total + div]))
4265 continue;
4266 isl_int_addmul(qp->div->row[i][1],
4267 qp->div->row[i][2 + total + div], v);
4268 isl_int_set_si(qp->div->row[i][2 + total + div], 0);
4269 }
4270
4271 cst = isl_upoly_rat_cst(qp->dim->ctx, v, qp->dim->ctx->one);
4272 qp = substitute_div(qp, div, cst);
4273
4274 return split_periods(set, qp, data);
4275 error:
4276 isl_set_free(set);
4277 isl_qpolynomial_free(qp);
4278 return -1;
4279 }
4280
4281 /* Split the domain "set" such that integer division "div"
4282 * has a fixed value (ranging from "min" to "max") on each slice
4283 * and add the results to data->res.
4284 */
4285 static int split_div(__isl_take isl_set *set,
4286 __isl_take isl_qpolynomial *qp, int div, isl_int min, isl_int max,
4287 struct isl_split_periods_data *data)
4288 {
4289 for (; isl_int_le(min, max); isl_int_add_ui(min, min, 1)) {
4290 isl_set *set_i = isl_set_copy(set);
4291 isl_qpolynomial *qp_i = isl_qpolynomial_copy(qp);
4292
4293 if (set_div(set_i, qp_i, div, min, data) < 0)
4294 goto error;
4295 }
4296 isl_set_free(set);
4297 isl_qpolynomial_free(qp);
4298 return 0;
4299 error:
4300 isl_set_free(set);
4301 isl_qpolynomial_free(qp);
4302 return -1;
4303 }
4304
4305 /* If "qp" refers to any integer division
4306 * that can only attain "max_periods" distinct values on "set"
4307 * then split the domain along those distinct values.
4308 * Add the results (or the original if no splitting occurs)
4309 * to data->res.
4310 */
4311 static int split_periods(__isl_take isl_set *set,
4312 __isl_take isl_qpolynomial *qp, void *user)
4313 {
4314 int i;
4315 isl_pw_qpolynomial *pwqp;
4316 struct isl_split_periods_data *data;
4317 isl_int min, max;
4318 int total;
4319 int r = 0;
4320
4321 data = (struct isl_split_periods_data *)user;
4322
4323 if (!set || !qp)
4324 goto error;
4325
4326 if (qp->div->n_row == 0) {
4327 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4328 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4329 return 0;
4330 }
4331
4332 isl_int_init(min);
4333 isl_int_init(max);
4334 total = isl_space_dim(qp->dim, isl_dim_all);
4335 for (i = 0; i < qp->div->n_row; ++i) {
4336 enum isl_lp_result lp_res;
4337
4338 if (isl_seq_first_non_zero(qp->div->row[i] + 2 + total,
4339 qp->div->n_row) != -1)
4340 continue;
4341
4342 lp_res = isl_set_solve_lp(set, 0, qp->div->row[i] + 1,
4343 set->ctx->one, &min, NULL, NULL);
4344 if (lp_res == isl_lp_error)
4345 goto error2;
4346 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4347 continue;
4348 isl_int_fdiv_q(min, min, qp->div->row[i][0]);
4349
4350 lp_res = isl_set_solve_lp(set, 1, qp->div->row[i] + 1,
4351 set->ctx->one, &max, NULL, NULL);
4352 if (lp_res == isl_lp_error)
4353 goto error2;
4354 if (lp_res == isl_lp_unbounded || lp_res == isl_lp_empty)
4355 continue;
4356 isl_int_fdiv_q(max, max, qp->div->row[i][0]);
4357
4358 isl_int_sub(max, max, min);
4359 if (isl_int_cmp_si(max, data->max_periods) < 0) {
4360 isl_int_add(max, max, min);
4361 break;
4362 }
4363 }
4364
4365 if (i < qp->div->n_row) {
4366 r = split_div(set, qp, i, min, max, data);
4367 } else {
4368 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4369 data->res = isl_pw_qpolynomial_add_disjoint(data->res, pwqp);
4370 }
4371
4372 isl_int_clear(max);
4373 isl_int_clear(min);
4374
4375 return r;
4376 error2:
4377 isl_int_clear(max);
4378 isl_int_clear(min);
4379 error:
4380 isl_set_free(set);
4381 isl_qpolynomial_free(qp);
4382 return -1;
4383 }
4384
4385 /* If any quasi-polynomial in pwqp refers to any integer division
4386 * that can only attain "max_periods" distinct values on its domain
4387 * then split the domain along those distinct values.
4388 */
4389 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_split_periods(
4390 __isl_take isl_pw_qpolynomial *pwqp, int max_periods)
4391 {
4392 struct isl_split_periods_data data;
4393
4394 data.max_periods = max_periods;
4395 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4396
4397 if (isl_pw_qpolynomial_foreach_piece(pwqp, &split_periods, &data) < 0)
4398 goto error;
4399
4400 isl_pw_qpolynomial_free(pwqp);
4401
4402 return data.res;
4403 error:
4404 isl_pw_qpolynomial_free(data.res);
4405 isl_pw_qpolynomial_free(pwqp);
4406 return NULL;
4407 }
4408
4409 /* Construct a piecewise quasipolynomial that is constant on the given
4410 * domain. In particular, it is
4411 * 0 if cst == 0
4412 * 1 if cst == 1
4413 * infinity if cst == -1
4414 */
4415 static __isl_give isl_pw_qpolynomial *constant_on_domain(
4416 __isl_take isl_basic_set *bset, int cst)
4417 {
4418 isl_space *dim;
4419 isl_qpolynomial *qp;
4420
4421 if (!bset)
4422 return NULL;
4423
4424 bset = isl_basic_set_params(bset);
4425 dim = isl_basic_set_get_space(bset);
4426 if (cst < 0)
4427 qp = isl_qpolynomial_infty_on_domain(dim);
4428 else if (cst == 0)
4429 qp = isl_qpolynomial_zero_on_domain(dim);
4430 else
4431 qp = isl_qpolynomial_one_on_domain(dim);
4432 return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset), qp);
4433 }
4434
4435 /* Factor bset, call fn on each of the factors and return the product.
4436 *
4437 * If no factors can be found, simply call fn on the input.
4438 * Otherwise, construct the factors based on the factorizer,
4439 * call fn on each factor and compute the product.
4440 */
4441 static __isl_give isl_pw_qpolynomial *compressed_multiplicative_call(
4442 __isl_take isl_basic_set *bset,
4443 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4444 {
4445 int i, n;
4446 isl_space *dim;
4447 isl_set *set;
4448 isl_factorizer *f;
4449 isl_qpolynomial *qp;
4450 isl_pw_qpolynomial *pwqp;
4451 unsigned nparam;
4452 unsigned nvar;
4453
4454 f = isl_basic_set_factorizer(bset);
4455 if (!f)
4456 goto error;
4457 if (f->n_group == 0) {
4458 isl_factorizer_free(f);
4459 return fn(bset);
4460 }
4461
4462 nparam = isl_basic_set_dim(bset, isl_dim_param);
4463 nvar = isl_basic_set_dim(bset, isl_dim_set);
4464
4465 dim = isl_basic_set_get_space(bset);
4466 dim = isl_space_domain(dim);
4467 set = isl_set_universe(isl_space_copy(dim));
4468 qp = isl_qpolynomial_one_on_domain(dim);
4469 pwqp = isl_pw_qpolynomial_alloc(set, qp);
4470
4471 bset = isl_morph_basic_set(isl_morph_copy(f->morph), bset);
4472
4473 for (i = 0, n = 0; i < f->n_group; ++i) {
4474 isl_basic_set *bset_i;
4475 isl_pw_qpolynomial *pwqp_i;
4476
4477 bset_i = isl_basic_set_copy(bset);
4478 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4479 nparam + n + f->len[i], nvar - n - f->len[i]);
4480 bset_i = isl_basic_set_drop_constraints_involving(bset_i,
4481 nparam, n);
4482 bset_i = isl_basic_set_drop(bset_i, isl_dim_set,
4483 n + f->len[i], nvar - n - f->len[i]);
4484 bset_i = isl_basic_set_drop(bset_i, isl_dim_set, 0, n);
4485
4486 pwqp_i = fn(bset_i);
4487 pwqp = isl_pw_qpolynomial_mul(pwqp, pwqp_i);
4488
4489 n += f->len[i];
4490 }
4491
4492 isl_basic_set_free(bset);
4493 isl_factorizer_free(f);
4494
4495 return pwqp;
4496 error:
4497 isl_basic_set_free(bset);
4498 return NULL;
4499 }
4500
4501 /* Factor bset, call fn on each of the factors and return the product.
4502 * The function is assumed to evaluate to zero on empty domains,
4503 * to one on zero-dimensional domains and to infinity on unbounded domains
4504 * and will not be called explicitly on zero-dimensional or unbounded domains.
4505 *
4506 * We first check for some special cases and remove all equalities.
4507 * Then we hand over control to compressed_multiplicative_call.
4508 */
4509 __isl_give isl_pw_qpolynomial *isl_basic_set_multiplicative_call(
4510 __isl_take isl_basic_set *bset,
4511 __isl_give isl_pw_qpolynomial *(*fn)(__isl_take isl_basic_set *bset))
4512 {
4513 int bounded;
4514 isl_morph *morph;
4515 isl_pw_qpolynomial *pwqp;
4516
4517 if (!bset)
4518 return NULL;
4519
4520 if (isl_basic_set_plain_is_empty(bset))
4521 return constant_on_domain(bset, 0);
4522
4523 if (isl_basic_set_dim(bset, isl_dim_set) == 0)
4524 return constant_on_domain(bset, 1);
4525
4526 bounded = isl_basic_set_is_bounded(bset);
4527 if (bounded < 0)
4528 goto error;
4529 if (!bounded)
4530 return constant_on_domain(bset, -1);
4531
4532 if (bset->n_eq == 0)
4533 return compressed_multiplicative_call(bset, fn);
4534
4535 morph = isl_basic_set_full_compression(bset);
4536 bset = isl_morph_basic_set(isl_morph_copy(morph), bset);
4537
4538 pwqp = compressed_multiplicative_call(bset, fn);
4539
4540 morph = isl_morph_dom_params(morph);
4541 morph = isl_morph_ran_params(morph);
4542 morph = isl_morph_inverse(morph);
4543
4544 pwqp = isl_pw_qpolynomial_morph_domain(pwqp, morph);
4545
4546 return pwqp;
4547 error:
4548 isl_basic_set_free(bset);
4549 return NULL;
4550 }
4551
4552 /* Drop all floors in "qp", turning each integer division [a/m] into
4553 * a rational division a/m. If "down" is set, then the integer division
4554 * is replaced by (a-(m-1))/m instead.
4555 */
4556 static __isl_give isl_qpolynomial *qp_drop_floors(
4557 __isl_take isl_qpolynomial *qp, int down)
4558 {
4559 int i;
4560 struct isl_upoly *s;
4561
4562 if (!qp)
4563 return NULL;
4564 if (qp->div->n_row == 0)
4565 return qp;
4566
4567 qp = isl_qpolynomial_cow(qp);
4568 if (!qp)
4569 return NULL;
4570
4571 for (i = qp->div->n_row - 1; i >= 0; --i) {
4572 if (down) {
4573 isl_int_sub(qp->div->row[i][1],
4574 qp->div->row[i][1], qp->div->row[i][0]);
4575 isl_int_add_ui(qp->div->row[i][1],
4576 qp->div->row[i][1], 1);
4577 }
4578 s = isl_upoly_from_affine(qp->dim->ctx, qp->div->row[i] + 1,
4579 qp->div->row[i][0], qp->div->n_col - 1);
4580 qp = substitute_div(qp, i, s);
4581 if (!qp)
4582 return NULL;
4583 }
4584
4585 return qp;
4586 }
4587
4588 /* Drop all floors in "pwqp", turning each integer division [a/m] into
4589 * a rational division a/m.
4590 */
4591 static __isl_give isl_pw_qpolynomial *pwqp_drop_floors(
4592 __isl_take isl_pw_qpolynomial *pwqp)
4593 {
4594 int i;
4595
4596 if (!pwqp)
4597 return NULL;
4598
4599 if (isl_pw_qpolynomial_is_zero(pwqp))
4600 return pwqp;
4601
4602 pwqp = isl_pw_qpolynomial_cow(pwqp);
4603 if (!pwqp)
4604 return NULL;
4605
4606 for (i = 0; i < pwqp->n; ++i) {
4607 pwqp->p[i].qp = qp_drop_floors(pwqp->p[i].qp, 0);
4608 if (!pwqp->p[i].qp)
4609 goto error;
4610 }
4611
4612 return pwqp;
4613 error:
4614 isl_pw_qpolynomial_free(pwqp);
4615 return NULL;
4616 }
4617
4618 /* Adjust all the integer divisions in "qp" such that they are at least
4619 * one over the given orthant (identified by "signs"). This ensures
4620 * that they will still be non-negative even after subtracting (m-1)/m.
4621 *
4622 * In particular, f is replaced by f' + v, changing f = [a/m]
4623 * to f' = [(a - m v)/m].
4624 * If the constant term k in a is smaller than m,
4625 * the constant term of v is set to floor(k/m) - 1.
4626 * For any other term, if the coefficient c and the variable x have
4627 * the same sign, then no changes are needed.
4628 * Otherwise, if the variable is positive (and c is negative),
4629 * then the coefficient of x in v is set to floor(c/m).
4630 * If the variable is negative (and c is positive),
4631 * then the coefficient of x in v is set to ceil(c/m).
4632 */
4633 static __isl_give isl_qpolynomial *make_divs_pos(__isl_take isl_qpolynomial *qp,
4634 int *signs)
4635 {
4636 int i, j;
4637 int total;
4638 isl_vec *v = NULL;
4639 struct isl_upoly *s;
4640
4641 qp = isl_qpolynomial_cow(qp);
4642 if (!qp)
4643 return NULL;
4644 qp->div = isl_mat_cow(qp->div);
4645 if (!qp->div)
4646 goto error;
4647
4648 total = isl_space_dim(qp->dim, isl_dim_all);
4649 v = isl_vec_alloc(qp->div->ctx, qp->div->n_col - 1);
4650
4651 for (i = 0; i < qp->div->n_row; ++i) {
4652 isl_int *row = qp->div->row[i];
4653 v = isl_vec_clr(v);
4654 if (!v)
4655 goto error;
4656 if (isl_int_lt(row[1], row[0])) {
4657 isl_int_fdiv_q(v->el[0], row[1], row[0]);
4658 isl_int_sub_ui(v->el[0], v->el[0], 1);
4659 isl_int_submul(row[1], row[0], v->el[0]);
4660 }
4661 for (j = 0; j < total; ++j) {
4662 if (isl_int_sgn(row[2 + j]) * signs[j] >= 0)
4663 continue;
4664 if (signs[j] < 0)
4665 isl_int_cdiv_q(v->el[1 + j], row[2 + j], row[0]);
4666 else
4667 isl_int_fdiv_q(v->el[1 + j], row[2 + j], row[0]);
4668 isl_int_submul(row[2 + j], row[0], v->el[1 + j]);
4669 }
4670 for (j = 0; j < i; ++j) {
4671 if (isl_int_sgn(row[2 + total + j]) >= 0)
4672 continue;
4673 isl_int_fdiv_q(v->el[1 + total + j],
4674 row[2 + total + j], row[0]);
4675 isl_int_submul(row[2 + total + j],
4676 row[0], v->el[1 + total + j]);
4677 }
4678 for (j = i + 1; j < qp->div->n_row; ++j) {
4679 if (isl_int_is_zero(qp->div->row[j][2 + total + i]))
4680 continue;
4681 isl_seq_combine(qp->div->row[j] + 1,
4682 qp->div->ctx->one, qp->div->row[j] + 1,
4683 qp->div->row[j][2 + total + i], v->el, v->size);
4684 }
4685 isl_int_set_si(v->el[1 + total + i], 1);
4686 s = isl_upoly_from_affine(qp->dim->ctx, v->el,
4687 qp->div->ctx->one, v->size);
4688 qp->upoly = isl_upoly_subs(qp->upoly, total + i, 1, &s);
4689 isl_upoly_free(s);
4690 if (!qp->upoly)
4691 goto error;
4692 }
4693
4694 isl_vec_free(v);
4695 return qp;
4696 error:
4697 isl_vec_free(v);
4698 isl_qpolynomial_free(qp);
4699 return NULL;
4700 }
4701
4702 struct isl_to_poly_data {
4703 int sign;
4704 isl_pw_qpolynomial *res;
4705 isl_qpolynomial *qp;
4706 };
4707
4708 /* Appoximate data->qp by a polynomial on the orthant identified by "signs".
4709 * We first make all integer divisions positive and then split the
4710 * quasipolynomials into terms with sign data->sign (the direction
4711 * of the requested approximation) and terms with the opposite sign.
4712 * In the first set of terms, each integer division [a/m] is
4713 * overapproximated by a/m, while in the second it is underapproximated
4714 * by (a-(m-1))/m.
4715 */
4716 static int to_polynomial_on_orthant(__isl_take isl_set *orthant, int *signs,
4717 void *user)
4718 {
4719 struct isl_to_poly_data *data = user;
4720 isl_pw_qpolynomial *t;
4721 isl_qpolynomial *qp, *up, *down;
4722
4723 qp = isl_qpolynomial_copy(data->qp);
4724 qp = make_divs_pos(qp, signs);
4725
4726 up = isl_qpolynomial_terms_of_sign(qp, signs, data->sign);
4727 up = qp_drop_floors(up, 0);
4728 down = isl_qpolynomial_terms_of_sign(qp, signs, -data->sign);
4729 down = qp_drop_floors(down, 1);
4730
4731 isl_qpolynomial_free(qp);
4732 qp = isl_qpolynomial_add(up, down);
4733
4734 t = isl_pw_qpolynomial_alloc(orthant, qp);
4735 data->res = isl_pw_qpolynomial_add_disjoint(data->res, t);
4736
4737 return 0;
4738 }
4739
4740 /* Approximate each quasipolynomial by a polynomial. If "sign" is positive,
4741 * the polynomial will be an overapproximation. If "sign" is negative,
4742 * it will be an underapproximation. If "sign" is zero, the approximation
4743 * will lie somewhere in between.
4744 *
4745 * In particular, is sign == 0, we simply drop the floors, turning
4746 * the integer divisions into rational divisions.
4747 * Otherwise, we split the domains into orthants, make all integer divisions
4748 * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
4749 * depending on the requested sign and the sign of the term in which
4750 * the integer division appears.
4751 */
4752 __isl_give isl_pw_qpolynomial *isl_pw_qpolynomial_to_polynomial(
4753 __isl_take isl_pw_qpolynomial *pwqp, int sign)
4754 {
4755 int i;
4756 struct isl_to_poly_data data;
4757
4758 if (sign == 0)
4759 return pwqp_drop_floors(pwqp);
4760
4761 if (!pwqp)
4762 return NULL;
4763
4764 data.sign = sign;
4765 data.res = isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp));
4766
4767 for (i = 0; i < pwqp->n; ++i) {
4768 if (pwqp->p[i].qp->div->n_row == 0) {
4769 isl_pw_qpolynomial *t;
4770 t = isl_pw_qpolynomial_alloc(
4771 isl_set_copy(pwqp->p[i].set),
4772 isl_qpolynomial_copy(pwqp->p[i].qp));
4773 data.res = isl_pw_qpolynomial_add_disjoint(data.res, t);
4774 continue;
4775 }
4776 data.qp = pwqp->p[i].qp;
4777 if (isl_set_foreach_orthant(pwqp->p[i].set,
4778 &to_polynomial_on_orthant, &data) < 0)
4779 goto error;
4780 }
4781
4782 isl_pw_qpolynomial_free(pwqp);
4783
4784 return data.res;
4785 error:
4786 isl_pw_qpolynomial_free(pwqp);
4787 isl_pw_qpolynomial_free(data.res);
4788 return NULL;
4789 }
4790
4791 static int poly_entry(void **entry, void *user)
4792 {
4793 int *sign = user;
4794 isl_pw_qpolynomial **pwqp = (isl_pw_qpolynomial **)entry;
4795
4796 *pwqp = isl_pw_qpolynomial_to_polynomial(*pwqp, *sign);
4797
4798 return *pwqp ? 0 : -1;
4799 }
4800
4801 __isl_give isl_union_pw_qpolynomial *isl_union_pw_qpolynomial_to_polynomial(
4802 __isl_take isl_union_pw_qpolynomial *upwqp, int sign)
4803 {
4804 upwqp = isl_union_pw_qpolynomial_cow(upwqp);
4805 if (!upwqp)
4806 return NULL;
4807
4808 if (isl_hash_table_foreach(upwqp->space->ctx, &upwqp->table,
4809 &poly_entry, &sign) < 0)
4810 goto error;
4811
4812 return upwqp;
4813 error:
4814 isl_union_pw_qpolynomial_free(upwqp);
4815 return NULL;
4816 }
4817
4818 __isl_give isl_basic_map *isl_basic_map_from_qpolynomial(
4819 __isl_take isl_qpolynomial *qp)
4820 {
4821 int i, k;
4822 isl_space *dim;
4823 isl_vec *aff = NULL;
4824 isl_basic_map *bmap = NULL;
4825 unsigned pos;
4826 unsigned n_div;
4827
4828 if (!qp)
4829 return NULL;
4830 if (!isl_upoly_is_affine(qp->upoly))
4831 isl_die(qp->dim->ctx, isl_error_invalid,
4832 "input quasi-polynomial not affine", goto error);
4833 aff = isl_qpolynomial_extract_affine(qp);
4834 if (!aff)
4835 goto error;
4836 dim = isl_qpolynomial_get_space(qp);
4837 pos = 1 + isl_space_offset(dim, isl_dim_out);
4838 n_div = qp->div->n_row;
4839 bmap = isl_basic_map_alloc_space(dim, n_div, 1, 2 * n_div);
4840
4841 for (i = 0; i < n_div; ++i) {
4842 k = isl_basic_map_alloc_div(bmap);
4843 if (k < 0)
4844 goto error;
4845 isl_seq_cpy(bmap->div[k], qp->div->row[i], qp->div->n_col);
4846 isl_int_set_si(bmap->div[k][qp->div->n_col], 0);
4847 if (isl_basic_map_add_div_constraints(bmap, k) < 0)
4848 goto error;
4849 }
4850 k = isl_basic_map_alloc_equality(bmap);
4851 if (k < 0)
4852 goto error;
4853 isl_int_neg(bmap->eq[k][pos], aff->el[0]);
4854 isl_seq_cpy(bmap->eq[k], aff->el + 1, pos);
4855 isl_seq_cpy(bmap->eq[k] + pos + 1, aff->el + 1 + pos, n_div);
4856
4857 isl_vec_free(aff);
4858 isl_qpolynomial_free(qp);
4859 bmap = isl_basic_map_finalize(bmap);
4860 return bmap;
4861 error:
4862 isl_vec_free(aff);
4863 isl_qpolynomial_free(qp);
4864 isl_basic_map_free(bmap);
4865 return NULL;
4866 }
Integer set library