change arguments of normalize methods
[barvinok.git] / piputil.c
blobaf0d8117c455d8911de1b2abd548252bbc8febb8
1 #include <polylib/polylibgmp.h>
2 #include <stdio.h> /* needed for piplib ! */
3 #include <piplib/piplibMP.h>
4 #include <assert.h>
6 static PipMatrix *poly2pip(Polyhedron *P, int exist, int nparam)
8 PipMatrix * matrix;
9 int i, j;
10 unsigned int nvar = P->Dimension - exist - nparam;
12 matrix = pip_matrix_alloc(P->NbConstraints, P->Dimension+2);
13 for (i = 0; i < P->NbConstraints; ++i) {
14 value_assign(matrix->p[i][0], P->Constraint[i][0]);
15 for (j = 0; j < exist; ++j)
16 value_assign(matrix->p[i][1+j], P->Constraint[i][1+nvar+j]);
17 for (j = 0; j < nvar; ++j)
18 value_assign(matrix->p[i][1+exist+j], P->Constraint[i][1+j]);
19 for (j = 1+exist+nvar; j < P->Dimension+2; ++j)
20 value_assign(matrix->p[i][j], P->Constraint[i][j]);
23 return matrix;
26 static int max_new(PipQuast *q, int max, int d, int *maxd)
28 PipNewparm *p;
30 for (p = q->newparm; p; p = p->next)
31 if (p->rank > max)
32 max = p->rank;
33 if (q->condition) {
34 if (++d > *maxd)
35 *maxd = d;
36 max = max_new(q->next_else, max, d, maxd);
37 max = max_new(q->next_then, max, d, maxd);
39 return max;
42 static int vectorpos2cpos(int v, int n, int e, int d, int j)
44 int np = d-v-n-e;
45 int l;
47 if (j < v)
48 l = 1+j;
49 else if (j < np+v)
50 l = 1+v+n+e+j-v;
51 else
52 l = 1+v+n+j-v-np;
54 return l;
58 * i: next row
59 * n: dimension of space in which minimum was taken
60 * i.e., number of existential vars in original problem
61 * v: number of variables in original problem
62 * e: number of extra existential variables
63 * d: total dimension
65 * Dimensions in the resulting domain are:
66 * v n e (d-v-n-e)
68 static void add_quast(Polyhedron**D, Matrix* M, PipQuast *q,
69 int i, int v, int n, int e, int d)
71 PipNewparm *p;
73 for (p = q->newparm; p; p = p->next) {
74 int j, l;
75 PipVector *vec = p->vector;
77 value_set_si(M->p[i][0], 1);
78 Vector_Set(M->p[i]+1, 0, d);
80 for (j = 0; j < vec->nb_elements-1; ++j) {
81 l = vectorpos2cpos(v, n, e, d, j);
82 value_assign(M->p[i][l], vec->the_vector[j]);
84 l = vectorpos2cpos(v, n, e, d, p->rank);
85 value_oppose(M->p[i][l], p->deno);
86 value_assign(M->p[i][1+d], vec->the_vector[j]);
87 ++i;
89 value_set_si(M->p[i][0], 1);
90 Vector_Set(M->p[i]+1, 0, d);
92 for (j = 0; j < vec->nb_elements-1; ++j) {
93 l = vectorpos2cpos(v, n, e, d, j);
94 value_oppose(M->p[i][l], vec->the_vector[j]);
96 l = vectorpos2cpos(v, n, e, d, p->rank);
97 value_assign(M->p[i][l], p->deno);
98 value_oppose(M->p[i][1+d], vec->the_vector[j]);
99 value_addto(M->p[i][1+d], M->p[i][1+d], p->deno);
100 value_decrement(M->p[i][1+d], M->p[i][1+d]);
101 ++i;
104 if (q->condition) {
105 int j;
106 value_set_si(M->p[i][0], 1);
107 Vector_Set(M->p[i]+1, 0, d);
109 for (j = 0; j < q->condition->nb_elements-1; ++j) {
110 int l = vectorpos2cpos(v, n, e, d, j);
111 value_assign(M->p[i][l], q->condition->the_vector[j]);
113 value_assign(M->p[i][1+d], q->condition->the_vector[j]);
114 add_quast(D, M, q->next_then, i+1, v, n, e, d);
116 for (j = 0; j < q->condition->nb_elements-1; ++j) {
117 int l = vectorpos2cpos(v, n, e, d, j);
118 value_oppose(M->p[i][l], q->condition->the_vector[j]);
120 value_oppose(M->p[i][1+d], q->condition->the_vector[j]);
121 value_decrement(M->p[i][1+d], M->p[i][1+d]);
122 add_quast(D, M, q->next_else, i+1, v, n, e, d);
123 } else if (q->list) {
124 PipList *l;
125 Matrix *C;
126 Polyhedron *P;
127 int j, k;
128 for (j = 0, l = q->list; l; ++j, l = l->next) {
129 Vector_Set(M->p[i], 0, d+1);
130 value_set_si(M->p[i][1+v+j], -1);
132 for (k = 0; k < l->vector->nb_elements-1; ++k) {
133 int ll = vectorpos2cpos(v, n, e, d, k);
134 value_assign(M->p[i][ll], l->vector->the_vector[k]);
136 value_assign(M->p[i][1+d], l->vector->the_vector[k]);
138 ++i;
140 C = Matrix_Alloc(i, 1+d+1);
141 Vector_Copy(M->p[0], C->p[0], i * (1+d+1));
142 P = Constraints2Polyhedron(C, 0);
143 Matrix_Free(C);
144 P->next = *D;
145 *D = P;
149 static Polyhedron *quast2poly(PipQuast *q, int nvar, int nset)
151 int nparam, nexist;
152 PipList* l;
153 int i, d, rows;
154 Matrix* M;
155 Polyhedron* D;
157 assert(nset != 0); /* required ? */
158 nparam = q->condition ? q->condition->nb_elements-1 :
159 q->list->vector->nb_elements-1;
160 d = 0;
161 nexist = max_new(q, nparam-1, 0, &d) - nparam+1;
162 rows = nparam + 2 * nexist + d + nset;
163 M = Matrix_Alloc(rows, 1+nset+nexist+nparam+1);
164 /* All parameters are/should be positive */
165 for (i = 0; i < nvar; ++i) {
166 value_set_si(M->p[i][0], 1);
167 value_set_si(M->p[i][1+i], 1);
169 for ( ; i < nparam; ++i) {
170 value_set_si(M->p[i][0], 1);
171 value_set_si(M->p[i][1+nvar+nset+nexist+i-nvar], 1);
173 D = 0;
174 add_quast(&D, M, q, nparam, nvar, nset, nexist, nset+nexist+nparam);
175 Matrix_Free(M);
176 return D;
179 Polyhedron *pip_lexmin(Polyhedron *P, int exist, int nparam)
181 PipOptions *options;
182 PipMatrix * domain;
183 unsigned int nvar = P->Dimension - exist - nparam;
184 PipMatrix *context = pip_matrix_alloc(0, nvar+nparam+2);
185 PipQuast *sol;
186 Polyhedron *min;
188 domain = poly2pip(P, exist, nparam);
190 pip_matrix_print(stderr, domain);
191 pip_matrix_print(stderr, context);
194 options = pip_options_init();
195 sol = pip_solve(domain, context, -1, options);
198 pip_quast_print(stderr, sol, 0);
201 min = quast2poly(sol, nvar, exist);
203 pip_quast_free(sol);
204 pip_matrix_free(context);
205 pip_matrix_free(domain);
206 pip_options_free(options);
208 return min;