3 #include <barvinok/evalue.h>
4 #include <barvinok/util.h>
5 #include <barvinok/barvinok.h>
9 #include "verif_ehrhart.h"
10 #include "remove_equalities.h"
12 #undef CS /* for Solaris 10 */
14 /* The input of this example program is the same as that of testehrhart
15 * in the PolyLib distribution, i.e., a polytope in combined
16 * data and parameter space, a context polytope in parameter space
17 * and (optionally) the names of the parameters.
18 * Both polytopes are in PolyLib notation.
21 #define PRINT_STATS (BV_OPT_LAST+1)
23 struct argp_option argp_options
[] = {
24 { "convert", 'c', 0, 0, "convert fractionals to periodics" },
25 { "floor", 'f', 0, 0, "convert fractionals to floorings" },
27 { "series", 's', 0, 0, "compute rational generating function" },
28 { "explicit", 'e', 0, 0, "convert rgf to psp" },
30 { "print-stats", PRINT_STATS
, 0, 0 },
42 struct verify_options verify
;
45 static error_t
parse_opt(int key
, char *arg
, struct argp_state
*state
)
47 struct arguments
*options
= (struct arguments
*) state
->input
;
51 state
->child_inputs
[0] = options
->verify
.barvinok
;
52 state
->child_inputs
[1] = &options
->verify
;
57 options
->function
= 0;
59 options
->print_stats
= 0;
62 options
->print_stats
= 1;
74 options
->function
= 1;
83 return ARGP_ERR_UNKNOWN
;
88 struct skewed_gen_fun
{
90 /* maps original space to space in which gf is defined */
92 /* equalities in the original space that need to be satisfied for
96 /* divisibilities in the original space that need to be satisfied for
101 skewed_gen_fun(gen_fun
*gf
, Matrix
*T
, Matrix
*eq
, Matrix
*div
) :
102 gf(gf
), T(T
), eq(eq
), div(div
) {}
113 void print(FILE *out
, unsigned int nparam
, char **param_name
) const;
114 operator evalue
*() const {
115 assert(T
== NULL
&& eq
== NULL
); /* other cases not supported for now */
118 void coefficient(Value
* params
, Value
* c
, barvinok_options
*options
) const;
121 void skewed_gen_fun::print(FILE *out
, unsigned int nparam
,
122 char **param_name
) const
124 fdostream
os(dup(fileno(out
)));
126 fprintf(out
, "T:\n");
127 Matrix_Print(out
, P_VALUE_FMT
, T
);
130 fprintf(out
, "eq:\n");
131 Matrix_Print(out
, P_VALUE_FMT
, eq
);
134 fprintf(out
, "div:\n");
135 Matrix_Print(out
, P_VALUE_FMT
, div
);
137 gf
->print(os
, nparam
, param_name
);
140 void skewed_gen_fun::coefficient(Value
* params
, Value
* c
,
141 barvinok_options
*options
) const
144 for (int i
= 0; i
< eq
->NbRows
; ++i
) {
145 Inner_Product(eq
->p
[i
]+1, params
, eq
->NbColumns
-2, eq
->p
[i
]);
146 if (value_notzero_p(eq
->p
[i
][0])) {
155 for (int i
= 0; i
< div
->NbRows
; ++i
) {
156 Inner_Product(div
->p
[i
], params
, div
->NbColumns
-1, &tmp
);
157 if (!mpz_divisible_p(tmp
, div
->p
[i
][div
->NbColumns
-1])) {
167 coeff
= gf
->coefficient(params
, options
);
169 Vector
*p2
= Vector_Alloc(T
->NbRows
);
170 Matrix_Vector_Product(T
, params
, p2
->p
);
171 if (value_notone_p(p2
->p
[T
->NbRows
-1]))
172 Vector_AntiScale(p2
->p
, p2
->p
, p2
->p
[T
->NbRows
-1], T
->NbRows
);
173 coeff
= gf
->coefficient(p2
->p
, options
);
180 static int check_series(Polyhedron
*S
, Polyhedron
*CS
, skewed_gen_fun
*gf
,
181 int nparam
, int pos
, Value
*z
, verify_options
*options
)
193 /* Computes the coefficient */
194 gf
->coefficient(&z
[S
->Dimension
-nparam
+1], &c
, options
->barvinok
);
196 /* if c=0 we may be out of context. */
197 /* scanning is useless in this case*/
199 if (options
->print_all
) {
201 value_print(stdout
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
202 for(k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
204 value_print(stdout
,VALUE_FMT
,z
[k
]);
207 value_print(stdout
,VALUE_FMT
,c
);
211 /* Manually count the number of points */
212 count_points(1,S
,z
,&tmp
);
213 if (options
->print_all
) {
214 printf(", count = ");
215 value_print(stdout
, P_VALUE_FMT
, tmp
);
219 if (value_ne(tmp
,c
)) {
222 fprintf(stderr
,"Error !\n");
223 fprintf(stderr
,"EP( ");
224 value_print(stderr
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
225 for (k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
226 fprintf(stderr
,", ");
227 value_print(stderr
,VALUE_FMT
,z
[k
]);
229 fprintf(stderr
," ) should be ");
230 value_print(stderr
,VALUE_FMT
,tmp
);
231 fprintf(stderr
,", while EP eval gives ");
232 value_print(stderr
,VALUE_FMT
,c
);
233 fprintf(stderr
,".\n");
234 if (!options
->continue_on_error
) {
235 value_clear(c
); value_clear(tmp
);
238 } else if (options
->print_all
)
242 !(lower_upper_bounds(1+pos
, CS
, &z
[S
->Dimension
-nparam
], &LB
, &UB
));
244 for (value_assign(tmp
,LB
); value_le(tmp
,UB
); value_increment(tmp
,tmp
)) {
245 if (!options
->print_all
) {
246 k
= VALUE_TO_INT(tmp
);
247 if(!pos
&& !(k
% options
->st
)) {
252 value_assign(z
[pos
+S
->Dimension
-nparam
+1],tmp
);
253 if (!check_series(S
, CS
->next
, gf
, nparam
, pos
+1, z
, options
)) {
254 value_clear(c
); value_clear(tmp
);
260 value_set_si(z
[pos
+S
->Dimension
-nparam
+1],0);
270 static int verify(Polyhedron
*P
, Polyhedron
**C
, evalue
*EP
, skewed_gen_fun
*gf
,
273 Polyhedron
*CC
, *PP
, *CS
, *S
;
278 /******* Compute true context *******/
279 CC
= align_context(*C
, P
->Dimension
, options
->verify
.barvinok
->MaxRays
);
280 PP
= DomainIntersection(P
, CC
, options
->verify
.barvinok
->MaxRays
);
282 C1
= Matrix_Alloc((*C
)->Dimension
+1, P
->Dimension
+1);
284 for (int i
= 0; i
< C1
->NbRows
; i
++)
285 for (int j
= 0; j
< C1
->NbColumns
; j
++)
286 if (i
== j
-P
->Dimension
+(*C
)->Dimension
)
287 value_set_si(C1
->p
[i
][j
], 1);
289 value_set_si(C1
->p
[i
][j
], 0);
290 CC
= Polyhedron_Image(PP
, C1
, options
->verify
.barvinok
->MaxRays
);
296 CS
= check_poly_context_scan(*C
, &options
->verify
);
298 p
= Vector_Alloc(P
->Dimension
+2);
299 value_set_si(p
->p
[P
->Dimension
+1], 1);
301 /* S = scanning list of polyhedra */
302 S
= Polyhedron_Scan(P
, *C
, options
->verify
.barvinok
->MaxRays
);
304 check_poly_init(*C
, &options
->verify
);
306 /******* CHECK NOW *********/
308 if (!options
->series
|| options
->function
) {
309 if (!check_poly(S
, CS
, EP
, 0, (*C
)->Dimension
, 0, p
->p
,
313 if (!check_series(S
, CS
, gf
, (*C
)->Dimension
, 0, p
->p
, &options
->verify
))
320 fprintf(stderr
,"Check failed !\n");
322 if (!options
->verify
.print_all
)
332 static void unimodular_complete(Matrix
*M
, int row
)
335 left_hermite(M
, &H
, &Q
, &U
);
338 for (int r
= row
; r
< M
->NbRows
; ++r
)
339 Vector_Copy(Q
->p
[r
], M
->p
[r
], M
->NbColumns
);
343 /* frees M and Minv */
344 static void apply_transformation(Polyhedron
**P
, Polyhedron
**C
,
345 bool free_P
, bool free_C
,
346 Matrix
*M
, Matrix
*Minv
, Matrix
**inv
,
347 barvinok_options
*options
)
352 M2
= align_matrix(M
, (*P
)->Dimension
+ 1);
354 *P
= Polyhedron_Preimage(*P
, M2
, options
->MaxRays
);
360 *C
= Polyhedron_Preimage(*C
, M
, options
->MaxRays
);
368 *inv
= Matrix_Alloc(Minv
->NbRows
, T
->NbColumns
);
369 Matrix_Product(Minv
, T
, *inv
);
376 static skewed_gen_fun
*series(Polyhedron
*P
, Polyhedron
* C
,
377 barvinok_options
*options
)
386 /* Compute true context */
387 C1
= Polyhedron_Project(P
, C
->Dimension
);
388 C2
= DomainIntersection(C
, C1
, options
->MaxRays
);
391 POL_ENSURE_VERTICES(C2
);
392 if (C2
->NbBid
!= 0) {
394 Matrix
*M
, *Minv
, *M2
;
396 if (C2
->NbEq
|| P
->NbEq
) {
397 /* We remove all equalities to be sure all lines are unit vectors */
399 remove_all_equalities(&PT
, &CT
, &CP
, NULL
, C2
->Dimension
,
406 inv
= left_inverse(CP
, &eq
);
412 div
= Matrix_Alloc(inv
->NbRows
-1, inv
->NbColumns
+1);
413 for (int i
= 0; i
< inv
->NbRows
-1; ++i
) {
414 Vector_Gcd(inv
->p
[i
], inv
->NbColumns
, &tmp
);
415 if (mpz_divisible_p(tmp
,
416 inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]))
418 Vector_Copy(inv
->p
[i
], div
->p
[d
], inv
->NbColumns
);
419 value_assign(div
->p
[d
][inv
->NbColumns
],
420 inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]);
432 POL_ENSURE_VERTICES(C2
);
434 /* Since we have "compressed" the parameters (in case there were
435 * any equalities), the result is independent of the coordinates in the
436 * coordinate subspace spanned by the lines. We can therefore assume
437 * these coordinates are zero and compute the inverse image of the map
438 * from a lower dimensional space that adds zeros in the appropriate
441 M
= Matrix_Alloc(C2
->Dimension
+1, C2
->Dimension
-C2
->NbBid
+1);
443 for (int i
= 0; i
< C2
->NbBid
; ++i
) {
444 int j
= First_Non_Zero(C2
->Ray
[i
]+1, C2
->Dimension
);
445 assert(First_Non_Zero(C2
->Ray
[i
]+1+j
+1, C2
->Dimension
-j
-1) == -1);
447 value_set_si(M
->p
[k
+i
][k
], 1);
449 for ( ; k
< C2
->Dimension
-C2
->NbBid
+1; k
++)
450 value_set_si(M
->p
[k
+C2
->NbBid
][k
], 1);
453 apply_transformation(&PT
, &C2
, PT
!= P
, C2
!= C
, M
, Minv
, &inv
, options
);
455 POL_ENSURE_VERTICES(C2
);
456 if (!Polyhedron_has_revlex_positive_rays(C2
, C2
->Dimension
)) {
460 Constraints
= Matrix_Alloc(C2
->NbConstraints
, C2
->Dimension
+1);
461 for (int i
= 0; i
< C2
->NbConstraints
; ++i
)
462 Vector_Copy(C2
->Constraint
[i
]+1, Constraints
->p
[i
], C2
->Dimension
);
463 left_hermite(Constraints
, &H
, &Q
, &U
);
465 for (int i
= 0; i
< C2
->Dimension
/2; ++i
)
466 Vector_Exchange(Q
->p
[i
], Q
->p
[C2
->Dimension
-1-i
], C2
->Dimension
);
469 Matrix
*M
= Matrix_Alloc(C2
->Dimension
+1, C2
->Dimension
+1);
471 int ok
= Matrix_Inverse(U
, M
);
475 apply_transformation(&PT
, &C2
, PT
!= P
, C2
!= C
, M
, Q
, &inv
, options
);
477 gf
= barvinok_series_with_options(PT
, C2
, options
);
481 return new skewed_gen_fun(gf
, inv
, eq
, div
);
484 int main(int argc
, char **argv
)
489 skewed_gen_fun
*gf
= NULL
;
491 int print_solution
= 1;
493 struct arguments options
;
494 static struct argp_child argp_children
[] = {
495 { &barvinok_argp
, 0, 0, 0 },
496 { &verify_argp
, 0, "verification", 1 },
499 static struct argp argp
= { argp_options
, parse_opt
, 0, 0, argp_children
};
500 struct barvinok_options
*bv_options
= barvinok_options_new_with_defaults();
502 options
.verify
.barvinok
= bv_options
;
503 argp_parse(&argp
, argc
, argv
, 0, 0, &options
);
506 A
= Constraints2Polyhedron(M
, bv_options
->MaxRays
);
509 C
= Constraints2Polyhedron(M
, bv_options
->MaxRays
);
511 param_name
= Read_ParamNames(stdin
, C
->Dimension
);
513 if (options
.verify
.verify
) {
514 verify_options_set_range(&options
.verify
, A
);
515 if (!options
.verbose
)
519 if (print_solution
) {
520 Polyhedron_Print(stdout
, P_VALUE_FMT
, A
);
521 Polyhedron_Print(stdout
, P_VALUE_FMT
, C
);
524 if (options
.series
) {
525 gf
= series(A
, C
, bv_options
);
526 if (print_solution
) {
527 gf
->print(stdout
, C
->Dimension
, param_name
);
530 if (options
.function
) {
533 print_evalue(stdout
, EP
, param_name
);
536 EP
= barvinok_enumerate_with_options(A
, C
, bv_options
);
538 print_evalue(stdout
, EP
, param_name
);
540 printf("\nSize: %d\n", evalue_size(EP
));
542 fprintf(stderr
, "WARNING: floor conversion not supported\n");
543 evalue_frac2floor2(EP
, 0);
544 print_evalue(stdout
, EP
, param_name
);
545 } else if (options
.convert
) {
546 evalue_mod2table(EP
, C
->Dimension
);
547 print_evalue(stdout
, EP
, param_name
);
549 printf("\nSize: %d\n", evalue_size(EP
));
553 if (options
.verify
.verify
) {
554 options
.verify
.params
= param_name
;
555 result
= verify(A
, &C
, EP
, gf
, &options
);
561 free_evalue_refs(EP
);
565 if (options
.print_stats
)
566 barvinok_stats_print(options
.verify
.barvinok
->stats
, stdout
);
568 Free_ParamNames(param_name
, C
->Dimension
);
571 barvinok_options_free(bv_options
);