2 * Copyright 2010 INRIA Saclay
4 * Use of this software is governed by the MIT license
6 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
13 #include <isl_ctx_private.h>
14 #include <isl_map_private.h>
15 #include <isl_factorization.h>
16 #include <isl_lp_private.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>
31 static unsigned pos(__isl_keep isl_space
*dim
, enum isl_dim_type type
)
34 case isl_dim_param
: return 0;
35 case isl_dim_in
: return dim
->nparam
;
36 case isl_dim_out
: return dim
->nparam
+ dim
->n_in
;
41 int isl_upoly_is_cst(__isl_keep
struct isl_upoly
*up
)
49 __isl_keep
struct isl_upoly_cst
*isl_upoly_as_cst(__isl_keep
struct isl_upoly
*up
)
54 isl_assert(up
->ctx
, up
->var
< 0, return NULL
);
56 return (struct isl_upoly_cst
*)up
;
59 __isl_keep
struct isl_upoly_rec
*isl_upoly_as_rec(__isl_keep
struct isl_upoly
*up
)
64 isl_assert(up
->ctx
, up
->var
>= 0, return NULL
);
66 return (struct isl_upoly_rec
*)up
;
69 int isl_upoly_is_equal(__isl_keep
struct isl_upoly
*up1
,
70 __isl_keep
struct isl_upoly
*up2
)
73 struct isl_upoly_rec
*rec1
, *rec2
;
79 if (up1
->var
!= up2
->var
)
81 if (isl_upoly_is_cst(up1
)) {
82 struct isl_upoly_cst
*cst1
, *cst2
;
83 cst1
= isl_upoly_as_cst(up1
);
84 cst2
= isl_upoly_as_cst(up2
);
87 return isl_int_eq(cst1
->n
, cst2
->n
) &&
88 isl_int_eq(cst1
->d
, cst2
->d
);
91 rec1
= isl_upoly_as_rec(up1
);
92 rec2
= isl_upoly_as_rec(up2
);
96 if (rec1
->n
!= rec2
->n
)
99 for (i
= 0; i
< rec1
->n
; ++i
) {
100 int eq
= isl_upoly_is_equal(rec1
->p
[i
], rec2
->p
[i
]);
108 int isl_upoly_is_zero(__isl_keep
struct isl_upoly
*up
)
110 struct isl_upoly_cst
*cst
;
114 if (!isl_upoly_is_cst(up
))
117 cst
= isl_upoly_as_cst(up
);
121 return isl_int_is_zero(cst
->n
) && isl_int_is_pos(cst
->d
);
124 int isl_upoly_sgn(__isl_keep
struct isl_upoly
*up
)
126 struct isl_upoly_cst
*cst
;
130 if (!isl_upoly_is_cst(up
))
133 cst
= isl_upoly_as_cst(up
);
137 return isl_int_sgn(cst
->n
);
140 int isl_upoly_is_nan(__isl_keep
struct isl_upoly
*up
)
142 struct isl_upoly_cst
*cst
;
146 if (!isl_upoly_is_cst(up
))
149 cst
= isl_upoly_as_cst(up
);
153 return isl_int_is_zero(cst
->n
) && isl_int_is_zero(cst
->d
);
156 int isl_upoly_is_infty(__isl_keep
struct isl_upoly
*up
)
158 struct isl_upoly_cst
*cst
;
162 if (!isl_upoly_is_cst(up
))
165 cst
= isl_upoly_as_cst(up
);
169 return isl_int_is_pos(cst
->n
) && isl_int_is_zero(cst
->d
);
172 int isl_upoly_is_neginfty(__isl_keep
struct isl_upoly
*up
)
174 struct isl_upoly_cst
*cst
;
178 if (!isl_upoly_is_cst(up
))
181 cst
= isl_upoly_as_cst(up
);
185 return isl_int_is_neg(cst
->n
) && isl_int_is_zero(cst
->d
);
188 int isl_upoly_is_one(__isl_keep
struct isl_upoly
*up
)
190 struct isl_upoly_cst
*cst
;
194 if (!isl_upoly_is_cst(up
))
197 cst
= isl_upoly_as_cst(up
);
201 return isl_int_eq(cst
->n
, cst
->d
) && isl_int_is_pos(cst
->d
);
204 int isl_upoly_is_negone(__isl_keep
struct isl_upoly
*up
)
206 struct isl_upoly_cst
*cst
;
210 if (!isl_upoly_is_cst(up
))
213 cst
= isl_upoly_as_cst(up
);
217 return isl_int_is_negone(cst
->n
) && isl_int_is_one(cst
->d
);
220 __isl_give
struct isl_upoly_cst
*isl_upoly_cst_alloc(struct isl_ctx
*ctx
)
222 struct isl_upoly_cst
*cst
;
224 cst
= isl_alloc_type(ctx
, struct isl_upoly_cst
);
233 isl_int_init(cst
->n
);
234 isl_int_init(cst
->d
);
239 __isl_give
struct isl_upoly
*isl_upoly_zero(struct isl_ctx
*ctx
)
241 struct isl_upoly_cst
*cst
;
243 cst
= isl_upoly_cst_alloc(ctx
);
247 isl_int_set_si(cst
->n
, 0);
248 isl_int_set_si(cst
->d
, 1);
253 __isl_give
struct isl_upoly
*isl_upoly_one(struct isl_ctx
*ctx
)
255 struct isl_upoly_cst
*cst
;
257 cst
= isl_upoly_cst_alloc(ctx
);
261 isl_int_set_si(cst
->n
, 1);
262 isl_int_set_si(cst
->d
, 1);
267 __isl_give
struct isl_upoly
*isl_upoly_infty(struct isl_ctx
*ctx
)
269 struct isl_upoly_cst
*cst
;
271 cst
= isl_upoly_cst_alloc(ctx
);
275 isl_int_set_si(cst
->n
, 1);
276 isl_int_set_si(cst
->d
, 0);
281 __isl_give
struct isl_upoly
*isl_upoly_neginfty(struct isl_ctx
*ctx
)
283 struct isl_upoly_cst
*cst
;
285 cst
= isl_upoly_cst_alloc(ctx
);
289 isl_int_set_si(cst
->n
, -1);
290 isl_int_set_si(cst
->d
, 0);
295 __isl_give
struct isl_upoly
*isl_upoly_nan(struct isl_ctx
*ctx
)
297 struct isl_upoly_cst
*cst
;
299 cst
= isl_upoly_cst_alloc(ctx
);
303 isl_int_set_si(cst
->n
, 0);
304 isl_int_set_si(cst
->d
, 0);
309 __isl_give
struct isl_upoly
*isl_upoly_rat_cst(struct isl_ctx
*ctx
,
310 isl_int n
, isl_int d
)
312 struct isl_upoly_cst
*cst
;
314 cst
= isl_upoly_cst_alloc(ctx
);
318 isl_int_set(cst
->n
, n
);
319 isl_int_set(cst
->d
, d
);
324 __isl_give
struct isl_upoly_rec
*isl_upoly_alloc_rec(struct isl_ctx
*ctx
,
327 struct isl_upoly_rec
*rec
;
329 isl_assert(ctx
, var
>= 0, return NULL
);
330 isl_assert(ctx
, size
>= 0, return NULL
);
331 rec
= isl_calloc(ctx
, struct isl_upoly_rec
,
332 sizeof(struct isl_upoly_rec
) +
333 size
* sizeof(struct isl_upoly
*));
348 __isl_give isl_qpolynomial
*isl_qpolynomial_reset_domain_space(
349 __isl_take isl_qpolynomial
*qp
, __isl_take isl_space
*dim
)
351 qp
= isl_qpolynomial_cow(qp
);
355 isl_space_free(qp
->dim
);
360 isl_qpolynomial_free(qp
);
365 /* Reset the space of "qp". This function is called from isl_pw_templ.c
366 * and doesn't know if the space of an element object is represented
367 * directly or through its domain. It therefore passes along both.
369 __isl_give isl_qpolynomial
*isl_qpolynomial_reset_space_and_domain(
370 __isl_take isl_qpolynomial
*qp
, __isl_take isl_space
*space
,
371 __isl_take isl_space
*domain
)
373 isl_space_free(space
);
374 return isl_qpolynomial_reset_domain_space(qp
, domain
);
377 isl_ctx
*isl_qpolynomial_get_ctx(__isl_keep isl_qpolynomial
*qp
)
379 return qp
? qp
->dim
->ctx
: NULL
;
382 __isl_give isl_space
*isl_qpolynomial_get_domain_space(
383 __isl_keep isl_qpolynomial
*qp
)
385 return qp
? isl_space_copy(qp
->dim
) : NULL
;
388 __isl_give isl_space
*isl_qpolynomial_get_space(__isl_keep isl_qpolynomial
*qp
)
393 space
= isl_space_copy(qp
->dim
);
394 space
= isl_space_from_domain(space
);
395 space
= isl_space_add_dims(space
, isl_dim_out
, 1);
399 /* Externally, an isl_qpolynomial has a map space, but internally, the
400 * ls field corresponds to the domain of that space.
402 unsigned isl_qpolynomial_dim(__isl_keep isl_qpolynomial
*qp
,
403 enum isl_dim_type type
)
407 if (type
== isl_dim_out
)
409 if (type
== isl_dim_in
)
411 return isl_space_dim(qp
->dim
, type
);
414 int isl_qpolynomial_is_zero(__isl_keep isl_qpolynomial
*qp
)
416 return qp
? isl_upoly_is_zero(qp
->upoly
) : -1;
419 int isl_qpolynomial_is_one(__isl_keep isl_qpolynomial
*qp
)
421 return qp
? isl_upoly_is_one(qp
->upoly
) : -1;
424 int isl_qpolynomial_is_nan(__isl_keep isl_qpolynomial
*qp
)
426 return qp
? isl_upoly_is_nan(qp
->upoly
) : -1;
429 int isl_qpolynomial_is_infty(__isl_keep isl_qpolynomial
*qp
)
431 return qp
? isl_upoly_is_infty(qp
->upoly
) : -1;
434 int isl_qpolynomial_is_neginfty(__isl_keep isl_qpolynomial
*qp
)
436 return qp
? isl_upoly_is_neginfty(qp
->upoly
) : -1;
439 int isl_qpolynomial_sgn(__isl_keep isl_qpolynomial
*qp
)
441 return qp
? isl_upoly_sgn(qp
->upoly
) : 0;
444 static void upoly_free_cst(__isl_take
struct isl_upoly_cst
*cst
)
446 isl_int_clear(cst
->n
);
447 isl_int_clear(cst
->d
);
450 static void upoly_free_rec(__isl_take
struct isl_upoly_rec
*rec
)
454 for (i
= 0; i
< rec
->n
; ++i
)
455 isl_upoly_free(rec
->p
[i
]);
458 __isl_give
struct isl_upoly
*isl_upoly_copy(__isl_keep
struct isl_upoly
*up
)
467 __isl_give
struct isl_upoly
*isl_upoly_dup_cst(__isl_keep
struct isl_upoly
*up
)
469 struct isl_upoly_cst
*cst
;
470 struct isl_upoly_cst
*dup
;
472 cst
= isl_upoly_as_cst(up
);
476 dup
= isl_upoly_as_cst(isl_upoly_zero(up
->ctx
));
479 isl_int_set(dup
->n
, cst
->n
);
480 isl_int_set(dup
->d
, cst
->d
);
485 __isl_give
struct isl_upoly
*isl_upoly_dup_rec(__isl_keep
struct isl_upoly
*up
)
488 struct isl_upoly_rec
*rec
;
489 struct isl_upoly_rec
*dup
;
491 rec
= isl_upoly_as_rec(up
);
495 dup
= isl_upoly_alloc_rec(up
->ctx
, up
->var
, rec
->n
);
499 for (i
= 0; i
< rec
->n
; ++i
) {
500 dup
->p
[i
] = isl_upoly_copy(rec
->p
[i
]);
508 isl_upoly_free(&dup
->up
);
512 __isl_give
struct isl_upoly
*isl_upoly_dup(__isl_keep
struct isl_upoly
*up
)
517 if (isl_upoly_is_cst(up
))
518 return isl_upoly_dup_cst(up
);
520 return isl_upoly_dup_rec(up
);
523 __isl_give
struct isl_upoly
*isl_upoly_cow(__isl_take
struct isl_upoly
*up
)
531 return isl_upoly_dup(up
);
534 void isl_upoly_free(__isl_take
struct isl_upoly
*up
)
543 upoly_free_cst((struct isl_upoly_cst
*)up
);
545 upoly_free_rec((struct isl_upoly_rec
*)up
);
547 isl_ctx_deref(up
->ctx
);
551 static void isl_upoly_cst_reduce(__isl_keep
struct isl_upoly_cst
*cst
)
556 isl_int_gcd(gcd
, cst
->n
, cst
->d
);
557 if (!isl_int_is_zero(gcd
) && !isl_int_is_one(gcd
)) {
558 isl_int_divexact(cst
->n
, cst
->n
, gcd
);
559 isl_int_divexact(cst
->d
, cst
->d
, gcd
);
564 __isl_give
struct isl_upoly
*isl_upoly_sum_cst(__isl_take
struct isl_upoly
*up1
,
565 __isl_take
struct isl_upoly
*up2
)
567 struct isl_upoly_cst
*cst1
;
568 struct isl_upoly_cst
*cst2
;
570 up1
= isl_upoly_cow(up1
);
574 cst1
= isl_upoly_as_cst(up1
);
575 cst2
= isl_upoly_as_cst(up2
);
577 if (isl_int_eq(cst1
->d
, cst2
->d
))
578 isl_int_add(cst1
->n
, cst1
->n
, cst2
->n
);
580 isl_int_mul(cst1
->n
, cst1
->n
, cst2
->d
);
581 isl_int_addmul(cst1
->n
, cst2
->n
, cst1
->d
);
582 isl_int_mul(cst1
->d
, cst1
->d
, cst2
->d
);
585 isl_upoly_cst_reduce(cst1
);
595 static __isl_give
struct isl_upoly
*replace_by_zero(
596 __isl_take
struct isl_upoly
*up
)
604 return isl_upoly_zero(ctx
);
607 static __isl_give
struct isl_upoly
*replace_by_constant_term(
608 __isl_take
struct isl_upoly
*up
)
610 struct isl_upoly_rec
*rec
;
611 struct isl_upoly
*cst
;
616 rec
= isl_upoly_as_rec(up
);
619 cst
= isl_upoly_copy(rec
->p
[0]);
627 __isl_give
struct isl_upoly
*isl_upoly_sum(__isl_take
struct isl_upoly
*up1
,
628 __isl_take
struct isl_upoly
*up2
)
631 struct isl_upoly_rec
*rec1
, *rec2
;
636 if (isl_upoly_is_nan(up1
)) {
641 if (isl_upoly_is_nan(up2
)) {
646 if (isl_upoly_is_zero(up1
)) {
651 if (isl_upoly_is_zero(up2
)) {
656 if (up1
->var
< up2
->var
)
657 return isl_upoly_sum(up2
, up1
);
659 if (up2
->var
< up1
->var
) {
660 struct isl_upoly_rec
*rec
;
661 if (isl_upoly_is_infty(up2
) || isl_upoly_is_neginfty(up2
)) {
665 up1
= isl_upoly_cow(up1
);
666 rec
= isl_upoly_as_rec(up1
);
669 rec
->p
[0] = isl_upoly_sum(rec
->p
[0], up2
);
671 up1
= replace_by_constant_term(up1
);
675 if (isl_upoly_is_cst(up1
))
676 return isl_upoly_sum_cst(up1
, up2
);
678 rec1
= isl_upoly_as_rec(up1
);
679 rec2
= isl_upoly_as_rec(up2
);
683 if (rec1
->n
< rec2
->n
)
684 return isl_upoly_sum(up2
, up1
);
686 up1
= isl_upoly_cow(up1
);
687 rec1
= isl_upoly_as_rec(up1
);
691 for (i
= rec2
->n
- 1; i
>= 0; --i
) {
692 rec1
->p
[i
] = isl_upoly_sum(rec1
->p
[i
],
693 isl_upoly_copy(rec2
->p
[i
]));
696 if (i
== rec1
->n
- 1 && isl_upoly_is_zero(rec1
->p
[i
])) {
697 isl_upoly_free(rec1
->p
[i
]);
703 up1
= replace_by_zero(up1
);
704 else if (rec1
->n
== 1)
705 up1
= replace_by_constant_term(up1
);
716 __isl_give
struct isl_upoly
*isl_upoly_cst_add_isl_int(
717 __isl_take
struct isl_upoly
*up
, isl_int v
)
719 struct isl_upoly_cst
*cst
;
721 up
= isl_upoly_cow(up
);
725 cst
= isl_upoly_as_cst(up
);
727 isl_int_addmul(cst
->n
, cst
->d
, v
);
732 __isl_give
struct isl_upoly
*isl_upoly_add_isl_int(
733 __isl_take
struct isl_upoly
*up
, isl_int v
)
735 struct isl_upoly_rec
*rec
;
740 if (isl_upoly_is_cst(up
))
741 return isl_upoly_cst_add_isl_int(up
, v
);
743 up
= isl_upoly_cow(up
);
744 rec
= isl_upoly_as_rec(up
);
748 rec
->p
[0] = isl_upoly_add_isl_int(rec
->p
[0], v
);
758 __isl_give
struct isl_upoly
*isl_upoly_cst_mul_isl_int(
759 __isl_take
struct isl_upoly
*up
, isl_int v
)
761 struct isl_upoly_cst
*cst
;
763 if (isl_upoly_is_zero(up
))
766 up
= isl_upoly_cow(up
);
770 cst
= isl_upoly_as_cst(up
);
772 isl_int_mul(cst
->n
, cst
->n
, v
);
777 __isl_give
struct isl_upoly
*isl_upoly_mul_isl_int(
778 __isl_take
struct isl_upoly
*up
, isl_int v
)
781 struct isl_upoly_rec
*rec
;
786 if (isl_upoly_is_cst(up
))
787 return isl_upoly_cst_mul_isl_int(up
, v
);
789 up
= isl_upoly_cow(up
);
790 rec
= isl_upoly_as_rec(up
);
794 for (i
= 0; i
< rec
->n
; ++i
) {
795 rec
->p
[i
] = isl_upoly_mul_isl_int(rec
->p
[i
], v
);
806 /* Multiply the constant polynomial "up" by "v".
808 static __isl_give
struct isl_upoly
*isl_upoly_cst_scale_val(
809 __isl_take
struct isl_upoly
*up
, __isl_keep isl_val
*v
)
811 struct isl_upoly_cst
*cst
;
813 if (isl_upoly_is_zero(up
))
816 up
= isl_upoly_cow(up
);
820 cst
= isl_upoly_as_cst(up
);
822 isl_int_mul(cst
->n
, cst
->n
, v
->n
);
823 isl_int_mul(cst
->d
, cst
->d
, v
->d
);
824 isl_upoly_cst_reduce(cst
);
829 /* Multiply the polynomial "up" by "v".
831 static __isl_give
struct isl_upoly
*isl_upoly_scale_val(
832 __isl_take
struct isl_upoly
*up
, __isl_keep isl_val
*v
)
835 struct isl_upoly_rec
*rec
;
840 if (isl_upoly_is_cst(up
))
841 return isl_upoly_cst_scale_val(up
, v
);
843 up
= isl_upoly_cow(up
);
844 rec
= isl_upoly_as_rec(up
);
848 for (i
= 0; i
< rec
->n
; ++i
) {
849 rec
->p
[i
] = isl_upoly_scale_val(rec
->p
[i
], v
);
860 __isl_give
struct isl_upoly
*isl_upoly_mul_cst(__isl_take
struct isl_upoly
*up1
,
861 __isl_take
struct isl_upoly
*up2
)
863 struct isl_upoly_cst
*cst1
;
864 struct isl_upoly_cst
*cst2
;
866 up1
= isl_upoly_cow(up1
);
870 cst1
= isl_upoly_as_cst(up1
);
871 cst2
= isl_upoly_as_cst(up2
);
873 isl_int_mul(cst1
->n
, cst1
->n
, cst2
->n
);
874 isl_int_mul(cst1
->d
, cst1
->d
, cst2
->d
);
876 isl_upoly_cst_reduce(cst1
);
886 __isl_give
struct isl_upoly
*isl_upoly_mul_rec(__isl_take
struct isl_upoly
*up1
,
887 __isl_take
struct isl_upoly
*up2
)
889 struct isl_upoly_rec
*rec1
;
890 struct isl_upoly_rec
*rec2
;
891 struct isl_upoly_rec
*res
= NULL
;
895 rec1
= isl_upoly_as_rec(up1
);
896 rec2
= isl_upoly_as_rec(up2
);
899 size
= rec1
->n
+ rec2
->n
- 1;
900 res
= isl_upoly_alloc_rec(up1
->ctx
, up1
->var
, size
);
904 for (i
= 0; i
< rec1
->n
; ++i
) {
905 res
->p
[i
] = isl_upoly_mul(isl_upoly_copy(rec2
->p
[0]),
906 isl_upoly_copy(rec1
->p
[i
]));
911 for (; i
< size
; ++i
) {
912 res
->p
[i
] = isl_upoly_zero(up1
->ctx
);
917 for (i
= 0; i
< rec1
->n
; ++i
) {
918 for (j
= 1; j
< rec2
->n
; ++j
) {
919 struct isl_upoly
*up
;
920 up
= isl_upoly_mul(isl_upoly_copy(rec2
->p
[j
]),
921 isl_upoly_copy(rec1
->p
[i
]));
922 res
->p
[i
+ j
] = isl_upoly_sum(res
->p
[i
+ j
], up
);
935 isl_upoly_free(&res
->up
);
939 __isl_give
struct isl_upoly
*isl_upoly_mul(__isl_take
struct isl_upoly
*up1
,
940 __isl_take
struct isl_upoly
*up2
)
945 if (isl_upoly_is_nan(up1
)) {
950 if (isl_upoly_is_nan(up2
)) {
955 if (isl_upoly_is_zero(up1
)) {
960 if (isl_upoly_is_zero(up2
)) {
965 if (isl_upoly_is_one(up1
)) {
970 if (isl_upoly_is_one(up2
)) {
975 if (up1
->var
< up2
->var
)
976 return isl_upoly_mul(up2
, up1
);
978 if (up2
->var
< up1
->var
) {
980 struct isl_upoly_rec
*rec
;
981 if (isl_upoly_is_infty(up2
) || isl_upoly_is_neginfty(up2
)) {
982 isl_ctx
*ctx
= up1
->ctx
;
985 return isl_upoly_nan(ctx
);
987 up1
= isl_upoly_cow(up1
);
988 rec
= isl_upoly_as_rec(up1
);
992 for (i
= 0; i
< rec
->n
; ++i
) {
993 rec
->p
[i
] = isl_upoly_mul(rec
->p
[i
],
994 isl_upoly_copy(up2
));
1002 if (isl_upoly_is_cst(up1
))
1003 return isl_upoly_mul_cst(up1
, up2
);
1005 return isl_upoly_mul_rec(up1
, up2
);
1007 isl_upoly_free(up1
);
1008 isl_upoly_free(up2
);
1012 __isl_give
struct isl_upoly
*isl_upoly_pow(__isl_take
struct isl_upoly
*up
,
1015 struct isl_upoly
*res
;
1023 res
= isl_upoly_copy(up
);
1025 res
= isl_upoly_one(up
->ctx
);
1027 while (power
>>= 1) {
1028 up
= isl_upoly_mul(up
, isl_upoly_copy(up
));
1030 res
= isl_upoly_mul(res
, isl_upoly_copy(up
));
1037 __isl_give isl_qpolynomial
*isl_qpolynomial_alloc(__isl_take isl_space
*dim
,
1038 unsigned n_div
, __isl_take
struct isl_upoly
*up
)
1040 struct isl_qpolynomial
*qp
= NULL
;
1046 if (!isl_space_is_set(dim
))
1047 isl_die(isl_space_get_ctx(dim
), isl_error_invalid
,
1048 "domain of polynomial should be a set", goto error
);
1050 total
= isl_space_dim(dim
, isl_dim_all
);
1052 qp
= isl_calloc_type(dim
->ctx
, struct isl_qpolynomial
);
1057 qp
->div
= isl_mat_alloc(dim
->ctx
, n_div
, 1 + 1 + total
+ n_div
);
1066 isl_space_free(dim
);
1068 isl_qpolynomial_free(qp
);
1072 __isl_give isl_qpolynomial
*isl_qpolynomial_copy(__isl_keep isl_qpolynomial
*qp
)
1081 __isl_give isl_qpolynomial
*isl_qpolynomial_dup(__isl_keep isl_qpolynomial
*qp
)
1083 struct isl_qpolynomial
*dup
;
1088 dup
= isl_qpolynomial_alloc(isl_space_copy(qp
->dim
), qp
->div
->n_row
,
1089 isl_upoly_copy(qp
->upoly
));
1092 isl_mat_free(dup
->div
);
1093 dup
->div
= isl_mat_copy(qp
->div
);
1099 isl_qpolynomial_free(dup
);
1103 __isl_give isl_qpolynomial
*isl_qpolynomial_cow(__isl_take isl_qpolynomial
*qp
)
1111 return isl_qpolynomial_dup(qp
);
1114 void *isl_qpolynomial_free(__isl_take isl_qpolynomial
*qp
)
1122 isl_space_free(qp
->dim
);
1123 isl_mat_free(qp
->div
);
1124 isl_upoly_free(qp
->upoly
);
1130 __isl_give
struct isl_upoly
*isl_upoly_var_pow(isl_ctx
*ctx
, int pos
, int power
)
1133 struct isl_upoly_rec
*rec
;
1134 struct isl_upoly_cst
*cst
;
1136 rec
= isl_upoly_alloc_rec(ctx
, pos
, 1 + power
);
1139 for (i
= 0; i
< 1 + power
; ++i
) {
1140 rec
->p
[i
] = isl_upoly_zero(ctx
);
1145 cst
= isl_upoly_as_cst(rec
->p
[power
]);
1146 isl_int_set_si(cst
->n
, 1);
1150 isl_upoly_free(&rec
->up
);
1154 /* r array maps original positions to new positions.
1156 static __isl_give
struct isl_upoly
*reorder(__isl_take
struct isl_upoly
*up
,
1160 struct isl_upoly_rec
*rec
;
1161 struct isl_upoly
*base
;
1162 struct isl_upoly
*res
;
1164 if (isl_upoly_is_cst(up
))
1167 rec
= isl_upoly_as_rec(up
);
1171 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
1173 base
= isl_upoly_var_pow(up
->ctx
, r
[up
->var
], 1);
1174 res
= reorder(isl_upoly_copy(rec
->p
[rec
->n
- 1]), r
);
1176 for (i
= rec
->n
- 2; i
>= 0; --i
) {
1177 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
1178 res
= isl_upoly_sum(res
, reorder(isl_upoly_copy(rec
->p
[i
]), r
));
1181 isl_upoly_free(base
);
1190 static int compatible_divs(__isl_keep isl_mat
*div1
, __isl_keep isl_mat
*div2
)
1195 isl_assert(div1
->ctx
, div1
->n_row
>= div2
->n_row
&&
1196 div1
->n_col
>= div2
->n_col
, return -1);
1198 if (div1
->n_row
== div2
->n_row
)
1199 return isl_mat_is_equal(div1
, div2
);
1201 n_row
= div1
->n_row
;
1202 n_col
= div1
->n_col
;
1203 div1
->n_row
= div2
->n_row
;
1204 div1
->n_col
= div2
->n_col
;
1206 equal
= isl_mat_is_equal(div1
, div2
);
1208 div1
->n_row
= n_row
;
1209 div1
->n_col
= n_col
;
1214 static int cmp_row(__isl_keep isl_mat
*div
, int i
, int j
)
1218 li
= isl_seq_last_non_zero(div
->row
[i
], div
->n_col
);
1219 lj
= isl_seq_last_non_zero(div
->row
[j
], div
->n_col
);
1224 return isl_seq_cmp(div
->row
[i
], div
->row
[j
], div
->n_col
);
1227 struct isl_div_sort_info
{
1232 static int div_sort_cmp(const void *p1
, const void *p2
)
1234 const struct isl_div_sort_info
*i1
, *i2
;
1235 i1
= (const struct isl_div_sort_info
*) p1
;
1236 i2
= (const struct isl_div_sort_info
*) p2
;
1238 return cmp_row(i1
->div
, i1
->row
, i2
->row
);
1241 /* Sort divs and remove duplicates.
1243 static __isl_give isl_qpolynomial
*sort_divs(__isl_take isl_qpolynomial
*qp
)
1248 struct isl_div_sort_info
*array
= NULL
;
1249 int *pos
= NULL
, *at
= NULL
;
1250 int *reordering
= NULL
;
1255 if (qp
->div
->n_row
<= 1)
1258 div_pos
= isl_space_dim(qp
->dim
, isl_dim_all
);
1260 array
= isl_alloc_array(qp
->div
->ctx
, struct isl_div_sort_info
,
1262 pos
= isl_alloc_array(qp
->div
->ctx
, int, qp
->div
->n_row
);
1263 at
= isl_alloc_array(qp
->div
->ctx
, int, qp
->div
->n_row
);
1264 len
= qp
->div
->n_col
- 2;
1265 reordering
= isl_alloc_array(qp
->div
->ctx
, int, len
);
1266 if (!array
|| !pos
|| !at
|| !reordering
)
1269 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
1270 array
[i
].div
= qp
->div
;
1276 qsort(array
, qp
->div
->n_row
, sizeof(struct isl_div_sort_info
),
1279 for (i
= 0; i
< div_pos
; ++i
)
1282 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
1283 if (pos
[array
[i
].row
] == i
)
1285 qp
->div
= isl_mat_swap_rows(qp
->div
, i
, pos
[array
[i
].row
]);
1286 pos
[at
[i
]] = pos
[array
[i
].row
];
1287 at
[pos
[array
[i
].row
]] = at
[i
];
1288 at
[i
] = array
[i
].row
;
1289 pos
[array
[i
].row
] = i
;
1293 for (i
= 0; i
< len
- div_pos
; ++i
) {
1295 isl_seq_eq(qp
->div
->row
[i
- skip
- 1],
1296 qp
->div
->row
[i
- skip
], qp
->div
->n_col
)) {
1297 qp
->div
= isl_mat_drop_rows(qp
->div
, i
- skip
, 1);
1298 isl_mat_col_add(qp
->div
, 2 + div_pos
+ i
- skip
- 1,
1299 2 + div_pos
+ i
- skip
);
1300 qp
->div
= isl_mat_drop_cols(qp
->div
,
1301 2 + div_pos
+ i
- skip
, 1);
1304 reordering
[div_pos
+ array
[i
].row
] = div_pos
+ i
- skip
;
1307 qp
->upoly
= reorder(qp
->upoly
, reordering
);
1309 if (!qp
->upoly
|| !qp
->div
)
1323 isl_qpolynomial_free(qp
);
1327 static __isl_give
struct isl_upoly
*expand(__isl_take
struct isl_upoly
*up
,
1328 int *exp
, int first
)
1331 struct isl_upoly_rec
*rec
;
1333 if (isl_upoly_is_cst(up
))
1336 if (up
->var
< first
)
1339 if (exp
[up
->var
- first
] == up
->var
- first
)
1342 up
= isl_upoly_cow(up
);
1346 up
->var
= exp
[up
->var
- first
] + first
;
1348 rec
= isl_upoly_as_rec(up
);
1352 for (i
= 0; i
< rec
->n
; ++i
) {
1353 rec
->p
[i
] = expand(rec
->p
[i
], exp
, first
);
1364 static __isl_give isl_qpolynomial
*with_merged_divs(
1365 __isl_give isl_qpolynomial
*(*fn
)(__isl_take isl_qpolynomial
*qp1
,
1366 __isl_take isl_qpolynomial
*qp2
),
1367 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
1371 isl_mat
*div
= NULL
;
1373 qp1
= isl_qpolynomial_cow(qp1
);
1374 qp2
= isl_qpolynomial_cow(qp2
);
1379 isl_assert(qp1
->div
->ctx
, qp1
->div
->n_row
>= qp2
->div
->n_row
&&
1380 qp1
->div
->n_col
>= qp2
->div
->n_col
, goto error
);
1382 exp1
= isl_alloc_array(qp1
->div
->ctx
, int, qp1
->div
->n_row
);
1383 exp2
= isl_alloc_array(qp2
->div
->ctx
, int, qp2
->div
->n_row
);
1387 div
= isl_merge_divs(qp1
->div
, qp2
->div
, exp1
, exp2
);
1391 isl_mat_free(qp1
->div
);
1392 qp1
->div
= isl_mat_copy(div
);
1393 isl_mat_free(qp2
->div
);
1394 qp2
->div
= isl_mat_copy(div
);
1396 qp1
->upoly
= expand(qp1
->upoly
, exp1
, div
->n_col
- div
->n_row
- 2);
1397 qp2
->upoly
= expand(qp2
->upoly
, exp2
, div
->n_col
- div
->n_row
- 2);
1399 if (!qp1
->upoly
|| !qp2
->upoly
)
1406 return fn(qp1
, qp2
);
1411 isl_qpolynomial_free(qp1
);
1412 isl_qpolynomial_free(qp2
);
1416 __isl_give isl_qpolynomial
*isl_qpolynomial_add(__isl_take isl_qpolynomial
*qp1
,
1417 __isl_take isl_qpolynomial
*qp2
)
1419 qp1
= isl_qpolynomial_cow(qp1
);
1424 if (qp1
->div
->n_row
< qp2
->div
->n_row
)
1425 return isl_qpolynomial_add(qp2
, qp1
);
1427 isl_assert(qp1
->dim
->ctx
, isl_space_is_equal(qp1
->dim
, qp2
->dim
), goto error
);
1428 if (!compatible_divs(qp1
->div
, qp2
->div
))
1429 return with_merged_divs(isl_qpolynomial_add
, qp1
, qp2
);
1431 qp1
->upoly
= isl_upoly_sum(qp1
->upoly
, isl_upoly_copy(qp2
->upoly
));
1435 isl_qpolynomial_free(qp2
);
1439 isl_qpolynomial_free(qp1
);
1440 isl_qpolynomial_free(qp2
);
1444 __isl_give isl_qpolynomial
*isl_qpolynomial_add_on_domain(
1445 __isl_keep isl_set
*dom
,
1446 __isl_take isl_qpolynomial
*qp1
,
1447 __isl_take isl_qpolynomial
*qp2
)
1449 qp1
= isl_qpolynomial_add(qp1
, qp2
);
1450 qp1
= isl_qpolynomial_gist(qp1
, isl_set_copy(dom
));
1454 __isl_give isl_qpolynomial
*isl_qpolynomial_sub(__isl_take isl_qpolynomial
*qp1
,
1455 __isl_take isl_qpolynomial
*qp2
)
1457 return isl_qpolynomial_add(qp1
, isl_qpolynomial_neg(qp2
));
1460 __isl_give isl_qpolynomial
*isl_qpolynomial_add_isl_int(
1461 __isl_take isl_qpolynomial
*qp
, isl_int v
)
1463 if (isl_int_is_zero(v
))
1466 qp
= isl_qpolynomial_cow(qp
);
1470 qp
->upoly
= isl_upoly_add_isl_int(qp
->upoly
, v
);
1476 isl_qpolynomial_free(qp
);
1481 __isl_give isl_qpolynomial
*isl_qpolynomial_neg(__isl_take isl_qpolynomial
*qp
)
1486 return isl_qpolynomial_mul_isl_int(qp
, qp
->dim
->ctx
->negone
);
1489 __isl_give isl_qpolynomial
*isl_qpolynomial_mul_isl_int(
1490 __isl_take isl_qpolynomial
*qp
, isl_int v
)
1492 if (isl_int_is_one(v
))
1495 if (qp
&& isl_int_is_zero(v
)) {
1496 isl_qpolynomial
*zero
;
1497 zero
= isl_qpolynomial_zero_on_domain(isl_space_copy(qp
->dim
));
1498 isl_qpolynomial_free(qp
);
1502 qp
= isl_qpolynomial_cow(qp
);
1506 qp
->upoly
= isl_upoly_mul_isl_int(qp
->upoly
, v
);
1512 isl_qpolynomial_free(qp
);
1516 __isl_give isl_qpolynomial
*isl_qpolynomial_scale(
1517 __isl_take isl_qpolynomial
*qp
, isl_int v
)
1519 return isl_qpolynomial_mul_isl_int(qp
, v
);
1522 /* Multiply "qp" by "v".
1524 __isl_give isl_qpolynomial
*isl_qpolynomial_scale_val(
1525 __isl_take isl_qpolynomial
*qp
, __isl_take isl_val
*v
)
1530 if (!isl_val_is_rat(v
))
1531 isl_die(isl_qpolynomial_get_ctx(qp
), isl_error_invalid
,
1532 "expecting rational factor", goto error
);
1534 if (isl_val_is_one(v
)) {
1539 if (isl_val_is_zero(v
)) {
1542 space
= isl_qpolynomial_get_domain_space(qp
);
1543 isl_qpolynomial_free(qp
);
1545 return isl_qpolynomial_zero_on_domain(space
);
1548 qp
= isl_qpolynomial_cow(qp
);
1552 qp
->upoly
= isl_upoly_scale_val(qp
->upoly
, v
);
1554 qp
= isl_qpolynomial_free(qp
);
1560 isl_qpolynomial_free(qp
);
1564 __isl_give isl_qpolynomial
*isl_qpolynomial_mul(__isl_take isl_qpolynomial
*qp1
,
1565 __isl_take isl_qpolynomial
*qp2
)
1567 qp1
= isl_qpolynomial_cow(qp1
);
1572 if (qp1
->div
->n_row
< qp2
->div
->n_row
)
1573 return isl_qpolynomial_mul(qp2
, qp1
);
1575 isl_assert(qp1
->dim
->ctx
, isl_space_is_equal(qp1
->dim
, qp2
->dim
), goto error
);
1576 if (!compatible_divs(qp1
->div
, qp2
->div
))
1577 return with_merged_divs(isl_qpolynomial_mul
, qp1
, qp2
);
1579 qp1
->upoly
= isl_upoly_mul(qp1
->upoly
, isl_upoly_copy(qp2
->upoly
));
1583 isl_qpolynomial_free(qp2
);
1587 isl_qpolynomial_free(qp1
);
1588 isl_qpolynomial_free(qp2
);
1592 __isl_give isl_qpolynomial
*isl_qpolynomial_pow(__isl_take isl_qpolynomial
*qp
,
1595 qp
= isl_qpolynomial_cow(qp
);
1600 qp
->upoly
= isl_upoly_pow(qp
->upoly
, power
);
1606 isl_qpolynomial_free(qp
);
1610 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_pow(
1611 __isl_take isl_pw_qpolynomial
*pwqp
, unsigned power
)
1618 pwqp
= isl_pw_qpolynomial_cow(pwqp
);
1622 for (i
= 0; i
< pwqp
->n
; ++i
) {
1623 pwqp
->p
[i
].qp
= isl_qpolynomial_pow(pwqp
->p
[i
].qp
, power
);
1625 return isl_pw_qpolynomial_free(pwqp
);
1631 __isl_give isl_qpolynomial
*isl_qpolynomial_zero_on_domain(
1632 __isl_take isl_space
*dim
)
1636 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_zero(dim
->ctx
));
1639 __isl_give isl_qpolynomial
*isl_qpolynomial_one_on_domain(
1640 __isl_take isl_space
*dim
)
1644 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_one(dim
->ctx
));
1647 __isl_give isl_qpolynomial
*isl_qpolynomial_infty_on_domain(
1648 __isl_take isl_space
*dim
)
1652 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_infty(dim
->ctx
));
1655 __isl_give isl_qpolynomial
*isl_qpolynomial_neginfty_on_domain(
1656 __isl_take isl_space
*dim
)
1660 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_neginfty(dim
->ctx
));
1663 __isl_give isl_qpolynomial
*isl_qpolynomial_nan_on_domain(
1664 __isl_take isl_space
*dim
)
1668 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_nan(dim
->ctx
));
1671 __isl_give isl_qpolynomial
*isl_qpolynomial_cst_on_domain(
1672 __isl_take isl_space
*dim
,
1675 struct isl_qpolynomial
*qp
;
1676 struct isl_upoly_cst
*cst
;
1681 qp
= isl_qpolynomial_alloc(dim
, 0, isl_upoly_zero(dim
->ctx
));
1685 cst
= isl_upoly_as_cst(qp
->upoly
);
1686 isl_int_set(cst
->n
, v
);
1691 int isl_qpolynomial_is_cst(__isl_keep isl_qpolynomial
*qp
,
1692 isl_int
*n
, isl_int
*d
)
1694 struct isl_upoly_cst
*cst
;
1699 if (!isl_upoly_is_cst(qp
->upoly
))
1702 cst
= isl_upoly_as_cst(qp
->upoly
);
1707 isl_int_set(*n
, cst
->n
);
1709 isl_int_set(*d
, cst
->d
);
1714 /* Return the constant term of "up".
1716 static __isl_give isl_val
*isl_upoly_get_constant_val(
1717 __isl_keep
struct isl_upoly
*up
)
1719 struct isl_upoly_cst
*cst
;
1724 while (!isl_upoly_is_cst(up
)) {
1725 struct isl_upoly_rec
*rec
;
1727 rec
= isl_upoly_as_rec(up
);
1733 cst
= isl_upoly_as_cst(up
);
1736 return isl_val_rat_from_isl_int(cst
->up
.ctx
, cst
->n
, cst
->d
);
1739 /* Return the constant term of "qp".
1741 __isl_give isl_val
*isl_qpolynomial_get_constant_val(
1742 __isl_keep isl_qpolynomial
*qp
)
1747 return isl_upoly_get_constant_val(qp
->upoly
);
1750 int isl_upoly_is_affine(__isl_keep
struct isl_upoly
*up
)
1753 struct isl_upoly_rec
*rec
;
1761 rec
= isl_upoly_as_rec(up
);
1768 isl_assert(up
->ctx
, rec
->n
> 1, return -1);
1770 is_cst
= isl_upoly_is_cst(rec
->p
[1]);
1776 return isl_upoly_is_affine(rec
->p
[0]);
1779 int isl_qpolynomial_is_affine(__isl_keep isl_qpolynomial
*qp
)
1784 if (qp
->div
->n_row
> 0)
1787 return isl_upoly_is_affine(qp
->upoly
);
1790 static void update_coeff(__isl_keep isl_vec
*aff
,
1791 __isl_keep
struct isl_upoly_cst
*cst
, int pos
)
1796 if (isl_int_is_zero(cst
->n
))
1801 isl_int_gcd(gcd
, cst
->d
, aff
->el
[0]);
1802 isl_int_divexact(f
, cst
->d
, gcd
);
1803 isl_int_divexact(gcd
, aff
->el
[0], gcd
);
1804 isl_seq_scale(aff
->el
, aff
->el
, f
, aff
->size
);
1805 isl_int_mul(aff
->el
[1 + pos
], gcd
, cst
->n
);
1810 int isl_upoly_update_affine(__isl_keep
struct isl_upoly
*up
,
1811 __isl_keep isl_vec
*aff
)
1813 struct isl_upoly_cst
*cst
;
1814 struct isl_upoly_rec
*rec
;
1820 struct isl_upoly_cst
*cst
;
1822 cst
= isl_upoly_as_cst(up
);
1825 update_coeff(aff
, cst
, 0);
1829 rec
= isl_upoly_as_rec(up
);
1832 isl_assert(up
->ctx
, rec
->n
== 2, return -1);
1834 cst
= isl_upoly_as_cst(rec
->p
[1]);
1837 update_coeff(aff
, cst
, 1 + up
->var
);
1839 return isl_upoly_update_affine(rec
->p
[0], aff
);
1842 __isl_give isl_vec
*isl_qpolynomial_extract_affine(
1843 __isl_keep isl_qpolynomial
*qp
)
1851 d
= isl_space_dim(qp
->dim
, isl_dim_all
);
1852 aff
= isl_vec_alloc(qp
->div
->ctx
, 2 + d
+ qp
->div
->n_row
);
1856 isl_seq_clr(aff
->el
+ 1, 1 + d
+ qp
->div
->n_row
);
1857 isl_int_set_si(aff
->el
[0], 1);
1859 if (isl_upoly_update_affine(qp
->upoly
, aff
) < 0)
1868 int isl_qpolynomial_plain_is_equal(__isl_keep isl_qpolynomial
*qp1
,
1869 __isl_keep isl_qpolynomial
*qp2
)
1876 equal
= isl_space_is_equal(qp1
->dim
, qp2
->dim
);
1877 if (equal
< 0 || !equal
)
1880 equal
= isl_mat_is_equal(qp1
->div
, qp2
->div
);
1881 if (equal
< 0 || !equal
)
1884 return isl_upoly_is_equal(qp1
->upoly
, qp2
->upoly
);
1887 static void upoly_update_den(__isl_keep
struct isl_upoly
*up
, isl_int
*d
)
1890 struct isl_upoly_rec
*rec
;
1892 if (isl_upoly_is_cst(up
)) {
1893 struct isl_upoly_cst
*cst
;
1894 cst
= isl_upoly_as_cst(up
);
1897 isl_int_lcm(*d
, *d
, cst
->d
);
1901 rec
= isl_upoly_as_rec(up
);
1905 for (i
= 0; i
< rec
->n
; ++i
)
1906 upoly_update_den(rec
->p
[i
], d
);
1909 void isl_qpolynomial_get_den(__isl_keep isl_qpolynomial
*qp
, isl_int
*d
)
1911 isl_int_set_si(*d
, 1);
1914 upoly_update_den(qp
->upoly
, d
);
1917 __isl_give isl_qpolynomial
*isl_qpolynomial_var_pow_on_domain(
1918 __isl_take isl_space
*dim
, int pos
, int power
)
1920 struct isl_ctx
*ctx
;
1927 return isl_qpolynomial_alloc(dim
, 0, isl_upoly_var_pow(ctx
, pos
, power
));
1930 __isl_give isl_qpolynomial
*isl_qpolynomial_var_on_domain(__isl_take isl_space
*dim
,
1931 enum isl_dim_type type
, unsigned pos
)
1936 isl_assert(dim
->ctx
, isl_space_dim(dim
, isl_dim_in
) == 0, goto error
);
1937 isl_assert(dim
->ctx
, pos
< isl_space_dim(dim
, type
), goto error
);
1939 if (type
== isl_dim_set
)
1940 pos
+= isl_space_dim(dim
, isl_dim_param
);
1942 return isl_qpolynomial_var_pow_on_domain(dim
, pos
, 1);
1944 isl_space_free(dim
);
1948 __isl_give
struct isl_upoly
*isl_upoly_subs(__isl_take
struct isl_upoly
*up
,
1949 unsigned first
, unsigned n
, __isl_keep
struct isl_upoly
**subs
)
1952 struct isl_upoly_rec
*rec
;
1953 struct isl_upoly
*base
, *res
;
1958 if (isl_upoly_is_cst(up
))
1961 if (up
->var
< first
)
1964 rec
= isl_upoly_as_rec(up
);
1968 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
1970 if (up
->var
>= first
+ n
)
1971 base
= isl_upoly_var_pow(up
->ctx
, up
->var
, 1);
1973 base
= isl_upoly_copy(subs
[up
->var
- first
]);
1975 res
= isl_upoly_subs(isl_upoly_copy(rec
->p
[rec
->n
- 1]), first
, n
, subs
);
1976 for (i
= rec
->n
- 2; i
>= 0; --i
) {
1977 struct isl_upoly
*t
;
1978 t
= isl_upoly_subs(isl_upoly_copy(rec
->p
[i
]), first
, n
, subs
);
1979 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
1980 res
= isl_upoly_sum(res
, t
);
1983 isl_upoly_free(base
);
1992 __isl_give
struct isl_upoly
*isl_upoly_from_affine(isl_ctx
*ctx
, isl_int
*f
,
1993 isl_int denom
, unsigned len
)
1996 struct isl_upoly
*up
;
1998 isl_assert(ctx
, len
>= 1, return NULL
);
2000 up
= isl_upoly_rat_cst(ctx
, f
[0], denom
);
2001 for (i
= 0; i
< len
- 1; ++i
) {
2002 struct isl_upoly
*t
;
2003 struct isl_upoly
*c
;
2005 if (isl_int_is_zero(f
[1 + i
]))
2008 c
= isl_upoly_rat_cst(ctx
, f
[1 + i
], denom
);
2009 t
= isl_upoly_var_pow(ctx
, i
, 1);
2010 t
= isl_upoly_mul(c
, t
);
2011 up
= isl_upoly_sum(up
, t
);
2017 /* Remove common factor of non-constant terms and denominator.
2019 static void normalize_div(__isl_keep isl_qpolynomial
*qp
, int div
)
2021 isl_ctx
*ctx
= qp
->div
->ctx
;
2022 unsigned total
= qp
->div
->n_col
- 2;
2024 isl_seq_gcd(qp
->div
->row
[div
] + 2, total
, &ctx
->normalize_gcd
);
2025 isl_int_gcd(ctx
->normalize_gcd
,
2026 ctx
->normalize_gcd
, qp
->div
->row
[div
][0]);
2027 if (isl_int_is_one(ctx
->normalize_gcd
))
2030 isl_seq_scale_down(qp
->div
->row
[div
] + 2, qp
->div
->row
[div
] + 2,
2031 ctx
->normalize_gcd
, total
);
2032 isl_int_divexact(qp
->div
->row
[div
][0], qp
->div
->row
[div
][0],
2033 ctx
->normalize_gcd
);
2034 isl_int_fdiv_q(qp
->div
->row
[div
][1], qp
->div
->row
[div
][1],
2035 ctx
->normalize_gcd
);
2038 /* Replace the integer division identified by "div" by the polynomial "s".
2039 * The integer division is assumed not to appear in the definition
2040 * of any other integer divisions.
2042 static __isl_give isl_qpolynomial
*substitute_div(
2043 __isl_take isl_qpolynomial
*qp
,
2044 int div
, __isl_take
struct isl_upoly
*s
)
2053 qp
= isl_qpolynomial_cow(qp
);
2057 total
= isl_space_dim(qp
->dim
, isl_dim_all
);
2058 qp
->upoly
= isl_upoly_subs(qp
->upoly
, total
+ div
, 1, &s
);
2062 reordering
= isl_alloc_array(qp
->dim
->ctx
, int, total
+ qp
->div
->n_row
);
2065 for (i
= 0; i
< total
+ div
; ++i
)
2067 for (i
= total
+ div
+ 1; i
< total
+ qp
->div
->n_row
; ++i
)
2068 reordering
[i
] = i
- 1;
2069 qp
->div
= isl_mat_drop_rows(qp
->div
, div
, 1);
2070 qp
->div
= isl_mat_drop_cols(qp
->div
, 2 + total
+ div
, 1);
2071 qp
->upoly
= reorder(qp
->upoly
, reordering
);
2074 if (!qp
->upoly
|| !qp
->div
)
2080 isl_qpolynomial_free(qp
);
2085 /* Replace all integer divisions [e/d] that turn out to not actually be integer
2086 * divisions because d is equal to 1 by their definition, i.e., e.
2088 static __isl_give isl_qpolynomial
*substitute_non_divs(
2089 __isl_take isl_qpolynomial
*qp
)
2093 struct isl_upoly
*s
;
2098 total
= isl_space_dim(qp
->dim
, isl_dim_all
);
2099 for (i
= 0; qp
&& i
< qp
->div
->n_row
; ++i
) {
2100 if (!isl_int_is_one(qp
->div
->row
[i
][0]))
2102 for (j
= i
+ 1; j
< qp
->div
->n_row
; ++j
) {
2103 if (isl_int_is_zero(qp
->div
->row
[j
][2 + total
+ i
]))
2105 isl_seq_combine(qp
->div
->row
[j
] + 1,
2106 qp
->div
->ctx
->one
, qp
->div
->row
[j
] + 1,
2107 qp
->div
->row
[j
][2 + total
+ i
],
2108 qp
->div
->row
[i
] + 1, 1 + total
+ i
);
2109 isl_int_set_si(qp
->div
->row
[j
][2 + total
+ i
], 0);
2110 normalize_div(qp
, j
);
2112 s
= isl_upoly_from_affine(qp
->dim
->ctx
, qp
->div
->row
[i
] + 1,
2113 qp
->div
->row
[i
][0], qp
->div
->n_col
- 1);
2114 qp
= substitute_div(qp
, i
, s
);
2121 /* Reduce the coefficients of div "div" to lie in the interval [0, d-1],
2122 * with d the denominator. When replacing the coefficient e of x by
2123 * d * frac(e/d) = e - d * floor(e/d), we are subtracting d * floor(e/d) * x
2124 * inside the division, so we need to add floor(e/d) * x outside.
2125 * That is, we replace q by q' + floor(e/d) * x and we therefore need
2126 * to adjust the coefficient of x in each later div that depends on the
2127 * current div "div" and also in the affine expression "aff"
2128 * (if it too depends on "div").
2130 static void reduce_div(__isl_keep isl_qpolynomial
*qp
, int div
,
2131 __isl_keep isl_vec
*aff
)
2135 unsigned total
= qp
->div
->n_col
- qp
->div
->n_row
- 2;
2138 for (i
= 0; i
< 1 + total
+ div
; ++i
) {
2139 if (isl_int_is_nonneg(qp
->div
->row
[div
][1 + i
]) &&
2140 isl_int_lt(qp
->div
->row
[div
][1 + i
], qp
->div
->row
[div
][0]))
2142 isl_int_fdiv_q(v
, qp
->div
->row
[div
][1 + i
], qp
->div
->row
[div
][0]);
2143 isl_int_fdiv_r(qp
->div
->row
[div
][1 + i
],
2144 qp
->div
->row
[div
][1 + i
], qp
->div
->row
[div
][0]);
2145 if (!isl_int_is_zero(aff
->el
[1 + total
+ div
]))
2146 isl_int_addmul(aff
->el
[i
], v
, aff
->el
[1 + total
+ div
]);
2147 for (j
= div
+ 1; j
< qp
->div
->n_row
; ++j
) {
2148 if (isl_int_is_zero(qp
->div
->row
[j
][2 + total
+ div
]))
2150 isl_int_addmul(qp
->div
->row
[j
][1 + i
],
2151 v
, qp
->div
->row
[j
][2 + total
+ div
]);
2157 /* Check if the last non-zero coefficient is bigger that half of the
2158 * denominator. If so, we will invert the div to further reduce the number
2159 * of distinct divs that may appear.
2160 * If the last non-zero coefficient is exactly half the denominator,
2161 * then we continue looking for earlier coefficients that are bigger
2162 * than half the denominator.
2164 static int needs_invert(__isl_keep isl_mat
*div
, int row
)
2169 for (i
= div
->n_col
- 1; i
>= 1; --i
) {
2170 if (isl_int_is_zero(div
->row
[row
][i
]))
2172 isl_int_mul_ui(div
->row
[row
][i
], div
->row
[row
][i
], 2);
2173 cmp
= isl_int_cmp(div
->row
[row
][i
], div
->row
[row
][0]);
2174 isl_int_divexact_ui(div
->row
[row
][i
], div
->row
[row
][i
], 2);
2184 /* Replace div "div" q = [e/d] by -[(-e+(d-1))/d].
2185 * We only invert the coefficients of e (and the coefficient of q in
2186 * later divs and in "aff"). After calling this function, the
2187 * coefficients of e should be reduced again.
2189 static void invert_div(__isl_keep isl_qpolynomial
*qp
, int div
,
2190 __isl_keep isl_vec
*aff
)
2192 unsigned total
= qp
->div
->n_col
- qp
->div
->n_row
- 2;
2194 isl_seq_neg(qp
->div
->row
[div
] + 1,
2195 qp
->div
->row
[div
] + 1, qp
->div
->n_col
- 1);
2196 isl_int_sub_ui(qp
->div
->row
[div
][1], qp
->div
->row
[div
][1], 1);
2197 isl_int_add(qp
->div
->row
[div
][1],
2198 qp
->div
->row
[div
][1], qp
->div
->row
[div
][0]);
2199 if (!isl_int_is_zero(aff
->el
[1 + total
+ div
]))
2200 isl_int_neg(aff
->el
[1 + total
+ div
], aff
->el
[1 + total
+ div
]);
2201 isl_mat_col_mul(qp
->div
, 2 + total
+ div
,
2202 qp
->div
->ctx
->negone
, 2 + total
+ div
);
2205 /* Assuming "qp" is a monomial, reduce all its divs to have coefficients
2206 * in the interval [0, d-1], with d the denominator and such that the
2207 * last non-zero coefficient that is not equal to d/2 is smaller than d/2.
2209 * After the reduction, some divs may have become redundant or identical,
2210 * so we call substitute_non_divs and sort_divs. If these functions
2211 * eliminate divs or merge two or more divs into one, the coefficients
2212 * of the enclosing divs may have to be reduced again, so we call
2213 * ourselves recursively if the number of divs decreases.
2215 static __isl_give isl_qpolynomial
*reduce_divs(__isl_take isl_qpolynomial
*qp
)
2218 isl_vec
*aff
= NULL
;
2219 struct isl_upoly
*s
;
2225 aff
= isl_vec_alloc(qp
->div
->ctx
, qp
->div
->n_col
- 1);
2226 aff
= isl_vec_clr(aff
);
2230 isl_int_set_si(aff
->el
[1 + qp
->upoly
->var
], 1);
2232 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
2233 normalize_div(qp
, i
);
2234 reduce_div(qp
, i
, aff
);
2235 if (needs_invert(qp
->div
, i
)) {
2236 invert_div(qp
, i
, aff
);
2237 reduce_div(qp
, i
, aff
);
2241 s
= isl_upoly_from_affine(qp
->div
->ctx
, aff
->el
,
2242 qp
->div
->ctx
->one
, aff
->size
);
2243 qp
->upoly
= isl_upoly_subs(qp
->upoly
, qp
->upoly
->var
, 1, &s
);
2250 n_div
= qp
->div
->n_row
;
2251 qp
= substitute_non_divs(qp
);
2253 if (qp
&& qp
->div
->n_row
< n_div
)
2254 return reduce_divs(qp
);
2258 isl_qpolynomial_free(qp
);
2263 __isl_give isl_qpolynomial
*isl_qpolynomial_rat_cst_on_domain(
2264 __isl_take isl_space
*dim
, const isl_int n
, const isl_int d
)
2266 struct isl_qpolynomial
*qp
;
2267 struct isl_upoly_cst
*cst
;
2272 qp
= isl_qpolynomial_alloc(dim
, 0, isl_upoly_zero(dim
->ctx
));
2276 cst
= isl_upoly_as_cst(qp
->upoly
);
2277 isl_int_set(cst
->n
, n
);
2278 isl_int_set(cst
->d
, d
);
2283 /* Return an isl_qpolynomial that is equal to "val" on domain space "domain".
2285 __isl_give isl_qpolynomial
*isl_qpolynomial_val_on_domain(
2286 __isl_take isl_space
*domain
, __isl_take isl_val
*val
)
2288 isl_qpolynomial
*qp
;
2289 struct isl_upoly_cst
*cst
;
2291 if (!domain
|| !val
)
2294 qp
= isl_qpolynomial_alloc(domain
, 0, isl_upoly_zero(domain
->ctx
));
2298 cst
= isl_upoly_as_cst(qp
->upoly
);
2299 isl_int_set(cst
->n
, val
->n
);
2300 isl_int_set(cst
->d
, val
->d
);
2305 isl_space_free(domain
);
2310 static int up_set_active(__isl_keep
struct isl_upoly
*up
, int *active
, int d
)
2312 struct isl_upoly_rec
*rec
;
2318 if (isl_upoly_is_cst(up
))
2322 active
[up
->var
] = 1;
2324 rec
= isl_upoly_as_rec(up
);
2325 for (i
= 0; i
< rec
->n
; ++i
)
2326 if (up_set_active(rec
->p
[i
], active
, d
) < 0)
2332 static int set_active(__isl_keep isl_qpolynomial
*qp
, int *active
)
2335 int d
= isl_space_dim(qp
->dim
, isl_dim_all
);
2340 for (i
= 0; i
< d
; ++i
)
2341 for (j
= 0; j
< qp
->div
->n_row
; ++j
) {
2342 if (isl_int_is_zero(qp
->div
->row
[j
][2 + i
]))
2348 return up_set_active(qp
->upoly
, active
, d
);
2351 int isl_qpolynomial_involves_dims(__isl_keep isl_qpolynomial
*qp
,
2352 enum isl_dim_type type
, unsigned first
, unsigned n
)
2363 isl_assert(qp
->dim
->ctx
,
2364 first
+ n
<= isl_qpolynomial_dim(qp
, type
), return -1);
2365 isl_assert(qp
->dim
->ctx
, type
== isl_dim_param
||
2366 type
== isl_dim_in
, return -1);
2368 active
= isl_calloc_array(qp
->dim
->ctx
, int,
2369 isl_space_dim(qp
->dim
, isl_dim_all
));
2370 if (set_active(qp
, active
) < 0)
2373 if (type
== isl_dim_in
)
2374 first
+= isl_space_dim(qp
->dim
, isl_dim_param
);
2375 for (i
= 0; i
< n
; ++i
)
2376 if (active
[first
+ i
]) {
2389 /* Remove divs that do not appear in the quasi-polynomial, nor in any
2390 * of the divs that do appear in the quasi-polynomial.
2392 static __isl_give isl_qpolynomial
*remove_redundant_divs(
2393 __isl_take isl_qpolynomial
*qp
)
2400 int *reordering
= NULL
;
2407 if (qp
->div
->n_row
== 0)
2410 d
= isl_space_dim(qp
->dim
, isl_dim_all
);
2411 len
= qp
->div
->n_col
- 2;
2412 ctx
= isl_qpolynomial_get_ctx(qp
);
2413 active
= isl_calloc_array(ctx
, int, len
);
2417 if (up_set_active(qp
->upoly
, active
, len
) < 0)
2420 for (i
= qp
->div
->n_row
- 1; i
>= 0; --i
) {
2421 if (!active
[d
+ i
]) {
2425 for (j
= 0; j
< i
; ++j
) {
2426 if (isl_int_is_zero(qp
->div
->row
[i
][2 + d
+ j
]))
2438 reordering
= isl_alloc_array(qp
->div
->ctx
, int, len
);
2442 for (i
= 0; i
< d
; ++i
)
2446 n_div
= qp
->div
->n_row
;
2447 for (i
= 0; i
< n_div
; ++i
) {
2448 if (!active
[d
+ i
]) {
2449 qp
->div
= isl_mat_drop_rows(qp
->div
, i
- skip
, 1);
2450 qp
->div
= isl_mat_drop_cols(qp
->div
,
2451 2 + d
+ i
- skip
, 1);
2454 reordering
[d
+ i
] = d
+ i
- skip
;
2457 qp
->upoly
= reorder(qp
->upoly
, reordering
);
2459 if (!qp
->upoly
|| !qp
->div
)
2469 isl_qpolynomial_free(qp
);
2473 __isl_give
struct isl_upoly
*isl_upoly_drop(__isl_take
struct isl_upoly
*up
,
2474 unsigned first
, unsigned n
)
2477 struct isl_upoly_rec
*rec
;
2481 if (n
== 0 || up
->var
< 0 || up
->var
< first
)
2483 if (up
->var
< first
+ n
) {
2484 up
= replace_by_constant_term(up
);
2485 return isl_upoly_drop(up
, first
, n
);
2487 up
= isl_upoly_cow(up
);
2491 rec
= isl_upoly_as_rec(up
);
2495 for (i
= 0; i
< rec
->n
; ++i
) {
2496 rec
->p
[i
] = isl_upoly_drop(rec
->p
[i
], first
, n
);
2507 __isl_give isl_qpolynomial
*isl_qpolynomial_set_dim_name(
2508 __isl_take isl_qpolynomial
*qp
,
2509 enum isl_dim_type type
, unsigned pos
, const char *s
)
2511 qp
= isl_qpolynomial_cow(qp
);
2514 qp
->dim
= isl_space_set_dim_name(qp
->dim
, type
, pos
, s
);
2519 isl_qpolynomial_free(qp
);
2523 __isl_give isl_qpolynomial
*isl_qpolynomial_drop_dims(
2524 __isl_take isl_qpolynomial
*qp
,
2525 enum isl_dim_type type
, unsigned first
, unsigned n
)
2529 if (type
== isl_dim_out
)
2530 isl_die(qp
->dim
->ctx
, isl_error_invalid
,
2531 "cannot drop output/set dimension",
2533 if (type
== isl_dim_in
)
2535 if (n
== 0 && !isl_space_is_named_or_nested(qp
->dim
, type
))
2538 qp
= isl_qpolynomial_cow(qp
);
2542 isl_assert(qp
->dim
->ctx
, first
+ n
<= isl_space_dim(qp
->dim
, type
),
2544 isl_assert(qp
->dim
->ctx
, type
== isl_dim_param
||
2545 type
== isl_dim_set
, goto error
);
2547 qp
->dim
= isl_space_drop_dims(qp
->dim
, type
, first
, n
);
2551 if (type
== isl_dim_set
)
2552 first
+= isl_space_dim(qp
->dim
, isl_dim_param
);
2554 qp
->div
= isl_mat_drop_cols(qp
->div
, 2 + first
, n
);
2558 qp
->upoly
= isl_upoly_drop(qp
->upoly
, first
, n
);
2564 isl_qpolynomial_free(qp
);
2568 /* Project the domain of the quasi-polynomial onto its parameter space.
2569 * The quasi-polynomial may not involve any of the domain dimensions.
2571 __isl_give isl_qpolynomial
*isl_qpolynomial_project_domain_on_params(
2572 __isl_take isl_qpolynomial
*qp
)
2578 n
= isl_qpolynomial_dim(qp
, isl_dim_in
);
2579 involves
= isl_qpolynomial_involves_dims(qp
, isl_dim_in
, 0, n
);
2581 return isl_qpolynomial_free(qp
);
2583 isl_die(isl_qpolynomial_get_ctx(qp
), isl_error_invalid
,
2584 "polynomial involves some of the domain dimensions",
2585 return isl_qpolynomial_free(qp
));
2586 qp
= isl_qpolynomial_drop_dims(qp
, isl_dim_in
, 0, n
);
2587 space
= isl_qpolynomial_get_domain_space(qp
);
2588 space
= isl_space_params(space
);
2589 qp
= isl_qpolynomial_reset_domain_space(qp
, space
);
2593 static __isl_give isl_qpolynomial
*isl_qpolynomial_substitute_equalities_lifted(
2594 __isl_take isl_qpolynomial
*qp
, __isl_take isl_basic_set
*eq
)
2600 struct isl_upoly
*up
;
2604 if (eq
->n_eq
== 0) {
2605 isl_basic_set_free(eq
);
2609 qp
= isl_qpolynomial_cow(qp
);
2612 qp
->div
= isl_mat_cow(qp
->div
);
2616 total
= 1 + isl_space_dim(eq
->dim
, isl_dim_all
);
2618 isl_int_init(denom
);
2619 for (i
= 0; i
< eq
->n_eq
; ++i
) {
2620 j
= isl_seq_last_non_zero(eq
->eq
[i
], total
+ n_div
);
2621 if (j
< 0 || j
== 0 || j
>= total
)
2624 for (k
= 0; k
< qp
->div
->n_row
; ++k
) {
2625 if (isl_int_is_zero(qp
->div
->row
[k
][1 + j
]))
2627 isl_seq_elim(qp
->div
->row
[k
] + 1, eq
->eq
[i
], j
, total
,
2628 &qp
->div
->row
[k
][0]);
2629 normalize_div(qp
, k
);
2632 if (isl_int_is_pos(eq
->eq
[i
][j
]))
2633 isl_seq_neg(eq
->eq
[i
], eq
->eq
[i
], total
);
2634 isl_int_abs(denom
, eq
->eq
[i
][j
]);
2635 isl_int_set_si(eq
->eq
[i
][j
], 0);
2637 up
= isl_upoly_from_affine(qp
->dim
->ctx
,
2638 eq
->eq
[i
], denom
, total
);
2639 qp
->upoly
= isl_upoly_subs(qp
->upoly
, j
- 1, 1, &up
);
2642 isl_int_clear(denom
);
2647 isl_basic_set_free(eq
);
2649 qp
= substitute_non_divs(qp
);
2654 isl_basic_set_free(eq
);
2655 isl_qpolynomial_free(qp
);
2659 /* Exploit the equalities in "eq" to simplify the quasi-polynomial.
2661 __isl_give isl_qpolynomial
*isl_qpolynomial_substitute_equalities(
2662 __isl_take isl_qpolynomial
*qp
, __isl_take isl_basic_set
*eq
)
2666 if (qp
->div
->n_row
> 0)
2667 eq
= isl_basic_set_add_dims(eq
, isl_dim_set
, qp
->div
->n_row
);
2668 return isl_qpolynomial_substitute_equalities_lifted(qp
, eq
);
2670 isl_basic_set_free(eq
);
2671 isl_qpolynomial_free(qp
);
2675 static __isl_give isl_basic_set
*add_div_constraints(
2676 __isl_take isl_basic_set
*bset
, __isl_take isl_mat
*div
)
2684 bset
= isl_basic_set_extend_constraints(bset
, 0, 2 * div
->n_row
);
2687 total
= isl_basic_set_total_dim(bset
);
2688 for (i
= 0; i
< div
->n_row
; ++i
)
2689 if (isl_basic_set_add_div_constraints_var(bset
,
2690 total
- div
->n_row
+ i
, div
->row
[i
]) < 0)
2697 isl_basic_set_free(bset
);
2701 /* Look for equalities among the variables shared by context and qp
2702 * and the integer divisions of qp, if any.
2703 * The equalities are then used to eliminate variables and/or integer
2704 * divisions from qp.
2706 __isl_give isl_qpolynomial
*isl_qpolynomial_gist(
2707 __isl_take isl_qpolynomial
*qp
, __isl_take isl_set
*context
)
2713 if (qp
->div
->n_row
> 0) {
2714 isl_basic_set
*bset
;
2715 context
= isl_set_add_dims(context
, isl_dim_set
,
2717 bset
= isl_basic_set_universe(isl_set_get_space(context
));
2718 bset
= add_div_constraints(bset
, isl_mat_copy(qp
->div
));
2719 context
= isl_set_intersect(context
,
2720 isl_set_from_basic_set(bset
));
2723 aff
= isl_set_affine_hull(context
);
2724 return isl_qpolynomial_substitute_equalities_lifted(qp
, aff
);
2726 isl_qpolynomial_free(qp
);
2727 isl_set_free(context
);
2731 __isl_give isl_qpolynomial
*isl_qpolynomial_gist_params(
2732 __isl_take isl_qpolynomial
*qp
, __isl_take isl_set
*context
)
2734 isl_space
*space
= isl_qpolynomial_get_domain_space(qp
);
2735 isl_set
*dom_context
= isl_set_universe(space
);
2736 dom_context
= isl_set_intersect_params(dom_context
, context
);
2737 return isl_qpolynomial_gist(qp
, dom_context
);
2740 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_from_qpolynomial(
2741 __isl_take isl_qpolynomial
*qp
)
2747 if (isl_qpolynomial_is_zero(qp
)) {
2748 isl_space
*dim
= isl_qpolynomial_get_space(qp
);
2749 isl_qpolynomial_free(qp
);
2750 return isl_pw_qpolynomial_zero(dim
);
2753 dom
= isl_set_universe(isl_qpolynomial_get_domain_space(qp
));
2754 return isl_pw_qpolynomial_alloc(dom
, qp
);
2758 #define PW isl_pw_qpolynomial
2760 #define EL isl_qpolynomial
2762 #define EL_IS_ZERO is_zero
2766 #define IS_ZERO is_zero
2769 #undef DEFAULT_IS_ZERO
2770 #define DEFAULT_IS_ZERO 1
2774 #include <isl_pw_templ.c>
2777 #define UNION isl_union_pw_qpolynomial
2779 #define PART isl_pw_qpolynomial
2781 #define PARTS pw_qpolynomial
2782 #define ALIGN_DOMAIN
2784 #include <isl_union_templ.c>
2786 int isl_pw_qpolynomial_is_one(__isl_keep isl_pw_qpolynomial
*pwqp
)
2794 if (!isl_set_plain_is_universe(pwqp
->p
[0].set
))
2797 return isl_qpolynomial_is_one(pwqp
->p
[0].qp
);
2800 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_add(
2801 __isl_take isl_pw_qpolynomial
*pwqp1
,
2802 __isl_take isl_pw_qpolynomial
*pwqp2
)
2804 return isl_pw_qpolynomial_union_add_(pwqp1
, pwqp2
);
2807 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_mul(
2808 __isl_take isl_pw_qpolynomial
*pwqp1
,
2809 __isl_take isl_pw_qpolynomial
*pwqp2
)
2812 struct isl_pw_qpolynomial
*res
;
2814 if (!pwqp1
|| !pwqp2
)
2817 isl_assert(pwqp1
->dim
->ctx
, isl_space_is_equal(pwqp1
->dim
, pwqp2
->dim
),
2820 if (isl_pw_qpolynomial_is_zero(pwqp1
)) {
2821 isl_pw_qpolynomial_free(pwqp2
);
2825 if (isl_pw_qpolynomial_is_zero(pwqp2
)) {
2826 isl_pw_qpolynomial_free(pwqp1
);
2830 if (isl_pw_qpolynomial_is_one(pwqp1
)) {
2831 isl_pw_qpolynomial_free(pwqp1
);
2835 if (isl_pw_qpolynomial_is_one(pwqp2
)) {
2836 isl_pw_qpolynomial_free(pwqp2
);
2840 n
= pwqp1
->n
* pwqp2
->n
;
2841 res
= isl_pw_qpolynomial_alloc_size(isl_space_copy(pwqp1
->dim
), n
);
2843 for (i
= 0; i
< pwqp1
->n
; ++i
) {
2844 for (j
= 0; j
< pwqp2
->n
; ++j
) {
2845 struct isl_set
*common
;
2846 struct isl_qpolynomial
*prod
;
2847 common
= isl_set_intersect(isl_set_copy(pwqp1
->p
[i
].set
),
2848 isl_set_copy(pwqp2
->p
[j
].set
));
2849 if (isl_set_plain_is_empty(common
)) {
2850 isl_set_free(common
);
2854 prod
= isl_qpolynomial_mul(
2855 isl_qpolynomial_copy(pwqp1
->p
[i
].qp
),
2856 isl_qpolynomial_copy(pwqp2
->p
[j
].qp
));
2858 res
= isl_pw_qpolynomial_add_piece(res
, common
, prod
);
2862 isl_pw_qpolynomial_free(pwqp1
);
2863 isl_pw_qpolynomial_free(pwqp2
);
2867 isl_pw_qpolynomial_free(pwqp1
);
2868 isl_pw_qpolynomial_free(pwqp2
);
2872 __isl_give
struct isl_upoly
*isl_upoly_eval(
2873 __isl_take
struct isl_upoly
*up
, __isl_take isl_vec
*vec
)
2876 struct isl_upoly_rec
*rec
;
2877 struct isl_upoly
*res
;
2878 struct isl_upoly
*base
;
2880 if (isl_upoly_is_cst(up
)) {
2885 rec
= isl_upoly_as_rec(up
);
2889 isl_assert(up
->ctx
, rec
->n
>= 1, goto error
);
2891 base
= isl_upoly_rat_cst(up
->ctx
, vec
->el
[1 + up
->var
], vec
->el
[0]);
2893 res
= isl_upoly_eval(isl_upoly_copy(rec
->p
[rec
->n
- 1]),
2896 for (i
= rec
->n
- 2; i
>= 0; --i
) {
2897 res
= isl_upoly_mul(res
, isl_upoly_copy(base
));
2898 res
= isl_upoly_sum(res
,
2899 isl_upoly_eval(isl_upoly_copy(rec
->p
[i
]),
2900 isl_vec_copy(vec
)));
2903 isl_upoly_free(base
);
2913 __isl_give isl_qpolynomial
*isl_qpolynomial_eval(
2914 __isl_take isl_qpolynomial
*qp
, __isl_take isl_point
*pnt
)
2917 struct isl_upoly
*up
;
2922 isl_assert(pnt
->dim
->ctx
, isl_space_is_equal(pnt
->dim
, qp
->dim
), goto error
);
2924 if (qp
->div
->n_row
== 0)
2925 ext
= isl_vec_copy(pnt
->vec
);
2928 unsigned dim
= isl_space_dim(qp
->dim
, isl_dim_all
);
2929 ext
= isl_vec_alloc(qp
->dim
->ctx
, 1 + dim
+ qp
->div
->n_row
);
2933 isl_seq_cpy(ext
->el
, pnt
->vec
->el
, pnt
->vec
->size
);
2934 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
2935 isl_seq_inner_product(qp
->div
->row
[i
] + 1, ext
->el
,
2936 1 + dim
+ i
, &ext
->el
[1+dim
+i
]);
2937 isl_int_fdiv_q(ext
->el
[1+dim
+i
], ext
->el
[1+dim
+i
],
2938 qp
->div
->row
[i
][0]);
2942 up
= isl_upoly_eval(isl_upoly_copy(qp
->upoly
), ext
);
2946 dim
= isl_space_copy(qp
->dim
);
2947 isl_qpolynomial_free(qp
);
2948 isl_point_free(pnt
);
2950 return isl_qpolynomial_alloc(dim
, 0, up
);
2952 isl_qpolynomial_free(qp
);
2953 isl_point_free(pnt
);
2957 int isl_upoly_cmp(__isl_keep
struct isl_upoly_cst
*cst1
,
2958 __isl_keep
struct isl_upoly_cst
*cst2
)
2963 isl_int_mul(t
, cst1
->n
, cst2
->d
);
2964 isl_int_submul(t
, cst2
->n
, cst1
->d
);
2965 cmp
= isl_int_sgn(t
);
2970 int isl_qpolynomial_le_cst(__isl_keep isl_qpolynomial
*qp1
,
2971 __isl_keep isl_qpolynomial
*qp2
)
2973 struct isl_upoly_cst
*cst1
, *cst2
;
2977 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), return -1);
2978 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), return -1);
2979 if (isl_qpolynomial_is_nan(qp1
))
2981 if (isl_qpolynomial_is_nan(qp2
))
2983 cst1
= isl_upoly_as_cst(qp1
->upoly
);
2984 cst2
= isl_upoly_as_cst(qp2
->upoly
);
2986 return isl_upoly_cmp(cst1
, cst2
) <= 0;
2989 __isl_give isl_qpolynomial
*isl_qpolynomial_min_cst(
2990 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
2992 struct isl_upoly_cst
*cst1
, *cst2
;
2997 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), goto error
);
2998 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), goto error
);
2999 cst1
= isl_upoly_as_cst(qp1
->upoly
);
3000 cst2
= isl_upoly_as_cst(qp2
->upoly
);
3001 cmp
= isl_upoly_cmp(cst1
, cst2
);
3004 isl_qpolynomial_free(qp2
);
3006 isl_qpolynomial_free(qp1
);
3011 isl_qpolynomial_free(qp1
);
3012 isl_qpolynomial_free(qp2
);
3016 __isl_give isl_qpolynomial
*isl_qpolynomial_max_cst(
3017 __isl_take isl_qpolynomial
*qp1
, __isl_take isl_qpolynomial
*qp2
)
3019 struct isl_upoly_cst
*cst1
, *cst2
;
3024 isl_assert(qp1
->dim
->ctx
, isl_upoly_is_cst(qp1
->upoly
), goto error
);
3025 isl_assert(qp2
->dim
->ctx
, isl_upoly_is_cst(qp2
->upoly
), goto error
);
3026 cst1
= isl_upoly_as_cst(qp1
->upoly
);
3027 cst2
= isl_upoly_as_cst(qp2
->upoly
);
3028 cmp
= isl_upoly_cmp(cst1
, cst2
);
3031 isl_qpolynomial_free(qp2
);
3033 isl_qpolynomial_free(qp1
);
3038 isl_qpolynomial_free(qp1
);
3039 isl_qpolynomial_free(qp2
);
3043 __isl_give isl_qpolynomial
*isl_qpolynomial_insert_dims(
3044 __isl_take isl_qpolynomial
*qp
, enum isl_dim_type type
,
3045 unsigned first
, unsigned n
)
3053 if (type
== isl_dim_out
)
3054 isl_die(qp
->div
->ctx
, isl_error_invalid
,
3055 "cannot insert output/set dimensions",
3057 if (type
== isl_dim_in
)
3059 if (n
== 0 && !isl_space_is_named_or_nested(qp
->dim
, type
))
3062 qp
= isl_qpolynomial_cow(qp
);
3066 isl_assert(qp
->div
->ctx
, first
<= isl_space_dim(qp
->dim
, type
),
3069 g_pos
= pos(qp
->dim
, type
) + first
;
3071 qp
->div
= isl_mat_insert_zero_cols(qp
->div
, 2 + g_pos
, n
);
3075 total
= qp
->div
->n_col
- 2;
3076 if (total
> g_pos
) {
3078 exp
= isl_alloc_array(qp
->div
->ctx
, int, total
- g_pos
);
3081 for (i
= 0; i
< total
- g_pos
; ++i
)
3083 qp
->upoly
= expand(qp
->upoly
, exp
, g_pos
);
3089 qp
->dim
= isl_space_insert_dims(qp
->dim
, type
, first
, n
);
3095 isl_qpolynomial_free(qp
);
3099 __isl_give isl_qpolynomial
*isl_qpolynomial_add_dims(
3100 __isl_take isl_qpolynomial
*qp
, enum isl_dim_type type
, unsigned n
)
3104 pos
= isl_qpolynomial_dim(qp
, type
);
3106 return isl_qpolynomial_insert_dims(qp
, type
, pos
, n
);
3109 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_add_dims(
3110 __isl_take isl_pw_qpolynomial
*pwqp
,
3111 enum isl_dim_type type
, unsigned n
)
3115 pos
= isl_pw_qpolynomial_dim(pwqp
, type
);
3117 return isl_pw_qpolynomial_insert_dims(pwqp
, type
, pos
, n
);
3120 static int *reordering_move(isl_ctx
*ctx
,
3121 unsigned len
, unsigned dst
, unsigned src
, unsigned n
)
3126 reordering
= isl_alloc_array(ctx
, int, len
);
3131 for (i
= 0; i
< dst
; ++i
)
3133 for (i
= 0; i
< n
; ++i
)
3134 reordering
[src
+ i
] = dst
+ i
;
3135 for (i
= 0; i
< src
- dst
; ++i
)
3136 reordering
[dst
+ i
] = dst
+ n
+ i
;
3137 for (i
= 0; i
< len
- src
- n
; ++i
)
3138 reordering
[src
+ n
+ i
] = src
+ n
+ i
;
3140 for (i
= 0; i
< src
; ++i
)
3142 for (i
= 0; i
< n
; ++i
)
3143 reordering
[src
+ i
] = dst
+ i
;
3144 for (i
= 0; i
< dst
- src
; ++i
)
3145 reordering
[src
+ n
+ i
] = src
+ i
;
3146 for (i
= 0; i
< len
- dst
- n
; ++i
)
3147 reordering
[dst
+ n
+ i
] = dst
+ n
+ i
;
3153 __isl_give isl_qpolynomial
*isl_qpolynomial_move_dims(
3154 __isl_take isl_qpolynomial
*qp
,
3155 enum isl_dim_type dst_type
, unsigned dst_pos
,
3156 enum isl_dim_type src_type
, unsigned src_pos
, unsigned n
)
3162 qp
= isl_qpolynomial_cow(qp
);
3166 if (dst_type
== isl_dim_out
|| src_type
== isl_dim_out
)
3167 isl_die(qp
->dim
->ctx
, isl_error_invalid
,
3168 "cannot move output/set dimension",
3170 if (dst_type
== isl_dim_in
)
3171 dst_type
= isl_dim_set
;
3172 if (src_type
== isl_dim_in
)
3173 src_type
= isl_dim_set
;
3175 isl_assert(qp
->dim
->ctx
, src_pos
+ n
<= isl_space_dim(qp
->dim
, src_type
),
3178 g_dst_pos
= pos(qp
->dim
, dst_type
) + dst_pos
;
3179 g_src_pos
= pos(qp
->dim
, src_type
) + src_pos
;
3180 if (dst_type
> src_type
)
3183 qp
->div
= isl_mat_move_cols(qp
->div
, 2 + g_dst_pos
, 2 + g_src_pos
, n
);
3190 reordering
= reordering_move(qp
->dim
->ctx
,
3191 qp
->div
->n_col
- 2, g_dst_pos
, g_src_pos
, n
);
3195 qp
->upoly
= reorder(qp
->upoly
, reordering
);
3200 qp
->dim
= isl_space_move_dims(qp
->dim
, dst_type
, dst_pos
, src_type
, src_pos
, n
);
3206 isl_qpolynomial_free(qp
);
3210 __isl_give isl_qpolynomial
*isl_qpolynomial_from_affine(__isl_take isl_space
*dim
,
3211 isl_int
*f
, isl_int denom
)
3213 struct isl_upoly
*up
;
3215 dim
= isl_space_domain(dim
);
3219 up
= isl_upoly_from_affine(dim
->ctx
, f
, denom
,
3220 1 + isl_space_dim(dim
, isl_dim_all
));
3222 return isl_qpolynomial_alloc(dim
, 0, up
);
3225 __isl_give isl_qpolynomial
*isl_qpolynomial_from_aff(__isl_take isl_aff
*aff
)
3228 struct isl_upoly
*up
;
3229 isl_qpolynomial
*qp
;
3234 ctx
= isl_aff_get_ctx(aff
);
3235 up
= isl_upoly_from_affine(ctx
, aff
->v
->el
+ 1, aff
->v
->el
[0],
3238 qp
= isl_qpolynomial_alloc(isl_aff_get_domain_space(aff
),
3239 aff
->ls
->div
->n_row
, up
);
3243 isl_mat_free(qp
->div
);
3244 qp
->div
= isl_mat_copy(aff
->ls
->div
);
3245 qp
->div
= isl_mat_cow(qp
->div
);
3250 qp
= reduce_divs(qp
);
3251 qp
= remove_redundant_divs(qp
);
3258 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_from_pw_aff(
3259 __isl_take isl_pw_aff
*pwaff
)
3262 isl_pw_qpolynomial
*pwqp
;
3267 pwqp
= isl_pw_qpolynomial_alloc_size(isl_pw_aff_get_space(pwaff
),
3270 for (i
= 0; i
< pwaff
->n
; ++i
) {
3272 isl_qpolynomial
*qp
;
3274 dom
= isl_set_copy(pwaff
->p
[i
].set
);
3275 qp
= isl_qpolynomial_from_aff(isl_aff_copy(pwaff
->p
[i
].aff
));
3276 pwqp
= isl_pw_qpolynomial_add_piece(pwqp
, dom
, qp
);
3279 isl_pw_aff_free(pwaff
);
3283 __isl_give isl_qpolynomial
*isl_qpolynomial_from_constraint(
3284 __isl_take isl_constraint
*c
, enum isl_dim_type type
, unsigned pos
)
3288 aff
= isl_constraint_get_bound(c
, type
, pos
);
3289 isl_constraint_free(c
);
3290 return isl_qpolynomial_from_aff(aff
);
3293 /* For each 0 <= i < "n", replace variable "first" + i of type "type"
3294 * in "qp" by subs[i].
3296 __isl_give isl_qpolynomial
*isl_qpolynomial_substitute(
3297 __isl_take isl_qpolynomial
*qp
,
3298 enum isl_dim_type type
, unsigned first
, unsigned n
,
3299 __isl_keep isl_qpolynomial
**subs
)
3302 struct isl_upoly
**ups
;
3307 qp
= isl_qpolynomial_cow(qp
);
3311 if (type
== isl_dim_out
)
3312 isl_die(qp
->dim
->ctx
, isl_error_invalid
,
3313 "cannot substitute output/set dimension",
3315 if (type
== isl_dim_in
)
3318 for (i
= 0; i
< n
; ++i
)
3322 isl_assert(qp
->dim
->ctx
, first
+ n
<= isl_space_dim(qp
->dim
, type
),
3325 for (i
= 0; i
< n
; ++i
)
3326 isl_assert(qp
->dim
->ctx
, isl_space_is_equal(qp
->dim
, subs
[i
]->dim
),
3329 isl_assert(qp
->dim
->ctx
, qp
->div
->n_row
== 0, goto error
);
3330 for (i
= 0; i
< n
; ++i
)
3331 isl_assert(qp
->dim
->ctx
, subs
[i
]->div
->n_row
== 0, goto error
);
3333 first
+= pos(qp
->dim
, type
);
3335 ups
= isl_alloc_array(qp
->dim
->ctx
, struct isl_upoly
*, n
);
3338 for (i
= 0; i
< n
; ++i
)
3339 ups
[i
] = subs
[i
]->upoly
;
3341 qp
->upoly
= isl_upoly_subs(qp
->upoly
, first
, n
, ups
);
3350 isl_qpolynomial_free(qp
);
3354 /* Extend "bset" with extra set dimensions for each integer division
3355 * in "qp" and then call "fn" with the extended bset and the polynomial
3356 * that results from replacing each of the integer divisions by the
3357 * corresponding extra set dimension.
3359 int isl_qpolynomial_as_polynomial_on_domain(__isl_keep isl_qpolynomial
*qp
,
3360 __isl_keep isl_basic_set
*bset
,
3361 int (*fn
)(__isl_take isl_basic_set
*bset
,
3362 __isl_take isl_qpolynomial
*poly
, void *user
), void *user
)
3366 isl_qpolynomial
*poly
;
3370 if (qp
->div
->n_row
== 0)
3371 return fn(isl_basic_set_copy(bset
), isl_qpolynomial_copy(qp
),
3374 div
= isl_mat_copy(qp
->div
);
3375 dim
= isl_space_copy(qp
->dim
);
3376 dim
= isl_space_add_dims(dim
, isl_dim_set
, qp
->div
->n_row
);
3377 poly
= isl_qpolynomial_alloc(dim
, 0, isl_upoly_copy(qp
->upoly
));
3378 bset
= isl_basic_set_copy(bset
);
3379 bset
= isl_basic_set_add_dims(bset
, isl_dim_set
, qp
->div
->n_row
);
3380 bset
= add_div_constraints(bset
, div
);
3382 return fn(bset
, poly
, user
);
3387 /* Return total degree in variables first (inclusive) up to last (exclusive).
3389 int isl_upoly_degree(__isl_keep
struct isl_upoly
*up
, int first
, int last
)
3393 struct isl_upoly_rec
*rec
;
3397 if (isl_upoly_is_zero(up
))
3399 if (isl_upoly_is_cst(up
) || up
->var
< first
)
3402 rec
= isl_upoly_as_rec(up
);
3406 for (i
= 0; i
< rec
->n
; ++i
) {
3409 if (isl_upoly_is_zero(rec
->p
[i
]))
3411 d
= isl_upoly_degree(rec
->p
[i
], first
, last
);
3421 /* Return total degree in set variables.
3423 int isl_qpolynomial_degree(__isl_keep isl_qpolynomial
*poly
)
3431 ovar
= isl_space_offset(poly
->dim
, isl_dim_set
);
3432 nvar
= isl_space_dim(poly
->dim
, isl_dim_set
);
3433 return isl_upoly_degree(poly
->upoly
, ovar
, ovar
+ nvar
);
3436 __isl_give
struct isl_upoly
*isl_upoly_coeff(__isl_keep
struct isl_upoly
*up
,
3437 unsigned pos
, int deg
)
3440 struct isl_upoly_rec
*rec
;
3445 if (isl_upoly_is_cst(up
) || up
->var
< pos
) {
3447 return isl_upoly_copy(up
);
3449 return isl_upoly_zero(up
->ctx
);
3452 rec
= isl_upoly_as_rec(up
);
3456 if (up
->var
== pos
) {
3458 return isl_upoly_copy(rec
->p
[deg
]);
3460 return isl_upoly_zero(up
->ctx
);
3463 up
= isl_upoly_copy(up
);
3464 up
= isl_upoly_cow(up
);
3465 rec
= isl_upoly_as_rec(up
);
3469 for (i
= 0; i
< rec
->n
; ++i
) {
3470 struct isl_upoly
*t
;
3471 t
= isl_upoly_coeff(rec
->p
[i
], pos
, deg
);
3474 isl_upoly_free(rec
->p
[i
]);
3484 /* Return coefficient of power "deg" of variable "t_pos" of type "type".
3486 __isl_give isl_qpolynomial
*isl_qpolynomial_coeff(
3487 __isl_keep isl_qpolynomial
*qp
,
3488 enum isl_dim_type type
, unsigned t_pos
, int deg
)
3491 struct isl_upoly
*up
;
3497 if (type
== isl_dim_out
)
3498 isl_die(qp
->div
->ctx
, isl_error_invalid
,
3499 "output/set dimension does not have a coefficient",
3501 if (type
== isl_dim_in
)
3504 isl_assert(qp
->div
->ctx
, t_pos
< isl_space_dim(qp
->dim
, type
),
3507 g_pos
= pos(qp
->dim
, type
) + t_pos
;
3508 up
= isl_upoly_coeff(qp
->upoly
, g_pos
, deg
);
3510 c
= isl_qpolynomial_alloc(isl_space_copy(qp
->dim
), qp
->div
->n_row
, up
);
3513 isl_mat_free(c
->div
);
3514 c
->div
= isl_mat_copy(qp
->div
);
3519 isl_qpolynomial_free(c
);
3523 /* Homogenize the polynomial in the variables first (inclusive) up to
3524 * last (exclusive) by inserting powers of variable first.
3525 * Variable first is assumed not to appear in the input.
3527 __isl_give
struct isl_upoly
*isl_upoly_homogenize(
3528 __isl_take
struct isl_upoly
*up
, int deg
, int target
,
3529 int first
, int last
)
3532 struct isl_upoly_rec
*rec
;
3536 if (isl_upoly_is_zero(up
))
3540 if (isl_upoly_is_cst(up
) || up
->var
< first
) {
3541 struct isl_upoly
*hom
;
3543 hom
= isl_upoly_var_pow(up
->ctx
, first
, target
- deg
);
3546 rec
= isl_upoly_as_rec(hom
);
3547 rec
->p
[target
- deg
] = isl_upoly_mul(rec
->p
[target
- deg
], up
);
3552 up
= isl_upoly_cow(up
);
3553 rec
= isl_upoly_as_rec(up
);
3557 for (i
= 0; i
< rec
->n
; ++i
) {
3558 if (isl_upoly_is_zero(rec
->p
[i
]))
3560 rec
->p
[i
] = isl_upoly_homogenize(rec
->p
[i
],
3561 up
->var
< last
? deg
+ i
: i
, target
,
3573 /* Homogenize the polynomial in the set variables by introducing
3574 * powers of an extra set variable at position 0.
3576 __isl_give isl_qpolynomial
*isl_qpolynomial_homogenize(
3577 __isl_take isl_qpolynomial
*poly
)
3581 int deg
= isl_qpolynomial_degree(poly
);
3586 poly
= isl_qpolynomial_insert_dims(poly
, isl_dim_in
, 0, 1);
3587 poly
= isl_qpolynomial_cow(poly
);
3591 ovar
= isl_space_offset(poly
->dim
, isl_dim_set
);
3592 nvar
= isl_space_dim(poly
->dim
, isl_dim_set
);
3593 poly
->upoly
= isl_upoly_homogenize(poly
->upoly
, 0, deg
,
3600 isl_qpolynomial_free(poly
);
3604 __isl_give isl_term
*isl_term_alloc(__isl_take isl_space
*dim
,
3605 __isl_take isl_mat
*div
)
3613 n
= isl_space_dim(dim
, isl_dim_all
) + div
->n_row
;
3615 term
= isl_calloc(dim
->ctx
, struct isl_term
,
3616 sizeof(struct isl_term
) + (n
- 1) * sizeof(int));
3623 isl_int_init(term
->n
);
3624 isl_int_init(term
->d
);
3628 isl_space_free(dim
);
3633 __isl_give isl_term
*isl_term_copy(__isl_keep isl_term
*term
)
3642 __isl_give isl_term
*isl_term_dup(__isl_keep isl_term
*term
)
3651 total
= isl_space_dim(term
->dim
, isl_dim_all
) + term
->div
->n_row
;
3653 dup
= isl_term_alloc(isl_space_copy(term
->dim
), isl_mat_copy(term
->div
));
3657 isl_int_set(dup
->n
, term
->n
);
3658 isl_int_set(dup
->d
, term
->d
);
3660 for (i
= 0; i
< total
; ++i
)
3661 dup
->pow
[i
] = term
->pow
[i
];
3666 __isl_give isl_term
*isl_term_cow(__isl_take isl_term
*term
)
3674 return isl_term_dup(term
);
3677 void isl_term_free(__isl_take isl_term
*term
)
3682 if (--term
->ref
> 0)
3685 isl_space_free(term
->dim
);
3686 isl_mat_free(term
->div
);
3687 isl_int_clear(term
->n
);
3688 isl_int_clear(term
->d
);
3692 unsigned isl_term_dim(__isl_keep isl_term
*term
, enum isl_dim_type type
)
3700 case isl_dim_out
: return isl_space_dim(term
->dim
, type
);
3701 case isl_dim_div
: return term
->div
->n_row
;
3702 case isl_dim_all
: return isl_space_dim(term
->dim
, isl_dim_all
) +
3708 isl_ctx
*isl_term_get_ctx(__isl_keep isl_term
*term
)
3710 return term
? term
->dim
->ctx
: NULL
;
3713 void isl_term_get_num(__isl_keep isl_term
*term
, isl_int
*n
)
3717 isl_int_set(*n
, term
->n
);
3720 void isl_term_get_den(__isl_keep isl_term
*term
, isl_int
*d
)
3724 isl_int_set(*d
, term
->d
);
3727 /* Return the coefficient of the term "term".
3729 __isl_give isl_val
*isl_term_get_coefficient_val(__isl_keep isl_term
*term
)
3734 return isl_val_rat_from_isl_int(isl_term_get_ctx(term
),
3738 int isl_term_get_exp(__isl_keep isl_term
*term
,
3739 enum isl_dim_type type
, unsigned pos
)
3744 isl_assert(term
->dim
->ctx
, pos
< isl_term_dim(term
, type
), return -1);
3746 if (type
>= isl_dim_set
)
3747 pos
+= isl_space_dim(term
->dim
, isl_dim_param
);
3748 if (type
>= isl_dim_div
)
3749 pos
+= isl_space_dim(term
->dim
, isl_dim_set
);
3751 return term
->pow
[pos
];
3754 __isl_give isl_aff
*isl_term_get_div(__isl_keep isl_term
*term
, unsigned pos
)
3756 isl_local_space
*ls
;
3762 isl_assert(term
->dim
->ctx
, pos
< isl_term_dim(term
, isl_dim_div
),
3765 ls
= isl_local_space_alloc_div(isl_space_copy(term
->dim
),
3766 isl_mat_copy(term
->div
));
3767 aff
= isl_aff_alloc(ls
);
3771 isl_seq_cpy(aff
->v
->el
, term
->div
->row
[pos
], aff
->v
->size
);
3773 aff
= isl_aff_normalize(aff
);
3778 __isl_give isl_term
*isl_upoly_foreach_term(__isl_keep
struct isl_upoly
*up
,
3779 int (*fn
)(__isl_take isl_term
*term
, void *user
),
3780 __isl_take isl_term
*term
, void *user
)
3783 struct isl_upoly_rec
*rec
;
3788 if (isl_upoly_is_zero(up
))
3791 isl_assert(up
->ctx
, !isl_upoly_is_nan(up
), goto error
);
3792 isl_assert(up
->ctx
, !isl_upoly_is_infty(up
), goto error
);
3793 isl_assert(up
->ctx
, !isl_upoly_is_neginfty(up
), goto error
);
3795 if (isl_upoly_is_cst(up
)) {
3796 struct isl_upoly_cst
*cst
;
3797 cst
= isl_upoly_as_cst(up
);
3800 term
= isl_term_cow(term
);
3803 isl_int_set(term
->n
, cst
->n
);
3804 isl_int_set(term
->d
, cst
->d
);
3805 if (fn(isl_term_copy(term
), user
) < 0)
3810 rec
= isl_upoly_as_rec(up
);
3814 for (i
= 0; i
< rec
->n
; ++i
) {
3815 term
= isl_term_cow(term
);
3818 term
->pow
[up
->var
] = i
;
3819 term
= isl_upoly_foreach_term(rec
->p
[i
], fn
, term
, user
);
3823 term
->pow
[up
->var
] = 0;
3827 isl_term_free(term
);
3831 int isl_qpolynomial_foreach_term(__isl_keep isl_qpolynomial
*qp
,
3832 int (*fn
)(__isl_take isl_term
*term
, void *user
), void *user
)
3839 term
= isl_term_alloc(isl_space_copy(qp
->dim
), isl_mat_copy(qp
->div
));
3843 term
= isl_upoly_foreach_term(qp
->upoly
, fn
, term
, user
);
3845 isl_term_free(term
);
3847 return term
? 0 : -1;
3850 __isl_give isl_qpolynomial
*isl_qpolynomial_from_term(__isl_take isl_term
*term
)
3852 struct isl_upoly
*up
;
3853 isl_qpolynomial
*qp
;
3859 n
= isl_space_dim(term
->dim
, isl_dim_all
) + term
->div
->n_row
;
3861 up
= isl_upoly_rat_cst(term
->dim
->ctx
, term
->n
, term
->d
);
3862 for (i
= 0; i
< n
; ++i
) {
3865 up
= isl_upoly_mul(up
,
3866 isl_upoly_var_pow(term
->dim
->ctx
, i
, term
->pow
[i
]));
3869 qp
= isl_qpolynomial_alloc(isl_space_copy(term
->dim
), term
->div
->n_row
, up
);
3872 isl_mat_free(qp
->div
);
3873 qp
->div
= isl_mat_copy(term
->div
);
3877 isl_term_free(term
);
3880 isl_qpolynomial_free(qp
);
3881 isl_term_free(term
);
3885 __isl_give isl_qpolynomial
*isl_qpolynomial_lift(__isl_take isl_qpolynomial
*qp
,
3886 __isl_take isl_space
*dim
)
3895 if (isl_space_is_equal(qp
->dim
, dim
)) {
3896 isl_space_free(dim
);
3900 qp
= isl_qpolynomial_cow(qp
);
3904 extra
= isl_space_dim(dim
, isl_dim_set
) -
3905 isl_space_dim(qp
->dim
, isl_dim_set
);
3906 total
= isl_space_dim(qp
->dim
, isl_dim_all
);
3907 if (qp
->div
->n_row
) {
3910 exp
= isl_alloc_array(qp
->div
->ctx
, int, qp
->div
->n_row
);
3913 for (i
= 0; i
< qp
->div
->n_row
; ++i
)
3915 qp
->upoly
= expand(qp
->upoly
, exp
, total
);
3920 qp
->div
= isl_mat_insert_cols(qp
->div
, 2 + total
, extra
);
3923 for (i
= 0; i
< qp
->div
->n_row
; ++i
)
3924 isl_seq_clr(qp
->div
->row
[i
] + 2 + total
, extra
);
3926 isl_space_free(qp
->dim
);
3931 isl_space_free(dim
);
3932 isl_qpolynomial_free(qp
);
3936 /* For each parameter or variable that does not appear in qp,
3937 * first eliminate the variable from all constraints and then set it to zero.
3939 static __isl_give isl_set
*fix_inactive(__isl_take isl_set
*set
,
3940 __isl_keep isl_qpolynomial
*qp
)
3951 d
= isl_space_dim(set
->dim
, isl_dim_all
);
3952 active
= isl_calloc_array(set
->ctx
, int, d
);
3953 if (set_active(qp
, active
) < 0)
3956 for (i
= 0; i
< d
; ++i
)
3965 nparam
= isl_space_dim(set
->dim
, isl_dim_param
);
3966 nvar
= isl_space_dim(set
->dim
, isl_dim_set
);
3967 for (i
= 0; i
< nparam
; ++i
) {
3970 set
= isl_set_eliminate(set
, isl_dim_param
, i
, 1);
3971 set
= isl_set_fix_si(set
, isl_dim_param
, i
, 0);
3973 for (i
= 0; i
< nvar
; ++i
) {
3974 if (active
[nparam
+ i
])
3976 set
= isl_set_eliminate(set
, isl_dim_set
, i
, 1);
3977 set
= isl_set_fix_si(set
, isl_dim_set
, i
, 0);
3989 struct isl_opt_data
{
3990 isl_qpolynomial
*qp
;
3992 isl_qpolynomial
*opt
;
3996 static int opt_fn(__isl_take isl_point
*pnt
, void *user
)
3998 struct isl_opt_data
*data
= (struct isl_opt_data
*)user
;
3999 isl_qpolynomial
*val
;
4001 val
= isl_qpolynomial_eval(isl_qpolynomial_copy(data
->qp
), pnt
);
4005 } else if (data
->max
) {
4006 data
->opt
= isl_qpolynomial_max_cst(data
->opt
, val
);
4008 data
->opt
= isl_qpolynomial_min_cst(data
->opt
, val
);
4014 __isl_give isl_qpolynomial
*isl_qpolynomial_opt_on_domain(
4015 __isl_take isl_qpolynomial
*qp
, __isl_take isl_set
*set
, int max
)
4017 struct isl_opt_data data
= { NULL
, 1, NULL
, max
};
4022 if (isl_upoly_is_cst(qp
->upoly
)) {
4027 set
= fix_inactive(set
, qp
);
4030 if (isl_set_foreach_point(set
, opt_fn
, &data
) < 0)
4034 isl_space
*space
= isl_qpolynomial_get_domain_space(qp
);
4035 data
.opt
= isl_qpolynomial_zero_on_domain(space
);
4039 isl_qpolynomial_free(qp
);
4043 isl_qpolynomial_free(qp
);
4044 isl_qpolynomial_free(data
.opt
);
4048 __isl_give isl_qpolynomial
*isl_qpolynomial_morph_domain(
4049 __isl_take isl_qpolynomial
*qp
, __isl_take isl_morph
*morph
)
4054 struct isl_upoly
**subs
;
4055 isl_mat
*mat
, *diag
;
4057 qp
= isl_qpolynomial_cow(qp
);
4062 isl_assert(ctx
, isl_space_is_equal(qp
->dim
, morph
->dom
->dim
), goto error
);
4064 n_sub
= morph
->inv
->n_row
- 1;
4065 if (morph
->inv
->n_row
!= morph
->inv
->n_col
)
4066 n_sub
+= qp
->div
->n_row
;
4067 subs
= isl_calloc_array(ctx
, struct isl_upoly
*, n_sub
);
4071 for (i
= 0; 1 + i
< morph
->inv
->n_row
; ++i
)
4072 subs
[i
] = isl_upoly_from_affine(ctx
, morph
->inv
->row
[1 + i
],
4073 morph
->inv
->row
[0][0], morph
->inv
->n_col
);
4074 if (morph
->inv
->n_row
!= morph
->inv
->n_col
)
4075 for (i
= 0; i
< qp
->div
->n_row
; ++i
)
4076 subs
[morph
->inv
->n_row
- 1 + i
] =
4077 isl_upoly_var_pow(ctx
, morph
->inv
->n_col
- 1 + i
, 1);
4079 qp
->upoly
= isl_upoly_subs(qp
->upoly
, 0, n_sub
, subs
);
4081 for (i
= 0; i
< n_sub
; ++i
)
4082 isl_upoly_free(subs
[i
]);
4085 diag
= isl_mat_diag(ctx
, 1, morph
->inv
->row
[0][0]);
4086 mat
= isl_mat_diagonal(diag
, isl_mat_copy(morph
->inv
));
4087 diag
= isl_mat_diag(ctx
, qp
->div
->n_row
, morph
->inv
->row
[0][0]);
4088 mat
= isl_mat_diagonal(mat
, diag
);
4089 qp
->div
= isl_mat_product(qp
->div
, mat
);
4090 isl_space_free(qp
->dim
);
4091 qp
->dim
= isl_space_copy(morph
->ran
->dim
);
4093 if (!qp
->upoly
|| !qp
->div
|| !qp
->dim
)
4096 isl_morph_free(morph
);
4100 isl_qpolynomial_free(qp
);
4101 isl_morph_free(morph
);
4105 static int neg_entry(void **entry
, void *user
)
4107 isl_pw_qpolynomial
**pwqp
= (isl_pw_qpolynomial
**)entry
;
4109 *pwqp
= isl_pw_qpolynomial_neg(*pwqp
);
4111 return *pwqp
? 0 : -1;
4114 __isl_give isl_union_pw_qpolynomial
*isl_union_pw_qpolynomial_neg(
4115 __isl_take isl_union_pw_qpolynomial
*upwqp
)
4117 upwqp
= isl_union_pw_qpolynomial_cow(upwqp
);
4121 if (isl_hash_table_foreach(upwqp
->dim
->ctx
, &upwqp
->table
,
4122 &neg_entry
, NULL
) < 0)
4127 isl_union_pw_qpolynomial_free(upwqp
);
4131 __isl_give isl_union_pw_qpolynomial
*isl_union_pw_qpolynomial_mul(
4132 __isl_take isl_union_pw_qpolynomial
*upwqp1
,
4133 __isl_take isl_union_pw_qpolynomial
*upwqp2
)
4135 return match_bin_op(upwqp1
, upwqp2
, &isl_pw_qpolynomial_mul
);
4138 /* Reorder the columns of the given div definitions according to the
4141 static __isl_give isl_mat
*reorder_divs(__isl_take isl_mat
*div
,
4142 __isl_take isl_reordering
*r
)
4151 extra
= isl_space_dim(r
->dim
, isl_dim_all
) + div
->n_row
- r
->len
;
4152 mat
= isl_mat_alloc(div
->ctx
, div
->n_row
, div
->n_col
+ extra
);
4156 for (i
= 0; i
< div
->n_row
; ++i
) {
4157 isl_seq_cpy(mat
->row
[i
], div
->row
[i
], 2);
4158 isl_seq_clr(mat
->row
[i
] + 2, mat
->n_col
- 2);
4159 for (j
= 0; j
< r
->len
; ++j
)
4160 isl_int_set(mat
->row
[i
][2 + r
->pos
[j
]],
4161 div
->row
[i
][2 + j
]);
4164 isl_reordering_free(r
);
4168 isl_reordering_free(r
);
4173 /* Reorder the dimension of "qp" according to the given reordering.
4175 __isl_give isl_qpolynomial
*isl_qpolynomial_realign_domain(
4176 __isl_take isl_qpolynomial
*qp
, __isl_take isl_reordering
*r
)
4178 qp
= isl_qpolynomial_cow(qp
);
4182 r
= isl_reordering_extend(r
, qp
->div
->n_row
);
4186 qp
->div
= reorder_divs(qp
->div
, isl_reordering_copy(r
));
4190 qp
->upoly
= reorder(qp
->upoly
, r
->pos
);
4194 qp
= isl_qpolynomial_reset_domain_space(qp
, isl_space_copy(r
->dim
));
4196 isl_reordering_free(r
);
4199 isl_qpolynomial_free(qp
);
4200 isl_reordering_free(r
);
4204 __isl_give isl_qpolynomial
*isl_qpolynomial_align_params(
4205 __isl_take isl_qpolynomial
*qp
, __isl_take isl_space
*model
)
4210 if (!isl_space_match(qp
->dim
, isl_dim_param
, model
, isl_dim_param
)) {
4211 isl_reordering
*exp
;
4213 model
= isl_space_drop_dims(model
, isl_dim_in
,
4214 0, isl_space_dim(model
, isl_dim_in
));
4215 model
= isl_space_drop_dims(model
, isl_dim_out
,
4216 0, isl_space_dim(model
, isl_dim_out
));
4217 exp
= isl_parameter_alignment_reordering(qp
->dim
, model
);
4218 exp
= isl_reordering_extend_space(exp
,
4219 isl_qpolynomial_get_domain_space(qp
));
4220 qp
= isl_qpolynomial_realign_domain(qp
, exp
);
4223 isl_space_free(model
);
4226 isl_space_free(model
);
4227 isl_qpolynomial_free(qp
);
4231 struct isl_split_periods_data
{
4233 isl_pw_qpolynomial
*res
;
4236 /* Create a slice where the integer division "div" has the fixed value "v".
4237 * In particular, if "div" refers to floor(f/m), then create a slice
4239 * m v <= f <= m v + (m - 1)
4244 * -f + m v + (m - 1) >= 0
4246 static __isl_give isl_set
*set_div_slice(__isl_take isl_space
*dim
,
4247 __isl_keep isl_qpolynomial
*qp
, int div
, isl_int v
)
4250 isl_basic_set
*bset
= NULL
;
4256 total
= isl_space_dim(dim
, isl_dim_all
);
4257 bset
= isl_basic_set_alloc_space(isl_space_copy(dim
), 0, 0, 2);
4259 k
= isl_basic_set_alloc_inequality(bset
);
4262 isl_seq_cpy(bset
->ineq
[k
], qp
->div
->row
[div
] + 1, 1 + total
);
4263 isl_int_submul(bset
->ineq
[k
][0], v
, qp
->div
->row
[div
][0]);
4265 k
= isl_basic_set_alloc_inequality(bset
);
4268 isl_seq_neg(bset
->ineq
[k
], qp
->div
->row
[div
] + 1, 1 + total
);
4269 isl_int_addmul(bset
->ineq
[k
][0], v
, qp
->div
->row
[div
][0]);
4270 isl_int_add(bset
->ineq
[k
][0], bset
->ineq
[k
][0], qp
->div
->row
[div
][0]);
4271 isl_int_sub_ui(bset
->ineq
[k
][0], bset
->ineq
[k
][0], 1);
4273 isl_space_free(dim
);
4274 return isl_set_from_basic_set(bset
);
4276 isl_basic_set_free(bset
);
4277 isl_space_free(dim
);
4281 static int split_periods(__isl_take isl_set
*set
,
4282 __isl_take isl_qpolynomial
*qp
, void *user
);
4284 /* Create a slice of the domain "set" such that integer division "div"
4285 * has the fixed value "v" and add the results to data->res,
4286 * replacing the integer division by "v" in "qp".
4288 static int set_div(__isl_take isl_set
*set
,
4289 __isl_take isl_qpolynomial
*qp
, int div
, isl_int v
,
4290 struct isl_split_periods_data
*data
)
4295 struct isl_upoly
*cst
;
4297 slice
= set_div_slice(isl_set_get_space(set
), qp
, div
, v
);
4298 set
= isl_set_intersect(set
, slice
);
4303 total
= isl_space_dim(qp
->dim
, isl_dim_all
);
4305 for (i
= div
+ 1; i
< qp
->div
->n_row
; ++i
) {
4306 if (isl_int_is_zero(qp
->div
->row
[i
][2 + total
+ div
]))
4308 isl_int_addmul(qp
->div
->row
[i
][1],
4309 qp
->div
->row
[i
][2 + total
+ div
], v
);
4310 isl_int_set_si(qp
->div
->row
[i
][2 + total
+ div
], 0);
4313 cst
= isl_upoly_rat_cst(qp
->dim
->ctx
, v
, qp
->dim
->ctx
->one
);
4314 qp
= substitute_div(qp
, div
, cst
);
4316 return split_periods(set
, qp
, data
);
4319 isl_qpolynomial_free(qp
);
4323 /* Split the domain "set" such that integer division "div"
4324 * has a fixed value (ranging from "min" to "max") on each slice
4325 * and add the results to data->res.
4327 static int split_div(__isl_take isl_set
*set
,
4328 __isl_take isl_qpolynomial
*qp
, int div
, isl_int min
, isl_int max
,
4329 struct isl_split_periods_data
*data
)
4331 for (; isl_int_le(min
, max
); isl_int_add_ui(min
, min
, 1)) {
4332 isl_set
*set_i
= isl_set_copy(set
);
4333 isl_qpolynomial
*qp_i
= isl_qpolynomial_copy(qp
);
4335 if (set_div(set_i
, qp_i
, div
, min
, data
) < 0)
4339 isl_qpolynomial_free(qp
);
4343 isl_qpolynomial_free(qp
);
4347 /* If "qp" refers to any integer division
4348 * that can only attain "max_periods" distinct values on "set"
4349 * then split the domain along those distinct values.
4350 * Add the results (or the original if no splitting occurs)
4353 static int split_periods(__isl_take isl_set
*set
,
4354 __isl_take isl_qpolynomial
*qp
, void *user
)
4357 isl_pw_qpolynomial
*pwqp
;
4358 struct isl_split_periods_data
*data
;
4363 data
= (struct isl_split_periods_data
*)user
;
4368 if (qp
->div
->n_row
== 0) {
4369 pwqp
= isl_pw_qpolynomial_alloc(set
, qp
);
4370 data
->res
= isl_pw_qpolynomial_add_disjoint(data
->res
, pwqp
);
4376 total
= isl_space_dim(qp
->dim
, isl_dim_all
);
4377 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
4378 enum isl_lp_result lp_res
;
4380 if (isl_seq_first_non_zero(qp
->div
->row
[i
] + 2 + total
,
4381 qp
->div
->n_row
) != -1)
4384 lp_res
= isl_set_solve_lp(set
, 0, qp
->div
->row
[i
] + 1,
4385 set
->ctx
->one
, &min
, NULL
, NULL
);
4386 if (lp_res
== isl_lp_error
)
4388 if (lp_res
== isl_lp_unbounded
|| lp_res
== isl_lp_empty
)
4390 isl_int_fdiv_q(min
, min
, qp
->div
->row
[i
][0]);
4392 lp_res
= isl_set_solve_lp(set
, 1, qp
->div
->row
[i
] + 1,
4393 set
->ctx
->one
, &max
, NULL
, NULL
);
4394 if (lp_res
== isl_lp_error
)
4396 if (lp_res
== isl_lp_unbounded
|| lp_res
== isl_lp_empty
)
4398 isl_int_fdiv_q(max
, max
, qp
->div
->row
[i
][0]);
4400 isl_int_sub(max
, max
, min
);
4401 if (isl_int_cmp_si(max
, data
->max_periods
) < 0) {
4402 isl_int_add(max
, max
, min
);
4407 if (i
< qp
->div
->n_row
) {
4408 r
= split_div(set
, qp
, i
, min
, max
, data
);
4410 pwqp
= isl_pw_qpolynomial_alloc(set
, qp
);
4411 data
->res
= isl_pw_qpolynomial_add_disjoint(data
->res
, pwqp
);
4423 isl_qpolynomial_free(qp
);
4427 /* If any quasi-polynomial in pwqp refers to any integer division
4428 * that can only attain "max_periods" distinct values on its domain
4429 * then split the domain along those distinct values.
4431 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_split_periods(
4432 __isl_take isl_pw_qpolynomial
*pwqp
, int max_periods
)
4434 struct isl_split_periods_data data
;
4436 data
.max_periods
= max_periods
;
4437 data
.res
= isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp
));
4439 if (isl_pw_qpolynomial_foreach_piece(pwqp
, &split_periods
, &data
) < 0)
4442 isl_pw_qpolynomial_free(pwqp
);
4446 isl_pw_qpolynomial_free(data
.res
);
4447 isl_pw_qpolynomial_free(pwqp
);
4451 /* Construct a piecewise quasipolynomial that is constant on the given
4452 * domain. In particular, it is
4455 * infinity if cst == -1
4457 static __isl_give isl_pw_qpolynomial
*constant_on_domain(
4458 __isl_take isl_basic_set
*bset
, int cst
)
4461 isl_qpolynomial
*qp
;
4466 bset
= isl_basic_set_params(bset
);
4467 dim
= isl_basic_set_get_space(bset
);
4469 qp
= isl_qpolynomial_infty_on_domain(dim
);
4471 qp
= isl_qpolynomial_zero_on_domain(dim
);
4473 qp
= isl_qpolynomial_one_on_domain(dim
);
4474 return isl_pw_qpolynomial_alloc(isl_set_from_basic_set(bset
), qp
);
4477 /* Factor bset, call fn on each of the factors and return the product.
4479 * If no factors can be found, simply call fn on the input.
4480 * Otherwise, construct the factors based on the factorizer,
4481 * call fn on each factor and compute the product.
4483 static __isl_give isl_pw_qpolynomial
*compressed_multiplicative_call(
4484 __isl_take isl_basic_set
*bset
,
4485 __isl_give isl_pw_qpolynomial
*(*fn
)(__isl_take isl_basic_set
*bset
))
4491 isl_qpolynomial
*qp
;
4492 isl_pw_qpolynomial
*pwqp
;
4496 f
= isl_basic_set_factorizer(bset
);
4499 if (f
->n_group
== 0) {
4500 isl_factorizer_free(f
);
4504 nparam
= isl_basic_set_dim(bset
, isl_dim_param
);
4505 nvar
= isl_basic_set_dim(bset
, isl_dim_set
);
4507 dim
= isl_basic_set_get_space(bset
);
4508 dim
= isl_space_domain(dim
);
4509 set
= isl_set_universe(isl_space_copy(dim
));
4510 qp
= isl_qpolynomial_one_on_domain(dim
);
4511 pwqp
= isl_pw_qpolynomial_alloc(set
, qp
);
4513 bset
= isl_morph_basic_set(isl_morph_copy(f
->morph
), bset
);
4515 for (i
= 0, n
= 0; i
< f
->n_group
; ++i
) {
4516 isl_basic_set
*bset_i
;
4517 isl_pw_qpolynomial
*pwqp_i
;
4519 bset_i
= isl_basic_set_copy(bset
);
4520 bset_i
= isl_basic_set_drop_constraints_involving(bset_i
,
4521 nparam
+ n
+ f
->len
[i
], nvar
- n
- f
->len
[i
]);
4522 bset_i
= isl_basic_set_drop_constraints_involving(bset_i
,
4524 bset_i
= isl_basic_set_drop(bset_i
, isl_dim_set
,
4525 n
+ f
->len
[i
], nvar
- n
- f
->len
[i
]);
4526 bset_i
= isl_basic_set_drop(bset_i
, isl_dim_set
, 0, n
);
4528 pwqp_i
= fn(bset_i
);
4529 pwqp
= isl_pw_qpolynomial_mul(pwqp
, pwqp_i
);
4534 isl_basic_set_free(bset
);
4535 isl_factorizer_free(f
);
4539 isl_basic_set_free(bset
);
4543 /* Factor bset, call fn on each of the factors and return the product.
4544 * The function is assumed to evaluate to zero on empty domains,
4545 * to one on zero-dimensional domains and to infinity on unbounded domains
4546 * and will not be called explicitly on zero-dimensional or unbounded domains.
4548 * We first check for some special cases and remove all equalities.
4549 * Then we hand over control to compressed_multiplicative_call.
4551 __isl_give isl_pw_qpolynomial
*isl_basic_set_multiplicative_call(
4552 __isl_take isl_basic_set
*bset
,
4553 __isl_give isl_pw_qpolynomial
*(*fn
)(__isl_take isl_basic_set
*bset
))
4557 isl_pw_qpolynomial
*pwqp
;
4562 if (isl_basic_set_plain_is_empty(bset
))
4563 return constant_on_domain(bset
, 0);
4565 if (isl_basic_set_dim(bset
, isl_dim_set
) == 0)
4566 return constant_on_domain(bset
, 1);
4568 bounded
= isl_basic_set_is_bounded(bset
);
4572 return constant_on_domain(bset
, -1);
4574 if (bset
->n_eq
== 0)
4575 return compressed_multiplicative_call(bset
, fn
);
4577 morph
= isl_basic_set_full_compression(bset
);
4578 bset
= isl_morph_basic_set(isl_morph_copy(morph
), bset
);
4580 pwqp
= compressed_multiplicative_call(bset
, fn
);
4582 morph
= isl_morph_dom_params(morph
);
4583 morph
= isl_morph_ran_params(morph
);
4584 morph
= isl_morph_inverse(morph
);
4586 pwqp
= isl_pw_qpolynomial_morph_domain(pwqp
, morph
);
4590 isl_basic_set_free(bset
);
4594 /* Drop all floors in "qp", turning each integer division [a/m] into
4595 * a rational division a/m. If "down" is set, then the integer division
4596 * is replaced by (a-(m-1))/m instead.
4598 static __isl_give isl_qpolynomial
*qp_drop_floors(
4599 __isl_take isl_qpolynomial
*qp
, int down
)
4602 struct isl_upoly
*s
;
4606 if (qp
->div
->n_row
== 0)
4609 qp
= isl_qpolynomial_cow(qp
);
4613 for (i
= qp
->div
->n_row
- 1; i
>= 0; --i
) {
4615 isl_int_sub(qp
->div
->row
[i
][1],
4616 qp
->div
->row
[i
][1], qp
->div
->row
[i
][0]);
4617 isl_int_add_ui(qp
->div
->row
[i
][1],
4618 qp
->div
->row
[i
][1], 1);
4620 s
= isl_upoly_from_affine(qp
->dim
->ctx
, qp
->div
->row
[i
] + 1,
4621 qp
->div
->row
[i
][0], qp
->div
->n_col
- 1);
4622 qp
= substitute_div(qp
, i
, s
);
4630 /* Drop all floors in "pwqp", turning each integer division [a/m] into
4631 * a rational division a/m.
4633 static __isl_give isl_pw_qpolynomial
*pwqp_drop_floors(
4634 __isl_take isl_pw_qpolynomial
*pwqp
)
4641 if (isl_pw_qpolynomial_is_zero(pwqp
))
4644 pwqp
= isl_pw_qpolynomial_cow(pwqp
);
4648 for (i
= 0; i
< pwqp
->n
; ++i
) {
4649 pwqp
->p
[i
].qp
= qp_drop_floors(pwqp
->p
[i
].qp
, 0);
4656 isl_pw_qpolynomial_free(pwqp
);
4660 /* Adjust all the integer divisions in "qp" such that they are at least
4661 * one over the given orthant (identified by "signs"). This ensures
4662 * that they will still be non-negative even after subtracting (m-1)/m.
4664 * In particular, f is replaced by f' + v, changing f = [a/m]
4665 * to f' = [(a - m v)/m].
4666 * If the constant term k in a is smaller than m,
4667 * the constant term of v is set to floor(k/m) - 1.
4668 * For any other term, if the coefficient c and the variable x have
4669 * the same sign, then no changes are needed.
4670 * Otherwise, if the variable is positive (and c is negative),
4671 * then the coefficient of x in v is set to floor(c/m).
4672 * If the variable is negative (and c is positive),
4673 * then the coefficient of x in v is set to ceil(c/m).
4675 static __isl_give isl_qpolynomial
*make_divs_pos(__isl_take isl_qpolynomial
*qp
,
4681 struct isl_upoly
*s
;
4683 qp
= isl_qpolynomial_cow(qp
);
4686 qp
->div
= isl_mat_cow(qp
->div
);
4690 total
= isl_space_dim(qp
->dim
, isl_dim_all
);
4691 v
= isl_vec_alloc(qp
->div
->ctx
, qp
->div
->n_col
- 1);
4693 for (i
= 0; i
< qp
->div
->n_row
; ++i
) {
4694 isl_int
*row
= qp
->div
->row
[i
];
4698 if (isl_int_lt(row
[1], row
[0])) {
4699 isl_int_fdiv_q(v
->el
[0], row
[1], row
[0]);
4700 isl_int_sub_ui(v
->el
[0], v
->el
[0], 1);
4701 isl_int_submul(row
[1], row
[0], v
->el
[0]);
4703 for (j
= 0; j
< total
; ++j
) {
4704 if (isl_int_sgn(row
[2 + j
]) * signs
[j
] >= 0)
4707 isl_int_cdiv_q(v
->el
[1 + j
], row
[2 + j
], row
[0]);
4709 isl_int_fdiv_q(v
->el
[1 + j
], row
[2 + j
], row
[0]);
4710 isl_int_submul(row
[2 + j
], row
[0], v
->el
[1 + j
]);
4712 for (j
= 0; j
< i
; ++j
) {
4713 if (isl_int_sgn(row
[2 + total
+ j
]) >= 0)
4715 isl_int_fdiv_q(v
->el
[1 + total
+ j
],
4716 row
[2 + total
+ j
], row
[0]);
4717 isl_int_submul(row
[2 + total
+ j
],
4718 row
[0], v
->el
[1 + total
+ j
]);
4720 for (j
= i
+ 1; j
< qp
->div
->n_row
; ++j
) {
4721 if (isl_int_is_zero(qp
->div
->row
[j
][2 + total
+ i
]))
4723 isl_seq_combine(qp
->div
->row
[j
] + 1,
4724 qp
->div
->ctx
->one
, qp
->div
->row
[j
] + 1,
4725 qp
->div
->row
[j
][2 + total
+ i
], v
->el
, v
->size
);
4727 isl_int_set_si(v
->el
[1 + total
+ i
], 1);
4728 s
= isl_upoly_from_affine(qp
->dim
->ctx
, v
->el
,
4729 qp
->div
->ctx
->one
, v
->size
);
4730 qp
->upoly
= isl_upoly_subs(qp
->upoly
, total
+ i
, 1, &s
);
4740 isl_qpolynomial_free(qp
);
4744 struct isl_to_poly_data
{
4746 isl_pw_qpolynomial
*res
;
4747 isl_qpolynomial
*qp
;
4750 /* Appoximate data->qp by a polynomial on the orthant identified by "signs".
4751 * We first make all integer divisions positive and then split the
4752 * quasipolynomials into terms with sign data->sign (the direction
4753 * of the requested approximation) and terms with the opposite sign.
4754 * In the first set of terms, each integer division [a/m] is
4755 * overapproximated by a/m, while in the second it is underapproximated
4758 static int to_polynomial_on_orthant(__isl_take isl_set
*orthant
, int *signs
,
4761 struct isl_to_poly_data
*data
= user
;
4762 isl_pw_qpolynomial
*t
;
4763 isl_qpolynomial
*qp
, *up
, *down
;
4765 qp
= isl_qpolynomial_copy(data
->qp
);
4766 qp
= make_divs_pos(qp
, signs
);
4768 up
= isl_qpolynomial_terms_of_sign(qp
, signs
, data
->sign
);
4769 up
= qp_drop_floors(up
, 0);
4770 down
= isl_qpolynomial_terms_of_sign(qp
, signs
, -data
->sign
);
4771 down
= qp_drop_floors(down
, 1);
4773 isl_qpolynomial_free(qp
);
4774 qp
= isl_qpolynomial_add(up
, down
);
4776 t
= isl_pw_qpolynomial_alloc(orthant
, qp
);
4777 data
->res
= isl_pw_qpolynomial_add_disjoint(data
->res
, t
);
4782 /* Approximate each quasipolynomial by a polynomial. If "sign" is positive,
4783 * the polynomial will be an overapproximation. If "sign" is negative,
4784 * it will be an underapproximation. If "sign" is zero, the approximation
4785 * will lie somewhere in between.
4787 * In particular, is sign == 0, we simply drop the floors, turning
4788 * the integer divisions into rational divisions.
4789 * Otherwise, we split the domains into orthants, make all integer divisions
4790 * positive and then approximate each [a/m] by either a/m or (a-(m-1))/m,
4791 * depending on the requested sign and the sign of the term in which
4792 * the integer division appears.
4794 __isl_give isl_pw_qpolynomial
*isl_pw_qpolynomial_to_polynomial(
4795 __isl_take isl_pw_qpolynomial
*pwqp
, int sign
)
4798 struct isl_to_poly_data data
;
4801 return pwqp_drop_floors(pwqp
);
4807 data
.res
= isl_pw_qpolynomial_zero(isl_pw_qpolynomial_get_space(pwqp
));
4809 for (i
= 0; i
< pwqp
->n
; ++i
) {
4810 if (pwqp
->p
[i
].qp
->div
->n_row
== 0) {
4811 isl_pw_qpolynomial
*t
;
4812 t
= isl_pw_qpolynomial_alloc(
4813 isl_set_copy(pwqp
->p
[i
].set
),
4814 isl_qpolynomial_copy(pwqp
->p
[i
].qp
));
4815 data
.res
= isl_pw_qpolynomial_add_disjoint(data
.res
, t
);
4818 data
.qp
= pwqp
->p
[i
].qp
;
4819 if (isl_set_foreach_orthant(pwqp
->p
[i
].set
,
4820 &to_polynomial_on_orthant
, &data
) < 0)
4824 isl_pw_qpolynomial_free(pwqp
);
4828 isl_pw_qpolynomial_free(pwqp
);
4829 isl_pw_qpolynomial_free(data
.res
);
4833 static int poly_entry(void **entry
, void *user
)
4836 isl_pw_qpolynomial
**pwqp
= (isl_pw_qpolynomial
**)entry
;
4838 *pwqp
= isl_pw_qpolynomial_to_polynomial(*pwqp
, *sign
);
4840 return *pwqp
? 0 : -1;
4843 __isl_give isl_union_pw_qpolynomial
*isl_union_pw_qpolynomial_to_polynomial(
4844 __isl_take isl_union_pw_qpolynomial
*upwqp
, int sign
)
4846 upwqp
= isl_union_pw_qpolynomial_cow(upwqp
);
4850 if (isl_hash_table_foreach(upwqp
->dim
->ctx
, &upwqp
->table
,
4851 &poly_entry
, &sign
) < 0)
4856 isl_union_pw_qpolynomial_free(upwqp
);
4860 __isl_give isl_basic_map
*isl_basic_map_from_qpolynomial(
4861 __isl_take isl_qpolynomial
*qp
)
4865 isl_vec
*aff
= NULL
;
4866 isl_basic_map
*bmap
= NULL
;
4872 if (!isl_upoly_is_affine(qp
->upoly
))
4873 isl_die(qp
->dim
->ctx
, isl_error_invalid
,
4874 "input quasi-polynomial not affine", goto error
);
4875 aff
= isl_qpolynomial_extract_affine(qp
);
4878 dim
= isl_qpolynomial_get_space(qp
);
4879 pos
= 1 + isl_space_offset(dim
, isl_dim_out
);
4880 n_div
= qp
->div
->n_row
;
4881 bmap
= isl_basic_map_alloc_space(dim
, n_div
, 1, 2 * n_div
);
4883 for (i
= 0; i
< n_div
; ++i
) {
4884 k
= isl_basic_map_alloc_div(bmap
);
4887 isl_seq_cpy(bmap
->div
[k
], qp
->div
->row
[i
], qp
->div
->n_col
);
4888 isl_int_set_si(bmap
->div
[k
][qp
->div
->n_col
], 0);
4889 if (isl_basic_map_add_div_constraints(bmap
, k
) < 0)
4892 k
= isl_basic_map_alloc_equality(bmap
);
4895 isl_int_neg(bmap
->eq
[k
][pos
], aff
->el
[0]);
4896 isl_seq_cpy(bmap
->eq
[k
], aff
->el
+ 1, pos
);
4897 isl_seq_cpy(bmap
->eq
[k
] + pos
+ 1, aff
->el
+ 1 + pos
, n_div
);
4900 isl_qpolynomial_free(qp
);
4901 bmap
= isl_basic_map_finalize(bmap
);
4905 isl_qpolynomial_free(qp
);
4906 isl_basic_map_free(bmap
);