a0bdab3c48d685e7ead4e4ff5930e6bc8e5418c8
[platform/upstream/isl.git] / isl_fold.c
1 /*
2  * Copyright 2010      INRIA Saclay
3  *
4  * Use of this software is governed by the GNU LGPLv2.1 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 <isl_polynomial_private.h>
12 #include <isl_point_private.h>
13 #include <isl_dim.h>
14 #include <isl_map_private.h>
15 #include <isl_lp.h>
16 #include <isl_seq.h>
17
18 static __isl_give isl_qpolynomial_fold *qpolynomial_fold_alloc(
19         enum isl_fold type, __isl_take isl_dim *dim, int n)
20 {
21         isl_qpolynomial_fold *fold;
22
23         if (!dim)
24                 goto error;
25
26         isl_assert(dim->ctx, n >= 0, goto error);
27         fold = isl_alloc(dim->ctx, struct isl_qpolynomial_fold,
28                         sizeof(struct isl_qpolynomial_fold) +
29                         (n - 1) * sizeof(struct isl_qpolynomial *));
30         if (!fold)
31                 goto error;
32
33         fold->ref = 1;
34         fold->size = n;
35         fold->n = 0;
36         fold->type = type;
37         fold->dim = dim;
38
39         return fold;
40 error:
41         isl_dim_free(dim);
42         return NULL;
43 }
44
45 int isl_qpolynomial_fold_involves_dims(__isl_keep isl_qpolynomial_fold *fold,
46         enum isl_dim_type type, unsigned first, unsigned n)
47 {
48         int i;
49
50         if (!fold)
51                 return -1;
52         if (fold->n == 0 || n == 0)
53                 return 0;
54
55         for (i = 0; i < fold->n; ++i) {
56                 int involves = isl_qpolynomial_involves_dims(fold->qp[i],
57                                                             type, first, n);
58                 if (involves < 0 || involves)
59                         return involves;
60         }
61         return 0;
62 }
63
64 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_drop_dims(
65         __isl_take isl_qpolynomial_fold *fold,
66         enum isl_dim_type type, unsigned first, unsigned n)
67 {
68         int i;
69
70         if (!fold)
71                 return NULL;
72         if (n == 0)
73                 return fold;
74
75         fold = isl_qpolynomial_fold_cow(fold);
76         if (!fold)
77                 return NULL;
78         fold->dim = isl_dim_drop(fold->dim, type, first, n);
79         if (!fold->dim)
80                 goto error;
81
82         for (i = 0; i < fold->n; ++i) {
83                 fold->qp[i] = isl_qpolynomial_drop_dims(fold->qp[i],
84                                                             type, first, n);
85                 if (!fold->qp[i])
86                         goto error;
87         }
88
89         return fold;
90 error:
91         isl_qpolynomial_fold_free(fold);
92         return NULL;
93 }
94
95 static int isl_qpolynomial_cst_sign(__isl_keep isl_qpolynomial *qp)
96 {
97         struct isl_upoly_cst *cst;
98
99         cst = isl_upoly_as_cst(qp->upoly);
100         if (!cst)
101                 return 0;
102
103         return isl_int_sgn(cst->n) < 0 ? -1 : 1;
104 }
105
106 static int isl_qpolynomial_aff_sign(__isl_keep isl_set *set,
107         __isl_keep isl_qpolynomial *qp)
108 {
109         enum isl_lp_result res;
110         isl_vec *aff;
111         isl_int opt;
112         int sgn = 0;
113
114         aff = isl_qpolynomial_extract_affine(qp);
115         if (!aff)
116                 return 0;
117
118         isl_int_init(opt);
119
120         res = isl_set_solve_lp(set, 0, aff->el + 1, aff->el[0],
121                                 &opt, NULL, NULL);
122         if (res == isl_lp_error)
123                 goto done;
124         if (res == isl_lp_empty ||
125             (res == isl_lp_ok && !isl_int_is_neg(opt))) {
126                 sgn = 1;
127                 goto done;
128         }
129
130         res = isl_set_solve_lp(set, 1, aff->el + 1, aff->el[0],
131                                 &opt, NULL, NULL);
132         if (res == isl_lp_ok && !isl_int_is_pos(opt))
133                 sgn = -1;
134
135 done:
136         isl_int_clear(opt);
137         isl_vec_free(aff);
138         return sgn;
139 }
140
141 /* Determine, if possible, the sign of the quasipolynomial "qp" on
142  * the domain "set".
143  *
144  * If qp is a constant, then the problem is trivial.
145  * If qp is linear, then we check if the minimum of the corresponding
146  * affine constraint is non-negative or if the maximum is non-positive.
147  *
148  * Otherwise, we check if the outermost variable "v" has a lower bound "l"
149  * in "set".  If so, we write qp(v,v') as
150  *
151  *      q(v,v') * (v - l) + r(v')
152  *
153  * if q(v,v') and r(v') have the same known sign, then the original
154  * quasipolynomial has the same sign as well.
155  *
156  * Return
157  *      -1 if qp <= 0
158  *       1 if qp >= 0
159  *       0 if unknown
160  */
161 static int isl_qpolynomial_sign(__isl_keep isl_set *set,
162         __isl_keep isl_qpolynomial *qp)
163 {
164         int d;
165         int i;
166         int is;
167         struct isl_upoly_rec *rec;
168         isl_vec *v;
169         isl_int l;
170         enum isl_lp_result res;
171         int sgn = 0;
172
173         is = isl_qpolynomial_is_cst(qp, NULL, NULL);
174         if (is < 0)
175                 return 0;
176         if (is)
177                 return isl_qpolynomial_cst_sign(qp);
178
179         is = isl_qpolynomial_is_affine(qp);
180         if (is < 0)
181                 return 0;
182         if (is)
183                 return isl_qpolynomial_aff_sign(set, qp);
184
185         if (qp->div->n_row > 0)
186                 return 0;
187
188         rec = isl_upoly_as_rec(qp->upoly);
189         if (!rec)
190                 return 0;
191
192         d = isl_dim_total(qp->dim);
193         v = isl_vec_alloc(set->ctx, 2 + d);
194         if (!v)
195                 return 0;
196
197         isl_seq_clr(v->el + 1, 1 + d);
198         isl_int_set_si(v->el[0], 1);
199         isl_int_set_si(v->el[2 + qp->upoly->var], 1);
200
201         isl_int_init(l);
202
203         res = isl_set_solve_lp(set, 0, v->el + 1, v->el[0], &l, NULL, NULL);
204         if (res == isl_lp_ok) {
205                 isl_qpolynomial *min;
206                 isl_qpolynomial *base;
207                 isl_qpolynomial *r, *q;
208                 isl_qpolynomial *t;
209
210                 min = isl_qpolynomial_cst(isl_dim_copy(qp->dim), l);
211                 base = isl_qpolynomial_pow(isl_dim_copy(qp->dim),
212                                                 qp->upoly->var, 1);
213
214                 r = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), 0,
215                                           isl_upoly_copy(rec->p[rec->n - 1]));
216                 q = isl_qpolynomial_copy(r);
217
218                 for (i = rec->n - 2; i >= 0; --i) {
219                         r = isl_qpolynomial_mul(r, isl_qpolynomial_copy(min));
220                         t = isl_qpolynomial_alloc(isl_dim_copy(qp->dim), 0,
221                                                   isl_upoly_copy(rec->p[i]));
222                         r = isl_qpolynomial_add(r, t);
223                         if (i == 0)
224                                 break;
225                         q = isl_qpolynomial_mul(q, isl_qpolynomial_copy(base));
226                         q = isl_qpolynomial_add(q, isl_qpolynomial_copy(r));
227                 }
228
229                 if (isl_qpolynomial_is_zero(q))
230                         sgn = isl_qpolynomial_sign(set, r);
231                 else if (isl_qpolynomial_is_zero(r))
232                         sgn = isl_qpolynomial_sign(set, q);
233                 else {
234                         int sgn_q, sgn_r;
235                         sgn_r = isl_qpolynomial_sign(set, r);
236                         sgn_q = isl_qpolynomial_sign(set, q);
237                         if (sgn_r == sgn_q)
238                                 sgn = sgn_r;
239                 }
240
241                 isl_qpolynomial_free(min);
242                 isl_qpolynomial_free(base);
243                 isl_qpolynomial_free(q);
244                 isl_qpolynomial_free(r);
245         }
246
247         isl_int_clear(l);
248
249         isl_vec_free(v);
250
251         return sgn;
252 }
253
254 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold_on_domain(
255         __isl_keep isl_set *set,
256         __isl_take isl_qpolynomial_fold *fold1,
257         __isl_take isl_qpolynomial_fold *fold2)
258 {
259         int i, j;
260         int n1;
261         struct isl_qpolynomial_fold *res = NULL;
262         int better;
263
264         if (!fold1 || !fold2)
265                 goto error;
266
267         isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
268         isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
269                         goto error);
270
271         better = fold1->type == isl_fold_max ? -1 : 1;
272
273         if (isl_qpolynomial_fold_is_empty(fold1)) {
274                 isl_qpolynomial_fold_free(fold1);
275                 return fold2;
276         }
277
278         if (isl_qpolynomial_fold_is_empty(fold2)) {
279                 isl_qpolynomial_fold_free(fold2);
280                 return fold1;
281         }
282
283         res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
284                                         fold1->n + fold2->n);
285         if (!res)
286                 goto error;
287
288         for (i = 0; i < fold1->n; ++i) {
289                 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
290                 if (!res->qp[res->n])
291                         goto error;
292                 res->n++;
293         }
294         n1 = res->n;
295
296         for (i = 0; i < fold2->n; ++i) {
297                 for (j = n1 - 1; j >= 0; --j) {
298                         isl_qpolynomial *d;
299                         int sgn;
300                         d = isl_qpolynomial_sub(
301                                 isl_qpolynomial_copy(res->qp[j]),
302                                 isl_qpolynomial_copy(fold2->qp[i]));
303                         sgn = isl_qpolynomial_sign(set, d);
304                         isl_qpolynomial_free(d);
305                         if (sgn == 0)
306                                 continue;
307                         if (sgn != better)
308                                 break;
309                         isl_qpolynomial_free(res->qp[j]);
310                         if (j != n1 - 1)
311                                 res->qp[j] = res->qp[n1 - 1];
312                         n1--;
313                         if (n1 != res->n - 1)
314                                 res->qp[n1] = res->qp[res->n - 1];
315                         res->n--;
316                 }
317                 if (j >= 0)
318                         continue;
319                 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
320                 if (!res->qp[res->n])
321                         goto error;
322                 res->n++;
323         }
324
325         isl_qpolynomial_fold_free(fold1);
326         isl_qpolynomial_fold_free(fold2);
327
328         return res;
329 error:
330         isl_qpolynomial_fold_free(res);
331         isl_qpolynomial_fold_free(fold1);
332         isl_qpolynomial_fold_free(fold2);
333         return NULL;
334 }
335
336 #undef PW
337 #define PW isl_pw_qpolynomial_fold
338 #undef EL
339 #define EL isl_qpolynomial_fold
340 #undef IS_ZERO
341 #define IS_ZERO is_empty
342 #undef FIELD
343 #define FIELD fold
344 #undef ADD
345 #define ADD(d,e1,e2)    isl_qpolynomial_fold_fold_on_domain(d,e1,e2) 
346
347 #include <isl_pw_templ.c>
348
349 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_empty(enum isl_fold type,
350         __isl_take isl_dim *dim)
351 {
352         return qpolynomial_fold_alloc(type, dim, 0);
353 }
354
355 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_alloc(
356         enum isl_fold type, __isl_take isl_qpolynomial *qp)
357 {
358         isl_qpolynomial_fold *fold;
359
360         if (!qp)
361                 return NULL;
362
363         fold = qpolynomial_fold_alloc(type, isl_dim_copy(qp->dim), 1);
364         if (!fold)
365                 goto error;
366
367         fold->qp[0] = qp;
368         fold->n++;
369
370         return fold;
371 error:
372         isl_qpolynomial_fold_free(fold);
373         isl_qpolynomial_free(qp);
374         return NULL;
375 }
376
377 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_copy(
378         __isl_keep isl_qpolynomial_fold *fold)
379 {
380         if (!fold)
381                 return NULL;
382
383         fold->ref++;
384         return fold;
385 }
386
387 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_dup(
388         __isl_keep isl_qpolynomial_fold *fold)
389 {
390         int i;
391         isl_qpolynomial_fold *dup;
392
393         if (!fold)
394                 return NULL;
395         dup = qpolynomial_fold_alloc(fold->type,
396                                         isl_dim_copy(fold->dim), fold->n);
397         if (!dup)
398                 return NULL;
399         
400         for (i = 0; i < fold->n; ++i) {
401                 dup->qp[i] = isl_qpolynomial_copy(fold->qp[i]);
402                 if (!dup->qp[i])
403                         goto error;
404         }
405
406         return dup;
407 error:
408         isl_qpolynomial_fold_free(dup);
409         return NULL;
410 }
411
412 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_cow(
413         __isl_take isl_qpolynomial_fold *fold)
414 {
415         if (!fold)
416                 return NULL;
417
418         if (fold->ref == 1)
419                 return fold;
420         fold->ref--;
421         return isl_qpolynomial_fold_dup(fold);
422 }
423
424 void isl_qpolynomial_fold_free(__isl_take isl_qpolynomial_fold *fold)
425 {
426         int i;
427
428         if (!fold)
429                 return;
430         if (--fold->ref > 0)
431                 return;
432
433         for (i = 0; i < fold->n; ++i)
434                 isl_qpolynomial_free(fold->qp[i]);
435         isl_dim_free(fold->dim);
436         free(fold);
437 }
438
439 int isl_qpolynomial_fold_is_empty(__isl_keep isl_qpolynomial_fold *fold)
440 {
441         if (!fold)
442                 return -1;
443
444         return fold->n == 0;
445 }
446
447 __isl_give isl_qpolynomial_fold *isl_qpolynomial_fold_fold(
448         __isl_take isl_qpolynomial_fold *fold1,
449         __isl_take isl_qpolynomial_fold *fold2)
450 {
451         int i;
452         struct isl_qpolynomial_fold *res = NULL;
453
454         if (!fold1 || !fold2)
455                 goto error;
456
457         isl_assert(fold1->dim->ctx, fold1->type == fold2->type, goto error);
458         isl_assert(fold1->dim->ctx, isl_dim_equal(fold1->dim, fold2->dim),
459                         goto error);
460
461         if (isl_qpolynomial_fold_is_empty(fold1)) {
462                 isl_qpolynomial_fold_free(fold1);
463                 return fold2;
464         }
465
466         if (isl_qpolynomial_fold_is_empty(fold2)) {
467                 isl_qpolynomial_fold_free(fold2);
468                 return fold1;
469         }
470
471         res = qpolynomial_fold_alloc(fold1->type, isl_dim_copy(fold1->dim),
472                                         fold1->n + fold2->n);
473         if (!res)
474                 goto error;
475
476         for (i = 0; i < fold1->n; ++i) {
477                 res->qp[res->n] = isl_qpolynomial_copy(fold1->qp[i]);
478                 if (!res->qp[res->n])
479                         goto error;
480                 res->n++;
481         }
482
483         for (i = 0; i < fold2->n; ++i) {
484                 res->qp[res->n] = isl_qpolynomial_copy(fold2->qp[i]);
485                 if (!res->qp[res->n])
486                         goto error;
487                 res->n++;
488         }
489
490         isl_qpolynomial_fold_free(fold1);
491         isl_qpolynomial_fold_free(fold2);
492
493         return res;
494 error:
495         isl_qpolynomial_fold_free(res);
496         isl_qpolynomial_fold_free(fold1);
497         isl_qpolynomial_fold_free(fold2);
498         return NULL;
499 }
500
501 __isl_give isl_pw_qpolynomial_fold *isl_pw_qpolynomial_fold_from_pw_qpolynomial(
502         enum isl_fold type, __isl_take isl_pw_qpolynomial *pwqp)
503 {
504         int i;
505         isl_pw_qpolynomial_fold *pwf;
506
507         if (!pwqp)
508                 return NULL;
509         
510         pwf = isl_pw_qpolynomial_fold_alloc_(isl_dim_copy(pwqp->dim), pwqp->n);
511
512         for (i = 0; i < pwqp->n; ++i)
513                 pwf = isl_pw_qpolynomial_fold_add_piece(pwf,
514                         isl_set_copy(pwqp->p[i].set),
515                         isl_qpolynomial_fold_alloc(type,
516                                 isl_qpolynomial_copy(pwqp->p[i].qp)));
517
518         isl_pw_qpolynomial_free(pwqp);
519
520         return pwf;
521 }
522
523 int isl_qpolynomial_fold_is_equal(__isl_keep isl_qpolynomial_fold *fold1,
524         __isl_keep isl_qpolynomial_fold *fold2)
525 {
526         int i;
527
528         if (!fold1 || !fold2)
529                 return -1;
530
531         if (fold1->n != fold2->n)
532                 return 0;
533
534         /* We probably want to sort the qps first... */
535         for (i = 0; i < fold1->n; ++i) {
536                 int eq = isl_qpolynomial_is_equal(fold1->qp[i], fold2->qp[i]);
537                 if (eq < 0 || !eq)
538                         return eq;
539         }
540
541         return 1;
542 }
543
544 __isl_give isl_qpolynomial *isl_qpolynomial_fold_eval(
545         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_point *pnt)
546 {
547         isl_qpolynomial *qp;
548
549         if (!fold || !pnt)
550                 goto error;
551         isl_assert(pnt->dim->ctx, isl_dim_equal(pnt->dim, fold->dim), goto error);
552         isl_assert(pnt->dim->ctx,
553                 fold->type == isl_fold_max || fold->type == isl_fold_min,
554                 goto error);
555
556         if (fold->n == 0)
557                 qp = isl_qpolynomial_zero(isl_dim_copy(fold->dim));
558         else {
559                 int i;
560                 qp = isl_qpolynomial_eval(isl_qpolynomial_copy(fold->qp[0]),
561                                                 isl_point_copy(pnt));
562                 for (i = 1; i < fold->n; ++i) {
563                         isl_qpolynomial *qp_i;
564                         qp_i = isl_qpolynomial_eval(
565                                             isl_qpolynomial_copy(fold->qp[i]),
566                                             isl_point_copy(pnt));
567                         if (fold->type == isl_fold_max)
568                                 qp = isl_qpolynomial_max_cst(qp, qp_i);
569                         else
570                                 qp = isl_qpolynomial_min_cst(qp, qp_i);
571                 }
572         }
573         isl_qpolynomial_fold_free(fold);
574         isl_point_free(pnt);
575
576         return qp;
577 error:
578         isl_qpolynomial_fold_free(fold);
579         isl_point_free(pnt);
580         return NULL;
581 }
582
583 size_t isl_pw_qpolynomial_fold_size(__isl_keep isl_pw_qpolynomial_fold *pwf)
584 {
585         int i;
586         size_t n = 0;
587
588         for (i = 0; i < pwf->n; ++i)
589                 n += pwf->p[i].fold->n;
590
591         return n;
592 }
593
594 __isl_give isl_qpolynomial *isl_qpolynomial_fold_opt_on_domain(
595         __isl_take isl_qpolynomial_fold *fold, __isl_take isl_set *set, int max)
596 {
597         int i;
598         isl_qpolynomial *opt;
599
600         if (!set || !fold)
601                 goto error;
602
603         if (fold->n == 0) {
604                 isl_dim *dim = isl_dim_copy(fold->dim);
605                 isl_set_free(set);
606                 isl_qpolynomial_fold_free(fold);
607                 return isl_qpolynomial_zero(dim);
608         }
609
610         opt = isl_qpolynomial_opt_on_domain(isl_qpolynomial_copy(fold->qp[0]),
611                                                 isl_set_copy(set), max);
612         for (i = 1; i < fold->n; ++i) {
613                 isl_qpolynomial *opt_i;
614                 opt_i = isl_qpolynomial_opt_on_domain(
615                                 isl_qpolynomial_copy(fold->qp[i]),
616                                 isl_set_copy(set), max);
617                 if (max)
618                         opt = isl_qpolynomial_max_cst(opt, opt_i);
619                 else
620                         opt = isl_qpolynomial_min_cst(opt, opt_i);
621         }
622
623         isl_set_free(set);
624         isl_qpolynomial_fold_free(fold);
625
626         return opt;
627 error:
628         isl_set_free(set);
629         isl_qpolynomial_fold_free(fold);
630         return NULL;
631 }
632
633 /* Check whether for each quasi-polynomial in "fold2" there is
634  * a quasi-polynomial in "fold1" that dominates it on "set".
635  */
636 static int qpolynomial_fold_covers_on_domain(__isl_keep isl_set *set,
637         __isl_keep isl_qpolynomial_fold *fold1,
638         __isl_keep isl_qpolynomial_fold *fold2)
639 {
640         int i, j;
641         int covers;
642
643         if (!set || !fold1 || !fold2)
644                 return -1;
645
646         covers = fold1->type == isl_fold_max ? 1 : -1;
647
648         for (i = 0; i < fold2->n; ++i) {
649                 for (j = 0; j < fold1->n; ++j) {
650                         isl_qpolynomial *d;
651                         int sgn;
652
653                         d = isl_qpolynomial_sub(
654                                 isl_qpolynomial_copy(fold1->qp[j]),
655                                 isl_qpolynomial_copy(fold2->qp[i]));
656                         sgn = isl_qpolynomial_sign(set, d);
657                         isl_qpolynomial_free(d);
658                         if (sgn == covers)
659                                 break;
660                 }
661                 if (j >= fold1->n)
662                         return 0;
663         }
664
665         return 1;
666 }
667
668 /* Check whether "pwf1" dominated "pwf2", i.e., the domain of "pwf1" contains
669  * that of "pwf2" and on each cell, the corresponding fold from pwf1 dominates
670  * that of pwf2.
671  */
672 int isl_pw_qpolynomial_fold_covers(__isl_keep isl_pw_qpolynomial_fold *pwf1,
673         __isl_keep isl_pw_qpolynomial_fold *pwf2)
674 {
675         int i, j;
676         isl_set *dom1, *dom2;
677         int is_subset;
678
679         if (!pwf1 || !pwf2)
680                 return -1;
681
682         if (pwf2->n == 0)
683                 return 1;
684         if (pwf1->n == 0)
685                 return 0;
686
687         dom1 = isl_pw_qpolynomial_fold_domain(isl_pw_qpolynomial_fold_copy(pwf1));
688         dom2 = isl_pw_qpolynomial_fold_domain(isl_pw_qpolynomial_fold_copy(pwf2));
689         is_subset = isl_set_is_subset(dom2, dom1);
690         isl_set_free(dom1);
691         isl_set_free(dom2);
692
693         if (is_subset < 0 || !is_subset)
694                 return is_subset;
695
696         for (i = 0; i < pwf2->n; ++i) {
697                 for (j = 0; j < pwf1->n; ++j) {
698                         int is_empty;
699                         isl_set *common;
700                         int covers;
701
702                         common = isl_set_intersect(isl_set_copy(pwf1->p[j].set),
703                                                    isl_set_copy(pwf2->p[i].set));
704                         is_empty = isl_set_is_empty(common);
705                         if (is_empty < 0 || is_empty) {
706                                 isl_set_free(common);
707                                 if (is_empty < 0)
708                                         return -1;
709                                 continue;
710                         }
711                         covers = qpolynomial_fold_covers_on_domain(common,
712                                         pwf1->p[j].fold, pwf2->p[i].fold);
713                         isl_set_free(common);
714                         if (covers < 0 || !covers)
715                                 return covers;
716                 }
717         }
718
719         return 1;
720 }