add isl_basic_set_dim_residue_class
[platform/upstream/isl.git] / isl_map.c
1 #include <string.h>
2 #include <strings.h>
3 #include "isl_ctx.h"
4 #include "isl_blk.h"
5 #include "isl_equalities.h"
6 #include "isl_lp.h"
7 #include "isl_seq.h"
8 #include "isl_set.h"
9 #include "isl_map.h"
10 #include "isl_map_private.h"
11 #include "isl_map_piplib.h"
12 #include "isl_sample.h"
13 #include "isl_vec.h"
14
15 static struct isl_basic_map *basic_map_init(struct isl_ctx *ctx,
16                 struct isl_basic_map *bmap,
17                 unsigned nparam, unsigned n_in, unsigned n_out, unsigned extra,
18                 unsigned n_eq, unsigned n_ineq)
19 {
20         int i;
21         size_t row_size = 1 + nparam + n_in + n_out + extra;
22
23         bmap->block = isl_blk_alloc(ctx, (n_eq + n_ineq) * row_size);
24         if (isl_blk_is_error(bmap->block)) {
25                 free(bmap);
26                 return NULL;
27         }
28
29         bmap->eq = isl_alloc_array(ctx, isl_int *, n_eq + n_ineq);
30         if (!bmap->eq) {
31                 isl_blk_free(ctx, bmap->block);
32                 free(bmap);
33                 return NULL;
34         }
35
36         if (extra == 0) {
37                 bmap->block2 = isl_blk_empty();
38                 bmap->div = NULL;
39         } else {
40                 bmap->block2 = isl_blk_alloc(ctx, extra * (1 + row_size));
41                 if (isl_blk_is_error(bmap->block2)) {
42                         free(bmap->eq);
43                         isl_blk_free(ctx, bmap->block);
44                         free(bmap);
45                         return NULL;
46                 }
47
48                 bmap->div = isl_alloc_array(ctx, isl_int *, extra);
49                 if (!bmap->div) {
50                         isl_blk_free(ctx, bmap->block2);
51                         free(bmap->eq);
52                         isl_blk_free(ctx, bmap->block);
53                         free(bmap);
54                         return NULL;
55                 }
56         }
57
58         for (i = 0; i < n_eq + n_ineq; ++i)
59                 bmap->eq[i] = bmap->block.data + i * row_size;
60
61         for (i = 0; i < extra; ++i)
62                 bmap->div[i] = bmap->block2.data + i * (1 + row_size);
63
64         bmap->ctx = ctx;
65         isl_ctx_ref(ctx);
66         bmap->ref = 1;
67         bmap->flags = 0;
68         bmap->c_size = n_eq + n_ineq;
69         bmap->ineq = bmap->eq + n_eq;
70         bmap->nparam = nparam;
71         bmap->n_in = n_in;
72         bmap->n_out = n_out;
73         bmap->extra = extra;
74         bmap->n_eq = 0;
75         bmap->n_ineq = 0;
76         bmap->n_div = 0;
77         bmap->sample = NULL;
78
79         return bmap;
80 }
81
82 struct isl_basic_set *isl_basic_set_alloc(struct isl_ctx *ctx,
83                 unsigned nparam, unsigned dim, unsigned extra,
84                 unsigned n_eq, unsigned n_ineq)
85 {
86         struct isl_basic_map *bmap;
87         bmap = isl_basic_map_alloc(ctx, nparam, 0, dim, extra, n_eq, n_ineq);
88         return (struct isl_basic_set *)bmap;
89 }
90
91 struct isl_basic_map *isl_basic_map_alloc(struct isl_ctx *ctx,
92                 unsigned nparam, unsigned in, unsigned out, unsigned extra,
93                 unsigned n_eq, unsigned n_ineq)
94 {
95         struct isl_basic_map *bmap;
96
97         bmap = isl_alloc_type(ctx, struct isl_basic_map);
98         if (!bmap)
99                 return NULL;
100
101         return basic_map_init(ctx, bmap,
102                                nparam, in, out, extra, n_eq, n_ineq);
103 }
104
105 static void dup_constraints(
106                 struct isl_basic_map *dst, struct isl_basic_map *src)
107 {
108         int i;
109         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
110
111         for (i = 0; i < src->n_eq; ++i) {
112                 int j = isl_basic_map_alloc_equality(dst);
113                 isl_seq_cpy(dst->eq[j], src->eq[i], 1+total);
114         }
115
116         for (i = 0; i < src->n_ineq; ++i) {
117                 int j = isl_basic_map_alloc_inequality(dst);
118                 isl_seq_cpy(dst->ineq[j], src->ineq[i], 1+total);
119         }
120
121         for (i = 0; i < src->n_div; ++i) {
122                 int j = isl_basic_map_alloc_div(dst);
123                 isl_seq_cpy(dst->div[j], src->div[i], 1+1+total);
124         }
125         F_SET(dst, ISL_BASIC_SET_FINAL);
126 }
127
128 struct isl_basic_map *isl_basic_map_dup(struct isl_basic_map *bmap)
129 {
130         struct isl_basic_map *dup;
131
132         if (!bmap)
133                 return NULL;
134         dup = isl_basic_map_alloc(bmap->ctx, bmap->nparam,
135                         bmap->n_in, bmap->n_out,
136                         bmap->n_div, bmap->n_eq, bmap->n_ineq);
137         if (!dup)
138                 return NULL;
139         dup->flags = bmap->flags;
140         dup_constraints(dup, bmap);
141         dup->sample = isl_vec_copy(bmap->ctx, bmap->sample);
142         return dup;
143 }
144
145 struct isl_basic_set *isl_basic_set_dup(struct isl_basic_set *bset)
146 {
147         struct isl_basic_map *dup;
148
149         dup = isl_basic_map_dup((struct isl_basic_map *)bset);
150         return (struct isl_basic_set *)dup;
151 }
152
153 struct isl_basic_set *isl_basic_set_copy(struct isl_basic_set *bset)
154 {
155         if (!bset)
156                 return NULL;
157
158         if (F_ISSET(bset, ISL_BASIC_SET_FINAL)) {
159                 bset->ref++;
160                 return bset;
161         }
162         return isl_basic_set_dup(bset);
163 }
164
165 struct isl_set *isl_set_copy(struct isl_set *set)
166 {
167         if (!set)
168                 return NULL;
169
170         set->ref++;
171         return set;
172 }
173
174 struct isl_basic_map *isl_basic_map_copy(struct isl_basic_map *bmap)
175 {
176         if (!bmap)
177                 return NULL;
178
179         if (F_ISSET(bmap, ISL_BASIC_SET_FINAL)) {
180                 bmap->ref++;
181                 return bmap;
182         }
183         return isl_basic_map_dup(bmap);
184 }
185
186 struct isl_map *isl_map_copy(struct isl_map *map)
187 {
188         if (!map)
189                 return NULL;
190
191         map->ref++;
192         return map;
193 }
194
195 void isl_basic_map_free(struct isl_basic_map *bmap)
196 {
197         if (!bmap)
198                 return;
199
200         if (--bmap->ref > 0)
201                 return;
202
203         isl_ctx_deref(bmap->ctx);
204         free(bmap->div);
205         isl_blk_free(bmap->ctx, bmap->block2);
206         free(bmap->eq);
207         isl_blk_free(bmap->ctx, bmap->block);
208         isl_vec_free(bmap->ctx, bmap->sample);
209         free(bmap);
210 }
211
212 void isl_basic_set_free(struct isl_basic_set *bset)
213 {
214         isl_basic_map_free((struct isl_basic_map *)bset);
215 }
216
217 int isl_basic_map_alloc_equality(struct isl_basic_map *bmap)
218 {
219         struct isl_ctx *ctx;
220         if (!bmap)
221                 return -1;
222         ctx = bmap->ctx;
223         isl_assert(ctx, bmap->n_eq + bmap->n_ineq < bmap->c_size, return -1);
224         isl_assert(ctx, bmap->eq + bmap->n_eq <= bmap->ineq, return -1);
225         if (bmap->eq + bmap->n_eq == bmap->ineq) {
226                 isl_int *t;
227                 int j = isl_basic_map_alloc_inequality(bmap);
228                 if (j < 0)
229                         return -1;
230                 t = bmap->ineq[0];
231                 bmap->ineq[0] = bmap->ineq[j];
232                 bmap->ineq[j] = t;
233                 bmap->n_ineq--;
234                 bmap->ineq++;
235         }
236         return bmap->n_eq++;
237 }
238
239 int isl_basic_set_alloc_equality(struct isl_basic_set *bset)
240 {
241         return isl_basic_map_alloc_equality((struct isl_basic_map *)bset);
242 }
243
244 int isl_basic_map_free_equality(struct isl_basic_map *bmap, unsigned n)
245 {
246         if (!bmap)
247                 return -1;
248         isl_assert(bmap->ctx, n <= bmap->n_eq, return -1);
249         bmap->n_eq -= n;
250         return 0;
251 }
252
253 int isl_basic_map_drop_equality(struct isl_basic_map *bmap, unsigned pos)
254 {
255         isl_int *t;
256         if (!bmap)
257                 return -1;
258         isl_assert(bmap->ctx, pos < bmap->n_eq, return -1);
259
260         if (pos != bmap->n_eq - 1) {
261                 t = bmap->eq[pos];
262                 bmap->eq[pos] = bmap->eq[bmap->n_eq - 1];
263                 bmap->eq[bmap->n_eq - 1] = t;
264         }
265         bmap->n_eq--;
266         return 0;
267 }
268
269 int isl_basic_set_drop_equality(struct isl_basic_set *bset, unsigned pos)
270 {
271         return isl_basic_map_drop_equality((struct isl_basic_map *)bset, pos);
272 }
273
274 void isl_basic_map_inequality_to_equality(
275                 struct isl_basic_map *bmap, unsigned pos)
276 {
277         isl_int *t;
278
279         t = bmap->ineq[pos];
280         bmap->ineq[pos] = bmap->ineq[0];
281         bmap->ineq[0] = bmap->eq[bmap->n_eq];
282         bmap->eq[bmap->n_eq] = t;
283         bmap->n_eq++;
284         bmap->n_ineq--;
285         bmap->ineq++;
286 }
287
288 int isl_basic_map_alloc_inequality(struct isl_basic_map *bmap)
289 {
290         struct isl_ctx *ctx;
291         if (!bmap)
292                 return -1;
293         ctx = bmap->ctx;
294         isl_assert(ctx, (bmap->ineq - bmap->eq) + bmap->n_ineq < bmap->c_size,
295                         return -1);
296         F_CLR(bmap, ISL_BASIC_MAP_NO_IMPLICIT);
297         F_CLR(bmap, ISL_BASIC_MAP_NO_REDUNDANT);
298         return bmap->n_ineq++;
299 }
300
301 int isl_basic_set_alloc_inequality(struct isl_basic_set *bset)
302 {
303         return isl_basic_map_alloc_inequality((struct isl_basic_map *)bset);
304 }
305
306 int isl_basic_map_free_inequality(struct isl_basic_map *bmap, unsigned n)
307 {
308         if (!bmap)
309                 return -1;
310         isl_assert(bmap->ctx, n <= bmap->n_ineq, return -1);
311         bmap->n_ineq -= n;
312         return 0;
313 }
314
315 int isl_basic_set_free_inequality(struct isl_basic_set *bset, unsigned n)
316 {
317         return isl_basic_map_free_inequality((struct isl_basic_map *)bset, n);
318 }
319
320 int isl_basic_map_drop_inequality(struct isl_basic_map *bmap, unsigned pos)
321 {
322         isl_int *t;
323         if (!bmap)
324                 return -1;
325         isl_assert(bmap->ctx, pos < bmap->n_ineq, return -1);
326
327         if (pos != bmap->n_ineq - 1) {
328                 t = bmap->ineq[pos];
329                 bmap->ineq[pos] = bmap->ineq[bmap->n_ineq - 1];
330                 bmap->ineq[bmap->n_ineq - 1] = t;
331         }
332         bmap->n_ineq--;
333         return 0;
334 }
335
336 int isl_basic_set_drop_inequality(struct isl_basic_set *bset, unsigned pos)
337 {
338         return isl_basic_map_drop_inequality((struct isl_basic_map *)bset, pos);
339 }
340
341 int isl_basic_map_alloc_div(struct isl_basic_map *bmap)
342 {
343         if (!bmap)
344                 return -1;
345         isl_assert(bmap->ctx, bmap->n_div < bmap->extra, return -1);
346         return bmap->n_div++;
347 }
348
349 int isl_basic_map_free_div(struct isl_basic_map *bmap, unsigned n)
350 {
351         if (!bmap)
352                 return -1;
353         isl_assert(bmap->ctx, n <= bmap->n_div, return -1);
354         bmap->n_div -= n;
355         return 0;
356 }
357
358 /* Copy constraint from src to dst, putting the vars of src at offset
359  * dim_off in dst and the divs of src at offset div_off in dst.
360  * If both sets are actually map, then dim_off applies to the input
361  * variables.
362  */
363 static void copy_constraint(struct isl_basic_map *dst_map, isl_int *dst,
364                             struct isl_basic_map *src_map, isl_int *src,
365                             unsigned in_off, unsigned out_off, unsigned div_off)
366 {
367         isl_int_set(dst[0], src[0]);
368         isl_seq_cpy(dst+1, src+1, isl_min(dst_map->nparam, src_map->nparam));
369         if (dst_map->nparam > src_map->nparam)
370                 isl_seq_clr(dst+1+src_map->nparam,
371                                 dst_map->nparam - src_map->nparam);
372         isl_seq_clr(dst+1+dst_map->nparam, in_off);
373         isl_seq_cpy(dst+1+dst_map->nparam+in_off,
374                     src+1+src_map->nparam,
375                     isl_min(dst_map->n_in-in_off, src_map->n_in));
376         if (dst_map->n_in-in_off > src_map->n_in)
377                 isl_seq_clr(dst+1+dst_map->nparam+in_off+src_map->n_in,
378                                 dst_map->n_in - in_off - src_map->n_in);
379         isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in, out_off);
380         isl_seq_cpy(dst+1+dst_map->nparam+dst_map->n_in+out_off,
381                     src+1+src_map->nparam+src_map->n_in,
382                     isl_min(dst_map->n_out-out_off, src_map->n_out));
383         if (dst_map->n_out-out_off > src_map->n_out)
384                 isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in+out_off+src_map->n_out,
385                                 dst_map->n_out - out_off - src_map->n_out);
386         isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out, div_off);
387         isl_seq_cpy(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out+div_off,
388                     src+1+src_map->nparam+src_map->n_in+src_map->n_out,
389                     isl_min(dst_map->extra-div_off, src_map->extra));
390         if (dst_map->extra-div_off > src_map->extra)
391                 isl_seq_clr(dst+1+dst_map->nparam+dst_map->n_in+dst_map->n_out+
392                                 div_off+src_map->extra,
393                                 dst_map->extra - div_off - src_map->extra);
394 }
395
396 static void copy_div(struct isl_basic_map *dst_map, isl_int *dst,
397                      struct isl_basic_map *src_map, isl_int *src,
398                      unsigned in_off, unsigned out_off, unsigned div_off)
399 {
400         isl_int_set(dst[0], src[0]);
401         copy_constraint(dst_map, dst+1, src_map, src+1, in_off, out_off, div_off);
402 }
403
404 static struct isl_basic_map *add_constraints(struct isl_basic_map *bmap1,
405                 struct isl_basic_map *bmap2, unsigned i_pos, unsigned o_pos)
406 {
407         int i;
408         unsigned div_off;
409
410         if (!bmap1 || !bmap2)
411                 goto error;
412
413         div_off = bmap1->n_div;
414
415         for (i = 0; i < bmap2->n_eq; ++i) {
416                 int i1 = isl_basic_map_alloc_equality(bmap1);
417                 if (i1 < 0)
418                         goto error;
419                 copy_constraint(bmap1, bmap1->eq[i1], bmap2, bmap2->eq[i],
420                                 i_pos, o_pos, div_off);
421         }
422
423         for (i = 0; i < bmap2->n_ineq; ++i) {
424                 int i1 = isl_basic_map_alloc_inequality(bmap1);
425                 if (i1 < 0)
426                         goto error;
427                 copy_constraint(bmap1, bmap1->ineq[i1], bmap2, bmap2->ineq[i],
428                                 i_pos, o_pos, div_off);
429         }
430
431         for (i = 0; i < bmap2->n_div; ++i) {
432                 int i1 = isl_basic_map_alloc_div(bmap1);
433                 if (i1 < 0)
434                         goto error;
435                 copy_div(bmap1, bmap1->div[i1], bmap2, bmap2->div[i],
436                          i_pos, o_pos, div_off);
437         }
438
439         isl_basic_map_free(bmap2);
440
441         return bmap1;
442
443 error:
444         isl_basic_map_free(bmap1);
445         isl_basic_map_free(bmap2);
446         return NULL;
447 }
448
449 static struct isl_basic_set *set_add_constraints(struct isl_basic_set *bset1,
450                 struct isl_basic_set *bset2, unsigned pos)
451 {
452         return (struct isl_basic_set *)
453                 add_constraints((struct isl_basic_map *)bset1,
454                                 (struct isl_basic_map *)bset2, 0, pos);
455 }
456
457 struct isl_basic_map *isl_basic_map_extend(struct isl_basic_map *base,
458                 unsigned nparam, unsigned n_in, unsigned n_out, unsigned extra,
459                 unsigned n_eq, unsigned n_ineq)
460 {
461         struct isl_basic_map *ext;
462         unsigned flags;
463         int dims_ok;
464
465         base = isl_basic_map_cow(base);
466         if (!base)
467                 goto error;
468
469         dims_ok = base && base->nparam == nparam &&
470                   base->n_in == n_in && base->n_out == n_out &&
471                   base->extra >= base->n_div + extra;
472
473         if (dims_ok && n_eq == 0 && n_ineq == 0)
474                 return base;
475
476         isl_assert(base->ctx, base->nparam <= nparam, goto error);
477         isl_assert(base->ctx, base->n_in <= n_in, goto error);
478         isl_assert(base->ctx, base->n_out <= n_out, goto error);
479         extra += base->extra;
480         n_eq += base->n_eq;
481         n_ineq += base->n_ineq;
482
483         ext = isl_basic_map_alloc(base->ctx, nparam, n_in, n_out,
484                                         extra, n_eq, n_ineq);
485         if (!base)
486                 return ext;
487         if (!ext)
488                 goto error;
489
490         flags = base->flags;
491         ext = add_constraints(ext, base, 0, 0);
492         if (ext) {
493                 ext->flags = flags;
494                 F_CLR(ext, ISL_BASIC_SET_FINAL);
495         }
496
497         return ext;
498
499 error:
500         isl_basic_map_free(base);
501         return NULL;
502 }
503
504 struct isl_basic_set *isl_basic_set_extend(struct isl_basic_set *base,
505                 unsigned nparam, unsigned dim, unsigned extra,
506                 unsigned n_eq, unsigned n_ineq)
507 {
508         return (struct isl_basic_set *)
509                 isl_basic_map_extend((struct isl_basic_map *)base,
510                                         nparam, 0, dim, extra, n_eq, n_ineq);
511 }
512
513 struct isl_basic_set *isl_basic_set_cow(struct isl_basic_set *bset)
514 {
515         return (struct isl_basic_set *)
516                 isl_basic_map_cow((struct isl_basic_map *)bset);
517 }
518
519 struct isl_basic_map *isl_basic_map_cow(struct isl_basic_map *bmap)
520 {
521         if (!bmap)
522                 return NULL;
523
524         if (bmap->ref > 1) {
525                 bmap->ref--;
526                 bmap = isl_basic_map_dup(bmap);
527         }
528         F_CLR(bmap, ISL_BASIC_SET_FINAL);
529         return bmap;
530 }
531
532 struct isl_set *isl_set_cow(struct isl_set *set)
533 {
534         if (!set)
535                 return NULL;
536
537         if (set->ref == 1)
538                 return set;
539         set->ref--;
540         return isl_set_dup(set);
541 }
542
543 struct isl_map *isl_map_cow(struct isl_map *map)
544 {
545         if (!map)
546                 return NULL;
547
548         if (map->ref == 1)
549                 return map;
550         map->ref--;
551         return isl_map_dup(map);
552 }
553
554 static void swap_vars(struct isl_blk blk, isl_int *a,
555                         unsigned a_len, unsigned b_len)
556 {
557         isl_seq_cpy(blk.data, a+a_len, b_len);
558         isl_seq_cpy(blk.data+b_len, a, a_len);
559         isl_seq_cpy(a, blk.data, b_len+a_len);
560 }
561
562 struct isl_basic_set *isl_basic_set_swap_vars(
563                 struct isl_basic_set *bset, unsigned n)
564 {
565         int i;
566         struct isl_blk blk;
567
568         if (!bset)
569                 goto error;
570
571         isl_assert(bset->ctx, n <= bset->dim, goto error);
572
573         if (n == bset->dim)
574                 return bset;
575
576         bset = isl_basic_set_cow(bset);
577         if (!bset)
578                 return NULL;
579
580         blk = isl_blk_alloc(bset->ctx, bset->dim);
581         if (isl_blk_is_error(blk))
582                 goto error;
583
584         for (i = 0; i < bset->n_eq; ++i)
585                 swap_vars(blk,
586                           bset->eq[i]+1+bset->nparam, n, bset->dim - n);
587
588         for (i = 0; i < bset->n_ineq; ++i)
589                 swap_vars(blk,
590                           bset->ineq[i]+1+bset->nparam, n, bset->dim - n);
591
592         for (i = 0; i < bset->n_div; ++i)
593                 swap_vars(blk,
594                           bset->div[i]+1+1+bset->nparam, n, bset->dim - n);
595
596         isl_blk_free(bset->ctx, blk);
597
598         return bset;
599
600 error:
601         isl_basic_set_free(bset);
602         return NULL;
603 }
604
605 struct isl_set *isl_set_swap_vars(struct isl_set *set, unsigned n)
606 {
607         int i;
608         set = isl_set_cow(set);
609         if (!set)
610                 return NULL;
611
612         for (i = 0; i < set->n; ++i) {
613                 set->p[i] = isl_basic_set_swap_vars(set->p[i], n);
614                 if (!set->p[i]) {
615                         isl_set_free(set);
616                         return NULL;
617                 }
618         }
619         return set;
620 }
621
622 static void constraint_drop_vars(isl_int *c, unsigned n, unsigned rem)
623 {
624         isl_seq_cpy(c, c + n, rem);
625         isl_seq_clr(c + rem, n);
626 }
627
628 struct isl_basic_set *isl_basic_set_drop_dims(
629                 struct isl_basic_set *bset, unsigned first, unsigned n)
630 {
631         int i;
632
633         if (!bset)
634                 goto error;
635
636         isl_assert(bset->ctx, first + n <= bset->dim, goto error);
637
638         if (n == 0)
639                 return bset;
640
641         bset = isl_basic_set_cow(bset);
642         if (!bset)
643                 return NULL;
644
645         for (i = 0; i < bset->n_eq; ++i)
646                 constraint_drop_vars(bset->eq[i]+1+bset->nparam+first, n,
647                                      (bset->dim-first-n)+bset->extra);
648
649         for (i = 0; i < bset->n_ineq; ++i)
650                 constraint_drop_vars(bset->ineq[i]+1+bset->nparam+first, n,
651                                      (bset->dim-first-n)+bset->extra);
652
653         for (i = 0; i < bset->n_div; ++i)
654                 constraint_drop_vars(bset->div[i]+1+1+bset->nparam+first, n,
655                                      (bset->dim-first-n)+bset->extra);
656
657         bset->dim -= n;
658         bset->extra += n;
659
660         return isl_basic_set_finalize(bset);
661 error:
662         isl_basic_set_free(bset);
663         return NULL;
664 }
665
666 static struct isl_set *isl_set_drop_dims(
667                 struct isl_set *set, unsigned first, unsigned n)
668 {
669         int i;
670
671         if (!set)
672                 goto error;
673
674         isl_assert(set->ctx, first + n <= set->dim, goto error);
675
676         if (n == 0)
677                 return set;
678         set = isl_set_cow(set);
679         if (!set)
680                 goto error;
681
682         for (i = 0; i < set->n; ++i) {
683                 set->p[i] = isl_basic_set_drop_dims(set->p[i], first, n);
684                 if (!set->p[i])
685                         goto error;
686         }
687         set->dim -= n;
688
689         return set;
690 error:
691         isl_set_free(set);
692         return NULL;
693 }
694
695 /*
696  * We don't cow, as the div is assumed to be redundant.
697  */
698 static struct isl_basic_map *isl_basic_map_drop_div(
699                 struct isl_basic_map *bmap, unsigned div)
700 {
701         int i;
702         unsigned pos;
703
704         if (!bmap)
705                 goto error;
706
707         pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
708
709         isl_assert(bmap->ctx, div < bmap->n_div, goto error);
710
711         for (i = 0; i < bmap->n_eq; ++i)
712                 constraint_drop_vars(bmap->eq[i]+pos, 1, bmap->extra-div-1);
713
714         for (i = 0; i < bmap->n_ineq; ++i) {
715                 if (!isl_int_is_zero(bmap->ineq[i][pos])) {
716                         isl_basic_map_drop_inequality(bmap, i);
717                         --i;
718                         continue;
719                 }
720                 constraint_drop_vars(bmap->ineq[i]+pos, 1, bmap->extra-div-1);
721         }
722
723         for (i = 0; i < bmap->n_div; ++i)
724                 constraint_drop_vars(bmap->div[i]+1+pos, 1, bmap->extra-div-1);
725
726         if (div != bmap->n_div - 1) {
727                 int j;
728                 isl_int *t = bmap->div[div];
729
730                 for (j = div; j < bmap->n_div - 1; ++j)
731                         bmap->div[j] = bmap->div[j+1];
732
733                 bmap->div[bmap->n_div - 1] = t;
734         }
735         isl_basic_map_free_div(bmap, 1);
736
737         return bmap;
738 error:
739         isl_basic_map_free(bmap);
740         return NULL;
741 }
742
743 static void eliminate_div(struct isl_basic_map *bmap, isl_int *eq, unsigned div)
744 {
745         int i;
746         unsigned pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
747         unsigned len;
748         len = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
749
750         for (i = 0; i < bmap->n_eq; ++i)
751                 if (bmap->eq[i] != eq)
752                         isl_seq_elim(bmap->eq[i], eq, pos, len, NULL);
753
754         for (i = 0; i < bmap->n_ineq; ++i)
755                 isl_seq_elim(bmap->ineq[i], eq, pos, len, NULL);
756
757         /* We need to be careful about circular definitions,
758          * so for now we just remove the definitions of other divs that
759          * depend on this div and (possibly) recompute them later.
760          */
761         for (i = 0; i < bmap->n_div; ++i)
762                 if (!isl_int_is_zero(bmap->div[i][0]) &&
763                     !isl_int_is_zero(bmap->div[i][1 + pos]))
764                         isl_seq_clr(bmap->div[i], 1 + len);
765
766         isl_basic_map_drop_div(bmap, div);
767 }
768
769 struct isl_basic_map *isl_basic_map_set_to_empty(struct isl_basic_map *bmap)
770 {
771         int i = 0;
772         unsigned total;
773         if (!bmap)
774                 goto error;
775         total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra;
776         isl_basic_map_free_div(bmap, bmap->n_div);
777         isl_basic_map_free_inequality(bmap, bmap->n_ineq);
778         if (bmap->n_eq > 0)
779                 isl_basic_map_free_equality(bmap, bmap->n_eq-1);
780         else {
781                 isl_basic_map_alloc_equality(bmap);
782                 if (i < 0)
783                         goto error;
784         }
785         isl_int_set_si(bmap->eq[i][0], 1);
786         isl_seq_clr(bmap->eq[i]+1, total);
787         F_SET(bmap, ISL_BASIC_MAP_EMPTY);
788         return isl_basic_map_finalize(bmap);
789 error:
790         isl_basic_map_free(bmap);
791         return NULL;
792 }
793
794 struct isl_basic_set *isl_basic_set_set_to_empty(struct isl_basic_set *bset)
795 {
796         return (struct isl_basic_set *)
797                 isl_basic_map_set_to_empty((struct isl_basic_map *)bset);
798 }
799
800 static void swap_equality(struct isl_basic_map *bmap, int a, int b)
801 {
802         isl_int *t = bmap->eq[a];
803         bmap->eq[a] = bmap->eq[b];
804         bmap->eq[b] = t;
805 }
806
807 static void swap_inequality(struct isl_basic_map *bmap, int a, int b)
808 {
809         if (a != b) {
810                 isl_int *t = bmap->ineq[a];
811                 bmap->ineq[a] = bmap->ineq[b];
812                 bmap->ineq[b] = t;
813         }
814 }
815
816 static void set_swap_inequality(struct isl_basic_set *bset, int a, int b)
817 {
818         swap_inequality((struct isl_basic_map *)bset, a, b);
819 }
820
821 static void swap_div(struct isl_basic_map *bmap, int a, int b)
822 {
823         int i;
824         unsigned off = bmap->nparam + bmap->n_in + bmap->n_out;
825         isl_int *t = bmap->div[a];
826         bmap->div[a] = bmap->div[b];
827         bmap->div[b] = t;
828
829         for (i = 0; i < bmap->n_eq; ++i)
830                 isl_int_swap(bmap->eq[i][1+off+a], bmap->eq[i][1+off+b]);
831
832         for (i = 0; i < bmap->n_ineq; ++i)
833                 isl_int_swap(bmap->ineq[i][1+off+a], bmap->ineq[i][1+off+b]);
834
835         for (i = 0; i < bmap->n_div; ++i)
836                 isl_int_swap(bmap->div[i][1+1+off+a], bmap->div[i][1+1+off+b]);
837 }
838
839 struct isl_basic_map *isl_basic_map_gauss(
840         struct isl_basic_map *bmap, int *progress)
841 {
842         int k;
843         int done;
844         int last_var;
845         unsigned total_var;
846         unsigned total;
847
848         if (!bmap)
849                 return NULL;
850
851         total_var = bmap->nparam + bmap->n_in + bmap->n_out;
852         total = total_var + bmap->n_div;
853
854         last_var = total - 1;
855         for (done = 0; done < bmap->n_eq; ++done) {
856                 for (; last_var >= 0; --last_var) {
857                         for (k = done; k < bmap->n_eq; ++k)
858                                 if (!isl_int_is_zero(bmap->eq[k][1+last_var]))
859                                         break;
860                         if (k < bmap->n_eq)
861                                 break;
862                 }
863                 if (last_var < 0)
864                         break;
865                 if (k != done)
866                         swap_equality(bmap, k, done);
867                 if (isl_int_is_neg(bmap->eq[done][1+last_var]))
868                         isl_seq_neg(bmap->eq[done], bmap->eq[done], 1+total);
869
870                 for (k = 0; k < bmap->n_eq; ++k) {
871                         if (k == done)
872                                 continue;
873                         if (isl_int_is_zero(bmap->eq[k][1+last_var]))
874                                 continue;
875                         if (progress)
876                                 *progress = 1;
877                         isl_seq_elim(bmap->eq[k], bmap->eq[done],
878                                         1+last_var, 1+total, NULL);
879                 }
880
881                 for (k = 0; k < bmap->n_ineq; ++k) {
882                         if (isl_int_is_zero(bmap->ineq[k][1+last_var]))
883                                 continue;
884                         if (progress)
885                                 *progress = 1;
886                         isl_seq_elim(bmap->ineq[k], bmap->eq[done],
887                                         1+last_var, 1+total, NULL);
888                 }
889
890                 for (k = 0; k < bmap->n_div; ++k) {
891                         if (isl_int_is_zero(bmap->div[k][0]))
892                                 continue;
893                         if (isl_int_is_zero(bmap->div[k][1+1+last_var]))
894                                 continue;
895                         if (progress)
896                                 *progress = 1;
897                         isl_seq_elim(bmap->div[k]+1, bmap->eq[done],
898                                         1+last_var, 1+total, &bmap->div[k][0]);
899                 }
900
901                 if (last_var >= total_var &&
902                     isl_int_is_zero(bmap->div[last_var - total_var][0])) {
903                         unsigned div = last_var - total_var;
904                         isl_seq_neg(bmap->div[div]+1, bmap->eq[done], 1+total);
905                         isl_int_set_si(bmap->div[div][1+1+last_var], 0);
906                         isl_int_set(bmap->div[div][0],
907                                     bmap->eq[done][1+last_var]);
908                 }
909         }
910         if (done == bmap->n_eq)
911                 return bmap;
912         for (k = done; k < bmap->n_eq; ++k) {
913                 if (isl_int_is_zero(bmap->eq[k][0]))
914                         continue;
915                 return isl_basic_map_set_to_empty(bmap);
916         }
917         isl_basic_map_free_equality(bmap, bmap->n_eq-done);
918         return bmap;
919 }
920
921 struct isl_basic_set *isl_basic_set_gauss(
922         struct isl_basic_set *bset, int *progress)
923 {
924         return (struct isl_basic_set*)isl_basic_map_gauss(
925                         (struct isl_basic_map *)bset, progress);
926 }
927
928 static unsigned int round_up(unsigned int v)
929 {
930         int old_v = v;
931
932         while (v) {
933                 old_v = v;
934                 v ^= v & -v;
935         }
936         return old_v << 1;
937 }
938
939 static int hash_index(int *index, unsigned int size, int bits,
940                         struct isl_basic_map *bmap, int k)
941 {
942         int h;
943         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
944         u_int32_t hash = isl_seq_hash(bmap->ineq[k]+1, total, bits);
945         for (h = hash; index[h]; h = (h+1) % size)
946                 if (k != index[h]-1 &&
947                     isl_seq_eq(bmap->ineq[k]+1,
948                                bmap->ineq[index[h]-1]+1, total))
949                         break;
950         return h;
951 }
952
953 static struct isl_basic_map *remove_duplicate_constraints(
954         struct isl_basic_map *bmap, int *progress)
955 {
956         unsigned int size;
957         int *index;
958         int k, l, h;
959         int bits;
960         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
961         isl_int sum;
962
963         if (bmap->n_ineq <= 1)
964                 return bmap;
965
966         size = round_up(4 * (bmap->n_ineq+1) / 3 - 1);
967         bits = ffs(size) - 1;
968         index = isl_alloc_array(ctx, int, size);
969         memset(index, 0, size * sizeof(int));
970         if (!index)
971                 return bmap;
972
973         index[isl_seq_hash(bmap->ineq[0]+1, total, bits)] = 1;
974         for (k = 1; k < bmap->n_ineq; ++k) {
975                 h = hash_index(index, size, bits, bmap, k);
976                 if (!index[h]) {
977                         index[h] = k+1;
978                         continue;
979                 }
980                 if (progress)
981                         *progress = 1;
982                 l = index[h] - 1;
983                 if (isl_int_lt(bmap->ineq[k][0], bmap->ineq[l][0]))
984                         swap_inequality(bmap, k, l);
985                 isl_basic_map_drop_inequality(bmap, k);
986                 --k;
987         }
988         isl_int_init(sum);
989         for (k = 0; k < bmap->n_ineq-1; ++k) {
990                 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
991                 h = hash_index(index, size, bits, bmap, k);
992                 isl_seq_neg(bmap->ineq[k]+1, bmap->ineq[k]+1, total);
993                 if (!index[h])
994                         continue;
995                 l = index[h] - 1;
996                 isl_int_add(sum, bmap->ineq[k][0], bmap->ineq[l][0]);
997                 if (isl_int_is_pos(sum))
998                         continue;
999                 if (isl_int_is_zero(sum)) {
1000                         /* We need to break out of the loop after these
1001                          * changes since the contents of the hash
1002                          * will no longer be valid.
1003                          * Plus, we probably we want to regauss first.
1004                          */
1005                         isl_basic_map_drop_inequality(bmap, l);
1006                         isl_basic_map_inequality_to_equality(bmap, k);
1007                 } else
1008                         bmap = isl_basic_map_set_to_empty(bmap);
1009                 break;
1010         }
1011         isl_int_clear(sum);
1012
1013         free(index);
1014         return bmap;
1015 }
1016
1017 static struct isl_basic_map *remove_duplicate_divs(
1018         struct isl_basic_map *bmap, int *progress)
1019 {
1020         unsigned int size;
1021         int *index;
1022         int k, l, h;
1023         int bits;
1024         struct isl_blk eq;
1025         unsigned total_var = bmap->nparam + bmap->n_in + bmap->n_out;
1026         unsigned total = total_var + bmap->n_div;
1027         struct isl_ctx *ctx;
1028
1029         if (bmap->n_div <= 1)
1030                 return bmap;
1031
1032         ctx = bmap->ctx;
1033         for (k = bmap->n_div - 1; k >= 0; --k)
1034                 if (!isl_int_is_zero(bmap->div[k][0]))
1035                         break;
1036         if (k <= 0)
1037                 return bmap;
1038
1039         size = round_up(4 * bmap->n_div / 3 - 1);
1040         bits = ffs(size) - 1;
1041         index = isl_alloc_array(ctx, int, size);
1042         memset(index, 0, size * sizeof(int));
1043         if (!index)
1044                 return bmap;
1045         eq = isl_blk_alloc(ctx, 1+total);
1046         if (isl_blk_is_error(eq))
1047                 goto out;
1048
1049         isl_seq_clr(eq.data, 1+total);
1050         index[isl_seq_hash(bmap->div[k], 2+total, bits)] = k + 1;
1051         for (--k; k >= 0; --k) {
1052                 u_int32_t hash;
1053
1054                 if (isl_int_is_zero(bmap->div[k][0]))
1055                         continue;
1056
1057                 hash = isl_seq_hash(bmap->div[k], 2+total, bits);
1058                 for (h = hash; index[h]; h = (h+1) % size)
1059                         if (isl_seq_eq(bmap->div[k],
1060                                        bmap->div[index[h]-1], 2+total))
1061                                 break;
1062                 if (index[h]) {
1063                         *progress = 1;
1064                         l = index[h] - 1;
1065                         isl_int_set_si(eq.data[1+total_var+k], -1);
1066                         isl_int_set_si(eq.data[1+total_var+l], 1);
1067                         eliminate_div(bmap, eq.data, l);
1068                         isl_int_set_si(eq.data[1+total_var+k], 0);
1069                         isl_int_set_si(eq.data[1+total_var+l], 0);
1070                 }
1071                 index[h] = k+1;
1072         }
1073
1074         isl_blk_free(ctx, eq);
1075 out:
1076         free(index);
1077         return bmap;
1078 }
1079
1080 /* Elimininate divs based on equalities
1081  */
1082 static struct isl_basic_map *eliminate_divs_eq(
1083                 struct isl_basic_map *bmap, int *progress)
1084 {
1085         int d;
1086         int i;
1087         int modified = 0;
1088         unsigned off;
1089
1090         if (!bmap)
1091                 return NULL;
1092
1093         off = 1 + bmap->nparam + bmap->n_in + bmap->n_out;
1094
1095         for (d = bmap->n_div - 1; d >= 0 ; --d) {
1096                 for (i = 0; i < bmap->n_eq; ++i) {
1097                         if (!isl_int_is_one(bmap->eq[i][off + d]) &&
1098                             !isl_int_is_negone(bmap->eq[i][off + d]))
1099                                 continue;
1100                         modified = 1;
1101                         *progress = 1;
1102                         eliminate_div(bmap, bmap->eq[i], d);
1103                         isl_basic_map_drop_equality(bmap, i);
1104                         break;
1105                 }
1106         }
1107         if (modified)
1108                 return eliminate_divs_eq(bmap, progress);
1109         return bmap;
1110 }
1111
1112 static struct isl_basic_map *normalize_constraints(struct isl_basic_map *bmap)
1113 {
1114         int i;
1115         isl_int gcd;
1116         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1117
1118         isl_int_init(gcd);
1119         for (i = bmap->n_eq - 1; i >= 0; --i) {
1120                 isl_seq_gcd(bmap->eq[i]+1, total, &gcd);
1121                 if (isl_int_is_zero(gcd)) {
1122                         if (!isl_int_is_zero(bmap->eq[i][0])) {
1123                                 bmap = isl_basic_map_set_to_empty(bmap);
1124                                 break;
1125                         }
1126                         isl_basic_map_drop_equality(bmap, i);
1127                         continue;
1128                 }
1129                 if (F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
1130                         isl_int_gcd(gcd, gcd, bmap->eq[i][0]);
1131                 if (isl_int_is_one(gcd))
1132                         continue;
1133                 if (!isl_int_is_divisible_by(bmap->eq[i][0], gcd)) {
1134                         bmap = isl_basic_map_set_to_empty(bmap);
1135                         break;
1136                 }
1137                 isl_seq_scale_down(bmap->eq[i], bmap->eq[i], gcd, 1+total);
1138         }
1139
1140         for (i = bmap->n_ineq - 1; i >= 0; --i) {
1141                 isl_seq_gcd(bmap->ineq[i]+1, total, &gcd);
1142                 if (isl_int_is_zero(gcd)) {
1143                         if (isl_int_is_neg(bmap->ineq[i][0])) {
1144                                 bmap = isl_basic_map_set_to_empty(bmap);
1145                                 break;
1146                         }
1147                         isl_basic_map_drop_inequality(bmap, i);
1148                         continue;
1149                 }
1150                 if (F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL))
1151                         isl_int_gcd(gcd, gcd, bmap->ineq[i][0]);
1152                 if (isl_int_is_one(gcd))
1153                         continue;
1154                 isl_int_fdiv_q(bmap->ineq[i][0], bmap->ineq[i][0], gcd);
1155                 isl_seq_scale_down(bmap->ineq[i]+1, bmap->ineq[i]+1, gcd, total);
1156         }
1157         isl_int_clear(gcd);
1158
1159         return bmap;
1160 }
1161
1162
1163 /* Eliminate the specified variables from the constraints using
1164  * Fourier-Motzkin.  The variables themselves are not removed.
1165  */
1166 struct isl_basic_map *isl_basic_map_eliminate_vars(
1167         struct isl_basic_map *bmap, int pos, unsigned n)
1168 {
1169         int d;
1170         int i, j, k;
1171         unsigned total;
1172         unsigned extra;
1173
1174         if (!bmap)
1175                 return NULL;
1176         total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1177         extra = bmap->extra - bmap->n_div;
1178
1179         bmap = isl_basic_map_cow(bmap);
1180         bmap = isl_basic_map_gauss(bmap, NULL);
1181         for (d = pos + n - 1; d >= pos; --d) {
1182                 int n_lower, n_upper;
1183                 if (!bmap)
1184                         return NULL;
1185                 for (i = 0; i < bmap->n_eq; ++i) {
1186                         if (isl_int_is_zero(bmap->eq[i][1+d]))
1187                                 continue;
1188                         isl_basic_map_drop_equality(bmap, i);
1189                         break;
1190                 }
1191                 if (i < bmap->n_eq)
1192                         continue;
1193                 n_lower = 0;
1194                 n_upper = 0;
1195                 for (i = 0; i < bmap->n_ineq; ++i) {
1196                         if (isl_int_is_pos(bmap->ineq[i][1+d]))
1197                                 n_lower++;
1198                         else if (isl_int_is_neg(bmap->ineq[i][1+d]))
1199                                 n_upper++;
1200                 }
1201                 bmap = isl_basic_map_extend(bmap,
1202                                 bmap->nparam, bmap->n_in, bmap->n_out, 0,
1203                                 0, n_lower * n_upper);
1204                 for (i = bmap->n_ineq - 1; i >= 0; --i) {
1205                         int last;
1206                         if (isl_int_is_zero(bmap->ineq[i][1+d]))
1207                                 continue;
1208                         last = -1;
1209                         for (j = 0; j < i; ++j) {
1210                                 if (isl_int_is_zero(bmap->ineq[j][1+d]))
1211                                         continue;
1212                                 last = j;
1213                                 if (isl_int_sgn(bmap->ineq[i][1+d]) ==
1214                                     isl_int_sgn(bmap->ineq[j][1+d]))
1215                                         continue;
1216                                 k = isl_basic_map_alloc_inequality(bmap);
1217                                 if (k < 0)
1218                                         goto error;
1219                                 isl_seq_cpy(bmap->ineq[k], bmap->ineq[i],
1220                                                 1+total);
1221                                 isl_seq_clr(bmap->ineq[k]+1+total, extra);
1222                                 isl_seq_elim(bmap->ineq[k], bmap->ineq[j],
1223                                                 1+d, 1+total, NULL);
1224                         }
1225                         isl_basic_map_drop_inequality(bmap, i);
1226                         i = last + 1;
1227                 }
1228                 if (n_lower > 0 && n_upper > 0) {
1229                         bmap = normalize_constraints(bmap);
1230                         bmap = remove_duplicate_constraints(bmap, NULL);
1231                         bmap = isl_basic_map_gauss(bmap, NULL);
1232                         bmap = isl_basic_map_convex_hull(bmap);
1233                         if (!bmap)
1234                                 goto error;
1235                         if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1236                                 break;
1237                 }
1238         }
1239         return bmap;
1240 error:
1241         isl_basic_map_free(bmap);
1242         return NULL;
1243 }
1244
1245 struct isl_basic_set *isl_basic_set_eliminate_vars(
1246         struct isl_basic_set *bset, unsigned pos, unsigned n)
1247 {
1248         return (struct isl_basic_set *)isl_basic_map_eliminate_vars(
1249                         (struct isl_basic_map *)bset, pos, n);
1250 }
1251
1252 /* Project out n dimensions starting at first using Fourier-Motzkin */
1253 struct isl_set *isl_set_remove_dims(struct isl_set *set,
1254         unsigned first, unsigned n)
1255 {
1256         int i;
1257
1258         if (n == 0)
1259                 return set;
1260
1261         set = isl_set_cow(set);
1262         if (!set)
1263                 return NULL;
1264         isl_assert(set->ctx, first+n <= set->dim, goto error);
1265         
1266         for (i = 0; i < set->n; ++i) {
1267                 set->p[i] = isl_basic_set_eliminate_vars(set->p[i],
1268                                                         set->nparam + first, n);
1269                 if (!set->p[i])
1270                         goto error;
1271         }
1272         set = isl_set_drop_dims(set, first, n);
1273         return set;
1274 error:
1275         isl_set_free(set);
1276         return NULL;
1277 }
1278
1279 /* Project out n dimensions starting at first using Fourier-Motzkin */
1280 struct isl_basic_set *isl_basic_set_remove_dims(struct isl_basic_set *bset,
1281         unsigned first, unsigned n)
1282 {
1283         bset = isl_basic_set_eliminate_vars(bset, bset->nparam + first, n);
1284         bset = isl_basic_set_drop_dims(bset, first, n);
1285         return bset;
1286 }
1287
1288 /* Elimininate divs based on inequalities
1289  */
1290 static struct isl_basic_map *eliminate_divs_ineq(
1291                 struct isl_basic_map *bmap, int *progress)
1292 {
1293         int d;
1294         int i;
1295         unsigned off;
1296         struct isl_ctx *ctx;
1297
1298         if (!bmap)
1299                 return NULL;
1300
1301         ctx = bmap->ctx;
1302         off = 1 + bmap->nparam + bmap->n_in + bmap->n_out;
1303
1304         for (d = bmap->n_div - 1; d >= 0 ; --d) {
1305                 for (i = 0; i < bmap->n_eq; ++i)
1306                         if (!isl_int_is_zero(bmap->eq[i][off + d]))
1307                                 break;
1308                 if (i < bmap->n_eq)
1309                         continue;
1310                 for (i = 0; i < bmap->n_ineq; ++i)
1311                         if (isl_int_abs_gt(bmap->ineq[i][off + d], ctx->one))
1312                                 break;
1313                 if (i < bmap->n_ineq)
1314                         continue;
1315                 *progress = 1;
1316                 bmap = isl_basic_map_eliminate_vars(bmap, (off-1)+d, 1);
1317                 if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
1318                         break;
1319                 bmap = isl_basic_map_drop_div(bmap, d);
1320                 if (!bmap)
1321                         break;
1322         }
1323         return bmap;
1324 }
1325
1326 struct isl_basic_map *isl_basic_map_simplify(struct isl_basic_map *bmap)
1327 {
1328         int progress = 1;
1329         if (!bmap)
1330                 return NULL;
1331         while (progress) {
1332                 progress = 0;
1333                 bmap = normalize_constraints(bmap);
1334                 bmap = eliminate_divs_eq(bmap, &progress);
1335                 bmap = eliminate_divs_ineq(bmap, &progress);
1336                 bmap = isl_basic_map_gauss(bmap, &progress);
1337                 bmap = remove_duplicate_divs(bmap, &progress);
1338                 bmap = remove_duplicate_constraints(bmap, &progress);
1339         }
1340         return bmap;
1341 }
1342
1343 struct isl_basic_set *isl_basic_set_simplify(struct isl_basic_set *bset)
1344 {
1345         return (struct isl_basic_set *)
1346                 isl_basic_map_simplify((struct isl_basic_map *)bset);
1347 }
1348
1349 static void dump_term(struct isl_basic_map *bmap,
1350                         isl_int c, int pos, FILE *out)
1351 {
1352         unsigned in = bmap->n_in;
1353         unsigned dim = bmap->n_in + bmap->n_out;
1354         if (!pos)
1355                 isl_int_print(out, c);
1356         else {
1357                 if (!isl_int_is_one(c))
1358                         isl_int_print(out, c);
1359                 if (pos < 1 + bmap->nparam)
1360                         fprintf(out, "p%d", pos - 1);
1361                 else if (pos < 1 + bmap->nparam + in)
1362                         fprintf(out, "i%d", pos - 1 - bmap->nparam);
1363                 else if (pos < 1 + bmap->nparam + dim)
1364                         fprintf(out, "o%d", pos - 1 - bmap->nparam - in);
1365                 else
1366                         fprintf(out, "e%d", pos - 1 - bmap->nparam - dim);
1367         }
1368 }
1369
1370 static void dump_constraint_sign(struct isl_basic_map *bmap, isl_int *c,
1371                                 int sign, FILE *out)
1372 {
1373         int i;
1374         int first;
1375         unsigned len = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1376         isl_int v;
1377
1378         isl_int_init(v);
1379         for (i = 0, first = 1; i < len; ++i) {
1380                 if (isl_int_sgn(c[i]) * sign <= 0)
1381                         continue;
1382                 if (!first)
1383                         fprintf(out, " + ");
1384                 first = 0;
1385                 isl_int_abs(v, c[i]);
1386                 dump_term(bmap, v, i, out);
1387         }
1388         isl_int_clear(v);
1389         if (first)
1390                 fprintf(out, "0");
1391 }
1392
1393 static void dump_constraint(struct isl_basic_map *bmap, isl_int *c,
1394                                 const char *op, FILE *out, int indent)
1395 {
1396         int i;
1397
1398         fprintf(out, "%*s", indent, "");
1399
1400         dump_constraint_sign(bmap, c, 1, out);
1401         fprintf(out, " %s ", op);
1402         dump_constraint_sign(bmap, c, -1, out);
1403
1404         fprintf(out, "\n");
1405
1406         for (i = bmap->n_div; i < bmap->extra; ++i) {
1407                 if (isl_int_is_zero(c[1+bmap->nparam+bmap->n_in+bmap->n_out+i]))
1408                         continue;
1409                 fprintf(out, "%*s", indent, "");
1410                 fprintf(out, "ERROR: unused div coefficient not zero\n");
1411         }
1412 }
1413
1414 static void dump_constraints(struct isl_basic_map *bmap,
1415                                 isl_int **c, unsigned n,
1416                                 const char *op, FILE *out, int indent)
1417 {
1418         int i;
1419
1420         for (i = 0; i < n; ++i)
1421                 dump_constraint(bmap, c[i], op, out, indent);
1422 }
1423
1424 static void dump_affine(struct isl_basic_map *bmap, isl_int *exp, FILE *out)
1425 {
1426         int j;
1427         int first = 1;
1428         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1429
1430         for (j = 0; j < 1 + total; ++j) {
1431                 if (isl_int_is_zero(exp[j]))
1432                         continue;
1433                 if (!first && isl_int_is_pos(exp[j]))
1434                         fprintf(out, "+");
1435                 dump_term(bmap, exp[j], j, out);
1436                 first = 0;
1437         }
1438 }
1439
1440 static void dump(struct isl_basic_map *bmap, FILE *out, int indent)
1441 {
1442         int i;
1443
1444         dump_constraints(bmap, bmap->eq, bmap->n_eq, "=", out, indent);
1445         dump_constraints(bmap, bmap->ineq, bmap->n_ineq, ">=", out, indent);
1446
1447         for (i = 0; i < bmap->n_div; ++i) {
1448                 fprintf(out, "%*s", indent, "");
1449                 fprintf(out, "e%d = [(", i);
1450                 dump_affine(bmap, bmap->div[i]+1, out);
1451                 fprintf(out, ")/");
1452                 isl_int_print(out, bmap->div[i][0]);
1453                 fprintf(out, "]\n");
1454         }
1455 }
1456
1457 void isl_basic_set_dump(struct isl_basic_set *bset, FILE *out, int indent)
1458 {
1459         if (!bset) {
1460                 fprintf(out, "null basic set\n");
1461                 return;
1462         }
1463
1464         fprintf(out, "%*s", indent, "");
1465         fprintf(out, "ref: %d, nparam: %d, dim: %d, extra: %d, flags: %x\n",
1466                         bset->ref, bset->nparam, bset->dim, bset->extra,
1467                         bset->flags);
1468         dump((struct isl_basic_map *)bset, out, indent);
1469 }
1470
1471 void isl_basic_map_dump(struct isl_basic_map *bmap, FILE *out, int indent)
1472 {
1473         if (!bmap) {
1474                 fprintf(out, "null basic map\n");
1475                 return;
1476         }
1477
1478         fprintf(out, "%*s", indent, "");
1479         fprintf(out, "ref: %d, nparam: %d, in: %d, out: %d, extra: %d, flags: %x\n",
1480                 bmap->ref,
1481                 bmap->nparam, bmap->n_in, bmap->n_out, bmap->extra, bmap->flags);
1482         dump(bmap, out, indent);
1483 }
1484
1485 int isl_inequality_negate(struct isl_basic_map *bmap, unsigned pos)
1486 {
1487         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
1488         if (!bmap)
1489                 return -1;
1490         isl_assert(bmap->ctx, pos < bmap->n_ineq, return -1);
1491         isl_seq_neg(bmap->ineq[pos], bmap->ineq[pos], 1 + total);
1492         isl_int_sub_ui(bmap->ineq[pos][0], bmap->ineq[pos][0], 1);
1493         return 0;
1494 }
1495
1496 struct isl_set *isl_set_alloc(struct isl_ctx *ctx,
1497                 unsigned nparam, unsigned dim, int n, unsigned flags)
1498 {
1499         struct isl_set *set;
1500
1501         isl_assert(ctx, n >= 0, return NULL);
1502         set = isl_alloc(ctx, struct isl_set,
1503                         sizeof(struct isl_set) +
1504                         n * sizeof(struct isl_basic_set *));
1505         if (!set)
1506                 return NULL;
1507
1508         set->ctx = ctx;
1509         isl_ctx_ref(ctx);
1510         set->ref = 1;
1511         set->size = n;
1512         set->n = 0;
1513         set->nparam = nparam;
1514         set->zero = 0;
1515         set->dim = dim;
1516         set->flags = flags;
1517         return set;
1518 }
1519
1520 struct isl_set *isl_set_dup(struct isl_set *set)
1521 {
1522         int i;
1523         struct isl_set *dup;
1524
1525         if (!set)
1526                 return NULL;
1527
1528         dup = isl_set_alloc(set->ctx, set->nparam, set->dim, set->n, set->flags);
1529         if (!dup)
1530                 return NULL;
1531         for (i = 0; i < set->n; ++i)
1532                 dup = isl_set_add(dup,
1533                                   isl_basic_set_dup(set->p[i]));
1534         return dup;
1535 }
1536
1537 struct isl_set *isl_set_from_basic_set(struct isl_basic_set *bset)
1538 {
1539         struct isl_set *set;
1540
1541         if (!bset)
1542                 return NULL;
1543
1544         set = isl_set_alloc(bset->ctx, bset->nparam, bset->dim, 1, ISL_MAP_DISJOINT);
1545         if (!set) {
1546                 isl_basic_set_free(bset);
1547                 return NULL;
1548         }
1549         return isl_set_add(set, bset);
1550 }
1551
1552 struct isl_map *isl_map_from_basic_map(struct isl_basic_map *bmap)
1553 {
1554         struct isl_map *map;
1555
1556         if (!bmap)
1557                 return NULL;
1558
1559         map = isl_map_alloc(bmap->ctx, bmap->nparam, bmap->n_in, bmap->n_out, 1,
1560                                 ISL_MAP_DISJOINT);
1561         if (!map) {
1562                 isl_basic_map_free(bmap);
1563                 return NULL;
1564         }
1565         return isl_map_add(map, bmap);
1566 }
1567
1568 struct isl_set *isl_set_add(struct isl_set *set, struct isl_basic_set *bset)
1569 {
1570         if (!bset || !set)
1571                 goto error;
1572         isl_assert(set->ctx, set->nparam == bset->nparam, goto error);
1573         isl_assert(set->ctx, set->dim == bset->dim, goto error);
1574         isl_assert(set->ctx, set->n < set->size, goto error);
1575         set->p[set->n] = bset;
1576         set->n++;
1577         return set;
1578 error:
1579         if (set)
1580                 isl_set_free(set);
1581         if (bset)
1582                 isl_basic_set_free(bset);
1583         return NULL;
1584 }
1585
1586 void isl_set_free(struct isl_set *set)
1587 {
1588         int i;
1589
1590         if (!set)
1591                 return;
1592
1593         if (--set->ref > 0)
1594                 return;
1595
1596         isl_ctx_deref(set->ctx);
1597         for (i = 0; i < set->n; ++i)
1598                 isl_basic_set_free(set->p[i]);
1599         free(set);
1600 }
1601
1602 void isl_set_dump(struct isl_set *set, FILE *out, int indent)
1603 {
1604         int i;
1605
1606         if (!set) {
1607                 fprintf(out, "null set\n");
1608                 return;
1609         }
1610
1611         fprintf(out, "%*s", indent, "");
1612         fprintf(out, "ref: %d, n: %d, nparam: %d, dim: %d, flags: %x\n",
1613                         set->ref, set->n, set->nparam, set->dim, set->flags);
1614         for (i = 0; i < set->n; ++i) {
1615                 fprintf(out, "%*s", indent, "");
1616                 fprintf(out, "basic set %d:\n", i);
1617                 isl_basic_set_dump(set->p[i], out, indent+4);
1618         }
1619 }
1620
1621 void isl_map_dump(struct isl_map *map, FILE *out, int indent)
1622 {
1623         int i;
1624
1625         if (!map) {
1626                 fprintf(out, "null map\n");
1627                 return;
1628         }
1629
1630         fprintf(out, "%*s", indent, "");
1631         fprintf(out, "ref: %d, n: %d, nparam: %d, in: %d, out: %d, flags: %x\n",
1632                         map->ref, map->n, map->nparam, map->n_in, map->n_out,
1633                         map->flags);
1634         for (i = 0; i < map->n; ++i) {
1635                 fprintf(out, "%*s", indent, "");
1636                 fprintf(out, "basic map %d:\n", i);
1637                 isl_basic_map_dump(map->p[i], out, indent+4);
1638         }
1639 }
1640
1641 struct isl_basic_map *isl_basic_map_intersect_domain(
1642                 struct isl_basic_map *bmap, struct isl_basic_set *bset)
1643 {
1644         struct isl_basic_map *bmap_domain;
1645
1646         if (!bmap || !bset)
1647                 goto error;
1648
1649         isl_assert(bset->ctx, bset->dim == bmap->n_in, goto error);
1650         isl_assert(bset->ctx, bset->nparam == bmap->nparam, goto error);
1651
1652         bmap = isl_basic_map_extend(bmap,
1653                         bset->nparam, bmap->n_in, bmap->n_out,
1654                         bset->n_div, bset->n_eq, bset->n_ineq);
1655         if (!bmap)
1656                 goto error;
1657         bmap_domain = isl_basic_map_from_basic_set(bset,
1658                                                 bset->dim, 0);
1659         bmap = add_constraints(bmap, bmap_domain, 0, 0);
1660
1661         bmap = isl_basic_map_simplify(bmap);
1662         return isl_basic_map_finalize(bmap);
1663 error:
1664         isl_basic_map_free(bmap);
1665         isl_basic_set_free(bset);
1666         return NULL;
1667 }
1668
1669 struct isl_basic_map *isl_basic_map_intersect_range(
1670                 struct isl_basic_map *bmap, struct isl_basic_set *bset)
1671 {
1672         struct isl_basic_map *bmap_range;
1673
1674         if (!bmap || !bset)
1675                 goto error;
1676
1677         isl_assert(bset->ctx, bset->dim == bmap->n_out, goto error);
1678         isl_assert(bset->ctx, bset->nparam == bmap->nparam, goto error);
1679
1680         bmap = isl_basic_map_extend(bmap,
1681                         bset->nparam, bmap->n_in, bmap->n_out,
1682                         bset->n_div, bset->n_eq, bset->n_ineq);
1683         if (!bmap)
1684                 goto error;
1685         bmap_range = isl_basic_map_from_basic_set(bset,
1686                                                 0, bset->dim);
1687         bmap = add_constraints(bmap, bmap_range, 0, 0);
1688
1689         bmap = isl_basic_map_simplify(bmap);
1690         return isl_basic_map_finalize(bmap);
1691 error:
1692         isl_basic_map_free(bmap);
1693         isl_basic_set_free(bset);
1694         return NULL;
1695 }
1696
1697 struct isl_basic_map *isl_basic_map_intersect(
1698                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
1699 {
1700         if (!bmap1 || !bmap2)
1701                 goto error;
1702
1703         isl_assert(bmap1->ctx, bmap1->nparam == bmap2->nparam, goto error);
1704         isl_assert(bmap1->ctx, bmap1->n_in == bmap2->n_in, goto error);
1705         isl_assert(bmap1->ctx, bmap1->n_out == bmap2->n_out, goto error);
1706
1707         bmap1 = isl_basic_map_extend(bmap1,
1708                         bmap1->nparam, bmap1->n_in, bmap1->n_out,
1709                         bmap2->n_div, bmap2->n_eq, bmap2->n_ineq);
1710         if (!bmap1)
1711                 goto error;
1712         bmap1 = add_constraints(bmap1, bmap2, 0, 0);
1713
1714         bmap1 = isl_basic_map_simplify(bmap1);
1715         return isl_basic_map_finalize(bmap1);
1716 error:
1717         isl_basic_map_free(bmap1);
1718         isl_basic_map_free(bmap2);
1719         return NULL;
1720 }
1721
1722 struct isl_basic_set *isl_basic_set_intersect(
1723                 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
1724 {
1725         return (struct isl_basic_set *)
1726                 isl_basic_map_intersect(
1727                         (struct isl_basic_map *)bset1,
1728                         (struct isl_basic_map *)bset2);
1729 }
1730
1731 struct isl_map *isl_map_intersect(struct isl_map *map1, struct isl_map *map2)
1732 {
1733         unsigned flags = 0;
1734         struct isl_map *result;
1735         int i, j;
1736
1737         if (!map1 || !map2)
1738                 goto error;
1739
1740         if (F_ISSET(map1, ISL_MAP_DISJOINT) &&
1741             F_ISSET(map2, ISL_MAP_DISJOINT))
1742                 FL_SET(flags, ISL_MAP_DISJOINT);
1743
1744         result = isl_map_alloc(map1->ctx, map1->nparam, map1->n_in, map1->n_out,
1745                                 map1->n * map2->n, flags);
1746         if (!result)
1747                 goto error;
1748         for (i = 0; i < map1->n; ++i)
1749                 for (j = 0; j < map2->n; ++j) {
1750                         struct isl_basic_map *part;
1751                         part = isl_basic_map_intersect(
1752                                     isl_basic_map_copy(map1->p[i]),
1753                                     isl_basic_map_copy(map2->p[j]));
1754                         if (isl_basic_map_is_empty(part))
1755                                 isl_basic_map_free(part);
1756                         else
1757                                 result = isl_map_add(result, part);
1758                         if (!result)
1759                                 goto error;
1760                 }
1761         isl_map_free(map1);
1762         isl_map_free(map2);
1763         return result;
1764 error:
1765         isl_map_free(map1);
1766         isl_map_free(map2);
1767         return NULL;
1768 }
1769
1770 struct isl_set *isl_set_intersect(struct isl_set *set1, struct isl_set *set2)
1771 {
1772         return (struct isl_set *)
1773                 isl_map_intersect((struct isl_map *)set1,
1774                                   (struct isl_map *)set2);
1775 }
1776
1777 struct isl_basic_map *isl_basic_map_reverse(struct isl_basic_map *bmap)
1778 {
1779         struct isl_basic_set *bset;
1780         unsigned in;
1781
1782         if (!bmap)
1783                 return NULL;
1784         bmap = isl_basic_map_cow(bmap);
1785         if (!bmap)
1786                 return NULL;
1787         in = bmap->n_in;
1788         bset = isl_basic_set_from_basic_map(bmap);
1789         bset = isl_basic_set_swap_vars(bset, in);
1790         return isl_basic_map_from_basic_set(bset, bset->dim-in, in);
1791 }
1792
1793 /* Turn final n dimensions into existentially quantified variables.
1794  */
1795 struct isl_basic_set *isl_basic_set_project_out(
1796                 struct isl_basic_set *bset, unsigned n, unsigned flags)
1797 {
1798         int i;
1799         size_t row_size;
1800         isl_int **new_div;
1801         isl_int *old;
1802
1803         if (!bset)
1804                 return NULL;
1805
1806         isl_assert(bset->ctx, n <= bset->dim, goto error);
1807
1808         if (n == 0)
1809                 return bset;
1810
1811         bset = isl_basic_set_cow(bset);
1812
1813         row_size = 1 + bset->nparam + bset->dim + bset->extra;
1814         old = bset->block2.data;
1815         bset->block2 = isl_blk_extend(bset->ctx, bset->block2,
1816                                         (bset->extra + n) * (1 + row_size));
1817         if (!bset->block2.data)
1818                 goto error;
1819         new_div = isl_alloc_array(ctx, isl_int *, bset->extra + n);
1820         if (!new_div)
1821                 goto error;
1822         for (i = 0; i < n; ++i) {
1823                 new_div[i] = bset->block2.data +
1824                                 (bset->extra + i) * (1 + row_size);
1825                 isl_seq_clr(new_div[i], 1 + row_size);
1826         }
1827         for (i = 0; i < bset->extra; ++i)
1828                 new_div[n + i] = bset->block2.data + (bset->div[i] - old);
1829         free(bset->div);
1830         bset->div = new_div;
1831         bset->n_div += n;
1832         bset->extra += n;
1833         bset->dim -= n;
1834         bset = isl_basic_set_simplify(bset);
1835         return isl_basic_set_finalize(bset);
1836 error:
1837         isl_basic_set_free(bset);
1838         return NULL;
1839 }
1840
1841 struct isl_basic_map *isl_basic_map_apply_range(
1842                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
1843 {
1844         struct isl_basic_set *bset;
1845         unsigned n_in, n_out;
1846
1847         if (!bmap1 || !bmap2)
1848                 goto error;
1849
1850         isl_assert(bmap->ctx, bmap1->n_out == bmap2->n_in, goto error);
1851         isl_assert(bmap->ctx, bmap1->nparam == bmap2->nparam, goto error);
1852
1853         n_in = bmap1->n_in;
1854         n_out = bmap2->n_out;
1855
1856         bmap2 = isl_basic_map_reverse(bmap2);
1857         if (!bmap2)
1858                 goto error;
1859         bmap1 = isl_basic_map_extend(bmap1, bmap1->nparam,
1860                         bmap1->n_in + bmap2->n_in, bmap1->n_out,
1861                         bmap2->extra, bmap2->n_eq, bmap2->n_ineq);
1862         if (!bmap1)
1863                 goto error;
1864         bmap1 = add_constraints(bmap1, bmap2, bmap1->n_in - bmap2->n_in, 0);
1865         bmap1 = isl_basic_map_simplify(bmap1);
1866         bset = isl_basic_set_from_basic_map(bmap1);
1867         bset = isl_basic_set_project_out(bset,
1868                                                 bset->dim - (n_in + n_out), 0);
1869         return isl_basic_map_from_basic_set(bset, n_in, n_out);
1870 error:
1871         isl_basic_map_free(bmap1);
1872         isl_basic_map_free(bmap2);
1873         return NULL;
1874 }
1875
1876 struct isl_basic_set *isl_basic_set_apply(
1877                 struct isl_basic_set *bset, struct isl_basic_map *bmap)
1878 {
1879         if (!bset || !bmap)
1880                 goto error;
1881
1882         isl_assert(bset->ctx, bset->dim == bmap->n_in, goto error);
1883
1884         return (struct isl_basic_set *)
1885                 isl_basic_map_apply_range((struct isl_basic_map *)bset, bmap);
1886 error:
1887         isl_basic_set_free(bset);
1888         isl_basic_map_free(bmap);
1889         return NULL;
1890 }
1891
1892 struct isl_basic_map *isl_basic_map_apply_domain(
1893                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
1894 {
1895         if (!bmap1 || !bmap2)
1896                 goto error;
1897
1898         isl_assert(ctx, bmap1->n_in == bmap2->n_in, goto error);
1899         isl_assert(ctx, bmap1->nparam == bmap2->nparam, goto error);
1900
1901         bmap1 = isl_basic_map_reverse(bmap1);
1902         bmap1 = isl_basic_map_apply_range(bmap1, bmap2);
1903         return isl_basic_map_reverse(bmap1);
1904 error:
1905         isl_basic_map_free(bmap1);
1906         isl_basic_map_free(bmap2);
1907         return NULL;
1908 }
1909
1910 static struct isl_basic_map *var_equal(struct isl_ctx *ctx,
1911                 struct isl_basic_map *bmap, unsigned pos)
1912 {
1913         int i;
1914         i = isl_basic_map_alloc_equality(bmap);
1915         if (i < 0)
1916                 goto error;
1917         isl_seq_clr(bmap->eq[i],
1918                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
1919         isl_int_set_si(bmap->eq[i][1+bmap->nparam+pos], -1);
1920         isl_int_set_si(bmap->eq[i][1+bmap->nparam+bmap->n_in+pos], 1);
1921         return isl_basic_map_finalize(bmap);
1922 error:
1923         isl_basic_map_free(bmap);
1924         return NULL;
1925 }
1926
1927 static struct isl_basic_map *var_more(struct isl_ctx *ctx,
1928                 struct isl_basic_map *bmap, unsigned pos)
1929 {
1930         int i;
1931         i = isl_basic_map_alloc_inequality(bmap);
1932         if (i < 0)
1933                 goto error;
1934         isl_seq_clr(bmap->ineq[i],
1935                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
1936         isl_int_set_si(bmap->ineq[i][0], -1);
1937         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+pos], -1);
1938         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+bmap->n_in+pos], 1);
1939         return isl_basic_map_finalize(bmap);
1940 error:
1941         isl_basic_map_free(bmap);
1942         return NULL;
1943 }
1944
1945 static struct isl_basic_map *var_less(struct isl_ctx *ctx,
1946                 struct isl_basic_map *bmap, unsigned pos)
1947 {
1948         int i;
1949         i = isl_basic_map_alloc_inequality(bmap);
1950         if (i < 0)
1951                 goto error;
1952         isl_seq_clr(bmap->ineq[i],
1953                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
1954         isl_int_set_si(bmap->ineq[i][0], -1);
1955         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+pos], 1);
1956         isl_int_set_si(bmap->ineq[i][1+bmap->nparam+bmap->n_in+pos], -1);
1957         return isl_basic_map_finalize(bmap);
1958 error:
1959         isl_basic_map_free(bmap);
1960         return NULL;
1961 }
1962
1963 struct isl_basic_map *isl_basic_map_equal(struct isl_ctx *ctx,
1964                 unsigned nparam, unsigned in, unsigned out, unsigned n_equal)
1965 {
1966         int i;
1967         struct isl_basic_map *bmap;
1968         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, n_equal, 0);
1969         if (!bmap)
1970                 return NULL;
1971         for (i = 0; i < n_equal && bmap; ++i)
1972                 bmap = var_equal(ctx, bmap, i);
1973         return isl_basic_map_finalize(bmap);
1974 }
1975
1976 struct isl_basic_map *isl_basic_map_less_at(struct isl_ctx *ctx,
1977                 unsigned nparam, unsigned in, unsigned out, unsigned pos)
1978 {
1979         int i;
1980         struct isl_basic_map *bmap;
1981         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, pos, 1);
1982         if (!bmap)
1983                 return NULL;
1984         for (i = 0; i < pos && bmap; ++i)
1985                 bmap = var_equal(ctx, bmap, i);
1986         if (bmap)
1987                 bmap = var_less(ctx, bmap, pos);
1988         return isl_basic_map_finalize(bmap);
1989 }
1990
1991 struct isl_basic_map *isl_basic_map_more_at(struct isl_ctx *ctx,
1992                 unsigned nparam, unsigned in, unsigned out, unsigned pos)
1993 {
1994         int i;
1995         struct isl_basic_map *bmap;
1996         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, pos, 1);
1997         if (!bmap)
1998                 return NULL;
1999         for (i = 0; i < pos && bmap; ++i)
2000                 bmap = var_equal(ctx, bmap, i);
2001         if (bmap)
2002                 bmap = var_more(ctx, bmap, pos);
2003         return isl_basic_map_finalize(bmap);
2004 }
2005
2006 struct isl_basic_map *isl_basic_map_from_basic_set(
2007                 struct isl_basic_set *bset, unsigned n_in, unsigned n_out)
2008 {
2009         struct isl_basic_map *bmap;
2010
2011         bset = isl_basic_set_cow(bset);
2012         if (!bset)
2013                 return NULL;
2014
2015         isl_assert(ctx, bset->dim == n_in + n_out, goto error);
2016         bmap = (struct isl_basic_map *) bset;
2017         bmap->n_in = n_in;
2018         bmap->n_out = n_out;
2019         return isl_basic_map_finalize(bmap);
2020 error:
2021         isl_basic_set_free(bset);
2022         return NULL;
2023 }
2024
2025 struct isl_basic_set *isl_basic_set_from_basic_map(struct isl_basic_map *bmap)
2026 {
2027         if (!bmap)
2028                 goto error;
2029         if (bmap->n_in == 0)
2030                 return (struct isl_basic_set *)bmap;
2031         bmap = isl_basic_map_cow(bmap);
2032         if (!bmap)
2033                 goto error;
2034         bmap->n_out += bmap->n_in;
2035         bmap->n_in = 0;
2036         bmap = isl_basic_map_finalize(bmap);
2037         return (struct isl_basic_set *)bmap;
2038 error:
2039         return NULL;
2040 }
2041
2042 /* For a div d = floor(f/m), add the constraints
2043  *
2044  *              f - m d >= 0
2045  *              -(f-(n-1)) + m d >= 0
2046  *
2047  * Note that the second constraint is the negation of
2048  *
2049  *              f - m d >= n
2050  */
2051 static int add_div_constraints(struct isl_basic_map *bmap, unsigned div)
2052 {
2053         int i, j;
2054         unsigned div_pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
2055         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra;
2056
2057         i = isl_basic_map_alloc_inequality(bmap);
2058         if (i < 0)
2059                 return -1;
2060         isl_seq_cpy(bmap->ineq[i], bmap->div[div]+1, 1+total);
2061         isl_int_neg(bmap->ineq[i][div_pos], bmap->div[div][0]);
2062
2063         j = isl_basic_map_alloc_inequality(bmap);
2064         if (j < 0)
2065                 return -1;
2066         isl_seq_neg(bmap->ineq[j], bmap->ineq[i], 1 + total);
2067         isl_int_add(bmap->ineq[j][0], bmap->ineq[j][0], bmap->ineq[j][div_pos]);
2068         isl_int_sub_ui(bmap->ineq[j][0], bmap->ineq[j][0], 1);
2069         return j;
2070 }
2071
2072 struct isl_basic_set *isl_basic_map_underlying_set(
2073                 struct isl_basic_map *bmap)
2074 {
2075         if (!bmap)
2076                 goto error;
2077         if (bmap->nparam == 0 && bmap->n_in == 0 && bmap->n_div == 0)
2078                 return (struct isl_basic_set *)bmap;
2079         bmap = isl_basic_map_cow(bmap);
2080         if (!bmap)
2081                 goto error;
2082         bmap->n_out += bmap->nparam + bmap->n_in + bmap->n_div;
2083         bmap->nparam = 0;
2084         bmap->n_in = 0;
2085         bmap->extra -= bmap->n_div;
2086         bmap->n_div = 0;
2087         bmap = isl_basic_map_finalize(bmap);
2088         return (struct isl_basic_set *)bmap;
2089 error:
2090         return NULL;
2091 }
2092
2093 struct isl_basic_map *isl_basic_map_overlying_set(
2094         struct isl_basic_set *bset, struct isl_basic_map *like)
2095 {
2096         struct isl_basic_map *bmap;
2097         struct isl_ctx *ctx;
2098         unsigned total;
2099         int i, k;
2100
2101         if (!bset || !like)
2102                 goto error;
2103         ctx = bset->ctx;
2104         isl_assert(ctx, bset->dim ==
2105                         like->nparam + like->n_in + like->n_out + like->n_div,
2106                         goto error);
2107         if (like->nparam == 0 && like->n_in == 0 && like->n_div == 0) {
2108                 isl_basic_map_free(like);
2109                 return (struct isl_basic_map *)bset;
2110         }
2111         bset = isl_basic_set_cow(bset);
2112         if (!bset)
2113                 goto error;
2114         total = bset->dim + bset->extra;
2115         bmap = (struct isl_basic_map *)bset;
2116         bmap->nparam = like->nparam;
2117         bmap->n_in = like->n_in;
2118         bmap->n_out = like->n_out;
2119         bmap->extra += like->n_div;
2120         if (bmap->extra) {
2121                 unsigned ltotal;
2122                 ltotal = total - bmap->extra + like->extra;
2123                 if (ltotal > total)
2124                         ltotal = total;
2125                 bmap->block2 = isl_blk_extend(ctx, bmap->block2,
2126                                         bmap->extra * (1 + 1 + total));
2127                 if (isl_blk_is_error(bmap->block2))
2128                         goto error;
2129                 bmap->div = isl_realloc_array(ctx, bmap->div, isl_int *,
2130                                                 bmap->extra);
2131                 if (!bmap->div)
2132                         goto error;
2133                 bmap = isl_basic_map_extend(bmap, bmap->nparam,
2134                         bmap->n_in, bmap->n_out, 0, 0, 2 * like->n_div);
2135                 for (i = 0; i < like->n_div; ++i) {
2136                         k = isl_basic_map_alloc_div(bmap);
2137                         if (k < 0)
2138                                 goto error;
2139                         isl_seq_cpy(bmap->div[k], like->div[i], 1 + 1 + ltotal);
2140                         isl_seq_clr(bmap->div[k]+1+1+ltotal, total - ltotal);
2141                         if (add_div_constraints(bmap, k) < 0)
2142                                 goto error;
2143                 }
2144         }
2145         isl_basic_map_free(like);
2146         bmap = isl_basic_map_finalize(bmap);
2147         return bmap;
2148 error:
2149         isl_basic_map_free(like);
2150         isl_basic_set_free(bset);
2151         return NULL;
2152 }
2153
2154 struct isl_basic_set *isl_basic_set_from_underlying_set(
2155         struct isl_basic_set *bset, struct isl_basic_set *like)
2156 {
2157         return (struct isl_basic_set *)
2158                 isl_basic_map_overlying_set(bset, (struct isl_basic_map *)like);
2159 }
2160
2161 struct isl_set *isl_set_from_underlying_set(
2162         struct isl_set *set, struct isl_basic_set *like)
2163 {
2164         int i;
2165
2166         if (!set || !like)
2167                 goto error;
2168         isl_assert(set->ctx, set->dim == like->nparam + like->dim + like->n_div,
2169                     goto error);
2170         if (like->nparam == 0 && like->n_div == 0) {
2171                 isl_basic_set_free(like);
2172                 return set;
2173         }
2174         set = isl_set_cow(set);
2175         if (!set)
2176                 goto error;
2177         for (i = 0; i < set->n; ++i) {
2178                 set->p[i] = isl_basic_set_from_underlying_set(set->p[i],
2179                                                       isl_basic_set_copy(like));
2180                 if (!set->p[i])
2181                         goto error;
2182         }
2183         set->nparam = like->nparam;
2184         set->dim = like->dim;
2185         isl_basic_set_free(like);
2186         return set;
2187 error:
2188         isl_basic_set_free(like);
2189         isl_set_free(set);
2190         return NULL;
2191 }
2192
2193 struct isl_set *isl_map_underlying_set(struct isl_map *map)
2194 {
2195         int i;
2196
2197         map = isl_map_align_divs(map);
2198         map = isl_map_cow(map);
2199         if (!map)
2200                 return NULL;
2201
2202         for (i = 0; i < map->n; ++i) {
2203                 map->p[i] = (struct isl_basic_map *)
2204                                 isl_basic_map_underlying_set(map->p[i]);
2205                 if (!map->p[i])
2206                         goto error;
2207         }
2208         if (map->n == 0)
2209                 map->n_out += map->nparam + map->n_in;
2210         else
2211                 map->n_out = map->p[0]->n_out;
2212         map->nparam = 0;
2213         map->n_in = 0;
2214         return (struct isl_set *)map;
2215 error:
2216         isl_map_free(map);
2217         return NULL;
2218 }
2219
2220 struct isl_set *isl_set_to_underlying_set(struct isl_set *set)
2221 {
2222         return (struct isl_set *)isl_map_underlying_set((struct isl_map *)set);
2223 }
2224
2225 struct isl_basic_set *isl_basic_map_domain(struct isl_basic_map *bmap)
2226 {
2227         struct isl_basic_set *domain;
2228         unsigned n_out;
2229         if (!bmap)
2230                 return NULL;
2231         n_out = bmap->n_out;
2232         domain = isl_basic_set_from_basic_map(bmap);
2233         return isl_basic_set_project_out(domain, n_out, 0);
2234 }
2235
2236 struct isl_basic_set *isl_basic_map_range(struct isl_basic_map *bmap)
2237 {
2238         return isl_basic_map_domain(isl_basic_map_reverse(bmap));
2239 }
2240
2241 struct isl_set *isl_map_range(struct isl_map *map)
2242 {
2243         int i;
2244         struct isl_set *set;
2245
2246         if (!map)
2247                 goto error;
2248         map = isl_map_cow(map);
2249         if (!map)
2250                 goto error;
2251
2252         set = (struct isl_set *) map;
2253         set->zero = 0;
2254         for (i = 0; i < map->n; ++i) {
2255                 set->p[i] = isl_basic_map_range(map->p[i]);
2256                 if (!set->p[i])
2257                         goto error;
2258         }
2259         F_CLR(set, ISL_MAP_DISJOINT);
2260         return set;
2261 error:
2262         isl_map_free(map);
2263         return NULL;
2264 }
2265
2266 struct isl_map *isl_map_from_set(struct isl_set *set,
2267                 unsigned n_in, unsigned n_out)
2268 {
2269         int i;
2270         struct isl_map *map = NULL;
2271
2272         if (!set)
2273                 return NULL;
2274         isl_assert(set->ctx, set->dim == n_in + n_out, goto error);
2275         set = isl_set_cow(set);
2276         if (!set)
2277                 return NULL;
2278         map = (struct isl_map *)set;
2279         for (i = 0; i < set->n; ++i) {
2280                 map->p[i] = isl_basic_map_from_basic_set(
2281                                 set->p[i], n_in, n_out);
2282                 if (!map->p[i])
2283                         goto error;
2284         }
2285         map->n_in = n_in;
2286         map->n_out = n_out;
2287         return map;
2288 error:
2289         isl_set_free(set);
2290         return NULL;
2291 }
2292
2293 struct isl_set *isl_set_from_map(struct isl_map *map)
2294 {
2295         int i;
2296         struct isl_set *set = NULL;
2297
2298         if (!map)
2299                 return NULL;
2300         map = isl_map_cow(map);
2301         if (!map)
2302                 return NULL;
2303         map->n_out += map->n_in;
2304         map->n_in = 0;
2305         set = (struct isl_set *)map;
2306         for (i = 0; i < map->n; ++i) {
2307                 set->p[i] = isl_basic_set_from_basic_map(map->p[i]);
2308                 if (!set->p[i])
2309                         goto error;
2310         }
2311         return set;
2312 error:
2313         isl_map_free(map);
2314         return NULL;
2315 }
2316
2317 struct isl_map *isl_map_alloc(struct isl_ctx *ctx,
2318                 unsigned nparam, unsigned in, unsigned out, int n,
2319                 unsigned flags)
2320 {
2321         struct isl_map *map;
2322
2323         map = isl_alloc(ctx, struct isl_map,
2324                         sizeof(struct isl_map) +
2325                         n * sizeof(struct isl_basic_map *));
2326         if (!map)
2327                 return NULL;
2328
2329         map->ctx = ctx;
2330         isl_ctx_ref(ctx);
2331         map->ref = 1;
2332         map->size = n;
2333         map->n = 0;
2334         map->nparam = nparam;
2335         map->n_in = in;
2336         map->n_out = out;
2337         map->flags = flags;
2338         return map;
2339 }
2340
2341 struct isl_basic_map *isl_basic_map_empty(struct isl_ctx *ctx,
2342                 unsigned nparam, unsigned in, unsigned out)
2343 {
2344         struct isl_basic_map *bmap;
2345         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, 1, 0);
2346         bmap = isl_basic_map_set_to_empty(bmap);
2347         return bmap;
2348 }
2349
2350 struct isl_basic_set *isl_basic_set_empty(struct isl_ctx *ctx,
2351                 unsigned nparam, unsigned dim)
2352 {
2353         struct isl_basic_set *bset;
2354         bset = isl_basic_set_alloc(ctx, nparam, dim, 0, 1, 0);
2355         bset = isl_basic_set_set_to_empty(bset);
2356         return bset;
2357 }
2358
2359 struct isl_basic_map *isl_basic_map_universe(struct isl_ctx *ctx,
2360                 unsigned nparam, unsigned in, unsigned out)
2361 {
2362         struct isl_basic_map *bmap;
2363         bmap = isl_basic_map_alloc(ctx, nparam, in, out, 0, 0, 0);
2364         return bmap;
2365 }
2366
2367 struct isl_basic_set *isl_basic_set_universe(struct isl_ctx *ctx,
2368                 unsigned nparam, unsigned dim)
2369 {
2370         struct isl_basic_set *bset;
2371         bset = isl_basic_set_alloc(ctx, nparam, dim, 0, 0, 0);
2372         return bset;
2373 }
2374
2375 struct isl_map *isl_map_empty(struct isl_ctx *ctx,
2376                 unsigned nparam, unsigned in, unsigned out)
2377 {
2378         return isl_map_alloc(ctx, nparam, in, out, 0, ISL_MAP_DISJOINT);
2379 }
2380
2381 struct isl_set *isl_set_empty(struct isl_ctx *ctx,
2382                 unsigned nparam, unsigned dim)
2383 {
2384         return isl_set_alloc(ctx, nparam, dim, 0, ISL_MAP_DISJOINT);
2385 }
2386
2387 struct isl_map *isl_map_dup(struct isl_map *map)
2388 {
2389         int i;
2390         struct isl_map *dup;
2391
2392         if (!map)
2393                 return NULL;
2394         dup = isl_map_alloc(map->ctx, map->nparam, map->n_in, map->n_out, map->n,
2395                                 map->flags);
2396         for (i = 0; i < map->n; ++i)
2397                 dup = isl_map_add(dup,
2398                                   isl_basic_map_dup(map->p[i]));
2399         return dup;
2400 }
2401
2402 struct isl_map *isl_map_add(struct isl_map *map, struct isl_basic_map *bmap)
2403 {
2404         if (!bmap || !map)
2405                 goto error;
2406         isl_assert(map->ctx, map->nparam == bmap->nparam, goto error);
2407         isl_assert(map->ctx, map->n_in == bmap->n_in, goto error);
2408         isl_assert(map->ctx, map->n_out == bmap->n_out, goto error);
2409         isl_assert(map->ctx, map->n < map->size, goto error);
2410         map->p[map->n] = bmap;
2411         map->n++;
2412         return map;
2413 error:
2414         if (map)
2415                 isl_map_free(map);
2416         if (bmap)
2417                 isl_basic_map_free(bmap);
2418         return NULL;
2419 }
2420
2421 void isl_map_free(struct isl_map *map)
2422 {
2423         int i;
2424
2425         if (!map)
2426                 return;
2427
2428         if (--map->ref > 0)
2429                 return;
2430
2431         isl_ctx_deref(map->ctx);
2432         for (i = 0; i < map->n; ++i)
2433                 isl_basic_map_free(map->p[i]);
2434         free(map);
2435 }
2436
2437 struct isl_map *isl_map_extend(struct isl_map *base,
2438                 unsigned nparam, unsigned n_in, unsigned n_out)
2439 {
2440         int i;
2441
2442         base = isl_map_cow(base);
2443         if (!base)
2444                 return NULL;
2445
2446         isl_assert(base->ctx, base->nparam <= nparam, goto error);
2447         isl_assert(base->ctx, base->n_in <= n_in, goto error);
2448         isl_assert(base->ctx, base->n_out <= n_out, goto error);
2449         base->nparam = nparam;
2450         base->n_in = n_in;
2451         base->n_out = n_out;
2452         for (i = 0; i < base->n; ++i) {
2453                 base->p[i] = isl_basic_map_extend(base->p[i],
2454                                 nparam, n_in, n_out, 0, 0, 0);
2455                 if (!base->p[i])
2456                         goto error;
2457         }
2458         return base;
2459 error:
2460         isl_map_free(base);
2461         return NULL;
2462 }
2463
2464 struct isl_set *isl_set_extend(struct isl_set *base,
2465                 unsigned nparam, unsigned dim)
2466 {
2467         return (struct isl_set *)isl_map_extend((struct isl_map *)base,
2468                                                         nparam, 0, dim);
2469 }
2470
2471 static struct isl_basic_map *isl_basic_map_fix_var(struct isl_basic_map *bmap,
2472                 unsigned var, int value)
2473 {
2474         int j;
2475
2476         bmap = isl_basic_map_cow(bmap);
2477         bmap = isl_basic_map_extend(bmap,
2478                         bmap->nparam, bmap->n_in, bmap->n_out, 0, 1, 0);
2479         j = isl_basic_map_alloc_equality(bmap);
2480         if (j < 0)
2481                 goto error;
2482         isl_seq_clr(bmap->eq[j],
2483                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
2484         isl_int_set_si(bmap->eq[j][1 + var], -1);
2485         isl_int_set_si(bmap->eq[j][0], value);
2486         bmap = isl_basic_map_simplify(bmap);
2487         return isl_basic_map_finalize(bmap);
2488 error:
2489         isl_basic_map_free(bmap);
2490         return NULL;
2491 }
2492
2493 struct isl_basic_map *isl_basic_map_fix_input_si(struct isl_basic_map *bmap,
2494                 unsigned input, int value)
2495 {
2496         if (!bmap)
2497                 return NULL;
2498         isl_assert(bmap->ctx, input < bmap->n_in, goto error);
2499         return isl_basic_map_fix_var(bmap, bmap->nparam + input, value);
2500 error:
2501         isl_basic_map_free(bmap);
2502         return NULL;
2503 }
2504
2505 struct isl_basic_set *isl_basic_set_fix_dim_si(struct isl_basic_set *bset,
2506                 unsigned dim, int value)
2507 {
2508         if (!bset)
2509                 return NULL;
2510         isl_assert(bset->ctx, dim < bset->dim, goto error);
2511         return (struct isl_basic_set *)
2512                 isl_basic_map_fix_var((struct isl_basic_map *)bset,
2513                                                 bset->nparam + dim, value);
2514 error:
2515         isl_basic_set_free(bset);
2516         return NULL;
2517 }
2518
2519 struct isl_map *isl_map_fix_input_si(struct isl_map *map,
2520                 unsigned input, int value)
2521 {
2522         int i;
2523
2524         map = isl_map_cow(map);
2525         if (!map)
2526                 return NULL;
2527
2528         isl_assert(ctx, input < map->n_in, goto error);
2529         for (i = 0; i < map->n; ++i) {
2530                 map->p[i] = isl_basic_map_fix_input_si(map->p[i], input, value);
2531                 if (!map->p[i])
2532                         goto error;
2533         }
2534         return map;
2535 error:
2536         isl_map_free(map);
2537         return NULL;
2538 }
2539
2540 struct isl_set *isl_set_fix_dim_si(struct isl_set *set, unsigned dim, int value)
2541 {
2542         int i;
2543
2544         set = isl_set_cow(set);
2545         if (!set)
2546                 return NULL;
2547
2548         isl_assert(set->ctx, dim < set->dim, goto error);
2549         for (i = 0; i < set->n; ++i) {
2550                 set->p[i] = isl_basic_set_fix_dim_si(set->p[i], dim, value);
2551                 if (!set->p[i])
2552                         goto error;
2553         }
2554         return set;
2555 error:
2556         isl_set_free(set);
2557         return NULL;
2558 }
2559
2560 struct isl_map *isl_map_reverse(struct isl_map *map)
2561 {
2562         int i;
2563         unsigned t;
2564
2565         map = isl_map_cow(map);
2566         if (!map)
2567                 return NULL;
2568
2569         t = map->n_in;
2570         map->n_in = map->n_out;
2571         map->n_out = t;
2572         for (i = 0; i < map->n; ++i) {
2573                 map->p[i] = isl_basic_map_reverse(map->p[i]);
2574                 if (!map->p[i])
2575                         goto error;
2576         }
2577         return map;
2578 error:
2579         isl_map_free(map);
2580         return NULL;
2581 }
2582
2583 struct isl_map *isl_basic_map_lexmax(
2584                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
2585                 struct isl_set **empty)
2586 {
2587         return isl_pip_basic_map_lexmax(bmap, dom, empty);
2588 }
2589
2590 struct isl_map *isl_basic_map_lexmin(
2591                 struct isl_basic_map *bmap, struct isl_basic_set *dom,
2592                 struct isl_set **empty)
2593 {
2594         return isl_pip_basic_map_lexmin(bmap, dom, empty);
2595 }
2596
2597 struct isl_set *isl_basic_set_lexmin(struct isl_basic_set *bset)
2598 {
2599         struct isl_basic_map *bmap = NULL;
2600         struct isl_basic_set *dom = NULL;
2601         struct isl_map *min;
2602
2603         if (!bset)
2604                 goto error;
2605         bmap = isl_basic_map_from_basic_set(bset, 0, bset->dim);
2606         if (!bmap)
2607                 goto error;
2608         dom = isl_basic_set_alloc(bmap->ctx, bmap->nparam, 0, 0, 0, 0);
2609         if (!dom)
2610                 goto error;
2611         min = isl_basic_map_lexmin(bmap, dom, NULL);
2612         return isl_map_range(min);
2613 error:
2614         isl_basic_map_free(bmap);
2615         return NULL;
2616 }
2617
2618 struct isl_map *isl_basic_map_compute_divs(struct isl_basic_map *bmap)
2619 {
2620         if (!bmap)
2621                 return NULL;
2622         if (bmap->n_div == 0)
2623                 return isl_map_from_basic_map(bmap);
2624         return isl_pip_basic_map_compute_divs(bmap);
2625 }
2626
2627 struct isl_map *isl_map_compute_divs(struct isl_map *map)
2628 {
2629         int i;
2630         struct isl_map *res;
2631
2632         if (!map)
2633                 return NULL;
2634         if (map->n == 0)
2635                 return map;
2636         res = isl_basic_map_compute_divs(isl_basic_map_copy(map->p[0]));
2637         for (i = 1 ; i < map->n; ++i) {
2638                 struct isl_map *r2;
2639                 r2 = isl_basic_map_compute_divs(isl_basic_map_copy(map->p[i]));
2640                 if (F_ISSET(map, ISL_MAP_DISJOINT))
2641                         res = isl_map_union_disjoint(res, r2);
2642                 else
2643                         res = isl_map_union(res, r2);
2644         }
2645         isl_map_free(map);
2646
2647         return res;
2648 }
2649
2650 struct isl_set *isl_basic_set_compute_divs(struct isl_basic_set *bset)
2651 {
2652         return (struct isl_set *)
2653                 isl_basic_map_compute_divs((struct isl_basic_map *)bset);
2654 }
2655
2656 struct isl_set *isl_set_compute_divs(struct isl_set *set)
2657 {
2658         return (struct isl_set *)
2659                 isl_map_compute_divs((struct isl_map *)set);
2660 }
2661
2662 struct isl_set *isl_map_domain(struct isl_map *map)
2663 {
2664         int i;
2665         struct isl_set *set;
2666
2667         if (!map)
2668                 goto error;
2669
2670         map = isl_map_cow(map);
2671         if (!map)
2672                 return NULL;
2673
2674         set = (struct isl_set *)map;
2675         set->dim = map->n_in;
2676         set->zero = 0;
2677         for (i = 0; i < map->n; ++i) {
2678                 set->p[i] = isl_basic_map_domain(map->p[i]);
2679                 if (!set->p[i])
2680                         goto error;
2681         }
2682         F_CLR(set, ISL_MAP_DISJOINT);
2683         return set;
2684 error:
2685         isl_map_free(map);
2686         return NULL;
2687 }
2688
2689 struct isl_map *isl_map_union_disjoint(
2690                         struct isl_map *map1, struct isl_map *map2)
2691 {
2692         int i;
2693         unsigned flags = 0;
2694         struct isl_map *map = NULL;
2695
2696         if (!map1 || !map2)
2697                 goto error;
2698
2699         if (map1->n == 0) {
2700                 isl_map_free(map1);
2701                 return map2;
2702         }
2703         if (map2->n == 0) {
2704                 isl_map_free(map2);
2705                 return map1;
2706         }
2707
2708         isl_assert(map1->ctx, map1->nparam == map2->nparam, goto error);
2709         isl_assert(map1->ctx, map1->n_in == map2->n_in, goto error);
2710         isl_assert(map1->ctx, map1->n_out == map2->n_out, goto error);
2711
2712         if (F_ISSET(map1, ISL_MAP_DISJOINT) &&
2713             F_ISSET(map2, ISL_MAP_DISJOINT))
2714                 FL_SET(flags, ISL_MAP_DISJOINT);
2715
2716         map = isl_map_alloc(map1->ctx, map1->nparam, map1->n_in, map1->n_out,
2717                                 map1->n + map2->n, flags);
2718         if (!map)
2719                 goto error;
2720         for (i = 0; i < map1->n; ++i) {
2721                 map = isl_map_add(map,
2722                                   isl_basic_map_copy(map1->p[i]));
2723                 if (!map)
2724                         goto error;
2725         }
2726         for (i = 0; i < map2->n; ++i) {
2727                 map = isl_map_add(map,
2728                                   isl_basic_map_copy(map2->p[i]));
2729                 if (!map)
2730                         goto error;
2731         }
2732         isl_map_free(map1);
2733         isl_map_free(map2);
2734         return map;
2735 error:
2736         isl_map_free(map);
2737         isl_map_free(map1);
2738         isl_map_free(map2);
2739         return NULL;
2740 }
2741
2742 struct isl_map *isl_map_union(struct isl_map *map1, struct isl_map *map2)
2743 {
2744         map1 = isl_map_union_disjoint(map1, map2);
2745         if (!map1)
2746                 return NULL;
2747         if (map1->n > 1)
2748                 F_CLR(map1, ISL_MAP_DISJOINT);
2749         return map1;
2750 }
2751
2752 struct isl_set *isl_set_union_disjoint(
2753                         struct isl_set *set1, struct isl_set *set2)
2754 {
2755         return (struct isl_set *)
2756                 isl_map_union_disjoint(
2757                         (struct isl_map *)set1, (struct isl_map *)set2);
2758 }
2759
2760 struct isl_set *isl_set_union(struct isl_set *set1, struct isl_set *set2)
2761 {
2762         return (struct isl_set *)
2763                 isl_map_union((struct isl_map *)set1, (struct isl_map *)set2);
2764 }
2765
2766 struct isl_map *isl_map_intersect_range(
2767                 struct isl_map *map, struct isl_set *set)
2768 {
2769         unsigned flags = 0;
2770         struct isl_map *result;
2771         int i, j;
2772
2773         if (!map || !set)
2774                 goto error;
2775
2776         if (F_ISSET(map, ISL_MAP_DISJOINT) &&
2777             F_ISSET(set, ISL_MAP_DISJOINT))
2778                 FL_SET(flags, ISL_MAP_DISJOINT);
2779
2780         result = isl_map_alloc(map->ctx, map->nparam, map->n_in, map->n_out,
2781                                 map->n * set->n, flags);
2782         if (!result)
2783                 goto error;
2784         for (i = 0; i < map->n; ++i)
2785                 for (j = 0; j < set->n; ++j) {
2786                         result = isl_map_add(result,
2787                             isl_basic_map_intersect_range(
2788                                 isl_basic_map_copy(map->p[i]),
2789                                 isl_basic_set_copy(set->p[j])));
2790                         if (!result)
2791                                 goto error;
2792                 }
2793         isl_map_free(map);
2794         isl_set_free(set);
2795         return result;
2796 error:
2797         isl_map_free(map);
2798         isl_set_free(set);
2799         return NULL;
2800 }
2801
2802 struct isl_map *isl_map_intersect_domain(
2803                 struct isl_map *map, struct isl_set *set)
2804 {
2805         return isl_map_reverse(
2806                 isl_map_intersect_range(isl_map_reverse(map), set));
2807 }
2808
2809 struct isl_map *isl_map_apply_domain(
2810                 struct isl_map *map1, struct isl_map *map2)
2811 {
2812         if (!map1 || !map2)
2813                 goto error;
2814         map1 = isl_map_reverse(map1);
2815         map1 = isl_map_apply_range(map1, map2);
2816         return isl_map_reverse(map1);
2817 error:
2818         isl_map_free(map1);
2819         isl_map_free(map2);
2820         return NULL;
2821 }
2822
2823 struct isl_map *isl_map_apply_range(
2824                 struct isl_map *map1, struct isl_map *map2)
2825 {
2826         struct isl_map *result;
2827         int i, j;
2828
2829         if (!map1 || !map2)
2830                 goto error;
2831
2832         isl_assert(map1->ctx, map1->nparam == map2->nparam, goto error);
2833         isl_assert(map1->ctx, map1->n_out == map2->n_in, goto error);
2834
2835         result = isl_map_alloc(map1->ctx, map1->nparam, map1->n_in, map2->n_out,
2836                                 map1->n * map2->n, 0);
2837         if (!result)
2838                 goto error;
2839         for (i = 0; i < map1->n; ++i)
2840                 for (j = 0; j < map2->n; ++j) {
2841                         result = isl_map_add(result,
2842                             isl_basic_map_apply_range(
2843                                 isl_basic_map_copy(map1->p[i]),
2844                                 isl_basic_map_copy(map2->p[j])));
2845                         if (!result)
2846                                 goto error;
2847                 }
2848         isl_map_free(map1);
2849         isl_map_free(map2);
2850         if (result && result->n <= 1)
2851                 F_SET(result, ISL_MAP_DISJOINT);
2852         return result;
2853 error:
2854         isl_map_free(map1);
2855         isl_map_free(map2);
2856         return NULL;
2857 }
2858
2859 /*
2860  * returns range - domain
2861  */
2862 struct isl_basic_set *isl_basic_map_deltas(struct isl_basic_map *bmap)
2863 {
2864         struct isl_basic_set *bset;
2865         unsigned dim;
2866         int i;
2867
2868         if (!bmap)
2869                 goto error;
2870         isl_assert(bmap->ctx, bmap->n_in == bmap->n_out, goto error);
2871         dim = bmap->n_in;
2872         bset = isl_basic_set_from_basic_map(bmap);
2873         bset = isl_basic_set_extend(bset, bset->nparam, 3*dim, 0,
2874                                         dim, 0);
2875         bset = isl_basic_set_swap_vars(bset, 2*dim);
2876         for (i = 0; i < dim; ++i) {
2877                 int j = isl_basic_map_alloc_equality(
2878                                             (struct isl_basic_map *)bset);
2879                 if (j < 0)
2880                         goto error;
2881                 isl_seq_clr(bset->eq[j],
2882                             1 + bset->nparam + bset->dim + bset->extra);
2883                 isl_int_set_si(bset->eq[j][1+bset->nparam+i], 1);
2884                 isl_int_set_si(bset->eq[j][1+bset->nparam+dim+i], 1);
2885                 isl_int_set_si(bset->eq[j][1+bset->nparam+2*dim+i], -1);
2886         }
2887         return isl_basic_set_project_out(bset, 2*dim, 0);
2888 error:
2889         isl_basic_map_free(bmap);
2890         return NULL;
2891 }
2892
2893 /*
2894  * returns range - domain
2895  */
2896 struct isl_set *isl_map_deltas(struct isl_map *map)
2897 {
2898         int i;
2899         struct isl_set *result;
2900
2901         if (!map)
2902                 return NULL;
2903
2904         isl_assert(map->ctx, map->n_in == map->n_out, goto error);
2905         result = isl_set_alloc(map->ctx, map->nparam, map->n_in, map->n, map->flags);
2906         if (!result)
2907                 goto error;
2908         for (i = 0; i < map->n; ++i)
2909                 result = isl_set_add(result,
2910                           isl_basic_map_deltas(isl_basic_map_copy(map->p[i])));
2911         isl_map_free(map);
2912         return result;
2913 error:
2914         isl_map_free(map);
2915         return NULL;
2916 }
2917
2918 /* If the only constraints a div d=floor(f/m)
2919  * appears in are its two defining constraints
2920  *
2921  *      f - m d >=0
2922  *      -(f - (m - 1)) + m d >= 0
2923  *
2924  * then it can safely be removed.
2925  */
2926 static int div_is_redundant(struct isl_basic_map *bmap, int div)
2927 {
2928         int i;
2929         unsigned pos = 1 + bmap->nparam + bmap->n_in + bmap->n_out + div;
2930
2931         for (i = 0; i < bmap->n_eq; ++i)
2932                 if (!isl_int_is_zero(bmap->eq[i][pos]))
2933                         return 0;
2934
2935         for (i = 0; i < bmap->n_ineq; ++i) {
2936                 if (isl_int_is_zero(bmap->ineq[i][pos]))
2937                         continue;
2938                 if (isl_int_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
2939                         int neg;
2940                         isl_int_sub(bmap->div[div][1],
2941                                         bmap->div[div][1], bmap->div[div][0]);
2942                         isl_int_add_ui(bmap->div[div][1], bmap->div[div][1], 1);
2943                         neg = isl_seq_is_neg(bmap->ineq[i], bmap->div[div]+1, pos);
2944                         isl_int_sub_ui(bmap->div[div][1], bmap->div[div][1], 1);
2945                         isl_int_add(bmap->div[div][1],
2946                                         bmap->div[div][1], bmap->div[div][0]);
2947                         if (!neg)
2948                                 return 0;
2949                         if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
2950                                                     bmap->n_div-div-1) != -1)
2951                                 return 0;
2952                 } else if (isl_int_abs_eq(bmap->ineq[i][pos], bmap->div[div][0])) {
2953                         if (!isl_seq_eq(bmap->ineq[i], bmap->div[div]+1, pos))
2954                                 return 0;
2955                         if (isl_seq_first_non_zero(bmap->ineq[i]+pos+1,
2956                                                     bmap->n_div-div-1) != -1)
2957                                 return 0;
2958                 } else
2959                         return 0;
2960         }
2961
2962         for (i = 0; i < bmap->n_div; ++i)
2963                 if (!isl_int_is_zero(bmap->div[i][1+pos]))
2964                         return 0;
2965
2966         return 1;
2967 }
2968
2969 /*
2970  * Remove divs that don't occur in any of the constraints or other divs.
2971  * These can arise when dropping some of the variables in a quast
2972  * returned by piplib.
2973  */
2974 static struct isl_basic_map *remove_redundant_divs(struct isl_basic_map *bmap)
2975 {
2976         int i;
2977
2978         if (!bmap)
2979                 return NULL;
2980
2981         for (i = bmap->n_div-1; i >= 0; --i) {
2982                 if (!div_is_redundant(bmap, i))
2983                         continue;
2984                 bmap = isl_basic_map_drop_div(bmap, i);
2985         }
2986         return bmap;
2987 }
2988
2989 struct isl_basic_map *isl_basic_map_finalize(struct isl_basic_map *bmap)
2990 {
2991         bmap = remove_redundant_divs(bmap);
2992         if (!bmap)
2993                 return NULL;
2994         F_SET(bmap, ISL_BASIC_SET_FINAL);
2995         return bmap;
2996 }
2997
2998 struct isl_basic_set *isl_basic_set_finalize(struct isl_basic_set *bset)
2999 {
3000         return (struct isl_basic_set *)
3001                 isl_basic_map_finalize((struct isl_basic_map *)bset);
3002 }
3003
3004 struct isl_set *isl_set_finalize(struct isl_set *set)
3005 {
3006         int i;
3007
3008         if (!set)
3009                 return NULL;
3010         for (i = 0; i < set->n; ++i) {
3011                 set->p[i] = isl_basic_set_finalize(set->p[i]);
3012                 if (!set->p[i])
3013                         goto error;
3014         }
3015         return set;
3016 error:
3017         isl_set_free(set);
3018         return NULL;
3019 }
3020
3021 struct isl_map *isl_map_finalize(struct isl_map *map)
3022 {
3023         int i;
3024
3025         if (!map)
3026                 return NULL;
3027         for (i = 0; i < map->n; ++i) {
3028                 map->p[i] = isl_basic_map_finalize(map->p[i]);
3029                 if (!map->p[i])
3030                         goto error;
3031         }
3032         return map;
3033 error:
3034         isl_map_free(map);
3035         return NULL;
3036 }
3037
3038 struct isl_basic_map *isl_basic_map_identity(struct isl_ctx *ctx,
3039                 unsigned nparam, unsigned dim)
3040 {
3041         struct isl_basic_map *bmap;
3042         int i;
3043
3044         bmap = isl_basic_map_alloc(ctx, nparam, dim, dim, 0, dim, 0);
3045         if (!bmap)
3046                 goto error;
3047
3048         for (i = 0; i < dim; ++i) {
3049                 int j = isl_basic_map_alloc_equality(bmap);
3050                 if (j < 0)
3051                         goto error;
3052                 isl_seq_clr(bmap->eq[j],
3053                     1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->extra);
3054                 isl_int_set_si(bmap->eq[j][1+nparam+i], 1);
3055                 isl_int_set_si(bmap->eq[j][1+nparam+dim+i], -1);
3056         }
3057         return isl_basic_map_finalize(bmap);
3058 error:
3059         isl_basic_map_free(bmap);
3060         return NULL;
3061 }
3062
3063 struct isl_map *isl_map_identity(struct isl_ctx *ctx,
3064                 unsigned nparam, unsigned dim)
3065 {
3066         struct isl_map *map = isl_map_alloc(ctx, nparam, dim, dim, 1,
3067                                                 ISL_MAP_DISJOINT);
3068         if (!map)
3069                 goto error;
3070         map = isl_map_add(map,
3071                 isl_basic_map_identity(ctx, nparam, dim));
3072         return map;
3073 error:
3074         isl_map_free(map);
3075         return NULL;
3076 }
3077
3078 int isl_set_is_equal(struct isl_set *set1, struct isl_set *set2)
3079 {
3080         return isl_map_is_equal((struct isl_map *)set1, (struct isl_map *)set2);
3081 }
3082
3083 int isl_set_is_subset(struct isl_set *set1, struct isl_set *set2)
3084 {
3085         return isl_map_is_subset(
3086                         (struct isl_map *)set1, (struct isl_map *)set2);
3087 }
3088
3089 int isl_basic_map_is_subset(
3090                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3091 {
3092         int is_subset;
3093         struct isl_map *map1;
3094         struct isl_map *map2;
3095
3096         if (!bmap1 || !bmap2)
3097                 return -1;
3098
3099         map1 = isl_map_from_basic_map(isl_basic_map_copy(bmap1));
3100         map2 = isl_map_from_basic_map(isl_basic_map_copy(bmap2));
3101
3102         is_subset = isl_map_is_subset(map1, map2);
3103
3104         isl_map_free(map1);
3105         isl_map_free(map2);
3106
3107         return is_subset;
3108 }
3109
3110 int isl_basic_map_is_equal(
3111                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3112 {
3113         int is_subset;
3114
3115         if (!bmap1 || !bmap2)
3116                 return -1;
3117         is_subset = isl_basic_map_is_subset(bmap1, bmap2);
3118         if (is_subset != 1)
3119                 return is_subset;
3120         is_subset = isl_basic_map_is_subset(bmap2, bmap1);
3121         return is_subset;
3122 }
3123
3124 int isl_basic_set_is_equal(
3125                 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
3126 {
3127         return isl_basic_map_is_equal(
3128                 (struct isl_basic_map *)bset1, (struct isl_basic_map *)bset2);
3129 }
3130
3131 int isl_map_is_empty(struct isl_map *map)
3132 {
3133         int i;
3134         int is_empty;
3135
3136         if (!map)
3137                 return -1;
3138         for (i = 0; i < map->n; ++i) {
3139                 is_empty = isl_basic_map_is_empty(map->p[i]);
3140                 if (is_empty < 0)
3141                         return -1;
3142                 if (!is_empty)
3143                         return 0;
3144         }
3145         return 1;
3146 }
3147
3148 int isl_set_is_empty(struct isl_set *set)
3149 {
3150         return isl_map_is_empty((struct isl_map *)set);
3151 }
3152
3153 int isl_map_is_subset(struct isl_map *map1, struct isl_map *map2)
3154 {
3155         int i;
3156         int is_subset = 0;
3157         struct isl_map *diff;
3158
3159         if (!map1 || !map2)
3160                 return -1;
3161
3162         if (isl_map_is_empty(map1))
3163                 return 1;
3164
3165         if (isl_map_is_empty(map2))
3166                 return 0;
3167
3168         diff = isl_map_subtract(isl_map_copy(map1), isl_map_copy(map2));
3169         if (!diff)
3170                 return -1;
3171
3172         is_subset = isl_map_is_empty(diff);
3173         isl_map_free(diff);
3174
3175         return is_subset;
3176 }
3177
3178 int isl_map_is_equal(struct isl_map *map1, struct isl_map *map2)
3179 {
3180         int is_subset;
3181
3182         if (!map1 || !map2)
3183                 return -1;
3184         is_subset = isl_map_is_subset(map1, map2);
3185         if (is_subset != 1)
3186                 return is_subset;
3187         is_subset = isl_map_is_subset(map2, map1);
3188         return is_subset;
3189 }
3190
3191 int isl_basic_map_is_strict_subset(
3192                 struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3193 {
3194         int is_subset;
3195
3196         if (!bmap1 || !bmap2)
3197                 return -1;
3198         is_subset = isl_basic_map_is_subset(bmap1, bmap2);
3199         if (is_subset != 1)
3200                 return is_subset;
3201         is_subset = isl_basic_map_is_subset(bmap2, bmap1);
3202         if (is_subset == -1)
3203                 return is_subset;
3204         return !is_subset;
3205 }
3206
3207 static int basic_map_contains(struct isl_basic_map *bmap, struct isl_vec *vec)
3208 {
3209         int i;
3210         unsigned total;
3211         isl_int s;
3212
3213         total = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3214         if (total != vec->size)
3215                 return -1;
3216
3217         isl_int_init(s);
3218
3219         for (i = 0; i < bmap->n_eq; ++i) {
3220                 isl_seq_inner_product(vec->block.data, bmap->eq[i], total, &s);
3221                 if (!isl_int_is_zero(s)) {
3222                         isl_int_clear(s);
3223                         return 0;
3224                 }
3225         }
3226
3227         for (i = 0; i < bmap->n_ineq; ++i) {
3228                 isl_seq_inner_product(vec->block.data, bmap->ineq[i], total, &s);
3229                 if (isl_int_is_neg(s)) {
3230                         isl_int_clear(s);
3231                         return 0;
3232                 }
3233         }
3234
3235         isl_int_clear(s);
3236
3237         return 1;
3238 }
3239
3240 int isl_basic_map_is_universe(struct isl_basic_map *bmap)
3241 {
3242         if (!bmap)
3243                 return -1;
3244         return bmap->n_eq == 0 && bmap->n_ineq == 0;
3245 }
3246
3247 int isl_basic_map_is_empty(struct isl_basic_map *bmap)
3248 {
3249         struct isl_basic_set *bset = NULL;
3250         struct isl_vec *sample = NULL;
3251         int empty;
3252         unsigned total;
3253
3254         if (!bmap)
3255                 return -1;
3256
3257         if (F_ISSET(bmap, ISL_BASIC_MAP_EMPTY))
3258                 return 1;
3259
3260         if (F_ISSET(bmap, ISL_BASIC_MAP_RATIONAL)) {
3261                 struct isl_basic_map *copy = isl_basic_map_copy(bmap);
3262                 copy = isl_basic_map_convex_hull(copy);
3263                 empty = F_ISSET(copy, ISL_BASIC_MAP_EMPTY);
3264                 isl_basic_map_free(copy);
3265                 return empty;
3266         }
3267
3268         total = 1 + bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3269         if (bmap->sample && bmap->sample->size == total) {
3270                 int contains = basic_map_contains(bmap, bmap->sample);
3271                 if (contains < 0)
3272                         return -1;
3273                 if (contains)
3274                         return 0;
3275         }
3276         bset = isl_basic_map_underlying_set(isl_basic_map_copy(bmap));
3277         if (!bset)
3278                 return -1;
3279         sample = isl_basic_set_sample(bset);
3280         if (!sample)
3281                 return -1;
3282         empty = sample->size == 0;
3283         if (bmap->sample)
3284                 isl_vec_free(bmap->ctx, bmap->sample);
3285         bmap->sample = sample;
3286
3287         return empty;
3288 }
3289
3290 struct isl_map *isl_basic_map_union(
3291         struct isl_basic_map *bmap1, struct isl_basic_map *bmap2)
3292 {
3293         struct isl_map *map;
3294         if (!bmap1 || !bmap2)
3295                 return NULL;
3296
3297         isl_assert(bmap1->ctx, bmap1->nparam == bmap2->nparam, goto error);
3298         isl_assert(bmap1->ctx, bmap1->n_in == bmap2->n_in, goto error);
3299         isl_assert(bmap1->ctx, bmap1->n_out == bmap2->n_out, goto error);
3300
3301         map = isl_map_alloc(bmap1->ctx, bmap1->nparam,
3302                                 bmap1->n_in, bmap1->n_out, 2, 0);
3303         if (!map)
3304                 goto error;
3305         map = isl_map_add(map, bmap1);
3306         map = isl_map_add(map, bmap2);
3307         return map;
3308 error:
3309         isl_basic_map_free(bmap1);
3310         isl_basic_map_free(bmap2);
3311         return NULL;
3312 }
3313
3314 struct isl_set *isl_basic_set_union(
3315                 struct isl_basic_set *bset1, struct isl_basic_set *bset2)
3316 {
3317         return (struct isl_set *)isl_basic_map_union(
3318                                             (struct isl_basic_map *)bset1,
3319                                             (struct isl_basic_map *)bset2);
3320 }
3321
3322 /* Order divs such that any div only depends on previous divs */
3323 static struct isl_basic_map *order_divs(struct isl_basic_map *bmap)
3324 {
3325         int i;
3326         unsigned off = bmap->nparam + bmap->n_in + bmap->n_out;
3327
3328         for (i = 0; i < bmap->n_div; ++i) {
3329                 int pos;
3330                 pos = isl_seq_first_non_zero(bmap->div[i]+1+1+off+i,
3331                                                             bmap->n_div-i);
3332                 if (pos == -1)
3333                         continue;
3334                 swap_div(bmap, i, pos);
3335                 --i;
3336         }
3337         return bmap;
3338 }
3339
3340 static int find_div(struct isl_basic_map *dst,
3341                         struct isl_basic_map *src, unsigned div)
3342 {
3343         int i;
3344
3345         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
3346
3347         for (i = 0; i < dst->n_div; ++i)
3348                 if (isl_seq_eq(dst->div[i], src->div[div], 1+1+total) &&
3349                     isl_seq_first_non_zero(dst->div[i]+1+1+total,
3350                                                 dst->n_div - src->n_div) == -1)
3351                         return i;
3352         return -1;
3353 }
3354
3355 struct isl_basic_map *isl_basic_map_align_divs(
3356                 struct isl_basic_map *dst, struct isl_basic_map *src)
3357 {
3358         int i;
3359         unsigned total = src->nparam + src->n_in + src->n_out + src->n_div;
3360
3361         if (src->n_div == 0)
3362                 return dst;
3363
3364         src = order_divs(src);
3365         dst = isl_basic_map_extend(dst, dst->nparam, dst->n_in,
3366                         dst->n_out, src->n_div, 0, 2 * src->n_div);
3367         if (!dst)
3368                 return NULL;
3369         for (i = 0; i < src->n_div; ++i) {
3370                 int j = find_div(dst, src, i);
3371                 if (j < 0) {
3372                         j = isl_basic_map_alloc_div(dst);
3373                         if (j < 0)
3374                                 goto error;
3375                         isl_seq_cpy(dst->div[j], src->div[i], 1+1+total);
3376                         isl_seq_clr(dst->div[j]+1+1+total,
3377                                                     dst->extra - src->n_div);
3378                         if (add_div_constraints(dst, j) < 0)
3379                                 goto error;
3380                 }
3381                 if (j != i)
3382                         swap_div(dst, i, j);
3383         }
3384         return dst;
3385 error:
3386         isl_basic_map_free(dst);
3387         return NULL;
3388 }
3389
3390 struct isl_map *isl_map_align_divs(struct isl_map *map)
3391 {
3392         int i;
3393
3394         map = isl_map_compute_divs(map);
3395         map = isl_map_cow(map);
3396         if (!map)
3397                 return NULL;
3398
3399         for (i = 1; i < map->n; ++i)
3400                 map->p[0] = isl_basic_map_align_divs(map->p[0], map->p[i]);
3401         for (i = 1; i < map->n; ++i)
3402                 map->p[i] = isl_basic_map_align_divs(map->p[i], map->p[0]);
3403
3404         return map;
3405 }
3406
3407 static struct isl_map *add_cut_constraint(struct isl_map *dst,
3408                 struct isl_basic_map *src, isl_int *c,
3409                 unsigned len, int oppose)
3410 {
3411         struct isl_basic_map *copy = NULL;
3412         int is_empty;
3413         int k;
3414         unsigned total;
3415
3416         copy = isl_basic_map_copy(src);
3417         copy = isl_basic_map_cow(copy);
3418         if (!copy)
3419                 goto error;
3420         copy = isl_basic_map_extend(copy,
3421                 copy->nparam, copy->n_in, copy->n_out, 0, 0, 1);
3422         k = isl_basic_map_alloc_inequality(copy);
3423         if (k < 0)
3424                 goto error;
3425         if (oppose)
3426                 isl_seq_neg(copy->ineq[k], c, len);
3427         else
3428                 isl_seq_cpy(copy->ineq[k], c, len);
3429         total = 1 + copy->nparam + copy->n_in + copy->n_out + copy->extra;
3430         isl_seq_clr(copy->ineq[k]+len, total - len);
3431         isl_inequality_negate(copy, k);
3432         copy = isl_basic_map_simplify(copy);
3433         copy = isl_basic_map_finalize(copy);
3434         is_empty = isl_basic_map_is_empty(copy);
3435         if (is_empty < 0)
3436                 goto error;
3437         if (!is_empty)
3438                 dst = isl_map_add(dst, copy);
3439         else
3440                 isl_basic_map_free(copy);
3441         return dst;
3442 error:
3443         isl_basic_map_free(copy);
3444         isl_map_free(dst);
3445         return NULL;
3446 }
3447
3448 static struct isl_map *subtract(struct isl_map *map, struct isl_basic_map *bmap)
3449 {
3450         int i, j, k;
3451         unsigned flags = 0;
3452         struct isl_map *rest = NULL;
3453         unsigned max;
3454         unsigned total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3455
3456         assert(bmap);
3457
3458         if (!map)
3459                 goto error;
3460
3461         if (F_ISSET(map, ISL_MAP_DISJOINT))
3462                 FL_SET(flags, ISL_MAP_DISJOINT);
3463
3464         max = map->n * (2 * bmap->n_eq + bmap->n_ineq);
3465         rest = isl_map_alloc(map->ctx, map->nparam, map->n_in, map->n_out,
3466                                 max, flags);
3467         if (!rest)
3468                 goto error;
3469
3470         for (i = 0; i < map->n; ++i) {
3471                 map->p[i] = isl_basic_map_align_divs(map->p[i], bmap);
3472                 if (!map->p[i])
3473                         goto error;
3474         }
3475
3476         for (j = 0; j < map->n; ++j)
3477                 map->p[j] = isl_basic_map_cow(map->p[j]);
3478
3479         for (i = 0; i < bmap->n_eq; ++i) {
3480                 for (j = 0; j < map->n; ++j) {
3481                         rest = add_cut_constraint(rest,
3482                                 map->p[j], bmap->eq[i], 1+total, 0);
3483                         if (!rest)
3484                                 goto error;
3485
3486                         rest = add_cut_constraint(rest,
3487                                 map->p[j], bmap->eq[i], 1+total, 1);
3488                         if (!rest)
3489                                 goto error;
3490
3491                         map->p[j] = isl_basic_map_extend(map->p[j],
3492                                 map->p[j]->nparam, map->p[j]->n_in,
3493                                 map->p[j]->n_out, 0, 1, 0);
3494                         if (!map->p[j])
3495                                 goto error;
3496                         k = isl_basic_map_alloc_equality(map->p[j]);
3497                         if (k < 0)
3498                                 goto error;
3499                         isl_seq_cpy(map->p[j]->eq[k], bmap->eq[i], 1+total);
3500                         isl_seq_clr(map->p[j]->eq[k]+1+total,
3501                                         map->p[j]->extra - bmap->n_div);
3502                 }
3503         }
3504
3505         for (i = 0; i < bmap->n_ineq; ++i) {
3506                 for (j = 0; j < map->n; ++j) {
3507                         rest = add_cut_constraint(rest,
3508                                 map->p[j], bmap->ineq[i], 1+total, 0);
3509                         if (!rest)
3510                                 goto error;
3511
3512                         map->p[j] = isl_basic_map_extend(map->p[j],
3513                                 map->p[j]->nparam, map->p[j]->n_in,
3514                                 map->p[j]->n_out, 0, 0, 1);
3515                         if (!map->p[j])
3516                                 goto error;
3517                         k = isl_basic_map_alloc_inequality(map->p[j]);
3518                         if (k < 0)
3519                                 goto error;
3520                         isl_seq_cpy(map->p[j]->ineq[k], bmap->ineq[i], 1+total);
3521                         isl_seq_clr(map->p[j]->ineq[k]+1+total,
3522                                         map->p[j]->extra - bmap->n_div);
3523                 }
3524         }
3525
3526         isl_map_free(map);
3527         return rest;
3528 error:
3529         isl_map_free(map);
3530         isl_map_free(rest);
3531         return NULL;
3532 }
3533
3534 struct isl_map *isl_map_subtract(struct isl_map *map1, struct isl_map *map2)
3535 {
3536         int i;
3537         if (!map1 || !map2)
3538                 goto error;
3539
3540         isl_assert(map1->ctx, map1->nparam == map2->nparam, goto error);
3541         isl_assert(map1->ctx, map1->n_in == map2->n_in, goto error);
3542         isl_assert(map1->ctx, map1->n_out == map2->n_out, goto error);
3543
3544         if (isl_map_is_empty(map2)) {
3545                 isl_map_free(map2);
3546                 return map1;
3547         }
3548
3549         map1 = isl_map_compute_divs(map1);
3550         map2 = isl_map_compute_divs(map2);
3551         if (!map1 || !map2)
3552                 goto error;
3553
3554         for (i = 0; map1 && i < map2->n; ++i)
3555                 map1 = subtract(map1, map2->p[i]);
3556
3557         isl_map_free(map2);
3558         return map1;
3559 error:
3560         isl_map_free(map1);
3561         isl_map_free(map2);
3562         return NULL;
3563 }
3564
3565 struct isl_set *isl_set_subtract(struct isl_set *set1, struct isl_set *set2)
3566 {
3567         return (struct isl_set *)
3568                 isl_map_subtract(
3569                         (struct isl_map *)set1, (struct isl_map *)set2);
3570 }
3571
3572 struct isl_set *isl_set_apply(struct isl_set *set, struct isl_map *map)
3573 {
3574         if (!set || !map)
3575                 goto error;
3576         isl_assert(set->ctx, set->dim == map->n_in, goto error);
3577         map = isl_map_intersect_domain(map, set);
3578         set = isl_map_range(map);
3579         return set;
3580 error:
3581         isl_set_free(set);
3582         isl_map_free(map);
3583         return NULL;
3584 }
3585
3586 /* There is no need to cow as removing empty parts doesn't change
3587  * the meaning of the set.
3588  */
3589 struct isl_map *isl_map_remove_empty_parts(struct isl_map *map)
3590 {
3591         int i;
3592
3593         if (!map)
3594                 return NULL;
3595
3596         for (i = map->n-1; i >= 0; --i) {
3597                 if (!F_ISSET(map->p[i], ISL_BASIC_MAP_EMPTY))
3598                         continue;
3599                 isl_basic_map_free(map->p[i]);
3600                 if (i != map->n-1)
3601                         map->p[i] = map->p[map->n-1];
3602                 map->n--;
3603         }
3604
3605         return map;
3606 }
3607
3608 struct isl_set *isl_set_remove_empty_parts(struct isl_set *set)
3609 {
3610         return (struct isl_set *)
3611                 isl_map_remove_empty_parts((struct isl_map *)set);
3612 }
3613
3614 struct isl_basic_set *isl_set_copy_basic_set(struct isl_set *set)
3615 {
3616         struct isl_basic_set *bset;
3617         if (!set || set->n == 0)
3618                 return NULL;
3619         bset = set->p[set->n-1];
3620         isl_assert(set->ctx, F_ISSET(bset, ISL_BASIC_SET_FINAL), return NULL);
3621         return isl_basic_set_copy(bset);
3622 }
3623
3624 struct isl_set *isl_set_drop_basic_set(struct isl_set *set,
3625                                                 struct isl_basic_set *bset)
3626 {
3627         int i;
3628
3629         if (!set || !bset)
3630                 goto error;
3631         for (i = set->n-1; i >= 0; --i) {
3632                 if (set->p[i] != bset)
3633                         continue;
3634                 set = isl_set_cow(set);
3635                 if (!set)
3636                         goto error;
3637                 isl_basic_set_free(set->p[i]);
3638                 if (i != set->n-1)
3639                         set->p[i] = set->p[set->n-1];
3640                 set->n--;
3641                 return set;
3642         }
3643         isl_basic_set_free(bset);
3644         return set;
3645 error:
3646         isl_set_free(set);
3647         isl_basic_set_free(bset);
3648         return NULL;
3649 }
3650
3651 /* Given two _disjoint_ basic sets bset1 and bset2, check whether
3652  * for any common value of the parameters and dimensions preceding dim
3653  * in both basic sets, the values of dimension pos in bset1 are
3654  * smaller or larger then those in bset2.
3655  *
3656  * Returns
3657  *       1 if bset1 follows bset2
3658  *      -1 if bset1 precedes bset2
3659  *       0 if bset1 and bset2 are incomparable
3660  *      -2 if some error occurred.
3661  */
3662 int isl_basic_set_compare_at(struct isl_basic_set *bset1,
3663         struct isl_basic_set *bset2, int pos)
3664 {
3665         struct isl_basic_map *bmap1 = NULL;
3666         struct isl_basic_map *bmap2 = NULL;
3667         struct isl_ctx *ctx;
3668         struct isl_vec *obj;
3669         unsigned total;
3670         isl_int num, den;
3671         enum isl_lp_result res;
3672         int cmp;
3673
3674         if (!bset1 || !bset2)
3675                 return -2;
3676
3677         bmap1 = isl_basic_map_from_basic_set(isl_basic_set_copy(bset1),
3678                                                 pos, bset1->dim - pos);
3679         bmap2 = isl_basic_map_from_basic_set(isl_basic_set_copy(bset2),
3680                                                 pos, bset2->dim - pos);
3681         if (!bmap1 || !bmap2)
3682                 goto error;
3683         bmap1 = isl_basic_map_extend(bmap1, bmap1->nparam,
3684                         bmap1->n_in, bmap1->n_out + bmap2->n_out,
3685                         bmap2->n_div, bmap2->n_eq, bmap2->n_ineq);
3686         if (!bmap1)
3687                 goto error;
3688         total = bmap1->nparam + bmap1->n_in + bmap1->n_out + bmap1->n_div;
3689         ctx = bmap1->ctx;
3690         obj = isl_vec_alloc(ctx, total);
3691         isl_seq_clr(obj->block.data, total);
3692         isl_int_set_si(obj->block.data[bmap1->nparam+bmap1->n_in], 1);
3693         isl_int_set_si(obj->block.data[bmap1->nparam+bmap1->n_in+
3694                                         bmap1->n_out-bmap2->n_out], -1);
3695         if (!obj)
3696                 goto error;
3697         bmap1 = add_constraints(bmap1, bmap2, 0, bmap1->n_out - bmap2->n_out);
3698         isl_int_init(num);
3699         isl_int_init(den);
3700         res = isl_solve_lp(bmap1, 0, obj->block.data, ctx->one, &num, &den);
3701         if (res == isl_lp_empty)
3702                 cmp = 0;
3703         else if (res == isl_lp_ok && isl_int_is_pos(num))
3704                 cmp = 1;
3705         else if ((res == isl_lp_ok && isl_int_is_neg(num)) ||
3706                   res == isl_lp_unbounded)
3707                 cmp = -1;
3708         else
3709                 cmp = -2;
3710         isl_int_clear(num);
3711         isl_int_clear(den);
3712         isl_basic_map_free(bmap1);
3713         isl_vec_free(ctx, obj);
3714         return cmp;
3715 error:
3716         isl_basic_map_free(bmap1);
3717         isl_basic_map_free(bmap2);
3718         return -2;
3719 }
3720
3721 static int isl_basic_map_fast_has_fixed_var(struct isl_basic_map *bmap,
3722         unsigned pos, isl_int *val)
3723 {
3724         int i;
3725         int d;
3726         unsigned total;
3727
3728         if (!bmap)
3729                 return -1;
3730         total = bmap->nparam + bmap->n_in + bmap->n_out + bmap->n_div;
3731         for (i = 0, d = total-1; i < bmap->n_eq && d+1 > pos; ++i) {
3732                 for (; d+1 > pos; --d)
3733                         if (!isl_int_is_zero(bmap->eq[i][1+d]))
3734                                 break;
3735                 if (d != pos)
3736                         continue;
3737                 if (isl_seq_first_non_zero(bmap->eq[i]+1, d) != -1)
3738                         return 0;
3739                 if (isl_seq_first_non_zero(bmap->eq[i]+1+d+1, total-d-1) != -1)
3740                         return 0;
3741                 if (!isl_int_is_one(bmap->eq[i][1+d]))
3742                         return 0;
3743                 if (val)
3744                         isl_int_neg(*val, bmap->eq[i][0]);
3745                 return 1;
3746         }
3747         return 0;
3748 }
3749
3750 static int isl_map_fast_has_fixed_var(struct isl_map *map,
3751         unsigned pos, isl_int *val)
3752 {
3753         int i;
3754         isl_int v;
3755         isl_int tmp;
3756         int fixed;
3757
3758         if (!map)
3759                 return -1;
3760         if (map->n == 0)
3761                 return 0;
3762         if (map->n == 1)
3763                 return isl_basic_map_fast_has_fixed_var(map->p[0], pos, val); 
3764         isl_int_init(v);
3765         isl_int_init(tmp);
3766         fixed = isl_basic_map_fast_has_fixed_var(map->p[0], pos, &v); 
3767         for (i = 1; fixed == 1 && i < map->n; ++i) {
3768                 fixed = isl_basic_map_fast_has_fixed_var(map->p[i], pos, &tmp); 
3769                 if (fixed == 1 && isl_int_ne(tmp, v))
3770                         fixed = 0;
3771         }
3772         if (val)
3773                 isl_int_set(*val, v);
3774         isl_int_clear(tmp);
3775         isl_int_clear(v);
3776         return fixed;
3777 }
3778
3779 static int isl_set_fast_has_fixed_var(struct isl_set *set, unsigned pos,
3780         isl_int *val)
3781 {
3782         return isl_map_fast_has_fixed_var((struct isl_map *)set, pos, val);
3783 }
3784
3785 /* Check if dimension dim has fixed value and if so and if val is not NULL,
3786  * then return this fixed value in *val.
3787  */
3788 int isl_set_fast_dim_is_fixed(struct isl_set *set, unsigned dim, isl_int *val)
3789 {
3790         return isl_set_fast_has_fixed_var(set, set->nparam + dim, val);
3791 }
3792
3793 /* Check if input variable in has fixed value and if so and if val is not NULL,
3794  * then return this fixed value in *val.
3795  */
3796 int isl_map_fast_input_is_fixed(struct isl_map *map, unsigned in, isl_int *val)
3797 {
3798         return isl_map_fast_has_fixed_var(map, map->nparam + in, val);
3799 }
3800
3801 /* Remove all information from bset that is redundant in the context
3802  * of context.  In particular, equalities that are linear combinations
3803  * of those in context are remobed.  Then the inequalities that are
3804  * redundant in the context of the equalities and inequalities of
3805  * context are removed.
3806  */
3807 static struct isl_basic_set *uset_gist(struct isl_basic_set *bset,
3808         struct isl_basic_set *context)
3809 {
3810         int i;
3811         isl_int opt;
3812         struct isl_basic_set *combined;
3813         struct isl_ctx *ctx;
3814         enum isl_lp_result res = isl_lp_ok;
3815
3816         if (!bset || !context)
3817                 goto error;
3818
3819         if (context->n_eq > 0) {
3820                 struct isl_mat *T;
3821                 struct isl_mat *T2;
3822                 struct isl_ctx *ctx = context->ctx;
3823                 context = isl_basic_set_remove_equalities(context, &T, &T2);
3824                 if (!context)
3825                         goto error;
3826                 bset = isl_basic_set_preimage(ctx, bset, T);
3827                 bset = uset_gist(bset, context);
3828                 bset = isl_basic_set_preimage(ctx, bset, T2);
3829                 return bset;
3830         }
3831         combined = isl_basic_set_extend(isl_basic_set_copy(bset),
3832                         0, bset->dim, 0, context->n_eq, context->n_ineq);
3833         context = set_add_constraints(combined, context, 0);
3834         if (!context)
3835                 goto error;
3836         ctx = context->ctx;
3837         isl_int_init(opt);
3838         for (i = bset->n_ineq-1; i >= 0; --i) {
3839                 set_swap_inequality(context, i, context->n_ineq-1);
3840                 context->n_ineq--;
3841                 res = isl_solve_lp((struct isl_basic_map *)context, 0,
3842                         context->ineq[context->n_ineq]+1, ctx->one, &opt, NULL);
3843                 context->n_ineq++;
3844                 set_swap_inequality(context, i, context->n_ineq-1);
3845                 if (res == isl_lp_unbounded)
3846                         continue;
3847                 if (res == isl_lp_error)
3848                         break;
3849                 if (res == isl_lp_empty) {
3850                         bset = isl_basic_set_set_to_empty(bset);
3851                         break;
3852                 }
3853                 isl_int_add(opt, opt, context->ineq[i][0]);
3854                 if (!isl_int_is_neg(opt)) {
3855                         isl_basic_set_drop_inequality(context, i);
3856                         isl_basic_set_drop_inequality(bset, i);
3857                 }
3858         }
3859         isl_int_clear(opt);
3860         if (res == isl_lp_error)
3861                 goto error;
3862         isl_basic_set_free(context);
3863         return bset;
3864 error:
3865         isl_basic_set_free(context);
3866         isl_basic_set_free(bset);
3867         return NULL;
3868 }
3869
3870 struct isl_basic_map *isl_basic_map_gist(struct isl_basic_map *bmap,
3871         struct isl_basic_map *context)
3872 {
3873         struct isl_basic_set *bset;
3874
3875         if (!bmap || !context)
3876                 goto error;
3877
3878         context = isl_basic_map_align_divs(context, bmap);
3879         bmap = isl_basic_map_align_divs(bmap, context);
3880
3881         bset = uset_gist(isl_basic_map_underlying_set(isl_basic_map_copy(bmap)),
3882                          isl_basic_map_underlying_set(context));
3883
3884         return isl_basic_map_overlying_set(bset, bmap);
3885 error:
3886         isl_basic_map_free(bmap);
3887         isl_basic_map_free(context);
3888         return NULL;
3889 }
3890
3891 struct isl_map *isl_map_gist(struct isl_map *map, struct isl_basic_map *context)
3892 {
3893         int i;
3894
3895         map = isl_map_cow(map);
3896         if (!map || !context)
3897                 return NULL;
3898         isl_assert(map->ctx, map->nparam == context->nparam, goto error);
3899         isl_assert(map->ctx, map->n_in == context->n_in, goto error);
3900         isl_assert(map->ctx, map->n_out == context->n_out, goto error);
3901         for (i = 0; i < map->n; ++i)
3902                 context = isl_basic_map_align_divs(context, map->p[i]);
3903         for (i = 0; i < map->n; ++i) {
3904                 map->p[i] = isl_basic_map_gist(map->p[i],
3905                                                 isl_basic_map_copy(context));
3906                 if (!map->p[i])
3907                         goto error;
3908         }
3909         isl_basic_map_free(context);
3910         return map;
3911 error:
3912         isl_map_free(map);
3913         isl_basic_map_free(context);
3914         return NULL;
3915 }
3916
3917 struct isl_set *isl_set_gist(struct isl_set *set, struct isl_basic_set *context)
3918 {
3919         return (struct isl_set *)isl_map_gist((struct isl_map *)set,
3920                                         (struct isl_basic_map *)context);
3921 }