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