5 #include <isl_set_polylib.h>
6 #include <barvinok/evalue.h>
7 #include <barvinok/util.h>
8 #include <barvinok/barvinok.h>
9 #include "barvinok_enumerate_options.h"
11 #include "verify_series.h"
12 #include "remove_equalities.h"
13 #include "evalue_convert.h"
14 #include "conversion.h"
15 #include "skewed_genfun.h"
17 #undef CS /* for Solaris 10 */
22 /* The input of this example program is the same as that of testehrhart
23 * in the PolyLib distribution, i.e., a polytope in combined
24 * data and parameter space, a context polytope in parameter space
25 * and (optionally) the names of the parameters.
26 * Both polytopes are in PolyLib notation.
29 struct verify_point_enum
{
30 struct verify_point_data vpd
;
32 isl_pw_qpolynomial
*pwqp
;
35 static isl_stat
verify_point(__isl_take isl_point
*pnt
, void *user
)
37 struct verify_point_enum
*vpe
= (struct verify_point_enum
*) user
;
42 int pa
= vpe
->vpd
.options
->barvinok
->approx
->approximation
;
44 FILE *out
= vpe
->vpd
.options
->print_all
? stdout
: stderr
;
48 set
= isl_set_copy(vpe
->set
);
49 nparam
= isl_set_dim(set
, isl_dim_param
);
50 for (i
= 0; i
< nparam
; ++i
) {
51 v
= isl_point_get_coordinate_val(pnt
, isl_dim_param
, i
);
52 set
= isl_set_fix_val(set
, isl_dim_param
, i
, v
);
55 v
= isl_set_count_val(set
);
57 n
= isl_pw_qpolynomial_eval(isl_pw_qpolynomial_copy(vpe
->pwqp
),
60 if (pa
== BV_APPROX_SIGN_LOWER
)
62 else if (pa
== BV_APPROX_SIGN_UPPER
)
67 if (pa
== BV_APPROX_SIGN_APPROX
)
68 /* just accept everything */
70 else if (pa
== BV_APPROX_SIGN_LOWER
)
71 ok
= isl_val_le(n
, v
);
72 else if (pa
== BV_APPROX_SIGN_UPPER
)
73 ok
= isl_val_ge(n
, v
);
75 ok
= isl_val_eq(n
, v
);
77 if (vpe
->vpd
.options
->print_all
|| !ok
) {
78 isl_ctx
*ctx
= isl_point_get_ctx(pnt
);
80 p
= isl_printer_to_file(ctx
, out
);
81 p
= isl_printer_print_str(p
, "EP(");
82 for (i
= 0; i
< nparam
; ++i
) {
84 p
= isl_printer_print_str(p
, ", ");
85 t
= isl_point_get_coordinate_val(pnt
, isl_dim_param
, i
);
86 p
= isl_printer_print_val(p
, t
);
89 p
= isl_printer_print_str(p
, ") = ");
90 p
= isl_printer_print_val(p
, n
);
91 p
= isl_printer_print_str(p
, ", count = ");
92 p
= isl_printer_print_val(p
, v
);
94 p
= isl_printer_print_str(p
, ". OK");
96 p
= isl_printer_print_str(p
, ". NOT OK");
97 p
= isl_printer_end_line(p
);
99 } else if ((vpe
->vpd
.n
% vpe
->vpd
.s
) == 0) {
116 if (vpe
->vpd
.options
->continue_on_error
)
119 return (vpe
->vpd
.n
>= 1 && ok
) ? isl_stat_ok
: isl_stat_error
;
122 static int verify_isl(Polyhedron
*P
, Polyhedron
*C
,
123 evalue
*EP
, const struct verify_options
*options
)
125 struct verify_point_enum vpe
= { { options
} };
127 isl_ctx
*ctx
= isl_ctx_alloc();
133 dim
= isl_space_set_alloc(ctx
, C
->Dimension
, P
->Dimension
- C
->Dimension
);
134 for (i
= 0; i
< C
->Dimension
; ++i
)
135 dim
= isl_space_set_dim_name(dim
, isl_dim_param
, i
, options
->params
[i
]);
136 set
= isl_set_new_from_polylib(P
, isl_space_copy(dim
));
137 dim
= isl_space_params(dim
);
138 set_C
= isl_set_new_from_polylib(C
, dim
);
139 set_C
= isl_set_intersect_params(isl_set_copy(set
), set_C
);
140 set_C
= isl_set_params(set_C
);
142 set_C
= verify_context_set_bounds(set_C
, options
);
144 r
= verify_point_data_init(&vpe
.vpd
, set_C
);
147 vpe
.pwqp
= isl_pw_qpolynomial_from_evalue(isl_set_get_space(set_C
), EP
);
149 isl_set_foreach_point(set_C
, verify_point
, &vpe
);
153 isl_pw_qpolynomial_free(vpe
.pwqp
);
159 verify_point_data_fini(&vpe
.vpd
);
164 static int verify(Polyhedron
*P
, Polyhedron
*C
, evalue
*EP
, skewed_gen_fun
*gf
,
165 struct enumerate_options
*options
)
171 if (!options
->series
|| options
->function
)
172 return verify_isl(P
, C
, EP
, options
->verify
);
174 CS
= check_poly_context_scan(P
, &C
, C
->Dimension
, options
->verify
);
176 p
= Vector_Alloc(P
->Dimension
+2);
177 value_set_si(p
->p
[P
->Dimension
+1], 1);
179 /* S = scanning list of polyhedra */
180 S
= Polyhedron_Scan(P
, C
, options
->verify
->barvinok
->MaxRays
);
182 check_poly_init(C
, options
->verify
);
184 /******* CHECK NOW *********/
186 if (!check_poly_gf(S
, CS
, gf
, 0, C
->Dimension
, 0, p
->p
,
193 fprintf(stderr
,"Check failed !\n");
195 if (!options
->verify
->print_all
)
207 /* frees M and Minv */
208 static void apply_transformation(Polyhedron
**P
, Polyhedron
**C
,
209 bool free_P
, bool free_C
,
210 Matrix
*M
, Matrix
*Minv
, Matrix
**inv
,
211 barvinok_options
*options
)
216 M2
= align_matrix(M
, (*P
)->Dimension
+ 1);
218 *P
= Polyhedron_Preimage(*P
, M2
, options
->MaxRays
);
224 *C
= Polyhedron_Preimage(*C
, M
, options
->MaxRays
);
232 *inv
= Matrix_Alloc(Minv
->NbRows
, T
->NbColumns
);
233 Matrix_Product(Minv
, T
, *inv
);
240 /* Since we have "compressed" the parameters (in case there were
241 * any equalities), the result is independent of the coordinates in the
242 * coordinate subspace spanned by the lines. We can therefore assume
243 * these coordinates are zero and compute the inverse image of the map
244 * from a lower dimensional space that adds zeros in the appropriate
247 static void remove_lines(Polyhedron
*C
, Matrix
**M
, Matrix
**Minv
)
249 Matrix
*L
= Matrix_Alloc(C
->Dimension
+1, C
->Dimension
+1);
250 for (int r
= 0; r
< C
->NbBid
; ++r
)
251 Vector_Copy(C
->Ray
[r
]+1, L
->p
[r
], C
->Dimension
);
252 unimodular_complete(L
, C
->NbBid
);
253 assert(value_one_p(L
->p
[C
->Dimension
][C
->Dimension
]));
254 assert(First_Non_Zero(L
->p
[C
->Dimension
], C
->Dimension
) == -1);
255 Matrix_Transposition(L
);
256 assert(First_Non_Zero(L
->p
[C
->Dimension
], C
->Dimension
) == -1);
258 *M
= Matrix_Alloc(C
->Dimension
+1, C
->Dimension
-C
->NbBid
+1);
259 for (int i
= 0; i
< C
->Dimension
+1; ++i
)
260 Vector_Copy(L
->p
[i
]+C
->NbBid
, (*M
)->p
[i
], C
->Dimension
-C
->NbBid
+1);
262 Matrix
*Linv
= Matrix_Alloc(C
->Dimension
+1, C
->Dimension
+1);
263 int ok
= Matrix_Inverse(L
, Linv
);
267 *Minv
= Matrix_Alloc(C
->Dimension
-C
->NbBid
+1, C
->Dimension
+1);
268 for (int i
= C
->NbBid
; i
< C
->Dimension
+1; ++i
)
269 Vector_AntiScale(Linv
->p
[i
], (*Minv
)->p
[i
-C
->NbBid
],
270 Linv
->p
[C
->Dimension
][C
->Dimension
], C
->Dimension
+1);
274 static skewed_gen_fun
*series(Polyhedron
*P
, Polyhedron
* C
,
275 barvinok_options
*options
)
284 /* Compute true context */
285 C1
= Polyhedron_Project(P
, C
->Dimension
);
286 C2
= DomainIntersection(C
, C1
, options
->MaxRays
);
289 POL_ENSURE_VERTICES(C2
);
290 if (C2
->NbBid
!= 0) {
292 if (C2
->NbEq
|| P
->NbEq
) {
293 /* We remove all equalities to be sure all lines are unit vectors */
295 remove_all_equalities(&PT
, &CT
, &CP
, NULL
, C2
->Dimension
,
302 inv
= left_inverse(CP
, &eq
);
308 div
= Matrix_Alloc(inv
->NbRows
-1, inv
->NbColumns
+1);
309 for (int i
= 0; i
< inv
->NbRows
-1; ++i
) {
310 Vector_Gcd(inv
->p
[i
], inv
->NbColumns
, &tmp
);
311 if (mpz_divisible_p(tmp
,
312 inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]))
314 Vector_Copy(inv
->p
[i
], div
->p
[d
], inv
->NbColumns
);
315 value_assign(div
->p
[d
][inv
->NbColumns
],
316 inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]);
328 POL_ENSURE_VERTICES(C2
);
332 remove_lines(C2
, &M
, &Minv
);
333 apply_transformation(&PT
, &C2
, PT
!= P
, C2
!= C
, M
, Minv
, &inv
,
337 POL_ENSURE_VERTICES(C2
);
338 if (!Polyhedron_has_revlex_positive_rays(C2
, C2
->Dimension
)) {
341 Constraints
= Matrix_Alloc(C2
->NbConstraints
, C2
->Dimension
+1);
342 for (int i
= 0; i
< C2
->NbConstraints
; ++i
)
343 Vector_Copy(C2
->Constraint
[i
]+1, Constraints
->p
[i
], C2
->Dimension
);
344 left_hermite(Constraints
, &H
, &Q
, &U
);
345 Matrix_Free(Constraints
);
347 for (int i
= 0; i
< C2
->Dimension
/2; ++i
)
348 Vector_Exchange(Q
->p
[i
], Q
->p
[C2
->Dimension
-1-i
], C2
->Dimension
);
351 Matrix
*M
= Matrix_Alloc(C2
->Dimension
+1, C2
->Dimension
+1);
353 int ok
= Matrix_Inverse(U
, M
);
357 apply_transformation(&PT
, &C2
, PT
!= P
, C2
!= C
, M
, Q
, &inv
, options
);
359 gf
= barvinok_series_with_options(PT
, C2
, options
);
363 return new skewed_gen_fun(gf
, inv
, eq
, div
);
366 int main(int argc
, char **argv
)
371 skewed_gen_fun
*gf
= NULL
;
372 const char **param_name
;
373 int print_solution
= 1;
375 struct enumerate_options
*options
= enumerate_options_new_with_defaults();
377 argc
= enumerate_options_parse(options
, argc
, argv
, ISL_ARG_ALL
);
381 A
= Constraints2Polyhedron(M
, options
->verify
->barvinok
->MaxRays
);
385 C
= Constraints2Polyhedron(M
, options
->verify
->barvinok
->MaxRays
);
387 assert(A
->Dimension
>= C
->Dimension
);
388 param_name
= Read_ParamNames(stdin
, C
->Dimension
);
390 if (options
->verify
->verify
) {
391 verify_options_set_range(options
->verify
, A
->Dimension
);
392 if (!options
->verify
->barvinok
->verbose
)
396 if (print_solution
&& options
->verify
->barvinok
->verbose
) {
397 Polyhedron_Print(stdout
, P_VALUE_FMT
, A
);
398 Polyhedron_Print(stdout
, P_VALUE_FMT
, C
);
401 if (options
->series
) {
402 gf
= series(A
, C
, options
->verify
->barvinok
);
403 if (print_solution
) {
404 gf
->print(cout
, C
->Dimension
, param_name
);
407 if (options
->function
) {
410 print_evalue(stdout
, EP
, param_name
);
413 EP
= barvinok_enumerate_with_options(A
, C
, options
->verify
->barvinok
);
415 if (evalue_convert(EP
, options
->convert
, options
->verify
->barvinok
->verbose
,
416 C
->Dimension
, param_name
))
419 printf("\nSize: %zd\n", evalue_size(EP
));
421 print_evalue(stdout
, EP
, param_name
);
424 if (options
->verify
->verify
) {
425 options
->verify
->params
= param_name
;
426 result
= verify(A
, C
, EP
, gf
, options
);
434 if (options
->verify
->barvinok
->print_stats
)
435 barvinok_stats_print(options
->verify
->barvinok
->stats
, stdout
);
437 Free_ParamNames(param_name
, C
->Dimension
);
440 enumerate_options_free(options
);