doc: add another paper refering to the library
[barvinok.git] / remove_equalities.c
blob378d77dae75144e6ea4c0a7422687f94a3c0ba15
1 #include <assert.h>
2 #include <barvinok/util.h>
3 #include "remove_equalities.h"
5 static void transform(Polyhedron **P, Polyhedron **C, Matrix *CP, int free,
6 unsigned MaxRays)
8 Polyhedron *Q = *P;
9 Polyhedron *D = *C;
10 Matrix *T;
12 T = align_matrix(CP, Q->Dimension+1);
13 *P = Polyhedron_Preimage(Q, T, MaxRays);
14 if (free)
15 Polyhedron_Free(Q);
16 Matrix_Free(T);
17 if (!D)
18 return;
19 *C = Polyhedron_Preimage(D, CP, MaxRays);
20 if (free)
21 Polyhedron_Free(D);
24 static Matrix *compose_transformations(Matrix *first, Matrix *second)
26 Matrix *combined;
28 if (!first)
29 return second;
31 combined = Matrix_Alloc(first->NbRows, second->NbColumns);
32 Matrix_Product(first, second, combined);
33 Matrix_Free(first);
34 Matrix_Free(second);
35 return combined;
38 static Polyhedron *replace_by_empty_polyhedron(Polyhedron *Q, int free)
40 unsigned dim = Q->Dimension;
41 if (free)
42 Polyhedron_Free(Q);
43 return Empty_Polyhedron(dim);
46 static int first_parameter_equality(Polyhedron *Q, unsigned nparam)
48 int i;
50 if (emptyQ(Q))
51 return Q->NbEq;
53 for (i = 0; i < Q->NbEq; ++i)
54 if (First_Non_Zero(Q->Constraint[i]+1, Q->Dimension-nparam) == -1)
55 break;
57 return i;
60 static void remove_parameter_equalities(Polyhedron **Q, Polyhedron **D,
61 Matrix **CP, unsigned *nparam,
62 int free, unsigned MaxRays)
64 int i;
66 /* We need to loop until we can't find any more equalities
67 * because the transformation may enable a simplification of the
68 * constraints resulting in new implicit equalities.
70 while ((i = first_parameter_equality(*Q, *nparam)) < (*Q)->NbEq) {
71 Matrix *M = Matrix_Alloc((*Q)->NbEq, 1+*nparam+1);
72 Matrix *T;
73 int n = 0;
75 for (; i < (*Q)->NbEq; ++i) {
76 if (First_Non_Zero((*Q)->Constraint[i]+1, (*Q)->Dimension-*nparam) == -1)
77 Vector_Copy((*Q)->Constraint[i]+1+(*Q)->Dimension-*nparam,
78 M->p[n++]+1, *nparam+1);
80 M->NbRows = n;
81 T = compress_variables(M, 0);
82 Matrix_Free(M);
83 if (!T) {
84 *Q = replace_by_empty_polyhedron(*Q, free);
85 break;
86 } else {
87 transform(Q, D, T, free, MaxRays);
88 *nparam = T->NbColumns-1;
89 *CP = compose_transformations(*CP, T);
90 free = 1;
95 static int is_translation(Matrix *M)
97 unsigned i, j;
99 if (M->NbRows != M->NbColumns)
100 return 0;
102 for (i = 0; i < M->NbRows; ++i)
103 for (j = 0; j < M->NbColumns-1; ++j)
104 if (i == j ? value_notone_p(M->p[i][j])
105 : value_notzero_p(M->p[i][j]))
106 return 0;
108 return 1;
111 /* Remove all equalities in P and the context C (if not NULL).
112 * Does not destroy P (or C).
113 * Returns transformation on the parameters in the Matrix pointed to by CPP
114 * (unless NULL) and transformation on the variables in the Matrix pointed to
115 * by CVP (unless NULL).
116 * Each of these transformation matrices maps the new parameters/variables
117 * back to the original ones.
118 * If it can be shown that there is no integer point in P, then
119 * an empty polyhedron will be returned.
121 int remove_all_equalities(Polyhedron **P, Polyhedron **C, Matrix **CPP, Matrix **CVP,
122 unsigned nparam, unsigned MaxRays)
124 Matrix *CV = NULL;
125 Matrix *CP = NULL;
126 Polyhedron *Q = *P;
127 Polyhedron *D = NULL;
128 Polyhedron *R;
129 int changed;
130 Matrix M;
132 if (C) {
133 D = *C;
134 assert(D->Dimension == nparam);
137 if (Q->NbEq == 0 && (!D || D->NbEq == 0))
138 return 0;
140 if (D && D->NbEq) {
141 Polyhedron_Matrix_View(D, &M, D->NbEq);
142 CP = compress_variables(&M, 0);
143 if (!CP) {
144 *C = replace_by_empty_polyhedron(D, D != *C);
145 return 1;
147 transform(&Q, &D, CP, Q != *P, MaxRays);
148 nparam = CP->NbColumns-1;
151 /* compress_parms doesn't like equalities that only involve parameters */
152 remove_parameter_equalities(&Q, &D, &CP, &nparam, Q != *P, MaxRays);
154 /* We need to loop until we can't find any more equalities
155 * because the transformation may enable a simplification of the
156 * constraints resulting in new implicit equalities.
158 while (!emptyQ2(Q) && Q->NbEq) {
159 Matrix *T;
161 do {
162 Polyhedron_Matrix_View(Q, &M, Q->NbEq);
163 T = compress_parms(&M, nparam);
164 if (!T) {
165 Q = replace_by_empty_polyhedron(Q, Q != *P);
166 } else if (is_translation(T)) {
167 Matrix_Free(T);
168 T = NULL;
169 } else {
170 transform(&Q, &D, T, Q != *P, MaxRays);
171 CP = compose_transformations(CP, T);
172 nparam = CP->NbColumns-1;
173 remove_parameter_equalities(&Q, &D, &CP, &nparam, Q != *P, MaxRays);
175 } while (!emptyQ(Q) && Q->NbEq && T);
177 if (emptyQ(Q) || !Q->NbEq)
178 break;
180 Polyhedron_Matrix_View(Q, &M, Q->NbEq);
181 T = compress_variables(&M, nparam);
182 if (!T) {
183 Q = replace_by_empty_polyhedron(Q, Q != *P);
184 } else if (isIdentity(T)) {
185 Matrix_Free(T);
186 } else {
187 R = Polyhedron_Preimage(Q, T, MaxRays);
188 CV = compose_transformations(CV, T);
189 if (Q != *P)
190 Polyhedron_Free(Q);
191 Q = R;
195 changed = *P != Q;
196 *P = Q;
197 if (C)
198 *C = D;
199 if (CPP)
200 *CPP = CP;
201 else if (CP)
202 Matrix_Free(CP);
203 if (CVP)
204 *CVP = CV;
205 else if (CV)
206 Matrix_Free(CV);
207 return changed;