isl_pw_qpolynomial_sum: handle existentials in wrapped domains
[barvinok.git] / normalization.c
blob22075be15974db9d267720f7af27aac02fee6541
1 #include <assert.h>
2 #include <barvinok/util.h>
3 #include "normalization.h"
5 static int is_unit_row(Value *row, int pos, int len)
7 if (!value_one_p(row[pos]) && !value_mone_p(row[pos]))
8 return 0;
9 return First_Non_Zero(row+pos+1, len-(pos+1)) == -1;
12 /* Transform the constraints C to "standard form".
13 * In particular, if C is described by
14 * A x + b(p) >= 0
15 * then this function returns a matrix H = A U, A = H Q, such
16 * that D x' = D Q x >= -b(p), with D a diagonal matrix with
17 * positive entries. The calling function can then construct
18 * the standard form H' x' - I s + b(p) = 0, with H' the rows of H
19 * that are not positive multiples of unit vectors
20 * (since those correspond to D x' >= -b(p)).
21 * The number of rows in H' is returned in *rows_p.
22 * Optionally returns the matrix that maps the new variables
23 * back to the old variables x = U x'.
24 * Note that the rows of H (and C) may be reordered by this function.
26 Matrix *standard_constraints(Matrix *C, unsigned nparam, int *rows_p,
27 Matrix **T)
29 unsigned dim = C->NbColumns - 2;
30 unsigned nvar = dim - nparam;
31 int i, j, d;
32 int rows;
33 Matrix *M;
34 Matrix *H, *U, *Q;
36 /* move constraints only involving parameters down
37 * and move unit vectors (if there are any) to the right place.
39 for (d = 0, j = C->NbRows; d < j; ++d) {
40 int index;
41 index = First_Non_Zero(C->p[d]+1, nvar);
42 if (index != -1) {
43 if (index != d &&
44 is_unit_row(C->p[d]+1, index, nvar)) {
45 Vector_Exchange(C->p[d], C->p[index], dim+2);
46 if (index > d)
47 --d;
49 continue;
51 while (d < --j && First_Non_Zero(C->p[j]+1, nvar) == -1)
53 if (d >= j)
54 break;
55 Vector_Exchange(C->p[d], C->p[j], dim+2);
57 M = Matrix_Alloc(d, nvar);
58 for (j = 0; j < d; ++j)
59 Vector_Copy(C->p[j]+1, M->p[j], nvar);
61 neg_left_hermite(M, &H, &Q, &U);
62 Matrix_Free(M);
63 Matrix_Free(Q);
64 if (T)
65 *T = U;
66 else
67 Matrix_Free(U);
69 /* Rearrange rows such that top of H is lower diagonal and
70 * compute the number of non (multiple of) unit-vector rows.
72 rows = H->NbRows-nvar;
73 for (i = 0; i < H->NbColumns; ++i) {
74 for (j = i; j < H->NbRows; ++j)
75 if (value_notzero_p(H->p[j][i]))
76 break;
77 if (j != i) {
78 Vector_Exchange(C->p[i], C->p[j], dim+2);
79 Vector_Exchange(H->p[i], H->p[j], H->NbColumns);
81 if (First_Non_Zero(H->p[i], i) != -1)
82 rows++;
84 if (rows_p)
85 *rows_p = rows;
87 return H;
91 * Transform the constraints of PE involving existentially
92 * quantified variables, placed in front in PE, to "standard form".
94 * The resulting matrix has the existential variables after
95 * the unquantified variables again.
97 static Matrix *exist_standard_constraints(Polyhedron *PE, unsigned nvar)
99 int i, i1, i2;
100 unsigned exist = PE->Dimension - nvar;
101 Matrix *M = Matrix_Alloc(PE->NbConstraints, 2+PE->Dimension);
102 Matrix *M2 = Matrix_Alloc(PE->NbConstraints, 2+PE->Dimension);
103 Matrix *S2;
104 Value c;
106 for (i = i1 = i2 = 0; i < PE->NbConstraints; ++i) {
107 if (First_Non_Zero(PE->Constraint[i]+1, exist) == -1) {
108 value_assign(M->p[i1][0], PE->Constraint[i][0]);
109 Vector_Copy(PE->Constraint[i]+1+exist, M->p[i1]+1, nvar);
110 value_assign(M->p[i1][1+PE->Dimension],
111 PE->Constraint[i][1+PE->Dimension]);
112 ++i1;
113 } else {
114 Vector_Copy(PE->Constraint[i], M2->p[i2], M2->NbColumns);
115 ++i2;
118 M->NbRows = i1;
119 M2->NbRows = i2;
121 S2 = standard_constraints(M2, nvar, NULL, NULL);
123 for (i = 0; i < M2->NbRows; ++i) {
124 value_assign(M->p[M->NbRows+i][0], M2->p[i][0]);
125 Vector_Copy(M2->p[i]+1+exist, M->p[M->NbRows+i]+1, nvar);
126 Vector_Copy(S2->p[i], M->p[M->NbRows+i]+1+nvar, exist);
127 value_assign(M->p[M->NbRows+i][1+PE->Dimension],
128 M2->p[i][1+PE->Dimension]);
131 /* Make sure "standard" constraints for existential variables
132 * have only non-positive coefficients for the variables.
134 assert(M2->NbRows >= exist);
135 value_init(c);
136 for (i = 0; i < exist; ++i) {
137 assert(value_pos_p(S2->p[i][i]));
138 for (i1 = 0; i1 < nvar; ++i1) {
139 if (value_negz_p(M->p[M->NbRows+i][1+i1]))
140 continue;
141 mpz_cdiv_q(c, M->p[M->NbRows+i][1+i1], S2->p[i][i]);
142 value_oppose(c, c);
143 for (i2 = i; i2 < M2->NbRows; ++i2)
144 value_addmul(M->p[M->NbRows+i2][1+i1], c, S2->p[i2][i]);
147 value_clear(c);
148 M->NbRows += M2->NbRows;
150 Matrix_Free(M2);
152 return M;
156 * For those variables x that may attain negative values,
157 * compute a translations t such that x' = x + t is sure
158 * to attain only non-negative values, i.e., if M is
160 * A x + b >= 0
162 * then the x' in
164 * A x' + b - A t >= 0
166 * are sure to attain only non-negative values.
168 * The first nvar shifts have already been set prior to
169 * calling this function.
171 static void compute_shifts(Matrix *M, unsigned nvar, Vector *shifts)
173 int c, r;
174 unsigned dim = M->NbColumns-2;
175 for (c = nvar, r = 0; c < dim; ++c) {
176 for ( ; r < M->NbRows; ++r) {
177 if (value_notzero_p(M->p[r][1+c]))
178 break;
180 assert(r < M->NbRows);
181 assert(value_pos_p(M->p[r][1+c]));
182 if (c > 0)
183 Inner_Product(M->p[r]+1, shifts->p+1, c, &shifts->p[1+c]);
184 value_subtract(shifts->p[1+c], M->p[r][1+dim], shifts->p[1+c]);
185 if (value_pos_p(shifts->p[1+c]))
186 mpz_cdiv_q(shifts->p[1+c], shifts->p[1+c], M->p[r][1+c]);
187 else
188 value_set_si(shifts->p[1+c], 0);
193 * Move the existential variables before the unquantified variables.
195 static Polyhedron *move_exists_in_front(Polyhedron *D, unsigned nvar,
196 unsigned MaxRays)
198 Polyhedron *P;
199 Polyhedron *DE = NULL;
200 Polyhedron **next = &DE;
202 for (P = D; P; P = P->next) {
203 int i;
204 unsigned exist = P->Dimension - nvar;
205 Matrix *M = Matrix_Alloc(P->NbConstraints, 2+P->Dimension);
207 for (i = 0; i < P->NbConstraints; ++i) {
208 value_assign(M->p[i][0], P->Constraint[i][0]);
209 Vector_Copy(P->Constraint[i]+1, M->p[i]+1+exist, nvar);
210 Vector_Copy(P->Constraint[i]+1+nvar, M->p[i]+1, exist);
211 value_assign(M->p[i][1+P->Dimension],
212 P->Constraint[i][1+P->Dimension]);
215 *next = Constraints2Polyhedron(M, MaxRays);
216 Matrix_Free(M);
218 next = &(*next)->next;
221 return DE;
225 * Compute the shifts s for the variables x, such that x + s >= 0.
226 * DE is a union of polyhedra with the existential variables first.
228 * We simply project out the existential variables, compute
229 * the convex hull and then compute the minimal coordinate
230 * in each direction of all the vertices.
231 * The shift is the opposite of this minimal coordinate or zero.
233 static Vector *compute_var_shifts(Polyhedron *DE, unsigned nvar,
234 unsigned MaxRays)
236 Polyhedron *PE;
237 Polyhedron *D = NULL;
238 Polyhedron **next = &D;
239 Polyhedron *C;
240 Vector *shifts;
241 int i, j;
242 Value tmp;
244 for (PE = DE; PE; PE = PE->next) {
245 *next = Polyhedron_Project(PE, nvar);
246 next = &(*next)->next;
248 C = DomainConvex(D, MaxRays);
249 Domain_Free(D);
250 POL_ENSURE_VERTICES(C);
252 shifts = Vector_Alloc(C->Dimension);
253 value_init(tmp);
254 for (i = 0; i < C->Dimension; ++i) {
255 for (j = 0; j < C->NbRays; ++j) {
256 assert(value_notzero_p(C->Ray[j][1+C->Dimension]));
257 if (value_posz_p(C->Ray[j][1+i]))
258 continue;
259 mpz_fdiv_q(tmp, C->Ray[j][1+i], C->Ray[j][1+C->Dimension]);
260 if (value_lt(tmp, shifts->p[i]))
261 value_assign(shifts->p[i], tmp);
263 value_oppose(shifts->p[i], shifts->p[i]);
265 value_clear(tmp);
266 Polyhedron_Free(C);
268 return shifts;
272 * Given a union of polyhedra, with possibly existentially
273 * quantified variables, but no parameters, apply a unimodular
274 * transformation with constant translation to the variables,
275 * and a unimodular transformation with a translation possibly
276 * depending on the variables to the existentially quantified variables,
277 * such that all variables and all existentially quantified variables
278 * attain only non-negative values.
280 * The projection onto the variables is assumed to be a polytope.
282 Polyhedron *skew_to_positive_orthant(Polyhedron *D, unsigned nvar,
283 unsigned MaxRays)
285 Polyhedron *PE;
286 Polyhedron *R = NULL;
287 Polyhedron **next = &R;
288 Vector *var_shifts;
289 Polyhedron *DE;
291 DE = move_exists_in_front(D, nvar, MaxRays);
292 var_shifts = compute_var_shifts(DE, nvar, MaxRays);
294 for (PE = DE; PE; PE = PE->next) {
295 Matrix *M, *C;
296 Vector *shifts, *offsets;
297 int i;
299 M = exist_standard_constraints(PE, nvar);
300 shifts = Vector_Alloc(M->NbColumns);
301 Vector_Copy(var_shifts->p, shifts->p+1, nvar);
302 compute_shifts(M, nvar, shifts);
304 offsets = Vector_Alloc(M->NbRows);
305 Matrix_Vector_Product(M, shifts->p, offsets->p);
306 Vector_Free(shifts);
308 C = Matrix_Alloc(M->NbRows, M->NbColumns);
309 for (i = 0; i < M->NbRows; ++i) {
310 value_assign(C->p[i][0], M->p[i][0]);
311 Vector_Copy(M->p[i]+1, C->p[i]+1, M->NbColumns-2);
312 value_subtract(C->p[i][M->NbColumns-1],
313 M->p[i][M->NbColumns-1], offsets->p[i]);
315 Vector_Free(offsets);
316 Matrix_Free(M);
318 *next = Constraints2Polyhedron(C, MaxRays);
319 Matrix_Free(C);
321 next = &(*next)->next;
323 Vector_Free(var_shifts);
324 Domain_Free(DE);
325 return R;