1 #include <assert.h>
2 #include <unistd.h>
3 #include <stdlib.h>
4 #include <gmp.h>
5 #include <isl/ctx.h>
6 #include <isl/val.h>
7 #include <isl/space.h>
8 #include <isl/point.h>
9 #include <isl/set.h>
10 #include <isl/polynomial.h>
11 #include <isl/printer.h>
12 #include <isl_set_polylib.h>
13 #include <barvinok/evalue.h>
14 #include <barvinok/util.h>
15 #include <barvinok/barvinok.h>
16 #include "barvinok_enumerate_options.h"
17 #include "verify.h"
18 #include "verify_series.h"
19 #include "remove_equalities.h"
20 #include "evalue_convert.h"
21 #include "conversion.h"
22 #include "skewed_genfun.h"
23
24 #undef CS /* for Solaris 10 */
25
26 using std::cout;
27 using std::endl;
28
29 /* The input of this example program is the same as that of testehrhart
30 * in the PolyLib distribution, i.e., a polytope in combined
31 * data and parameter space, a context polytope in parameter space
32 * and (optionally) the names of the parameters.
33 * Both polytopes are in PolyLib notation.
34 */
35
36 struct verify_point_enum {
37 struct verify_point_data vpd;
38 isl_set *set;
39 isl_pw_qpolynomial *pwqp;
40 };
41
verify_point(__isl_take isl_point * pnt,void * user)42 static isl_stat verify_point(__isl_take isl_point *pnt, void *user)
43 {
44 struct verify_point_enum *vpe = (struct verify_point_enum *) user;
45 isl_set *set;
46 int i;
47 unsigned nparam;
48 isl_val *v, *n, *t;
49 int pa = vpe->vpd.options->barvinok->approx->approximation;
50 int ok;
51 FILE *out = vpe->vpd.options->print_all ? stdout : stderr;
52
53 vpe->vpd.n--;
54
55 set = isl_set_copy(vpe->set);
56 nparam = isl_set_dim(set, isl_dim_param);
57 for (i = 0; i < nparam; ++i) {
58 v = isl_point_get_coordinate_val(pnt, isl_dim_param, i);
59 set = isl_set_fix_val(set, isl_dim_param, i, v);
60 }
61
62 v = isl_set_count_val(set);
63
64 n = isl_pw_qpolynomial_eval(isl_pw_qpolynomial_copy(vpe->pwqp),
65 isl_point_copy(pnt));
66
67 if (pa == BV_APPROX_SIGN_LOWER)
68 n = isl_val_ceil(n);
69 else if (pa == BV_APPROX_SIGN_UPPER)
70 n = isl_val_floor(n);
71 else
72 n = isl_val_trunc(n);
73
74 if (pa == BV_APPROX_SIGN_APPROX)
75 /* just accept everything */
76 ok = 1;
77 else if (pa == BV_APPROX_SIGN_LOWER)
78 ok = isl_val_le(n, v);
79 else if (pa == BV_APPROX_SIGN_UPPER)
80 ok = isl_val_ge(n, v);
81 else
82 ok = isl_val_eq(n, v);
83
84 if (vpe->vpd.options->print_all || !ok) {
85 isl_ctx *ctx = isl_point_get_ctx(pnt);
86 isl_printer *p;
87 p = isl_printer_to_file(ctx, out);
88 p = isl_printer_print_str(p, "EP(");
89 for (i = 0; i < nparam; ++i) {
90 if (i)
91 p = isl_printer_print_str(p, ", ");
92 t = isl_point_get_coordinate_val(pnt, isl_dim_param, i);
93 p = isl_printer_print_val(p, t);
94 isl_val_free(t);
95 }
96 p = isl_printer_print_str(p, ") = ");
97 p = isl_printer_print_val(p, n);
98 p = isl_printer_print_str(p, ", count = ");
99 p = isl_printer_print_val(p, v);
100 if (ok)
101 p = isl_printer_print_str(p, ". OK");
102 else
103 p = isl_printer_print_str(p, ". NOT OK");
104 p = isl_printer_end_line(p);
105 isl_printer_free(p);
106 } else if ((vpe->vpd.n % vpe->vpd.s) == 0) {
107 printf("o");
108 fflush(stdout);
109 }
110
111 isl_set_free(set);
112 isl_val_free(v);
113 isl_val_free(n);
114 isl_point_free(pnt);
115
116 if (!ok)
117 vpe->vpd.error = 1;
118
119 if (vpe->vpd.options->continue_on_error)
120 ok = 1;
121
122 return (vpe->vpd.n >= 1 && ok) ? isl_stat_ok : isl_stat_error;
123 }
124
verify_isl(Polyhedron * P,Polyhedron * C,evalue * EP,const struct verify_options * options)125 static int verify_isl(Polyhedron *P, Polyhedron *C,
126 evalue *EP, const struct verify_options *options)
127 {
128 struct verify_point_enum vpe = { { options } };
129 int i;
130 isl_ctx *ctx = isl_ctx_alloc();
131 isl_space *dim;
132 isl_set *set;
133 isl_set *set_C;
134 int r;
135
136 dim = isl_space_set_alloc(ctx, C->Dimension, P->Dimension - C->Dimension);
137 for (i = 0; i < C->Dimension; ++i)
138 dim = isl_space_set_dim_name(dim, isl_dim_param, i, options->params[i]);
139 set = isl_set_new_from_polylib(P, isl_space_copy(dim));
140 dim = isl_space_params(dim);
141 set_C = isl_set_new_from_polylib(C, dim);
142 set_C = isl_set_intersect_params(isl_set_copy(set), set_C);
143 set_C = isl_set_params(set_C);
144
145 set_C = verify_context_set_bounds(set_C, options);
146
147 r = verify_point_data_init(&vpe.vpd, set_C);
148
149 vpe.set = set;
150 vpe.pwqp = isl_pw_qpolynomial_from_evalue(isl_set_get_space(set_C), EP);
151 if (r == 0)
152 isl_set_foreach_point(set_C, verify_point, &vpe);
153 if (vpe.vpd.error)
154 r = -1;
155
156 isl_pw_qpolynomial_free(vpe.pwqp);
157 isl_set_free(set);
158 isl_set_free(set_C);
159
160 isl_ctx_free(ctx);
161
162 verify_point_data_fini(&vpe.vpd);
163
164 return r;
165 }
166
verify(Polyhedron * P,Polyhedron * C,evalue * EP,skewed_gen_fun * gf,struct enumerate_options * options)167 static int verify(Polyhedron *P, Polyhedron *C, evalue *EP, skewed_gen_fun *gf,
168 struct enumerate_options *options)
169 {
170 Polyhedron *CS, *S;
171 Vector *p;
172 int result = 0;
173
174 if (!options->series || options->function)
175 return verify_isl(P, C, EP, options->verify);
176
177 CS = check_poly_context_scan(P, &C, C->Dimension, options->verify);
178
179 p = Vector_Alloc(P->Dimension+2);
180 value_set_si(p->p[P->Dimension+1], 1);
181
182 /* S = scanning list of polyhedra */
183 S = Polyhedron_Scan(P, C, options->verify->barvinok->MaxRays);
184
185 check_poly_init(C, options->verify);
186
187 /******* CHECK NOW *********/
188 if (S) {
189 if (!check_poly_gf(S, CS, gf, 0, C->Dimension, 0, p->p,
190 options->verify))
191 result = -1;
192 Domain_Free(S);
193 }
194
195 if (result == -1)
196 fprintf(stderr,"Check failed !\n");
197
198 if (!options->verify->print_all)
199 printf( "\n" );
200
201 Vector_Free(p);
202 if (CS) {
203 Domain_Free(CS);
204 Domain_Free(C);
205 }
206
207 return result;
208 }
209
210 /* frees M and Minv */
apply_transformation(Polyhedron ** P,Polyhedron ** C,bool free_P,bool free_C,Matrix * M,Matrix * Minv,Matrix ** inv,barvinok_options * options)211 static void apply_transformation(Polyhedron **P, Polyhedron **C,
212 bool free_P, bool free_C,
213 Matrix *M, Matrix *Minv, Matrix **inv,
214 barvinok_options *options)
215 {
216 Polyhedron *T;
217 Matrix *M2;
218
219 M2 = align_matrix(M, (*P)->Dimension + 1);
220 T = *P;
221 *P = Polyhedron_Preimage(*P, M2, options->MaxRays);
222 if (free_P)
223 Polyhedron_Free(T);
224 Matrix_Free(M2);
225
226 T = *C;
227 *C = Polyhedron_Preimage(*C, M, options->MaxRays);
228 if (free_C)
229 Polyhedron_Free(T);
230
231 Matrix_Free(M);
232
233 if (*inv) {
234 Matrix *T = *inv;
235 *inv = Matrix_Alloc(Minv->NbRows, T->NbColumns);
236 Matrix_Product(Minv, T, *inv);
237 Matrix_Free(T);
238 Matrix_Free(Minv);
239 } else
240 *inv = Minv;
241 }
242
243 /* Since we have "compressed" the parameters (in case there were
244 * any equalities), the result is independent of the coordinates in the
245 * coordinate subspace spanned by the lines. We can therefore assume
246 * these coordinates are zero and compute the inverse image of the map
247 * from a lower dimensional space that adds zeros in the appropriate
248 * places.
249 */
remove_lines(Polyhedron * C,Matrix ** M,Matrix ** Minv)250 static void remove_lines(Polyhedron *C, Matrix **M, Matrix **Minv)
251 {
252 Matrix *L = Matrix_Alloc(C->Dimension+1, C->Dimension+1);
253 for (int r = 0; r < C->NbBid; ++r)
254 Vector_Copy(C->Ray[r]+1, L->p[r], C->Dimension);
255 unimodular_complete(L, C->NbBid);
256 assert(value_one_p(L->p[C->Dimension][C->Dimension]));
257 assert(First_Non_Zero(L->p[C->Dimension], C->Dimension) == -1);
258 Matrix_Transposition(L);
259 assert(First_Non_Zero(L->p[C->Dimension], C->Dimension) == -1);
260
261 *M = Matrix_Alloc(C->Dimension+1, C->Dimension-C->NbBid+1);
262 for (int i = 0; i < C->Dimension+1; ++i)
263 Vector_Copy(L->p[i]+C->NbBid, (*M)->p[i], C->Dimension-C->NbBid+1);
264
265 Matrix *Linv = Matrix_Alloc(C->Dimension+1, C->Dimension+1);
266 int ok = Matrix_Inverse(L, Linv);
267 assert(ok);
268 Matrix_Free(L);
269
270 *Minv = Matrix_Alloc(C->Dimension-C->NbBid+1, C->Dimension+1);
271 for (int i = C->NbBid; i < C->Dimension+1; ++i)
272 Vector_AntiScale(Linv->p[i], (*Minv)->p[i-C->NbBid],
273 Linv->p[C->Dimension][C->Dimension], C->Dimension+1);
274 Matrix_Free(Linv);
275 }
276
series(Polyhedron * P,Polyhedron * C,barvinok_options * options)277 static skewed_gen_fun *series(Polyhedron *P, Polyhedron* C,
278 barvinok_options *options)
279 {
280 Polyhedron *C1, *C2;
281 gen_fun *gf;
282 Matrix *inv = NULL;
283 Matrix *eq = NULL;
284 Matrix *div = NULL;
285 Polyhedron *PT = P;
286
287 /* Compute true context */
288 C1 = Polyhedron_Project(P, C->Dimension);
289 C2 = DomainIntersection(C, C1, options->MaxRays);
290 Polyhedron_Free(C1);
291
292 POL_ENSURE_VERTICES(C2);
293 if (C2->NbBid != 0) {
294 Matrix *CP;
295 if (C2->NbEq || P->NbEq) {
296 /* We remove all equalities to be sure all lines are unit vectors */
297 Polyhedron *CT = C2;
298 remove_all_equalities(&PT, &CT, &CP, NULL, C2->Dimension,
299 options->MaxRays);
300 if (CT != C2) {
301 Polyhedron_Free(C2);
302 C2 = CT;
303 }
304 if (CP) {
305 inv = left_inverse(CP, &eq);
306 Matrix_Free(CP);
307
308 int d = 0;
309 Value tmp;
310 value_init(tmp);
311 div = Matrix_Alloc(inv->NbRows-1, inv->NbColumns+1);
312 for (int i = 0; i < inv->NbRows-1; ++i) {
313 Vector_Gcd(inv->p[i], inv->NbColumns, &tmp);
314 if (mpz_divisible_p(tmp,
315 inv->p[inv->NbRows-1][inv->NbColumns-1]))
316 continue;
317 Vector_Copy(inv->p[i], div->p[d], inv->NbColumns);
318 value_assign(div->p[d][inv->NbColumns],
319 inv->p[inv->NbRows-1][inv->NbColumns-1]);
320 ++d;
321 }
322 value_clear(tmp);
323
324 if (!d) {
325 Matrix_Free(div);
326 div = NULL;
327 } else
328 div->NbRows = d;
329 }
330 }
331 POL_ENSURE_VERTICES(C2);
332
333 if (C2->NbBid) {
334 Matrix *M, *Minv;
335 remove_lines(C2, &M, &Minv);
336 apply_transformation(&PT, &C2, PT != P, C2 != C, M, Minv, &inv,
337 options);
338 }
339 }
340 POL_ENSURE_VERTICES(C2);
341 if (!Polyhedron_has_revlex_positive_rays(C2, C2->Dimension)) {
342 Matrix *Constraints;
343 Matrix *H, *Q, *U;
344 Constraints = Matrix_Alloc(C2->NbConstraints, C2->Dimension+1);
345 for (int i = 0; i < C2->NbConstraints; ++i)
346 Vector_Copy(C2->Constraint[i]+1, Constraints->p[i], C2->Dimension);
347 left_hermite(Constraints, &H, &Q, &U);
348 Matrix_Free(Constraints);
349 /* flip rows of Q */
350 for (int i = 0; i < C2->Dimension/2; ++i)
351 Vector_Exchange(Q->p[i], Q->p[C2->Dimension-1-i], C2->Dimension);
352 Matrix_Free(H);
353 Matrix_Free(U);
354 Matrix *M = Matrix_Alloc(C2->Dimension+1, C2->Dimension+1);
355 U = Matrix_Copy(Q);
356 int ok = Matrix_Inverse(U, M);
357 assert(ok);
358 Matrix_Free(U);
359
360 apply_transformation(&PT, &C2, PT != P, C2 != C, M, Q, &inv, options);
361 }
362 gf = barvinok_series_with_options(PT, C2, options);
363 Polyhedron_Free(C2);
364 if (PT != P)
365 Polyhedron_Free(PT);
366 return new skewed_gen_fun(gf, inv, eq, div);
367 }
368
main(int argc,char ** argv)369 int main(int argc, char **argv)
370 {
371 Polyhedron *A, *C;
372 Matrix *M;
373 evalue *EP = NULL;
374 skewed_gen_fun *gf = NULL;
375 const char **param_name;
376 int print_solution = 1;
377 int result = 0;
378 struct enumerate_options *options = enumerate_options_new_with_defaults();
379
380 argc = enumerate_options_parse(options, argc, argv, ISL_ARG_ALL);
381
382 M = Matrix_Read();
383 assert(M);
384 A = Constraints2Polyhedron(M, options->verify->barvinok->MaxRays);
385 Matrix_Free(M);
386 M = Matrix_Read();
387 assert(M);
388 C = Constraints2Polyhedron(M, options->verify->barvinok->MaxRays);
389 Matrix_Free(M);
390 assert(A->Dimension >= C->Dimension);
391 param_name = Read_ParamNames(stdin, C->Dimension);
392
393 if (options->verify->verify) {
394 verify_options_set_range(options->verify, A->Dimension);
395 if (!options->verify->barvinok->verbose)
396 print_solution = 0;
397 }
398
399 if (print_solution && options->verify->barvinok->verbose) {
400 Polyhedron_Print(stdout, P_VALUE_FMT, A);
401 Polyhedron_Print(stdout, P_VALUE_FMT, C);
402 }
403
404 if (options->series) {
405 gf = series(A, C, options->verify->barvinok);
406 if (print_solution) {
407 gf->print(cout, C->Dimension, param_name);
408 puts("");
409 }
410 if (options->function) {
411 EP = *gf;
412 if (print_solution)
413 print_evalue(stdout, EP, param_name);
414 }
415 } else {
416 EP = barvinok_enumerate_with_options(A, C, options->verify->barvinok);
417 assert(EP);
418 if (evalue_convert(EP, options->convert, options->verify->barvinok->verbose,
419 C->Dimension, param_name))
420 print_solution = 0;
421 if (options->size)
422 printf("\nSize: %zd\n", evalue_size(EP));
423 if (print_solution)
424 print_evalue(stdout, EP, param_name);
425 }
426
427 if (options->verify->verify) {
428 options->verify->params = param_name;
429 result = verify(A, C, EP, gf, options);
430 }
431
432 if (gf)
433 delete gf;
434 if (EP)
435 evalue_free(EP);
436
437 if (options->verify->barvinok->print_stats)
438 barvinok_stats_print(options->verify->barvinok->stats, stdout);
439
440 Free_ParamNames(param_name, C->Dimension);
441 Polyhedron_Free(A);
442 Polyhedron_Free(C);
443 enumerate_options_free(options);
444 return result;
445 }
446