add isl_set_indicator_function
[platform/upstream/isl.git] / isl_aff.c
1 /*
2  * Copyright 2011      INRIA Saclay
3  * Copyright 2011      Sven Verdoolaege
4  * Copyright 2012      Ecole Normale Superieure
5  *
6  * Use of this software is governed by the GNU LGPLv2.1 license
7  *
8  * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
9  * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
10  * 91893 Orsay, France
11  * and Ecole Normale Superieure, 45 rue d’Ulm, 75230 Paris, France
12  */
13
14 #include <isl_ctx_private.h>
15 #include <isl_map_private.h>
16 #include <isl_aff_private.h>
17 #include <isl_space_private.h>
18 #include <isl_local_space_private.h>
19 #include <isl_mat_private.h>
20 #include <isl_list_private.h>
21 #include <isl/constraint.h>
22 #include <isl/seq.h>
23 #include <isl/set.h>
24 #include <isl_config.h>
25
26 __isl_give isl_aff *isl_aff_alloc_vec(__isl_take isl_local_space *ls,
27         __isl_take isl_vec *v)
28 {
29         isl_aff *aff;
30
31         if (!ls || !v)
32                 goto error;
33
34         aff = isl_calloc_type(v->ctx, struct isl_aff);
35         if (!aff)
36                 goto error;
37
38         aff->ref = 1;
39         aff->ls = ls;
40         aff->v = v;
41
42         return aff;
43 error:
44         isl_local_space_free(ls);
45         isl_vec_free(v);
46         return NULL;
47 }
48
49 __isl_give isl_aff *isl_aff_alloc(__isl_take isl_local_space *ls)
50 {
51         isl_ctx *ctx;
52         isl_vec *v;
53         unsigned total;
54
55         if (!ls)
56                 return NULL;
57
58         ctx = isl_local_space_get_ctx(ls);
59         if (!isl_local_space_divs_known(ls))
60                 isl_die(ctx, isl_error_invalid, "local space has unknown divs",
61                         goto error);
62         if (!isl_local_space_is_set(ls))
63                 isl_die(ctx, isl_error_invalid,
64                         "domain of affine expression should be a set",
65                         goto error);
66
67         total = isl_local_space_dim(ls, isl_dim_all);
68         v = isl_vec_alloc(ctx, 1 + 1 + total);
69         return isl_aff_alloc_vec(ls, v);
70 error:
71         isl_local_space_free(ls);
72         return NULL;
73 }
74
75 __isl_give isl_aff *isl_aff_zero_on_domain(__isl_take isl_local_space *ls)
76 {
77         isl_aff *aff;
78
79         aff = isl_aff_alloc(ls);
80         if (!aff)
81                 return NULL;
82
83         isl_int_set_si(aff->v->el[0], 1);
84         isl_seq_clr(aff->v->el + 1, aff->v->size - 1);
85
86         return aff;
87 }
88
89 __isl_give isl_aff *isl_aff_copy(__isl_keep isl_aff *aff)
90 {
91         if (!aff)
92                 return NULL;
93
94         aff->ref++;
95         return aff;
96 }
97
98 __isl_give isl_aff *isl_aff_dup(__isl_keep isl_aff *aff)
99 {
100         if (!aff)
101                 return NULL;
102
103         return isl_aff_alloc_vec(isl_local_space_copy(aff->ls),
104                                  isl_vec_copy(aff->v));
105 }
106
107 __isl_give isl_aff *isl_aff_cow(__isl_take isl_aff *aff)
108 {
109         if (!aff)
110                 return NULL;
111
112         if (aff->ref == 1)
113                 return aff;
114         aff->ref--;
115         return isl_aff_dup(aff);
116 }
117
118 void *isl_aff_free(__isl_take isl_aff *aff)
119 {
120         if (!aff)
121                 return NULL;
122
123         if (--aff->ref > 0)
124                 return NULL;
125
126         isl_local_space_free(aff->ls);
127         isl_vec_free(aff->v);
128
129         free(aff);
130
131         return NULL;
132 }
133
134 isl_ctx *isl_aff_get_ctx(__isl_keep isl_aff *aff)
135 {
136         return aff ? isl_local_space_get_ctx(aff->ls) : NULL;
137 }
138
139 /* Externally, an isl_aff has a map space, but internally, the
140  * ls field corresponds to the domain of that space.
141  */
142 int isl_aff_dim(__isl_keep isl_aff *aff, enum isl_dim_type type)
143 {
144         if (!aff)
145                 return 0;
146         if (type == isl_dim_out)
147                 return 1;
148         if (type == isl_dim_in)
149                 type = isl_dim_set;
150         return isl_local_space_dim(aff->ls, type);
151 }
152
153 __isl_give isl_space *isl_aff_get_domain_space(__isl_keep isl_aff *aff)
154 {
155         return aff ? isl_local_space_get_space(aff->ls) : NULL;
156 }
157
158 __isl_give isl_space *isl_aff_get_space(__isl_keep isl_aff *aff)
159 {
160         isl_space *space;
161         if (!aff)
162                 return NULL;
163         space = isl_local_space_get_space(aff->ls);
164         space = isl_space_from_domain(space);
165         space = isl_space_add_dims(space, isl_dim_out, 1);
166         return space;
167 }
168
169 __isl_give isl_local_space *isl_aff_get_domain_local_space(
170         __isl_keep isl_aff *aff)
171 {
172         return aff ? isl_local_space_copy(aff->ls) : NULL;
173 }
174
175 __isl_give isl_local_space *isl_aff_get_local_space(__isl_keep isl_aff *aff)
176 {
177         isl_local_space *ls;
178         if (!aff)
179                 return NULL;
180         ls = isl_local_space_copy(aff->ls);
181         ls = isl_local_space_from_domain(ls);
182         ls = isl_local_space_add_dims(ls, isl_dim_out, 1);
183         return ls;
184 }
185
186 /* Externally, an isl_aff has a map space, but internally, the
187  * ls field corresponds to the domain of that space.
188  */
189 const char *isl_aff_get_dim_name(__isl_keep isl_aff *aff,
190         enum isl_dim_type type, unsigned pos)
191 {
192         if (!aff)
193                 return NULL;
194         if (type == isl_dim_out)
195                 return NULL;
196         if (type == isl_dim_in)
197                 type = isl_dim_set;
198         return isl_local_space_get_dim_name(aff->ls, type, pos);
199 }
200
201 __isl_give isl_aff *isl_aff_reset_domain_space(__isl_take isl_aff *aff,
202         __isl_take isl_space *dim)
203 {
204         aff = isl_aff_cow(aff);
205         if (!aff || !dim)
206                 goto error;
207
208         aff->ls = isl_local_space_reset_space(aff->ls, dim);
209         if (!aff->ls)
210                 return isl_aff_free(aff);
211
212         return aff;
213 error:
214         isl_aff_free(aff);
215         isl_space_free(dim);
216         return NULL;
217 }
218
219 /* Reset the space of "aff".  This function is called from isl_pw_templ.c
220  * and doesn't know if the space of an element object is represented
221  * directly or through its domain.  It therefore passes along both.
222  */
223 __isl_give isl_aff *isl_aff_reset_space_and_domain(__isl_take isl_aff *aff,
224         __isl_take isl_space *space, __isl_take isl_space *domain)
225 {
226         isl_space_free(space);
227         return isl_aff_reset_domain_space(aff, domain);
228 }
229
230 /* Reorder the coefficients of the affine expression based
231  * on the given reodering.
232  * The reordering r is assumed to have been extended with the local
233  * variables.
234  */
235 static __isl_give isl_vec *vec_reorder(__isl_take isl_vec *vec,
236         __isl_take isl_reordering *r, int n_div)
237 {
238         isl_vec *res;
239         int i;
240
241         if (!vec || !r)
242                 goto error;
243
244         res = isl_vec_alloc(vec->ctx,
245                             2 + isl_space_dim(r->dim, isl_dim_all) + n_div);
246         isl_seq_cpy(res->el, vec->el, 2);
247         isl_seq_clr(res->el + 2, res->size - 2);
248         for (i = 0; i < r->len; ++i)
249                 isl_int_set(res->el[2 + r->pos[i]], vec->el[2 + i]);
250
251         isl_reordering_free(r);
252         isl_vec_free(vec);
253         return res;
254 error:
255         isl_vec_free(vec);
256         isl_reordering_free(r);
257         return NULL;
258 }
259
260 /* Reorder the dimensions of the domain of "aff" according
261  * to the given reordering.
262  */
263 __isl_give isl_aff *isl_aff_realign_domain(__isl_take isl_aff *aff,
264         __isl_take isl_reordering *r)
265 {
266         aff = isl_aff_cow(aff);
267         if (!aff)
268                 goto error;
269
270         r = isl_reordering_extend(r, aff->ls->div->n_row);
271         aff->v = vec_reorder(aff->v, isl_reordering_copy(r),
272                                 aff->ls->div->n_row);
273         aff->ls = isl_local_space_realign(aff->ls, r);
274
275         if (!aff->v || !aff->ls)
276                 return isl_aff_free(aff);
277
278         return aff;
279 error:
280         isl_aff_free(aff);
281         isl_reordering_free(r);
282         return NULL;
283 }
284
285 __isl_give isl_aff *isl_aff_align_params(__isl_take isl_aff *aff,
286         __isl_take isl_space *model)
287 {
288         if (!aff || !model)
289                 goto error;
290
291         if (!isl_space_match(aff->ls->dim, isl_dim_param,
292                              model, isl_dim_param)) {
293                 isl_reordering *exp;
294
295                 model = isl_space_drop_dims(model, isl_dim_in,
296                                         0, isl_space_dim(model, isl_dim_in));
297                 model = isl_space_drop_dims(model, isl_dim_out,
298                                         0, isl_space_dim(model, isl_dim_out));
299                 exp = isl_parameter_alignment_reordering(aff->ls->dim, model);
300                 exp = isl_reordering_extend_space(exp,
301                                         isl_aff_get_domain_space(aff));
302                 aff = isl_aff_realign_domain(aff, exp);
303         }
304
305         isl_space_free(model);
306         return aff;
307 error:
308         isl_space_free(model);
309         isl_aff_free(aff);
310         return NULL;
311 }
312
313 int isl_aff_plain_is_zero(__isl_keep isl_aff *aff)
314 {
315         if (!aff)
316                 return -1;
317
318         return isl_seq_first_non_zero(aff->v->el + 1, aff->v->size - 1) < 0;
319 }
320
321 int isl_aff_plain_is_equal(__isl_keep isl_aff *aff1, __isl_keep isl_aff *aff2)
322 {
323         int equal;
324
325         if (!aff1 || !aff2)
326                 return -1;
327
328         equal = isl_local_space_is_equal(aff1->ls, aff2->ls);
329         if (equal < 0 || !equal)
330                 return equal;
331
332         return isl_vec_is_equal(aff1->v, aff2->v);
333 }
334
335 int isl_aff_get_denominator(__isl_keep isl_aff *aff, isl_int *v)
336 {
337         if (!aff)
338                 return -1;
339         isl_int_set(*v, aff->v->el[0]);
340         return 0;
341 }
342
343 int isl_aff_get_constant(__isl_keep isl_aff *aff, isl_int *v)
344 {
345         if (!aff)
346                 return -1;
347         isl_int_set(*v, aff->v->el[1]);
348         return 0;
349 }
350
351 int isl_aff_get_coefficient(__isl_keep isl_aff *aff,
352         enum isl_dim_type type, int pos, isl_int *v)
353 {
354         if (!aff)
355                 return -1;
356
357         if (type == isl_dim_out)
358                 isl_die(aff->v->ctx, isl_error_invalid,
359                         "output/set dimension does not have a coefficient",
360                         return -1);
361         if (type == isl_dim_in)
362                 type = isl_dim_set;
363
364         if (pos >= isl_local_space_dim(aff->ls, type))
365                 isl_die(aff->v->ctx, isl_error_invalid,
366                         "position out of bounds", return -1);
367
368         pos += isl_local_space_offset(aff->ls, type);
369         isl_int_set(*v, aff->v->el[1 + pos]);
370
371         return 0;
372 }
373
374 __isl_give isl_aff *isl_aff_set_denominator(__isl_take isl_aff *aff, isl_int v)
375 {
376         aff = isl_aff_cow(aff);
377         if (!aff)
378                 return NULL;
379
380         aff->v = isl_vec_cow(aff->v);
381         if (!aff->v)
382                 return isl_aff_free(aff);
383
384         isl_int_set(aff->v->el[0], v);
385
386         return aff;
387 }
388
389 __isl_give isl_aff *isl_aff_set_constant(__isl_take isl_aff *aff, isl_int v)
390 {
391         aff = isl_aff_cow(aff);
392         if (!aff)
393                 return NULL;
394
395         aff->v = isl_vec_cow(aff->v);
396         if (!aff->v)
397                 return isl_aff_free(aff);
398
399         isl_int_set(aff->v->el[1], v);
400
401         return aff;
402 }
403
404 __isl_give isl_aff *isl_aff_add_constant(__isl_take isl_aff *aff, isl_int v)
405 {
406         if (isl_int_is_zero(v))
407                 return aff;
408
409         aff = isl_aff_cow(aff);
410         if (!aff)
411                 return NULL;
412
413         aff->v = isl_vec_cow(aff->v);
414         if (!aff->v)
415                 return isl_aff_free(aff);
416
417         isl_int_addmul(aff->v->el[1], aff->v->el[0], v);
418
419         return aff;
420 }
421
422 __isl_give isl_aff *isl_aff_add_constant_si(__isl_take isl_aff *aff, int v)
423 {
424         isl_int t;
425
426         isl_int_init(t);
427         isl_int_set_si(t, v);
428         aff = isl_aff_add_constant(aff, t);
429         isl_int_clear(t);
430
431         return aff;
432 }
433
434 __isl_give isl_aff *isl_aff_set_constant_si(__isl_take isl_aff *aff, int v)
435 {
436         aff = isl_aff_cow(aff);
437         if (!aff)
438                 return NULL;
439
440         aff->v = isl_vec_cow(aff->v);
441         if (!aff->v)
442                 return isl_aff_free(aff);
443
444         isl_int_set_si(aff->v->el[1], v);
445
446         return aff;
447 }
448
449 __isl_give isl_aff *isl_aff_set_coefficient(__isl_take isl_aff *aff,
450         enum isl_dim_type type, int pos, isl_int v)
451 {
452         if (!aff)
453                 return NULL;
454
455         if (type == isl_dim_out)
456                 isl_die(aff->v->ctx, isl_error_invalid,
457                         "output/set dimension does not have a coefficient",
458                         return isl_aff_free(aff));
459         if (type == isl_dim_in)
460                 type = isl_dim_set;
461
462         if (pos >= isl_local_space_dim(aff->ls, type))
463                 isl_die(aff->v->ctx, isl_error_invalid,
464                         "position out of bounds", return isl_aff_free(aff));
465
466         aff = isl_aff_cow(aff);
467         if (!aff)
468                 return NULL;
469
470         aff->v = isl_vec_cow(aff->v);
471         if (!aff->v)
472                 return isl_aff_free(aff);
473
474         pos += isl_local_space_offset(aff->ls, type);
475         isl_int_set(aff->v->el[1 + pos], v);
476
477         return aff;
478 }
479
480 __isl_give isl_aff *isl_aff_set_coefficient_si(__isl_take isl_aff *aff,
481         enum isl_dim_type type, int pos, int v)
482 {
483         if (!aff)
484                 return NULL;
485
486         if (type == isl_dim_out)
487                 isl_die(aff->v->ctx, isl_error_invalid,
488                         "output/set dimension does not have a coefficient",
489                         return isl_aff_free(aff));
490         if (type == isl_dim_in)
491                 type = isl_dim_set;
492
493         if (pos >= isl_local_space_dim(aff->ls, type))
494                 isl_die(aff->v->ctx, isl_error_invalid,
495                         "position out of bounds", return isl_aff_free(aff));
496
497         aff = isl_aff_cow(aff);
498         if (!aff)
499                 return NULL;
500
501         aff->v = isl_vec_cow(aff->v);
502         if (!aff->v)
503                 return isl_aff_free(aff);
504
505         pos += isl_local_space_offset(aff->ls, type);
506         isl_int_set_si(aff->v->el[1 + pos], v);
507
508         return aff;
509 }
510
511 __isl_give isl_aff *isl_aff_add_coefficient(__isl_take isl_aff *aff,
512         enum isl_dim_type type, int pos, isl_int v)
513 {
514         if (!aff)
515                 return NULL;
516
517         if (type == isl_dim_out)
518                 isl_die(aff->v->ctx, isl_error_invalid,
519                         "output/set dimension does not have a coefficient",
520                         return isl_aff_free(aff));
521         if (type == isl_dim_in)
522                 type = isl_dim_set;
523
524         if (pos >= isl_local_space_dim(aff->ls, type))
525                 isl_die(aff->v->ctx, isl_error_invalid,
526                         "position out of bounds", return isl_aff_free(aff));
527
528         aff = isl_aff_cow(aff);
529         if (!aff)
530                 return NULL;
531
532         aff->v = isl_vec_cow(aff->v);
533         if (!aff->v)
534                 return isl_aff_free(aff);
535
536         pos += isl_local_space_offset(aff->ls, type);
537         isl_int_addmul(aff->v->el[1 + pos], aff->v->el[0], v);
538
539         return aff;
540 }
541
542 __isl_give isl_aff *isl_aff_add_coefficient_si(__isl_take isl_aff *aff,
543         enum isl_dim_type type, int pos, int v)
544 {
545         isl_int t;
546
547         isl_int_init(t);
548         isl_int_set_si(t, v);
549         aff = isl_aff_add_coefficient(aff, type, pos, t);
550         isl_int_clear(t);
551
552         return aff;
553 }
554
555 __isl_give isl_aff *isl_aff_get_div(__isl_keep isl_aff *aff, int pos)
556 {
557         if (!aff)
558                 return NULL;
559
560         return isl_local_space_get_div(aff->ls, pos);
561 }
562
563 __isl_give isl_aff *isl_aff_neg(__isl_take isl_aff *aff)
564 {
565         aff = isl_aff_cow(aff);
566         if (!aff)
567                 return NULL;
568         aff->v = isl_vec_cow(aff->v);
569         if (!aff->v)
570                 return isl_aff_free(aff);
571
572         isl_seq_neg(aff->v->el + 1, aff->v->el + 1, aff->v->size - 1);
573
574         return aff;
575 }
576
577 /* Remove divs from the local space that do not appear in the affine
578  * expression.
579  * We currently only remove divs at the end.
580  * Some intermediate divs may also not appear directly in the affine
581  * expression, but we would also need to check that no other divs are
582  * defined in terms of them.
583  */
584 __isl_give isl_aff *isl_aff_remove_unused_divs( __isl_take isl_aff *aff)
585 {
586         int pos;
587         int off;
588         int n;
589
590         if (!aff)
591                 return NULL;
592
593         n = isl_local_space_dim(aff->ls, isl_dim_div);
594         off = isl_local_space_offset(aff->ls, isl_dim_div);
595
596         pos = isl_seq_last_non_zero(aff->v->el + 1 + off, n) + 1;
597         if (pos == n)
598                 return aff;
599
600         aff = isl_aff_cow(aff);
601         if (!aff)
602                 return NULL;
603
604         aff->ls = isl_local_space_drop_dims(aff->ls, isl_dim_div, pos, n - pos);
605         aff->v = isl_vec_drop_els(aff->v, 1 + off + pos, n - pos);
606         if (!aff->ls || !aff->v)
607                 return isl_aff_free(aff);
608
609         return aff;
610 }
611
612 __isl_give isl_aff *isl_aff_normalize(__isl_take isl_aff *aff)
613 {
614         if (!aff)
615                 return NULL;
616         aff->v = isl_vec_normalize(aff->v);
617         if (!aff->v)
618                 return isl_aff_free(aff);
619         aff = isl_aff_remove_unused_divs(aff);
620         return aff;
621 }
622
623 /* Given f, return floor(f).
624  * If f is an integer expression, then just return f.
625  * Otherwise, if f = g/m, write g = q m + r,
626  * create a new div d = [r/m] and return the expression q + d.
627  * The coefficients in r are taken to lie between -m/2 and m/2.
628  */
629 __isl_give isl_aff *isl_aff_floor(__isl_take isl_aff *aff)
630 {
631         int i;
632         int size;
633         isl_ctx *ctx;
634         isl_vec *div;
635
636         if (!aff)
637                 return NULL;
638
639         if (isl_int_is_one(aff->v->el[0]))
640                 return aff;
641
642         aff = isl_aff_cow(aff);
643         if (!aff)
644                 return NULL;
645
646         aff->v = isl_vec_cow(aff->v);
647         div = isl_vec_copy(aff->v);
648         div = isl_vec_cow(div);
649         if (!div)
650                 return isl_aff_free(aff);
651
652         ctx = isl_aff_get_ctx(aff);
653         isl_int_fdiv_q(aff->v->el[0], aff->v->el[0], ctx->two);
654         for (i = 1; i < aff->v->size; ++i) {
655                 isl_int_fdiv_r(div->el[i], div->el[i], div->el[0]);
656                 isl_int_fdiv_q(aff->v->el[i], aff->v->el[i], div->el[0]);
657                 if (isl_int_gt(div->el[i], aff->v->el[0])) {
658                         isl_int_sub(div->el[i], div->el[i], div->el[0]);
659                         isl_int_add_ui(aff->v->el[i], aff->v->el[i], 1);
660                 }
661         }
662
663         aff->ls = isl_local_space_add_div(aff->ls, div);
664         if (!aff->ls)
665                 return isl_aff_free(aff);
666
667         size = aff->v->size;
668         aff->v = isl_vec_extend(aff->v, size + 1);
669         if (!aff->v)
670                 return isl_aff_free(aff);
671         isl_int_set_si(aff->v->el[0], 1);
672         isl_int_set_si(aff->v->el[size], 1);
673
674         return aff;
675 }
676
677 /* Compute
678  *
679  *      aff mod m = aff - m * floor(aff/m)
680  */
681 __isl_give isl_aff *isl_aff_mod(__isl_take isl_aff *aff, isl_int m)
682 {
683         isl_aff *res;
684
685         res = isl_aff_copy(aff);
686         aff = isl_aff_scale_down(aff, m);
687         aff = isl_aff_floor(aff);
688         aff = isl_aff_scale(aff, m);
689         res = isl_aff_sub(res, aff);
690
691         return res;
692 }
693
694 /* Compute
695  *
696  *      pwaff mod m = pwaff - m * floor(pwaff/m)
697  */
698 __isl_give isl_pw_aff *isl_pw_aff_mod(__isl_take isl_pw_aff *pwaff, isl_int m)
699 {
700         isl_pw_aff *res;
701
702         res = isl_pw_aff_copy(pwaff);
703         pwaff = isl_pw_aff_scale_down(pwaff, m);
704         pwaff = isl_pw_aff_floor(pwaff);
705         pwaff = isl_pw_aff_scale(pwaff, m);
706         res = isl_pw_aff_sub(res, pwaff);
707
708         return res;
709 }
710
711 /* Given f, return ceil(f).
712  * If f is an integer expression, then just return f.
713  * Otherwise, create a new div d = [-f] and return the expression -d.
714  */
715 __isl_give isl_aff *isl_aff_ceil(__isl_take isl_aff *aff)
716 {
717         if (!aff)
718                 return NULL;
719
720         if (isl_int_is_one(aff->v->el[0]))
721                 return aff;
722
723         aff = isl_aff_neg(aff);
724         aff = isl_aff_floor(aff);
725         aff = isl_aff_neg(aff);
726
727         return aff;
728 }
729
730 /* Apply the expansion computed by isl_merge_divs.
731  * The expansion itself is given by "exp" while the resulting
732  * list of divs is given by "div".
733  */
734 __isl_give isl_aff *isl_aff_expand_divs( __isl_take isl_aff *aff,
735         __isl_take isl_mat *div, int *exp)
736 {
737         int i, j;
738         int old_n_div;
739         int new_n_div;
740         int offset;
741
742         aff = isl_aff_cow(aff);
743         if (!aff || !div)
744                 goto error;
745
746         old_n_div = isl_local_space_dim(aff->ls, isl_dim_div);
747         new_n_div = isl_mat_rows(div);
748         if (new_n_div < old_n_div)
749                 isl_die(isl_mat_get_ctx(div), isl_error_invalid,
750                         "not an expansion", goto error);
751
752         aff->v = isl_vec_extend(aff->v, aff->v->size + new_n_div - old_n_div);
753         if (!aff->v)
754                 goto error;
755
756         offset = 1 + isl_local_space_offset(aff->ls, isl_dim_div);
757         j = old_n_div - 1;
758         for (i = new_n_div - 1; i >= 0; --i) {
759                 if (j >= 0 && exp[j] == i) {
760                         if (i != j)
761                                 isl_int_swap(aff->v->el[offset + i],
762                                              aff->v->el[offset + j]);
763                         j--;
764                 } else
765                         isl_int_set_si(aff->v->el[offset + i], 0);
766         }
767
768         aff->ls = isl_local_space_replace_divs(aff->ls, isl_mat_copy(div));
769         if (!aff->ls)
770                 goto error;
771         isl_mat_free(div);
772         return aff;
773 error:
774         isl_aff_free(aff);
775         isl_mat_free(div);
776         return NULL;
777 }
778
779 /* Add two affine expressions that live in the same local space.
780  */
781 static __isl_give isl_aff *add_expanded(__isl_take isl_aff *aff1,
782         __isl_take isl_aff *aff2)
783 {
784         isl_int gcd, f;
785
786         aff1 = isl_aff_cow(aff1);
787         if (!aff1 || !aff2)
788                 goto error;
789
790         aff1->v = isl_vec_cow(aff1->v);
791         if (!aff1->v)
792                 goto error;
793
794         isl_int_init(gcd);
795         isl_int_init(f);
796         isl_int_gcd(gcd, aff1->v->el[0], aff2->v->el[0]);
797         isl_int_divexact(f, aff2->v->el[0], gcd);
798         isl_seq_scale(aff1->v->el + 1, aff1->v->el + 1, f, aff1->v->size - 1);
799         isl_int_divexact(f, aff1->v->el[0], gcd);
800         isl_seq_addmul(aff1->v->el + 1, f, aff2->v->el + 1, aff1->v->size - 1);
801         isl_int_divexact(f, aff2->v->el[0], gcd);
802         isl_int_mul(aff1->v->el[0], aff1->v->el[0], f);
803         isl_int_clear(f);
804         isl_int_clear(gcd);
805
806         isl_aff_free(aff2);
807         return aff1;
808 error:
809         isl_aff_free(aff1);
810         isl_aff_free(aff2);
811         return NULL;
812 }
813
814 __isl_give isl_aff *isl_aff_add(__isl_take isl_aff *aff1,
815         __isl_take isl_aff *aff2)
816 {
817         isl_ctx *ctx;
818         int *exp1 = NULL;
819         int *exp2 = NULL;
820         isl_mat *div;
821
822         if (!aff1 || !aff2)
823                 goto error;
824
825         ctx = isl_aff_get_ctx(aff1);
826         if (!isl_space_is_equal(aff1->ls->dim, aff2->ls->dim))
827                 isl_die(ctx, isl_error_invalid,
828                         "spaces don't match", goto error);
829
830         if (aff1->ls->div->n_row == 0 && aff2->ls->div->n_row == 0)
831                 return add_expanded(aff1, aff2);
832
833         exp1 = isl_alloc_array(ctx, int, aff1->ls->div->n_row);
834         exp2 = isl_alloc_array(ctx, int, aff2->ls->div->n_row);
835         if (!exp1 || !exp2)
836                 goto error;
837
838         div = isl_merge_divs(aff1->ls->div, aff2->ls->div, exp1, exp2);
839         aff1 = isl_aff_expand_divs(aff1, isl_mat_copy(div), exp1);
840         aff2 = isl_aff_expand_divs(aff2, div, exp2);
841         free(exp1);
842         free(exp2);
843
844         return add_expanded(aff1, aff2);
845 error:
846         free(exp1);
847         free(exp2);
848         isl_aff_free(aff1);
849         isl_aff_free(aff2);
850         return NULL;
851 }
852
853 __isl_give isl_aff *isl_aff_sub(__isl_take isl_aff *aff1,
854         __isl_take isl_aff *aff2)
855 {
856         return isl_aff_add(aff1, isl_aff_neg(aff2));
857 }
858
859 __isl_give isl_aff *isl_aff_scale(__isl_take isl_aff *aff, isl_int f)
860 {
861         isl_int gcd;
862
863         if (isl_int_is_one(f))
864                 return aff;
865
866         aff = isl_aff_cow(aff);
867         if (!aff)
868                 return NULL;
869         aff->v = isl_vec_cow(aff->v);
870         if (!aff->v)
871                 return isl_aff_free(aff);
872
873         isl_int_init(gcd);
874         isl_int_gcd(gcd, aff->v->el[0], f);
875         isl_int_divexact(aff->v->el[0], aff->v->el[0], gcd);
876         isl_int_divexact(gcd, f, gcd);
877         isl_seq_scale(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
878         isl_int_clear(gcd);
879
880         return aff;
881 }
882
883 __isl_give isl_aff *isl_aff_scale_down(__isl_take isl_aff *aff, isl_int f)
884 {
885         isl_int gcd;
886
887         if (isl_int_is_one(f))
888                 return aff;
889
890         aff = isl_aff_cow(aff);
891         if (!aff)
892                 return NULL;
893         aff->v = isl_vec_cow(aff->v);
894         if (!aff->v)
895                 return isl_aff_free(aff);
896
897         isl_int_init(gcd);
898         isl_seq_gcd(aff->v->el + 1, aff->v->size - 1, &gcd);
899         isl_int_gcd(gcd, gcd, f);
900         isl_seq_scale_down(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
901         isl_int_divexact(gcd, f, gcd);
902         isl_int_mul(aff->v->el[0], aff->v->el[0], gcd);
903         isl_int_clear(gcd);
904
905         return aff;
906 }
907
908 __isl_give isl_aff *isl_aff_scale_down_ui(__isl_take isl_aff *aff, unsigned f)
909 {
910         isl_int v;
911
912         if (f == 1)
913                 return aff;
914
915         isl_int_init(v);
916         isl_int_set_ui(v, f);
917         aff = isl_aff_scale_down(aff, v);
918         isl_int_clear(v);
919
920         return aff;
921 }
922
923 __isl_give isl_aff *isl_aff_set_dim_name(__isl_take isl_aff *aff,
924         enum isl_dim_type type, unsigned pos, const char *s)
925 {
926         aff = isl_aff_cow(aff);
927         if (!aff)
928                 return NULL;
929         if (type == isl_dim_out)
930                 isl_die(aff->v->ctx, isl_error_invalid,
931                         "cannot set name of output/set dimension",
932                         return isl_aff_free(aff));
933         if (type == isl_dim_in)
934                 type = isl_dim_set;
935         aff->ls = isl_local_space_set_dim_name(aff->ls, type, pos, s);
936         if (!aff->ls)
937                 return isl_aff_free(aff);
938
939         return aff;
940 }
941
942 __isl_give isl_aff *isl_aff_set_dim_id(__isl_take isl_aff *aff,
943         enum isl_dim_type type, unsigned pos, __isl_take isl_id *id)
944 {
945         aff = isl_aff_cow(aff);
946         if (!aff)
947                 return isl_id_free(id);
948         if (type == isl_dim_out)
949                 isl_die(aff->v->ctx, isl_error_invalid,
950                         "cannot set name of output/set dimension",
951                         goto error);
952         if (type == isl_dim_in)
953                 type = isl_dim_set;
954         aff->ls = isl_local_space_set_dim_id(aff->ls, type, pos, id);
955         if (!aff->ls)
956                 return isl_aff_free(aff);
957
958         return aff;
959 error:
960         isl_id_free(id);
961         isl_aff_free(aff);
962         return NULL;
963 }
964
965 /* Exploit the equalities in "eq" to simplify the affine expression
966  * and the expressions of the integer divisions in the local space.
967  * The integer divisions in this local space are assumed to appear
968  * as regular dimensions in "eq".
969  */
970 static __isl_give isl_aff *isl_aff_substitute_equalities_lifted(
971         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
972 {
973         int i, j;
974         unsigned total;
975         unsigned n_div;
976
977         if (!eq)
978                 goto error;
979         if (eq->n_eq == 0) {
980                 isl_basic_set_free(eq);
981                 return aff;
982         }
983
984         aff = isl_aff_cow(aff);
985         if (!aff)
986                 goto error;
987
988         aff->ls = isl_local_space_substitute_equalities(aff->ls,
989                                                         isl_basic_set_copy(eq));
990         if (!aff->ls)
991                 goto error;
992
993         total = 1 + isl_space_dim(eq->dim, isl_dim_all);
994         n_div = eq->n_div;
995         for (i = 0; i < eq->n_eq; ++i) {
996                 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
997                 if (j < 0 || j == 0 || j >= total)
998                         continue;
999
1000                 isl_seq_elim(aff->v->el + 1, eq->eq[i], j, total,
1001                                 &aff->v->el[0]);
1002         }
1003
1004         isl_basic_set_free(eq);
1005         aff = isl_aff_normalize(aff);
1006         return aff;
1007 error:
1008         isl_basic_set_free(eq);
1009         isl_aff_free(aff);
1010         return NULL;
1011 }
1012
1013 /* Exploit the equalities in "eq" to simplify the affine expression
1014  * and the expressions of the integer divisions in the local space.
1015  */
1016 static __isl_give isl_aff *isl_aff_substitute_equalities(
1017         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1018 {
1019         int n_div;
1020
1021         if (!aff || !eq)
1022                 goto error;
1023         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1024         if (n_div > 0)
1025                 eq = isl_basic_set_add(eq, isl_dim_set, n_div);
1026         return isl_aff_substitute_equalities_lifted(aff, eq);
1027 error:
1028         isl_basic_set_free(eq);
1029         isl_aff_free(aff);
1030         return NULL;
1031 }
1032
1033 /* Look for equalities among the variables shared by context and aff
1034  * and the integer divisions of aff, if any.
1035  * The equalities are then used to eliminate coefficients and/or integer
1036  * divisions from aff.
1037  */
1038 __isl_give isl_aff *isl_aff_gist(__isl_take isl_aff *aff,
1039         __isl_take isl_set *context)
1040 {
1041         isl_basic_set *hull;
1042         int n_div;
1043
1044         if (!aff)
1045                 goto error;
1046         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1047         if (n_div > 0) {
1048                 isl_basic_set *bset;
1049                 isl_local_space *ls;
1050                 context = isl_set_add_dims(context, isl_dim_set, n_div);
1051                 ls = isl_aff_get_domain_local_space(aff);
1052                 bset = isl_basic_set_from_local_space(ls);
1053                 bset = isl_basic_set_lift(bset);
1054                 bset = isl_basic_set_flatten(bset);
1055                 context = isl_set_intersect(context,
1056                                             isl_set_from_basic_set(bset));
1057         }
1058
1059         hull = isl_set_affine_hull(context);
1060         return isl_aff_substitute_equalities_lifted(aff, hull);
1061 error:
1062         isl_aff_free(aff);
1063         isl_set_free(context);
1064         return NULL;
1065 }
1066
1067 __isl_give isl_aff *isl_aff_gist_params(__isl_take isl_aff *aff,
1068         __isl_take isl_set *context)
1069 {
1070         isl_set *dom_context = isl_set_universe(isl_aff_get_domain_space(aff));
1071         dom_context = isl_set_intersect_params(dom_context, context);
1072         return isl_aff_gist(aff, dom_context);
1073 }
1074
1075 /* Return a basic set containing those elements in the space
1076  * of aff where it is non-negative.
1077  */
1078 __isl_give isl_basic_set *isl_aff_nonneg_basic_set(__isl_take isl_aff *aff)
1079 {
1080         isl_constraint *ineq;
1081         isl_basic_set *bset;
1082
1083         ineq = isl_inequality_from_aff(aff);
1084
1085         bset = isl_basic_set_from_constraint(ineq);
1086         bset = isl_basic_set_simplify(bset);
1087         return bset;
1088 }
1089
1090 /* Return a basic set containing those elements in the space
1091  * of aff where it is zero.
1092  */
1093 __isl_give isl_basic_set *isl_aff_zero_basic_set(__isl_take isl_aff *aff)
1094 {
1095         isl_constraint *ineq;
1096         isl_basic_set *bset;
1097
1098         ineq = isl_equality_from_aff(aff);
1099
1100         bset = isl_basic_set_from_constraint(ineq);
1101         bset = isl_basic_set_simplify(bset);
1102         return bset;
1103 }
1104
1105 /* Return a basic set containing those elements in the shared space
1106  * of aff1 and aff2 where aff1 is greater than or equal to aff2.
1107  */
1108 __isl_give isl_basic_set *isl_aff_ge_basic_set(__isl_take isl_aff *aff1,
1109         __isl_take isl_aff *aff2)
1110 {
1111         aff1 = isl_aff_sub(aff1, aff2);
1112
1113         return isl_aff_nonneg_basic_set(aff1);
1114 }
1115
1116 /* Return a basic set containing those elements in the shared space
1117  * of aff1 and aff2 where aff1 is smaller than or equal to aff2.
1118  */
1119 __isl_give isl_basic_set *isl_aff_le_basic_set(__isl_take isl_aff *aff1,
1120         __isl_take isl_aff *aff2)
1121 {
1122         return isl_aff_ge_basic_set(aff2, aff1);
1123 }
1124
1125 __isl_give isl_aff *isl_aff_add_on_domain(__isl_keep isl_set *dom,
1126         __isl_take isl_aff *aff1, __isl_take isl_aff *aff2)
1127 {
1128         aff1 = isl_aff_add(aff1, aff2);
1129         aff1 = isl_aff_gist(aff1, isl_set_copy(dom));
1130         return aff1;
1131 }
1132
1133 int isl_aff_is_empty(__isl_keep isl_aff *aff)
1134 {
1135         if (!aff)
1136                 return -1;
1137
1138         return 0;
1139 }
1140
1141 /* Check whether the given affine expression has non-zero coefficient
1142  * for any dimension in the given range or if any of these dimensions
1143  * appear with non-zero coefficients in any of the integer divisions
1144  * involved in the affine expression.
1145  */
1146 int isl_aff_involves_dims(__isl_keep isl_aff *aff,
1147         enum isl_dim_type type, unsigned first, unsigned n)
1148 {
1149         int i;
1150         isl_ctx *ctx;
1151         int *active = NULL;
1152         int involves = 0;
1153
1154         if (!aff)
1155                 return -1;
1156         if (n == 0)
1157                 return 0;
1158
1159         ctx = isl_aff_get_ctx(aff);
1160         if (first + n > isl_aff_dim(aff, type))
1161                 isl_die(ctx, isl_error_invalid,
1162                         "range out of bounds", return -1);
1163
1164         active = isl_local_space_get_active(aff->ls, aff->v->el + 2);
1165         if (!active)
1166                 goto error;
1167
1168         first += isl_local_space_offset(aff->ls, type) - 1;
1169         for (i = 0; i < n; ++i)
1170                 if (active[first + i]) {
1171                         involves = 1;
1172                         break;
1173                 }
1174
1175         free(active);
1176
1177         return involves;
1178 error:
1179         free(active);
1180         return -1;
1181 }
1182
1183 __isl_give isl_aff *isl_aff_drop_dims(__isl_take isl_aff *aff,
1184         enum isl_dim_type type, unsigned first, unsigned n)
1185 {
1186         isl_ctx *ctx;
1187
1188         if (!aff)
1189                 return NULL;
1190         if (type == isl_dim_out)
1191                 isl_die(aff->v->ctx, isl_error_invalid,
1192                         "cannot drop output/set dimension",
1193                         return isl_aff_free(aff));
1194         if (type == isl_dim_in)
1195                 type = isl_dim_set;
1196         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1197                 return aff;
1198
1199         ctx = isl_aff_get_ctx(aff);
1200         if (first + n > isl_local_space_dim(aff->ls, type))
1201                 isl_die(ctx, isl_error_invalid, "range out of bounds",
1202                         return isl_aff_free(aff));
1203
1204         aff = isl_aff_cow(aff);
1205         if (!aff)
1206                 return NULL;
1207
1208         aff->ls = isl_local_space_drop_dims(aff->ls, type, first, n);
1209         if (!aff->ls)
1210                 return isl_aff_free(aff);
1211
1212         first += 1 + isl_local_space_offset(aff->ls, type);
1213         aff->v = isl_vec_drop_els(aff->v, first, n);
1214         if (!aff->v)
1215                 return isl_aff_free(aff);
1216
1217         return aff;
1218 }
1219
1220 /* Project the domain of the affine expression onto its parameter space.
1221  * The affine expression may not involve any of the domain dimensions.
1222  */
1223 __isl_give isl_aff *isl_aff_project_domain_on_params(__isl_take isl_aff *aff)
1224 {
1225         isl_space *space;
1226         unsigned n;
1227         int involves;
1228
1229         n = isl_aff_dim(aff, isl_dim_in);
1230         involves = isl_aff_involves_dims(aff, isl_dim_in, 0, n);
1231         if (involves < 0)
1232                 return isl_aff_free(aff);
1233         if (involves)
1234                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1235                     "affine expression involves some of the domain dimensions",
1236                     return isl_aff_free(aff));
1237         aff = isl_aff_drop_dims(aff, isl_dim_in, 0, n);
1238         space = isl_aff_get_domain_space(aff);
1239         space = isl_space_params(space);
1240         aff = isl_aff_reset_domain_space(aff, space);
1241         return aff;
1242 }
1243
1244 __isl_give isl_aff *isl_aff_insert_dims(__isl_take isl_aff *aff,
1245         enum isl_dim_type type, unsigned first, unsigned n)
1246 {
1247         isl_ctx *ctx;
1248
1249         if (!aff)
1250                 return NULL;
1251         if (type == isl_dim_out)
1252                 isl_die(aff->v->ctx, isl_error_invalid,
1253                         "cannot insert output/set dimensions",
1254                         return isl_aff_free(aff));
1255         if (type == isl_dim_in)
1256                 type = isl_dim_set;
1257         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1258                 return aff;
1259
1260         ctx = isl_aff_get_ctx(aff);
1261         if (first > isl_local_space_dim(aff->ls, type))
1262                 isl_die(ctx, isl_error_invalid, "position out of bounds",
1263                         return isl_aff_free(aff));
1264
1265         aff = isl_aff_cow(aff);
1266         if (!aff)
1267                 return NULL;
1268
1269         aff->ls = isl_local_space_insert_dims(aff->ls, type, first, n);
1270         if (!aff->ls)
1271                 return isl_aff_free(aff);
1272
1273         first += 1 + isl_local_space_offset(aff->ls, type);
1274         aff->v = isl_vec_insert_zero_els(aff->v, first, n);
1275         if (!aff->v)
1276                 return isl_aff_free(aff);
1277
1278         return aff;
1279 }
1280
1281 __isl_give isl_aff *isl_aff_add_dims(__isl_take isl_aff *aff,
1282         enum isl_dim_type type, unsigned n)
1283 {
1284         unsigned pos;
1285
1286         pos = isl_aff_dim(aff, type);
1287
1288         return isl_aff_insert_dims(aff, type, pos, n);
1289 }
1290
1291 __isl_give isl_pw_aff *isl_pw_aff_add_dims(__isl_take isl_pw_aff *pwaff,
1292         enum isl_dim_type type, unsigned n)
1293 {
1294         unsigned pos;
1295
1296         pos = isl_pw_aff_dim(pwaff, type);
1297
1298         return isl_pw_aff_insert_dims(pwaff, type, pos, n);
1299 }
1300
1301 __isl_give isl_pw_aff *isl_pw_aff_from_aff(__isl_take isl_aff *aff)
1302 {
1303         isl_set *dom = isl_set_universe(isl_aff_get_domain_space(aff));
1304         return isl_pw_aff_alloc(dom, aff);
1305 }
1306
1307 #undef PW
1308 #define PW isl_pw_aff
1309 #undef EL
1310 #define EL isl_aff
1311 #undef EL_IS_ZERO
1312 #define EL_IS_ZERO is_empty
1313 #undef ZERO
1314 #define ZERO empty
1315 #undef IS_ZERO
1316 #define IS_ZERO is_empty
1317 #undef FIELD
1318 #define FIELD aff
1319 #undef DEFAULT_IS_ZERO
1320 #define DEFAULT_IS_ZERO 0
1321
1322 #define NO_EVAL
1323 #define NO_OPT
1324 #define NO_MOVE_DIMS
1325 #define NO_LIFT
1326 #define NO_MORPH
1327
1328 #include <isl_pw_templ.c>
1329
1330 static __isl_give isl_set *align_params_pw_pw_set_and(
1331         __isl_take isl_pw_aff *pwaff1, __isl_take isl_pw_aff *pwaff2,
1332         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1333                                     __isl_take isl_pw_aff *pwaff2))
1334 {
1335         if (!pwaff1 || !pwaff2)
1336                 goto error;
1337         if (isl_space_match(pwaff1->dim, isl_dim_param,
1338                           pwaff2->dim, isl_dim_param))
1339                 return fn(pwaff1, pwaff2);
1340         if (!isl_space_has_named_params(pwaff1->dim) ||
1341             !isl_space_has_named_params(pwaff2->dim))
1342                 isl_die(isl_pw_aff_get_ctx(pwaff1), isl_error_invalid,
1343                         "unaligned unnamed parameters", goto error);
1344         pwaff1 = isl_pw_aff_align_params(pwaff1, isl_pw_aff_get_space(pwaff2));
1345         pwaff2 = isl_pw_aff_align_params(pwaff2, isl_pw_aff_get_space(pwaff1));
1346         return fn(pwaff1, pwaff2);
1347 error:
1348         isl_pw_aff_free(pwaff1);
1349         isl_pw_aff_free(pwaff2);
1350         return NULL;
1351 }
1352
1353 /* Compute a piecewise quasi-affine expression with a domain that
1354  * is the union of those of pwaff1 and pwaff2 and such that on each
1355  * cell, the quasi-affine expression is the better (according to cmp)
1356  * of those of pwaff1 and pwaff2.  If only one of pwaff1 or pwaff2
1357  * is defined on a given cell, then the associated expression
1358  * is the defined one.
1359  */
1360 static __isl_give isl_pw_aff *pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1361         __isl_take isl_pw_aff *pwaff2,
1362         __isl_give isl_basic_set *(*cmp)(__isl_take isl_aff *aff1,
1363                                         __isl_take isl_aff *aff2))
1364 {
1365         int i, j, n;
1366         isl_pw_aff *res;
1367         isl_ctx *ctx;
1368         isl_set *set;
1369
1370         if (!pwaff1 || !pwaff2)
1371                 goto error;
1372
1373         ctx = isl_space_get_ctx(pwaff1->dim);
1374         if (!isl_space_is_equal(pwaff1->dim, pwaff2->dim))
1375                 isl_die(ctx, isl_error_invalid,
1376                         "arguments should live in same space", goto error);
1377
1378         if (isl_pw_aff_is_empty(pwaff1)) {
1379                 isl_pw_aff_free(pwaff1);
1380                 return pwaff2;
1381         }
1382
1383         if (isl_pw_aff_is_empty(pwaff2)) {
1384                 isl_pw_aff_free(pwaff2);
1385                 return pwaff1;
1386         }
1387
1388         n = 2 * (pwaff1->n + 1) * (pwaff2->n + 1);
1389         res = isl_pw_aff_alloc_size(isl_space_copy(pwaff1->dim), n);
1390
1391         for (i = 0; i < pwaff1->n; ++i) {
1392                 set = isl_set_copy(pwaff1->p[i].set);
1393                 for (j = 0; j < pwaff2->n; ++j) {
1394                         struct isl_set *common;
1395                         isl_set *better;
1396
1397                         common = isl_set_intersect(
1398                                         isl_set_copy(pwaff1->p[i].set),
1399                                         isl_set_copy(pwaff2->p[j].set));
1400                         better = isl_set_from_basic_set(cmp(
1401                                         isl_aff_copy(pwaff2->p[j].aff),
1402                                         isl_aff_copy(pwaff1->p[i].aff)));
1403                         better = isl_set_intersect(common, better);
1404                         if (isl_set_plain_is_empty(better)) {
1405                                 isl_set_free(better);
1406                                 continue;
1407                         }
1408                         set = isl_set_subtract(set, isl_set_copy(better));
1409
1410                         res = isl_pw_aff_add_piece(res, better,
1411                                                 isl_aff_copy(pwaff2->p[j].aff));
1412                 }
1413                 res = isl_pw_aff_add_piece(res, set,
1414                                                 isl_aff_copy(pwaff1->p[i].aff));
1415         }
1416
1417         for (j = 0; j < pwaff2->n; ++j) {
1418                 set = isl_set_copy(pwaff2->p[j].set);
1419                 for (i = 0; i < pwaff1->n; ++i)
1420                         set = isl_set_subtract(set,
1421                                         isl_set_copy(pwaff1->p[i].set));
1422                 res = isl_pw_aff_add_piece(res, set,
1423                                                 isl_aff_copy(pwaff2->p[j].aff));
1424         }
1425
1426         isl_pw_aff_free(pwaff1);
1427         isl_pw_aff_free(pwaff2);
1428
1429         return res;
1430 error:
1431         isl_pw_aff_free(pwaff1);
1432         isl_pw_aff_free(pwaff2);
1433         return NULL;
1434 }
1435
1436 /* Compute a piecewise quasi-affine expression with a domain that
1437  * is the union of those of pwaff1 and pwaff2 and such that on each
1438  * cell, the quasi-affine expression is the maximum of those of pwaff1
1439  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1440  * cell, then the associated expression is the defined one.
1441  */
1442 static __isl_give isl_pw_aff *pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1443         __isl_take isl_pw_aff *pwaff2)
1444 {
1445         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_ge_basic_set);
1446 }
1447
1448 __isl_give isl_pw_aff *isl_pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1449         __isl_take isl_pw_aff *pwaff2)
1450 {
1451         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1452                                                         &pw_aff_union_max);
1453 }
1454
1455 /* Compute a piecewise quasi-affine expression with a domain that
1456  * is the union of those of pwaff1 and pwaff2 and such that on each
1457  * cell, the quasi-affine expression is the minimum of those of pwaff1
1458  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1459  * cell, then the associated expression is the defined one.
1460  */
1461 static __isl_give isl_pw_aff *pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1462         __isl_take isl_pw_aff *pwaff2)
1463 {
1464         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_le_basic_set);
1465 }
1466
1467 __isl_give isl_pw_aff *isl_pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1468         __isl_take isl_pw_aff *pwaff2)
1469 {
1470         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1471                                                         &pw_aff_union_min);
1472 }
1473
1474 __isl_give isl_pw_aff *isl_pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1475         __isl_take isl_pw_aff *pwaff2, int max)
1476 {
1477         if (max)
1478                 return isl_pw_aff_union_max(pwaff1, pwaff2);
1479         else
1480                 return isl_pw_aff_union_min(pwaff1, pwaff2);
1481 }
1482
1483 /* Construct a map with as domain the domain of pwaff and
1484  * one-dimensional range corresponding to the affine expressions.
1485  */
1486 static __isl_give isl_map *map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1487 {
1488         int i;
1489         isl_space *dim;
1490         isl_map *map;
1491
1492         if (!pwaff)
1493                 return NULL;
1494
1495         dim = isl_pw_aff_get_space(pwaff);
1496         map = isl_map_empty(dim);
1497
1498         for (i = 0; i < pwaff->n; ++i) {
1499                 isl_basic_map *bmap;
1500                 isl_map *map_i;
1501
1502                 bmap = isl_basic_map_from_aff(isl_aff_copy(pwaff->p[i].aff));
1503                 map_i = isl_map_from_basic_map(bmap);
1504                 map_i = isl_map_intersect_domain(map_i,
1505                                                 isl_set_copy(pwaff->p[i].set));
1506                 map = isl_map_union_disjoint(map, map_i);
1507         }
1508
1509         isl_pw_aff_free(pwaff);
1510
1511         return map;
1512 }
1513
1514 /* Construct a map with as domain the domain of pwaff and
1515  * one-dimensional range corresponding to the affine expressions.
1516  */
1517 __isl_give isl_map *isl_map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1518 {
1519         if (!pwaff)
1520                 return NULL;
1521         if (isl_space_is_set(pwaff->dim))
1522                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1523                         "space of input is not a map",
1524                         return isl_pw_aff_free(pwaff));
1525         return map_from_pw_aff(pwaff);
1526 }
1527
1528 /* Construct a one-dimensional set with as parameter domain
1529  * the domain of pwaff and the single set dimension
1530  * corresponding to the affine expressions.
1531  */
1532 __isl_give isl_set *isl_set_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1533 {
1534         if (!pwaff)
1535                 return NULL;
1536         if (!isl_space_is_set(pwaff->dim))
1537                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1538                         "space of input is not a set",
1539                         return isl_pw_aff_free(pwaff));
1540         return map_from_pw_aff(pwaff);
1541 }
1542
1543 /* Return a set containing those elements in the domain
1544  * of pwaff where it is non-negative.
1545  */
1546 __isl_give isl_set *isl_pw_aff_nonneg_set(__isl_take isl_pw_aff *pwaff)
1547 {
1548         int i;
1549         isl_set *set;
1550
1551         if (!pwaff)
1552                 return NULL;
1553
1554         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1555
1556         for (i = 0; i < pwaff->n; ++i) {
1557                 isl_basic_set *bset;
1558                 isl_set *set_i;
1559
1560                 bset = isl_aff_nonneg_basic_set(isl_aff_copy(pwaff->p[i].aff));
1561                 set_i = isl_set_from_basic_set(bset);
1562                 set_i = isl_set_intersect(set_i, isl_set_copy(pwaff->p[i].set));
1563                 set = isl_set_union_disjoint(set, set_i);
1564         }
1565
1566         isl_pw_aff_free(pwaff);
1567
1568         return set;
1569 }
1570
1571 /* Return a set containing those elements in the domain
1572  * of pwaff where it is zero (if complement is 0) or not zero
1573  * (if complement is 1).
1574  */
1575 static __isl_give isl_set *pw_aff_zero_set(__isl_take isl_pw_aff *pwaff,
1576         int complement)
1577 {
1578         int i;
1579         isl_set *set;
1580
1581         if (!pwaff)
1582                 return NULL;
1583
1584         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1585
1586         for (i = 0; i < pwaff->n; ++i) {
1587                 isl_basic_set *bset;
1588                 isl_set *set_i, *zero;
1589
1590                 bset = isl_aff_zero_basic_set(isl_aff_copy(pwaff->p[i].aff));
1591                 zero = isl_set_from_basic_set(bset);
1592                 set_i = isl_set_copy(pwaff->p[i].set);
1593                 if (complement)
1594                         set_i = isl_set_subtract(set_i, zero);
1595                 else
1596                         set_i = isl_set_intersect(set_i, zero);
1597                 set = isl_set_union_disjoint(set, set_i);
1598         }
1599
1600         isl_pw_aff_free(pwaff);
1601
1602         return set;
1603 }
1604
1605 /* Return a set containing those elements in the domain
1606  * of pwaff where it is zero.
1607  */
1608 __isl_give isl_set *isl_pw_aff_zero_set(__isl_take isl_pw_aff *pwaff)
1609 {
1610         return pw_aff_zero_set(pwaff, 0);
1611 }
1612
1613 /* Return a set containing those elements in the domain
1614  * of pwaff where it is not zero.
1615  */
1616 __isl_give isl_set *isl_pw_aff_non_zero_set(__isl_take isl_pw_aff *pwaff)
1617 {
1618         return pw_aff_zero_set(pwaff, 1);
1619 }
1620
1621 /* Return a set containing those elements in the shared domain
1622  * of pwaff1 and pwaff2 where pwaff1 is greater than (or equal) to pwaff2.
1623  *
1624  * We compute the difference on the shared domain and then construct
1625  * the set of values where this difference is non-negative.
1626  * If strict is set, we first subtract 1 from the difference.
1627  * If equal is set, we only return the elements where pwaff1 and pwaff2
1628  * are equal.
1629  */
1630 static __isl_give isl_set *pw_aff_gte_set(__isl_take isl_pw_aff *pwaff1,
1631         __isl_take isl_pw_aff *pwaff2, int strict, int equal)
1632 {
1633         isl_set *set1, *set2;
1634
1635         set1 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff1));
1636         set2 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff2));
1637         set1 = isl_set_intersect(set1, set2);
1638         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, isl_set_copy(set1));
1639         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, isl_set_copy(set1));
1640         pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_neg(pwaff2));
1641
1642         if (strict) {
1643                 isl_space *dim = isl_set_get_space(set1);
1644                 isl_aff *aff;
1645                 aff = isl_aff_zero_on_domain(isl_local_space_from_space(dim));
1646                 aff = isl_aff_add_constant_si(aff, -1);
1647                 pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_alloc(set1, aff));
1648         } else
1649                 isl_set_free(set1);
1650
1651         if (equal)
1652                 return isl_pw_aff_zero_set(pwaff1);
1653         return isl_pw_aff_nonneg_set(pwaff1);
1654 }
1655
1656 /* Return a set containing those elements in the shared domain
1657  * of pwaff1 and pwaff2 where pwaff1 is equal to pwaff2.
1658  */
1659 static __isl_give isl_set *pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1660         __isl_take isl_pw_aff *pwaff2)
1661 {
1662         return pw_aff_gte_set(pwaff1, pwaff2, 0, 1);
1663 }
1664
1665 __isl_give isl_set *isl_pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1666         __isl_take isl_pw_aff *pwaff2)
1667 {
1668         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_eq_set);
1669 }
1670
1671 /* Return a set containing those elements in the shared domain
1672  * of pwaff1 and pwaff2 where pwaff1 is greater than or equal to pwaff2.
1673  */
1674 static __isl_give isl_set *pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1675         __isl_take isl_pw_aff *pwaff2)
1676 {
1677         return pw_aff_gte_set(pwaff1, pwaff2, 0, 0);
1678 }
1679
1680 __isl_give isl_set *isl_pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1681         __isl_take isl_pw_aff *pwaff2)
1682 {
1683         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ge_set);
1684 }
1685
1686 /* Return a set containing those elements in the shared domain
1687  * of pwaff1 and pwaff2 where pwaff1 is strictly greater than pwaff2.
1688  */
1689 static __isl_give isl_set *pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1690         __isl_take isl_pw_aff *pwaff2)
1691 {
1692         return pw_aff_gte_set(pwaff1, pwaff2, 1, 0);
1693 }
1694
1695 __isl_give isl_set *isl_pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1696         __isl_take isl_pw_aff *pwaff2)
1697 {
1698         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_gt_set);
1699 }
1700
1701 __isl_give isl_set *isl_pw_aff_le_set(__isl_take isl_pw_aff *pwaff1,
1702         __isl_take isl_pw_aff *pwaff2)
1703 {
1704         return isl_pw_aff_ge_set(pwaff2, pwaff1);
1705 }
1706
1707 __isl_give isl_set *isl_pw_aff_lt_set(__isl_take isl_pw_aff *pwaff1,
1708         __isl_take isl_pw_aff *pwaff2)
1709 {
1710         return isl_pw_aff_gt_set(pwaff2, pwaff1);
1711 }
1712
1713 /* Return a set containing those elements in the shared domain
1714  * of the elements of list1 and list2 where each element in list1
1715  * has the relation specified by "fn" with each element in list2.
1716  */
1717 static __isl_give isl_set *pw_aff_list_set(__isl_take isl_pw_aff_list *list1,
1718         __isl_take isl_pw_aff_list *list2,
1719         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1720                                     __isl_take isl_pw_aff *pwaff2))
1721 {
1722         int i, j;
1723         isl_ctx *ctx;
1724         isl_set *set;
1725
1726         if (!list1 || !list2)
1727                 goto error;
1728
1729         ctx = isl_pw_aff_list_get_ctx(list1);
1730         if (list1->n < 1 || list2->n < 1)
1731                 isl_die(ctx, isl_error_invalid,
1732                         "list should contain at least one element", goto error);
1733
1734         set = isl_set_universe(isl_pw_aff_get_domain_space(list1->p[0]));
1735         for (i = 0; i < list1->n; ++i)
1736                 for (j = 0; j < list2->n; ++j) {
1737                         isl_set *set_ij;
1738
1739                         set_ij = fn(isl_pw_aff_copy(list1->p[i]),
1740                                     isl_pw_aff_copy(list2->p[j]));
1741                         set = isl_set_intersect(set, set_ij);
1742                 }
1743
1744         isl_pw_aff_list_free(list1);
1745         isl_pw_aff_list_free(list2);
1746         return set;
1747 error:
1748         isl_pw_aff_list_free(list1);
1749         isl_pw_aff_list_free(list2);
1750         return NULL;
1751 }
1752
1753 /* Return a set containing those elements in the shared domain
1754  * of the elements of list1 and list2 where each element in list1
1755  * is equal to each element in list2.
1756  */
1757 __isl_give isl_set *isl_pw_aff_list_eq_set(__isl_take isl_pw_aff_list *list1,
1758         __isl_take isl_pw_aff_list *list2)
1759 {
1760         return pw_aff_list_set(list1, list2, &isl_pw_aff_eq_set);
1761 }
1762
1763 __isl_give isl_set *isl_pw_aff_list_ne_set(__isl_take isl_pw_aff_list *list1,
1764         __isl_take isl_pw_aff_list *list2)
1765 {
1766         return pw_aff_list_set(list1, list2, &isl_pw_aff_ne_set);
1767 }
1768
1769 /* Return a set containing those elements in the shared domain
1770  * of the elements of list1 and list2 where each element in list1
1771  * is less than or equal to each element in list2.
1772  */
1773 __isl_give isl_set *isl_pw_aff_list_le_set(__isl_take isl_pw_aff_list *list1,
1774         __isl_take isl_pw_aff_list *list2)
1775 {
1776         return pw_aff_list_set(list1, list2, &isl_pw_aff_le_set);
1777 }
1778
1779 __isl_give isl_set *isl_pw_aff_list_lt_set(__isl_take isl_pw_aff_list *list1,
1780         __isl_take isl_pw_aff_list *list2)
1781 {
1782         return pw_aff_list_set(list1, list2, &isl_pw_aff_lt_set);
1783 }
1784
1785 __isl_give isl_set *isl_pw_aff_list_ge_set(__isl_take isl_pw_aff_list *list1,
1786         __isl_take isl_pw_aff_list *list2)
1787 {
1788         return pw_aff_list_set(list1, list2, &isl_pw_aff_ge_set);
1789 }
1790
1791 __isl_give isl_set *isl_pw_aff_list_gt_set(__isl_take isl_pw_aff_list *list1,
1792         __isl_take isl_pw_aff_list *list2)
1793 {
1794         return pw_aff_list_set(list1, list2, &isl_pw_aff_gt_set);
1795 }
1796
1797
1798 /* Return a set containing those elements in the shared domain
1799  * of pwaff1 and pwaff2 where pwaff1 is not equal to pwaff2.
1800  */
1801 static __isl_give isl_set *pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
1802         __isl_take isl_pw_aff *pwaff2)
1803 {
1804         isl_set *set_lt, *set_gt;
1805
1806         set_lt = isl_pw_aff_lt_set(isl_pw_aff_copy(pwaff1),
1807                                    isl_pw_aff_copy(pwaff2));
1808         set_gt = isl_pw_aff_gt_set(pwaff1, pwaff2);
1809         return isl_set_union_disjoint(set_lt, set_gt);
1810 }
1811
1812 __isl_give isl_set *isl_pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
1813         __isl_take isl_pw_aff *pwaff2)
1814 {
1815         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ne_set);
1816 }
1817
1818 __isl_give isl_pw_aff *isl_pw_aff_scale_down(__isl_take isl_pw_aff *pwaff,
1819         isl_int v)
1820 {
1821         int i;
1822
1823         if (isl_int_is_one(v))
1824                 return pwaff;
1825         if (!isl_int_is_pos(v))
1826                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1827                         "factor needs to be positive",
1828                         return isl_pw_aff_free(pwaff));
1829         pwaff = isl_pw_aff_cow(pwaff);
1830         if (!pwaff)
1831                 return NULL;
1832         if (pwaff->n == 0)
1833                 return pwaff;
1834
1835         for (i = 0; i < pwaff->n; ++i) {
1836                 pwaff->p[i].aff = isl_aff_scale_down(pwaff->p[i].aff, v);
1837                 if (!pwaff->p[i].aff)
1838                         return isl_pw_aff_free(pwaff);
1839         }
1840
1841         return pwaff;
1842 }
1843
1844 __isl_give isl_pw_aff *isl_pw_aff_floor(__isl_take isl_pw_aff *pwaff)
1845 {
1846         int i;
1847
1848         pwaff = isl_pw_aff_cow(pwaff);
1849         if (!pwaff)
1850                 return NULL;
1851         if (pwaff->n == 0)
1852                 return pwaff;
1853
1854         for (i = 0; i < pwaff->n; ++i) {
1855                 pwaff->p[i].aff = isl_aff_floor(pwaff->p[i].aff);
1856                 if (!pwaff->p[i].aff)
1857                         return isl_pw_aff_free(pwaff);
1858         }
1859
1860         return pwaff;
1861 }
1862
1863 __isl_give isl_pw_aff *isl_pw_aff_ceil(__isl_take isl_pw_aff *pwaff)
1864 {
1865         int i;
1866
1867         pwaff = isl_pw_aff_cow(pwaff);
1868         if (!pwaff)
1869                 return NULL;
1870         if (pwaff->n == 0)
1871                 return pwaff;
1872
1873         for (i = 0; i < pwaff->n; ++i) {
1874                 pwaff->p[i].aff = isl_aff_ceil(pwaff->p[i].aff);
1875                 if (!pwaff->p[i].aff)
1876                         return isl_pw_aff_free(pwaff);
1877         }
1878
1879         return pwaff;
1880 }
1881
1882 /* Return an affine expression that is equal to pwaff_true for elements
1883  * in "cond" and to pwaff_false for elements not in "cond".
1884  * That is, return cond ? pwaff_true : pwaff_false;
1885  */
1886 __isl_give isl_pw_aff *isl_pw_aff_cond(__isl_take isl_set *cond,
1887         __isl_take isl_pw_aff *pwaff_true, __isl_take isl_pw_aff *pwaff_false)
1888 {
1889         isl_set *comp;
1890
1891         comp = isl_set_complement(isl_set_copy(cond));
1892         pwaff_true = isl_pw_aff_intersect_domain(pwaff_true, cond);
1893         pwaff_false = isl_pw_aff_intersect_domain(pwaff_false, comp);
1894
1895         return isl_pw_aff_add_disjoint(pwaff_true, pwaff_false);
1896 }
1897
1898 int isl_aff_is_cst(__isl_keep isl_aff *aff)
1899 {
1900         if (!aff)
1901                 return -1;
1902
1903         return isl_seq_first_non_zero(aff->v->el + 2, aff->v->size - 2) == -1;
1904 }
1905
1906 /* Check whether pwaff is a piecewise constant.
1907  */
1908 int isl_pw_aff_is_cst(__isl_keep isl_pw_aff *pwaff)
1909 {
1910         int i;
1911
1912         if (!pwaff)
1913                 return -1;
1914
1915         for (i = 0; i < pwaff->n; ++i) {
1916                 int is_cst = isl_aff_is_cst(pwaff->p[i].aff);
1917                 if (is_cst < 0 || !is_cst)
1918                         return is_cst;
1919         }
1920
1921         return 1;
1922 }
1923
1924 __isl_give isl_aff *isl_aff_mul(__isl_take isl_aff *aff1,
1925         __isl_take isl_aff *aff2)
1926 {
1927         if (!isl_aff_is_cst(aff2) && isl_aff_is_cst(aff1))
1928                 return isl_aff_mul(aff2, aff1);
1929
1930         if (!isl_aff_is_cst(aff2))
1931                 isl_die(isl_aff_get_ctx(aff1), isl_error_invalid,
1932                         "at least one affine expression should be constant",
1933                         goto error);
1934
1935         aff1 = isl_aff_cow(aff1);
1936         if (!aff1 || !aff2)
1937                 goto error;
1938
1939         aff1 = isl_aff_scale(aff1, aff2->v->el[1]);
1940         aff1 = isl_aff_scale_down(aff1, aff2->v->el[0]);
1941
1942         isl_aff_free(aff2);
1943         return aff1;
1944 error:
1945         isl_aff_free(aff1);
1946         isl_aff_free(aff2);
1947         return NULL;
1948 }
1949
1950 static __isl_give isl_pw_aff *pw_aff_add(__isl_take isl_pw_aff *pwaff1,
1951         __isl_take isl_pw_aff *pwaff2)
1952 {
1953         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_add);
1954 }
1955
1956 __isl_give isl_pw_aff *isl_pw_aff_add(__isl_take isl_pw_aff *pwaff1,
1957         __isl_take isl_pw_aff *pwaff2)
1958 {
1959         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_add);
1960 }
1961
1962 __isl_give isl_pw_aff *isl_pw_aff_union_add(__isl_take isl_pw_aff *pwaff1,
1963         __isl_take isl_pw_aff *pwaff2)
1964 {
1965         return isl_pw_aff_union_add_(pwaff1, pwaff2);
1966 }
1967
1968 static __isl_give isl_pw_aff *pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
1969         __isl_take isl_pw_aff *pwaff2)
1970 {
1971         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_mul);
1972 }
1973
1974 __isl_give isl_pw_aff *isl_pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
1975         __isl_take isl_pw_aff *pwaff2)
1976 {
1977         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_mul);
1978 }
1979
1980 static __isl_give isl_pw_aff *pw_aff_min(__isl_take isl_pw_aff *pwaff1,
1981         __isl_take isl_pw_aff *pwaff2)
1982 {
1983         isl_set *le;
1984
1985         le = isl_pw_aff_le_set(isl_pw_aff_copy(pwaff1),
1986                                 isl_pw_aff_copy(pwaff2));
1987         return isl_pw_aff_cond(le, pwaff1, pwaff2);
1988 }
1989
1990 __isl_give isl_pw_aff *isl_pw_aff_min(__isl_take isl_pw_aff *pwaff1,
1991         __isl_take isl_pw_aff *pwaff2)
1992 {
1993         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_min);
1994 }
1995
1996 static __isl_give isl_pw_aff *pw_aff_max(__isl_take isl_pw_aff *pwaff1,
1997         __isl_take isl_pw_aff *pwaff2)
1998 {
1999         isl_set *le;
2000
2001         le = isl_pw_aff_ge_set(isl_pw_aff_copy(pwaff1),
2002                                 isl_pw_aff_copy(pwaff2));
2003         return isl_pw_aff_cond(le, pwaff1, pwaff2);
2004 }
2005
2006 __isl_give isl_pw_aff *isl_pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2007         __isl_take isl_pw_aff *pwaff2)
2008 {
2009         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_max);
2010 }
2011
2012 static __isl_give isl_pw_aff *pw_aff_list_reduce(
2013         __isl_take isl_pw_aff_list *list,
2014         __isl_give isl_pw_aff *(*fn)(__isl_take isl_pw_aff *pwaff1,
2015                                         __isl_take isl_pw_aff *pwaff2))
2016 {
2017         int i;
2018         isl_ctx *ctx;
2019         isl_pw_aff *res;
2020
2021         if (!list)
2022                 return NULL;
2023
2024         ctx = isl_pw_aff_list_get_ctx(list);
2025         if (list->n < 1)
2026                 isl_die(ctx, isl_error_invalid,
2027                         "list should contain at least one element",
2028                         return isl_pw_aff_list_free(list));
2029
2030         res = isl_pw_aff_copy(list->p[0]);
2031         for (i = 1; i < list->n; ++i)
2032                 res = fn(res, isl_pw_aff_copy(list->p[i]));
2033
2034         isl_pw_aff_list_free(list);
2035         return res;
2036 }
2037
2038 /* Return an isl_pw_aff that maps each element in the intersection of the
2039  * domains of the elements of list to the minimal corresponding affine
2040  * expression.
2041  */
2042 __isl_give isl_pw_aff *isl_pw_aff_list_min(__isl_take isl_pw_aff_list *list)
2043 {
2044         return pw_aff_list_reduce(list, &isl_pw_aff_min);
2045 }
2046
2047 /* Return an isl_pw_aff that maps each element in the intersection of the
2048  * domains of the elements of list to the maximal corresponding affine
2049  * expression.
2050  */
2051 __isl_give isl_pw_aff *isl_pw_aff_list_max(__isl_take isl_pw_aff_list *list)
2052 {
2053         return pw_aff_list_reduce(list, &isl_pw_aff_max);
2054 }
2055
2056 #undef BASE
2057 #define BASE aff
2058
2059 #include <isl_multi_templ.c>
2060
2061 __isl_give isl_multi_aff *isl_multi_aff_add(__isl_take isl_multi_aff *maff1,
2062         __isl_take isl_multi_aff *maff2)
2063 {
2064         int i;
2065         isl_ctx *ctx;
2066
2067         maff1 = isl_multi_aff_cow(maff1);
2068         if (!maff1 || !maff2)
2069                 goto error;
2070
2071         ctx = isl_multi_aff_get_ctx(maff1);
2072         if (!isl_space_is_equal(maff1->space, maff2->space))
2073                 isl_die(ctx, isl_error_invalid,
2074                         "spaces don't match", goto error);
2075
2076         for (i = 0; i < maff1->n; ++i) {
2077                 maff1->p[i] = isl_aff_add(maff1->p[i],
2078                                             isl_aff_copy(maff2->p[i]));
2079                 if (!maff1->p[i])
2080                         goto error;
2081         }
2082
2083         isl_multi_aff_free(maff2);
2084         return maff1;
2085 error:
2086         isl_multi_aff_free(maff1);
2087         isl_multi_aff_free(maff2);
2088         return NULL;
2089 }
2090
2091 /* Exploit the equalities in "eq" to simplify the affine expressions.
2092  */
2093 static __isl_give isl_multi_aff *isl_multi_aff_substitute_equalities(
2094         __isl_take isl_multi_aff *maff, __isl_take isl_basic_set *eq)
2095 {
2096         int i;
2097
2098         maff = isl_multi_aff_cow(maff);
2099         if (!maff || !eq)
2100                 goto error;
2101
2102         for (i = 0; i < maff->n; ++i) {
2103                 maff->p[i] = isl_aff_substitute_equalities(maff->p[i],
2104                                                     isl_basic_set_copy(eq));
2105                 if (!maff->p[i])
2106                         goto error;
2107         }
2108
2109         isl_basic_set_free(eq);
2110         return maff;
2111 error:
2112         isl_basic_set_free(eq);
2113         isl_multi_aff_free(maff);
2114         return NULL;
2115 }
2116
2117 __isl_give isl_multi_aff *isl_multi_aff_scale(__isl_take isl_multi_aff *maff,
2118         isl_int f)
2119 {
2120         int i;
2121
2122         maff = isl_multi_aff_cow(maff);
2123         if (!maff)
2124                 return NULL;
2125
2126         for (i = 0; i < maff->n; ++i) {
2127                 maff->p[i] = isl_aff_scale(maff->p[i], f);
2128                 if (!maff->p[i])
2129                         return isl_multi_aff_free(maff);
2130         }
2131
2132         return maff;
2133 }
2134
2135 __isl_give isl_multi_aff *isl_multi_aff_add_on_domain(__isl_keep isl_set *dom,
2136         __isl_take isl_multi_aff *maff1, __isl_take isl_multi_aff *maff2)
2137 {
2138         maff1 = isl_multi_aff_add(maff1, maff2);
2139         maff1 = isl_multi_aff_gist(maff1, isl_set_copy(dom));
2140         return maff1;
2141 }
2142
2143 int isl_multi_aff_is_empty(__isl_keep isl_multi_aff *maff)
2144 {
2145         if (!maff)
2146                 return -1;
2147
2148         return 0;
2149 }
2150
2151 int isl_multi_aff_plain_is_equal(__isl_keep isl_multi_aff *maff1,
2152         __isl_keep isl_multi_aff *maff2)
2153 {
2154         int i;
2155         int equal;
2156
2157         if (!maff1 || !maff2)
2158                 return -1;
2159         if (maff1->n != maff2->n)
2160                 return 0;
2161         equal = isl_space_is_equal(maff1->space, maff2->space);
2162         if (equal < 0 || !equal)
2163                 return equal;
2164
2165         for (i = 0; i < maff1->n; ++i) {
2166                 equal = isl_aff_plain_is_equal(maff1->p[i], maff2->p[i]);
2167                 if (equal < 0 || !equal)
2168                         return equal;
2169         }
2170
2171         return 1;
2172 }
2173
2174 __isl_give isl_multi_aff *isl_multi_aff_set_dim_name(
2175         __isl_take isl_multi_aff *maff,
2176         enum isl_dim_type type, unsigned pos, const char *s)
2177 {
2178         int i;
2179
2180         maff = isl_multi_aff_cow(maff);
2181         if (!maff)
2182                 return NULL;
2183
2184         maff->space = isl_space_set_dim_name(maff->space, type, pos, s);
2185         if (!maff->space)
2186                 return isl_multi_aff_free(maff);
2187         for (i = 0; i < maff->n; ++i) {
2188                 maff->p[i] = isl_aff_set_dim_name(maff->p[i], type, pos, s);
2189                 if (!maff->p[i])
2190                         return isl_multi_aff_free(maff);
2191         }
2192
2193         return maff;
2194 }
2195
2196 __isl_give isl_multi_aff *isl_multi_aff_drop_dims(__isl_take isl_multi_aff *maff,
2197         enum isl_dim_type type, unsigned first, unsigned n)
2198 {
2199         int i;
2200
2201         maff = isl_multi_aff_cow(maff);
2202         if (!maff)
2203                 return NULL;
2204
2205         maff->space = isl_space_drop_dims(maff->space, type, first, n);
2206         if (!maff->space)
2207                 return isl_multi_aff_free(maff);
2208
2209         if (type == isl_dim_out) {
2210                 for (i = 0; i < n; ++i)
2211                         isl_aff_free(maff->p[first + i]);
2212                 for (i = first; i + n < maff->n; ++i)
2213                         maff->p[i] = maff->p[i + n];
2214                 maff->n -= n;
2215                 return maff;
2216         }
2217
2218         for (i = 0; i < maff->n; ++i) {
2219                 maff->p[i] = isl_aff_drop_dims(maff->p[i], type, first, n);
2220                 if (!maff->p[i])
2221                         return isl_multi_aff_free(maff);
2222         }
2223
2224         return maff;
2225 }
2226
2227 #undef PW
2228 #define PW isl_pw_multi_aff
2229 #undef EL
2230 #define EL isl_multi_aff
2231 #undef EL_IS_ZERO
2232 #define EL_IS_ZERO is_empty
2233 #undef ZERO
2234 #define ZERO empty
2235 #undef IS_ZERO
2236 #define IS_ZERO is_empty
2237 #undef FIELD
2238 #define FIELD maff
2239 #undef DEFAULT_IS_ZERO
2240 #define DEFAULT_IS_ZERO 0
2241
2242 #define NO_NEG
2243 #define NO_EVAL
2244 #define NO_OPT
2245 #define NO_INVOLVES_DIMS
2246 #define NO_MOVE_DIMS
2247 #define NO_INSERT_DIMS
2248 #define NO_LIFT
2249 #define NO_MORPH
2250
2251 #include <isl_pw_templ.c>
2252
2253 static __isl_give isl_pw_multi_aff *pw_multi_aff_add(
2254         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2255 {
2256         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
2257                                                 &isl_multi_aff_add);
2258 }
2259
2260 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_add(
2261         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2262 {
2263         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2264                                                 &pw_multi_aff_add);
2265 }
2266
2267 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_add(
2268         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2269 {
2270         return isl_pw_multi_aff_union_add_(pma1, pma2);
2271 }
2272
2273 /* Construct a map mapping the domain the piecewise multi-affine expression
2274  * to its range, with each dimension in the range equated to the
2275  * corresponding affine expression on its cell.
2276  */
2277 __isl_give isl_map *isl_map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
2278 {
2279         int i;
2280         isl_map *map;
2281
2282         if (!pma)
2283                 return NULL;
2284
2285         map = isl_map_empty(isl_pw_multi_aff_get_space(pma));
2286
2287         for (i = 0; i < pma->n; ++i) {
2288                 isl_multi_aff *maff;
2289                 isl_basic_map *bmap;
2290                 isl_map *map_i;
2291
2292                 maff = isl_multi_aff_copy(pma->p[i].maff);
2293                 bmap = isl_basic_map_from_multi_aff(maff);
2294                 map_i = isl_map_from_basic_map(bmap);
2295                 map_i = isl_map_intersect_domain(map_i,
2296                                                 isl_set_copy(pma->p[i].set));
2297                 map = isl_map_union_disjoint(map, map_i);
2298         }
2299
2300         isl_pw_multi_aff_free(pma);
2301         return map;
2302 }
2303
2304 __isl_give isl_set *isl_set_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
2305 {
2306         if (!isl_space_is_set(pma->dim))
2307                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
2308                         "isl_pw_multi_aff cannot be converted into an isl_set",
2309                         return isl_pw_multi_aff_free(pma));
2310
2311         return isl_map_from_pw_multi_aff(pma);
2312 }
2313
2314 /* Given a basic map with a single output dimension that is defined
2315  * in terms of the parameters and input dimensions using an equality,
2316  * extract an isl_aff that expresses the output dimension in terms
2317  * of the parameters and input dimensions.
2318  *
2319  * Since some applications expect the result of isl_pw_multi_aff_from_map
2320  * to only contain integer affine expressions, we compute the floor
2321  * of the expression before returning.
2322  *
2323  * This function shares some similarities with
2324  * isl_basic_map_has_defining_equality and isl_constraint_get_bound.
2325  */
2326 static __isl_give isl_aff *extract_isl_aff_from_basic_map(
2327         __isl_take isl_basic_map *bmap)
2328 {
2329         int i;
2330         unsigned offset;
2331         unsigned total;
2332         isl_local_space *ls;
2333         isl_aff *aff;
2334
2335         if (!bmap)
2336                 return NULL;
2337         if (isl_basic_map_dim(bmap, isl_dim_out) != 1)
2338                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
2339                         "basic map should have a single output dimension",
2340                         goto error);
2341         offset = isl_basic_map_offset(bmap, isl_dim_out);
2342         total = isl_basic_map_total_dim(bmap);
2343         for (i = 0; i < bmap->n_eq; ++i) {
2344                 if (isl_int_is_zero(bmap->eq[i][offset]))
2345                         continue;
2346                 if (isl_seq_first_non_zero(bmap->eq[i] + offset + 1,
2347                                            1 + total - (offset + 1)) != -1)
2348                         continue;
2349                 break;
2350         }
2351         if (i >= bmap->n_eq)
2352                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
2353                         "unable to find suitable equality", goto error);
2354         ls = isl_basic_map_get_local_space(bmap);
2355         aff = isl_aff_alloc(isl_local_space_domain(ls));
2356         if (!aff)
2357                 goto error;
2358         if (isl_int_is_neg(bmap->eq[i][offset]))
2359                 isl_seq_cpy(aff->v->el + 1, bmap->eq[i], offset);
2360         else
2361                 isl_seq_neg(aff->v->el + 1, bmap->eq[i], offset);
2362         isl_seq_clr(aff->v->el + 1 + offset, aff->v->size - (1 + offset));
2363         isl_int_abs(aff->v->el[0], bmap->eq[i][offset]);
2364         isl_basic_map_free(bmap);
2365
2366         aff = isl_aff_remove_unused_divs(aff);
2367         aff = isl_aff_floor(aff);
2368         return aff;
2369 error:
2370         isl_basic_map_free(bmap);
2371         return NULL;
2372 }
2373
2374 /* Given a basic map where each output dimension is defined
2375  * in terms of the parameters and input dimensions using an equality,
2376  * extract an isl_multi_aff that expresses the output dimensions in terms
2377  * of the parameters and input dimensions.
2378  */
2379 static __isl_give isl_multi_aff *extract_isl_multi_aff_from_basic_map(
2380         __isl_take isl_basic_map *bmap)
2381 {
2382         int i;
2383         unsigned n_out;
2384         isl_multi_aff *ma;
2385
2386         if (!bmap)
2387                 return NULL;
2388
2389         ma = isl_multi_aff_alloc(isl_basic_map_get_space(bmap));
2390         n_out = isl_basic_map_dim(bmap, isl_dim_out);
2391
2392         for (i = 0; i < n_out; ++i) {
2393                 isl_basic_map *bmap_i;
2394                 isl_aff *aff;
2395
2396                 bmap_i = isl_basic_map_copy(bmap);
2397                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out,
2398                                                         i + 1, n_out - (1 + i));
2399                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out, 0, i);
2400                 aff = extract_isl_aff_from_basic_map(bmap_i);
2401                 ma = isl_multi_aff_set_aff(ma, i, aff);
2402         }
2403
2404         isl_basic_map_free(bmap);
2405
2406         return ma;
2407 }
2408
2409 /* Create an isl_pw_multi_aff that is equivalent to
2410  * isl_map_intersect_domain(isl_map_from_basic_map(bmap), domain).
2411  * The given basic map is such that each output dimension is defined
2412  * in terms of the parameters and input dimensions using an equality.
2413  */
2414 static __isl_give isl_pw_multi_aff *plain_pw_multi_aff_from_map(
2415         __isl_take isl_set *domain, __isl_take isl_basic_map *bmap)
2416 {
2417         isl_multi_aff *ma;
2418
2419         ma = extract_isl_multi_aff_from_basic_map(bmap);
2420         return isl_pw_multi_aff_alloc(domain, ma);
2421 }
2422
2423 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
2424  * This obivously only works if the input "map" is single-valued.
2425  * If so, we compute the lexicographic minimum of the image in the form
2426  * of an isl_pw_multi_aff.  Since the image is unique, it is equal
2427  * to its lexicographic minimum.
2428  * If the input is not single-valued, we produce an error.
2429  *
2430  * As a special case, we first check if all output dimensions are uniquely
2431  * defined in terms of the parameters and input dimensions over the entire
2432  * domain.  If so, we extract the desired isl_pw_multi_aff directly
2433  * from the affine hull of "map" and its domain.
2434  */
2435 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_map(__isl_take isl_map *map)
2436 {
2437         int i;
2438         int sv;
2439         isl_pw_multi_aff *pma;
2440         isl_basic_map *hull;
2441
2442         if (!map)
2443                 return NULL;
2444
2445         hull = isl_map_affine_hull(isl_map_copy(map));
2446         sv = isl_basic_map_plain_is_single_valued(hull);
2447         if (sv >= 0 && sv)
2448                 return plain_pw_multi_aff_from_map(isl_map_domain(map), hull);
2449         isl_basic_map_free(hull);
2450         if (sv < 0)
2451                 goto error;
2452
2453         sv = isl_map_is_single_valued(map);
2454         if (sv < 0)
2455                 goto error;
2456         if (!sv)
2457                 isl_die(isl_map_get_ctx(map), isl_error_invalid,
2458                         "map is not single-valued", goto error);
2459         map = isl_map_make_disjoint(map);
2460         if (!map)
2461                 return NULL;
2462
2463         pma = isl_pw_multi_aff_empty(isl_map_get_space(map));
2464
2465         for (i = 0; i < map->n; ++i) {
2466                 isl_pw_multi_aff *pma_i;
2467                 isl_basic_map *bmap;
2468                 bmap = isl_basic_map_copy(map->p[i]);
2469                 pma_i = isl_basic_map_lexmin_pw_multi_aff(bmap);
2470                 pma = isl_pw_multi_aff_add_disjoint(pma, pma_i);
2471         }
2472
2473         isl_map_free(map);
2474         return pma;
2475 error:
2476         isl_map_free(map);
2477         return NULL;
2478 }
2479
2480 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_set(__isl_take isl_set *set)
2481 {
2482         return isl_pw_multi_aff_from_map(set);
2483 }
2484
2485 /* Return the piecewise affine expression "set ? 1 : 0".
2486  */
2487 __isl_give isl_pw_aff *isl_set_indicator_function(__isl_take isl_set *set)
2488 {
2489         isl_pw_aff *pa;
2490         isl_space *space = isl_set_get_space(set);
2491         isl_local_space *ls = isl_local_space_from_space(space);
2492         isl_aff *zero = isl_aff_zero_on_domain(isl_local_space_copy(ls));
2493         isl_aff *one = isl_aff_zero_on_domain(ls);
2494
2495         one = isl_aff_add_constant_si(one, 1);
2496         pa = isl_pw_aff_alloc(isl_set_copy(set), one);
2497         set = isl_set_complement(set);
2498         pa = isl_pw_aff_add_disjoint(pa, isl_pw_aff_alloc(set, zero));
2499
2500         return pa;
2501 }
2502
2503 /* Plug in "subs" for dimension "type", "pos" of "aff".
2504  *
2505  * Let i be the dimension to replace and let "subs" be of the form
2506  *
2507  *      f/d
2508  *
2509  * and "aff" of the form
2510  *
2511  *      (a i + g)/m
2512  *
2513  * The result is
2514  *
2515  *      floor((a f + d g')/(m d))
2516  *
2517  * where g' is the result of plugging in "subs" in each of the integer
2518  * divisions in g.
2519  */
2520 __isl_give isl_aff *isl_aff_substitute(__isl_take isl_aff *aff,
2521         enum isl_dim_type type, unsigned pos, __isl_keep isl_aff *subs)
2522 {
2523         isl_ctx *ctx;
2524         isl_int v;
2525
2526         aff = isl_aff_cow(aff);
2527         if (!aff || !subs)
2528                 return isl_aff_free(aff);
2529
2530         ctx = isl_aff_get_ctx(aff);
2531         if (!isl_space_is_equal(aff->ls->dim, subs->ls->dim))
2532                 isl_die(ctx, isl_error_invalid,
2533                         "spaces don't match", return isl_aff_free(aff));
2534         if (isl_local_space_dim(subs->ls, isl_dim_div) != 0)
2535                 isl_die(ctx, isl_error_unsupported,
2536                         "cannot handle divs yet", return isl_aff_free(aff));
2537
2538         aff->ls = isl_local_space_substitute(aff->ls, type, pos, subs);
2539         if (!aff->ls)
2540                 return isl_aff_free(aff);
2541
2542         aff->v = isl_vec_cow(aff->v);
2543         if (!aff->v)
2544                 return isl_aff_free(aff);
2545
2546         pos += isl_local_space_offset(aff->ls, type);
2547
2548         isl_int_init(v);
2549         isl_int_set(v, aff->v->el[1 + pos]);
2550         isl_int_set_si(aff->v->el[1 + pos], 0);
2551         isl_seq_combine(aff->v->el + 1, subs->v->el[0], aff->v->el + 1,
2552                         v, subs->v->el + 1, subs->v->size - 1);
2553         isl_int_mul(aff->v->el[0], aff->v->el[0], subs->v->el[0]);
2554         isl_int_clear(v);
2555
2556         return aff;
2557 }
2558
2559 /* Plug in "subs" for dimension "type", "pos" in each of the affine
2560  * expressions in "maff".
2561  */
2562 __isl_give isl_multi_aff *isl_multi_aff_substitute(
2563         __isl_take isl_multi_aff *maff, enum isl_dim_type type, unsigned pos,
2564         __isl_keep isl_aff *subs)
2565 {
2566         int i;
2567
2568         maff = isl_multi_aff_cow(maff);
2569         if (!maff || !subs)
2570                 return isl_multi_aff_free(maff);
2571
2572         if (type == isl_dim_in)
2573                 type = isl_dim_set;
2574
2575         for (i = 0; i < maff->n; ++i) {
2576                 maff->p[i] = isl_aff_substitute(maff->p[i], type, pos, subs);
2577                 if (!maff->p[i])
2578                         return isl_multi_aff_free(maff);
2579         }
2580
2581         return maff;
2582 }
2583
2584 /* Plug in "subs" for dimension "type", "pos" of "pma".
2585  *
2586  * pma is of the form
2587  *
2588  *      A_i(v) -> M_i(v)
2589  *
2590  * while subs is of the form
2591  *
2592  *      v' = B_j(v) -> S_j
2593  *
2594  * Each pair i,j such that C_ij = A_i \cap B_i is non-empty
2595  * has a contribution in the result, in particular
2596  *
2597  *      C_ij(S_j) -> M_i(S_j)
2598  *
2599  * Note that plugging in S_j in C_ij may also result in an empty set
2600  * and this contribution should simply be discarded.
2601  */
2602 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_substitute(
2603         __isl_take isl_pw_multi_aff *pma, enum isl_dim_type type, unsigned pos,
2604         __isl_keep isl_pw_aff *subs)
2605 {
2606         int i, j, n;
2607         isl_pw_multi_aff *res;
2608
2609         if (!pma || !subs)
2610                 return isl_pw_multi_aff_free(pma);
2611
2612         n = pma->n * subs->n;
2613         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma->dim), n);
2614
2615         for (i = 0; i < pma->n; ++i) {
2616                 for (j = 0; j < subs->n; ++j) {
2617                         isl_set *common;
2618                         isl_multi_aff *res_ij;
2619                         common = isl_set_intersect(
2620                                         isl_set_copy(pma->p[i].set),
2621                                         isl_set_copy(subs->p[j].set));
2622                         common = isl_set_substitute(common,
2623                                         type, pos, subs->p[j].aff);
2624                         if (isl_set_plain_is_empty(common)) {
2625                                 isl_set_free(common);
2626                                 continue;
2627                         }
2628
2629                         res_ij = isl_multi_aff_substitute(
2630                                         isl_multi_aff_copy(pma->p[i].maff),
2631                                         type, pos, subs->p[j].aff);
2632
2633                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
2634                 }
2635         }
2636
2637         isl_pw_multi_aff_free(pma);
2638         return res;
2639 }
2640
2641 /* Extend the local space of "dst" to include the divs
2642  * in the local space of "src".
2643  */
2644 __isl_give isl_aff *isl_aff_align_divs(__isl_take isl_aff *dst,
2645         __isl_keep isl_aff *src)
2646 {
2647         isl_ctx *ctx;
2648         int *exp1 = NULL;
2649         int *exp2 = NULL;
2650         isl_mat *div;
2651
2652         if (!src || !dst)
2653                 return isl_aff_free(dst);
2654
2655         ctx = isl_aff_get_ctx(src);
2656         if (!isl_space_is_equal(src->ls->dim, dst->ls->dim))
2657                 isl_die(ctx, isl_error_invalid,
2658                         "spaces don't match", goto error);
2659
2660         if (src->ls->div->n_row == 0)
2661                 return dst;
2662
2663         exp1 = isl_alloc_array(ctx, int, src->ls->div->n_row);
2664         exp2 = isl_alloc_array(ctx, int, dst->ls->div->n_row);
2665         if (!exp1 || !exp2)
2666                 goto error;
2667
2668         div = isl_merge_divs(src->ls->div, dst->ls->div, exp1, exp2);
2669         dst = isl_aff_expand_divs(dst, div, exp2);
2670         free(exp1);
2671         free(exp2);
2672
2673         return dst;
2674 error:
2675         free(exp1);
2676         free(exp2);
2677         return isl_aff_free(dst);
2678 }
2679
2680 /* Adjust the local spaces of the affine expressions in "maff"
2681  * such that they all have the save divs.
2682  */
2683 __isl_give isl_multi_aff *isl_multi_aff_align_divs(
2684         __isl_take isl_multi_aff *maff)
2685 {
2686         int i;
2687
2688         if (!maff)
2689                 return NULL;
2690         if (maff->n == 0)
2691                 return maff;
2692         maff = isl_multi_aff_cow(maff);
2693         if (!maff)
2694                 return NULL;
2695
2696         for (i = 1; i < maff->n; ++i)
2697                 maff->p[0] = isl_aff_align_divs(maff->p[0], maff->p[i]);
2698         for (i = 1; i < maff->n; ++i) {
2699                 maff->p[i] = isl_aff_align_divs(maff->p[i], maff->p[0]);
2700                 if (!maff->p[i])
2701                         return isl_multi_aff_free(maff);
2702         }
2703
2704         return maff;
2705 }
2706
2707 __isl_give isl_aff *isl_aff_lift(__isl_take isl_aff *aff)
2708 {
2709         aff = isl_aff_cow(aff);
2710         if (!aff)
2711                 return NULL;
2712
2713         aff->ls = isl_local_space_lift(aff->ls);
2714         if (!aff->ls)
2715                 return isl_aff_free(aff);
2716
2717         return aff;
2718 }
2719
2720 /* Lift "maff" to a space with extra dimensions such that the result
2721  * has no more existentially quantified variables.
2722  * If "ls" is not NULL, then *ls is assigned the local space that lies
2723  * at the basis of the lifting applied to "maff".
2724  */
2725 __isl_give isl_multi_aff *isl_multi_aff_lift(__isl_take isl_multi_aff *maff,
2726         __isl_give isl_local_space **ls)
2727 {
2728         int i;
2729         isl_space *space;
2730         unsigned n_div;
2731
2732         if (ls)
2733                 *ls = NULL;
2734
2735         if (!maff)
2736                 return NULL;
2737
2738         if (maff->n == 0) {
2739                 if (ls) {
2740                         isl_space *space = isl_multi_aff_get_domain_space(maff);
2741                         *ls = isl_local_space_from_space(space);
2742                         if (!*ls)
2743                                 return isl_multi_aff_free(maff);
2744                 }
2745                 return maff;
2746         }
2747
2748         maff = isl_multi_aff_cow(maff);
2749         maff = isl_multi_aff_align_divs(maff);
2750         if (!maff)
2751                 return NULL;
2752
2753         n_div = isl_aff_dim(maff->p[0], isl_dim_div);
2754         space = isl_multi_aff_get_space(maff);
2755         space = isl_space_lift(isl_space_domain(space), n_div);
2756         space = isl_space_extend_domain_with_range(space,
2757                                                 isl_multi_aff_get_space(maff));
2758         if (!space)
2759                 return isl_multi_aff_free(maff);
2760         isl_space_free(maff->space);
2761         maff->space = space;
2762
2763         if (ls) {
2764                 *ls = isl_aff_get_domain_local_space(maff->p[0]);
2765                 if (!*ls)
2766                         return isl_multi_aff_free(maff);
2767         }
2768
2769         for (i = 0; i < maff->n; ++i) {
2770                 maff->p[i] = isl_aff_lift(maff->p[i]);
2771                 if (!maff->p[i])
2772                         goto error;
2773         }
2774
2775         return maff;
2776 error:
2777         if (ls)
2778                 isl_local_space_free(*ls);
2779         return isl_multi_aff_free(maff);
2780 }
2781
2782
2783 /* Extract an isl_pw_aff corresponding to output dimension "pos" of "pma".
2784  */
2785 __isl_give isl_pw_aff *isl_pw_multi_aff_get_pw_aff(
2786         __isl_keep isl_pw_multi_aff *pma, int pos)
2787 {
2788         int i;
2789         int n_out;
2790         isl_space *space;
2791         isl_pw_aff *pa;
2792
2793         if (!pma)
2794                 return NULL;
2795
2796         n_out = isl_pw_multi_aff_dim(pma, isl_dim_out);
2797         if (pos < 0 || pos >= n_out)
2798                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
2799                         "index out of bounds", return NULL);
2800
2801         space = isl_pw_multi_aff_get_space(pma);
2802         space = isl_space_drop_dims(space, isl_dim_out,
2803                                     pos + 1, n_out - pos - 1);
2804         space = isl_space_drop_dims(space, isl_dim_out, 0, pos);
2805
2806         pa = isl_pw_aff_alloc_size(space, pma->n);
2807         for (i = 0; i < pma->n; ++i) {
2808                 isl_aff *aff;
2809                 aff = isl_multi_aff_get_aff(pma->p[i].maff, pos);
2810                 pa = isl_pw_aff_add_piece(pa, isl_set_copy(pma->p[i].set), aff);
2811         }
2812
2813         return pa;
2814 }