Coverage Report

Created: 2017-03-27 23:01

/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/llvm/tools/polly/lib/External/isl/isl_transitive_closure.c
Line
Count
Source (jump to first uncovered line)
1
/*
2
 * Copyright 2010      INRIA Saclay
3
 *
4
 * Use of this software is governed by the MIT license
5
 *
6
 * Written by Sven Verdoolaege, INRIA Saclay - Ile-de-France,
7
 * Parc Club Orsay Universite, ZAC des vignes, 4 rue Jacques Monod,
8
 * 91893 Orsay, France 
9
 */
10
11
#include <isl_ctx_private.h>
12
#include <isl_map_private.h>
13
#include <isl/map.h>
14
#include <isl_seq.h>
15
#include <isl_space_private.h>
16
#include <isl_lp_private.h>
17
#include <isl/union_map.h>
18
#include <isl_mat_private.h>
19
#include <isl_vec_private.h>
20
#include <isl_options_private.h>
21
#include <isl_tarjan.h>
22
23
int isl_map_is_transitively_closed(__isl_keep isl_map *map)
24
96
{
25
96
  isl_map *map2;
26
96
  int closed;
27
96
28
96
  map2 = isl_map_apply_range(isl_map_copy(map), isl_map_copy(map));
29
96
  closed = isl_map_is_subset(map2, map);
30
96
  isl_map_free(map2);
31
96
32
96
  return closed;
33
96
}
34
35
int isl_union_map_is_transitively_closed(__isl_keep isl_union_map *umap)
36
53
{
37
53
  isl_union_map *umap2;
38
53
  int closed;
39
53
40
53
  umap2 = isl_union_map_apply_range(isl_union_map_copy(umap),
41
53
            isl_union_map_copy(umap));
42
53
  closed = isl_union_map_is_subset(umap2, umap);
43
53
  isl_union_map_free(umap2);
44
53
45
53
  return closed;
46
53
}
47
 
48
/* Given a map that represents a path with the length of the path
49
 * encoded as the difference between the last output coordindate
50
 * and the last input coordinate, set this length to either
51
 * exactly "length" (if "exactly" is set) or at least "length"
52
 * (if "exactly" is not set).
53
 */
54
static __isl_give isl_map *set_path_length(__isl_take isl_map *map,
55
  int exactly, int length)
56
291
{
57
291
  isl_space *dim;
58
291
  struct isl_basic_map *bmap;
59
291
  unsigned d;
60
291
  unsigned nparam;
61
291
  int k;
62
291
  isl_int *c;
63
291
64
291
  if (!map)
65
0
    return NULL;
66
291
67
291
  dim = isl_map_get_space(map);
68
291
  d = isl_space_dim(dim, isl_dim_in);
69
291
  nparam = isl_space_dim(dim, isl_dim_param);
70
291
  bmap = isl_basic_map_alloc_space(dim, 0, 1, 1);
71
291
  if (
exactly291
)
{6
72
6
    k = isl_basic_map_alloc_equality(bmap);
73
6
    if (k < 0)
74
0
      goto error;
75
6
    c = bmap->eq[k];
76
285
  } else {
77
285
    k = isl_basic_map_alloc_inequality(bmap);
78
285
    if (k < 0)
79
0
      goto error;
80
285
    c = bmap->ineq[k];
81
285
  }
82
291
  isl_seq_clr(c, 1 + isl_basic_map_total_dim(bmap));
83
291
  isl_int_set_si(c[0], -length);
84
291
  isl_int_set_si(c[1 + nparam + d - 1], -1);
85
291
  isl_int_set_si(c[1 + nparam + d + d - 1], 1);
86
291
87
291
  bmap = isl_basic_map_finalize(bmap);
88
291
  map = isl_map_intersect(map, isl_map_from_basic_map(bmap));
89
291
90
291
  return map;
91
0
error:
92
0
  isl_basic_map_free(bmap);
93
0
  isl_map_free(map);
94
0
  return NULL;
95
291
}
96
97
/* Check whether the overapproximation of the power of "map" is exactly
98
 * the power of "map".  Let R be "map" and A_k the overapproximation.
99
 * The approximation is exact if
100
 *
101
 *  A_1 = R
102
 *  A_k = A_{k-1} \circ R     k >= 2
103
 *
104
 * Since A_k is known to be an overapproximation, we only need to check
105
 *
106
 *  A_1 \subset R
107
 *  A_k \subset A_{k-1} \circ R   k >= 2
108
 *
109
 * In practice, "app" has an extra input and output coordinate
110
 * to encode the length of the path.  So, we first need to add
111
 * this coordinate to "map" and set the length of the path to
112
 * one.
113
 */
114
static int check_power_exactness(__isl_take isl_map *map,
115
  __isl_take isl_map *app)
116
3
{
117
3
  int exact;
118
3
  isl_map *app_1;
119
3
  isl_map *app_2;
120
3
121
3
  map = isl_map_add_dims(map, isl_dim_in, 1);
122
3
  map = isl_map_add_dims(map, isl_dim_out, 1);
123
3
  map = set_path_length(map, 1, 1);
124
3
125
3
  app_1 = set_path_length(isl_map_copy(app), 1, 1);
126
3
127
3
  exact = isl_map_is_subset(app_1, map);
128
3
  isl_map_free(app_1);
129
3
130
3
  if (
!exact || 3
exact < 03
)
{0
131
0
    isl_map_free(app);
132
0
    isl_map_free(map);
133
0
    return exact;
134
0
  }
135
3
136
3
  app_1 = set_path_length(isl_map_copy(app), 0, 1);
137
3
  app_2 = set_path_length(app, 0, 2);
138
3
  app_1 = isl_map_apply_range(map, app_1);
139
3
140
3
  exact = isl_map_is_subset(app_2, app_1);
141
3
142
3
  isl_map_free(app_1);
143
3
  isl_map_free(app_2);
144
3
145
3
  return exact;
146
3
}
147
148
/* Check whether the overapproximation of the power of "map" is exactly
149
 * the power of "map", possibly after projecting out the power (if "project"
150
 * is set).
151
 *
152
 * If "project" is set and if "steps" can only result in acyclic paths,
153
 * then we check
154
 *
155
 *  A = R \cup (A \circ R)
156
 *
157
 * where A is the overapproximation with the power projected out, i.e.,
158
 * an overapproximation of the transitive closure.
159
 * More specifically, since A is known to be an overapproximation, we check
160
 *
161
 *  A \subset R \cup (A \circ R)
162
 *
163
 * Otherwise, we check if the power is exact.
164
 *
165
 * Note that "app" has an extra input and output coordinate to encode
166
 * the length of the part.  If we are only interested in the transitive
167
 * closure, then we can simply project out these coordinates first.
168
 */
169
static int check_exactness(__isl_take isl_map *map, __isl_take isl_map *app,
170
  int project)
171
105
{
172
105
  isl_map *test;
173
105
  int exact;
174
105
  unsigned d;
175
105
176
105
  if (!project)
177
3
    return check_power_exactness(map, app);
178
105
179
102
  d = isl_map_dim(map, isl_dim_in);
180
102
  app = set_path_length(app, 0, 1);
181
102
  app = isl_map_project_out(app, isl_dim_in, d, 1);
182
102
  app = isl_map_project_out(app, isl_dim_out, d, 1);
183
102
184
102
  app = isl_map_reset_space(app, isl_map_get_space(map));
185
102
186
102
  test = isl_map_apply_range(isl_map_copy(map), isl_map_copy(app));
187
102
  test = isl_map_union(test, isl_map_copy(map));
188
102
189
102
  exact = isl_map_is_subset(app, test);
190
102
191
102
  isl_map_free(app);
192
102
  isl_map_free(test);
193
102
194
102
  isl_map_free(map);
195
102
196
102
  return exact;
197
105
}
198
199
/*
200
 * The transitive closure implementation is based on the paper
201
 * "Computing the Transitive Closure of a Union of Affine Integer
202
 * Tuple Relations" by Anna Beletska, Denis Barthou, Wlodzimierz Bielecki and
203
 * Albert Cohen.
204
 */
205
206
/* Given a set of n offsets v_i (the rows of "steps"), construct a relation
207
 * of the given dimension specification (Z^{n+1} -> Z^{n+1})
208
 * that maps an element x to any element that can be reached
209
 * by taking a non-negative number of steps along any of
210
 * the extended offsets v'_i = [v_i 1].
211
 * That is, construct
212
 *
213
 * { [x] -> [y] : exists k_i >= 0, y = x + \sum_i k_i v'_i }
214
 *
215
 * For any element in this relation, the number of steps taken
216
 * is equal to the difference in the final coordinates.
217
 */
218
static __isl_give isl_map *path_along_steps(__isl_take isl_space *dim,
219
  __isl_keep isl_mat *steps)
220
170
{
221
170
  int i, j, k;
222
170
  struct isl_basic_map *path = NULL;
223
170
  unsigned d;
224
170
  unsigned n;
225
170
  unsigned nparam;
226
170
227
170
  if (
!dim || 170
!steps170
)
228
0
    goto error;
229
170
230
170
  d = isl_space_dim(dim, isl_dim_in);
231
170
  n = steps->n_row;
232
170
  nparam = isl_space_dim(dim, isl_dim_param);
233
170
234
170
  path = isl_basic_map_alloc_space(isl_space_copy(dim), n, d, n);
235
170
236
362
  for (i = 0; 
i < n362
;
++i192
)
{192
237
192
    k = isl_basic_map_alloc_div(path);
238
192
    if (k < 0)
239
0
      goto error;
240
192
    
isl_assert192
(steps->ctx, i == k, goto error);192
241
192
    
isl_int_set_si192
(path->div[k][0], 0);192
242
192
  }
243
170
244
846
  
for (i = 0; 170
i < d846
;
++i676
)
{676
245
676
    k = isl_basic_map_alloc_equality(path);
246
676
    if (k < 0)
247
0
      goto error;
248
676
    isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
249
676
    isl_int_set_si(path->eq[k][1 + nparam + i], 1);
250
676
    isl_int_set_si(path->eq[k][1 + nparam + d + i], -1);
251
676
    if (i == d - 1)
252
362
      
for (j = 0; 170
j < n362
;
++j192
)
253
192
        isl_int_set_si(path->eq[k][1 + nparam + 2 * d + j], 1);
254
676
    else
255
1.07k
      
for (j = 0; 506
j < n1.07k
;
++j573
)
256
573
        isl_int_set(path->eq[k][1 + nparam + 2 * d + j],
257
676
              steps->row[j][i]);
258
676
  }
259
170
260
362
  
for (i = 0; 170
i < n362
;
++i192
)
{192
261
192
    k = isl_basic_map_alloc_inequality(path);
262
192
    if (k < 0)
263
0
      goto error;
264
192
    isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
265
192
    isl_int_set_si(path->ineq[k][1 + nparam + 2 * d + i], 1);
266
192
  }
267
170
268
170
  isl_space_free(dim);
269
170
270
170
  path = isl_basic_map_simplify(path);
271
170
  path = isl_basic_map_finalize(path);
272
170
  return isl_map_from_basic_map(path);
273
0
error:
274
0
  isl_space_free(dim);
275
0
  isl_basic_map_free(path);
276
0
  return NULL;
277
170
}
278
279
120
#define IMPURE    0
280
98
#define PURE_PARAM  1
281
350
#define PURE_VAR  2
282
44
#define MIXED   3
283
284
/* Check whether the parametric constant term of constraint c is never
285
 * positive in "bset".
286
 */
287
static isl_bool parametric_constant_never_positive(
288
  __isl_keep isl_basic_set *bset, isl_int *c, int *div_purity)
289
14
{
290
14
  unsigned d;
291
14
  unsigned n_div;
292
14
  unsigned nparam;
293
14
  int i;
294
14
  int k;
295
14
  isl_bool empty;
296
14
297
14
  n_div = isl_basic_set_dim(bset, isl_dim_div);
298
14
  d = isl_basic_set_dim(bset, isl_dim_set);
299
14
  nparam = isl_basic_set_dim(bset, isl_dim_param);
300
14
301
14
  bset = isl_basic_set_copy(bset);
302
14
  bset = isl_basic_set_cow(bset);
303
14
  bset = isl_basic_set_extend_constraints(bset, 0, 1);
304
14
  k = isl_basic_set_alloc_inequality(bset);
305
14
  if (k < 0)
306
0
    goto error;
307
14
  isl_seq_clr(bset->ineq[k], 1 + isl_basic_set_total_dim(bset));
308
14
  isl_seq_cpy(bset->ineq[k], c, 1 + nparam);
309
27
  for (i = 0; 
i < n_div27
;
++i13
)
{13
310
13
    if (
div_purity[i] != 13
PURE_PARAM13
)
311
1
      continue;
312
12
    
isl_int_set12
(bset->ineq[k][1 + nparam + d + i],12
313
12
          c[1 + nparam + d + i]);
314
12
  }
315
14
  isl_int_sub_ui(bset->ineq[k][0], bset->ineq[k][0], 1);
316
14
  empty = isl_basic_set_is_empty(bset);
317
14
  isl_basic_set_free(bset);
318
14
319
14
  return empty;
320
0
error:
321
0
  isl_basic_set_free(bset);
322
0
  return isl_bool_error;
323
14
}
324
325
/* Return PURE_PARAM if only the coefficients of the parameters are non-zero.
326
 * Return PURE_VAR if only the coefficients of the set variables are non-zero.
327
 * Return MIXED if only the coefficients of the parameters and the set
328
 *  variables are non-zero and if moreover the parametric constant
329
 *  can never attain positive values.
330
 * Return IMPURE otherwise.
331
 */
332
static int purity(__isl_keep isl_basic_set *bset, isl_int *c, int *div_purity,
333
  int eq)
334
69
{
335
69
  unsigned d;
336
69
  unsigned n_div;
337
69
  unsigned nparam;
338
69
  isl_bool empty;
339
69
  int i;
340
69
  int p = 0, v = 0;
341
69
342
69
  n_div = isl_basic_set_dim(bset, isl_dim_div);
343
69
  d = isl_basic_set_dim(bset, isl_dim_set);
344
69
  nparam = isl_basic_set_dim(bset, isl_dim_param);
345
69
346
157
  for (i = 0; 
i < n_div157
;
++i88
)
{88
347
88
    if (isl_int_is_zero(c[1 + nparam + d + i]))
348
75
      continue;
349
13
    switch (div_purity[i]) {
350
12
    
case 12
PURE_PARAM12
: p = 1; break;
351
1
    
case 1
PURE_VAR1
: v = 1; break;
352
0
    
default: return 0
IMPURE0
;
353
13
    }
354
13
  }
355
69
  
if (69
!p && 69
isl_seq_first_non_zero(c + 1, nparam) == -157
)
356
40
    
return 40
PURE_VAR40
;
357
29
  
if (29
!v && 29
isl_seq_first_non_zero(c + 1 + nparam, d) == -129
)
358
15
    
return 15
PURE_PARAM15
;
359
29
360
14
  empty = parametric_constant_never_positive(bset, c, div_purity);
361
14
  if (
eq && 14
empty >= 00
&&
!empty0
)
{0
362
0
    isl_seq_neg(c, c, 1 + nparam + d + n_div);
363
0
    empty = parametric_constant_never_positive(bset, c, div_purity);
364
0
  }
365
14
366
14
  return empty < 0 ? 
-10
:
empty ? 14
MIXED1
:
IMPURE13
;
367
29
}
368
369
/* Return an array of integers indicating the type of each div in bset.
370
 * If the div is (recursively) defined in terms of only the parameters,
371
 * then the type is PURE_PARAM.
372
 * If the div is (recursively) defined in terms of only the set variables,
373
 * then the type is PURE_VAR.
374
 * Otherwise, the type is IMPURE.
375
 */
376
static __isl_give int *get_div_purity(__isl_keep isl_basic_set *bset)
377
10
{
378
10
  int i, j;
379
10
  int *div_purity;
380
10
  unsigned d;
381
10
  unsigned n_div;
382
10
  unsigned nparam;
383
10
384
10
  if (!bset)
385
0
    return NULL;
386
10
387
10
  n_div = isl_basic_set_dim(bset, isl_dim_div);
388
10
  d = isl_basic_set_dim(bset, isl_dim_set);
389
10
  nparam = isl_basic_set_dim(bset, isl_dim_param);
390
10
391
10
  div_purity = isl_alloc_array(bset->ctx, int, n_div);
392
10
  if (
n_div && 10
!div_purity4
)
393
0
    return NULL;
394
10
395
17
  
for (i = 0; 10
i < bset->n_div17
;
++i7
)
{7
396
7
    int p = 0, v = 0;
397
7
    if (
isl_int_is_zero7
(bset->div[i][0]))
{0
398
0
      div_purity[i] = IMPURE;
399
0
      continue;
400
0
    }
401
7
    
if (7
isl_seq_first_non_zero(bset->div[i] + 2, nparam) != -17
)
402
6
      p = 1;
403
7
    if (isl_seq_first_non_zero(bset->div[i] + 2 + nparam, d) != -1)
404
1
      v = 1;
405
10
    for (j = 0; 
j < i10
;
++j3
)
{3
406
3
      if (isl_int_is_zero(bset->div[i][2 + nparam + d + j]))
407
3
        continue;
408
0
      switch (div_purity[j]) {
409
0
      
case 0
PURE_PARAM0
: p = 1; break;
410
0
      
case 0
PURE_VAR0
: v = 1; break;
411
0
      default: p = v = 1; break;
412
0
      }
413
0
    }
414
7
    
div_purity[i] = v ? 7
p ? 1
IMPURE0
:
PURE_VAR1
:
PURE_PARAM6
;
415
7
  }
416
10
417
10
  return div_purity;
418
10
}
419
420
/* Given a path with the as yet unconstrained length at position "pos",
421
 * check if setting the length to zero results in only the identity
422
 * mapping.
423
 */
424
static int empty_path_is_identity(__isl_keep isl_basic_map *path, unsigned pos)
425
10
{
426
10
  isl_basic_map *test = NULL;
427
10
  isl_basic_map *id = NULL;
428
10
  int k;
429
10
  int is_id;
430
10
431
10
  test = isl_basic_map_copy(path);
432
10
  test = isl_basic_map_extend_constraints(test, 1, 0);
433
10
  k = isl_basic_map_alloc_equality(test);
434
10
  if (k < 0)
435
0
    goto error;
436
10
  isl_seq_clr(test->eq[k], 1 + isl_basic_map_total_dim(test));
437
10
  isl_int_set_si(test->eq[k][pos], 1);
438
10
  test = isl_basic_map_gauss(test, NULL);
439
10
  id = isl_basic_map_identity(isl_basic_map_get_space(path));
440
10
  is_id = isl_basic_map_is_equal(test, id);
441
10
  isl_basic_map_free(test);
442
10
  isl_basic_map_free(id);
443
10
  return is_id;
444
0
error:
445
0
  isl_basic_map_free(test);
446
0
  return -1;
447
10
}
448
449
/* If any of the constraints is found to be impure then this function
450
 * sets *impurity to 1.
451
 *
452
 * If impurity is NULL then we are dealing with a non-parametric set
453
 * and so the constraints are obviously PURE_VAR.
454
 */
455
static __isl_give isl_basic_map *add_delta_constraints(
456
  __isl_take isl_basic_map *path,
457
  __isl_keep isl_basic_set *delta, unsigned off, unsigned nparam,
458
  unsigned d, int *div_purity, int eq, int *impurity)
459
38
{
460
38
  int i, k;
461
19
  int n = eq ? 
delta->n_eq19
:
delta->n_ineq19
;
462
19
  isl_int **delta_c = eq ? 
delta->eq19
:
delta->ineq19
;
463
38
  unsigned n_div;
464
38
465
38
  n_div = isl_basic_set_dim(delta, isl_dim_div);
466
38
467
145
  for (i = 0; 
i < n145
;
++i107
)
{107
468
107
    isl_int *path_c;
469
107
    int p = PURE_VAR;
470
107
    if (impurity)
471
69
      p = purity(delta, delta_c[i], div_purity, eq);
472
107
    if (p < 0)
473
0
      goto error;
474
107
    
if (107
p != 107
PURE_VAR107
&&
p != 29
PURE_PARAM29
&&
!*impurity14
)
475
9
      *impurity = 1;
476
107
    if (
p == 107
IMPURE107
)
477
13
      continue;
478
94
    
if (94
eq && 94
p != 43
MIXED43
)
{43
479
43
      k = isl_basic_map_alloc_equality(path);
480
43
      if (k < 0)
481
0
        goto error;
482
43
      path_c = path->eq[k];
483
51
    } else {
484
51
      k = isl_basic_map_alloc_inequality(path);
485
51
      if (k < 0)
486
0
        goto error;
487
51
      path_c = path->ineq[k];
488
51
    }
489
94
    isl_seq_clr(path_c, 1 + isl_basic_map_total_dim(path));
490
94
    if (
p == 94
PURE_VAR94
)
{78
491
78
      isl_seq_cpy(path_c + off,
492
78
            delta_c[i] + 1 + nparam, d);
493
78
      isl_int_set(path_c[off + d], delta_c[i][0]);
494
16
    } else 
if (16
p == 16
PURE_PARAM16
)
{15
495
15
      isl_seq_cpy(path_c, delta_c[i], 1 + nparam);
496
1
    } else {
497
1
      isl_seq_cpy(path_c + off,
498
1
            delta_c[i] + 1 + nparam, d);
499
1
      isl_seq_cpy(path_c, delta_c[i], 1 + nparam);
500
1
    }
501
94
    isl_seq_cpy(path_c + off - n_div,
502
94
          delta_c[i] + 1 + nparam + d, n_div);
503
94
  }
504
38
505
38
  return path;
506
0
error:
507
0
  isl_basic_map_free(path);
508
0
  return NULL;
509
38
}
510
511
/* Given a set of offsets "delta", construct a relation of the
512
 * given dimension specification (Z^{n+1} -> Z^{n+1}) that
513
 * is an overapproximation of the relations that
514
 * maps an element x to any element that can be reached
515
 * by taking a non-negative number of steps along any of
516
 * the elements in "delta".
517
 * That is, construct an approximation of
518
 *
519
 *  { [x] -> [y] : exists f \in \delta, k \in Z :
520
 *          y = x + k [f, 1] and k >= 0 }
521
 *
522
 * For any element in this relation, the number of steps taken
523
 * is equal to the difference in the final coordinates.
524
 *
525
 * In particular, let delta be defined as
526
 *
527
 *  \delta = [p] -> { [x] : A x + a >= 0 and B p + b >= 0 and
528
 *        C x + C'p + c >= 0 and
529
 *        D x + D'p + d >= 0 }
530
 *
531
 * where the constraints C x + C'p + c >= 0 are such that the parametric
532
 * constant term of each constraint j, "C_j x + C'_j p + c_j",
533
 * can never attain positive values, then the relation is constructed as
534
 *
535
 *  { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
536
 *      A f + k a >= 0 and B p + b >= 0 and
537
 *      C f + C'p + c >= 0 and k >= 1 }
538
 *  union { [x] -> [x] }
539
 *
540
 * If the zero-length paths happen to correspond exactly to the identity
541
 * mapping, then we return
542
 *
543
 *  { [x] -> [y] : exists [f, k] \in Z^{n+1} : y = x + f and
544
 *      A f + k a >= 0 and B p + b >= 0 and
545
 *      C f + C'p + c >= 0 and k >= 0 }
546
 *
547
 * instead.
548
 *
549
 * Existentially quantified variables in \delta are handled by
550
 * classifying them as independent of the parameters, purely
551
 * parameter dependent and others.  Constraints containing
552
 * any of the other existentially quantified variables are removed.
553
 * This is safe, but leads to an additional overapproximation.
554
 *
555
 * If there are any impure constraints, then we also eliminate
556
 * the parameters from \delta, resulting in a set
557
 *
558
 *  \delta' = { [x] : E x + e >= 0 }
559
 *
560
 * and add the constraints
561
 *
562
 *      E f + k e >= 0
563
 *
564
 * to the constructed relation.
565
 */
566
static __isl_give isl_map *path_along_delta(__isl_take isl_space *dim,
567
  __isl_take isl_basic_set *delta)
568
10
{
569
10
  isl_basic_map *path = NULL;
570
10
  unsigned d;
571
10
  unsigned n_div;
572
10
  unsigned nparam;
573
10
  unsigned off;
574
10
  int i, k;
575
10
  int is_id;
576
10
  int *div_purity = NULL;
577
10
  int impurity = 0;
578
10
579
10
  if (!delta)
580
0
    goto error;
581
10
  n_div = isl_basic_set_dim(delta, isl_dim_div);
582
10
  d = isl_basic_set_dim(delta, isl_dim_set);
583
10
  nparam = isl_basic_set_dim(delta, isl_dim_param);
584
10
  path = isl_basic_map_alloc_space(isl_space_copy(dim), n_div + d + 1,
585
10
      d + 1 + delta->n_eq, delta->n_eq + delta->n_ineq + 1);
586
10
  off = 1 + nparam + 2 * (d + 1) + n_div;
587
10
588
61
  for (i = 0; 
i < n_div + d + 161
;
++i51
)
{51
589
51
    k = isl_basic_map_alloc_div(path);
590
51
    if (k < 0)
591
0
      goto error;
592
51
    
isl_int_set_si51
(path->div[k][0], 0);51
593
51
  }
594
10
595
54
  
for (i = 0; 10
i < d + 154
;
++i44
)
{44
596
44
    k = isl_basic_map_alloc_equality(path);
597
44
    if (k < 0)
598
0
      goto error;
599
44
    isl_seq_clr(path->eq[k], 1 + isl_basic_map_total_dim(path));
600
44
    isl_int_set_si(path->eq[k][1 + nparam + i], 1);
601
44
    isl_int_set_si(path->eq[k][1 + nparam + d + 1 + i], -1);
602
44
    isl_int_set_si(path->eq[k][off + i], 1);
603
44
  }
604
10
605
10
  div_purity = get_div_purity(delta);
606
10
  if (
n_div && 10
!div_purity4
)
607
0
    goto error;
608
10
609
10
  path = add_delta_constraints(path, delta, off, nparam, d,
610
10
             div_purity, 1, &impurity);
611
10
  path = add_delta_constraints(path, delta, off, nparam, d,
612
10
             div_purity, 0, &impurity);
613
10
  if (
impurity10
)
{9
614
9
    isl_space *dim = isl_basic_set_get_space(delta);
615
9
    delta = isl_basic_set_project_out(delta,
616
9
              isl_dim_param, 0, nparam);
617
9
    delta = isl_basic_set_add_dims(delta, isl_dim_param, nparam);
618
9
    delta = isl_basic_set_reset_space(delta, dim);
619
9
    if (!delta)
620
0
      goto error;
621
9
    path = isl_basic_map_extend_constraints(path, delta->n_eq,
622
9
              delta->n_ineq + 1);
623
9
    path = add_delta_constraints(path, delta, off, nparam, d,
624
9
               NULL, 1, NULL);
625
9
    path = add_delta_constraints(path, delta, off, nparam, d,
626
9
               NULL, 0, NULL);
627
9
    path = isl_basic_map_gauss(path, NULL);
628
9
  }
629
10
630
10
  is_id = empty_path_is_identity(path, off + d);
631
10
  if (is_id < 0)
632
0
    goto error;
633
10
634
10
  k = isl_basic_map_alloc_inequality(path);
635
10
  if (k < 0)
636
0
    goto error;
637
10
  isl_seq_clr(path->ineq[k], 1 + isl_basic_map_total_dim(path));
638
10
  if (!is_id)
639
8
    isl_int_set_si(path->ineq[k][0], -1);
640
10
  isl_int_set_si(path->ineq[k][off + d], 1);
641
10
      
642
10
  free(div_purity);
643
10
  isl_basic_set_free(delta);
644
10
  path = isl_basic_map_finalize(path);
645
10
  if (
is_id10
)
{2
646
2
    isl_space_free(dim);
647
2
    return isl_map_from_basic_map(path);
648
2
  }
649
8
  return isl_basic_map_union(path, isl_basic_map_identity(dim));
650
0
error:
651
0
  free(div_purity);
652
0
  isl_space_free(dim);
653
0
  isl_basic_set_free(delta);
654
0
  isl_basic_map_free(path);
655
0
  return NULL;
656
10
}
657
658
/* Given a dimension specification Z^{n+1} -> Z^{n+1} and a parameter "param",
659
 * construct a map that equates the parameter to the difference
660
 * in the final coordinates and imposes that this difference is positive.
661
 * That is, construct
662
 *
663
 *  { [x,x_s] -> [y,y_s] : k = y_s - x_s > 0 }
664
 */
665
static __isl_give isl_map *equate_parameter_to_length(__isl_take isl_space *dim,
666
  unsigned param)
667
3
{
668
3
  struct isl_basic_map *bmap;
669
3
  unsigned d;
670
3
  unsigned nparam;
671
3
  int k;
672
3
673
3
  d = isl_space_dim(dim, isl_dim_in);
674
3
  nparam = isl_space_dim(dim, isl_dim_param);
675
3
  bmap = isl_basic_map_alloc_space(dim, 0, 1, 1);
676
3
  k = isl_basic_map_alloc_equality(bmap);
677
3
  if (k < 0)
678
0
    goto error;
679
3
  isl_seq_clr(bmap->eq[k], 1 + isl_basic_map_total_dim(bmap));
680
3
  isl_int_set_si(bmap->eq[k][1 + param], -1);
681
3
  isl_int_set_si(bmap->eq[k][1 + nparam + d - 1], -1);
682
3
  isl_int_set_si(bmap->eq[k][1 + nparam + d + d - 1], 1);
683
3
684
3
  k = isl_basic_map_alloc_inequality(bmap);
685
3
  if (k < 0)
686
0
    goto error;
687
3
  isl_seq_clr(bmap->ineq[k], 1 + isl_basic_map_total_dim(bmap));
688
3
  isl_int_set_si(bmap->ineq[k][1 + param], 1);
689
3
  isl_int_set_si(bmap->ineq[k][0], -1);
690
3
691
3
  bmap = isl_basic_map_finalize(bmap);
692
3
  return isl_map_from_basic_map(bmap);
693
0
error:
694
0
  isl_basic_map_free(bmap);
695
0
  return NULL;
696
3
}
697
698
/* Check whether "path" is acyclic, where the last coordinates of domain
699
 * and range of path encode the number of steps taken.
700
 * That is, check whether
701
 *
702
 *  { d | d = y - x and (x,y) in path }
703
 *
704
 * does not contain any element with positive last coordinate (positive length)
705
 * and zero remaining coordinates (cycle).
706
 */
707
static int is_acyclic(__isl_take isl_map *path)
708
102
{
709
102
  int i;
710
102
  int acyclic;
711
102
  unsigned dim;
712
102
  struct isl_set *delta;
713
102
714
102
  delta = isl_map_deltas(path);
715
102
  dim = isl_set_dim(delta, isl_dim_set);
716
513
  for (i = 0; 
i < dim513
;
++i411
)
{411
717
411
    if (i == dim -1)
718
102
      delta = isl_set_lower_bound_si(delta, isl_dim_set, i, 1);
719
411
    else
720
309
      delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
721
411
  }
722
102
723
102
  acyclic = isl_set_is_empty(delta);
724
102
  isl_set_free(delta);
725
102
726
102
  return acyclic;
727
102
}
728
729
/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
730
 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
731
 * construct a map that is an overapproximation of the map
732
 * that takes an element from the space D \times Z to another
733
 * element from the same space, such that the first n coordinates of the
734
 * difference between them is a sum of differences between images
735
 * and pre-images in one of the R_i and such that the last coordinate
736
 * is equal to the number of steps taken.
737
 * That is, let
738
 *
739
 *  \Delta_i = { y - x | (x, y) in R_i }
740
 *
741
 * then the constructed map is an overapproximation of
742
 *
743
 *  { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
744
 *        d = (\sum_i k_i \delta_i, \sum_i k_i) }
745
 *
746
 * The elements of the singleton \Delta_i's are collected as the
747
 * rows of the steps matrix.  For all these \Delta_i's together,
748
 * a single path is constructed.
749
 * For each of the other \Delta_i's, we compute an overapproximation
750
 * of the paths along elements of \Delta_i.
751
 * Since each of these paths performs an addition, composition is
752
 * symmetric and we can simply compose all resulting paths in any order.
753
 */
754
static __isl_give isl_map *construct_extended_path(__isl_take isl_space *dim,
755
  __isl_keep isl_map *map, int *project)
756
179
{
757
179
  struct isl_mat *steps = NULL;
758
179
  struct isl_map *path = NULL;
759
179
  unsigned d;
760
179
  int i, j, n;
761
179
762
179
  if (!map)
763
0
    goto error;
764
179
765
179
  d = isl_map_dim(map, isl_dim_in);
766
179
767
179
  path = isl_map_identity(isl_space_copy(dim));
768
179
769
179
  steps = isl_mat_alloc(map->ctx, map->n, d);
770
179
  if (!steps)
771
0
    goto error;
772
179
773
179
  n = 0;
774
381
  for (i = 0; 
i < map->n381
;
++i202
)
{202
775
202
    struct isl_basic_set *delta;
776
202
777
202
    delta = isl_basic_map_deltas(isl_basic_map_copy(map->p[i]));
778
202
779
780
    for (j = 0; 
j < d780
;
++j578
)
{588
780
588
      isl_bool fixed;
781
588
782
588
      fixed = isl_basic_set_plain_dim_is_fixed(delta, j,
783
588
                  &steps->row[n][j]);
784
588
      if (
fixed < 0588
)
{0
785
0
        isl_basic_set_free(delta);
786
0
        goto error;
787
0
      }
788
588
      
if (588
!fixed588
)
789
10
        break;
790
588
    }
791
202
792
202
793
202
    
if (202
j < d202
)
{10
794
10
      path = isl_map_apply_range(path,
795
10
        path_along_delta(isl_space_copy(dim), delta));
796
10
      path = isl_map_coalesce(path);
797
192
    } else {
798
192
      isl_basic_set_free(delta);
799
192
      ++n;
800
192
    }
801
202
  }
802
179
803
179
  
if (179
n > 0179
)
{170
804
170
    steps->n_row = n;
805
170
    path = isl_map_apply_range(path,
806
170
        path_along_steps(isl_space_copy(dim), steps));
807
170
  }
808
179
809
179
  if (
project && 179
*project105
)
{102
810
102
    *project = is_acyclic(isl_map_copy(path));
811
102
    if (*project < 0)
812
0
      goto error;
813
102
  }
814
179
815
179
  isl_space_free(dim);
816
179
  isl_mat_free(steps);
817
179
  return path;
818
0
error:
819
0
  isl_space_free(dim);
820
0
  isl_mat_free(steps);
821
0
  isl_map_free(path);
822
0
  return NULL;
823
179
}
824
825
static isl_bool isl_set_overlaps(__isl_keep isl_set *set1,
826
  __isl_keep isl_set *set2)
827
379
{
828
379
  isl_set *i;
829
379
  isl_bool no_overlap;
830
379
831
379
  if (
!set1 || 379
!set2379
)
832
0
    return isl_bool_error;
833
379
834
379
  
if (379
!isl_space_tuple_is_equal(set1->dim, isl_dim_set,379
835
379
          set2->dim, isl_dim_set))
836
0
    return isl_bool_false;
837
379
838
379
  i = isl_set_intersect(isl_set_copy(set1), isl_set_copy(set2));
839
379
  no_overlap = isl_set_is_empty(i);
840
379
  isl_set_free(i);
841
379
842
379
  return isl_bool_not(no_overlap);
843
379
}
844
845
/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
846
 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
847
 * construct a map that is an overapproximation of the map
848
 * that takes an element from the dom R \times Z to an
849
 * element from ran R \times Z, such that the first n coordinates of the
850
 * difference between them is a sum of differences between images
851
 * and pre-images in one of the R_i and such that the last coordinate
852
 * is equal to the number of steps taken.
853
 * That is, let
854
 *
855
 *  \Delta_i = { y - x | (x, y) in R_i }
856
 *
857
 * then the constructed map is an overapproximation of
858
 *
859
 *  { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
860
 *        d = (\sum_i k_i \delta_i, \sum_i k_i) and
861
 *        x in dom R and x + d in ran R and
862
 *        \sum_i k_i >= 1 }
863
 */
864
static __isl_give isl_map *construct_component(__isl_take isl_space *dim,
865
  __isl_keep isl_map *map, int *exact, int project)
866
93
{
867
93
  struct isl_set *domain = NULL;
868
93
  struct isl_set *range = NULL;
869
93
  struct isl_map *app = NULL;
870
93
  struct isl_map *path = NULL;
871
93
  isl_bool overlaps;
872
93
873
93
  domain = isl_map_domain(isl_map_copy(map));
874
93
  domain = isl_set_coalesce(domain);
875
93
  range = isl_map_range(isl_map_copy(map));
876
93
  range = isl_set_coalesce(range);
877
93
  overlaps = isl_set_overlaps(domain, range);
878
93
  if (
overlaps < 0 || 93
!overlaps93
)
{0
879
0
    isl_set_free(domain);
880
0
    isl_set_free(range);
881
0
    isl_space_free(dim);
882
0
883
0
    if (overlaps < 0)
884
0
      map = NULL;
885
0
    map = isl_map_copy(map);
886
0
    map = isl_map_add_dims(map, isl_dim_in, 1);
887
0
    map = isl_map_add_dims(map, isl_dim_out, 1);
888
0
    map = set_path_length(map, 1, 1);
889
0
    return map;
890
0
  }
891
93
  app = isl_map_from_domain_and_range(domain, range);
892
93
  app = isl_map_add_dims(app, isl_dim_in, 1);
893
93
  app = isl_map_add_dims(app, isl_dim_out, 1);
894
93
895
93
  path = construct_extended_path(isl_space_copy(dim), map,
896
19
          exact && 
*exact19
?
&project19
: NULL);
897
93
  app = isl_map_intersect(app, path);
898
93
899
93
  if (
exact && 93
*exact19
&&
900
19
      (*exact = check_exactness(isl_map_copy(map), isl_map_copy(app),
901
19
              project)) < 0)
902
0
    goto error;
903
93
904
93
  isl_space_free(dim);
905
93
  app = set_path_length(app, 0, 1);
906
93
  return app;
907
0
error:
908
0
  isl_space_free(dim);
909
0
  isl_map_free(app);
910
0
  return NULL;
911
93
}
912
913
/* Call construct_component and, if "project" is set, project out
914
 * the final coordinates.
915
 */
916
static __isl_give isl_map *construct_projected_component(
917
  __isl_take isl_space *dim,
918
  __isl_keep isl_map *map, int *exact, int project)
919
93
{
920
93
  isl_map *app;
921
93
  unsigned d;
922
93
923
93
  if (!dim)
924
0
    return NULL;
925
93
  d = isl_space_dim(dim, isl_dim_in);
926
93
927
93
  app = construct_component(dim, map, exact, project);
928
93
  if (
project93
)
{90
929
90
    app = isl_map_project_out(app, isl_dim_in, d - 1, 1);
930
90
    app = isl_map_project_out(app, isl_dim_out, d - 1, 1);
931
90
  }
932
93
  return app;
933
93
}
934
935
/* Compute an extended version, i.e., with path lengths, of
936
 * an overapproximation of the transitive closure of "bmap"
937
 * with path lengths greater than or equal to zero and with
938
 * domain and range equal to "dom".
939
 */
940
static __isl_give isl_map *q_closure(__isl_take isl_space *dim,
941
  __isl_take isl_set *dom, __isl_keep isl_basic_map *bmap, int *exact)
942
86
{
943
86
  int project = 1;
944
86
  isl_map *path;
945
86
  isl_map *map;
946
86
  isl_map *app;
947
86
948
86
  dom = isl_set_add_dims(dom, isl_dim_set, 1);
949
86
  app = isl_map_from_domain_and_range(dom, isl_set_copy(dom));
950
86
  map = isl_map_from_basic_map(isl_basic_map_copy(bmap));
951
86
  path = construct_extended_path(dim, map, &project);
952
86
  app = isl_map_intersect(app, path);
953
86
954
86
  if ((*exact = check_exactness(map, isl_map_copy(app), project)) < 0)
955
0
    goto error;
956
86
957
86
  return app;
958
0
error:
959
0
  isl_map_free(app);
960
0
  return NULL;
961
86
}
962
963
/* Check whether qc has any elements of length at least one
964
 * with domain and/or range outside of dom and ran.
965
 */
966
static int has_spurious_elements(__isl_keep isl_map *qc,
967
  __isl_keep isl_set *dom, __isl_keep isl_set *ran)
968
84
{
969
84
  isl_set *s;
970
84
  int subset;
971
84
  unsigned d;
972
84
973
84
  if (
!qc || 84
!dom84
||
!ran84
)
974
0
    return -1;
975
84
976
84
  d = isl_map_dim(qc, isl_dim_in);
977
84
978
84
  qc = isl_map_copy(qc);
979
84
  qc = set_path_length(qc, 0, 1);
980
84
  qc = isl_map_project_out(qc, isl_dim_in, d - 1, 1);
981
84
  qc = isl_map_project_out(qc, isl_dim_out, d - 1, 1);
982
84
983
84
  s = isl_map_domain(isl_map_copy(qc));
984
84
  subset = isl_set_is_subset(s, dom);
985
84
  isl_set_free(s);
986
84
  if (subset < 0)
987
0
    goto error;
988
84
  
if (84
!subset84
)
{77
989
77
    isl_map_free(qc);
990
77
    return 1;
991
77
  }
992
84
993
7
  s = isl_map_range(qc);
994
7
  subset = isl_set_is_subset(s, ran);
995
7
  isl_set_free(s);
996
7
997
7
  return subset < 0 ? 
-10
:
!subset7
;
998
0
error:
999
0
  isl_map_free(qc);
1000
0
  return -1;
1001
84
}
1002
1003
87
#define LEFT  2
1004
87
#define RIGHT 1
1005
1006
/* For each basic map in "map", except i, check whether it combines
1007
 * with the transitive closure that is reflexive on C combines
1008
 * to the left and to the right.
1009
 *
1010
 * In particular, if
1011
 *
1012
 *  dom map_j \subseteq C
1013
 *
1014
 * then right[j] is set to 1.  Otherwise, if
1015
 *
1016
 *  ran map_i \cap dom map_j = \emptyset
1017
 *
1018
 * then right[j] is set to 0.  Otherwise, composing to the right
1019
 * is impossible.
1020
 *
1021
 * Similar, for composing to the left, we have if
1022
 *
1023
 *  ran map_j \subseteq C
1024
 *
1025
 * then left[j] is set to 1.  Otherwise, if
1026
 *
1027
 *  dom map_i \cap ran map_j = \emptyset
1028
 *
1029
 * then left[j] is set to 0.  Otherwise, composing to the left
1030
 * is impossible.
1031
 *
1032
 * The return value is or'd with LEFT if composing to the left
1033
 * is possible and with RIGHT if composing to the right is possible.
1034
 */
1035
static int composability(__isl_keep isl_set *C, int i,
1036
  isl_set **dom, isl_set **ran, int *left, int *right,
1037
  __isl_keep isl_map *map)
1038
42
{
1039
42
  int j;
1040
42
  int ok;
1041
42
1042
42
  ok = 
LEFT42
|
RIGHT42
;
1043
126
  for (j = 0; 
j < map->n && 126
ok84
;
++j84
)
{84
1044
84
    isl_bool overlaps, subset;
1045
84
    if (j == i)
1046
42
      continue;
1047
84
1048
42
    
if (42
ok & 42
RIGHT42
)
{42
1049
42
      if (!dom[j])
1050
0
        dom[j] = isl_set_from_basic_set(
1051
0
          isl_basic_map_domain(
1052
0
            isl_basic_map_copy(map->p[j])));
1053
42
      if (!dom[j])
1054
0
        return -1;
1055
42
      overlaps = isl_set_overlaps(ran[i], dom[j]);
1056
42
      if (overlaps < 0)
1057
0
        return -1;
1058
42
      
if (42
!overlaps42
)
1059
0
        right[j] = 0;
1060
42
      else {
1061
42
        subset = isl_set_is_subset(dom[j], C);
1062
42
        if (subset < 0)
1063
0
          return -1;
1064
42
        
if (42
subset42
)
1065
41
          right[j] = 1;
1066
42
        else
1067
1
          
ok &= ~1
RIGHT1
;
1068
42
      }
1069
42
    }
1070
42
1071
42
    
if (42
ok & 42
LEFT42
)
{42
1072
42
      if (!ran[j])
1073
0
        ran[j] = isl_set_from_basic_set(
1074
0
          isl_basic_map_range(
1075
0
            isl_basic_map_copy(map->p[j])));
1076
42
      if (!ran[j])
1077
0
        return -1;
1078
42
      overlaps = isl_set_overlaps(dom[i], ran[j]);
1079
42
      if (overlaps < 0)
1080
0
        return -1;
1081
42
      
if (42
!overlaps42
)
1082
0
        left[j] = 0;
1083
42
      else {
1084
42
        subset = isl_set_is_subset(ran[j], C);
1085
42
        if (subset < 0)
1086
0
          return -1;
1087
42
        
if (42
subset42
)
1088
41
          left[j] = 1;
1089
42
        else
1090
1
          
ok &= ~1
LEFT1
;
1091
42
      }
1092
42
    }
1093
42
  }
1094
42
1095
42
  return ok;
1096
42
}
1097
1098
static __isl_give isl_map *anonymize(__isl_take isl_map *map)
1099
71
{
1100
71
  map = isl_map_reset(map, isl_dim_in);
1101
71
  map = isl_map_reset(map, isl_dim_out);
1102
71
  return map;
1103
71
}
1104
1105
/* Return a map that is a union of the basic maps in "map", except i,
1106
 * composed to left and right with qc based on the entries of "left"
1107
 * and "right".
1108
 */
1109
static __isl_give isl_map *compose(__isl_keep isl_map *map, int i,
1110
  __isl_take isl_map *qc, int *left, int *right)
1111
7
{
1112
7
  int j;
1113
7
  isl_map *comp;
1114
7
1115
7
  comp = isl_map_empty(isl_map_get_space(map));
1116
21
  for (j = 0; 
j < map->n21
;
++j14
)
{14
1117
14
    isl_map *map_j;
1118
14
1119
14
    if (j == i)
1120
7
      continue;
1121
14
1122
7
    map_j = isl_map_from_basic_map(isl_basic_map_copy(map->p[j]));
1123
7
    map_j = anonymize(map_j);
1124
7
    if (
left && 7
left[j]7
)
1125
7
      map_j = isl_map_apply_range(map_j, isl_map_copy(qc));
1126
7
    if (
right && 7
right[j]7
)
1127
7
      map_j = isl_map_apply_range(isl_map_copy(qc), map_j);
1128
7
    comp = isl_map_union(comp, map_j);
1129
7
  }
1130
7
1131
7
  comp = isl_map_compute_divs(comp);
1132
7
  comp = isl_map_coalesce(comp);
1133
7
1134
7
  isl_map_free(qc);
1135
7
1136
7
  return comp;
1137
7
}
1138
1139
/* Compute the transitive closure of "map" incrementally by
1140
 * computing
1141
 *
1142
 *  map_i^+ \cup qc^+
1143
 *
1144
 * or
1145
 *
1146
 *  map_i^+ \cup ((id \cup map_i^) \circ qc^+)
1147
 *
1148
 * or
1149
 *
1150
 *  map_i^+ \cup (qc^+ \circ (id \cup map_i^))
1151
 *
1152
 * depending on whether left or right are NULL.
1153
 */
1154
static __isl_give isl_map *compute_incremental(
1155
  __isl_take isl_space *dim, __isl_keep isl_map *map,
1156
  int i, __isl_take isl_map *qc, int *left, int *right, int *exact)
1157
3
{
1158
3
  isl_map *map_i;
1159
3
  isl_map *tc;
1160
3
  isl_map *rtc = NULL;
1161
3
1162
3
  if (!map)
1163
0
    goto error;
1164
3
  
isl_assert3
(map->ctx, left || right, goto error);3
1165
3
1166
3
  map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));
1167
3
  tc = construct_projected_component(isl_space_copy(dim), map_i,
1168
3
            exact, 1);
1169
3
  isl_map_free(map_i);
1170
3
1171
3
  if (*exact)
1172
3
    qc = isl_map_transitive_closure(qc, exact);
1173
3
1174
3
  if (
!*exact3
)
{0
1175
0
    isl_space_free(dim);
1176
0
    isl_map_free(tc);
1177
0
    isl_map_free(qc);
1178
0
    return isl_map_universe(isl_map_get_space(map));
1179
0
  }
1180
3
1181
3
  
if (3
!left || 3
!right3
)
1182
0
    rtc = isl_map_union(isl_map_copy(tc),
1183
0
            isl_map_identity(isl_map_get_space(tc)));
1184
3
  if (!right)
1185
0
    qc = isl_map_apply_range(rtc, qc);
1186
3
  if (!left)
1187
0
    qc = isl_map_apply_range(qc, rtc);
1188
3
  qc = isl_map_union(tc, qc);
1189
3
1190
3
  isl_space_free(dim);
1191
3
1192
3
  return qc;
1193
0
error:
1194
0
  isl_space_free(dim);
1195
0
  isl_map_free(qc);
1196
0
  return NULL;
1197
3
}
1198
1199
/* Given a map "map", try to find a basic map such that
1200
 * map^+ can be computed as
1201
 *
1202
 * map^+ = map_i^+ \cup
1203
 *    \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+
1204
 *
1205
 * with C the simple hull of the domain and range of the input map.
1206
 * map_i^ \cup Id_C is computed by allowing the path lengths to be zero
1207
 * and by intersecting domain and range with C.
1208
 * Of course, we need to check that this is actually equal to map_i^ \cup Id_C.
1209
 * Also, we only use the incremental computation if all the transitive
1210
 * closures are exact and if the number of basic maps in the union,
1211
 * after computing the integer divisions, is smaller than the number
1212
 * of basic maps in the input map.
1213
 */
1214
static int incemental_on_entire_domain(__isl_keep isl_space *dim,
1215
  __isl_keep isl_map *map,
1216
  isl_set **dom, isl_set **ran, int *left, int *right,
1217
  __isl_give isl_map **res)
1218
24
{
1219
24
  int i;
1220
24
  isl_set *C;
1221
24
  unsigned d;
1222
24
1223
24
  *res = NULL;
1224
24
1225
24
  C = isl_set_union(isl_map_domain(isl_map_copy(map)),
1226
24
        isl_map_range(isl_map_copy(map)));
1227
24
  C = isl_set_from_basic_set(isl_set_simple_hull(C));
1228
24
  if (!C)
1229
0
    return -1;
1230
24
  
if (24
C->n != 124
)
{0
1231
0
    isl_set_free(C);
1232
0
    return 0;
1233
0
  }
1234
24
1235
24
  d = isl_map_dim(map, isl_dim_in);
1236
24
1237
66
  for (i = 0; 
i < map->n66
;
++i42
)
{45
1238
45
    isl_map *qc;
1239
45
    int exact_i, spurious;
1240
45
    int j;
1241
45
    dom[i] = isl_set_from_basic_set(isl_basic_map_domain(
1242
45
          isl_basic_map_copy(map->p[i])));
1243
45
    ran[i] = isl_set_from_basic_set(isl_basic_map_range(
1244
45
          isl_basic_map_copy(map->p[i])));
1245
45
    qc = q_closure(isl_space_copy(dim), isl_set_copy(C),
1246
45
        map->p[i], &exact_i);
1247
45
    if (!qc)
1248
0
      goto error;
1249
45
    
if (45
!exact_i45
)
{1
1250
1
      isl_map_free(qc);
1251
1
      continue;
1252
1
    }
1253
44
    spurious = has_spurious_elements(qc, dom[i], ran[i]);
1254
44
    if (
spurious44
)
{39
1255
39
      isl_map_free(qc);
1256
39
      if (spurious < 0)
1257
0
        goto error;
1258
39
      continue;
1259
39
    }
1260
5
    qc = isl_map_project_out(qc, isl_dim_in, d, 1);
1261
5
    qc = isl_map_project_out(qc, isl_dim_out, d, 1);
1262
5
    qc = isl_map_compute_divs(qc);
1263
15
    for (j = 0; 
j < map->n15
;
++j10
)
1264
10
      left[j] = right[j] = 1;
1265
5
    qc = compose(map, i, qc, left, right);
1266
5
    if (!qc)
1267
0
      goto error;
1268
5
    
if (5
qc->n >= map->n5
)
{2
1269
2
      isl_map_free(qc);
1270
2
      continue;
1271
2
    }
1272
3
    *res = compute_incremental(isl_space_copy(dim), map, i, qc,
1273
3
        left, right, &exact_i);
1274
3
    if (!*res)
1275
0
      goto error;
1276
3
    
if (3
exact_i3
)
1277
3
      break;
1278
0
    isl_map_free(*res);
1279
0
    *res = NULL;
1280
0
  }
1281
24
1282
24
  isl_set_free(C);
1283
24
1284
24
  return *res != NULL;
1285
0
error:
1286
0
  isl_set_free(C);
1287
0
  return -1;
1288
24
}
1289
1290
/* Try and compute the transitive closure of "map" as
1291
 *
1292
 * map^+ = map_i^+ \cup
1293
 *    \bigcup_j ((map_i^+ \cup Id_C)^+ \circ map_j \circ (map_i^+ \cup Id_C))^+
1294
 *
1295
 * with C either the simple hull of the domain and range of the entire
1296
 * map or the simple hull of domain and range of map_i.
1297
 */
1298
static __isl_give isl_map *incremental_closure(__isl_take isl_space *dim,
1299
  __isl_keep isl_map *map, int *exact, int project)
1300
93
{
1301
93
  int i;
1302
93
  isl_set **dom = NULL;
1303
93
  isl_set **ran = NULL;
1304
93
  int *left = NULL;
1305
93
  int *right = NULL;
1306
93
  isl_set *C;
1307
93
  unsigned d;
1308
93
  isl_map *res = NULL;
1309
93
1310
93
  if (!project)
1311
3
    return construct_projected_component(dim, map, exact, project);
1312
93
1313
90
  
if (90
!map90
)
1314
0
    goto error;
1315
90
  
if (90
map->n <= 190
)
1316
66
    return construct_projected_component(dim, map, exact, project);
1317
90
1318
24
  d = isl_map_dim(map, isl_dim_in);
1319
24
1320
24
  dom = isl_calloc_array(map->ctx, isl_set *, map->n);
1321
24
  ran = isl_calloc_array(map->ctx, isl_set *, map->n);
1322
24
  left = isl_calloc_array(map->ctx, int, map->n);
1323
24
  right = isl_calloc_array(map->ctx, int, map->n);
1324
24
  if (
!ran || 24
!dom24
||
!left24
||
!right24
)
1325
0
    goto error;
1326
24
1327
24
  
if (24
incemental_on_entire_domain(dim, map, dom, ran, left, right, &res) < 024
)
1328
0
    goto error;
1329
24
1330
66
  
for (i = 0; 24
!res && 66
i < map->n63
;
++i42
)
{42
1331
42
    isl_map *qc;
1332
42
    int exact_i, spurious, comp;
1333
42
    if (!dom[i])
1334
0
      dom[i] = isl_set_from_basic_set(
1335
0
          isl_basic_map_domain(
1336
0
            isl_basic_map_copy(map->p[i])));
1337
42
    if (!dom[i])
1338
0
      goto error;
1339
42
    
if (42
!ran[i]42
)
1340
0
      ran[i] = isl_set_from_basic_set(
1341
0
          isl_basic_map_range(
1342
0
            isl_basic_map_copy(map->p[i])));
1343
42
    if (!ran[i])
1344
0
      goto error;
1345
42
    C = isl_set_union(isl_set_copy(dom[i]),
1346
42
              isl_set_copy(ran[i]));
1347
42
    C = isl_set_from_basic_set(isl_set_simple_hull(C));
1348
42
    if (!C)
1349
0
      goto error;
1350
42
    
if (42
C->n != 142
)
{0
1351
0
      isl_set_free(C);
1352
0
      continue;
1353
0
    }
1354
42
    comp = composability(C, i, dom, ran, left, right, map);
1355
42
    if (
!comp || 42
comp < 041
)
{1
1356
1
      isl_set_free(C);
1357
1
      if (comp < 0)
1358
0
        goto error;
1359
1
      continue;
1360
1
    }
1361
41
    qc = q_closure(isl_space_copy(dim), C, map->p[i], &exact_i);
1362
41
    if (!qc)
1363
0
      goto error;
1364
41
    
if (41
!exact_i41
)
{1
1365
1
      isl_map_free(qc);
1366
1
      continue;
1367
1
    }
1368
40
    spurious = has_spurious_elements(qc, dom[i], ran[i]);
1369
40
    if (
spurious40
)
{38
1370
38
      isl_map_free(qc);
1371
38
      if (spurious < 0)
1372
0
        goto error;
1373
38
      continue;
1374
38
    }
1375
2
    qc = isl_map_project_out(qc, isl_dim_in, d, 1);
1376
2
    qc = isl_map_project_out(qc, isl_dim_out, d, 1);
1377
2
    qc = isl_map_compute_divs(qc);
1378
2
    qc = compose(map, i, qc, (comp & 
LEFT2
) ?
left2
: NULL,
1379
2
        (comp & 
RIGHT2
) ?
right2
: NULL);
1380
2
    if (!qc)
1381
0
      goto error;
1382
2
    
if (2
qc->n >= map->n2
)
{2
1383
2
      isl_map_free(qc);
1384
2
      continue;
1385
2
    }
1386
0
    res = compute_incremental(isl_space_copy(dim), map, i, qc,
1387
0
        (comp & 
LEFT0
) ?
left0
: NULL,
1388
0
        (comp & 
RIGHT0
) ?
right0
: NULL, &exact_i);
1389
0
    if (!res)
1390
0
      goto error;
1391
0
    
if (0
exact_i0
)
1392
0
      break;
1393
0
    isl_map_free(res);
1394
0
    res = NULL;
1395
0
  }
1396
24
1397
72
  
for (i = 0; 24
i < map->n72
;
++i48
)
{48
1398
48
    isl_set_free(dom[i]);
1399
48
    isl_set_free(ran[i]);
1400
48
  }
1401
24
  free(dom);
1402
24
  free(ran);
1403
24
  free(left);
1404
24
  free(right);
1405
24
1406
24
  if (
res24
)
{3
1407
3
    isl_space_free(dim);
1408
3
    return res;
1409
3
  }
1410
24
1411
21
  return construct_projected_component(dim, map, exact, project);
1412
0
error:
1413
0
  if (dom)
1414
0
    
for (i = 0; 0
i < map->n0
;
++i0
)
1415
0
      isl_set_free(dom[i]);
1416
0
  free(dom);
1417
0
  if (ran)
1418
0
    
for (i = 0; 0
i < map->n0
;
++i0
)
1419
0
      isl_set_free(ran[i]);
1420
0
  free(ran);
1421
0
  free(left);
1422
0
  free(right);
1423
0
  isl_space_free(dim);
1424
0
  return NULL;
1425
24
}
1426
1427
/* Given an array of sets "set", add "dom" at position "pos"
1428
 * and search for elements at earlier positions that overlap with "dom".
1429
 * If any can be found, then merge all of them, together with "dom", into
1430
 * a single set and assign the union to the first in the array,
1431
 * which becomes the new group leader for all groups involved in the merge.
1432
 * During the search, we only consider group leaders, i.e., those with
1433
 * group[i] = i, as the other sets have already been combined
1434
 * with one of the group leaders.
1435
 */
1436
static int merge(isl_set **set, int *group, __isl_take isl_set *dom, int pos)
1437
298
{
1438
298
  int i;
1439
298
1440
298
  group[pos] = pos;
1441
298
  set[pos] = isl_set_copy(dom);
1442
298
1443
635
  for (i = pos - 1; 
i >= 0635
;
--i337
)
{337
1444
337
    isl_bool o;
1445
337
1446
337
    if (group[i] != i)
1447
135
      continue;
1448
337
1449
202
    o = isl_set_overlaps(set[i], dom);
1450
202
    if (o < 0)
1451
0
      goto error;
1452
202
    
if (202
!o202
)
1453
9
      continue;
1454
202
1455
193
    set[i] = isl_set_union(set[i], set[group[pos]]);
1456
193
    set[group[pos]] = NULL;
1457
193
    if (!set[i])
1458
0
      goto error;
1459
193
    group[group[pos]] = i;
1460
193
    group[pos] = i;
1461
193
  }
1462
298
1463
298
  isl_set_free(dom);
1464
298
  return 0;
1465
0
error:
1466
0
  isl_set_free(dom);
1467
0
  return -1;
1468
298
}
1469
1470
/* Replace each entry in the n by n grid of maps by the cross product
1471
 * with the relation { [i] -> [i + 1] }.
1472
 */
1473
static int add_length(__isl_keep isl_map *map, isl_map ***grid, int n)
1474
1
{
1475
1
  int i, j, k;
1476
1
  isl_space *dim;
1477
1
  isl_basic_map *bstep;
1478
1
  isl_map *step;
1479
1
  unsigned nparam;
1480
1
1481
1
  if (!map)
1482
0
    return -1;
1483
1
1484
1
  dim = isl_map_get_space(map);
1485
1
  nparam = isl_space_dim(dim, isl_dim_param);
1486
1
  dim = isl_space_drop_dims(dim, isl_dim_in, 0, isl_space_dim(dim, isl_dim_in));
1487
1
  dim = isl_space_drop_dims(dim, isl_dim_out, 0, isl_space_dim(dim, isl_dim_out));
1488
1
  dim = isl_space_add_dims(dim, isl_dim_in, 1);
1489
1
  dim = isl_space_add_dims(dim, isl_dim_out, 1);
1490
1
  bstep = isl_basic_map_alloc_space(dim, 0, 1, 0);
1491
1
  k = isl_basic_map_alloc_equality(bstep);
1492
1
  if (
k < 01
)
{0
1493
0
    isl_basic_map_free(bstep);
1494
0
    return -1;
1495
0
  }
1496
1
  isl_seq_clr(bstep->eq[k], 1 + isl_basic_map_total_dim(bstep));
1497
1
  isl_int_set_si(bstep->eq[k][0], 1);
1498
1
  isl_int_set_si(bstep->eq[k][1 + nparam], 1);
1499
1
  isl_int_set_si(bstep->eq[k][1 + nparam + 1], -1);
1500
1
  bstep = isl_basic_map_finalize(bstep);
1501
1
  step = isl_map_from_basic_map(bstep);
1502
1
1503
3
  for (i = 0; 
i < n3
;
++i2
)
1504
6
    
for (j = 0; 2
j < n6
;
++j4
)
1505
4
      grid[i][j] = isl_map_product(grid[i][j],
1506
4
                 isl_map_copy(step));
1507
1
1508
1
  isl_map_free(step);
1509
1
1510
1
  return 0;
1511
1
}
1512
1513
/* The core of the Floyd-Warshall algorithm.
1514
 * Updates the given n x x matrix of relations in place.
1515
 *
1516
 * The algorithm iterates over all vertices.  In each step, the whole
1517
 * matrix is updated to include all paths that go to the current vertex,
1518
 * possibly stay there a while (including passing through earlier vertices)
1519
 * and then come back.  At the start of each iteration, the diagonal
1520
 * element corresponding to the current vertex is replaced by its
1521
 * transitive closure to account for all indirect paths that stay
1522
 * in the current vertex.
1523
 */
1524
static void floyd_warshall_iterate(isl_map ***grid, int n, int *exact)
1525
76
{
1526
76
  int r, p, q;
1527
76
1528
155
  for (r = 0; 
r < n155
;
++r79
)
{79
1529
79
    int r_exact;
1530
79
    grid[r][r] = isl_map_transitive_closure(grid[r][r],
1531
4
        (exact && 
*exact4
) ?
&r_exact4
: NULL);
1532
79
    if (
exact && 79
*exact4
&&
!r_exact4
)
1533
0
      *exact = 0;
1534
79
1535
164
    for (p = 0; 
p < n164
;
++p85
)
1536
182
      
for (q = 0; 85
q < n182
;
++q97
)
{97
1537
97
        isl_map *loop;
1538
97
        if (
p == r && 97
q == r85
)
1539
79
          continue;
1540
18
        loop = isl_map_apply_range(
1541
18
            isl_map_copy(grid[p][r]),
1542
18
            isl_map_copy(grid[r][q]));
1543
18
        grid[p][q] = isl_map_union(grid[p][q], loop);
1544
18
        loop = isl_map_apply_range(
1545
18
            isl_map_copy(grid[p][r]),
1546
18
          isl_map_apply_range(
1547
18
            isl_map_copy(grid[r][r]),
1548
18
            isl_map_copy(grid[r][q])));
1549
18
        grid[p][q] = isl_map_union(grid[p][q], loop);
1550
18
        grid[p][q] = isl_map_coalesce(grid[p][q]);
1551
18
      }
1552
79
  }
1553
76
}
1554
1555
/* Given a partition of the domains and ranges of the basic maps in "map",
1556
 * apply the Floyd-Warshall algorithm with the elements in the partition
1557
 * as vertices.
1558
 *
1559
 * In particular, there are "n" elements in the partition and "group" is
1560
 * an array of length 2 * map->n with entries in [0,n-1].
1561
 *
1562
 * We first construct a matrix of relations based on the partition information,
1563
 * apply Floyd-Warshall on this matrix of relations and then take the
1564
 * union of all entries in the matrix as the final result.
1565
 *
1566
 * If we are actually computing the power instead of the transitive closure,
1567
 * i.e., when "project" is not set, then the result should have the
1568
 * path lengths encoded as the difference between an extra pair of
1569
 * coordinates.  We therefore apply the nested transitive closures
1570
 * to relations that include these lengths.  In particular, we replace
1571
 * the input relation by the cross product with the unit length relation
1572
 * { [i] -> [i + 1] }.
1573
 */
1574
static __isl_give isl_map *floyd_warshall_with_groups(__isl_take isl_space *dim,
1575
  __isl_keep isl_map *map, int *exact, int project, int *group, int n)
1576
28
{
1577
28
  int i, j, k;
1578
28
  isl_map ***grid = NULL;
1579
28
  isl_map *app;
1580
28
1581
28
  if (!map)
1582
0
    goto error;
1583
28
1584
28
  
if (28
n == 128
)
{26
1585
26
    free(group);
1586
26
    return incremental_closure(dim, map, exact, project);
1587
26
  }
1588
28
1589
2
  
grid = 2
isl_calloc_array2
(map->ctx, isl_map **, n);
1590
2
  if (!grid)
1591
0
    goto error;
1592
6
  
for (i = 0; 2
i < n6
;
++i4
)
{4
1593
4
    grid[i] = isl_calloc_array(map->ctx, isl_map *, n);
1594
4
    if (!grid[i])
1595
0
      goto error;
1596
12
    
for (j = 0; 4
j < n12
;
++j8
)
1597
8
      grid[i][j] = isl_map_empty(isl_map_get_space(map));
1598
4
  }
1599
2
1600
6
  
for (k = 0; 2
k < map->n6
;
++k4
)
{4
1601
4
    i = group[2 * k];
1602
4
    j = group[2 * k + 1];
1603
4
    grid[i][j] = isl_map_union(grid[i][j],
1604
4
        isl_map_from_basic_map(
1605
4
          isl_basic_map_copy(map->p[k])));
1606
4
  }
1607
2
1608
2
  if (
!project && 2
add_length(map, grid, n) < 01
)
1609
0
    goto error;
1610
2
1611
2
  floyd_warshall_iterate(grid, n, exact);
1612
2
1613
2
  app = isl_map_empty(isl_map_get_space(grid[0][0]));
1614
2
1615
6
  for (i = 0; 
i < n6
;
++i4
)
{4
1616
12
    for (j = 0; 
j < n12
;
++j8
)
1617
8
      app = isl_map_union(app, grid[i][j]);
1618
4
    free(grid[i]);
1619
4
  }
1620
2
  free(grid);
1621
2
1622
2
  free(group);
1623
2
  isl_space_free(dim);
1624
2
1625
2
  return app;
1626
0
error:
1627
0
  if (grid)
1628
0
    
for (i = 0; 0
i < n0
;
++i0
)
{0
1629
0
      if (!grid[i])
1630
0
        continue;
1631
0
      
for (j = 0; 0
j < n0
;
++j0
)
1632
0
        isl_map_free(grid[i][j]);
1633
0
      free(grid[i]);
1634
0
    }
1635
0
  free(grid);
1636
0
  free(group);
1637
0
  isl_space_free(dim);
1638
0
  return NULL;
1639
2
}
1640
1641
/* Partition the domains and ranges of the n basic relations in list
1642
 * into disjoint cells.
1643
 *
1644
 * To find the partition, we simply consider all of the domains
1645
 * and ranges in turn and combine those that overlap.
1646
 * "set" contains the partition elements and "group" indicates
1647
 * to which partition element a given domain or range belongs.
1648
 * The domain of basic map i corresponds to element 2 * i in these arrays,
1649
 * while the domain corresponds to element 2 * i + 1.
1650
 * During the construction group[k] is either equal to k,
1651
 * in which case set[k] contains the union of all the domains and
1652
 * ranges in the corresponding group, or is equal to some l < k,
1653
 * with l another domain or range in the same group.
1654
 */
1655
static int *setup_groups(isl_ctx *ctx, __isl_keep isl_basic_map **list, int n,
1656
  isl_set ***set, int *n_group)
1657
102
{
1658
102
  int i;
1659
102
  int *group = NULL;
1660
102
  int g;
1661
102
1662
102
  *set = isl_calloc_array(ctx, isl_set *, 2 * n);
1663
102
  group = isl_alloc_array(ctx, int, 2 * n);
1664
102
1665
102
  if (
!*set || 102
!group102
)
1666
0
    goto error;
1667
102
1668
251
  
for (i = 0; 102
i < n251
;
++i149
)
{149
1669
149
    isl_set *dom;
1670
149
    dom = isl_set_from_basic_set(isl_basic_map_domain(
1671
149
        isl_basic_map_copy(list[i])));
1672
149
    if (merge(*set, group, dom, 2 * i) < 0)
1673
0
      goto error;
1674
149
    dom = isl_set_from_basic_set(isl_basic_map_range(
1675
149
        isl_basic_map_copy(list[i])));
1676
149
    if (merge(*set, group, dom, 2 * i + 1) < 0)
1677
0
      goto error;
1678
149
  }
1679
102
1680
102
  g = 0;
1681
400
  for (i = 0; 
i < 2 * n400
;
++i298
)
1682
298
    
if (298
group[i] == i298
)
{105
1683
105
      if (
g != i105
)
{0
1684
0
        (*set)[g] = (*set)[i];
1685
0
        (*set)[i] = NULL;
1686
0
      }
1687
105
      group[i] = g++;
1688
105
    } else
1689
193
      group[i] = group[group[i]];
1690
102
1691
102
  *n_group = g;
1692
102
1693
102
  return group;
1694
0
error:
1695
0
  if (
*set0
)
{0
1696
0
    for (i = 0; 
i < 2 * n0
;
++i0
)
1697
0
      isl_set_free((*set)[i]);
1698
0
    free(*set);
1699
0
    *set = NULL;
1700
0
  }
1701
0
  free(group);
1702
0
  return NULL;
1703
102
}
1704
1705
/* Check if the domains and ranges of the basic maps in "map" can
1706
 * be partitioned, and if so, apply Floyd-Warshall on the elements
1707
 * of the partition.  Note that we also apply this algorithm
1708
 * if we want to compute the power, i.e., when "project" is not set.
1709
 * However, the results are unlikely to be exact since the recursive
1710
 * calls inside the Floyd-Warshall algorithm typically result in
1711
 * non-linear path lengths quite quickly.
1712
 */
1713
static __isl_give isl_map *floyd_warshall(__isl_take isl_space *dim,
1714
  __isl_keep isl_map *map, int *exact, int project)
1715
95
{
1716
95
  int i;
1717
95
  isl_set **set = NULL;
1718
95
  int *group = NULL;
1719
95
  int n;
1720
95
1721
95
  if (!map)
1722
0
    goto error;
1723
95
  
if (95
map->n <= 195
)
1724
67
    return incremental_closure(dim, map, exact, project);
1725
95
1726
28
  group = setup_groups(map->ctx, map->p, map->n, &set, &n);
1727
28
  if (!group)
1728
0
    goto error;
1729
28
1730
140
  
for (i = 0; 28
i < 2 * map->n140
;
++i112
)
1731
112
    isl_set_free(set[i]);
1732
28
1733
28
  free(set);
1734
28
1735
28
  return floyd_warshall_with_groups(dim, map, exact, project, group, n);
1736
0
error:
1737
0
  isl_space_free(dim);
1738
0
  return NULL;
1739
28
}
1740
1741
/* Structure for representing the nodes of the graph of which
1742
 * strongly connected components are being computed.
1743
 *
1744
 * list contains the actual nodes
1745
 * check_closed is set if we may have used the fact that
1746
 * a pair of basic maps can be interchanged
1747
 */
1748
struct isl_tc_follows_data {
1749
  isl_basic_map **list;
1750
  int check_closed;
1751
};
1752
1753
/* Check whether in the computation of the transitive closure
1754
 * "list[i]" (R_1) should follow (or be part of the same component as)
1755
 * "list[j]" (R_2).
1756
 *
1757
 * That is check whether
1758
 *
1759
 *  R_1 \circ R_2
1760
 *
1761
 * is a subset of
1762
 *
1763
 *  R_2 \circ R_1
1764
 *
1765
 * If so, then there is no reason for R_1 to immediately follow R_2
1766
 * in any path.
1767
 *
1768
 * *check_closed is set if the subset relation holds while
1769
 * R_1 \circ R_2 is not empty.
1770
 */
1771
static isl_bool basic_map_follows(int i, int j, void *user)
1772
480
{
1773
480
  struct isl_tc_follows_data *data = user;
1774
480
  struct isl_map *map12 = NULL;
1775
480
  struct isl_map *map21 = NULL;
1776
480
  isl_bool subset;
1777
480
1778
480
  if (!isl_space_tuple_is_equal(data->list[i]->dim, isl_dim_in,
1779
480
            data->list[j]->dim, isl_dim_out))
1780
377
    return isl_bool_false;
1781
480
1782
103
  map21 = isl_map_from_basic_map(
1783
103
      isl_basic_map_apply_range(
1784
103
        isl_basic_map_copy(data->list[j]),
1785
103
        isl_basic_map_copy(data->list[i])));
1786
103
  subset = isl_map_is_empty(map21);
1787
103
  if (subset < 0)
1788
0
    goto error;
1789
103
  
if (103
subset103
)
{2
1790
2
    isl_map_free(map21);
1791
2
    return isl_bool_false;
1792
2
  }
1793
103
1794
101
  
if (101
!isl_space_tuple_is_equal(data->list[i]->dim, isl_dim_in,101
1795
101
            data->list[i]->dim, isl_dim_out) ||
1796
101
      !isl_space_tuple_is_equal(data->list[j]->dim, isl_dim_in,
1797
0
            data->list[j]->dim, isl_dim_out)) {
1798
0
    isl_map_free(map21);
1799
0
    return isl_bool_true;
1800
0
  }
1801
101
1802
101
  map12 = isl_map_from_basic_map(
1803
101
      isl_basic_map_apply_range(
1804
101
        isl_basic_map_copy(data->list[i]),
1805
101
        isl_basic_map_copy(data->list[j])));
1806
101
1807
101
  subset = isl_map_is_subset(map21, map12);
1808
101
1809
101
  isl_map_free(map12);
1810
101
  isl_map_free(map21);
1811
101
1812
101
  if (subset)
1813
6
    data->check_closed = 1;
1814
101
1815
101
  return subset < 0 ? 
isl_bool_error0
:
!subset101
;
1816
0
error:
1817
0
  isl_map_free(map21);
1818
0
  return isl_bool_error;
1819
101
}
1820
1821
/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D
1822
 * and a dimension specification (Z^{n+1} -> Z^{n+1}),
1823
 * construct a map that is an overapproximation of the map
1824
 * that takes an element from the dom R \times Z to an
1825
 * element from ran R \times Z, such that the first n coordinates of the
1826
 * difference between them is a sum of differences between images
1827
 * and pre-images in one of the R_i and such that the last coordinate
1828
 * is equal to the number of steps taken.
1829
 * If "project" is set, then these final coordinates are not included,
1830
 * i.e., a relation of type Z^n -> Z^n is returned.
1831
 * That is, let
1832
 *
1833
 *  \Delta_i = { y - x | (x, y) in R_i }
1834
 *
1835
 * then the constructed map is an overapproximation of
1836
 *
1837
 *  { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1838
 *        d = (\sum_i k_i \delta_i, \sum_i k_i) and
1839
 *        x in dom R and x + d in ran R }
1840
 *
1841
 * or
1842
 *
1843
 *  { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1844
 *        d = (\sum_i k_i \delta_i) and
1845
 *        x in dom R and x + d in ran R }
1846
 *
1847
 * if "project" is set.
1848
 *
1849
 * We first split the map into strongly connected components, perform
1850
 * the above on each component and then join the results in the correct
1851
 * order, at each join also taking in the union of both arguments
1852
 * to allow for paths that do not go through one of the two arguments.
1853
 */
1854
static __isl_give isl_map *construct_power_components(__isl_take isl_space *dim,
1855
  __isl_keep isl_map *map, int *exact, int project)
1856
91
{
1857
91
  int i, n, c;
1858
91
  struct isl_map *path = NULL;
1859
91
  struct isl_tc_follows_data data;
1860
91
  struct isl_tarjan_graph *g = NULL;
1861
91
  int *orig_exact;
1862
91
  int local_exact;
1863
91
1864
91
  if (!map)
1865
0
    goto error;
1866
91
  
if (91
map->n <= 191
)
1867
61
    return floyd_warshall(dim, map, exact, project);
1868
91
1869
30
  data.list = map->p;
1870
30
  data.check_closed = 0;
1871
30
  g = isl_tarjan_graph_init(map->ctx, map->n, &basic_map_follows, &data);
1872
30
  if (!g)
1873
0
    goto error;
1874
30
1875
30
  orig_exact = exact;
1876
30
  if (
data.check_closed && 30
!exact4
)
1877
0
    exact = &local_exact;
1878
30
1879
30
  c = 0;
1880
30
  i = 0;
1881
30
  n = map->n;
1882
30
  if (project)
1883
27
    path = isl_map_empty(isl_map_get_space(map));
1884
30
  else
1885
3
    path = isl_map_empty(isl_space_copy(dim));
1886
30
  path = anonymize(path);
1887
64
  while (
n64
)
{34
1888
34
    struct isl_map *comp;
1889
34
    isl_map *path_comp, *path_comb;
1890
34
    comp = isl_map_alloc_space(isl_map_get_space(map), n, 0);
1891
96
    while (
g->order[i] != -196
)
{62
1892
62
      comp = isl_map_add_basic_map(comp,
1893
62
            isl_basic_map_copy(map->p[g->order[i]]));
1894
62
      --n;
1895
62
      ++i;
1896
62
    }
1897
34
    path_comp = floyd_warshall(isl_space_copy(dim),
1898
34
            comp, exact, project);
1899
34
    path_comp = anonymize(path_comp);
1900
34
    path_comb = isl_map_apply_range(isl_map_copy(path),
1901
34
            isl_map_copy(path_comp));
1902
34
    path = isl_map_union(path, path_comp);
1903
34
    path = isl_map_union(path, path_comb);
1904
34
    isl_map_free(comp);
1905
34
    ++i;
1906
34
    ++c;
1907
34
  }
1908
30
1909
30
  if (
c > 1 && 30
data.check_closed4
&&
!*exact4
)
{0
1910
0
    int closed;
1911
0
1912
0
    closed = isl_map_is_transitively_closed(path);
1913
0
    if (closed < 0)
1914
0
      goto error;
1915
0
    
if (0
!closed0
)
{0
1916
0
      isl_tarjan_graph_free(g);
1917
0
      isl_map_free(path);
1918
0
      return floyd_warshall(dim, map, orig_exact, project);
1919
0
    }
1920
0
  }
1921
30
1922
30
  isl_tarjan_graph_free(g);
1923
30
  isl_space_free(dim);
1924
30
1925
30
  return path;
1926
0
error:
1927
0
  isl_tarjan_graph_free(g);
1928
0
  isl_space_free(dim);
1929
0
  isl_map_free(path);
1930
0
  return NULL;
1931
30
}
1932
1933
/* Given a union of basic maps R = \cup_i R_i \subseteq D \times D,
1934
 * construct a map that is an overapproximation of the map
1935
 * that takes an element from the space D to another
1936
 * element from the same space, such that the difference between
1937
 * them is a strictly positive sum of differences between images
1938
 * and pre-images in one of the R_i.
1939
 * The number of differences in the sum is equated to parameter "param".
1940
 * That is, let
1941
 *
1942
 *  \Delta_i = { y - x | (x, y) in R_i }
1943
 *
1944
 * then the constructed map is an overapproximation of
1945
 *
1946
 *  { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1947
 *        d = \sum_i k_i \delta_i and k = \sum_i k_i > 0 }
1948
 * or
1949
 *
1950
 *  { (x) -> (x + d) | \exists k_i >= 0, \delta_i \in \Delta_i :
1951
 *        d = \sum_i k_i \delta_i and \sum_i k_i > 0 }
1952
 *
1953
 * if "project" is set.
1954
 *
1955
 * If "project" is not set, then
1956
 * we construct an extended mapping with an extra coordinate
1957
 * that indicates the number of steps taken.  In particular,
1958
 * the difference in the last coordinate is equal to the number
1959
 * of steps taken to move from a domain element to the corresponding
1960
 * image element(s).
1961
 */
1962
static __isl_give isl_map *construct_power(__isl_keep isl_map *map,
1963
  int *exact, int project)
1964
91
{
1965
91
  struct isl_map *app = NULL;
1966
91
  isl_space *dim = NULL;
1967
91
1968
91
  if (!map)
1969
0
    return NULL;
1970
91
1971
91
  dim = isl_map_get_space(map);
1972
91
1973
91
  dim = isl_space_add_dims(dim, isl_dim_in, 1);
1974
91
  dim = isl_space_add_dims(dim, isl_dim_out, 1);
1975
91
1976
91
  app = construct_power_components(isl_space_copy(dim), map,
1977
91
          exact, project);
1978
91
1979
91
  isl_space_free(dim);
1980
91
1981
91
  return app;
1982
91
}
1983
1984
/* Compute the positive powers of "map", or an overapproximation.
1985
 * If the result is exact, then *exact is set to 1.
1986
 *
1987
 * If project is set, then we are actually interested in the transitive
1988
 * closure, so we can use a more relaxed exactness check.
1989
 * The lengths of the paths are also projected out instead of being
1990
 * encoded as the difference between an extra pair of final coordinates.
1991
 */
1992
static __isl_give isl_map *map_power(__isl_take isl_map *map,
1993
  int *exact, int project)
1994
91
{
1995
91
  struct isl_map *app = NULL;
1996
91
1997
91
  if (exact)
1998
17
    *exact = 1;
1999
91
2000
91
  if (!map)
2001
0
    return NULL;
2002
91
2003
91
  
isl_assert91
(map->ctx,91
2004
91
    isl_map_dim(map, isl_dim_in) == isl_map_dim(map, isl_dim_out),
2005
91
    goto error);
2006
91
2007
91
  app = construct_power(map, exact, project);
2008
91
2009
91
  isl_map_free(map);
2010
91
  return app;
2011
0
error:
2012
0
  isl_map_free(map);
2013
0
  isl_map_free(app);
2014
0
  return NULL;
2015
91
}
2016
2017
/* Compute the positive powers of "map", or an overapproximation.
2018
 * The result maps the exponent to a nested copy of the corresponding power.
2019
 * If the result is exact, then *exact is set to 1.
2020
 * map_power constructs an extended relation with the path lengths
2021
 * encoded as the difference between the final coordinates.
2022
 * In the final step, this difference is equated to an extra parameter
2023
 * and made positive.  The extra coordinates are subsequently projected out
2024
 * and the parameter is turned into the domain of the result.
2025
 */
2026
__isl_give isl_map *isl_map_power(__isl_take isl_map *map, int *exact)
2027
3
{
2028
3
  isl_space *target_dim;
2029
3
  isl_space *dim;
2030
3
  isl_map *diff;
2031
3
  unsigned d;
2032
3
  unsigned param;
2033
3
2034
3
  if (!map)
2035
0
    return NULL;
2036
3
2037
3
  d = isl_map_dim(map, isl_dim_in);
2038
3
  param = isl_map_dim(map, isl_dim_param);
2039
3
2040
3
  map = isl_map_compute_divs(map);
2041
3
  map = isl_map_coalesce(map);
2042
3
2043
3
  if (
isl_map_plain_is_empty(map)3
)
{0
2044
0
    map = isl_map_from_range(isl_map_wrap(map));
2045
0
    map = isl_map_add_dims(map, isl_dim_in, 1);
2046
0
    map = isl_map_set_dim_name(map, isl_dim_in, 0, "k");
2047
0
    return map;
2048
0
  }
2049
3
2050
3
  target_dim = isl_map_get_space(map);
2051
3
  target_dim = isl_space_from_range(isl_space_wrap(target_dim));
2052
3
  target_dim = isl_space_add_dims(target_dim, isl_dim_in, 1);
2053
3
  target_dim = isl_space_set_dim_name(target_dim, isl_dim_in, 0, "k");
2054
3
2055
3
  map = map_power(map, exact, 0);
2056
3
2057
3
  map = isl_map_add_dims(map, isl_dim_param, 1);
2058
3
  dim = isl_map_get_space(map);
2059
3
  diff = equate_parameter_to_length(dim, param);
2060
3
  map = isl_map_intersect(map, diff);
2061
3
  map = isl_map_project_out(map, isl_dim_in, d, 1);
2062
3
  map = isl_map_project_out(map, isl_dim_out, d, 1);
2063
3
  map = isl_map_from_range(isl_map_wrap(map));
2064
3
  map = isl_map_move_dims(map, isl_dim_in, 0, isl_dim_param, param, 1);
2065
3
2066
3
  map = isl_map_reset_space(map, target_dim);
2067
3
2068
3
  return map;
2069
3
}
2070
2071
/* Compute a relation that maps each element in the range of the input
2072
 * relation to the lengths of all paths composed of edges in the input
2073
 * relation that end up in the given range element.
2074
 * The result may be an overapproximation, in which case *exact is set to 0.
2075
 * The resulting relation is very similar to the power relation.
2076
 * The difference are that the domain has been projected out, the
2077
 * range has become the domain and the exponent is the range instead
2078
 * of a parameter.
2079
 */
2080
__isl_give isl_map *isl_map_reaching_path_lengths(__isl_take isl_map *map,
2081
  int *exact)
2082
0
{
2083
0
  isl_space *dim;
2084
0
  isl_map *diff;
2085
0
  unsigned d;
2086
0
  unsigned param;
2087
0
2088
0
  if (!map)
2089
0
    return NULL;
2090
0
2091
0
  d = isl_map_dim(map, isl_dim_in);
2092
0
  param = isl_map_dim(map, isl_dim_param);
2093
0
2094
0
  map = isl_map_compute_divs(map);
2095
0
  map = isl_map_coalesce(map);
2096
0
2097
0
  if (
isl_map_plain_is_empty(map)0
)
{0
2098
0
    if (exact)
2099
0
      *exact = 1;
2100
0
    map = isl_map_project_out(map, isl_dim_out, 0, d);
2101
0
    map = isl_map_add_dims(map, isl_dim_out, 1);
2102
0
    return map;
2103
0
  }
2104
0
2105
0
  map = map_power(map, exact, 0);
2106
0
2107
0
  map = isl_map_add_dims(map, isl_dim_param, 1);
2108
0
  dim = isl_map_get_space(map);
2109
0
  diff = equate_parameter_to_length(dim, param);
2110
0
  map = isl_map_intersect(map, diff);
2111
0
  map = isl_map_project_out(map, isl_dim_in, 0, d + 1);
2112
0
  map = isl_map_project_out(map, isl_dim_out, d, 1);
2113
0
  map = isl_map_reverse(map);
2114
0
  map = isl_map_move_dims(map, isl_dim_out, 0, isl_dim_param, param, 1);
2115
0
2116
0
  return map;
2117
0
}
2118
2119
/* Check whether equality i of bset is a pure stride constraint
2120
 * on a single dimension, i.e., of the form
2121
 *
2122
 *  v = k e
2123
 *
2124
 * with k a constant and e an existentially quantified variable.
2125
 */
2126
static isl_bool is_eq_stride(__isl_keep isl_basic_set *bset, int i)
2127
0
{
2128
0
  unsigned nparam;
2129
0
  unsigned d;
2130
0
  unsigned n_div;
2131
0
  int pos1;
2132
0
  int pos2;
2133
0
2134
0
  if (!bset)
2135
0
    return isl_bool_error;
2136
0
2137
0
  
if (0
!0
isl_int_is_zero0
(bset->eq[i][0]))
2138
0
    return isl_bool_false;
2139
0
2140
0
  nparam = isl_basic_set_dim(bset, isl_dim_param);
2141
0
  d = isl_basic_set_dim(bset, isl_dim_set);
2142
0
  n_div = isl_basic_set_dim(bset, isl_dim_div);
2143
0
2144
0
  if (isl_seq_first_non_zero(bset->eq[i] + 1, nparam) != -1)
2145
0
    return isl_bool_false;
2146
0
  pos1 = isl_seq_first_non_zero(bset->eq[i] + 1 + nparam, d);
2147
0
  if (pos1 == -1)
2148
0
    return isl_bool_false;
2149
0
  
if (0
isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + pos1 + 1, 0
2150
0
          d - pos1 - 1) != -1)
2151
0
    return isl_bool_false;
2152
0
2153
0
  pos2 = isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + d, n_div);
2154
0
  if (pos2 == -1)
2155
0
    return isl_bool_false;
2156
0
  
if (0
isl_seq_first_non_zero(bset->eq[i] + 1 + nparam + d + pos2 + 1,0
2157
0
           n_div - pos2 - 1) != -1)
2158
0
    return isl_bool_false;
2159
0
  
if (0
!0
isl_int_is_one0
(bset->eq[i][1 + nparam + pos1]) &&
2160
0
      
!0
isl_int_is_negone0
(bset->eq[i][1 + nparam + pos1]))
2161
0
    return isl_bool_false;
2162
0
2163
0
  return isl_bool_true;
2164
0
}
2165
2166
/* Given a map, compute the smallest superset of this map that is of the form
2167
 *
2168
 *  { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }
2169
 *
2170
 * (where p ranges over the (non-parametric) dimensions),
2171
 * compute the transitive closure of this map, i.e.,
2172
 *
2173
 *  { i -> j : exists k > 0:
2174
 *    k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
2175
 *
2176
 * and intersect domain and range of this transitive closure with
2177
 * the given domain and range.
2178
 *
2179
 * If with_id is set, then try to include as much of the identity mapping
2180
 * as possible, by computing
2181
 *
2182
 *  { i -> j : exists k >= 0:
2183
 *    k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
2184
 *
2185
 * instead (i.e., allow k = 0).
2186
 *
2187
 * In practice, we compute the difference set
2188
 *
2189
 *  delta  = { j - i | i -> j in map },
2190
 *
2191
 * look for stride constraint on the individual dimensions and compute
2192
 * (constant) lower and upper bounds for each individual dimension,
2193
 * adding a constraint for each bound not equal to infinity.
2194
 */
2195
static __isl_give isl_map *box_closure_on_domain(__isl_take isl_map *map,
2196
  __isl_take isl_set *dom, __isl_take isl_set *ran, int with_id)
2197
0
{
2198
0
  int i;
2199
0
  int k;
2200
0
  unsigned d;
2201
0
  unsigned nparam;
2202
0
  unsigned total;
2203
0
  isl_space *dim;
2204
0
  isl_set *delta;
2205
0
  isl_map *app = NULL;
2206
0
  isl_basic_set *aff = NULL;
2207
0
  isl_basic_map *bmap = NULL;
2208
0
  isl_vec *obj = NULL;
2209
0
  isl_int opt;
2210
0
2211
0
  isl_int_init(opt);
2212
0
2213
0
  delta = isl_map_deltas(isl_map_copy(map));
2214
0
2215
0
  aff = isl_set_affine_hull(isl_set_copy(delta));
2216
0
  if (!aff)
2217
0
    goto error;
2218
0
  dim = isl_map_get_space(map);
2219
0
  d = isl_space_dim(dim, isl_dim_in);
2220
0
  nparam = isl_space_dim(dim, isl_dim_param);
2221
0
  total = isl_space_dim(dim, isl_dim_all);
2222
0
  bmap = isl_basic_map_alloc_space(dim,
2223
0
          aff->n_div + 1, aff->n_div, 2 * d + 1);
2224
0
  for (i = 0; 
i < aff->n_div + 10
;
++i0
)
{0
2225
0
    k = isl_basic_map_alloc_div(bmap);
2226
0
    if (k < 0)
2227
0
      goto error;
2228
0
    
isl_int_set_si0
(bmap->div[k][0], 0);0
2229
0
  }
2230
0
  
for (i = 0; 0
i < aff->n_eq0
;
++i0
)
{0
2231
0
    if (!is_eq_stride(aff, i))
2232
0
      continue;
2233
0
    k = isl_basic_map_alloc_equality(bmap);
2234
0
    if (k < 0)
2235
0
      goto error;
2236
0
    isl_seq_clr(bmap->eq[k], 1 + nparam);
2237
0
    isl_seq_cpy(bmap->eq[k] + 1 + nparam + d,
2238
0
        aff->eq[i] + 1 + nparam, d);
2239
0
    isl_seq_neg(bmap->eq[k] + 1 + nparam,
2240
0
        aff->eq[i] + 1 + nparam, d);
2241
0
    isl_seq_cpy(bmap->eq[k] + 1 + nparam + 2 * d,
2242
0
        aff->eq[i] + 1 + nparam + d, aff->n_div);
2243
0
    isl_int_set_si(bmap->eq[k][1 + total + aff->n_div], 0);
2244
0
  }
2245
0
  obj = isl_vec_alloc(map->ctx, 1 + nparam + d);
2246
0
  if (!obj)
2247
0
    goto error;
2248
0
  isl_seq_clr(obj->el, 1 + nparam + d);
2249
0
  for (i = 0; 
i < d0
;
++ i0
)
{0
2250
0
    enum isl_lp_result res;
2251
0
2252
0
    isl_int_set_si(obj->el[1 + nparam + i], 1);
2253
0
2254
0
    res = isl_set_solve_lp(delta, 0, obj->el, map->ctx->one, &opt,
2255
0
          NULL, NULL);
2256
0
    if (res == isl_lp_error)
2257
0
      goto error;
2258
0
    
if (0
res == isl_lp_ok0
)
{0
2259
0
      k = isl_basic_map_alloc_inequality(bmap);
2260
0
      if (k < 0)
2261
0
        goto error;
2262
0
      isl_seq_clr(bmap->ineq[k],
2263
0
          1 + nparam + 2 * d + bmap->n_div);
2264
0
      isl_int_set_si(bmap->ineq[k][1 + nparam + i], -1);
2265
0
      isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], 1);
2266
0
      isl_int_neg(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);
2267
0
    }
2268
0
2269
0
    res = isl_set_solve_lp(delta, 1, obj->el, map->ctx->one, &opt,
2270
0
          NULL, NULL);
2271
0
    if (res == isl_lp_error)
2272
0
      goto error;
2273
0
    
if (0
res == isl_lp_ok0
)
{0
2274
0
      k = isl_basic_map_alloc_inequality(bmap);
2275
0
      if (k < 0)
2276
0
        goto error;
2277
0
      isl_seq_clr(bmap->ineq[k],
2278
0
          1 + nparam + 2 * d + bmap->n_div);
2279
0
      isl_int_set_si(bmap->ineq[k][1 + nparam + i], 1);
2280
0
      isl_int_set_si(bmap->ineq[k][1 + nparam + d + i], -1);
2281
0
      isl_int_set(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], opt);
2282
0
    }
2283
0
2284
0
    
isl_int_set_si0
(obj->el[1 + nparam + i], 0);0
2285
0
  }
2286
0
  k = isl_basic_map_alloc_inequality(bmap);
2287
0
  if (k < 0)
2288
0
    goto error;
2289
0
  isl_seq_clr(bmap->ineq[k],
2290
0
      1 + nparam + 2 * d + bmap->n_div);
2291
0
  if (!with_id)
2292
0
    isl_int_set_si(bmap->ineq[k][0], -1);
2293
0
  isl_int_set_si(bmap->ineq[k][1 + nparam + 2 * d + aff->n_div], 1);
2294
0
2295
0
  app = isl_map_from_domain_and_range(dom, ran);
2296
0
2297
0
  isl_vec_free(obj);
2298
0
  isl_basic_set_free(aff);
2299
0
  isl_map_free(map);
2300
0
  bmap = isl_basic_map_finalize(bmap);
2301
0
  isl_set_free(delta);
2302
0
  isl_int_clear(opt);
2303
0
2304
0
  map = isl_map_from_basic_map(bmap);
2305
0
  map = isl_map_intersect(map, app);
2306
0
2307
0
  return map;
2308
0
error:
2309
0
  isl_vec_free(obj);
2310
0
  isl_basic_map_free(bmap);
2311
0
  isl_basic_set_free(aff);
2312
0
  isl_set_free(dom);
2313
0
  isl_set_free(ran);
2314
0
  isl_map_free(map);
2315
0
  isl_set_free(delta);
2316
0
  isl_int_clear(opt);
2317
0
  return NULL;
2318
0
}
2319
2320
/* Given a map, compute the smallest superset of this map that is of the form
2321
 *
2322
 *  { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }
2323
 *
2324
 * (where p ranges over the (non-parametric) dimensions),
2325
 * compute the transitive closure of this map, i.e.,
2326
 *
2327
 *  { i -> j : exists k > 0:
2328
 *    k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
2329
 *
2330
 * and intersect domain and range of this transitive closure with
2331
 * domain and range of the original map.
2332
 */
2333
static __isl_give isl_map *box_closure(__isl_take isl_map *map)
2334
0
{
2335
0
  isl_set *domain;
2336
0
  isl_set *range;
2337
0
2338
0
  domain = isl_map_domain(isl_map_copy(map));
2339
0
  domain = isl_set_coalesce(domain);
2340
0
  range = isl_map_range(isl_map_copy(map));
2341
0
  range = isl_set_coalesce(range);
2342
0
2343
0
  return box_closure_on_domain(map, domain, range, 0);
2344
0
}
2345
2346
/* Given a map, compute the smallest superset of this map that is of the form
2347
 *
2348
 *  { i -> j : L <= j - i <= U and exists a_p: j_p - i_p = M_p a_p }
2349
 *
2350
 * (where p ranges over the (non-parametric) dimensions),
2351
 * compute the transitive and partially reflexive closure of this map, i.e.,
2352
 *
2353
 *  { i -> j : exists k >= 0:
2354
 *    k L <= j - i <= k U and exists a: j_p - i_p = M_p a_p }
2355
 *
2356
 * and intersect domain and range of this transitive closure with
2357
 * the given domain.
2358
 */
2359
static __isl_give isl_map *box_closure_with_identity(__isl_take isl_map *map,
2360
  __isl_take isl_set *dom)
2361
0
{
2362
0
  return box_closure_on_domain(map, dom, isl_set_copy(dom), 1);
2363
0
}
2364
2365
/* Check whether app is the transitive closure of map.
2366
 * In particular, check that app is acyclic and, if so,
2367
 * check that
2368
 *
2369
 *  app \subset (map \cup (map \circ app))
2370
 */
2371
static int check_exactness_omega(__isl_keep isl_map *map,
2372
  __isl_keep isl_map *app)
2373
0
{
2374
0
  isl_set *delta;
2375
0
  int i;
2376
0
  int is_empty, is_exact;
2377
0
  unsigned d;
2378
0
  isl_map *test;
2379
0
2380
0
  delta = isl_map_deltas(isl_map_copy(app));
2381
0
  d = isl_set_dim(delta, isl_dim_set);
2382
0
  for (i = 0; 
i < d0
;
++i0
)
2383
0
    delta = isl_set_fix_si(delta, isl_dim_set, i, 0);
2384
0
  is_empty = isl_set_is_empty(delta);
2385
0
  isl_set_free(delta);
2386
0
  if (is_empty < 0)
2387
0
    return -1;
2388
0
  
if (0
!is_empty0
)
2389
0
    return 0;
2390
0
2391
0
  test = isl_map_apply_range(isl_map_copy(app), isl_map_copy(map));
2392
0
  test = isl_map_union(test, isl_map_copy(map));
2393
0
  is_exact = isl_map_is_subset(app, test);
2394
0
  isl_map_free(test);
2395
0
2396
0
  return is_exact;
2397
0
}
2398
2399
/* Check if basic map M_i can be combined with all the other
2400
 * basic maps such that
2401
 *
2402
 *  (\cup_j M_j)^+
2403
 *
2404
 * can be computed as
2405
 *
2406
 *  M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+
2407
 *
2408
 * In particular, check if we can compute a compact representation
2409
 * of
2410
 *
2411
 *    M_i^* \circ M_j \circ M_i^*
2412
 *
2413
 * for each j != i.
2414
 * Let M_i^? be an extension of M_i^+ that allows paths
2415
 * of length zero, i.e., the result of box_closure(., 1).
2416
 * The criterion, as proposed by Kelly et al., is that
2417
 * id = M_i^? - M_i^+ can be represented as a basic map
2418
 * and that
2419
 *
2420
 *  id \circ M_j \circ id = M_j
2421
 *
2422
 * for each j != i.
2423
 *
2424
 * If this function returns 1, then tc and qc are set to
2425
 * M_i^+ and M_i^?, respectively.
2426
 */
2427
static int can_be_split_off(__isl_keep isl_map *map, int i,
2428
  __isl_give isl_map **tc, __isl_give isl_map **qc)
2429
0
{
2430
0
  isl_map *map_i, *id = NULL;
2431
0
  int j = -1;
2432
0
  isl_set *C;
2433
0
2434
0
  *tc = NULL;
2435
0
  *qc = NULL;
2436
0
2437
0
  C = isl_set_union(isl_map_domain(isl_map_copy(map)),
2438
0
        isl_map_range(isl_map_copy(map)));
2439
0
  C = isl_set_from_basic_set(isl_set_simple_hull(C));
2440
0
  if (!C)
2441
0
    goto error;
2442
0
2443
0
  map_i = isl_map_from_basic_map(isl_basic_map_copy(map->p[i]));
2444
0
  *tc = box_closure(isl_map_copy(map_i));
2445
0
  *qc = box_closure_with_identity(map_i, C);
2446
0
  id = isl_map_subtract(isl_map_copy(*qc), isl_map_copy(*tc));
2447
0
2448
0
  if (
!id || 0
!*qc0
)
2449
0
    goto error;
2450
0
  
if (0
id->n != 1 || 0
(*qc)->n != 10
)
2451
0
    goto done;
2452
0
2453
0
  
for (j = 0; 0
j < map->n0
;
++j0
)
{0
2454
0
    isl_map *map_j, *test;
2455
0
    int is_ok;
2456
0
2457
0
    if (i == j)
2458
0
      continue;
2459
0
    map_j = isl_map_from_basic_map(
2460
0
          isl_basic_map_copy(map->p[j]));
2461
0
    test = isl_map_apply_range(isl_map_copy(id),
2462
0
            isl_map_copy(map_j));
2463
0
    test = isl_map_apply_range(test, isl_map_copy(id));
2464
0
    is_ok = isl_map_is_equal(test, map_j);
2465
0
    isl_map_free(map_j);
2466
0
    isl_map_free(test);
2467
0
    if (is_ok < 0)
2468
0
      goto error;
2469
0
    
if (0
!is_ok0
)
2470
0
      break;
2471
0
  }
2472
0
2473
0
done:
2474
0
  isl_map_free(id);
2475
0
  if (j == map->n)
2476
0
    return 1;
2477
0
2478
0
  isl_map_free(*qc);
2479
0
  isl_map_free(*tc);
2480
0
  *qc = NULL;
2481
0
  *tc = NULL;
2482
0
2483
0
  return 0;
2484
0
error:
2485
0
  isl_map_free(id);
2486
0
  isl_map_free(*qc);
2487
0
  isl_map_free(*tc);
2488
0
  *qc = NULL;
2489
0
  *tc = NULL;
2490
0
  return -1;
2491
0
}
2492
2493
static __isl_give isl_map *box_closure_with_check(__isl_take isl_map *map,
2494
  int *exact)
2495
0
{
2496
0
  isl_map *app;
2497
0
2498
0
  app = box_closure(isl_map_copy(map));
2499
0
  if (exact)
2500
0
    *exact = check_exactness_omega(map, app);
2501
0
2502
0
  isl_map_free(map);
2503
0
  return app;
2504
0
}
2505
2506
/* Compute an overapproximation of the transitive closure of "map"
2507
 * using a variation of the algorithm from
2508
 * "Transitive Closure of Infinite Graphs and its Applications"
2509
 * by Kelly et al.
2510
 *
2511
 * We first check whether we can can split of any basic map M_i and
2512
 * compute
2513
 *
2514
 *  (\cup_j M_j)^+
2515
 *
2516
 * as
2517
 *
2518
 *  M_i \cup (\cup_{j \ne i} M_i^* \circ M_j \circ M_i^*)^+
2519
 *
2520
 * using a recursive call on the remaining map.
2521
 *
2522
 * If not, we simply call box_closure on the whole map.
2523
 */
2524
static __isl_give isl_map *transitive_closure_omega(__isl_take isl_map *map,
2525
  int *exact)
2526
0
{
2527
0
  int i, j;
2528
0
  int exact_i;
2529
0
  isl_map *app;
2530
0
2531
0
  if (!map)
2532
0
    return NULL;
2533
0
  
if (0
map->n == 10
)
2534
0
    return box_closure_with_check(map, exact);
2535
0
2536
0
  
for (i = 0; 0
i < map->n0
;
++i0
)
{0
2537
0
    int ok;
2538
0
    isl_map *qc, *tc;
2539
0
    ok = can_be_split_off(map, i, &tc, &qc);
2540
0
    if (ok < 0)
2541
0
      goto error;
2542
0
    
if (0
!ok0
)
2543
0
      continue;
2544
0
2545
0
    app = isl_map_alloc_space(isl_map_get_space(map), map->n - 1, 0);
2546
0
2547
0
    for (j = 0; 
j < map->n0
;
++j0
)
{0
2548
0
      if (j == i)
2549
0
        continue;
2550
0
      app = isl_map_add_basic_map(app,
2551
0
            isl_basic_map_copy(map->p[j]));
2552
0
    }
2553
0
2554
0
    app = isl_map_apply_range(isl_map_copy(qc), app);
2555
0
    app = isl_map_apply_range(app, qc);
2556
0
2557
0
    app = isl_map_union(tc, transitive_closure_omega(app, NULL));
2558
0
    exact_i = check_exactness_omega(map, app);
2559
0
    if (
exact_i == 10
)
{0
2560
0
      if (exact)
2561
0
        *exact = exact_i;
2562
0
      isl_map_free(map);
2563
0
      return app;
2564
0
    }
2565
0
    isl_map_free(app);
2566
0
    if (exact_i < 0)
2567
0
      goto error;
2568
0
  }
2569
0
2570
0
  return box_closure_with_check(map, exact);
2571
0
error:
2572
0
  isl_map_free(map);
2573
0
  return NULL;
2574
0
}
2575
2576
/* Compute the transitive closure  of "map", or an overapproximation.
2577
 * If the result is exact, then *exact is set to 1.
2578
 * Simply use map_power to compute the powers of map, but tell
2579
 * it to project out the lengths of the paths instead of equating
2580
 * the length to a parameter.
2581
 */
2582
__isl_give isl_map *isl_map_transitive_closure(__isl_take isl_map *map,
2583
  int *exact)
2584
96
{
2585
96
  isl_space *target_dim;
2586
96
  int closed;
2587
96
2588
96
  if (!map)
2589
0
    goto error;
2590
96
2591
96
  
if (96
map->ctx->opt->closure == 96
ISL_CLOSURE_BOX96
)
2592
0
    return transitive_closure_omega(map, exact);
2593
96
2594
96
  map = isl_map_compute_divs(map);
2595
96
  map = isl_map_coalesce(map);
2596
96
  closed = isl_map_is_transitively_closed(map);
2597
96
  if (closed < 0)
2598
0
    goto error;
2599
96
  
if (96
closed96
)
{8
2600
8
    if (exact)
2601
6
      *exact = 1;
2602
8
    return map;
2603
8
  }
2604
96
2605
88
  target_dim = isl_map_get_space(map);
2606
88
  map = map_power(map, exact, 1);
2607
88
  map = isl_map_reset_space(map, target_dim);
2608
88
2609
88
  return map;
2610
0
error:
2611
0
  isl_map_free(map);
2612
0
  return NULL;
2613
96
}
2614
2615
static isl_stat inc_count(__isl_take isl_map *map, void *user)
2616
146
{
2617
146
  int *n = user;
2618
146
2619
146
  *n += map->n;
2620
146
2621
146
  isl_map_free(map);
2622
146
2623
146
  return isl_stat_ok;
2624
146
}
2625
2626
static isl_stat collect_basic_map(__isl_take isl_map *map, void *user)
2627
106
{
2628
106
  int i;
2629
106
  isl_basic_map ***next = user;
2630
106
2631
252
  for (i = 0; 
i < map->n252
;
++i146
)
{146
2632
146
    **next = isl_basic_map_copy(map->p[i]);
2633
146
    if (!**next)
2634
0
      goto error;
2635
146
    (*next)++;
2636
146
  }
2637
106
2638
106
  isl_map_free(map);
2639
106
  return isl_stat_ok;
2640
0
error:
2641
0
  isl_map_free(map);
2642
0
  return isl_stat_error;
2643
106
}
2644
2645
/* Perform Floyd-Warshall on the given list of basic relations.
2646
 * The basic relations may live in different dimensions,
2647
 * but basic relations that get assigned to the diagonal of the
2648
 * grid have domains and ranges of the same dimension and so
2649
 * the standard algorithm can be used because the nested transitive
2650
 * closures are only applied to diagonal elements and because all
2651
 * compositions are peformed on relations with compatible domains and ranges.
2652
 */
2653
static __isl_give isl_union_map *union_floyd_warshall_on_list(isl_ctx *ctx,
2654
  __isl_keep isl_basic_map **list, int n, int *exact)
2655
74
{
2656
74
  int i, j, k;
2657
74
  int n_group;
2658
74
  int *group = NULL;
2659
74
  isl_set **set = NULL;
2660
74
  isl_map ***grid = NULL;
2661
74
  isl_union_map *app;
2662
74
2663
74
  group = setup_groups(ctx, list, n, &set, &n_group);
2664
74
  if (!group)
2665
0
    goto error;
2666
74
2667
74
  
grid = 74
isl_calloc_array74
(ctx, isl_map **, n_group);
2668
74
  if (!grid)
2669
0
    goto error;
2670
149
  
for (i = 0; 74
i < n_group149
;
++i75
)
{75
2671
75
    grid[i] = isl_calloc_array(ctx, isl_map *, n_group);
2672
75
    if (!grid[i])
2673
0
      goto error;
2674
152
    
for (j = 0; 75
j < n_group152
;
++j77
)
{77
2675
77
      isl_space *dim1, *dim2, *dim;
2676
77
      dim1 = isl_space_reverse(isl_set_get_space(set[i]));
2677
77
      dim2 = isl_set_get_space(set[j]);
2678
77
      dim = isl_space_join(dim1, dim2);
2679
77
      grid[i][j] = isl_map_empty(dim);
2680
77
    }
2681
75
  }
2682
74
2683
167
  
for (k = 0; 74
k < n167
;
++k93
)
{93
2684
93
    i = group[2 * k];
2685
93
    j = group[2 * k + 1];
2686
93
    grid[i][j] = isl_map_union(grid[i][j],
2687
93
        isl_map_from_basic_map(
2688
93
          isl_basic_map_copy(list[k])));
2689
93
  }
2690
74
  
2691
74
  floyd_warshall_iterate(grid, n_group, exact);
2692
74
2693
74
  app = isl_union_map_empty(isl_map_get_space(grid[0][0]));
2694
74
2695
149
  for (i = 0; 
i < n_group149
;
++i75
)
{75
2696
152
    for (j = 0; 
j < n_group152
;
++j77
)
2697
77
      app = isl_union_map_add_map(app, grid[i][j]);
2698
75
    free(grid[i]);
2699
75
  }
2700
74
  free(grid);
2701
74
2702
260
  for (i = 0; 
i < 2 * n260
;
++i186
)
2703
186
    isl_set_free(set[i]);
2704
74
  free(set);
2705
74
2706
74
  free(group);
2707
74
  return app;
2708
0
error:
2709
0
  if (grid)
2710
0
    
for (i = 0; 0
i < n_group0
;
++i0
)
{0
2711
0
      if (!grid[i])
2712
0
        continue;
2713
0
      
for (j = 0; 0
j < n_group0
;
++j0
)
2714
0
        isl_map_free(grid[i][j]);
2715
0
      free(grid[i]);
2716
0
    }
2717
0
  free(grid);
2718
0
  if (
set0
)
{0
2719
0
    for (i = 0; 
i < 2 * n0
;
++i0
)
2720
0
      isl_set_free(set[i]);
2721
0
    free(set);
2722
0
  }
2723
0
  free(group);
2724
0
  return NULL;
2725
74
}
2726
2727
/* Perform Floyd-Warshall on the given union relation.
2728
 * The implementation is very similar to that for non-unions.
2729
 * The main difference is that it is applied unconditionally.
2730
 * We first extract a list of basic maps from the union map
2731
 * and then perform the algorithm on this list.
2732
 */
2733
static __isl_give isl_union_map *union_floyd_warshall(
2734
  __isl_take isl_union_map *umap, int *exact)
2735
74
{
2736
74
  int i, n;
2737
74
  isl_ctx *ctx;
2738
74
  isl_basic_map **list = NULL;
2739
74
  isl_basic_map **next;
2740
74
  isl_union_map *res;
2741
74
2742
74
  n = 0;
2743
74
  if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
2744
0
    goto error;
2745
74
2746
74
  ctx = isl_union_map_get_ctx(umap);
2747
74
  list = isl_calloc_array(ctx, isl_basic_map *, n);
2748
74
  if (!list)
2749
0
    goto error;
2750
74
2751
74
  next = list;
2752
74
  if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
2753
0
    goto error;
2754
74
2755
74
  res = union_floyd_warshall_on_list(ctx, list, n, exact);
2756
74
2757
74
  if (
list74
)
{74
2758
167
    for (i = 0; 
i < n167
;
++i93
)
2759
93
      isl_basic_map_free(list[i]);
2760
74
    free(list);
2761
74
  }
2762
74
2763
74
  isl_union_map_free(umap);
2764
74
  return res;
2765
0
error:
2766
0
  if (
list0
)
{0
2767
0
    for (i = 0; 
i < n0
;
++i0
)
2768
0
      isl_basic_map_free(list[i]);
2769
0
    free(list);
2770
0
  }
2771
0
  isl_union_map_free(umap);
2772
0
  return NULL;
2773
74
}
2774
2775
/* Decompose the give union relation into strongly connected components.
2776
 * The implementation is essentially the same as that of
2777
 * construct_power_components with the major difference that all
2778
 * operations are performed on union maps.
2779
 */
2780
static __isl_give isl_union_map *union_components(
2781
  __isl_take isl_union_map *umap, int *exact)
2782
52
{
2783
52
  int i;
2784
52
  int n;
2785
52
  isl_ctx *ctx;
2786
52
  isl_basic_map **list = NULL;
2787
52
  isl_basic_map **next;
2788
52
  isl_union_map *path = NULL;
2789
52
  struct isl_tc_follows_data data;
2790
52
  struct isl_tarjan_graph *g = NULL;
2791
52
  int c, l;
2792
52
  int recheck = 0;
2793
52
2794
52
  n = 0;
2795
52
  if (isl_union_map_foreach_map(umap, inc_count, &n) < 0)
2796
0
    goto error;
2797
52
2798
52
  
if (52
n == 052
)
2799
0
    return umap;
2800
52
  
if (52
n <= 152
)
2801
40
    return union_floyd_warshall(umap, exact);
2802
52
2803
12
  ctx = isl_union_map_get_ctx(umap);
2804
12
  list = isl_calloc_array(ctx, isl_basic_map *, n);
2805
12
  if (!list)
2806
0
    goto error;
2807
12
2808
12
  next = list;
2809
12
  if (isl_union_map_foreach_map(umap, collect_basic_map, &next) < 0)
2810
0
    goto error;
2811
12
2812
12
  data.list = list;
2813
12
  data.check_closed = 0;
2814
12
  g = isl_tarjan_graph_init(ctx, n, &basic_map_follows, &data);
2815
12
  if (!g)
2816
0
    goto error;
2817
12
2818
12
  c = 0;
2819
12
  i = 0;
2820
12
  l = n;
2821
12
  path = isl_union_map_empty(isl_union_map_get_space(umap));
2822
46
  while (
l46
)
{34
2823
34
    isl_union_map *comp;
2824
34
    isl_union_map *path_comp, *path_comb;
2825
34
    comp = isl_union_map_empty(isl_union_map_get_space(umap));
2826
87
    while (
g->order[i] != -187
)
{53
2827
53
      comp = isl_union_map_add_map(comp,
2828
53
            isl_map_from_basic_map(
2829
53
          isl_basic_map_copy(list[g->order[i]])));
2830
53
      --l;
2831
53
      ++i;
2832
53
    }
2833
34
    path_comp = union_floyd_warshall(comp, exact);
2834
34
    path_comb = isl_union_map_apply_range(isl_union_map_copy(path),
2835
34
            isl_union_map_copy(path_comp));
2836
34
    path = isl_union_map_union(path, path_comp);
2837
34
    path = isl_union_map_union(path, path_comb);
2838
34
    ++i;
2839
34
    ++c;
2840
34
  }
2841
12
2842
12
  if (
c > 1 && 12
data.check_closed7
&&
!*exact0
)
{0
2843
0
    int closed;
2844
0
2845
0
    closed = isl_union_map_is_transitively_closed(path);
2846
0
    if (closed < 0)
2847
0
      goto error;
2848
0
    recheck = !closed;
2849
0
  }
2850
12
2851
12
  isl_tarjan_graph_free(g);
2852
12
2853
65
  for (i = 0; 
i < n65
;
++i53
)
2854
53
    isl_basic_map_free(list[i]);
2855
12
  free(list);
2856
12
2857
12
  if (
recheck12
)
{0
2858
0
    isl_union_map_free(path);
2859
0
    return union_floyd_warshall(umap, exact);
2860
0
  }
2861
12
2862
12
  isl_union_map_free(umap);
2863
12
2864
12
  return path;
2865
0
error:
2866
0
  isl_tarjan_graph_free(g);
2867
0
  if (
list0
)
{0
2868
0
    for (i = 0; 
i < n0
;
++i0
)
2869
0
      isl_basic_map_free(list[i]);
2870
0
    free(list);
2871
0
  }
2872
0
  isl_union_map_free(umap);
2873
0
  isl_union_map_free(path);
2874
0
  return NULL;
2875
12
}
2876
2877
/* Compute the transitive closure  of "umap", or an overapproximation.
2878
 * If the result is exact, then *exact is set to 1.
2879
 */
2880
__isl_give isl_union_map *isl_union_map_transitive_closure(
2881
  __isl_take isl_union_map *umap, int *exact)
2882
53
{
2883
53
  int closed;
2884
53
2885
53
  if (!umap)
2886
0
    return NULL;
2887
53
2888
53
  
if (53
exact53
)
2889
0
    *exact = 1;
2890
53
2891
53
  umap = isl_union_map_compute_divs(umap);
2892
53
  umap = isl_union_map_coalesce(umap);
2893
53
  closed = isl_union_map_is_transitively_closed(umap);
2894
53
  if (closed < 0)
2895
0
    goto error;
2896
53
  
if (53
closed53
)
2897
1
    return umap;
2898
52
  umap = union_components(umap, exact);
2899
52
  return umap;
2900
0
error:
2901
0
  isl_union_map_free(umap);
2902
0
  return NULL;
2903
53
}
2904
2905
struct isl_union_power {
2906
  isl_union_map *pow;
2907
  int *exact;
2908
};
2909
2910
static isl_stat power(__isl_take isl_map *map, void *user)
2911
0
{
2912
0
  struct isl_union_power *up = user;
2913
0
2914
0
  map = isl_map_power(map, up->exact);
2915
0
  up->pow = isl_union_map_from_map(map);
2916
0
2917
0
  return isl_stat_error;
2918
0
}
2919
2920
/* Construct a map [x] -> [x+1], with parameters prescribed by "dim".
2921
 */
2922
static __isl_give isl_union_map *increment(__isl_take isl_space *dim)
2923
0
{
2924
0
  int k;
2925
0
  isl_basic_map *bmap;
2926
0
2927
0
  dim = isl_space_add_dims(dim, isl_dim_in, 1);
2928
0
  dim = isl_space_add_dims(dim, isl_dim_out, 1);
2929
0
  bmap = isl_basic_map_alloc_space(dim, 0, 1, 0);
2930
0
  k = isl_basic_map_alloc_equality(bmap);
2931
0
  if (k < 0)
2932
0
    goto error;
2933
0
  isl_seq_clr(bmap->eq[k], isl_basic_map_total_dim(bmap));
2934
0
  isl_int_set_si(bmap->eq[k][0], 1);
2935
0
  isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_in)], 1);
2936
0
  isl_int_set_si(bmap->eq[k][isl_basic_map_offset(bmap, isl_dim_out)], -1);
2937
0
  return isl_union_map_from_map(isl_map_from_basic_map(bmap));
2938
0
error:
2939
0
  isl_basic_map_free(bmap);
2940
0
  return NULL;
2941
0
}
2942
2943
/* Construct a map [[x]->[y]] -> [y-x], with parameters prescribed by "dim".
2944
 */
2945
static __isl_give isl_union_map *deltas_map(__isl_take isl_space *dim)
2946
0
{
2947
0
  isl_basic_map *bmap;
2948
0
2949
0
  dim = isl_space_add_dims(dim, isl_dim_in, 1);
2950
0
  dim = isl_space_add_dims(dim, isl_dim_out, 1);
2951
0
  bmap = isl_basic_map_universe(dim);
2952
0
  bmap = isl_basic_map_deltas_map(bmap);
2953
0
2954
0
  return isl_union_map_from_map(isl_map_from_basic_map(bmap));
2955
0
}
2956
2957
/* Compute the positive powers of "map", or an overapproximation.
2958
 * The result maps the exponent to a nested copy of the corresponding power.
2959
 * If the result is exact, then *exact is set to 1.
2960
 */
2961
__isl_give isl_union_map *isl_union_map_power(__isl_take isl_union_map *umap,
2962
  int *exact)
2963
0
{
2964
0
  int n;
2965
0
  isl_union_map *inc;
2966
0
  isl_union_map *dm;
2967
0
2968
0
  if (!umap)
2969
0
    return NULL;
2970
0
  n = isl_union_map_n_map(umap);
2971
0
  if (n == 0)
2972
0
    return umap;
2973
0
  
if (0
n == 10
)
{0
2974
0
    struct isl_union_power up = { NULL, exact };
2975
0
    isl_union_map_foreach_map(umap, &power, &up);
2976
0
    isl_union_map_free(umap);
2977
0
    return up.pow;
2978
0
  }
2979
0
  inc = increment(isl_union_map_get_space(umap));
2980
0
  umap = isl_union_map_product(inc, umap);
2981
0
  umap = isl_union_map_transitive_closure(umap, exact);
2982
0
  umap = isl_union_map_zip(umap);
2983
0
  dm = deltas_map(isl_union_map_get_space(umap));
2984
0
  umap = isl_union_map_apply_domain(umap, dm);
2985
0
  
2986
0
  return umap;
2987
0
}
2988
2989
#undef TYPE
2990
1
#define TYPE isl_map
2991
#include "isl_power_templ.c"
2992
2993
#undef TYPE
2994
0
#define TYPE isl_union_map
2995
#include "isl_power_templ.c"