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 /* Mark the domains of "pwaff" as rational.
2513  */
2514 __isl_give isl_pw_aff *isl_pw_aff_set_rational(__isl_take isl_pw_aff *pwaff)
2515 {
2516         int i;
2517
2518         pwaff = isl_pw_aff_cow(pwaff);
2519         if (!pwaff)
2520                 return NULL;
2521         if (pwaff->n == 0)
2522                 return pwaff;
2523
2524         for (i = 0; i < pwaff->n; ++i) {
2525                 pwaff->p[i].set = isl_set_set_rational(pwaff->p[i].set);
2526                 if (!pwaff->p[i].set)
2527                         return isl_pw_aff_free(pwaff);
2528         }
2529
2530         return pwaff;
2531 }
2532
2533 /* Mark the domains of the elements of "list" as rational.
2534  */
2535 __isl_give isl_pw_aff_list *isl_pw_aff_list_set_rational(
2536         __isl_take isl_pw_aff_list *list)
2537 {
2538         int i;
2539
2540         if (!list)
2541                 return NULL;
2542         if (list->n == 0)
2543                 return list;
2544
2545         for (i = 0; i < list->n; ++i) {
2546                 isl_pw_aff *pa;
2547
2548                 pa = isl_pw_aff_list_get_pw_aff(list, i);
2549                 pa = isl_pw_aff_set_rational(pa);
2550                 list = isl_pw_aff_list_set_pw_aff(list, i, pa);
2551         }
2552
2553         return list;
2554 }
2555
2556 #undef BASE
2557 #define BASE aff
2558
2559 #include <isl_multi_templ.c>
2560
2561 /* Create an isl_pw_multi_aff with the given isl_multi_aff on a universe
2562  * domain.
2563  */
2564 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_multi_aff(
2565         __isl_take isl_multi_aff *ma)
2566 {
2567         isl_set *dom = isl_set_universe(isl_multi_aff_get_domain_space(ma));
2568         return isl_pw_multi_aff_alloc(dom, ma);
2569 }
2570
2571 /* Create a piecewise multi-affine expression in the given space that maps each
2572  * input dimension to the corresponding output dimension.
2573  */
2574 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_identity(
2575         __isl_take isl_space *space)
2576 {
2577         return isl_pw_multi_aff_from_multi_aff(isl_multi_aff_identity(space));
2578 }
2579
2580 __isl_give isl_multi_aff *isl_multi_aff_add(__isl_take isl_multi_aff *maff1,
2581         __isl_take isl_multi_aff *maff2)
2582 {
2583         int i;
2584         isl_ctx *ctx;
2585
2586         maff1 = isl_multi_aff_cow(maff1);
2587         if (!maff1 || !maff2)
2588                 goto error;
2589
2590         ctx = isl_multi_aff_get_ctx(maff1);
2591         if (!isl_space_is_equal(maff1->space, maff2->space))
2592                 isl_die(ctx, isl_error_invalid,
2593                         "spaces don't match", goto error);
2594
2595         for (i = 0; i < maff1->n; ++i) {
2596                 maff1->p[i] = isl_aff_add(maff1->p[i],
2597                                             isl_aff_copy(maff2->p[i]));
2598                 if (!maff1->p[i])
2599                         goto error;
2600         }
2601
2602         isl_multi_aff_free(maff2);
2603         return maff1;
2604 error:
2605         isl_multi_aff_free(maff1);
2606         isl_multi_aff_free(maff2);
2607         return NULL;
2608 }
2609
2610 /* Given two multi-affine expressions A -> B and C -> D,
2611  * construct a multi-affine expression [A -> C] -> [B -> D].
2612  */
2613 __isl_give isl_multi_aff *isl_multi_aff_product(
2614         __isl_take isl_multi_aff *ma1, __isl_take isl_multi_aff *ma2)
2615 {
2616         int i;
2617         isl_aff *aff;
2618         isl_space *space;
2619         isl_multi_aff *res;
2620         int in1, in2, out1, out2;
2621
2622         in1 = isl_multi_aff_dim(ma1, isl_dim_in);
2623         in2 = isl_multi_aff_dim(ma2, isl_dim_in);
2624         out1 = isl_multi_aff_dim(ma1, isl_dim_out);
2625         out2 = isl_multi_aff_dim(ma2, isl_dim_out);
2626         space = isl_space_product(isl_multi_aff_get_space(ma1),
2627                                   isl_multi_aff_get_space(ma2));
2628         res = isl_multi_aff_alloc(isl_space_copy(space));
2629         space = isl_space_domain(space);
2630
2631         for (i = 0; i < out1; ++i) {
2632                 aff = isl_multi_aff_get_aff(ma1, i);
2633                 aff = isl_aff_insert_dims(aff, isl_dim_in, in1, in2);
2634                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2635                 res = isl_multi_aff_set_aff(res, i, aff);
2636         }
2637
2638         for (i = 0; i < out2; ++i) {
2639                 aff = isl_multi_aff_get_aff(ma2, i);
2640                 aff = isl_aff_insert_dims(aff, isl_dim_in, 0, in1);
2641                 aff = isl_aff_reset_domain_space(aff, isl_space_copy(space));
2642                 res = isl_multi_aff_set_aff(res, out1 + i, aff);
2643         }
2644
2645         isl_space_free(space);
2646         isl_multi_aff_free(ma1);
2647         isl_multi_aff_free(ma2);
2648         return res;
2649 }
2650
2651 /* Exploit the equalities in "eq" to simplify the affine expressions.
2652  */
2653 static __isl_give isl_multi_aff *isl_multi_aff_substitute_equalities(
2654         __isl_take isl_multi_aff *maff, __isl_take isl_basic_set *eq)
2655 {
2656         int i;
2657
2658         maff = isl_multi_aff_cow(maff);
2659         if (!maff || !eq)
2660                 goto error;
2661
2662         for (i = 0; i < maff->n; ++i) {
2663                 maff->p[i] = isl_aff_substitute_equalities(maff->p[i],
2664                                                     isl_basic_set_copy(eq));
2665                 if (!maff->p[i])
2666                         goto error;
2667         }
2668
2669         isl_basic_set_free(eq);
2670         return maff;
2671 error:
2672         isl_basic_set_free(eq);
2673         isl_multi_aff_free(maff);
2674         return NULL;
2675 }
2676
2677 __isl_give isl_multi_aff *isl_multi_aff_scale(__isl_take isl_multi_aff *maff,
2678         isl_int f)
2679 {
2680         int i;
2681
2682         maff = isl_multi_aff_cow(maff);
2683         if (!maff)
2684                 return NULL;
2685
2686         for (i = 0; i < maff->n; ++i) {
2687                 maff->p[i] = isl_aff_scale(maff->p[i], f);
2688                 if (!maff->p[i])
2689                         return isl_multi_aff_free(maff);
2690         }
2691
2692         return maff;
2693 }
2694
2695 __isl_give isl_multi_aff *isl_multi_aff_add_on_domain(__isl_keep isl_set *dom,
2696         __isl_take isl_multi_aff *maff1, __isl_take isl_multi_aff *maff2)
2697 {
2698         maff1 = isl_multi_aff_add(maff1, maff2);
2699         maff1 = isl_multi_aff_gist(maff1, isl_set_copy(dom));
2700         return maff1;
2701 }
2702
2703 int isl_multi_aff_is_empty(__isl_keep isl_multi_aff *maff)
2704 {
2705         if (!maff)
2706                 return -1;
2707
2708         return 0;
2709 }
2710
2711 int isl_multi_aff_plain_is_equal(__isl_keep isl_multi_aff *maff1,
2712         __isl_keep isl_multi_aff *maff2)
2713 {
2714         int i;
2715         int equal;
2716
2717         if (!maff1 || !maff2)
2718                 return -1;
2719         if (maff1->n != maff2->n)
2720                 return 0;
2721         equal = isl_space_is_equal(maff1->space, maff2->space);
2722         if (equal < 0 || !equal)
2723                 return equal;
2724
2725         for (i = 0; i < maff1->n; ++i) {
2726                 equal = isl_aff_plain_is_equal(maff1->p[i], maff2->p[i]);
2727                 if (equal < 0 || !equal)
2728                         return equal;
2729         }
2730
2731         return 1;
2732 }
2733
2734 /* Return the set of domain elements where "ma1" is lexicographically
2735  * smaller than or equal to "ma2".
2736  */
2737 __isl_give isl_set *isl_multi_aff_lex_le_set(__isl_take isl_multi_aff *ma1,
2738         __isl_take isl_multi_aff *ma2)
2739 {
2740         return isl_multi_aff_lex_ge_set(ma2, ma1);
2741 }
2742
2743 /* Return the set of domain elements where "ma1" is lexicographically
2744  * greater than or equal to "ma2".
2745  */
2746 __isl_give isl_set *isl_multi_aff_lex_ge_set(__isl_take isl_multi_aff *ma1,
2747         __isl_take isl_multi_aff *ma2)
2748 {
2749         isl_space *space;
2750         isl_map *map1, *map2;
2751         isl_map *map, *ge;
2752
2753         map1 = isl_map_from_multi_aff(ma1);
2754         map2 = isl_map_from_multi_aff(ma2);
2755         map = isl_map_range_product(map1, map2);
2756         space = isl_space_range(isl_map_get_space(map));
2757         space = isl_space_domain(isl_space_unwrap(space));
2758         ge = isl_map_lex_ge(space);
2759         map = isl_map_intersect_range(map, isl_map_wrap(ge));
2760
2761         return isl_map_domain(map);
2762 }
2763
2764 #undef PW
2765 #define PW isl_pw_multi_aff
2766 #undef EL
2767 #define EL isl_multi_aff
2768 #undef EL_IS_ZERO
2769 #define EL_IS_ZERO is_empty
2770 #undef ZERO
2771 #define ZERO empty
2772 #undef IS_ZERO
2773 #define IS_ZERO is_empty
2774 #undef FIELD
2775 #define FIELD maff
2776 #undef DEFAULT_IS_ZERO
2777 #define DEFAULT_IS_ZERO 0
2778
2779 #define NO_NEG
2780 #define NO_EVAL
2781 #define NO_OPT
2782 #define NO_INVOLVES_DIMS
2783 #define NO_MOVE_DIMS
2784 #define NO_INSERT_DIMS
2785 #define NO_LIFT
2786 #define NO_MORPH
2787
2788 #include <isl_pw_templ.c>
2789
2790 #undef UNION
2791 #define UNION isl_union_pw_multi_aff
2792 #undef PART
2793 #define PART isl_pw_multi_aff
2794 #undef PARTS
2795 #define PARTS pw_multi_aff
2796 #define ALIGN_DOMAIN
2797
2798 #define NO_EVAL
2799
2800 #include <isl_union_templ.c>
2801
2802 /* Given a function "cmp" that returns the set of elements where
2803  * "ma1" is "better" than "ma2", return the intersection of this
2804  * set with "dom1" and "dom2".
2805  */
2806 static __isl_give isl_set *shared_and_better(__isl_keep isl_set *dom1,
2807         __isl_keep isl_set *dom2, __isl_keep isl_multi_aff *ma1,
2808         __isl_keep isl_multi_aff *ma2,
2809         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2810                                     __isl_take isl_multi_aff *ma2))
2811 {
2812         isl_set *common;
2813         isl_set *better;
2814         int is_empty;
2815
2816         common = isl_set_intersect(isl_set_copy(dom1), isl_set_copy(dom2));
2817         is_empty = isl_set_plain_is_empty(common);
2818         if (is_empty >= 0 && is_empty)
2819                 return common;
2820         if (is_empty < 0)
2821                 return isl_set_free(common);
2822         better = cmp(isl_multi_aff_copy(ma1), isl_multi_aff_copy(ma2));
2823         better = isl_set_intersect(common, better);
2824
2825         return better;
2826 }
2827
2828 /* Given a function "cmp" that returns the set of elements where
2829  * "ma1" is "better" than "ma2", return a piecewise multi affine
2830  * expression defined on the union of the definition domains
2831  * of "pma1" and "pma2" that maps to the "best" of "pma1" and
2832  * "pma2" on each cell.  If only one of the two input functions
2833  * is defined on a given cell, then it is considered the best.
2834  */
2835 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_opt(
2836         __isl_take isl_pw_multi_aff *pma1,
2837         __isl_take isl_pw_multi_aff *pma2,
2838         __isl_give isl_set *(*cmp)(__isl_take isl_multi_aff *ma1,
2839                                     __isl_take isl_multi_aff *ma2))
2840 {
2841         int i, j, n;
2842         isl_pw_multi_aff *res = NULL;
2843         isl_ctx *ctx;
2844         isl_set *set = NULL;
2845
2846         if (!pma1 || !pma2)
2847                 goto error;
2848
2849         ctx = isl_space_get_ctx(pma1->dim);
2850         if (!isl_space_is_equal(pma1->dim, pma2->dim))
2851                 isl_die(ctx, isl_error_invalid,
2852                         "arguments should live in the same space", goto error);
2853
2854         if (isl_pw_multi_aff_is_empty(pma1)) {
2855                 isl_pw_multi_aff_free(pma1);
2856                 return pma2;
2857         }
2858
2859         if (isl_pw_multi_aff_is_empty(pma2)) {
2860                 isl_pw_multi_aff_free(pma2);
2861                 return pma1;
2862         }
2863
2864         n = 2 * (pma1->n + 1) * (pma2->n + 1);
2865         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma1->dim), n);
2866
2867         for (i = 0; i < pma1->n; ++i) {
2868                 set = isl_set_copy(pma1->p[i].set);
2869                 for (j = 0; j < pma2->n; ++j) {
2870                         isl_set *better;
2871                         int is_empty;
2872
2873                         better = shared_and_better(pma2->p[j].set,
2874                                         pma1->p[i].set, pma2->p[j].maff,
2875                                         pma1->p[i].maff, cmp);
2876                         is_empty = isl_set_plain_is_empty(better);
2877                         if (is_empty < 0 || is_empty) {
2878                                 isl_set_free(better);
2879                                 if (is_empty < 0)
2880                                         goto error;
2881                                 continue;
2882                         }
2883                         set = isl_set_subtract(set, isl_set_copy(better));
2884
2885                         res = isl_pw_multi_aff_add_piece(res, better,
2886                                         isl_multi_aff_copy(pma2->p[j].maff));
2887                 }
2888                 res = isl_pw_multi_aff_add_piece(res, set,
2889                                         isl_multi_aff_copy(pma1->p[i].maff));
2890         }
2891
2892         for (j = 0; j < pma2->n; ++j) {
2893                 set = isl_set_copy(pma2->p[j].set);
2894                 for (i = 0; i < pma1->n; ++i)
2895                         set = isl_set_subtract(set,
2896                                         isl_set_copy(pma1->p[i].set));
2897                 res = isl_pw_multi_aff_add_piece(res, set,
2898                                         isl_multi_aff_copy(pma2->p[j].maff));
2899         }
2900
2901         isl_pw_multi_aff_free(pma1);
2902         isl_pw_multi_aff_free(pma2);
2903
2904         return res;
2905 error:
2906         isl_pw_multi_aff_free(pma1);
2907         isl_pw_multi_aff_free(pma2);
2908         isl_set_free(set);
2909         return isl_pw_multi_aff_free(res);
2910 }
2911
2912 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmax(
2913         __isl_take isl_pw_multi_aff *pma1,
2914         __isl_take isl_pw_multi_aff *pma2)
2915 {
2916         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_ge_set);
2917 }
2918
2919 /* Given two piecewise multi affine expressions, return a piecewise
2920  * multi-affine expression defined on the union of the definition domains
2921  * of the inputs that is equal to the lexicographic maximum of the two
2922  * inputs on each cell.  If only one of the two inputs is defined on
2923  * a given cell, then it is considered to be the maximum.
2924  */
2925 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmax(
2926         __isl_take isl_pw_multi_aff *pma1,
2927         __isl_take isl_pw_multi_aff *pma2)
2928 {
2929         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2930                                                     &pw_multi_aff_union_lexmax);
2931 }
2932
2933 static __isl_give isl_pw_multi_aff *pw_multi_aff_union_lexmin(
2934         __isl_take isl_pw_multi_aff *pma1,
2935         __isl_take isl_pw_multi_aff *pma2)
2936 {
2937         return pw_multi_aff_union_opt(pma1, pma2, &isl_multi_aff_lex_le_set);
2938 }
2939
2940 /* Given two piecewise multi affine expressions, return a piecewise
2941  * multi-affine expression defined on the union of the definition domains
2942  * of the inputs that is equal to the lexicographic minimum of the two
2943  * inputs on each cell.  If only one of the two inputs is defined on
2944  * a given cell, then it is considered to be the minimum.
2945  */
2946 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_lexmin(
2947         __isl_take isl_pw_multi_aff *pma1,
2948         __isl_take isl_pw_multi_aff *pma2)
2949 {
2950         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2951                                                     &pw_multi_aff_union_lexmin);
2952 }
2953
2954 static __isl_give isl_pw_multi_aff *pw_multi_aff_add(
2955         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2956 {
2957         return isl_pw_multi_aff_on_shared_domain(pma1, pma2,
2958                                                 &isl_multi_aff_add);
2959 }
2960
2961 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_add(
2962         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2963 {
2964         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
2965                                                 &pw_multi_aff_add);
2966 }
2967
2968 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_union_add(
2969         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2970 {
2971         return isl_pw_multi_aff_union_add_(pma1, pma2);
2972 }
2973
2974 /* Given two piecewise multi-affine expressions A -> B and C -> D,
2975  * construct a piecewise multi-affine expression [A -> C] -> [B -> D].
2976  */
2977 static __isl_give isl_pw_multi_aff *pw_multi_aff_product(
2978         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
2979 {
2980         int i, j, n;
2981         isl_space *space;
2982         isl_pw_multi_aff *res;
2983
2984         if (!pma1 || !pma2)
2985                 goto error;
2986
2987         n = pma1->n * pma2->n;
2988         space = isl_space_product(isl_space_copy(pma1->dim),
2989                                   isl_space_copy(pma2->dim));
2990         res = isl_pw_multi_aff_alloc_size(space, n);
2991
2992         for (i = 0; i < pma1->n; ++i) {
2993                 for (j = 0; j < pma2->n; ++j) {
2994                         isl_set *domain;
2995                         isl_multi_aff *ma;
2996
2997                         domain = isl_set_product(isl_set_copy(pma1->p[i].set),
2998                                                  isl_set_copy(pma2->p[j].set));
2999                         ma = isl_multi_aff_product(
3000                                         isl_multi_aff_copy(pma1->p[i].maff),
3001                                         isl_multi_aff_copy(pma2->p[i].maff));
3002                         res = isl_pw_multi_aff_add_piece(res, domain, ma);
3003                 }
3004         }
3005
3006         isl_pw_multi_aff_free(pma1);
3007         isl_pw_multi_aff_free(pma2);
3008         return res;
3009 error:
3010         isl_pw_multi_aff_free(pma1);
3011         isl_pw_multi_aff_free(pma2);
3012         return NULL;
3013 }
3014
3015 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_product(
3016         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3017 {
3018         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3019                                                 &pw_multi_aff_product);
3020 }
3021
3022 /* Construct a map mapping the domain of the piecewise multi-affine expression
3023  * to its range, with each dimension in the range equated to the
3024  * corresponding affine expression on its cell.
3025  */
3026 __isl_give isl_map *isl_map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3027 {
3028         int i;
3029         isl_map *map;
3030
3031         if (!pma)
3032                 return NULL;
3033
3034         map = isl_map_empty(isl_pw_multi_aff_get_space(pma));
3035
3036         for (i = 0; i < pma->n; ++i) {
3037                 isl_multi_aff *maff;
3038                 isl_basic_map *bmap;
3039                 isl_map *map_i;
3040
3041                 maff = isl_multi_aff_copy(pma->p[i].maff);
3042                 bmap = isl_basic_map_from_multi_aff(maff);
3043                 map_i = isl_map_from_basic_map(bmap);
3044                 map_i = isl_map_intersect_domain(map_i,
3045                                                 isl_set_copy(pma->p[i].set));
3046                 map = isl_map_union_disjoint(map, map_i);
3047         }
3048
3049         isl_pw_multi_aff_free(pma);
3050         return map;
3051 }
3052
3053 __isl_give isl_set *isl_set_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma)
3054 {
3055         if (!pma)
3056                 return NULL;
3057
3058         if (!isl_space_is_set(pma->dim))
3059                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3060                         "isl_pw_multi_aff cannot be converted into an isl_set",
3061                         return isl_pw_multi_aff_free(pma));
3062
3063         return isl_map_from_pw_multi_aff(pma);
3064 }
3065
3066 /* Given a basic map with a single output dimension that is defined
3067  * in terms of the parameters and input dimensions using an equality,
3068  * extract an isl_aff that expresses the output dimension in terms
3069  * of the parameters and input dimensions.
3070  *
3071  * Since some applications expect the result of isl_pw_multi_aff_from_map
3072  * to only contain integer affine expressions, we compute the floor
3073  * of the expression before returning.
3074  *
3075  * This function shares some similarities with
3076  * isl_basic_map_has_defining_equality and isl_constraint_get_bound.
3077  */
3078 static __isl_give isl_aff *extract_isl_aff_from_basic_map(
3079         __isl_take isl_basic_map *bmap)
3080 {
3081         int i;
3082         unsigned offset;
3083         unsigned total;
3084         isl_local_space *ls;
3085         isl_aff *aff;
3086
3087         if (!bmap)
3088                 return NULL;
3089         if (isl_basic_map_dim(bmap, isl_dim_out) != 1)
3090                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3091                         "basic map should have a single output dimension",
3092                         goto error);
3093         offset = isl_basic_map_offset(bmap, isl_dim_out);
3094         total = isl_basic_map_total_dim(bmap);
3095         for (i = 0; i < bmap->n_eq; ++i) {
3096                 if (isl_int_is_zero(bmap->eq[i][offset]))
3097                         continue;
3098                 if (isl_seq_first_non_zero(bmap->eq[i] + offset + 1,
3099                                            1 + total - (offset + 1)) != -1)
3100                         continue;
3101                 break;
3102         }
3103         if (i >= bmap->n_eq)
3104                 isl_die(isl_basic_map_get_ctx(bmap), isl_error_invalid,
3105                         "unable to find suitable equality", goto error);
3106         ls = isl_basic_map_get_local_space(bmap);
3107         aff = isl_aff_alloc(isl_local_space_domain(ls));
3108         if (!aff)
3109                 goto error;
3110         if (isl_int_is_neg(bmap->eq[i][offset]))
3111                 isl_seq_cpy(aff->v->el + 1, bmap->eq[i], offset);
3112         else
3113                 isl_seq_neg(aff->v->el + 1, bmap->eq[i], offset);
3114         isl_seq_clr(aff->v->el + 1 + offset, aff->v->size - (1 + offset));
3115         isl_int_abs(aff->v->el[0], bmap->eq[i][offset]);
3116         isl_basic_map_free(bmap);
3117
3118         aff = isl_aff_remove_unused_divs(aff);
3119         aff = isl_aff_floor(aff);
3120         return aff;
3121 error:
3122         isl_basic_map_free(bmap);
3123         return NULL;
3124 }
3125
3126 /* Given a basic map where each output dimension is defined
3127  * in terms of the parameters and input dimensions using an equality,
3128  * extract an isl_multi_aff that expresses the output dimensions in terms
3129  * of the parameters and input dimensions.
3130  */
3131 static __isl_give isl_multi_aff *extract_isl_multi_aff_from_basic_map(
3132         __isl_take isl_basic_map *bmap)
3133 {
3134         int i;
3135         unsigned n_out;
3136         isl_multi_aff *ma;
3137
3138         if (!bmap)
3139                 return NULL;
3140
3141         ma = isl_multi_aff_alloc(isl_basic_map_get_space(bmap));
3142         n_out = isl_basic_map_dim(bmap, isl_dim_out);
3143
3144         for (i = 0; i < n_out; ++i) {
3145                 isl_basic_map *bmap_i;
3146                 isl_aff *aff;
3147
3148                 bmap_i = isl_basic_map_copy(bmap);
3149                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out,
3150                                                         i + 1, n_out - (1 + i));
3151                 bmap_i = isl_basic_map_project_out(bmap_i, isl_dim_out, 0, i);
3152                 aff = extract_isl_aff_from_basic_map(bmap_i);
3153                 ma = isl_multi_aff_set_aff(ma, i, aff);
3154         }
3155
3156         isl_basic_map_free(bmap);
3157
3158         return ma;
3159 }
3160
3161 /* Create an isl_pw_multi_aff that is equivalent to
3162  * isl_map_intersect_domain(isl_map_from_basic_map(bmap), domain).
3163  * The given basic map is such that each output dimension is defined
3164  * in terms of the parameters and input dimensions using an equality.
3165  */
3166 static __isl_give isl_pw_multi_aff *plain_pw_multi_aff_from_map(
3167         __isl_take isl_set *domain, __isl_take isl_basic_map *bmap)
3168 {
3169         isl_multi_aff *ma;
3170
3171         ma = extract_isl_multi_aff_from_basic_map(bmap);
3172         return isl_pw_multi_aff_alloc(domain, ma);
3173 }
3174
3175 /* Try and create an isl_pw_multi_aff that is equivalent to the given isl_map.
3176  * This obviously only works if the input "map" is single-valued.
3177  * If so, we compute the lexicographic minimum of the image in the form
3178  * of an isl_pw_multi_aff.  Since the image is unique, it is equal
3179  * to its lexicographic minimum.
3180  * If the input is not single-valued, we produce an error.
3181  *
3182  * As a special case, we first check if all output dimensions are uniquely
3183  * defined in terms of the parameters and input dimensions over the entire
3184  * domain.  If so, we extract the desired isl_pw_multi_aff directly
3185  * from the affine hull of "map" and its domain.
3186  */
3187 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_map(__isl_take isl_map *map)
3188 {
3189         int i;
3190         int sv;
3191         isl_pw_multi_aff *pma;
3192         isl_basic_map *hull;
3193
3194         if (!map)
3195                 return NULL;
3196
3197         hull = isl_map_affine_hull(isl_map_copy(map));
3198         sv = isl_basic_map_plain_is_single_valued(hull);
3199         if (sv >= 0 && sv)
3200                 return plain_pw_multi_aff_from_map(isl_map_domain(map), hull);
3201         isl_basic_map_free(hull);
3202         if (sv < 0)
3203                 goto error;
3204
3205         sv = isl_map_is_single_valued(map);
3206         if (sv < 0)
3207                 goto error;
3208         if (!sv)
3209                 isl_die(isl_map_get_ctx(map), isl_error_invalid,
3210                         "map is not single-valued", goto error);
3211         map = isl_map_make_disjoint(map);
3212         if (!map)
3213                 return NULL;
3214
3215         pma = isl_pw_multi_aff_empty(isl_map_get_space(map));
3216
3217         for (i = 0; i < map->n; ++i) {
3218                 isl_pw_multi_aff *pma_i;
3219                 isl_basic_map *bmap;
3220                 bmap = isl_basic_map_copy(map->p[i]);
3221                 pma_i = isl_basic_map_lexmin_pw_multi_aff(bmap);
3222                 pma = isl_pw_multi_aff_add_disjoint(pma, pma_i);
3223         }
3224
3225         isl_map_free(map);
3226         return pma;
3227 error:
3228         isl_map_free(map);
3229         return NULL;
3230 }
3231
3232 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_set(__isl_take isl_set *set)
3233 {
3234         return isl_pw_multi_aff_from_map(set);
3235 }
3236
3237 /* Return the piecewise affine expression "set ? 1 : 0".
3238  */
3239 __isl_give isl_pw_aff *isl_set_indicator_function(__isl_take isl_set *set)
3240 {
3241         isl_pw_aff *pa;
3242         isl_space *space = isl_set_get_space(set);
3243         isl_local_space *ls = isl_local_space_from_space(space);
3244         isl_aff *zero = isl_aff_zero_on_domain(isl_local_space_copy(ls));
3245         isl_aff *one = isl_aff_zero_on_domain(ls);
3246
3247         one = isl_aff_add_constant_si(one, 1);
3248         pa = isl_pw_aff_alloc(isl_set_copy(set), one);
3249         set = isl_set_complement(set);
3250         pa = isl_pw_aff_add_disjoint(pa, isl_pw_aff_alloc(set, zero));
3251
3252         return pa;
3253 }
3254
3255 /* Plug in "subs" for dimension "type", "pos" of "aff".
3256  *
3257  * Let i be the dimension to replace and let "subs" be of the form
3258  *
3259  *      f/d
3260  *
3261  * and "aff" of the form
3262  *
3263  *      (a i + g)/m
3264  *
3265  * The result is
3266  *
3267  *      (a f + d g')/(m d)
3268  *
3269  * where g' is the result of plugging in "subs" in each of the integer
3270  * divisions in g.
3271  */
3272 __isl_give isl_aff *isl_aff_substitute(__isl_take isl_aff *aff,
3273         enum isl_dim_type type, unsigned pos, __isl_keep isl_aff *subs)
3274 {
3275         isl_ctx *ctx;
3276         isl_int v;
3277
3278         aff = isl_aff_cow(aff);
3279         if (!aff || !subs)
3280                 return isl_aff_free(aff);
3281
3282         ctx = isl_aff_get_ctx(aff);
3283         if (!isl_space_is_equal(aff->ls->dim, subs->ls->dim))
3284                 isl_die(ctx, isl_error_invalid,
3285                         "spaces don't match", return isl_aff_free(aff));
3286         if (isl_local_space_dim(subs->ls, isl_dim_div) != 0)
3287                 isl_die(ctx, isl_error_unsupported,
3288                         "cannot handle divs yet", return isl_aff_free(aff));
3289
3290         aff->ls = isl_local_space_substitute(aff->ls, type, pos, subs);
3291         if (!aff->ls)
3292                 return isl_aff_free(aff);
3293
3294         aff->v = isl_vec_cow(aff->v);
3295         if (!aff->v)
3296                 return isl_aff_free(aff);
3297
3298         pos += isl_local_space_offset(aff->ls, type);
3299
3300         isl_int_init(v);
3301         isl_seq_substitute(aff->v->el, pos, subs->v->el,
3302                             aff->v->size, subs->v->size, v);
3303         isl_int_clear(v);
3304
3305         return aff;
3306 }
3307
3308 /* Plug in "subs" for dimension "type", "pos" in each of the affine
3309  * expressions in "maff".
3310  */
3311 __isl_give isl_multi_aff *isl_multi_aff_substitute(
3312         __isl_take isl_multi_aff *maff, enum isl_dim_type type, unsigned pos,
3313         __isl_keep isl_aff *subs)
3314 {
3315         int i;
3316
3317         maff = isl_multi_aff_cow(maff);
3318         if (!maff || !subs)
3319                 return isl_multi_aff_free(maff);
3320
3321         if (type == isl_dim_in)
3322                 type = isl_dim_set;
3323
3324         for (i = 0; i < maff->n; ++i) {
3325                 maff->p[i] = isl_aff_substitute(maff->p[i], type, pos, subs);
3326                 if (!maff->p[i])
3327                         return isl_multi_aff_free(maff);
3328         }
3329
3330         return maff;
3331 }
3332
3333 /* Plug in "subs" for dimension "type", "pos" of "pma".
3334  *
3335  * pma is of the form
3336  *
3337  *      A_i(v) -> M_i(v)
3338  *
3339  * while subs is of the form
3340  *
3341  *      v' = B_j(v) -> S_j
3342  *
3343  * Each pair i,j such that C_ij = A_i \cap B_i is non-empty
3344  * has a contribution in the result, in particular
3345  *
3346  *      C_ij(S_j) -> M_i(S_j)
3347  *
3348  * Note that plugging in S_j in C_ij may also result in an empty set
3349  * and this contribution should simply be discarded.
3350  */
3351 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_substitute(
3352         __isl_take isl_pw_multi_aff *pma, enum isl_dim_type type, unsigned pos,
3353         __isl_keep isl_pw_aff *subs)
3354 {
3355         int i, j, n;
3356         isl_pw_multi_aff *res;
3357
3358         if (!pma || !subs)
3359                 return isl_pw_multi_aff_free(pma);
3360
3361         n = pma->n * subs->n;
3362         res = isl_pw_multi_aff_alloc_size(isl_space_copy(pma->dim), n);
3363
3364         for (i = 0; i < pma->n; ++i) {
3365                 for (j = 0; j < subs->n; ++j) {
3366                         isl_set *common;
3367                         isl_multi_aff *res_ij;
3368                         common = isl_set_intersect(
3369                                         isl_set_copy(pma->p[i].set),
3370                                         isl_set_copy(subs->p[j].set));
3371                         common = isl_set_substitute(common,
3372                                         type, pos, subs->p[j].aff);
3373                         if (isl_set_plain_is_empty(common)) {
3374                                 isl_set_free(common);
3375                                 continue;
3376                         }
3377
3378                         res_ij = isl_multi_aff_substitute(
3379                                         isl_multi_aff_copy(pma->p[i].maff),
3380                                         type, pos, subs->p[j].aff);
3381
3382                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3383                 }
3384         }
3385
3386         isl_pw_multi_aff_free(pma);
3387         return res;
3388 }
3389
3390 /* Extend the local space of "dst" to include the divs
3391  * in the local space of "src".
3392  */
3393 __isl_give isl_aff *isl_aff_align_divs(__isl_take isl_aff *dst,
3394         __isl_keep isl_aff *src)
3395 {
3396         isl_ctx *ctx;
3397         int *exp1 = NULL;
3398         int *exp2 = NULL;
3399         isl_mat *div;
3400
3401         if (!src || !dst)
3402                 return isl_aff_free(dst);
3403
3404         ctx = isl_aff_get_ctx(src);
3405         if (!isl_space_is_equal(src->ls->dim, dst->ls->dim))
3406                 isl_die(ctx, isl_error_invalid,
3407                         "spaces don't match", goto error);
3408
3409         if (src->ls->div->n_row == 0)
3410                 return dst;
3411
3412         exp1 = isl_alloc_array(ctx, int, src->ls->div->n_row);
3413         exp2 = isl_alloc_array(ctx, int, dst->ls->div->n_row);
3414         if (!exp1 || !exp2)
3415                 goto error;
3416
3417         div = isl_merge_divs(src->ls->div, dst->ls->div, exp1, exp2);
3418         dst = isl_aff_expand_divs(dst, div, exp2);
3419         free(exp1);
3420         free(exp2);
3421
3422         return dst;
3423 error:
3424         free(exp1);
3425         free(exp2);
3426         return isl_aff_free(dst);
3427 }
3428
3429 /* Adjust the local spaces of the affine expressions in "maff"
3430  * such that they all have the save divs.
3431  */
3432 __isl_give isl_multi_aff *isl_multi_aff_align_divs(
3433         __isl_take isl_multi_aff *maff)
3434 {
3435         int i;
3436
3437         if (!maff)
3438                 return NULL;
3439         if (maff->n == 0)
3440                 return maff;
3441         maff = isl_multi_aff_cow(maff);
3442         if (!maff)
3443                 return NULL;
3444
3445         for (i = 1; i < maff->n; ++i)
3446                 maff->p[0] = isl_aff_align_divs(maff->p[0], maff->p[i]);
3447         for (i = 1; i < maff->n; ++i) {
3448                 maff->p[i] = isl_aff_align_divs(maff->p[i], maff->p[0]);
3449                 if (!maff->p[i])
3450                         return isl_multi_aff_free(maff);
3451         }
3452
3453         return maff;
3454 }
3455
3456 __isl_give isl_aff *isl_aff_lift(__isl_take isl_aff *aff)
3457 {
3458         aff = isl_aff_cow(aff);
3459         if (!aff)
3460                 return NULL;
3461
3462         aff->ls = isl_local_space_lift(aff->ls);
3463         if (!aff->ls)
3464                 return isl_aff_free(aff);
3465
3466         return aff;
3467 }
3468
3469 /* Lift "maff" to a space with extra dimensions such that the result
3470  * has no more existentially quantified variables.
3471  * If "ls" is not NULL, then *ls is assigned the local space that lies
3472  * at the basis of the lifting applied to "maff".
3473  */
3474 __isl_give isl_multi_aff *isl_multi_aff_lift(__isl_take isl_multi_aff *maff,
3475         __isl_give isl_local_space **ls)
3476 {
3477         int i;
3478         isl_space *space;
3479         unsigned n_div;
3480
3481         if (ls)
3482                 *ls = NULL;
3483
3484         if (!maff)
3485                 return NULL;
3486
3487         if (maff->n == 0) {
3488                 if (ls) {
3489                         isl_space *space = isl_multi_aff_get_domain_space(maff);
3490                         *ls = isl_local_space_from_space(space);
3491                         if (!*ls)
3492                                 return isl_multi_aff_free(maff);
3493                 }
3494                 return maff;
3495         }
3496
3497         maff = isl_multi_aff_cow(maff);
3498         maff = isl_multi_aff_align_divs(maff);
3499         if (!maff)
3500                 return NULL;
3501
3502         n_div = isl_aff_dim(maff->p[0], isl_dim_div);
3503         space = isl_multi_aff_get_space(maff);
3504         space = isl_space_lift(isl_space_domain(space), n_div);
3505         space = isl_space_extend_domain_with_range(space,
3506                                                 isl_multi_aff_get_space(maff));
3507         if (!space)
3508                 return isl_multi_aff_free(maff);
3509         isl_space_free(maff->space);
3510         maff->space = space;
3511
3512         if (ls) {
3513                 *ls = isl_aff_get_domain_local_space(maff->p[0]);
3514                 if (!*ls)
3515                         return isl_multi_aff_free(maff);
3516         }
3517
3518         for (i = 0; i < maff->n; ++i) {
3519                 maff->p[i] = isl_aff_lift(maff->p[i]);
3520                 if (!maff->p[i])
3521                         goto error;
3522         }
3523
3524         return maff;
3525 error:
3526         if (ls)
3527                 isl_local_space_free(*ls);
3528         return isl_multi_aff_free(maff);
3529 }
3530
3531
3532 /* Extract an isl_pw_aff corresponding to output dimension "pos" of "pma".
3533  */
3534 __isl_give isl_pw_aff *isl_pw_multi_aff_get_pw_aff(
3535         __isl_keep isl_pw_multi_aff *pma, int pos)
3536 {
3537         int i;
3538         int n_out;
3539         isl_space *space;
3540         isl_pw_aff *pa;
3541
3542         if (!pma)
3543                 return NULL;
3544
3545         n_out = isl_pw_multi_aff_dim(pma, isl_dim_out);
3546         if (pos < 0 || pos >= n_out)
3547                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3548                         "index out of bounds", return NULL);
3549
3550         space = isl_pw_multi_aff_get_space(pma);
3551         space = isl_space_drop_dims(space, isl_dim_out,
3552                                     pos + 1, n_out - pos - 1);
3553         space = isl_space_drop_dims(space, isl_dim_out, 0, pos);
3554
3555         pa = isl_pw_aff_alloc_size(space, pma->n);
3556         for (i = 0; i < pma->n; ++i) {
3557                 isl_aff *aff;
3558                 aff = isl_multi_aff_get_aff(pma->p[i].maff, pos);
3559                 pa = isl_pw_aff_add_piece(pa, isl_set_copy(pma->p[i].set), aff);
3560         }
3561
3562         return pa;
3563 }
3564
3565 /* Return an isl_pw_multi_aff with the given "set" as domain and
3566  * an unnamed zero-dimensional range.
3567  */
3568 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_from_domain(
3569         __isl_take isl_set *set)
3570 {
3571         isl_multi_aff *ma;
3572         isl_space *space;
3573
3574         space = isl_set_get_space(set);
3575         space = isl_space_from_domain(space);
3576         ma = isl_multi_aff_zero(space);
3577         return isl_pw_multi_aff_alloc(set, ma);
3578 }
3579
3580 /* Add an isl_pw_multi_aff with the given "set" as domain and
3581  * an unnamed zero-dimensional range to *user.
3582  */
3583 static int add_pw_multi_aff_from_domain(__isl_take isl_set *set, void *user)
3584 {
3585         isl_union_pw_multi_aff **upma = user;
3586         isl_pw_multi_aff *pma;
3587
3588         pma = isl_pw_multi_aff_from_domain(set);
3589         *upma = isl_union_pw_multi_aff_add_pw_multi_aff(*upma, pma);
3590
3591         return 0;
3592 }
3593
3594 /* Return an isl_union_pw_multi_aff with the given "uset" as domain and
3595  * an unnamed zero-dimensional range.
3596  */
3597 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_from_domain(
3598         __isl_take isl_union_set *uset)
3599 {
3600         isl_space *space;
3601         isl_union_pw_multi_aff *upma;
3602
3603         if (!uset)
3604                 return NULL;
3605
3606         space = isl_union_set_get_space(uset);
3607         upma = isl_union_pw_multi_aff_empty(space);
3608
3609         if (isl_union_set_foreach_set(uset,
3610                                     &add_pw_multi_aff_from_domain, &upma) < 0)
3611                 goto error;
3612
3613         isl_union_set_free(uset);
3614         return upma;
3615 error:
3616         isl_union_set_free(uset);
3617         isl_union_pw_multi_aff_free(upma);
3618         return NULL;
3619 }
3620
3621 /* Convert "pma" to an isl_map and add it to *umap.
3622  */
3623 static int map_from_pw_multi_aff(__isl_take isl_pw_multi_aff *pma, void *user)
3624 {
3625         isl_union_map **umap = user;
3626         isl_map *map;
3627
3628         map = isl_map_from_pw_multi_aff(pma);
3629         *umap = isl_union_map_add_map(*umap, map);
3630
3631         return 0;
3632 }
3633
3634 /* Construct a union map mapping the domain of the union
3635  * piecewise multi-affine expression to its range, with each dimension
3636  * in the range equated to the corresponding affine expression on its cell.
3637  */
3638 __isl_give isl_union_map *isl_union_map_from_union_pw_multi_aff(
3639         __isl_take isl_union_pw_multi_aff *upma)
3640 {
3641         isl_space *space;
3642         isl_union_map *umap;
3643
3644         if (!upma)
3645                 return NULL;
3646
3647         space = isl_union_pw_multi_aff_get_space(upma);
3648         umap = isl_union_map_empty(space);
3649
3650         if (isl_union_pw_multi_aff_foreach_pw_multi_aff(upma,
3651                                         &map_from_pw_multi_aff, &umap) < 0)
3652                 goto error;
3653
3654         isl_union_pw_multi_aff_free(upma);
3655         return umap;
3656 error:
3657         isl_union_pw_multi_aff_free(upma);
3658         isl_union_map_free(umap);
3659         return NULL;
3660 }
3661
3662 /* Local data for bin_entry and the callback "fn".
3663  */
3664 struct isl_union_pw_multi_aff_bin_data {
3665         isl_union_pw_multi_aff *upma2;
3666         isl_union_pw_multi_aff *res;
3667         isl_pw_multi_aff *pma;
3668         int (*fn)(void **entry, void *user);
3669 };
3670
3671 /* Given an isl_pw_multi_aff from upma1, store it in data->pma
3672  * and call data->fn for each isl_pw_multi_aff in data->upma2.
3673  */
3674 static int bin_entry(void **entry, void *user)
3675 {
3676         struct isl_union_pw_multi_aff_bin_data *data = user;
3677         isl_pw_multi_aff *pma = *entry;
3678
3679         data->pma = pma;
3680         if (isl_hash_table_foreach(data->upma2->dim->ctx, &data->upma2->table,
3681                                    data->fn, data) < 0)
3682                 return -1;
3683
3684         return 0;
3685 }
3686
3687 /* Call "fn" on each pair of isl_pw_multi_affs in "upma1" and "upma2".
3688  * The isl_pw_multi_aff from upma1 is stored in data->pma (where data is
3689  * passed as user field) and the isl_pw_multi_aff from upma2 is available
3690  * as *entry.  The callback should adjust data->res if desired.
3691  */
3692 static __isl_give isl_union_pw_multi_aff *bin_op(
3693         __isl_take isl_union_pw_multi_aff *upma1,
3694         __isl_take isl_union_pw_multi_aff *upma2,
3695         int (*fn)(void **entry, void *user))
3696 {
3697         isl_space *space;
3698         struct isl_union_pw_multi_aff_bin_data data = { NULL, NULL, NULL, fn };
3699
3700         space = isl_union_pw_multi_aff_get_space(upma2);
3701         upma1 = isl_union_pw_multi_aff_align_params(upma1, space);
3702         space = isl_union_pw_multi_aff_get_space(upma1);
3703         upma2 = isl_union_pw_multi_aff_align_params(upma2, space);
3704
3705         if (!upma1 || !upma2)
3706                 goto error;
3707
3708         data.upma2 = upma2;
3709         data.res = isl_union_pw_multi_aff_alloc(isl_space_copy(upma1->dim),
3710                                        upma1->table.n);
3711         if (isl_hash_table_foreach(upma1->dim->ctx, &upma1->table,
3712                                    &bin_entry, &data) < 0)
3713                 goto error;
3714
3715         isl_union_pw_multi_aff_free(upma1);
3716         isl_union_pw_multi_aff_free(upma2);
3717         return data.res;
3718 error:
3719         isl_union_pw_multi_aff_free(upma1);
3720         isl_union_pw_multi_aff_free(upma2);
3721         isl_union_pw_multi_aff_free(data.res);
3722         return NULL;
3723 }
3724
3725 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
3726  * construct an isl_pw_multi_aff (A * C) -> [B -> D].
3727  */
3728 static __isl_give isl_pw_multi_aff *pw_multi_aff_range_product(
3729         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3730 {
3731         isl_space *space;
3732
3733         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
3734                                         isl_pw_multi_aff_get_space(pma2));
3735         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
3736                                             &isl_multi_aff_range_product);
3737 }
3738
3739 /* Given two isl_pw_multi_affs A -> B and C -> D,
3740  * construct an isl_pw_multi_aff (A * C) -> [B -> D].
3741  */
3742 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_range_product(
3743         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3744 {
3745         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3746                                             &pw_multi_aff_range_product);
3747 }
3748
3749 /* Given two aligned isl_pw_multi_affs A -> B and C -> D,
3750  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3751  */
3752 static __isl_give isl_pw_multi_aff *pw_multi_aff_flat_range_product(
3753         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3754 {
3755         isl_space *space;
3756
3757         space = isl_space_range_product(isl_pw_multi_aff_get_space(pma1),
3758                                         isl_pw_multi_aff_get_space(pma2));
3759         space = isl_space_flatten_range(space);
3760         return isl_pw_multi_aff_on_shared_domain_in(pma1, pma2, space,
3761                                             &isl_multi_aff_flat_range_product);
3762 }
3763
3764 /* Given two isl_pw_multi_affs A -> B and C -> D,
3765  * construct an isl_pw_multi_aff (A * C) -> (B, D).
3766  */
3767 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_flat_range_product(
3768         __isl_take isl_pw_multi_aff *pma1, __isl_take isl_pw_multi_aff *pma2)
3769 {
3770         return isl_pw_multi_aff_align_params_pw_pw_and(pma1, pma2,
3771                                             &pw_multi_aff_flat_range_product);
3772 }
3773
3774 /* If data->pma and *entry have the same domain space, then compute
3775  * their flat range product and the result to data->res.
3776  */
3777 static int flat_range_product_entry(void **entry, void *user)
3778 {
3779         struct isl_union_pw_multi_aff_bin_data *data = user;
3780         isl_pw_multi_aff *pma2 = *entry;
3781
3782         if (!isl_space_tuple_match(data->pma->dim, isl_dim_in,
3783                                  pma2->dim, isl_dim_in))
3784                 return 0;
3785
3786         pma2 = isl_pw_multi_aff_flat_range_product(
3787                                         isl_pw_multi_aff_copy(data->pma),
3788                                         isl_pw_multi_aff_copy(pma2));
3789
3790         data->res = isl_union_pw_multi_aff_add_pw_multi_aff(data->res, pma2);
3791
3792         return 0;
3793 }
3794
3795 /* Given two isl_union_pw_multi_affs A -> B and C -> D,
3796  * construct an isl_union_pw_multi_aff (A * C) -> (B, D).
3797  */
3798 __isl_give isl_union_pw_multi_aff *isl_union_pw_multi_aff_flat_range_product(
3799         __isl_take isl_union_pw_multi_aff *upma1,
3800         __isl_take isl_union_pw_multi_aff *upma2)
3801 {
3802         return bin_op(upma1, upma2, &flat_range_product_entry);
3803 }
3804
3805 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3806  * The parameters are assumed to have been aligned.
3807  *
3808  * The implementation essentially performs an isl_pw_*_on_shared_domain,
3809  * except that it works on two different isl_pw_* types.
3810  */
3811 static __isl_give isl_pw_multi_aff *pw_multi_aff_set_pw_aff(
3812         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3813         __isl_take isl_pw_aff *pa)
3814 {
3815         int i, j, n;
3816         isl_pw_multi_aff *res = NULL;
3817
3818         if (!pma || !pa)
3819                 goto error;
3820
3821         if (!isl_space_tuple_match(pma->dim, isl_dim_in, pa->dim, isl_dim_in))
3822                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3823                         "domains don't match", goto error);
3824         if (pos >= isl_pw_multi_aff_dim(pma, isl_dim_out))
3825                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3826                         "index out of bounds", goto error);
3827
3828         n = pma->n * pa->n;
3829         res = isl_pw_multi_aff_alloc_size(isl_pw_multi_aff_get_space(pma), n);
3830
3831         for (i = 0; i < pma->n; ++i) {
3832                 for (j = 0; j < pa->n; ++j) {
3833                         isl_set *common;
3834                         isl_multi_aff *res_ij;
3835                         int empty;
3836
3837                         common = isl_set_intersect(isl_set_copy(pma->p[i].set),
3838                                                    isl_set_copy(pa->p[j].set));
3839                         empty = isl_set_plain_is_empty(common);
3840                         if (empty < 0 || empty) {
3841                                 isl_set_free(common);
3842                                 if (empty < 0)
3843                                         goto error;
3844                                 continue;
3845                         }
3846
3847                         res_ij = isl_multi_aff_set_aff(
3848                                         isl_multi_aff_copy(pma->p[i].maff), pos,
3849                                         isl_aff_copy(pa->p[j].aff));
3850                         res_ij = isl_multi_aff_gist(res_ij,
3851                                         isl_set_copy(common));
3852
3853                         res = isl_pw_multi_aff_add_piece(res, common, res_ij);
3854                 }
3855         }
3856
3857         isl_pw_multi_aff_free(pma);
3858         isl_pw_aff_free(pa);
3859         return res;
3860 error:
3861         isl_pw_multi_aff_free(pma);
3862         isl_pw_aff_free(pa);
3863         return isl_pw_multi_aff_free(res);
3864 }
3865
3866 /* Replace the affine expressions at position "pos" in "pma" by "pa".
3867  */
3868 __isl_give isl_pw_multi_aff *isl_pw_multi_aff_set_pw_aff(
3869         __isl_take isl_pw_multi_aff *pma, unsigned pos,
3870         __isl_take isl_pw_aff *pa)
3871 {
3872         if (!pma || !pa)
3873                 goto error;
3874         if (isl_space_match(pma->dim, isl_dim_param, pa->dim, isl_dim_param))
3875                 return pw_multi_aff_set_pw_aff(pma, pos, pa);
3876         if (!isl_space_has_named_params(pma->dim) ||
3877             !isl_space_has_named_params(pa->dim))
3878                 isl_die(isl_pw_multi_aff_get_ctx(pma), isl_error_invalid,
3879                         "unaligned unnamed parameters", goto error);
3880         pma = isl_pw_multi_aff_align_params(pma, isl_pw_aff_get_space(pa));
3881         pa = isl_pw_aff_align_params(pa, isl_pw_multi_aff_get_space(pma));
3882         return pw_multi_aff_set_pw_aff(pma, pos, pa);
3883 error:
3884         isl_pw_multi_aff_free(pma);
3885         isl_pw_aff_free(pa);
3886         return NULL;
3887 }
3888
3889 #undef BASE
3890 #define BASE pw_aff
3891
3892 #include <isl_multi_templ.c>