generalize isl_multi_aff_zero
[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 MIT 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 #define ISL_DIM_H
16 #include <isl_map_private.h>
17 #include <isl_union_map_private.h>
18 #include <isl_aff_private.h>
19 #include <isl_space_private.h>
20 #include <isl_local_space_private.h>
21 #include <isl_mat_private.h>
22 #include <isl_list_private.h>
23 #include <isl/constraint.h>
24 #include <isl/seq.h>
25 #include <isl/set.h>
26 #include <isl_config.h>
27
28 __isl_give isl_aff *isl_aff_alloc_vec(__isl_take isl_local_space *ls,
29         __isl_take isl_vec *v)
30 {
31         isl_aff *aff;
32
33         if (!ls || !v)
34                 goto error;
35
36         aff = isl_calloc_type(v->ctx, struct isl_aff);
37         if (!aff)
38                 goto error;
39
40         aff->ref = 1;
41         aff->ls = ls;
42         aff->v = v;
43
44         return aff;
45 error:
46         isl_local_space_free(ls);
47         isl_vec_free(v);
48         return NULL;
49 }
50
51 __isl_give isl_aff *isl_aff_alloc(__isl_take isl_local_space *ls)
52 {
53         isl_ctx *ctx;
54         isl_vec *v;
55         unsigned total;
56
57         if (!ls)
58                 return NULL;
59
60         ctx = isl_local_space_get_ctx(ls);
61         if (!isl_local_space_divs_known(ls))
62                 isl_die(ctx, isl_error_invalid, "local space has unknown divs",
63                         goto error);
64         if (!isl_local_space_is_set(ls))
65                 isl_die(ctx, isl_error_invalid,
66                         "domain of affine expression should be a set",
67                         goto error);
68
69         total = isl_local_space_dim(ls, isl_dim_all);
70         v = isl_vec_alloc(ctx, 1 + 1 + total);
71         return isl_aff_alloc_vec(ls, v);
72 error:
73         isl_local_space_free(ls);
74         return NULL;
75 }
76
77 __isl_give isl_aff *isl_aff_zero_on_domain(__isl_take isl_local_space *ls)
78 {
79         isl_aff *aff;
80
81         aff = isl_aff_alloc(ls);
82         if (!aff)
83                 return NULL;
84
85         isl_int_set_si(aff->v->el[0], 1);
86         isl_seq_clr(aff->v->el + 1, aff->v->size - 1);
87
88         return aff;
89 }
90
91 /* Return a piecewise affine expression defined on the specified domain
92  * that is equal to zero.
93  */
94 __isl_give isl_pw_aff *isl_pw_aff_zero_on_domain(__isl_take isl_local_space *ls)
95 {
96         return isl_pw_aff_from_aff(isl_aff_zero_on_domain(ls));
97 }
98
99 /* Return an affine expression that is equal to the specified dimension
100  * in "ls".
101  */
102 __isl_give isl_aff *isl_aff_var_on_domain(__isl_take isl_local_space *ls,
103         enum isl_dim_type type, unsigned pos)
104 {
105         isl_space *space;
106         isl_aff *aff;
107
108         if (!ls)
109                 return NULL;
110
111         space = isl_local_space_get_space(ls);
112         if (!space)
113                 goto error;
114         if (isl_space_is_map(space))
115                 isl_die(isl_space_get_ctx(space), isl_error_invalid,
116                         "expecting (parameter) set space", goto error);
117         if (pos >= isl_local_space_dim(ls, type))
118                 isl_die(isl_space_get_ctx(space), isl_error_invalid,
119                         "position out of bounds", goto error);
120
121         isl_space_free(space);
122         aff = isl_aff_alloc(ls);
123         if (!aff)
124                 return NULL;
125
126         pos += isl_local_space_offset(aff->ls, type);
127
128         isl_int_set_si(aff->v->el[0], 1);
129         isl_seq_clr(aff->v->el + 1, aff->v->size - 1);
130         isl_int_set_si(aff->v->el[1 + pos], 1);
131
132         return aff;
133 error:
134         isl_local_space_free(ls);
135         isl_space_free(space);
136         return NULL;
137 }
138
139 /* Return a piecewise affine expression that is equal to
140  * the specified dimension in "ls".
141  */
142 __isl_give isl_pw_aff *isl_pw_aff_var_on_domain(__isl_take isl_local_space *ls,
143         enum isl_dim_type type, unsigned pos)
144 {
145         return isl_pw_aff_from_aff(isl_aff_var_on_domain(ls, type, pos));
146 }
147
148 __isl_give isl_aff *isl_aff_copy(__isl_keep isl_aff *aff)
149 {
150         if (!aff)
151                 return NULL;
152
153         aff->ref++;
154         return aff;
155 }
156
157 __isl_give isl_aff *isl_aff_dup(__isl_keep isl_aff *aff)
158 {
159         if (!aff)
160                 return NULL;
161
162         return isl_aff_alloc_vec(isl_local_space_copy(aff->ls),
163                                  isl_vec_copy(aff->v));
164 }
165
166 __isl_give isl_aff *isl_aff_cow(__isl_take isl_aff *aff)
167 {
168         if (!aff)
169                 return NULL;
170
171         if (aff->ref == 1)
172                 return aff;
173         aff->ref--;
174         return isl_aff_dup(aff);
175 }
176
177 void *isl_aff_free(__isl_take isl_aff *aff)
178 {
179         if (!aff)
180                 return NULL;
181
182         if (--aff->ref > 0)
183                 return NULL;
184
185         isl_local_space_free(aff->ls);
186         isl_vec_free(aff->v);
187
188         free(aff);
189
190         return NULL;
191 }
192
193 isl_ctx *isl_aff_get_ctx(__isl_keep isl_aff *aff)
194 {
195         return aff ? isl_local_space_get_ctx(aff->ls) : NULL;
196 }
197
198 /* Externally, an isl_aff has a map space, but internally, the
199  * ls field corresponds to the domain of that space.
200  */
201 int isl_aff_dim(__isl_keep isl_aff *aff, enum isl_dim_type type)
202 {
203         if (!aff)
204                 return 0;
205         if (type == isl_dim_out)
206                 return 1;
207         if (type == isl_dim_in)
208                 type = isl_dim_set;
209         return isl_local_space_dim(aff->ls, type);
210 }
211
212 __isl_give isl_space *isl_aff_get_domain_space(__isl_keep isl_aff *aff)
213 {
214         return aff ? isl_local_space_get_space(aff->ls) : NULL;
215 }
216
217 __isl_give isl_space *isl_aff_get_space(__isl_keep isl_aff *aff)
218 {
219         isl_space *space;
220         if (!aff)
221                 return NULL;
222         space = isl_local_space_get_space(aff->ls);
223         space = isl_space_from_domain(space);
224         space = isl_space_add_dims(space, isl_dim_out, 1);
225         return space;
226 }
227
228 __isl_give isl_local_space *isl_aff_get_domain_local_space(
229         __isl_keep isl_aff *aff)
230 {
231         return aff ? isl_local_space_copy(aff->ls) : NULL;
232 }
233
234 __isl_give isl_local_space *isl_aff_get_local_space(__isl_keep isl_aff *aff)
235 {
236         isl_local_space *ls;
237         if (!aff)
238                 return NULL;
239         ls = isl_local_space_copy(aff->ls);
240         ls = isl_local_space_from_domain(ls);
241         ls = isl_local_space_add_dims(ls, isl_dim_out, 1);
242         return ls;
243 }
244
245 /* Externally, an isl_aff has a map space, but internally, the
246  * ls field corresponds to the domain of that space.
247  */
248 const char *isl_aff_get_dim_name(__isl_keep isl_aff *aff,
249         enum isl_dim_type type, unsigned pos)
250 {
251         if (!aff)
252                 return NULL;
253         if (type == isl_dim_out)
254                 return NULL;
255         if (type == isl_dim_in)
256                 type = isl_dim_set;
257         return isl_local_space_get_dim_name(aff->ls, type, pos);
258 }
259
260 __isl_give isl_aff *isl_aff_reset_domain_space(__isl_take isl_aff *aff,
261         __isl_take isl_space *dim)
262 {
263         aff = isl_aff_cow(aff);
264         if (!aff || !dim)
265                 goto error;
266
267         aff->ls = isl_local_space_reset_space(aff->ls, dim);
268         if (!aff->ls)
269                 return isl_aff_free(aff);
270
271         return aff;
272 error:
273         isl_aff_free(aff);
274         isl_space_free(dim);
275         return NULL;
276 }
277
278 /* Reset the space of "aff".  This function is called from isl_pw_templ.c
279  * and doesn't know if the space of an element object is represented
280  * directly or through its domain.  It therefore passes along both.
281  */
282 __isl_give isl_aff *isl_aff_reset_space_and_domain(__isl_take isl_aff *aff,
283         __isl_take isl_space *space, __isl_take isl_space *domain)
284 {
285         isl_space_free(space);
286         return isl_aff_reset_domain_space(aff, domain);
287 }
288
289 /* Reorder the coefficients of the affine expression based
290  * on the given reodering.
291  * The reordering r is assumed to have been extended with the local
292  * variables.
293  */
294 static __isl_give isl_vec *vec_reorder(__isl_take isl_vec *vec,
295         __isl_take isl_reordering *r, int n_div)
296 {
297         isl_vec *res;
298         int i;
299
300         if (!vec || !r)
301                 goto error;
302
303         res = isl_vec_alloc(vec->ctx,
304                             2 + isl_space_dim(r->dim, isl_dim_all) + n_div);
305         isl_seq_cpy(res->el, vec->el, 2);
306         isl_seq_clr(res->el + 2, res->size - 2);
307         for (i = 0; i < r->len; ++i)
308                 isl_int_set(res->el[2 + r->pos[i]], vec->el[2 + i]);
309
310         isl_reordering_free(r);
311         isl_vec_free(vec);
312         return res;
313 error:
314         isl_vec_free(vec);
315         isl_reordering_free(r);
316         return NULL;
317 }
318
319 /* Reorder the dimensions of the domain of "aff" according
320  * to the given reordering.
321  */
322 __isl_give isl_aff *isl_aff_realign_domain(__isl_take isl_aff *aff,
323         __isl_take isl_reordering *r)
324 {
325         aff = isl_aff_cow(aff);
326         if (!aff)
327                 goto error;
328
329         r = isl_reordering_extend(r, aff->ls->div->n_row);
330         aff->v = vec_reorder(aff->v, isl_reordering_copy(r),
331                                 aff->ls->div->n_row);
332         aff->ls = isl_local_space_realign(aff->ls, r);
333
334         if (!aff->v || !aff->ls)
335                 return isl_aff_free(aff);
336
337         return aff;
338 error:
339         isl_aff_free(aff);
340         isl_reordering_free(r);
341         return NULL;
342 }
343
344 __isl_give isl_aff *isl_aff_align_params(__isl_take isl_aff *aff,
345         __isl_take isl_space *model)
346 {
347         if (!aff || !model)
348                 goto error;
349
350         if (!isl_space_match(aff->ls->dim, isl_dim_param,
351                              model, isl_dim_param)) {
352                 isl_reordering *exp;
353
354                 model = isl_space_drop_dims(model, isl_dim_in,
355                                         0, isl_space_dim(model, isl_dim_in));
356                 model = isl_space_drop_dims(model, isl_dim_out,
357                                         0, isl_space_dim(model, isl_dim_out));
358                 exp = isl_parameter_alignment_reordering(aff->ls->dim, model);
359                 exp = isl_reordering_extend_space(exp,
360                                         isl_aff_get_domain_space(aff));
361                 aff = isl_aff_realign_domain(aff, exp);
362         }
363
364         isl_space_free(model);
365         return aff;
366 error:
367         isl_space_free(model);
368         isl_aff_free(aff);
369         return NULL;
370 }
371
372 int isl_aff_plain_is_zero(__isl_keep isl_aff *aff)
373 {
374         if (!aff)
375                 return -1;
376
377         return isl_seq_first_non_zero(aff->v->el + 1, aff->v->size - 1) < 0;
378 }
379
380 int isl_aff_plain_is_equal(__isl_keep isl_aff *aff1, __isl_keep isl_aff *aff2)
381 {
382         int equal;
383
384         if (!aff1 || !aff2)
385                 return -1;
386
387         equal = isl_local_space_is_equal(aff1->ls, aff2->ls);
388         if (equal < 0 || !equal)
389                 return equal;
390
391         return isl_vec_is_equal(aff1->v, aff2->v);
392 }
393
394 int isl_aff_get_denominator(__isl_keep isl_aff *aff, isl_int *v)
395 {
396         if (!aff)
397                 return -1;
398         isl_int_set(*v, aff->v->el[0]);
399         return 0;
400 }
401
402 int isl_aff_get_constant(__isl_keep isl_aff *aff, isl_int *v)
403 {
404         if (!aff)
405                 return -1;
406         isl_int_set(*v, aff->v->el[1]);
407         return 0;
408 }
409
410 int isl_aff_get_coefficient(__isl_keep isl_aff *aff,
411         enum isl_dim_type type, int pos, isl_int *v)
412 {
413         if (!aff)
414                 return -1;
415
416         if (type == isl_dim_out)
417                 isl_die(aff->v->ctx, isl_error_invalid,
418                         "output/set dimension does not have a coefficient",
419                         return -1);
420         if (type == isl_dim_in)
421                 type = isl_dim_set;
422
423         if (pos >= isl_local_space_dim(aff->ls, type))
424                 isl_die(aff->v->ctx, isl_error_invalid,
425                         "position out of bounds", return -1);
426
427         pos += isl_local_space_offset(aff->ls, type);
428         isl_int_set(*v, aff->v->el[1 + pos]);
429
430         return 0;
431 }
432
433 __isl_give isl_aff *isl_aff_set_denominator(__isl_take isl_aff *aff, isl_int v)
434 {
435         aff = isl_aff_cow(aff);
436         if (!aff)
437                 return NULL;
438
439         aff->v = isl_vec_cow(aff->v);
440         if (!aff->v)
441                 return isl_aff_free(aff);
442
443         isl_int_set(aff->v->el[0], v);
444
445         return aff;
446 }
447
448 __isl_give isl_aff *isl_aff_set_constant(__isl_take isl_aff *aff, isl_int v)
449 {
450         aff = isl_aff_cow(aff);
451         if (!aff)
452                 return NULL;
453
454         aff->v = isl_vec_cow(aff->v);
455         if (!aff->v)
456                 return isl_aff_free(aff);
457
458         isl_int_set(aff->v->el[1], v);
459
460         return aff;
461 }
462
463 __isl_give isl_aff *isl_aff_add_constant(__isl_take isl_aff *aff, isl_int v)
464 {
465         if (isl_int_is_zero(v))
466                 return aff;
467
468         aff = isl_aff_cow(aff);
469         if (!aff)
470                 return NULL;
471
472         aff->v = isl_vec_cow(aff->v);
473         if (!aff->v)
474                 return isl_aff_free(aff);
475
476         isl_int_addmul(aff->v->el[1], aff->v->el[0], v);
477
478         return aff;
479 }
480
481 __isl_give isl_aff *isl_aff_add_constant_si(__isl_take isl_aff *aff, int v)
482 {
483         isl_int t;
484
485         isl_int_init(t);
486         isl_int_set_si(t, v);
487         aff = isl_aff_add_constant(aff, t);
488         isl_int_clear(t);
489
490         return aff;
491 }
492
493 /* Add "v" to the numerator of the constant term of "aff".
494  */
495 __isl_give isl_aff *isl_aff_add_constant_num(__isl_take isl_aff *aff, isl_int v)
496 {
497         if (isl_int_is_zero(v))
498                 return aff;
499
500         aff = isl_aff_cow(aff);
501         if (!aff)
502                 return NULL;
503
504         aff->v = isl_vec_cow(aff->v);
505         if (!aff->v)
506                 return isl_aff_free(aff);
507
508         isl_int_add(aff->v->el[1], aff->v->el[1], v);
509
510         return aff;
511 }
512
513 /* Add "v" to the numerator of the constant term of "aff".
514  */
515 __isl_give isl_aff *isl_aff_add_constant_num_si(__isl_take isl_aff *aff, int v)
516 {
517         isl_int t;
518
519         if (v == 0)
520                 return aff;
521
522         isl_int_init(t);
523         isl_int_set_si(t, v);
524         aff = isl_aff_add_constant_num(aff, t);
525         isl_int_clear(t);
526
527         return aff;
528 }
529
530 __isl_give isl_aff *isl_aff_set_constant_si(__isl_take isl_aff *aff, int v)
531 {
532         aff = isl_aff_cow(aff);
533         if (!aff)
534                 return NULL;
535
536         aff->v = isl_vec_cow(aff->v);
537         if (!aff->v)
538                 return isl_aff_free(aff);
539
540         isl_int_set_si(aff->v->el[1], v);
541
542         return aff;
543 }
544
545 __isl_give isl_aff *isl_aff_set_coefficient(__isl_take isl_aff *aff,
546         enum isl_dim_type type, int pos, isl_int v)
547 {
548         if (!aff)
549                 return NULL;
550
551         if (type == isl_dim_out)
552                 isl_die(aff->v->ctx, isl_error_invalid,
553                         "output/set dimension does not have a coefficient",
554                         return isl_aff_free(aff));
555         if (type == isl_dim_in)
556                 type = isl_dim_set;
557
558         if (pos >= isl_local_space_dim(aff->ls, type))
559                 isl_die(aff->v->ctx, isl_error_invalid,
560                         "position out of bounds", return isl_aff_free(aff));
561
562         aff = isl_aff_cow(aff);
563         if (!aff)
564                 return NULL;
565
566         aff->v = isl_vec_cow(aff->v);
567         if (!aff->v)
568                 return isl_aff_free(aff);
569
570         pos += isl_local_space_offset(aff->ls, type);
571         isl_int_set(aff->v->el[1 + pos], v);
572
573         return aff;
574 }
575
576 __isl_give isl_aff *isl_aff_set_coefficient_si(__isl_take isl_aff *aff,
577         enum isl_dim_type type, int pos, int v)
578 {
579         if (!aff)
580                 return NULL;
581
582         if (type == isl_dim_out)
583                 isl_die(aff->v->ctx, isl_error_invalid,
584                         "output/set dimension does not have a coefficient",
585                         return isl_aff_free(aff));
586         if (type == isl_dim_in)
587                 type = isl_dim_set;
588
589         if (pos >= isl_local_space_dim(aff->ls, type))
590                 isl_die(aff->v->ctx, isl_error_invalid,
591                         "position out of bounds", return isl_aff_free(aff));
592
593         aff = isl_aff_cow(aff);
594         if (!aff)
595                 return NULL;
596
597         aff->v = isl_vec_cow(aff->v);
598         if (!aff->v)
599                 return isl_aff_free(aff);
600
601         pos += isl_local_space_offset(aff->ls, type);
602         isl_int_set_si(aff->v->el[1 + pos], v);
603
604         return aff;
605 }
606
607 __isl_give isl_aff *isl_aff_add_coefficient(__isl_take isl_aff *aff,
608         enum isl_dim_type type, int pos, isl_int v)
609 {
610         if (!aff)
611                 return NULL;
612
613         if (type == isl_dim_out)
614                 isl_die(aff->v->ctx, isl_error_invalid,
615                         "output/set dimension does not have a coefficient",
616                         return isl_aff_free(aff));
617         if (type == isl_dim_in)
618                 type = isl_dim_set;
619
620         if (pos >= isl_local_space_dim(aff->ls, type))
621                 isl_die(aff->v->ctx, isl_error_invalid,
622                         "position out of bounds", return isl_aff_free(aff));
623
624         aff = isl_aff_cow(aff);
625         if (!aff)
626                 return NULL;
627
628         aff->v = isl_vec_cow(aff->v);
629         if (!aff->v)
630                 return isl_aff_free(aff);
631
632         pos += isl_local_space_offset(aff->ls, type);
633         isl_int_addmul(aff->v->el[1 + pos], aff->v->el[0], v);
634
635         return aff;
636 }
637
638 __isl_give isl_aff *isl_aff_add_coefficient_si(__isl_take isl_aff *aff,
639         enum isl_dim_type type, int pos, int v)
640 {
641         isl_int t;
642
643         isl_int_init(t);
644         isl_int_set_si(t, v);
645         aff = isl_aff_add_coefficient(aff, type, pos, t);
646         isl_int_clear(t);
647
648         return aff;
649 }
650
651 __isl_give isl_aff *isl_aff_get_div(__isl_keep isl_aff *aff, int pos)
652 {
653         if (!aff)
654                 return NULL;
655
656         return isl_local_space_get_div(aff->ls, pos);
657 }
658
659 __isl_give isl_aff *isl_aff_neg(__isl_take isl_aff *aff)
660 {
661         aff = isl_aff_cow(aff);
662         if (!aff)
663                 return NULL;
664         aff->v = isl_vec_cow(aff->v);
665         if (!aff->v)
666                 return isl_aff_free(aff);
667
668         isl_seq_neg(aff->v->el + 1, aff->v->el + 1, aff->v->size - 1);
669
670         return aff;
671 }
672
673 /* Remove divs from the local space that do not appear in the affine
674  * expression.
675  * We currently only remove divs at the end.
676  * Some intermediate divs may also not appear directly in the affine
677  * expression, but we would also need to check that no other divs are
678  * defined in terms of them.
679  */
680 __isl_give isl_aff *isl_aff_remove_unused_divs( __isl_take isl_aff *aff)
681 {
682         int pos;
683         int off;
684         int n;
685
686         if (!aff)
687                 return NULL;
688
689         n = isl_local_space_dim(aff->ls, isl_dim_div);
690         off = isl_local_space_offset(aff->ls, isl_dim_div);
691
692         pos = isl_seq_last_non_zero(aff->v->el + 1 + off, n) + 1;
693         if (pos == n)
694                 return aff;
695
696         aff = isl_aff_cow(aff);
697         if (!aff)
698                 return NULL;
699
700         aff->ls = isl_local_space_drop_dims(aff->ls, isl_dim_div, pos, n - pos);
701         aff->v = isl_vec_drop_els(aff->v, 1 + off + pos, n - pos);
702         if (!aff->ls || !aff->v)
703                 return isl_aff_free(aff);
704
705         return aff;
706 }
707
708 /* Given two affine expressions "p" of length p_len (including the
709  * denominator and the constant term) and "subs" of length subs_len,
710  * plug in "subs" for the variable at position "pos".
711  * The variables of "subs" and "p" are assumed to match up to subs_len,
712  * but "p" may have additional variables.
713  * "v" is an initialized isl_int that can be used internally.
714  *
715  * In particular, if "p" represents the expression
716  *
717  *      (a i + g)/m
718  *
719  * with i the variable at position "pos" and "subs" represents the expression
720  *
721  *      f/d
722  *
723  * then the result represents the expression
724  *
725  *      (a f + d g)/(m d)
726  *
727  */
728 void isl_seq_substitute(isl_int *p, int pos, isl_int *subs,
729         int p_len, int subs_len, isl_int v)
730 {
731         isl_int_set(v, p[1 + pos]);
732         isl_int_set_si(p[1 + pos], 0);
733         isl_seq_combine(p + 1, subs[0], p + 1, v, subs + 1, subs_len - 1);
734         isl_seq_scale(p + subs_len, p + subs_len, subs[0], p_len - subs_len);
735         isl_int_mul(p[0], p[0], subs[0]);
736 }
737
738 /* Look for any divs in the aff->ls with a denominator equal to one
739  * and plug them into the affine expression and any subsequent divs
740  * that may reference the div.
741  */
742 static __isl_give isl_aff *plug_in_integral_divs(__isl_take isl_aff *aff)
743 {
744         int i, n;
745         int len;
746         isl_int v;
747         isl_vec *vec;
748         isl_local_space *ls;
749         unsigned pos;
750
751         if (!aff)
752                 return NULL;
753
754         n = isl_local_space_dim(aff->ls, isl_dim_div);
755         len = aff->v->size;
756         for (i = 0; i < n; ++i) {
757                 if (!isl_int_is_one(aff->ls->div->row[i][0]))
758                         continue;
759                 ls = isl_local_space_copy(aff->ls);
760                 ls = isl_local_space_substitute_seq(ls, isl_dim_div, i,
761                                             aff->ls->div->row[i], len, i + 1);
762                 vec = isl_vec_copy(aff->v);
763                 vec = isl_vec_cow(vec);
764                 if (!ls || !vec)
765                         goto error;
766
767                 isl_int_init(v);
768
769                 pos = isl_local_space_offset(aff->ls, isl_dim_div) + i;
770                 isl_seq_substitute(vec->el, pos, aff->ls->div->row[i],
771                                         len, len, v);
772
773                 isl_int_clear(v);
774
775                 isl_vec_free(aff->v);
776                 aff->v = vec;
777                 isl_local_space_free(aff->ls);
778                 aff->ls = ls;
779         }
780
781         return aff;
782 error:
783         isl_vec_free(vec);
784         isl_local_space_free(ls);
785         return isl_aff_free(aff);
786 }
787
788 /* Swap divs "a" and "b" in "aff", which is assumed to be non-NULL.
789  *
790  * Even though this function is only called on isl_affs with a single
791  * reference, we are careful to only change aff->v and aff->ls together.
792  */
793 static __isl_give isl_aff *swap_div(__isl_take isl_aff *aff, int a, int b)
794 {
795         unsigned off = isl_local_space_offset(aff->ls, isl_dim_div);
796         isl_local_space *ls;
797         isl_vec *v;
798
799         ls = isl_local_space_copy(aff->ls);
800         ls = isl_local_space_swap_div(ls, a, b);
801         v = isl_vec_copy(aff->v);
802         v = isl_vec_cow(v);
803         if (!ls || !v)
804                 goto error;
805
806         isl_int_swap(v->el[1 + off + a], v->el[1 + off + b]);
807         isl_vec_free(aff->v);
808         aff->v = v;
809         isl_local_space_free(aff->ls);
810         aff->ls = ls;
811
812         return aff;
813 error:
814         isl_vec_free(v);
815         isl_local_space_free(ls);
816         return isl_aff_free(aff);
817 }
818
819 /* Merge divs "a" and "b" in "aff", which is assumed to be non-NULL.
820  *
821  * We currently do not actually remove div "b", but simply add its
822  * coefficient to that of "a" and then zero it out.
823  */
824 static __isl_give isl_aff *merge_divs(__isl_take isl_aff *aff, int a, int b)
825 {
826         unsigned off = isl_local_space_offset(aff->ls, isl_dim_div);
827
828         if (isl_int_is_zero(aff->v->el[1 + off + b]))
829                 return aff;
830
831         aff->v = isl_vec_cow(aff->v);
832         if (!aff->v)
833                 return isl_aff_free(aff);
834
835         isl_int_add(aff->v->el[1 + off + a],
836                     aff->v->el[1 + off + a], aff->v->el[1 + off + b]);
837         isl_int_set_si(aff->v->el[1 + off + b], 0);
838
839         return aff;
840 }
841
842 /* Sort the divs in the local space of "aff" according to
843  * the comparison function "cmp_row" in isl_local_space.c,
844  * combining the coefficients of identical divs.
845  *
846  * Reordering divs does not change the semantics of "aff",
847  * so there is no need to call isl_aff_cow.
848  * Moreover, this function is currently only called on isl_affs
849  * with a single reference.
850  */
851 static __isl_give isl_aff *sort_divs(__isl_take isl_aff *aff)
852 {
853         int i, j, n;
854         unsigned off;
855
856         if (!aff)
857                 return NULL;
858
859         off = isl_local_space_offset(aff->ls, isl_dim_div);
860         n = isl_aff_dim(aff, isl_dim_div);
861         for (i = 1; i < n; ++i) {
862                 for (j = i - 1; j >= 0; --j) {
863                         int cmp = isl_mat_cmp_div(aff->ls->div, j, j + 1);
864                         if (cmp < 0)
865                                 break;
866                         if (cmp == 0)
867                                 aff = merge_divs(aff, j, j + 1);
868                         else
869                                 aff = swap_div(aff, j, j + 1);
870                         if (!aff)
871                                 return NULL;
872                 }
873         }
874
875         return aff;
876 }
877
878 /* Normalize the representation of "aff".
879  *
880  * This function should only be called of "new" isl_affs, i.e.,
881  * with only a single reference.  We therefore do not need to
882  * worry about affecting other instances.
883  */
884 __isl_give isl_aff *isl_aff_normalize(__isl_take isl_aff *aff)
885 {
886         if (!aff)
887                 return NULL;
888         aff->v = isl_vec_normalize(aff->v);
889         if (!aff->v)
890                 return isl_aff_free(aff);
891         aff = plug_in_integral_divs(aff);
892         aff = sort_divs(aff);
893         aff = isl_aff_remove_unused_divs(aff);
894         return aff;
895 }
896
897 /* Given f, return floor(f).
898  * If f is an integer expression, then just return f.
899  * If f is a constant, then return the constant floor(f).
900  * Otherwise, if f = g/m, write g = q m + r,
901  * create a new div d = [r/m] and return the expression q + d.
902  * The coefficients in r are taken to lie between -m/2 and m/2.
903  */
904 __isl_give isl_aff *isl_aff_floor(__isl_take isl_aff *aff)
905 {
906         int i;
907         int size;
908         isl_ctx *ctx;
909         isl_vec *div;
910
911         if (!aff)
912                 return NULL;
913
914         if (isl_int_is_one(aff->v->el[0]))
915                 return aff;
916
917         aff = isl_aff_cow(aff);
918         if (!aff)
919                 return NULL;
920
921         aff->v = isl_vec_cow(aff->v);
922         if (!aff->v)
923                 return isl_aff_free(aff);
924
925         if (isl_aff_is_cst(aff)) {
926                 isl_int_fdiv_q(aff->v->el[1], aff->v->el[1], aff->v->el[0]);
927                 isl_int_set_si(aff->v->el[0], 1);
928                 return aff;
929         }
930
931         div = isl_vec_copy(aff->v);
932         div = isl_vec_cow(div);
933         if (!div)
934                 return isl_aff_free(aff);
935
936         ctx = isl_aff_get_ctx(aff);
937         isl_int_fdiv_q(aff->v->el[0], aff->v->el[0], ctx->two);
938         for (i = 1; i < aff->v->size; ++i) {
939                 isl_int_fdiv_r(div->el[i], div->el[i], div->el[0]);
940                 isl_int_fdiv_q(aff->v->el[i], aff->v->el[i], div->el[0]);
941                 if (isl_int_gt(div->el[i], aff->v->el[0])) {
942                         isl_int_sub(div->el[i], div->el[i], div->el[0]);
943                         isl_int_add_ui(aff->v->el[i], aff->v->el[i], 1);
944                 }
945         }
946
947         aff->ls = isl_local_space_add_div(aff->ls, div);
948         if (!aff->ls)
949                 return isl_aff_free(aff);
950
951         size = aff->v->size;
952         aff->v = isl_vec_extend(aff->v, size + 1);
953         if (!aff->v)
954                 return isl_aff_free(aff);
955         isl_int_set_si(aff->v->el[0], 1);
956         isl_int_set_si(aff->v->el[size], 1);
957
958         return aff;
959 }
960
961 /* Compute
962  *
963  *      aff mod m = aff - m * floor(aff/m)
964  */
965 __isl_give isl_aff *isl_aff_mod(__isl_take isl_aff *aff, isl_int m)
966 {
967         isl_aff *res;
968
969         res = isl_aff_copy(aff);
970         aff = isl_aff_scale_down(aff, m);
971         aff = isl_aff_floor(aff);
972         aff = isl_aff_scale(aff, m);
973         res = isl_aff_sub(res, aff);
974
975         return res;
976 }
977
978 /* Compute
979  *
980  *      pwaff mod m = pwaff - m * floor(pwaff/m)
981  */
982 __isl_give isl_pw_aff *isl_pw_aff_mod(__isl_take isl_pw_aff *pwaff, isl_int m)
983 {
984         isl_pw_aff *res;
985
986         res = isl_pw_aff_copy(pwaff);
987         pwaff = isl_pw_aff_scale_down(pwaff, m);
988         pwaff = isl_pw_aff_floor(pwaff);
989         pwaff = isl_pw_aff_scale(pwaff, m);
990         res = isl_pw_aff_sub(res, pwaff);
991
992         return res;
993 }
994
995 /* Given f, return ceil(f).
996  * If f is an integer expression, then just return f.
997  * Otherwise, create a new div d = [-f] and return the expression -d.
998  */
999 __isl_give isl_aff *isl_aff_ceil(__isl_take isl_aff *aff)
1000 {
1001         if (!aff)
1002                 return NULL;
1003
1004         if (isl_int_is_one(aff->v->el[0]))
1005                 return aff;
1006
1007         aff = isl_aff_neg(aff);
1008         aff = isl_aff_floor(aff);
1009         aff = isl_aff_neg(aff);
1010
1011         return aff;
1012 }
1013
1014 /* Apply the expansion computed by isl_merge_divs.
1015  * The expansion itself is given by "exp" while the resulting
1016  * list of divs is given by "div".
1017  */
1018 __isl_give isl_aff *isl_aff_expand_divs( __isl_take isl_aff *aff,
1019         __isl_take isl_mat *div, int *exp)
1020 {
1021         int i, j;
1022         int old_n_div;
1023         int new_n_div;
1024         int offset;
1025
1026         aff = isl_aff_cow(aff);
1027         if (!aff || !div)
1028                 goto error;
1029
1030         old_n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1031         new_n_div = isl_mat_rows(div);
1032         if (new_n_div < old_n_div)
1033                 isl_die(isl_mat_get_ctx(div), isl_error_invalid,
1034                         "not an expansion", goto error);
1035
1036         aff->v = isl_vec_extend(aff->v, aff->v->size + new_n_div - old_n_div);
1037         if (!aff->v)
1038                 goto error;
1039
1040         offset = 1 + isl_local_space_offset(aff->ls, isl_dim_div);
1041         j = old_n_div - 1;
1042         for (i = new_n_div - 1; i >= 0; --i) {
1043                 if (j >= 0 && exp[j] == i) {
1044                         if (i != j)
1045                                 isl_int_swap(aff->v->el[offset + i],
1046                                              aff->v->el[offset + j]);
1047                         j--;
1048                 } else
1049                         isl_int_set_si(aff->v->el[offset + i], 0);
1050         }
1051
1052         aff->ls = isl_local_space_replace_divs(aff->ls, isl_mat_copy(div));
1053         if (!aff->ls)
1054                 goto error;
1055         isl_mat_free(div);
1056         return aff;
1057 error:
1058         isl_aff_free(aff);
1059         isl_mat_free(div);
1060         return NULL;
1061 }
1062
1063 /* Add two affine expressions that live in the same local space.
1064  */
1065 static __isl_give isl_aff *add_expanded(__isl_take isl_aff *aff1,
1066         __isl_take isl_aff *aff2)
1067 {
1068         isl_int gcd, f;
1069
1070         aff1 = isl_aff_cow(aff1);
1071         if (!aff1 || !aff2)
1072                 goto error;
1073
1074         aff1->v = isl_vec_cow(aff1->v);
1075         if (!aff1->v)
1076                 goto error;
1077
1078         isl_int_init(gcd);
1079         isl_int_init(f);
1080         isl_int_gcd(gcd, aff1->v->el[0], aff2->v->el[0]);
1081         isl_int_divexact(f, aff2->v->el[0], gcd);
1082         isl_seq_scale(aff1->v->el + 1, aff1->v->el + 1, f, aff1->v->size - 1);
1083         isl_int_divexact(f, aff1->v->el[0], gcd);
1084         isl_seq_addmul(aff1->v->el + 1, f, aff2->v->el + 1, aff1->v->size - 1);
1085         isl_int_divexact(f, aff2->v->el[0], gcd);
1086         isl_int_mul(aff1->v->el[0], aff1->v->el[0], f);
1087         isl_int_clear(f);
1088         isl_int_clear(gcd);
1089
1090         isl_aff_free(aff2);
1091         return aff1;
1092 error:
1093         isl_aff_free(aff1);
1094         isl_aff_free(aff2);
1095         return NULL;
1096 }
1097
1098 __isl_give isl_aff *isl_aff_add(__isl_take isl_aff *aff1,
1099         __isl_take isl_aff *aff2)
1100 {
1101         isl_ctx *ctx;
1102         int *exp1 = NULL;
1103         int *exp2 = NULL;
1104         isl_mat *div;
1105
1106         if (!aff1 || !aff2)
1107                 goto error;
1108
1109         ctx = isl_aff_get_ctx(aff1);
1110         if (!isl_space_is_equal(aff1->ls->dim, aff2->ls->dim))
1111                 isl_die(ctx, isl_error_invalid,
1112                         "spaces don't match", goto error);
1113
1114         if (aff1->ls->div->n_row == 0 && aff2->ls->div->n_row == 0)
1115                 return add_expanded(aff1, aff2);
1116
1117         exp1 = isl_alloc_array(ctx, int, aff1->ls->div->n_row);
1118         exp2 = isl_alloc_array(ctx, int, aff2->ls->div->n_row);
1119         if (!exp1 || !exp2)
1120                 goto error;
1121
1122         div = isl_merge_divs(aff1->ls->div, aff2->ls->div, exp1, exp2);
1123         aff1 = isl_aff_expand_divs(aff1, isl_mat_copy(div), exp1);
1124         aff2 = isl_aff_expand_divs(aff2, div, exp2);
1125         free(exp1);
1126         free(exp2);
1127
1128         return add_expanded(aff1, aff2);
1129 error:
1130         free(exp1);
1131         free(exp2);
1132         isl_aff_free(aff1);
1133         isl_aff_free(aff2);
1134         return NULL;
1135 }
1136
1137 __isl_give isl_aff *isl_aff_sub(__isl_take isl_aff *aff1,
1138         __isl_take isl_aff *aff2)
1139 {
1140         return isl_aff_add(aff1, isl_aff_neg(aff2));
1141 }
1142
1143 __isl_give isl_aff *isl_aff_scale(__isl_take isl_aff *aff, isl_int f)
1144 {
1145         isl_int gcd;
1146
1147         if (isl_int_is_one(f))
1148                 return aff;
1149
1150         aff = isl_aff_cow(aff);
1151         if (!aff)
1152                 return NULL;
1153         aff->v = isl_vec_cow(aff->v);
1154         if (!aff->v)
1155                 return isl_aff_free(aff);
1156
1157         isl_int_init(gcd);
1158         isl_int_gcd(gcd, aff->v->el[0], f);
1159         isl_int_divexact(aff->v->el[0], aff->v->el[0], gcd);
1160         isl_int_divexact(gcd, f, gcd);
1161         isl_seq_scale(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
1162         isl_int_clear(gcd);
1163
1164         return aff;
1165 }
1166
1167 __isl_give isl_aff *isl_aff_scale_down(__isl_take isl_aff *aff, isl_int f)
1168 {
1169         isl_int gcd;
1170
1171         if (isl_int_is_one(f))
1172                 return aff;
1173
1174         aff = isl_aff_cow(aff);
1175         if (!aff)
1176                 return NULL;
1177
1178         if (isl_int_is_zero(f))
1179                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1180                         "cannot scale down by zero", return isl_aff_free(aff));
1181
1182         aff->v = isl_vec_cow(aff->v);
1183         if (!aff->v)
1184                 return isl_aff_free(aff);
1185
1186         isl_int_init(gcd);
1187         isl_seq_gcd(aff->v->el + 1, aff->v->size - 1, &gcd);
1188         isl_int_gcd(gcd, gcd, f);
1189         isl_seq_scale_down(aff->v->el + 1, aff->v->el + 1, gcd, aff->v->size - 1);
1190         isl_int_divexact(gcd, f, gcd);
1191         isl_int_mul(aff->v->el[0], aff->v->el[0], gcd);
1192         isl_int_clear(gcd);
1193
1194         return aff;
1195 }
1196
1197 __isl_give isl_aff *isl_aff_scale_down_ui(__isl_take isl_aff *aff, unsigned f)
1198 {
1199         isl_int v;
1200
1201         if (f == 1)
1202                 return aff;
1203
1204         isl_int_init(v);
1205         isl_int_set_ui(v, f);
1206         aff = isl_aff_scale_down(aff, v);
1207         isl_int_clear(v);
1208
1209         return aff;
1210 }
1211
1212 __isl_give isl_aff *isl_aff_set_dim_name(__isl_take isl_aff *aff,
1213         enum isl_dim_type type, unsigned pos, const char *s)
1214 {
1215         aff = isl_aff_cow(aff);
1216         if (!aff)
1217                 return NULL;
1218         if (type == isl_dim_out)
1219                 isl_die(aff->v->ctx, isl_error_invalid,
1220                         "cannot set name of output/set dimension",
1221                         return isl_aff_free(aff));
1222         if (type == isl_dim_in)
1223                 type = isl_dim_set;
1224         aff->ls = isl_local_space_set_dim_name(aff->ls, type, pos, s);
1225         if (!aff->ls)
1226                 return isl_aff_free(aff);
1227
1228         return aff;
1229 }
1230
1231 __isl_give isl_aff *isl_aff_set_dim_id(__isl_take isl_aff *aff,
1232         enum isl_dim_type type, unsigned pos, __isl_take isl_id *id)
1233 {
1234         aff = isl_aff_cow(aff);
1235         if (!aff)
1236                 return isl_id_free(id);
1237         if (type == isl_dim_out)
1238                 isl_die(aff->v->ctx, isl_error_invalid,
1239                         "cannot set name of output/set dimension",
1240                         goto error);
1241         if (type == isl_dim_in)
1242                 type = isl_dim_set;
1243         aff->ls = isl_local_space_set_dim_id(aff->ls, type, pos, id);
1244         if (!aff->ls)
1245                 return isl_aff_free(aff);
1246
1247         return aff;
1248 error:
1249         isl_id_free(id);
1250         isl_aff_free(aff);
1251         return NULL;
1252 }
1253
1254 /* Exploit the equalities in "eq" to simplify the affine expression
1255  * and the expressions of the integer divisions in the local space.
1256  * The integer divisions in this local space are assumed to appear
1257  * as regular dimensions in "eq".
1258  */
1259 static __isl_give isl_aff *isl_aff_substitute_equalities_lifted(
1260         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1261 {
1262         int i, j;
1263         unsigned total;
1264         unsigned n_div;
1265
1266         if (!eq)
1267                 goto error;
1268         if (eq->n_eq == 0) {
1269                 isl_basic_set_free(eq);
1270                 return aff;
1271         }
1272
1273         aff = isl_aff_cow(aff);
1274         if (!aff)
1275                 goto error;
1276
1277         aff->ls = isl_local_space_substitute_equalities(aff->ls,
1278                                                         isl_basic_set_copy(eq));
1279         if (!aff->ls)
1280                 goto error;
1281
1282         total = 1 + isl_space_dim(eq->dim, isl_dim_all);
1283         n_div = eq->n_div;
1284         for (i = 0; i < eq->n_eq; ++i) {
1285                 j = isl_seq_last_non_zero(eq->eq[i], total + n_div);
1286                 if (j < 0 || j == 0 || j >= total)
1287                         continue;
1288
1289                 isl_seq_elim(aff->v->el + 1, eq->eq[i], j, total,
1290                                 &aff->v->el[0]);
1291         }
1292
1293         isl_basic_set_free(eq);
1294         aff = isl_aff_normalize(aff);
1295         return aff;
1296 error:
1297         isl_basic_set_free(eq);
1298         isl_aff_free(aff);
1299         return NULL;
1300 }
1301
1302 /* Exploit the equalities in "eq" to simplify the affine expression
1303  * and the expressions of the integer divisions in the local space.
1304  */
1305 static __isl_give isl_aff *isl_aff_substitute_equalities(
1306         __isl_take isl_aff *aff, __isl_take isl_basic_set *eq)
1307 {
1308         int n_div;
1309
1310         if (!aff || !eq)
1311                 goto error;
1312         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1313         if (n_div > 0)
1314                 eq = isl_basic_set_add(eq, isl_dim_set, n_div);
1315         return isl_aff_substitute_equalities_lifted(aff, eq);
1316 error:
1317         isl_basic_set_free(eq);
1318         isl_aff_free(aff);
1319         return NULL;
1320 }
1321
1322 /* Look for equalities among the variables shared by context and aff
1323  * and the integer divisions of aff, if any.
1324  * The equalities are then used to eliminate coefficients and/or integer
1325  * divisions from aff.
1326  */
1327 __isl_give isl_aff *isl_aff_gist(__isl_take isl_aff *aff,
1328         __isl_take isl_set *context)
1329 {
1330         isl_basic_set *hull;
1331         int n_div;
1332
1333         if (!aff)
1334                 goto error;
1335         n_div = isl_local_space_dim(aff->ls, isl_dim_div);
1336         if (n_div > 0) {
1337                 isl_basic_set *bset;
1338                 isl_local_space *ls;
1339                 context = isl_set_add_dims(context, isl_dim_set, n_div);
1340                 ls = isl_aff_get_domain_local_space(aff);
1341                 bset = isl_basic_set_from_local_space(ls);
1342                 bset = isl_basic_set_lift(bset);
1343                 bset = isl_basic_set_flatten(bset);
1344                 context = isl_set_intersect(context,
1345                                             isl_set_from_basic_set(bset));
1346         }
1347
1348         hull = isl_set_affine_hull(context);
1349         return isl_aff_substitute_equalities_lifted(aff, hull);
1350 error:
1351         isl_aff_free(aff);
1352         isl_set_free(context);
1353         return NULL;
1354 }
1355
1356 __isl_give isl_aff *isl_aff_gist_params(__isl_take isl_aff *aff,
1357         __isl_take isl_set *context)
1358 {
1359         isl_set *dom_context = isl_set_universe(isl_aff_get_domain_space(aff));
1360         dom_context = isl_set_intersect_params(dom_context, context);
1361         return isl_aff_gist(aff, dom_context);
1362 }
1363
1364 /* Return a basic set containing those elements in the space
1365  * of aff where it is non-negative.
1366  */
1367 __isl_give isl_basic_set *isl_aff_nonneg_basic_set(__isl_take isl_aff *aff)
1368 {
1369         isl_constraint *ineq;
1370         isl_basic_set *bset;
1371
1372         ineq = isl_inequality_from_aff(aff);
1373
1374         bset = isl_basic_set_from_constraint(ineq);
1375         bset = isl_basic_set_simplify(bset);
1376         return bset;
1377 }
1378
1379 /* Return a basic set containing those elements in the domain space
1380  * of aff where it is negative.
1381  */
1382 __isl_give isl_basic_set *isl_aff_neg_basic_set(__isl_take isl_aff *aff)
1383 {
1384         aff = isl_aff_neg(aff);
1385         aff = isl_aff_add_constant_num_si(aff, -1);
1386         return isl_aff_nonneg_basic_set(aff);
1387 }
1388
1389 /* Return a basic set containing those elements in the space
1390  * of aff where it is zero.
1391  */
1392 __isl_give isl_basic_set *isl_aff_zero_basic_set(__isl_take isl_aff *aff)
1393 {
1394         isl_constraint *ineq;
1395         isl_basic_set *bset;
1396
1397         ineq = isl_equality_from_aff(aff);
1398
1399         bset = isl_basic_set_from_constraint(ineq);
1400         bset = isl_basic_set_simplify(bset);
1401         return bset;
1402 }
1403
1404 /* Return a basic set containing those elements in the shared space
1405  * of aff1 and aff2 where aff1 is greater than or equal to aff2.
1406  */
1407 __isl_give isl_basic_set *isl_aff_ge_basic_set(__isl_take isl_aff *aff1,
1408         __isl_take isl_aff *aff2)
1409 {
1410         aff1 = isl_aff_sub(aff1, aff2);
1411
1412         return isl_aff_nonneg_basic_set(aff1);
1413 }
1414
1415 /* Return a basic set containing those elements in the shared space
1416  * of aff1 and aff2 where aff1 is smaller than or equal to aff2.
1417  */
1418 __isl_give isl_basic_set *isl_aff_le_basic_set(__isl_take isl_aff *aff1,
1419         __isl_take isl_aff *aff2)
1420 {
1421         return isl_aff_ge_basic_set(aff2, aff1);
1422 }
1423
1424 __isl_give isl_aff *isl_aff_add_on_domain(__isl_keep isl_set *dom,
1425         __isl_take isl_aff *aff1, __isl_take isl_aff *aff2)
1426 {
1427         aff1 = isl_aff_add(aff1, aff2);
1428         aff1 = isl_aff_gist(aff1, isl_set_copy(dom));
1429         return aff1;
1430 }
1431
1432 int isl_aff_is_empty(__isl_keep isl_aff *aff)
1433 {
1434         if (!aff)
1435                 return -1;
1436
1437         return 0;
1438 }
1439
1440 /* Check whether the given affine expression has non-zero coefficient
1441  * for any dimension in the given range or if any of these dimensions
1442  * appear with non-zero coefficients in any of the integer divisions
1443  * involved in the affine expression.
1444  */
1445 int isl_aff_involves_dims(__isl_keep isl_aff *aff,
1446         enum isl_dim_type type, unsigned first, unsigned n)
1447 {
1448         int i;
1449         isl_ctx *ctx;
1450         int *active = NULL;
1451         int involves = 0;
1452
1453         if (!aff)
1454                 return -1;
1455         if (n == 0)
1456                 return 0;
1457
1458         ctx = isl_aff_get_ctx(aff);
1459         if (first + n > isl_aff_dim(aff, type))
1460                 isl_die(ctx, isl_error_invalid,
1461                         "range out of bounds", return -1);
1462
1463         active = isl_local_space_get_active(aff->ls, aff->v->el + 2);
1464         if (!active)
1465                 goto error;
1466
1467         first += isl_local_space_offset(aff->ls, type) - 1;
1468         for (i = 0; i < n; ++i)
1469                 if (active[first + i]) {
1470                         involves = 1;
1471                         break;
1472                 }
1473
1474         free(active);
1475
1476         return involves;
1477 error:
1478         free(active);
1479         return -1;
1480 }
1481
1482 __isl_give isl_aff *isl_aff_drop_dims(__isl_take isl_aff *aff,
1483         enum isl_dim_type type, unsigned first, unsigned n)
1484 {
1485         isl_ctx *ctx;
1486
1487         if (!aff)
1488                 return NULL;
1489         if (type == isl_dim_out)
1490                 isl_die(aff->v->ctx, isl_error_invalid,
1491                         "cannot drop output/set dimension",
1492                         return isl_aff_free(aff));
1493         if (type == isl_dim_in)
1494                 type = isl_dim_set;
1495         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1496                 return aff;
1497
1498         ctx = isl_aff_get_ctx(aff);
1499         if (first + n > isl_local_space_dim(aff->ls, type))
1500                 isl_die(ctx, isl_error_invalid, "range out of bounds",
1501                         return isl_aff_free(aff));
1502
1503         aff = isl_aff_cow(aff);
1504         if (!aff)
1505                 return NULL;
1506
1507         aff->ls = isl_local_space_drop_dims(aff->ls, type, first, n);
1508         if (!aff->ls)
1509                 return isl_aff_free(aff);
1510
1511         first += 1 + isl_local_space_offset(aff->ls, type);
1512         aff->v = isl_vec_drop_els(aff->v, first, n);
1513         if (!aff->v)
1514                 return isl_aff_free(aff);
1515
1516         return aff;
1517 }
1518
1519 /* Project the domain of the affine expression onto its parameter space.
1520  * The affine expression may not involve any of the domain dimensions.
1521  */
1522 __isl_give isl_aff *isl_aff_project_domain_on_params(__isl_take isl_aff *aff)
1523 {
1524         isl_space *space;
1525         unsigned n;
1526         int involves;
1527
1528         n = isl_aff_dim(aff, isl_dim_in);
1529         involves = isl_aff_involves_dims(aff, isl_dim_in, 0, n);
1530         if (involves < 0)
1531                 return isl_aff_free(aff);
1532         if (involves)
1533                 isl_die(isl_aff_get_ctx(aff), isl_error_invalid,
1534                     "affine expression involves some of the domain dimensions",
1535                     return isl_aff_free(aff));
1536         aff = isl_aff_drop_dims(aff, isl_dim_in, 0, n);
1537         space = isl_aff_get_domain_space(aff);
1538         space = isl_space_params(space);
1539         aff = isl_aff_reset_domain_space(aff, space);
1540         return aff;
1541 }
1542
1543 __isl_give isl_aff *isl_aff_insert_dims(__isl_take isl_aff *aff,
1544         enum isl_dim_type type, unsigned first, unsigned n)
1545 {
1546         isl_ctx *ctx;
1547
1548         if (!aff)
1549                 return NULL;
1550         if (type == isl_dim_out)
1551                 isl_die(aff->v->ctx, isl_error_invalid,
1552                         "cannot insert output/set dimensions",
1553                         return isl_aff_free(aff));
1554         if (type == isl_dim_in)
1555                 type = isl_dim_set;
1556         if (n == 0 && !isl_local_space_is_named_or_nested(aff->ls, type))
1557                 return aff;
1558
1559         ctx = isl_aff_get_ctx(aff);
1560         if (first > isl_local_space_dim(aff->ls, type))
1561                 isl_die(ctx, isl_error_invalid, "position out of bounds",
1562                         return isl_aff_free(aff));
1563
1564         aff = isl_aff_cow(aff);
1565         if (!aff)
1566                 return NULL;
1567
1568         aff->ls = isl_local_space_insert_dims(aff->ls, type, first, n);
1569         if (!aff->ls)
1570                 return isl_aff_free(aff);
1571
1572         first += 1 + isl_local_space_offset(aff->ls, type);
1573         aff->v = isl_vec_insert_zero_els(aff->v, first, n);
1574         if (!aff->v)
1575                 return isl_aff_free(aff);
1576
1577         return aff;
1578 }
1579
1580 __isl_give isl_aff *isl_aff_add_dims(__isl_take isl_aff *aff,
1581         enum isl_dim_type type, unsigned n)
1582 {
1583         unsigned pos;
1584
1585         pos = isl_aff_dim(aff, type);
1586
1587         return isl_aff_insert_dims(aff, type, pos, n);
1588 }
1589
1590 __isl_give isl_pw_aff *isl_pw_aff_add_dims(__isl_take isl_pw_aff *pwaff,
1591         enum isl_dim_type type, unsigned n)
1592 {
1593         unsigned pos;
1594
1595         pos = isl_pw_aff_dim(pwaff, type);
1596
1597         return isl_pw_aff_insert_dims(pwaff, type, pos, n);
1598 }
1599
1600 __isl_give isl_pw_aff *isl_pw_aff_from_aff(__isl_take isl_aff *aff)
1601 {
1602         isl_set *dom = isl_set_universe(isl_aff_get_domain_space(aff));
1603         return isl_pw_aff_alloc(dom, aff);
1604 }
1605
1606 #undef PW
1607 #define PW isl_pw_aff
1608 #undef EL
1609 #define EL isl_aff
1610 #undef EL_IS_ZERO
1611 #define EL_IS_ZERO is_empty
1612 #undef ZERO
1613 #define ZERO empty
1614 #undef IS_ZERO
1615 #define IS_ZERO is_empty
1616 #undef FIELD
1617 #define FIELD aff
1618 #undef DEFAULT_IS_ZERO
1619 #define DEFAULT_IS_ZERO 0
1620
1621 #define NO_EVAL
1622 #define NO_OPT
1623 #define NO_MOVE_DIMS
1624 #define NO_LIFT
1625 #define NO_MORPH
1626
1627 #include <isl_pw_templ.c>
1628
1629 static __isl_give isl_set *align_params_pw_pw_set_and(
1630         __isl_take isl_pw_aff *pwaff1, __isl_take isl_pw_aff *pwaff2,
1631         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
1632                                     __isl_take isl_pw_aff *pwaff2))
1633 {
1634         if (!pwaff1 || !pwaff2)
1635                 goto error;
1636         if (isl_space_match(pwaff1->dim, isl_dim_param,
1637                           pwaff2->dim, isl_dim_param))
1638                 return fn(pwaff1, pwaff2);
1639         if (!isl_space_has_named_params(pwaff1->dim) ||
1640             !isl_space_has_named_params(pwaff2->dim))
1641                 isl_die(isl_pw_aff_get_ctx(pwaff1), isl_error_invalid,
1642                         "unaligned unnamed parameters", goto error);
1643         pwaff1 = isl_pw_aff_align_params(pwaff1, isl_pw_aff_get_space(pwaff2));
1644         pwaff2 = isl_pw_aff_align_params(pwaff2, isl_pw_aff_get_space(pwaff1));
1645         return fn(pwaff1, pwaff2);
1646 error:
1647         isl_pw_aff_free(pwaff1);
1648         isl_pw_aff_free(pwaff2);
1649         return NULL;
1650 }
1651
1652 /* Compute a piecewise quasi-affine expression with a domain that
1653  * is the union of those of pwaff1 and pwaff2 and such that on each
1654  * cell, the quasi-affine expression is the better (according to cmp)
1655  * of those of pwaff1 and pwaff2.  If only one of pwaff1 or pwaff2
1656  * is defined on a given cell, then the associated expression
1657  * is the defined one.
1658  */
1659 static __isl_give isl_pw_aff *pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1660         __isl_take isl_pw_aff *pwaff2,
1661         __isl_give isl_basic_set *(*cmp)(__isl_take isl_aff *aff1,
1662                                         __isl_take isl_aff *aff2))
1663 {
1664         int i, j, n;
1665         isl_pw_aff *res;
1666         isl_ctx *ctx;
1667         isl_set *set;
1668
1669         if (!pwaff1 || !pwaff2)
1670                 goto error;
1671
1672         ctx = isl_space_get_ctx(pwaff1->dim);
1673         if (!isl_space_is_equal(pwaff1->dim, pwaff2->dim))
1674                 isl_die(ctx, isl_error_invalid,
1675                         "arguments should live in same space", goto error);
1676
1677         if (isl_pw_aff_is_empty(pwaff1)) {
1678                 isl_pw_aff_free(pwaff1);
1679                 return pwaff2;
1680         }
1681
1682         if (isl_pw_aff_is_empty(pwaff2)) {
1683                 isl_pw_aff_free(pwaff2);
1684                 return pwaff1;
1685         }
1686
1687         n = 2 * (pwaff1->n + 1) * (pwaff2->n + 1);
1688         res = isl_pw_aff_alloc_size(isl_space_copy(pwaff1->dim), n);
1689
1690         for (i = 0; i < pwaff1->n; ++i) {
1691                 set = isl_set_copy(pwaff1->p[i].set);
1692                 for (j = 0; j < pwaff2->n; ++j) {
1693                         struct isl_set *common;
1694                         isl_set *better;
1695
1696                         common = isl_set_intersect(
1697                                         isl_set_copy(pwaff1->p[i].set),
1698                                         isl_set_copy(pwaff2->p[j].set));
1699                         better = isl_set_from_basic_set(cmp(
1700                                         isl_aff_copy(pwaff2->p[j].aff),
1701                                         isl_aff_copy(pwaff1->p[i].aff)));
1702                         better = isl_set_intersect(common, better);
1703                         if (isl_set_plain_is_empty(better)) {
1704                                 isl_set_free(better);
1705                                 continue;
1706                         }
1707                         set = isl_set_subtract(set, isl_set_copy(better));
1708
1709                         res = isl_pw_aff_add_piece(res, better,
1710                                                 isl_aff_copy(pwaff2->p[j].aff));
1711                 }
1712                 res = isl_pw_aff_add_piece(res, set,
1713                                                 isl_aff_copy(pwaff1->p[i].aff));
1714         }
1715
1716         for (j = 0; j < pwaff2->n; ++j) {
1717                 set = isl_set_copy(pwaff2->p[j].set);
1718                 for (i = 0; i < pwaff1->n; ++i)
1719                         set = isl_set_subtract(set,
1720                                         isl_set_copy(pwaff1->p[i].set));
1721                 res = isl_pw_aff_add_piece(res, set,
1722                                                 isl_aff_copy(pwaff2->p[j].aff));
1723         }
1724
1725         isl_pw_aff_free(pwaff1);
1726         isl_pw_aff_free(pwaff2);
1727
1728         return res;
1729 error:
1730         isl_pw_aff_free(pwaff1);
1731         isl_pw_aff_free(pwaff2);
1732         return NULL;
1733 }
1734
1735 /* Compute a piecewise quasi-affine expression with a domain that
1736  * is the union of those of pwaff1 and pwaff2 and such that on each
1737  * cell, the quasi-affine expression is the maximum of those of pwaff1
1738  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1739  * cell, then the associated expression is the defined one.
1740  */
1741 static __isl_give isl_pw_aff *pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1742         __isl_take isl_pw_aff *pwaff2)
1743 {
1744         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_ge_basic_set);
1745 }
1746
1747 __isl_give isl_pw_aff *isl_pw_aff_union_max(__isl_take isl_pw_aff *pwaff1,
1748         __isl_take isl_pw_aff *pwaff2)
1749 {
1750         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1751                                                         &pw_aff_union_max);
1752 }
1753
1754 /* Compute a piecewise quasi-affine expression with a domain that
1755  * is the union of those of pwaff1 and pwaff2 and such that on each
1756  * cell, the quasi-affine expression is the minimum of those of pwaff1
1757  * and pwaff2.  If only one of pwaff1 or pwaff2 is defined on a given
1758  * cell, then the associated expression is the defined one.
1759  */
1760 static __isl_give isl_pw_aff *pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1761         __isl_take isl_pw_aff *pwaff2)
1762 {
1763         return pw_aff_union_opt(pwaff1, pwaff2, &isl_aff_le_basic_set);
1764 }
1765
1766 __isl_give isl_pw_aff *isl_pw_aff_union_min(__isl_take isl_pw_aff *pwaff1,
1767         __isl_take isl_pw_aff *pwaff2)
1768 {
1769         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2,
1770                                                         &pw_aff_union_min);
1771 }
1772
1773 __isl_give isl_pw_aff *isl_pw_aff_union_opt(__isl_take isl_pw_aff *pwaff1,
1774         __isl_take isl_pw_aff *pwaff2, int max)
1775 {
1776         if (max)
1777                 return isl_pw_aff_union_max(pwaff1, pwaff2);
1778         else
1779                 return isl_pw_aff_union_min(pwaff1, pwaff2);
1780 }
1781
1782 /* Construct a map with as domain the domain of pwaff and
1783  * one-dimensional range corresponding to the affine expressions.
1784  */
1785 static __isl_give isl_map *map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1786 {
1787         int i;
1788         isl_space *dim;
1789         isl_map *map;
1790
1791         if (!pwaff)
1792                 return NULL;
1793
1794         dim = isl_pw_aff_get_space(pwaff);
1795         map = isl_map_empty(dim);
1796
1797         for (i = 0; i < pwaff->n; ++i) {
1798                 isl_basic_map *bmap;
1799                 isl_map *map_i;
1800
1801                 bmap = isl_basic_map_from_aff(isl_aff_copy(pwaff->p[i].aff));
1802                 map_i = isl_map_from_basic_map(bmap);
1803                 map_i = isl_map_intersect_domain(map_i,
1804                                                 isl_set_copy(pwaff->p[i].set));
1805                 map = isl_map_union_disjoint(map, map_i);
1806         }
1807
1808         isl_pw_aff_free(pwaff);
1809
1810         return map;
1811 }
1812
1813 /* Construct a map with as domain the domain of pwaff and
1814  * one-dimensional range corresponding to the affine expressions.
1815  */
1816 __isl_give isl_map *isl_map_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1817 {
1818         if (!pwaff)
1819                 return NULL;
1820         if (isl_space_is_set(pwaff->dim))
1821                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1822                         "space of input is not a map",
1823                         return isl_pw_aff_free(pwaff));
1824         return map_from_pw_aff(pwaff);
1825 }
1826
1827 /* Construct a one-dimensional set with as parameter domain
1828  * the domain of pwaff and the single set dimension
1829  * corresponding to the affine expressions.
1830  */
1831 __isl_give isl_set *isl_set_from_pw_aff(__isl_take isl_pw_aff *pwaff)
1832 {
1833         if (!pwaff)
1834                 return NULL;
1835         if (!isl_space_is_set(pwaff->dim))
1836                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
1837                         "space of input is not a set",
1838                         return isl_pw_aff_free(pwaff));
1839         return map_from_pw_aff(pwaff);
1840 }
1841
1842 /* Return a set containing those elements in the domain
1843  * of pwaff where it is non-negative.
1844  */
1845 __isl_give isl_set *isl_pw_aff_nonneg_set(__isl_take isl_pw_aff *pwaff)
1846 {
1847         int i;
1848         isl_set *set;
1849
1850         if (!pwaff)
1851                 return NULL;
1852
1853         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1854
1855         for (i = 0; i < pwaff->n; ++i) {
1856                 isl_basic_set *bset;
1857                 isl_set *set_i;
1858
1859                 bset = isl_aff_nonneg_basic_set(isl_aff_copy(pwaff->p[i].aff));
1860                 set_i = isl_set_from_basic_set(bset);
1861                 set_i = isl_set_intersect(set_i, isl_set_copy(pwaff->p[i].set));
1862                 set = isl_set_union_disjoint(set, set_i);
1863         }
1864
1865         isl_pw_aff_free(pwaff);
1866
1867         return set;
1868 }
1869
1870 /* Return a set containing those elements in the domain
1871  * of pwaff where it is zero (if complement is 0) or not zero
1872  * (if complement is 1).
1873  */
1874 static __isl_give isl_set *pw_aff_zero_set(__isl_take isl_pw_aff *pwaff,
1875         int complement)
1876 {
1877         int i;
1878         isl_set *set;
1879
1880         if (!pwaff)
1881                 return NULL;
1882
1883         set = isl_set_empty(isl_pw_aff_get_domain_space(pwaff));
1884
1885         for (i = 0; i < pwaff->n; ++i) {
1886                 isl_basic_set *bset;
1887                 isl_set *set_i, *zero;
1888
1889                 bset = isl_aff_zero_basic_set(isl_aff_copy(pwaff->p[i].aff));
1890                 zero = isl_set_from_basic_set(bset);
1891                 set_i = isl_set_copy(pwaff->p[i].set);
1892                 if (complement)
1893                         set_i = isl_set_subtract(set_i, zero);
1894                 else
1895                         set_i = isl_set_intersect(set_i, zero);
1896                 set = isl_set_union_disjoint(set, set_i);
1897         }
1898
1899         isl_pw_aff_free(pwaff);
1900
1901         return set;
1902 }
1903
1904 /* Return a set containing those elements in the domain
1905  * of pwaff where it is zero.
1906  */
1907 __isl_give isl_set *isl_pw_aff_zero_set(__isl_take isl_pw_aff *pwaff)
1908 {
1909         return pw_aff_zero_set(pwaff, 0);
1910 }
1911
1912 /* Return a set containing those elements in the domain
1913  * of pwaff where it is not zero.
1914  */
1915 __isl_give isl_set *isl_pw_aff_non_zero_set(__isl_take isl_pw_aff *pwaff)
1916 {
1917         return pw_aff_zero_set(pwaff, 1);
1918 }
1919
1920 /* Return a set containing those elements in the shared domain
1921  * of pwaff1 and pwaff2 where pwaff1 is greater than (or equal) to pwaff2.
1922  *
1923  * We compute the difference on the shared domain and then construct
1924  * the set of values where this difference is non-negative.
1925  * If strict is set, we first subtract 1 from the difference.
1926  * If equal is set, we only return the elements where pwaff1 and pwaff2
1927  * are equal.
1928  */
1929 static __isl_give isl_set *pw_aff_gte_set(__isl_take isl_pw_aff *pwaff1,
1930         __isl_take isl_pw_aff *pwaff2, int strict, int equal)
1931 {
1932         isl_set *set1, *set2;
1933
1934         set1 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff1));
1935         set2 = isl_pw_aff_domain(isl_pw_aff_copy(pwaff2));
1936         set1 = isl_set_intersect(set1, set2);
1937         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, isl_set_copy(set1));
1938         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, isl_set_copy(set1));
1939         pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_neg(pwaff2));
1940
1941         if (strict) {
1942                 isl_space *dim = isl_set_get_space(set1);
1943                 isl_aff *aff;
1944                 aff = isl_aff_zero_on_domain(isl_local_space_from_space(dim));
1945                 aff = isl_aff_add_constant_si(aff, -1);
1946                 pwaff1 = isl_pw_aff_add(pwaff1, isl_pw_aff_alloc(set1, aff));
1947         } else
1948                 isl_set_free(set1);
1949
1950         if (equal)
1951                 return isl_pw_aff_zero_set(pwaff1);
1952         return isl_pw_aff_nonneg_set(pwaff1);
1953 }
1954
1955 /* Return a set containing those elements in the shared domain
1956  * of pwaff1 and pwaff2 where pwaff1 is equal to pwaff2.
1957  */
1958 static __isl_give isl_set *pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1959         __isl_take isl_pw_aff *pwaff2)
1960 {
1961         return pw_aff_gte_set(pwaff1, pwaff2, 0, 1);
1962 }
1963
1964 __isl_give isl_set *isl_pw_aff_eq_set(__isl_take isl_pw_aff *pwaff1,
1965         __isl_take isl_pw_aff *pwaff2)
1966 {
1967         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_eq_set);
1968 }
1969
1970 /* Return a set containing those elements in the shared domain
1971  * of pwaff1 and pwaff2 where pwaff1 is greater than or equal to pwaff2.
1972  */
1973 static __isl_give isl_set *pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1974         __isl_take isl_pw_aff *pwaff2)
1975 {
1976         return pw_aff_gte_set(pwaff1, pwaff2, 0, 0);
1977 }
1978
1979 __isl_give isl_set *isl_pw_aff_ge_set(__isl_take isl_pw_aff *pwaff1,
1980         __isl_take isl_pw_aff *pwaff2)
1981 {
1982         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ge_set);
1983 }
1984
1985 /* Return a set containing those elements in the shared domain
1986  * of pwaff1 and pwaff2 where pwaff1 is strictly greater than pwaff2.
1987  */
1988 static __isl_give isl_set *pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1989         __isl_take isl_pw_aff *pwaff2)
1990 {
1991         return pw_aff_gte_set(pwaff1, pwaff2, 1, 0);
1992 }
1993
1994 __isl_give isl_set *isl_pw_aff_gt_set(__isl_take isl_pw_aff *pwaff1,
1995         __isl_take isl_pw_aff *pwaff2)
1996 {
1997         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_gt_set);
1998 }
1999
2000 __isl_give isl_set *isl_pw_aff_le_set(__isl_take isl_pw_aff *pwaff1,
2001         __isl_take isl_pw_aff *pwaff2)
2002 {
2003         return isl_pw_aff_ge_set(pwaff2, pwaff1);
2004 }
2005
2006 __isl_give isl_set *isl_pw_aff_lt_set(__isl_take isl_pw_aff *pwaff1,
2007         __isl_take isl_pw_aff *pwaff2)
2008 {
2009         return isl_pw_aff_gt_set(pwaff2, pwaff1);
2010 }
2011
2012 /* Return a set containing those elements in the shared domain
2013  * of the elements of list1 and list2 where each element in list1
2014  * has the relation specified by "fn" with each element in list2.
2015  */
2016 static __isl_give isl_set *pw_aff_list_set(__isl_take isl_pw_aff_list *list1,
2017         __isl_take isl_pw_aff_list *list2,
2018         __isl_give isl_set *(*fn)(__isl_take isl_pw_aff *pwaff1,
2019                                     __isl_take isl_pw_aff *pwaff2))
2020 {
2021         int i, j;
2022         isl_ctx *ctx;
2023         isl_set *set;
2024
2025         if (!list1 || !list2)
2026                 goto error;
2027
2028         ctx = isl_pw_aff_list_get_ctx(list1);
2029         if (list1->n < 1 || list2->n < 1)
2030                 isl_die(ctx, isl_error_invalid,
2031                         "list should contain at least one element", goto error);
2032
2033         set = isl_set_universe(isl_pw_aff_get_domain_space(list1->p[0]));
2034         for (i = 0; i < list1->n; ++i)
2035                 for (j = 0; j < list2->n; ++j) {
2036                         isl_set *set_ij;
2037
2038                         set_ij = fn(isl_pw_aff_copy(list1->p[i]),
2039                                     isl_pw_aff_copy(list2->p[j]));
2040                         set = isl_set_intersect(set, set_ij);
2041                 }
2042
2043         isl_pw_aff_list_free(list1);
2044         isl_pw_aff_list_free(list2);
2045         return set;
2046 error:
2047         isl_pw_aff_list_free(list1);
2048         isl_pw_aff_list_free(list2);
2049         return NULL;
2050 }
2051
2052 /* Return a set containing those elements in the shared domain
2053  * of the elements of list1 and list2 where each element in list1
2054  * is equal to each element in list2.
2055  */
2056 __isl_give isl_set *isl_pw_aff_list_eq_set(__isl_take isl_pw_aff_list *list1,
2057         __isl_take isl_pw_aff_list *list2)
2058 {
2059         return pw_aff_list_set(list1, list2, &isl_pw_aff_eq_set);
2060 }
2061
2062 __isl_give isl_set *isl_pw_aff_list_ne_set(__isl_take isl_pw_aff_list *list1,
2063         __isl_take isl_pw_aff_list *list2)
2064 {
2065         return pw_aff_list_set(list1, list2, &isl_pw_aff_ne_set);
2066 }
2067
2068 /* Return a set containing those elements in the shared domain
2069  * of the elements of list1 and list2 where each element in list1
2070  * is less than or equal to each element in list2.
2071  */
2072 __isl_give isl_set *isl_pw_aff_list_le_set(__isl_take isl_pw_aff_list *list1,
2073         __isl_take isl_pw_aff_list *list2)
2074 {
2075         return pw_aff_list_set(list1, list2, &isl_pw_aff_le_set);
2076 }
2077
2078 __isl_give isl_set *isl_pw_aff_list_lt_set(__isl_take isl_pw_aff_list *list1,
2079         __isl_take isl_pw_aff_list *list2)
2080 {
2081         return pw_aff_list_set(list1, list2, &isl_pw_aff_lt_set);
2082 }
2083
2084 __isl_give isl_set *isl_pw_aff_list_ge_set(__isl_take isl_pw_aff_list *list1,
2085         __isl_take isl_pw_aff_list *list2)
2086 {
2087         return pw_aff_list_set(list1, list2, &isl_pw_aff_ge_set);
2088 }
2089
2090 __isl_give isl_set *isl_pw_aff_list_gt_set(__isl_take isl_pw_aff_list *list1,
2091         __isl_take isl_pw_aff_list *list2)
2092 {
2093         return pw_aff_list_set(list1, list2, &isl_pw_aff_gt_set);
2094 }
2095
2096
2097 /* Return a set containing those elements in the shared domain
2098  * of pwaff1 and pwaff2 where pwaff1 is not equal to pwaff2.
2099  */
2100 static __isl_give isl_set *pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2101         __isl_take isl_pw_aff *pwaff2)
2102 {
2103         isl_set *set_lt, *set_gt;
2104
2105         set_lt = isl_pw_aff_lt_set(isl_pw_aff_copy(pwaff1),
2106                                    isl_pw_aff_copy(pwaff2));
2107         set_gt = isl_pw_aff_gt_set(pwaff1, pwaff2);
2108         return isl_set_union_disjoint(set_lt, set_gt);
2109 }
2110
2111 __isl_give isl_set *isl_pw_aff_ne_set(__isl_take isl_pw_aff *pwaff1,
2112         __isl_take isl_pw_aff *pwaff2)
2113 {
2114         return align_params_pw_pw_set_and(pwaff1, pwaff2, &pw_aff_ne_set);
2115 }
2116
2117 __isl_give isl_pw_aff *isl_pw_aff_scale_down(__isl_take isl_pw_aff *pwaff,
2118         isl_int v)
2119 {
2120         int i;
2121
2122         if (isl_int_is_one(v))
2123                 return pwaff;
2124         if (!isl_int_is_pos(v))
2125                 isl_die(isl_pw_aff_get_ctx(pwaff), isl_error_invalid,
2126                         "factor needs to be positive",
2127                         return isl_pw_aff_free(pwaff));
2128         pwaff = isl_pw_aff_cow(pwaff);
2129         if (!pwaff)
2130                 return NULL;
2131         if (pwaff->n == 0)
2132                 return pwaff;
2133
2134         for (i = 0; i < pwaff->n; ++i) {
2135                 pwaff->p[i].aff = isl_aff_scale_down(pwaff->p[i].aff, v);
2136                 if (!pwaff->p[i].aff)
2137                         return isl_pw_aff_free(pwaff);
2138         }
2139
2140         return pwaff;
2141 }
2142
2143 __isl_give isl_pw_aff *isl_pw_aff_floor(__isl_take isl_pw_aff *pwaff)
2144 {
2145         int i;
2146
2147         pwaff = isl_pw_aff_cow(pwaff);
2148         if (!pwaff)
2149                 return NULL;
2150         if (pwaff->n == 0)
2151                 return pwaff;
2152
2153         for (i = 0; i < pwaff->n; ++i) {
2154                 pwaff->p[i].aff = isl_aff_floor(pwaff->p[i].aff);
2155                 if (!pwaff->p[i].aff)
2156                         return isl_pw_aff_free(pwaff);
2157         }
2158
2159         return pwaff;
2160 }
2161
2162 __isl_give isl_pw_aff *isl_pw_aff_ceil(__isl_take isl_pw_aff *pwaff)
2163 {
2164         int i;
2165
2166         pwaff = isl_pw_aff_cow(pwaff);
2167         if (!pwaff)
2168                 return NULL;
2169         if (pwaff->n == 0)
2170                 return pwaff;
2171
2172         for (i = 0; i < pwaff->n; ++i) {
2173                 pwaff->p[i].aff = isl_aff_ceil(pwaff->p[i].aff);
2174                 if (!pwaff->p[i].aff)
2175                         return isl_pw_aff_free(pwaff);
2176         }
2177
2178         return pwaff;
2179 }
2180
2181 /* Assuming that "cond1" and "cond2" are disjoint,
2182  * return an affine expression that is equal to pwaff1 on cond1
2183  * and to pwaff2 on cond2.
2184  */
2185 static __isl_give isl_pw_aff *isl_pw_aff_select(
2186         __isl_take isl_set *cond1, __isl_take isl_pw_aff *pwaff1,
2187         __isl_take isl_set *cond2, __isl_take isl_pw_aff *pwaff2)
2188 {
2189         pwaff1 = isl_pw_aff_intersect_domain(pwaff1, cond1);
2190         pwaff2 = isl_pw_aff_intersect_domain(pwaff2, cond2);
2191
2192         return isl_pw_aff_add_disjoint(pwaff1, pwaff2);
2193 }
2194
2195 /* Return an affine expression that is equal to pwaff_true for elements
2196  * where "cond" is non-zero and to pwaff_false for elements where "cond"
2197  * is zero.
2198  * That is, return cond ? pwaff_true : pwaff_false;
2199  */
2200 __isl_give isl_pw_aff *isl_pw_aff_cond(__isl_take isl_pw_aff *cond,
2201         __isl_take isl_pw_aff *pwaff_true, __isl_take isl_pw_aff *pwaff_false)
2202 {
2203         isl_set *cond_true, *cond_false;
2204
2205         cond_true = isl_pw_aff_non_zero_set(isl_pw_aff_copy(cond));
2206         cond_false = isl_pw_aff_zero_set(cond);
2207         return isl_pw_aff_select(cond_true, pwaff_true,
2208                                  cond_false, pwaff_false);
2209 }
2210
2211 int isl_aff_is_cst(__isl_keep isl_aff *aff)
2212 {
2213         if (!aff)
2214                 return -1;
2215
2216         return isl_seq_first_non_zero(aff->v->el + 2, aff->v->size - 2) == -1;
2217 }
2218
2219 /* Check whether pwaff is a piecewise constant.
2220  */
2221 int isl_pw_aff_is_cst(__isl_keep isl_pw_aff *pwaff)
2222 {
2223         int i;
2224
2225         if (!pwaff)
2226                 return -1;
2227
2228         for (i = 0; i < pwaff->n; ++i) {
2229                 int is_cst = isl_aff_is_cst(pwaff->p[i].aff);
2230                 if (is_cst < 0 || !is_cst)
2231                         return is_cst;
2232         }
2233
2234         return 1;
2235 }
2236
2237 __isl_give isl_aff *isl_aff_mul(__isl_take isl_aff *aff1,
2238         __isl_take isl_aff *aff2)
2239 {
2240         if (!isl_aff_is_cst(aff2) && isl_aff_is_cst(aff1))
2241                 return isl_aff_mul(aff2, aff1);
2242
2243         if (!isl_aff_is_cst(aff2))
2244                 isl_die(isl_aff_get_ctx(aff1), isl_error_invalid,
2245                         "at least one affine expression should be constant",
2246                         goto error);
2247
2248         aff1 = isl_aff_cow(aff1);
2249         if (!aff1 || !aff2)
2250                 goto error;
2251
2252         aff1 = isl_aff_scale(aff1, aff2->v->el[1]);
2253         aff1 = isl_aff_scale_down(aff1, aff2->v->el[0]);
2254
2255         isl_aff_free(aff2);
2256         return aff1;
2257 error:
2258         isl_aff_free(aff1);
2259         isl_aff_free(aff2);
2260         return NULL;
2261 }
2262
2263 /* Divide "aff1" by "aff2", assuming "aff2" is a piecewise constant.
2264  */
2265 __isl_give isl_aff *isl_aff_div(__isl_take isl_aff *aff1,
2266         __isl_take isl_aff *aff2)
2267 {
2268         int is_cst;
2269         int neg;
2270
2271         is_cst = isl_aff_is_cst(aff2);
2272         if (is_cst < 0)
2273                 goto error;
2274         if (!is_cst)
2275                 isl_die(isl_aff_get_ctx(aff2), isl_error_invalid,
2276                         "second argument should be a constant", goto error);
2277
2278         if (!aff2)
2279                 goto error;
2280
2281         neg = isl_int_is_neg(aff2->v->el[1]);
2282         if (neg) {
2283                 isl_int_neg(aff2->v->el[0], aff2->v->el[0]);
2284                 isl_int_neg(aff2->v->el[1], aff2->v->el[1]);
2285         }
2286
2287         aff1 = isl_aff_scale(aff1, aff2->v->el[0]);
2288         aff1 = isl_aff_scale_down(aff1, aff2->v->el[1]);
2289
2290         if (neg) {
2291                 isl_int_neg(aff2->v->el[0], aff2->v->el[0]);
2292                 isl_int_neg(aff2->v->el[1], aff2->v->el[1]);
2293         }
2294
2295         isl_aff_free(aff2);
2296         return aff1;
2297 error:
2298         isl_aff_free(aff1);
2299         isl_aff_free(aff2);
2300         return NULL;
2301 }
2302
2303 static __isl_give isl_pw_aff *pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2304         __isl_take isl_pw_aff *pwaff2)
2305 {
2306         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_add);
2307 }
2308
2309 __isl_give isl_pw_aff *isl_pw_aff_add(__isl_take isl_pw_aff *pwaff1,
2310         __isl_take isl_pw_aff *pwaff2)
2311 {
2312         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_add);
2313 }
2314
2315 __isl_give isl_pw_aff *isl_pw_aff_union_add(__isl_take isl_pw_aff *pwaff1,
2316         __isl_take isl_pw_aff *pwaff2)
2317 {
2318         return isl_pw_aff_union_add_(pwaff1, pwaff2);
2319 }
2320
2321 static __isl_give isl_pw_aff *pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2322         __isl_take isl_pw_aff *pwaff2)
2323 {
2324         return isl_pw_aff_on_shared_domain(pwaff1, pwaff2, &isl_aff_mul);
2325 }
2326
2327 __isl_give isl_pw_aff *isl_pw_aff_mul(__isl_take isl_pw_aff *pwaff1,
2328         __isl_take isl_pw_aff *pwaff2)
2329 {
2330         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_mul);
2331 }
2332
2333 static __isl_give isl_pw_aff *pw_aff_div(__isl_take isl_pw_aff *pa1,
2334         __isl_take isl_pw_aff *pa2)
2335 {
2336         return isl_pw_aff_on_shared_domain(pa1, pa2, &isl_aff_div);
2337 }
2338
2339 /* Divide "pa1" by "pa2", assuming "pa2" is a piecewise constant.
2340  */
2341 __isl_give isl_pw_aff *isl_pw_aff_div(__isl_take isl_pw_aff *pa1,
2342         __isl_take isl_pw_aff *pa2)
2343 {
2344         int is_cst;
2345
2346         is_cst = isl_pw_aff_is_cst(pa2);
2347         if (is_cst < 0)
2348                 goto error;
2349         if (!is_cst)
2350                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2351                         "second argument should be a piecewise constant",
2352                         goto error);
2353         return isl_pw_aff_align_params_pw_pw_and(pa1, pa2, &pw_aff_div);
2354 error:
2355         isl_pw_aff_free(pa1);
2356         isl_pw_aff_free(pa2);
2357         return NULL;
2358 }
2359
2360 /* Compute the quotient of the integer division of "pa1" by "pa2"
2361  * with rounding towards zero.
2362  * "pa2" is assumed to be a piecewise constant.
2363  *
2364  * In particular, return
2365  *
2366  *      pa1 >= 0 ? floor(pa1/pa2) : ceil(pa1/pa2)
2367  *
2368  */
2369 __isl_give isl_pw_aff *isl_pw_aff_tdiv_q(__isl_take isl_pw_aff *pa1,
2370         __isl_take isl_pw_aff *pa2)
2371 {
2372         int is_cst;
2373         isl_set *cond;
2374         isl_pw_aff *f, *c;
2375
2376         is_cst = isl_pw_aff_is_cst(pa2);
2377         if (is_cst < 0)
2378                 goto error;
2379         if (!is_cst)
2380                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2381                         "second argument should be a piecewise constant",
2382                         goto error);
2383
2384         pa1 = isl_pw_aff_div(pa1, pa2);
2385
2386         cond = isl_pw_aff_nonneg_set(isl_pw_aff_copy(pa1));
2387         f = isl_pw_aff_floor(isl_pw_aff_copy(pa1));
2388         c = isl_pw_aff_ceil(pa1);
2389         return isl_pw_aff_cond(isl_set_indicator_function(cond), f, c);
2390 error:
2391         isl_pw_aff_free(pa1);
2392         isl_pw_aff_free(pa2);
2393         return NULL;
2394 }
2395
2396 /* Compute the remainder of the integer division of "pa1" by "pa2"
2397  * with rounding towards zero.
2398  * "pa2" is assumed to be a piecewise constant.
2399  *
2400  * In particular, return
2401  *
2402  *      pa1 - pa2 * (pa1 >= 0 ? floor(pa1/pa2) : ceil(pa1/pa2))
2403  *
2404  */
2405 __isl_give isl_pw_aff *isl_pw_aff_tdiv_r(__isl_take isl_pw_aff *pa1,
2406         __isl_take isl_pw_aff *pa2)
2407 {
2408         int is_cst;
2409         isl_pw_aff *res;
2410
2411         is_cst = isl_pw_aff_is_cst(pa2);
2412         if (is_cst < 0)
2413                 goto error;
2414         if (!is_cst)
2415                 isl_die(isl_pw_aff_get_ctx(pa2), isl_error_invalid,
2416                         "second argument should be a piecewise constant",
2417                         goto error);
2418         res = isl_pw_aff_tdiv_q(isl_pw_aff_copy(pa1), isl_pw_aff_copy(pa2));
2419         res = isl_pw_aff_mul(pa2, res);
2420         res = isl_pw_aff_sub(pa1, res);
2421         return res;
2422 error:
2423         isl_pw_aff_free(pa1);
2424         isl_pw_aff_free(pa2);
2425         return NULL;
2426 }
2427
2428 static __isl_give isl_pw_aff *pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2429         __isl_take isl_pw_aff *pwaff2)
2430 {
2431         isl_set *le;
2432         isl_set *dom;
2433
2434         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2435                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2436         le = isl_pw_aff_le_set(isl_pw_aff_copy(pwaff1),
2437                                 isl_pw_aff_copy(pwaff2));
2438         dom = isl_set_subtract(dom, isl_set_copy(le));
2439         return isl_pw_aff_select(le, pwaff1, dom, pwaff2);
2440 }
2441
2442 __isl_give isl_pw_aff *isl_pw_aff_min(__isl_take isl_pw_aff *pwaff1,
2443         __isl_take isl_pw_aff *pwaff2)
2444 {
2445         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_min);
2446 }
2447
2448 static __isl_give isl_pw_aff *pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2449         __isl_take isl_pw_aff *pwaff2)
2450 {
2451         isl_set *ge;
2452         isl_set *dom;
2453
2454         dom = isl_set_intersect(isl_pw_aff_domain(isl_pw_aff_copy(pwaff1)),
2455                                 isl_pw_aff_domain(isl_pw_aff_copy(pwaff2)));
2456         ge = isl_pw_aff_ge_set(isl_pw_aff_copy(pwaff1),
2457                                 isl_pw_aff_copy(pwaff2));
2458         dom = isl_set_subtract(dom, isl_set_copy(ge));
2459         return isl_pw_aff_select(ge, pwaff1, dom, pwaff2);
2460 }
2461
2462 __isl_give isl_pw_aff *isl_pw_aff_max(__isl_take isl_pw_aff *pwaff1,
2463         __isl_take isl_pw_aff *pwaff2)
2464 {
2465         return isl_pw_aff_align_params_pw_pw_and(pwaff1, pwaff2, &pw_aff_max);
2466 }
2467
2468 static __isl_give isl_pw_aff *pw_aff_list_reduce(
2469         __isl_take isl_pw_aff_list *list,
2470         __isl_give isl_pw_aff *(*fn)(__isl_take isl_pw_aff *pwaff1,
2471                                         __isl_take isl_pw_aff *pwaff2))
2472 {
2473         int i;
2474         isl_ctx *ctx;
2475         isl_pw_aff *res;
2476
2477         if (!list)
2478                 return NULL;
2479
2480         ctx = isl_pw_aff_list_get_ctx(list);
2481         if (list->n < 1)
2482                 isl_die(ctx, isl_error_invalid,
2483                         "list should contain at least one element",
2484                         return isl_pw_aff_list_free(list));
2485
2486         res = isl_pw_aff_copy(list->p[0]);
2487         for (i = 1; i < list->n; ++i)
2488                 res = fn(res, isl_pw_aff_copy(list->p[i]));
2489
2490         isl_pw_aff_list_free(list);
2491         return res;
2492 }
2493
2494 /* Return an isl_pw_aff that maps each element in the intersection of the
2495  * domains of the elements of list to the minimal corresponding affine
2496  * expression.
2497  */
2498 __isl_give isl_pw_aff *isl_pw_aff_list_min(__isl_take isl_pw_aff_list *list)
2499 {
2500         return pw_aff_list_reduce(list, &isl_pw_aff_min);
2501 }
2502
2503 /* Return an isl_pw_aff that maps each element in the intersection of the
2504  * domains of the elements of list to the maximal corresponding affine
2505  * expression.
2506  */
2507 __isl_give isl_pw_aff *isl_pw_aff_list_max(__isl_take isl_pw_aff_list *list)
2508 {
2509         return pw_aff_list_reduce(list, &isl_pw_aff_max);
2510 }
2511
2512 #undef BASE
2513 #define BASE aff
2514
2515 #include <isl_multi_templ.c>
2516
2517 /* Create an isl_pw_multi_aff with the given isl_multi_aff on a universe
2518  * domain.
2519  */
2520 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_multi_aff(
2521         __isl_take isl_multi_aff *ma)
2522 {
2523         isl_set *dom = isl_set_universe(isl_multi_aff_get_domain_space(ma));
2524         return isl_pw_multi_aff_alloc(dom, ma);
2525 }
2526
2527 /* Create a piecewise multi-affine expression in the given space that maps each
2528  * input dimension to the corresponding output dimension.
2529  */
2530 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_identity(
2531         __isl_take isl_space *space)
2532 {
2533         return isl_pw_multi_aff_from_multi_aff(isl_multi_aff_identity(space));
2534 }
2535
2536 __isl_give isl_multi_aff *isl_multi_aff_add(__isl_take isl_multi_aff *maff1,
2537         __isl_take isl_multi_aff *maff2)
2538 {
2539         int i;
2540         isl_ctx *ctx;
2541
2542         maff1 = isl_multi_aff_cow(maff1);
2543         if (!maff1 || !maff2)
2544                 goto error;
2545
2546         ctx = isl_multi_aff_get_ctx(maff1);
2547         if (!isl_space_is_equal(maff1->space, maff2->space))
2548                 isl_die(ctx, isl_error_invalid,
2549                         "spaces don't match", goto error);
2550
2551         for (i = 0; i < maff1->n; ++i) {
2552                 maff1->p[i] = isl_aff_add(maff1->p[i],
2553                                             isl_aff_copy(maff2->p[i]));
2554                 if (!maff1->p[i])
2555                         goto error;
2556         }
2557
2558         isl_multi_aff_free(maff2);
2559         return maff1;
2560 error:
2561         isl_multi_aff_free(maff1);
2562         isl_multi_aff_free(maff2);
2563         return NULL;
2564 }
2565
2566 /* Given two multi-affine expressions A -> B and C -> D,
2567  * construct a multi-affine expression [A -> C] -> [B -> D].
2568  */
2569 __isl_give isl_multi_aff *isl_multi_aff_product(
2570         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
2571 {
2572         int i;
2573         isl_aff *aff;
2574         isl_space *space;
2575         isl_multi_aff *res;
2576         int in1, in2, out1, out2;
2577
2578         in1 = isl_multi_aff_dim(ma1, isl_dim_in);
2579         in2 = isl_multi_aff_dim(ma2, isl_dim_in);
2580         out1 = isl_multi_aff_dim(ma1, isl_dim_out);
2581         out2 = isl_multi_aff_dim(ma2, isl_dim_out);
2582         space = isl_space_product(isl_multi_aff_get_space(ma1),
2583                                   isl_multi_aff_get_space(ma2));
2584         res = isl_multi_aff_alloc(isl_space_copy(space));
2585         space = isl_space_domain(space);
2586
2587         for (i = 0; i < out1; ++i) {
2588                 aff = isl_multi_aff_get_aff(ma1, i);
2589                 aff = isl_aff_insert_dims(aff, isl_dim_in, in1, in2);
2590                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2591                 res = isl_multi_aff_set_aff(res, i, aff);
2592         }
2593
2594         for (i = 0; i < out2; ++i) {
2595                 aff = isl_multi_aff_get_aff(ma2, i);
2596                 aff = isl_aff_insert_dims(aff, isl_dim_in, 0, in1);
2597                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2598                 res = isl_multi_aff_set_aff(res, out1 + i, aff);
2599         }
2600
2601         isl_space_free(space);
2602         isl_multi_aff_free(ma1);
2603         isl_multi_aff_free(ma2);
2604         return res;
2605 }
2606
2607 /* Exploit the equalities in "eq" to simplify the affine expressions.
2608  */
2609 static __isl_give isl_multi_aff *isl_multi_aff_substitute_equalities(
2610         __isl_take isl_multi_aff *maff, __isl_take isl_basic_set *eq)
2611 {
2612         int i;
2613
2614         maff = isl_multi_aff_cow(maff);
2615         if (!maff || !eq)
2616                 goto error;
2617
2618         for (i = 0; i < maff->n; ++i) {
2619                 maff->p[i] = isl_aff_substitute_equalities(maff->p[i],
2620                                                     isl_basic_set_copy(eq));
2621                 if (!maff->p[i])
2622                         goto error;
2623         }
2624
2625         isl_basic_set_free(eq);
2626         return maff;
2627 error:
2628         isl_basic_set_free(eq);
2629         isl_multi_aff_free(maff);
2630         return NULL;
2631 }
2632
2633 __isl_give isl_multi_aff *isl_multi_aff_scale(__isl_take isl_multi_aff *maff,
2634         isl_int f)
2635 {
2636         int i;
2637
2638         maff = isl_multi_aff_cow(maff);
2639         if (!maff)
2640                 return NULL;
2641
2642         for (i = 0; i < maff->n; ++i) {
2643                 maff->p[i] = isl_aff_scale(maff->p[i], f);
2644                 if (!maff->p[i])
2645                         return isl_multi_aff_free(maff);
2646         }
2647
2648         return maff;
2649 }
2650
2651 __isl_give isl_multi_aff *isl_multi_aff_add_on_domain(__isl_keep isl_set *dom,
2652         __isl_take isl_multi_aff *maff1, __isl_take isl_multi_aff *maff2)
2653 {
2654         maff1 = isl_multi_aff_add(maff1, maff2);
2655         maff1 = isl_multi_aff_gist(maff1, isl_set_copy(dom));
2656         return maff1;
2657 }
2658
2659 int isl_multi_aff_is_empty(__isl_keep isl_multi_aff *maff)
2660 {
2661         if (!maff)
2662                 return -1;
2663
2664         return 0;
2665 }
2666
2667 int isl_multi_aff_plain_is_equal(__isl_keep isl_multi_aff *maff1,
2668         __isl_keep isl_multi_aff *maff2)
2669 {
2670         int i;
2671         int equal;
2672
2673         if (!maff1 || !maff2)
2674                 return -1;
2675         if (maff1->n != maff2->n)
2676                 return 0;
2677         equal = isl_space_is_equal(maff1->space, maff2->space);
2678         if (equal < 0 || !equal)
2679                 return equal;
2680
2681         for (i = 0; i < maff1->n; ++i) {
2682                 equal = isl_aff_plain_is_equal(maff1->p[i], maff2->p[i]);
2683                 if (equal < 0 || !equal)
2684                         return equal;
2685         }
2686
2687         return 1;
2688 }
2689
2690 __isl_give isl_multi_aff *isl_multi_aff_set_dim_name(
2691         __isl_take isl_multi_aff *maff,
2692         enum isl_dim_type type, unsigned pos, const char *s)
2693 {
2694         int i;
2695
2696         maff = isl_multi_aff_cow(maff);
2697         if (!maff)
2698                 return NULL;
2699
2700         maff->space = isl_space_set_dim_name(maff->space, type, pos, s);
2701         if (!maff->space)
2702                 return isl_multi_aff_free(maff);
2703
2704         if (type == isl_dim_out)
2705                 return maff;
2706         for (i = 0; i < maff->n; ++i) {
2707                 maff->p[i] = isl_aff_set_dim_name(maff->p[i], type, pos, s);
2708                 if (!maff->p[i])
2709                         return isl_multi_aff_free(maff);
2710         }
2711
2712         return maff;
2713 }
2714
2715 __isl_give isl_multi_aff *isl_multi_aff_drop_dims(__isl_take isl_multi_aff *maff,
2716         enum isl_dim_type type, unsigned first, unsigned n)
2717 {
2718         int i;
2719
2720         maff = isl_multi_aff_cow(maff);
2721         if (!maff)
2722                 return NULL;
2723
2724         maff->space = isl_space_drop_dims(maff->space, type, first, n);
2725         if (!maff->space)
2726                 return isl_multi_aff_free(maff);
2727
2728         if (type == isl_dim_out) {
2729                 for (i = 0; i < n; ++i)
2730                         isl_aff_free(maff->p[first + i]);
2731                 for (i = first; i + n < maff->n; ++i)
2732                         maff->p[i] = maff->p[i + n];
2733                 maff->n -= n;
2734                 return maff;
2735         }
2736
2737         for (i = 0; i < maff->n; ++i) {
2738                 maff->p[i] = isl_aff_drop_dims(maff->p[i], type, first, n);
2739                 if (!maff->p[i])
2740                         return isl_multi_aff_free(maff);
2741         }
2742
2743         return maff;
2744 }
2745
2746 /* Return the set of domain elements where "ma1" is lexicographically
2747  * smaller than or equal to "ma2".
2748  */
2749 __isl_give isl_set *isl_multi_aff_lex_le_set(__isl_take isl_multi_aff *ma1,
2750         __isl_take isl_multi_aff *ma2)
2751 {
2752         return isl_multi_aff_lex_ge_set(ma2, ma1);
2753 }
2754
2755 /* Return the set of domain elements where "ma1" is lexicographically
2756  * greater than or equal to "ma2".
2757  */
2758 __isl_give isl_set *isl_multi_aff_lex_ge_set(__isl_take isl_multi_aff *ma1,
2759         __isl_take isl_multi_aff *ma2)
2760 {
2761         isl_space *space;
2762         isl_map *map1, *map2;
2763         isl_map *map, *ge;
2764
2765         map1 = isl_map_from_multi_aff(ma1);
2766         map2 = isl_map_from_multi_aff(ma2);
2767         map = isl_map_range_product(map1, map2);
2768         space = isl_space_range(isl_map_get_space(map));
2769         space = isl_space_domain(isl_space_unwrap(space));
2770         ge = isl_map_lex_ge(space);
2771         map = isl_map_intersect_range(map, isl_map_wrap(ge));
2772
2773         return isl_map_domain(map);
2774 }
2775
2776 #undef PW
2777 #define PW isl_pw_multi_aff
2778 #undef EL
2779 #define EL isl_multi_aff
2780 #undef EL_IS_ZERO
2781 #define EL_IS_ZERO is_empty
2782 #undef ZERO
2783 #define ZERO empty
2784 #undef IS_ZERO
2785 #define IS_ZERO is_empty
2786 #undef FIELD
2787 #define FIELD maff
2788 #undef DEFAULT_IS_ZERO
2789 #define DEFAULT_IS_ZERO 0
2790
2791 #define NO_NEG
2792 #define NO_EVAL
2793 #define NO_OPT
2794 #define NO_INVOLVES_DIMS
2795 #define NO_MOVE_DIMS
2796 #define NO_INSERT_DIMS
2797 #define NO_LIFT
2798 #define NO_MORPH
2799
2800 #include <isl_pw_templ.c>
2801
2802 #undef UNION
2803 #define UNION isl_union_pw_multi_aff
2804 #undef PART
2805 #define PART isl_pw_multi_aff
2806 #undef PARTS
2807 #define PARTS pw_multi_aff
2808 #define ALIGN_DOMAIN
2809
2810 #define NO_EVAL
2811
2812 #include <isl_union_templ.c>
2813
2814 /* Given a function "cmp" that returns the set of elements where
2815  * "ma1" is "better" than "ma2", return the intersection of this
2816  * set with "dom1" and "dom2".
2817  */
2818 static __isl_give isl_set *shared_and_better(__isl_keep isl_set *dom1,
2819         __isl_keep isl_set *dom2, __isl_keep isl_multi_aff *ma1,
2820         __isl_keep isl_multi_aff *ma2,
2821         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2822                                     __isl_take isl_multi_aff *ma2))
2823 {
2824         isl_set *common;
2825         isl_set *better;
2826         int is_empty;
2827
2828         common = isl_set_intersect(isl_set_copy(dom1), isl_set_copy(dom2));
2829         is_empty = isl_set_plain_is_empty(common);
2830         if (is_empty >= 0 && is_empty)
2831                 return common;
2832         if (is_empty < 0)
2833                 return isl_set_free(common);
2834         better = cmp(isl_multi_aff_copy(ma1), isl_multi_aff_copy(ma2));
2835         better = isl_set_intersect(common, better);
2836
2837         return better;
2838 }
2839
2840 /* Given a function "cmp" that returns the set of elements where
2841  * "ma1" is "better" than "ma2", return a piecewise multi affine
2842  * expression defined on the union of the definition domains
2843  * of "pma1" and "pma2" that maps to the "best" of "pma1" and
2844  * "pma2" on each cell.  If only one of the two input functions
2845  * is defined on a given cell, then it is considered the best.
2846  */
2847 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_opt(
2848         __isl_take isl_pw_multi_aff *pma1,
2849         __isl_take isl_pw_multi_aff *pma2,
2850         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2851                                     __isl_take isl_multi_aff *ma2))
2852 {
2853         int i, j, n;
2854         isl_pw_multi_aff *res = NULL;
2855         isl_ctx *ctx;
2856         isl_set *set = NULL;
2857
2858         if (!pma1 || !pma2)
2859                 goto error;
2860
2861         ctx = isl_space_get_ctx(pma1->dim);
2862         if (!isl_space_is_equal(pma1->dim, pma2->dim))
2863                 isl_die(ctx, isl_error_invalid,
2864                         "arguments should live in the same space", goto error);
2865
2866         if (isl_pw_multi_aff_is_empty(pma1)) {
2867                 isl_pw_multi_aff_free(pma1);
2868                 return pma2;
2869         }
2870
2871         if (isl_pw_multi_aff_is_empty(pma2)) {
2872                 isl_pw_multi_aff_free(pma2);
2873                 return pma1;
2874         }
2875
2876         n = 2 * (pma1->n + 1) * (pma2->n + 1);
2877         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma1->dim), n);
2878
2879         for (i = 0; i < pma1->n; ++i) {
2880                 set = isl_set_copy(pma1->p[i].set);
2881                 for (j = 0; j < pma2->n; ++j) {
2882                         isl_set *better;
2883                         int is_empty;
2884
2885                         better = shared_and_better(pma2->p[j].set,
2886                                         pma1->p[i].set, pma2->p[j].maff,
2887                                         pma1->p[i].maff, cmp);
2888                         is_empty = isl_set_plain_is_empty(better);
2889                         if (is_empty < 0 || is_empty) {
2890                                 isl_set_free(better);
2891                                 if (is_empty < 0)
2892                                         goto error;
2893                                 continue;
2894                         }
2895                         set = isl_set_subtract(set, isl_set_copy(better));
2896
2897                         res = isl_pw_multi_aff_add_piece(res, better,
2898                                         isl_multi_aff_copy(pma2->p[j].maff));
2899                 }
2900                 res = isl_pw_multi_aff_add_piece(res, set,
2901                                         isl_multi_aff_copy(pma1->p[i].maff));
2902         }
2903
2904         for (j = 0; j < pma2->n; ++j) {
2905                 set = isl_set_copy(pma2->p[j].set);
2906                 for (i = 0; i < pma1->n; ++i)
2907                         set = isl_set_subtract(set,
2908                                         isl_set_copy(pma1->p[i].set));
2909                 res = isl_pw_multi_aff_add_piece(res, set,
2910                                         isl_multi_aff_copy(pma2->p[j].maff));
2911         }
2912
2913         isl_pw_multi_aff_free(pma1);
2914         isl_pw_multi_aff_free(pma2);
2915
2916         return res;
2917 error:
2918         isl_pw_multi_aff_free(pma1);
2919         isl_pw_multi_aff_free(pma2);
2920         isl_set_free(set);
2921         return isl_pw_multi_aff_free(res);
2922 }
2923
2924 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmax(
2925         __isl_take isl_pw_multi_aff *pma1,
2926         __isl_take isl_pw_multi_aff *pma2)
2927 {
2928         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_ge_set);
2929 }
2930
2931 /* Given two piecewise multi affine expressions, return a piecewise
2932  * multi-affine expression defined on the union of the definition domains
2933  * of the inputs that is equal to the lexicographic maximum of the two
2934  * inputs on each cell.  If only one of the two inputs is defined on
2935  * a given cell, then it is considered to be the maximum.
2936  */
2937 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmax(
2938         __isl_take isl_pw_multi_aff *pma1,
2939         __isl_take isl_pw_multi_aff *pma2)
2940 {
2941         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2942                                                     &pw_multi_aff_union_lexmax);
2943 }
2944
2945 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmin(
2946         __isl_take isl_pw_multi_aff *pma1,
2947         __isl_take isl_pw_multi_aff *pma2)
2948 {
2949         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_le_set);
2950 }
2951
2952 /* Given two piecewise multi affine expressions, return a piecewise
2953  * multi-affine expression defined on the union of the definition domains
2954  * of the inputs that is equal to the lexicographic minimum of the two
2955  * inputs on each cell.  If only one of the two inputs is defined on
2956  * a given cell, then it is considered to be the minimum.
2957  */
2958 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmin(
2959         __isl_take isl_pw_multi_aff *pma1,
2960         __isl_take isl_pw_multi_aff *pma2)
2961 {
2962         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2963                                                     &pw_multi_aff_union_lexmin);
2964 }
2965
2966 static __isl_give isl_pw_multi_aff *pw_multi_aff_add(
2967         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2968 {
2969         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
2970                                                 &isl_multi_aff_add);
2971 }
2972
2973 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_add(
2974         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2975 {
2976         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2977                                                 &pw_multi_aff_add);
2978 }
2979
2980 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_add(
2981         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2982 {
2983         return isl_pw_multi_aff_union_add_(pma1, pma2);
2984 }
2985
2986 /* Given two piecewise multi-affine expressions A -> B and C -> D,
2987  * construct a piecewise multi-affine expression [A -> C] -> [B -> D].
2988  */
2989 static __isl_give isl_pw_multi_aff *pw_multi_aff_product(
2990         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2991 {
2992         int i, j, n;
2993         isl_space *space;
2994         isl_pw_multi_aff *res;
2995
2996         if (!pma1 || !pma2)
2997                 goto error;
2998
2999         n = pma1->n * pma2->n;
3000         space = isl_space_product(isl_space_copy(pma1->dim),
3001                                   isl_space_copy(pma2->dim));
3002         res = isl_pw_multi_aff_alloc_size(space, n);
3003
3004         for (i = 0; i < pma1->n; ++i) {
3005                 for (j = 0; j < pma2->n; ++j) {
3006                         isl_set *domain;
3007                         isl_multi_aff *ma;
3008
3009                         domain = isl_set_product(isl_set_copy(pma1->p[i].set),
3010                                                  isl_set_copy(pma2->p[j].set));
3011                         ma = isl_multi_aff_product(
3012                                         isl_multi_aff_copy(pma1->p[i].maff),
3013                                         isl_multi_aff_copy(pma2->p[i].maff));
3014                         res = isl_pw_multi_aff_add_piece(res, domain, ma);
3015                 }
3016         }
3017
3018         isl_pw_multi_aff_free(pma1);
3019         isl_pw_multi_aff_free(pma2);
3020         return res;
3021 error:
3022         isl_pw_multi_aff_free(pma1);
3023         isl_pw_multi_aff_free(pma2);
3024         return NULL;
3025 }
3026
3027 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_product(
3028         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3029 {
3030         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3031                                                 &pw_multi_aff_product);
3032 }
3033
3034 /* Construct a map mapping the domain of the piecewise multi-affine expression
3035  * to its range, with each dimension in the range equated to the
3036  * corresponding affine expression on its cell.
3037  */
3038 __isl_give isl_map *isl_map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3039 {
3040         int i;
3041         isl_map *map;
3042
3043         if (!pma)
3044                 return NULL;
3045
3046         map = isl_map_empty(isl_pw_multi_aff_get_space(pma));
3047
3048         for (i = 0; i < pma->n; ++i) {
3049                 isl_multi_aff *maff;
3050                 isl_basic_map *bmap;
3051                 isl_map *map_i;
3052
3053                 maff = isl_multi_aff_copy(pma->p[i].maff);
3054                 bmap = isl_basic_map_from_multi_aff(maff);
3055                 map_i = isl_map_from_basic_map(bmap);
3056                 map_i = isl_map_intersect_domain(map_i,
3057                                                 isl_set_copy(pma->p[i].set));
3058                 map = isl_map_union_disjoint(map, map_i);
3059         }
3060
3061         isl_pw_multi_aff_free(pma);
3062         return map;
3063 }
3064
3065 __isl_give isl_set *isl_set_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3066 {
3067         if (!pma)
3068                 return NULL;
3069
3070         if (!isl_space_is_set(pma->dim))
3071                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3072                         "isl_pw_multi_aff cannot be converted into an isl_set",
3073                         return isl_pw_multi_aff_free(pma));
3074
3075         return isl_map_from_pw_multi_aff(pma);
3076 }
3077
3078 /* Given a basic map with a single output dimension that is defined
3079  * in terms of the parameters and input dimensions using an equality,
3080  * extract an isl_aff that expresses the output dimension in terms
3081  * of the parameters and input dimensions.
3082  *
3083  * Since some applications expect the result of isl_pw_multi_aff_from_map
3084  * to only contain integer affine expressions, we compute the floor
3085  * of the expression before returning.
3086  *
3087  * This function shares some similarities with
3088  * isl_basic_map_has_defining_equality and isl_constraint_get_bound.
3089  */
3090 static __isl_give isl_aff *extract_isl_aff_from_basic_map(
3091         __isl_take isl_basic_map *bmap)
3092 {
3093         int i;
3094         unsigned offset;
3095         unsigned total;
3096         isl_local_space *ls;
3097         isl_aff *aff;
3098
3099         if (!bmap)
3100                 return NULL;
3101         if (isl_basic_map_dim(bmap, isl_dim_out) != 1)
3102                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3103                         "basic map should have a single output dimension",
3104                         goto error);
3105         offset = isl_basic_map_offset(bmap, isl_dim_out);
3106         total = isl_basic_map_total_dim(bmap);
3107         for (i = 0; i < bmap->n_eq; ++i) {
3108                 if (isl_int_is_zero(bmap->eq[i][offset]))
3109                         continue;
3110                 if (isl_seq_first_non_zero(bmap->eq[i] + offset + 1,
3111                                            1 + total - (offset + 1)) != -1)
3112                         continue;
3113                 break;
3114         }
3115         if (i >= bmap->n_eq)
3116                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3117                         "unable to find suitable equality", goto error);
3118         ls = isl_basic_map_get_local_space(bmap);
3119         aff = isl_aff_alloc(isl_local_space_domain(ls));
3120         if (!aff)
3121                 goto error;
3122         if (isl_int_is_neg(bmap->eq[i][offset]))
3123                 isl_seq_cpy(aff->v->el + 1, bmap->eq[i], offset);
3124         else
3125                 isl_seq_neg(aff->v->el + 1, bmap->eq[i], offset);
3126         isl_seq_clr(aff->v->el + 1 + offset, aff->v->size - (1 + offset));
3127         isl_int_abs(aff->v->el[0], bmap->eq[i][offset]);
3128         isl_basic_map_free(bmap);
3129
3130         aff = isl_aff_remove_unused_divs(aff);
3131         aff = isl_aff_floor(aff);
3132         return aff;
3133 error:
3134         isl_basic_map_free(bmap);
3135         return NULL;
3136 }
3137
3138 /* Given a basic map where each output dimension is defined
3139  * in terms of the parameters and input dimensions using an equality,
3140  * extract an isl_multi_aff that expresses the output dimensions in terms
3141  * of the parameters and input dimensions.
3142  */
3143 static __isl_give isl_multi_aff *extract_isl_multi_aff_from_basic_map(
3144         __isl_take isl_basic_map *bmap)
3145 {
3146         int i;
3147         unsigned n_out;
3148         isl_multi_aff *ma;
3149
3150         if (!bmap)
3151                 return NULL;
3152
3153         ma = isl_multi_aff_alloc(isl_basic_map_get_space(bmap));
3154         n_out = isl_basic_map_dim(bmap, isl_dim_out);
3155
3156         for (i = 0; i < n_out; ++i) {
3157                 isl_basic_map *bmap_i;
3158                 isl_aff *aff;
3159
3160                 bmap_i = isl_basic_map_copy(bmap);
3161                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out,
3162                                                         i + 1, n_out - (1 + i));
3163                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out, 0, i);
3164                 aff = extract_isl_aff_from_basic_map(bmap_i);
3165                 ma = isl_multi_aff_set_aff(ma, i, aff);
3166         }
3167
3168         isl_basic_map_free(bmap);
3169
3170         return ma;
3171 }
3172
3173 /* Create an isl_pw_multi_aff that is equivalent to
3174  * isl_map_intersect_domain(isl_map_from_basic_map(bmap), domain).
3175  * The given basic map is such that each output dimension is defined
3176  * in terms of the parameters and input dimensions using an equality.
3177  */
3178 static __isl_give isl_pw_multi_aff *plain_pw_multi_aff_from_map(
3179         __isl_take isl_set *domain, __isl_take isl_basic_map *bmap)
3180 {
3181         isl_multi_aff *ma;
3182
3183         ma = extract_isl_multi_aff_from_basic_map(bmap);
3184         return isl_pw_multi_aff_alloc(domain, ma);
3185 }
3186
3187 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3188  * This obviously only works if the input "map" is single-valued.
3189  * If so, we compute the lexicographic minimum of the image in the form
3190  * of an isl_pw_multi_aff.  Since the image is unique, it is equal
3191  * to its lexicographic minimum.
3192  * If the input is not single-valued, we produce an error.
3193  *
3194  * As a special case, we first check if all output dimensions are uniquely
3195  * defined in terms of the parameters and input dimensions over the entire
3196  * domain.  If so, we extract the desired isl_pw_multi_aff directly
3197  * from the affine hull of "map" and its domain.
3198  */
3199 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_map(__isl_take isl_map *map)
3200 {
3201         int i;
3202         int sv;
3203         isl_pw_multi_aff *pma;
3204         isl_basic_map *hull;
3205
3206         if (!map)
3207                 return NULL;
3208
3209         hull = isl_map_affine_hull(isl_map_copy(map));
3210         sv = isl_basic_map_plain_is_single_valued(hull);
3211         if (sv >= 0 && sv)
3212                 return plain_pw_multi_aff_from_map(isl_map_domain(map), hull);
3213         isl_basic_map_free(hull);
3214         if (sv < 0)
3215                 goto error;
3216
3217         sv = isl_map_is_single_valued(map);
3218         if (sv < 0)
3219                 goto error;
3220         if (!sv)
3221                 isl_die(isl_map_get_ctx(map), isl_error_invalid,
3222                         "map is not single-valued", goto error);
3223         map = isl_map_make_disjoint(map);
3224         if (!map)
3225                 return NULL;
3226
3227         pma = isl_pw_multi_aff_empty(isl_map_get_space(map));
3228
3229         for (i = 0; i < map->n; ++i) {
3230                 isl_pw_multi_aff *pma_i;
3231                 isl_basic_map *bmap;
3232                 bmap = isl_basic_map_copy(map->p[i]);
3233                 pma_i = isl_basic_map_lexmin_pw_multi_aff(bmap);
3234                 pma = isl_pw_multi_aff_add_disjoint(pma, pma_i);
3235         }
3236
3237         isl_map_free(map);
3238         return pma;
3239 error:
3240         isl_map_free(map);
3241         return NULL;
3242 }
3243
3244 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_set(__isl_take isl_set *set)
3245 {
3246         return isl_pw_multi_aff_from_map(set);
3247 }
3248
3249 /* Return the piecewise affine expression "set ? 1 : 0".
3250  */
3251 __isl_give isl_pw_aff *isl_set_indicator_function(__isl_take isl_set *set)
3252 {
3253         isl_pw_aff *pa;
3254         isl_space *space = isl_set_get_space(set);
3255         isl_local_space *ls = isl_local_space_from_space(space);
3256         isl_aff *zero = isl_aff_zero_on_domain(isl_local_space_copy(ls));
3257         isl_aff *one = isl_aff_zero_on_domain(ls);
3258
3259         one = isl_aff_add_constant_si(one, 1);
3260         pa = isl_pw_aff_alloc(isl_set_copy(set), one);
3261         set = isl_set_complement(set);
3262         pa = isl_pw_aff_add_disjoint(pa, isl_pw_aff_alloc(set, zero));
3263
3264         return pa;
3265 }
3266
3267 /* Plug in "subs" for dimension "type", "pos" of "aff".
3268  *
3269  * Let i be the dimension to replace and let "subs" be of the form
3270  *
3271  *      f/d
3272  *
3273  * and "aff" of the form
3274  *
3275  *      (a i + g)/m
3276  *
3277  * The result is
3278  *
3279  *      (a f + d g')/(m d)
3280  *
3281  * where g' is the result of plugging in "subs" in each of the integer
3282  * divisions in g.
3283  */
3284 __isl_give isl_aff *isl_aff_substitute(__isl_take isl_aff *aff,
3285         enum isl_dim_type type, unsigned pos, __isl_keep isl_aff *subs)
3286 {
3287         isl_ctx *ctx;
3288         isl_int v;
3289
3290         aff = isl_aff_cow(aff);
3291         if (!aff || !subs)
3292                 return isl_aff_free(aff);
3293
3294         ctx = isl_aff_get_ctx(aff);
3295         if (!isl_space_is_equal(aff->ls->dim, subs->ls->dim))
3296                 isl_die(ctx, isl_error_invalid,
3297                         "spaces don't match", return isl_aff_free(aff));
3298         if (isl_local_space_dim(subs->ls, isl_dim_div) != 0)
3299                 isl_die(ctx, isl_error_unsupported,
3300                         "cannot handle divs yet", return isl_aff_free(aff));
3301
3302         aff->ls = isl_local_space_substitute(aff->ls, type, pos, subs);
3303         if (!aff->ls)
3304                 return isl_aff_free(aff);
3305
3306         aff->v = isl_vec_cow(aff->v);
3307         if (!aff->v)
3308                 return isl_aff_free(aff);
3309
3310         pos += isl_local_space_offset(aff->ls, type);
3311
3312         isl_int_init(v);
3313         isl_seq_substitute(aff->v->el, pos, subs->v->el,
3314                             aff->v->size, subs->v->size, v);
3315         isl_int_clear(v);
3316
3317         return aff;
3318 }
3319
3320 /* Plug in "subs" for dimension "type", "pos" in each of the affine
3321  * expressions in "maff".
3322  */
3323 __isl_give isl_multi_aff *isl_multi_aff_substitute(
3324         __isl_take isl_multi_aff *maff, enum isl_dim_type type, unsigned pos,
3325         __isl_keep isl_aff *subs)
3326 {
3327         int i;
3328
3329         maff = isl_multi_aff_cow(maff);
3330         if (!maff || !subs)
3331                 return isl_multi_aff_free(maff);
3332
3333         if (type == isl_dim_in)
3334                 type = isl_dim_set;
3335
3336         for (i = 0; i < maff->n; ++i) {
3337                 maff->p[i] = isl_aff_substitute(maff->p[i], type, pos, subs);
3338                 if (!maff->p[i])
3339                         return isl_multi_aff_free(maff);
3340         }
3341
3342         return maff;
3343 }
3344
3345 /* Plug in "subs" for dimension "type", "pos" of "pma".
3346  *
3347  * pma is of the form
3348  *
3349  *      A_i(v) -> M_i(v)
3350  *
3351  * while subs is of the form
3352  *
3353  *      v' = B_j(v) -> S_j
3354  *
3355  * Each pair i,j such that C_ij = A_i \cap B_i is non-empty
3356  * has a contribution in the result, in particular
3357  *
3358  *      C_ij(S_j) -> M_i(S_j)
3359  *
3360  * Note that plugging in S_j in C_ij may also result in an empty set
3361  * and this contribution should simply be discarded.
3362  */
3363 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_substitute(
3364         __isl_take isl_pw_multi_aff *pma, enum isl_dim_type type, unsigned pos,
3365         __isl_keep isl_pw_aff *subs)
3366 {
3367         int i, j, n;
3368         isl_pw_multi_aff *res;
3369
3370         if (!pma || !subs)
3371                 return isl_pw_multi_aff_free(pma);
3372
3373         n = pma->n * subs->n;
3374         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma->dim), n);
3375
3376         for (i = 0; i < pma->n; ++i) {
3377                 for (j = 0; j < subs->n; ++j) {
3378                         isl_set *common;
3379                         isl_multi_aff *res_ij;
3380                         common = isl_set_intersect(
3381                                         isl_set_copy(pma->p[i].set),
3382                                         isl_set_copy(subs->p[j].set));
3383                         common = isl_set_substitute(common,
3384                                         type, pos, subs->p[j].aff);
3385                         if (isl_set_plain_is_empty(common)) {
3386                                 isl_set_free(common);
3387                                 continue;
3388                         }
3389
3390                         res_ij = isl_multi_aff_substitute(
3391                                         isl_multi_aff_copy(pma->p[i].maff),
3392                                         type, pos, subs->p[j].aff);
3393
3394                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3395                 }
3396         }
3397
3398         isl_pw_multi_aff_free(pma);
3399         return res;
3400 }
3401
3402 /* Extend the local space of "dst" to include the divs
3403  * in the local space of "src".
3404  */
3405 __isl_give isl_aff *isl_aff_align_divs(__isl_take isl_aff *dst,
3406         __isl_keep isl_aff *src)
3407 {
3408         isl_ctx *ctx;
3409         int *exp1 = NULL;
3410         int *exp2 = NULL;
3411         isl_mat *div;
3412
3413         if (!src || !dst)
3414                 return isl_aff_free(dst);
3415
3416         ctx = isl_aff_get_ctx(src);
3417         if (!isl_space_is_equal(src->ls->dim, dst->ls->dim))
3418                 isl_die(ctx, isl_error_invalid,
3419                         "spaces don't match", goto error);
3420
3421         if (src->ls->div->n_row == 0)
3422                 return dst;
3423
3424         exp1 = isl_alloc_array(ctx, int, src->ls->div->n_row);
3425         exp2 = isl_alloc_array(ctx, int, dst->ls->div->n_row);
3426         if (!exp1 || !exp2)
3427                 goto error;
3428
3429         div = isl_merge_divs(src->ls->div, dst->ls->div, exp1, exp2);
3430         dst = isl_aff_expand_divs(dst, div, exp2);
3431         free(exp1);
3432         free(exp2);
3433
3434         return dst;
3435 error:
3436         free(exp1);
3437         free(exp2);
3438         return isl_aff_free(dst);
3439 }
3440
3441 /* Adjust the local spaces of the affine expressions in "maff"
3442  * such that they all have the save divs.
3443  */
3444 __isl_give isl_multi_aff *isl_multi_aff_align_divs(
3445         __isl_take isl_multi_aff *maff)
3446 {
3447         int i;
3448
3449         if (!maff)
3450                 return NULL;
3451         if (maff->n == 0)
3452                 return maff;
3453         maff = isl_multi_aff_cow(maff);
3454         if (!maff)
3455                 return NULL;
3456
3457         for (i = 1; i < maff->n; ++i)
3458                 maff->p[0] = isl_aff_align_divs(maff->p[0], maff->p[i]);
3459         for (i = 1; i < maff->n; ++i) {
3460                 maff->p[i] = isl_aff_align_divs(maff->p[i], maff->p[0]);
3461                 if (!maff->p[i])
3462                         return isl_multi_aff_free(maff);
3463         }
3464
3465         return maff;
3466 }
3467
3468 __isl_give isl_aff *isl_aff_lift(__isl_take isl_aff *aff)
3469 {
3470         aff = isl_aff_cow(aff);
3471         if (!aff)
3472                 return NULL;
3473
3474         aff->ls = isl_local_space_lift(aff->ls);
3475         if (!aff->ls)
3476                 return isl_aff_free(aff);
3477
3478         return aff;
3479 }
3480
3481 /* Lift "maff" to a space with extra dimensions such that the result
3482  * has no more existentially quantified variables.
3483  * If "ls" is not NULL, then *ls is assigned the local space that lies
3484  * at the basis of the lifting applied to "maff".
3485  */
3486 __isl_give isl_multi_aff *isl_multi_aff_lift(__isl_take isl_multi_aff *maff,
3487         __isl_give isl_local_space **ls)
3488 {
3489         int i;
3490         isl_space *space;
3491         unsigned n_div;
3492
3493         if (ls)
3494                 *ls = NULL;
3495
3496         if (!maff)
3497                 return NULL;
3498
3499         if (maff->n == 0) {
3500                 if (ls) {
3501                         isl_space *space = isl_multi_aff_get_domain_space(maff);
3502                         *ls = isl_local_space_from_space(space);
3503                         if (!*ls)
3504                                 return isl_multi_aff_free(maff);
3505                 }
3506                 return maff;
3507         }
3508
3509         maff = isl_multi_aff_cow(maff);
3510         maff = isl_multi_aff_align_divs(maff);
3511         if (!maff)
3512                 return NULL;
3513
3514         n_div = isl_aff_dim(maff->p[0], isl_dim_div);
3515         space = isl_multi_aff_get_space(maff);
3516         space = isl_space_lift(isl_space_domain(space), n_div);
3517         space = isl_space_extend_domain_with_range(space,
3518                                                 isl_multi_aff_get_space(maff));
3519         if (!space)
3520                 return isl_multi_aff_free(maff);
3521         isl_space_free(maff->space);
3522         maff->space = space;
3523
3524         if (ls) {
3525                 *ls = isl_aff_get_domain_local_space(maff->p[0]);
3526                 if (!*ls)
3527                         return isl_multi_aff_free(maff);
3528         }
3529
3530         for (i = 0; i < maff->n; ++i) {
3531                 maff->p[i] = isl_aff_lift(maff->p[i]);
3532                 if (!maff->p[i])
3533                         goto error;
3534         }
3535
3536         return maff;
3537 error:
3538         if (ls)
3539                 isl_local_space_free(*ls);
3540         return isl_multi_aff_free(maff);
3541 }
3542
3543
3544 /* Extract an isl_pw_aff corresponding to output dimension "pos" of "pma".
3545  */
3546 __isl_give isl_pw_aff *isl_pw_multi_aff_get_pw_aff(
3547         __isl_keep isl_pw_multi_aff *pma, int pos)
3548 {
3549         int i;
3550         int n_out;
3551         isl_space *space;
3552         isl_pw_aff *pa;
3553
3554         if (!pma)
3555                 return NULL;
3556
3557         n_out = isl_pw_multi_aff_dim(pma, isl_dim_out);
3558         if (pos < 0 || pos >= n_out)
3559                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3560                         "index out of bounds", return NULL);
3561
3562         space = isl_pw_multi_aff_get_space(pma);
3563         space = isl_space_drop_dims(space, isl_dim_out,
3564                                     pos + 1, n_out - pos - 1);
3565         space = isl_space_drop_dims(space, isl_dim_out, 0, pos);
3566
3567         pa = isl_pw_aff_alloc_size(space, pma->n);
3568         for (i = 0; i < pma->n; ++i) {
3569                 isl_aff *aff;
3570                 aff = isl_multi_aff_get_aff(pma->p[i].maff, pos);
3571                 pa = isl_pw_aff_add_piece(pa, isl_set_copy(pma->p[i].set), aff);
3572         }
3573
3574         return pa;
3575 }
3576
3577 /* Return an isl_pw_multi_aff with the given "set" as domain and
3578  * an unnamed zero-dimensional range.
3579  */
3580 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_domain(
3581         __isl_take isl_set *set)
3582 {
3583         isl_multi_aff *ma;
3584         isl_space *space;
3585
3586         space = isl_set_get_space(set);
3587         space = isl_space_from_domain(space);
3588         ma = isl_multi_aff_zero(space);
3589         return isl_pw_multi_aff_alloc(set, ma);
3590 }
3591
3592 /* Add an isl_pw_multi_aff with the given "set" as domain and
3593  * an unnamed zero-dimensional range to *user.
3594  */
3595 static int add_pw_multi_aff_from_domain(__isl_take isl_set *set, void *user)
3596 {
3597         isl_union_pw_multi_aff **upma = user;
3598         isl_pw_multi_aff *pma;
3599
3600         pma = isl_pw_multi_aff_from_domain(set);
3601         *upma = isl_union_pw_multi_aff_add_pw_multi_aff(*upma, pma);
3602
3603         return 0;
3604 }
3605
3606 /* Return an isl_union_pw_multi_aff with the given "uset" as domain and
3607  * an unnamed zero-dimensional range.
3608  */
3609 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_from_domain(
3610         __isl_take isl_union_set *uset)
3611 {
3612         isl_space *space;
3613         isl_union_pw_multi_aff *upma;
3614
3615         if (!uset)
3616                 return NULL;
3617
3618         space = isl_union_set_get_space(uset);
3619         upma = isl_union_pw_multi_aff_empty(space);
3620
3621         if (isl_union_set_foreach_set(uset,
3622                                     &add_pw_multi_aff_from_domain, &upma) < 0)
3623                 goto error;
3624
3625         isl_union_set_free(uset);
3626         return upma;
3627 error:
3628         isl_union_set_free(uset);
3629         isl_union_pw_multi_aff_free(upma);
3630         return NULL;
3631 }
3632
3633 /* Convert "pma" to an isl_map and add it to *umap.
3634  */
3635 static int map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma, void *user)
3636 {
3637         isl_union_map **umap = user;
3638         isl_map *map;
3639
3640         map = isl_map_from_pw_multi_aff(pma);
3641         *umap = isl_union_map_add_map(*umap, map);
3642
3643         return 0;
3644 }
3645
3646 /* Construct a union map mapping the domain of the union
3647  * piecewise multi-affine expression to its range, with each dimension
3648  * in the range equated to the corresponding affine expression on its cell.
3649  */
3650 __isl_give isl_union_map *isl_union_map_from_union_pw_multi_aff(
3651         __isl_take isl_union_pw_multi_aff *upma)
3652 {
3653         isl_space *space;
3654         isl_union_map *umap;
3655
3656         if (!upma)
3657                 return NULL;
3658
3659         space = isl_union_pw_multi_aff_get_space(upma);
3660         umap = isl_union_map_empty(space);
3661
3662         if (isl_union_pw_multi_aff_foreach_pw_multi_aff(upma,
3663                                         &map_from_pw_multi_aff, &umap) < 0)
3664                 goto error;
3665
3666         isl_union_pw_multi_aff_free(upma);
3667         return umap;
3668 error:
3669         isl_union_pw_multi_aff_free(upma);
3670         isl_union_map_free(umap);
3671         return NULL;
3672 }
3673
3674 /* Local data for bin_entry and the callback "fn".
3675  */
3676 struct isl_union_pw_multi_aff_bin_data {
3677         isl_union_pw_multi_aff *upma2;
3678         isl_union_pw_multi_aff *res;
3679         isl_pw_multi_aff *pma;
3680         int (*fn)(void **entry, void *user);
3681 };
3682
3683 /* Given an isl_pw_multi_aff from upma1, store it in data->pma
3684  * and call data->fn for each isl_pw_multi_aff in data->upma2.
3685  */
3686 static int bin_entry(void **entry, void *user)
3687 {
3688         struct isl_union_pw_multi_aff_bin_data *data = user;
3689         isl_pw_multi_aff *pma = *entry;
3690
3691         data->pma = pma;
3692         if (isl_hash_table_foreach(data->upma2->dim->ctx, &data->upma2->table,
3693                                    data->fn, data) < 0)
3694                 return -1;
3695
3696         return 0;
3697 }
3698
3699 /* Call "fn" on each pair of isl_pw_multi_affs in "upma1" and "upma2".
3700  * The isl_pw_multi_aff from upma1 is stored in data->pma (where data is
3701  * passed as user field) and the isl_pw_multi_aff from upma2 is available
3702  * as *entry.  The callback should adjust data->res if desired.
3703  */
3704 static __isl_give isl_union_pw_multi_aff *bin_op(
3705         __isl_take isl_union_pw_multi_aff *upma1,
3706         __isl_take isl_union_pw_multi_aff *upma2,
3707         int (*fn)(void **entry, void *user))
3708 {
3709         isl_space *space;
3710         struct isl_union_pw_multi_aff_bin_data data = { NULL, NULL, NULL, fn };
3711
3712         space = isl_union_pw_multi_aff_get_space(upma2);
3713         upma1 = isl_union_pw_multi_aff_align_params(upma1, space);
3714         space = isl_union_pw_multi_aff_get_space(upma1);
3715         upma2 = isl_union_pw_multi_aff_align_params(upma2, space);
3716
3717         if (!upma1 || !upma2)
3718                 goto error;
3719
3720         data.upma2 = upma2;
3721         data.res = isl_union_pw_multi_aff_alloc(isl_space_copy(upma1->dim),
3722                                        upma1->table.n);
3723         if (isl_hash_table_foreach(upma1->dim->ctx, &upma1->table,
3724                                    &bin_entry, &data) < 0)
3725                 goto error;
3726
3727         isl_union_pw_multi_aff_free(upma1);
3728         isl_union_pw_multi_aff_free(upma2);
3729         return data.res;
3730 error:
3731         isl_union_pw_multi_aff_free(upma1);
3732         isl_union_pw_multi_aff_free(upma2);
3733         isl_union_pw_multi_aff_free(data.res);
3734         return NULL;
3735 }
3736
3737 /* Given two isl_multi_affs A -> B and C -> D,
3738  * construct an isl_multi_aff (A * C) -> (B, D).
3739  */
3740 __isl_give isl_multi_aff *isl_multi_aff_flat_range_product(
3741         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
3742 {
3743         int i, n1, n2;
3744         isl_aff *aff;
3745         isl_space *space;
3746         isl_multi_aff *res;
3747
3748         if (!ma1 || !ma2)
3749                 goto error;
3750
3751         space = isl_space_range_product(isl_multi_aff_get_space(ma1),
3752                                         isl_multi_aff_get_space(ma2));
3753         space = isl_space_flatten_range(space);
3754         res = isl_multi_aff_alloc(space);
3755
3756         n1 = isl_multi_aff_dim(ma1, isl_dim_out);
3757         n2 = isl_multi_aff_dim(ma2, isl_dim_out);
3758
3759         for (i = 0; i < n1; ++i) {
3760                 aff = isl_multi_aff_get_aff(ma1, i);
3761                 res = isl_multi_aff_set_aff(res, i, aff);
3762         }
3763
3764         for (i = 0; i < n2; ++i) {
3765                 aff = isl_multi_aff_get_aff(ma2, i);
3766                 res = isl_multi_aff_set_aff(res, n1 + i, aff);
3767         }
3768
3769         isl_multi_aff_free(ma1);
3770         isl_multi_aff_free(ma2);
3771         return res;
3772 error:
3773         isl_multi_aff_free(ma1);
3774         isl_multi_aff_free(ma2);
3775         return NULL;
3776 }
3777
3778 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
3779  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3780  */
3781 static __isl_give isl_pw_multi_aff *pw_multi_aff_flat_range_product(
3782         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3783 {
3784         isl_space *space;
3785
3786         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
3787                                         isl_pw_multi_aff_get_space(pma2));
3788         space = isl_space_flatten_range(space);
3789         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
3790                                             &isl_multi_aff_flat_range_product);
3791 }
3792
3793 /* Given two isl_pw_multi_affs A -> B and C -> D,
3794  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3795  */
3796 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_flat_range_product(
3797         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3798 {
3799         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3800                                             &pw_multi_aff_flat_range_product);
3801 }
3802
3803 /* If data->pma and *entry have the same domain space, then compute
3804  * their flat range product and the result to data->res.
3805  */
3806 static int flat_range_product_entry(void **entry, void *user)
3807 {
3808         struct isl_union_pw_multi_aff_bin_data *data = user;
3809         isl_pw_multi_aff *pma2 = *entry;
3810
3811         if (!isl_space_tuple_match(data->pma->dim, isl_dim_in,
3812                                  pma2->dim, isl_dim_in))
3813                 return 0;
3814
3815         pma2 = isl_pw_multi_aff_flat_range_product(
3816                                         isl_pw_multi_aff_copy(data->pma),
3817                                         isl_pw_multi_aff_copy(pma2));
3818
3819         data->res = isl_union_pw_multi_aff_add_pw_multi_aff(data->res, pma2);
3820
3821         return 0;
3822 }
3823
3824 /* Given two isl_union_pw_multi_affs A -> B and C -> D,
3825  * construct an isl_union_pw_multi_aff (A * C) -> (B, D).
3826  */
3827 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_flat_range_product(
3828         __isl_take isl_union_pw_multi_aff *upma1,
3829         __isl_take isl_union_pw_multi_aff *upma2)
3830 {
3831         return bin_op(upma1, upma2, &flat_range_product_entry);
3832 }
3833
3834 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3835  * The parameters are assumed to have been aligned.
3836  *
3837  * The implementation essentially performs an isl_pw_*_on_shared_domain,
3838  * except that it works on two different isl_pw_* types.
3839  */
3840 static __isl_give isl_pw_multi_aff *pw_multi_aff_set_pw_aff(
3841         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3842         __isl_take isl_pw_aff *pa)
3843 {
3844         int i, j, n;
3845         isl_pw_multi_aff *res = NULL;
3846
3847         if (!pma || !pa)
3848                 goto error;
3849
3850         if (!isl_space_tuple_match(pma->dim, isl_dim_in, pa->dim, isl_dim_in))
3851                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3852                         "domains don't match", goto error);
3853         if (pos >= isl_pw_multi_aff_dim(pma, isl_dim_out))
3854                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3855                         "index out of bounds", goto error);
3856
3857         n = pma->n * pa->n;
3858         res = isl_pw_multi_aff_alloc_size(isl_pw_multi_aff_get_space(pma), n);
3859
3860         for (i = 0; i < pma->n; ++i) {
3861                 for (j = 0; j < pa->n; ++j) {
3862                         isl_set *common;
3863                         isl_multi_aff *res_ij;
3864                         int empty;
3865
3866                         common = isl_set_intersect(isl_set_copy(pma->p[i].set),
3867                                                    isl_set_copy(pa->p[j].set));
3868                         empty = isl_set_plain_is_empty(common);
3869                         if (empty < 0 || empty) {
3870                                 isl_set_free(common);
3871                                 if (empty < 0)
3872                                         goto error;
3873                                 continue;
3874                         }
3875
3876                         res_ij = isl_multi_aff_set_aff(
3877                                         isl_multi_aff_copy(pma->p[i].maff), pos,
3878                                         isl_aff_copy(pa->p[j].aff));
3879                         res_ij = isl_multi_aff_gist(res_ij,
3880                                         isl_set_copy(common));
3881
3882                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3883                 }
3884         }
3885
3886         isl_pw_multi_aff_free(pma);
3887         isl_pw_aff_free(pa);
3888         return res;
3889 error:
3890         isl_pw_multi_aff_free(pma);
3891         isl_pw_aff_free(pa);
3892         return isl_pw_multi_aff_free(res);
3893 }
3894
3895 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3896  */
3897 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_set_pw_aff(
3898         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3899         __isl_take isl_pw_aff *pa)
3900 {
3901         if (!pma || !pa)
3902                 goto error;
3903         if (isl_space_match(pma->dim, isl_dim_param, pa->dim, isl_dim_param))
3904                 return pw_multi_aff_set_pw_aff(pma, pos, pa);
3905         if (!isl_space_has_named_params(pma->dim) ||
3906             !isl_space_has_named_params(pa->dim))
3907                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3908                         "unaligned unnamed parameters", goto error);
3909         pma = isl_pw_multi_aff_align_params(pma, isl_pw_aff_get_space(pa));
3910         pa = isl_pw_aff_align_params(pa, isl_pw_multi_aff_get_space(pma));
3911         return pw_multi_aff_set_pw_aff(pma, pos, pa);
3912 error:
3913         isl_pw_multi_aff_free(pma);
3914         isl_pw_aff_free(pa);
3915         return NULL;
3916 }