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 },
35 struct barvinok_options
*barvinok
;
43 struct verify_options verify
;
46 static error_t
parse_opt(int key
, char *arg
, struct argp_state
*state
)
48 struct arguments
*options
= (struct arguments
*) state
->input
;
52 state
->child_inputs
[0] = options
->barvinok
;
53 state
->child_inputs
[1] = &options
->verify
;
58 options
->function
= 0;
60 options
->print_stats
= 0;
63 options
->print_stats
= 1;
75 options
->function
= 1;
84 return ARGP_ERR_UNKNOWN
;
89 struct skewed_gen_fun
{
91 /* maps original space to space in which gf is defined */
93 /* equalities in the original space that need to be satisfied for
97 /* divisibilities in the original space that need to be satisfied for
102 skewed_gen_fun(gen_fun
*gf
, Matrix
*T
, Matrix
*eq
, Matrix
*div
) :
103 gf(gf
), T(T
), eq(eq
), div(div
) {}
114 void print(FILE *out
, unsigned int nparam
, char **param_name
) const;
115 operator evalue
*() const {
116 assert(T
== NULL
&& eq
== NULL
); /* other cases not supported for now */
119 void coefficient(Value
* params
, Value
* c
, barvinok_options
*options
) const;
122 void skewed_gen_fun::print(FILE *out
, unsigned int nparam
,
123 char **param_name
) const
125 fdostream
os(dup(fileno(out
)));
127 fprintf(out
, "T:\n");
128 Matrix_Print(out
, P_VALUE_FMT
, T
);
131 fprintf(out
, "eq:\n");
132 Matrix_Print(out
, P_VALUE_FMT
, eq
);
135 fprintf(out
, "div:\n");
136 Matrix_Print(out
, P_VALUE_FMT
, div
);
138 gf
->print(os
, nparam
, param_name
);
141 void skewed_gen_fun::coefficient(Value
* params
, Value
* c
,
142 barvinok_options
*options
) const
145 for (int i
= 0; i
< eq
->NbRows
; ++i
) {
146 Inner_Product(eq
->p
[i
]+1, params
, eq
->NbColumns
-2, eq
->p
[i
]);
147 if (value_notzero_p(eq
->p
[i
][0])) {
156 for (int i
= 0; i
< div
->NbRows
; ++i
) {
157 Inner_Product(div
->p
[i
], params
, div
->NbColumns
-1, &tmp
);
158 if (!mpz_divisible_p(tmp
, div
->p
[i
][div
->NbColumns
-1])) {
168 coeff
= gf
->coefficient(params
, options
);
170 Vector
*p2
= Vector_Alloc(T
->NbRows
);
171 Matrix_Vector_Product(T
, params
, p2
->p
);
172 if (value_notone_p(p2
->p
[T
->NbRows
-1]))
173 Vector_AntiScale(p2
->p
, p2
->p
, p2
->p
[T
->NbRows
-1], T
->NbRows
);
174 coeff
= gf
->coefficient(p2
->p
, options
);
181 static int check_series(Polyhedron
*S
, Polyhedron
*CS
, skewed_gen_fun
*gf
,
182 int nparam
, int pos
, Value
*z
, int print_all
,
183 barvinok_options
*options
)
195 /* Computes the coefficient */
196 gf
->coefficient(&z
[S
->Dimension
-nparam
+1], &c
, options
);
198 /* if c=0 we may be out of context. */
199 /* scanning is useless in this case*/
203 value_print(stdout
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
204 for(k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
206 value_print(stdout
,VALUE_FMT
,z
[k
]);
209 value_print(stdout
,VALUE_FMT
,c
);
213 /* Manually count the number of points */
214 count_points(1,S
,z
,&tmp
);
216 printf(", count = ");
217 value_print(stdout
, P_VALUE_FMT
, tmp
);
221 if (value_ne(tmp
,c
)) {
224 fprintf(stderr
,"Error !\n");
225 fprintf(stderr
,"EP( ");
226 value_print(stderr
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
227 for (k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
228 fprintf(stderr
,", ");
229 value_print(stderr
,VALUE_FMT
,z
[k
]);
231 fprintf(stderr
," ) should be ");
232 value_print(stderr
,VALUE_FMT
,tmp
);
233 fprintf(stderr
,", while EP eval gives ");
234 value_print(stderr
,VALUE_FMT
,c
);
235 fprintf(stderr
,".\n");
236 #ifndef DONT_BREAK_ON_ERROR
237 value_clear(c
); value_clear(tmp
);
240 } else if (print_all
)
244 !(lower_upper_bounds(1+pos
, CS
, &z
[S
->Dimension
-nparam
], &LB
, &UB
));
246 for (value_assign(tmp
,LB
); value_le(tmp
,UB
); value_increment(tmp
,tmp
)) {
248 k
= VALUE_TO_INT(tmp
);
249 if(!pos
&& !(k
%st
)) {
254 value_assign(z
[pos
+S
->Dimension
-nparam
+1],tmp
);
255 if (!check_series(S
, CS
->next
, gf
, nparam
, pos
+1, z
, print_all
,
257 value_clear(c
); value_clear(tmp
);
263 value_set_si(z
[pos
+S
->Dimension
-nparam
+1],0);
273 static int verify(Polyhedron
*P
, Polyhedron
**C
, Enumeration
*en
, skewed_gen_fun
*gf
,
276 Polyhedron
*CC
, *PP
, *CS
, *S
, *U
;
281 /******* Compute true context *******/
282 CC
= align_context(*C
, P
->Dimension
, options
->barvinok
->MaxRays
);
283 PP
= DomainIntersection(P
, CC
, options
->barvinok
->MaxRays
);
285 C1
= Matrix_Alloc((*C
)->Dimension
+1, P
->Dimension
+1);
287 for (int i
= 0; i
< C1
->NbRows
; i
++)
288 for (int j
= 0; j
< C1
->NbColumns
; j
++)
289 if (i
== j
-P
->Dimension
+(*C
)->Dimension
)
290 value_set_si(C1
->p
[i
][j
], 1);
292 value_set_si(C1
->p
[i
][j
], 0);
293 CC
= Polyhedron_Image(PP
, C1
, options
->barvinok
->MaxRays
);
299 /* Intersect context with range */
300 if ((*C
)->Dimension
> 0) {
301 MM
= Matrix_Alloc(2*(*C
)->Dimension
, (*C
)->Dimension
+2);
302 for (int i
= 0; i
< (*C
)->Dimension
; ++i
) {
303 value_set_si(MM
->p
[2*i
][0], 1);
304 value_set_si(MM
->p
[2*i
][1+i
], 1);
305 value_set_si(MM
->p
[2*i
][1+(*C
)->Dimension
], -options
->verify
.m
);
306 value_set_si(MM
->p
[2*i
+1][0], 1);
307 value_set_si(MM
->p
[2*i
+1][1+i
], -1);
308 value_set_si(MM
->p
[2*i
+1][1+(*C
)->Dimension
], options
->verify
.M
);
310 CC
= AddConstraints(MM
->p
[0], 2*(*C
)->Dimension
, *C
,
311 options
->barvinok
->MaxRays
);
312 U
= Universe_Polyhedron(0);
313 CS
= Polyhedron_Scan(CC
, U
, options
->barvinok
->MaxRays
);
320 p
= Vector_Alloc(P
->Dimension
+2);
321 value_set_si(p
->p
[P
->Dimension
+1], 1);
323 /* S = scanning list of polyhedra */
324 S
= Polyhedron_Scan(P
, *C
, options
->barvinok
->MaxRays
);
326 if (!options
->verify
.print_all
)
327 if ((*C
)->Dimension
> 0) {
328 int d
= options
->verify
.M
- options
->verify
.m
;
333 for (int i
= options
->verify
.m
; i
<= options
->verify
.M
; i
+= st
)
339 /******* CHECK NOW *********/
341 if (!options
->series
|| options
->function
) {
342 if (!check_poly(S
, CS
, en
, (*C
)->Dimension
, 0, p
->p
,
343 options
->verify
.print_all
))
346 if (!check_series(S
, CS
, gf
, (*C
)->Dimension
, 0, p
->p
,
347 options
->verify
.print_all
, options
->barvinok
))
354 fprintf(stderr
,"Check failed !\n");
356 if (!options
->verify
.print_all
)
366 static void unimodular_complete(Matrix
*M
, int row
)
369 left_hermite(M
, &H
, &Q
, &U
);
372 for (int r
= row
; r
< M
->NbRows
; ++r
)
373 Vector_Copy(Q
->p
[r
], M
->p
[r
], M
->NbColumns
);
377 /* frees M and Minv */
378 static void apply_transformation(Polyhedron
**P
, Polyhedron
**C
,
379 bool free_P
, bool free_C
,
380 Matrix
*M
, Matrix
*Minv
, Matrix
**inv
,
381 barvinok_options
*options
)
386 M2
= align_matrix(M
, (*P
)->Dimension
+ 1);
388 *P
= Polyhedron_Preimage(*P
, M2
, options
->MaxRays
);
394 *C
= Polyhedron_Preimage(*C
, M
, options
->MaxRays
);
402 *inv
= Matrix_Alloc(Minv
->NbRows
, T
->NbColumns
);
403 Matrix_Product(Minv
, T
, *inv
);
410 static skewed_gen_fun
*series(Polyhedron
*P
, Polyhedron
* C
,
411 barvinok_options
*options
)
420 /* Compute true context */
421 C1
= Polyhedron_Project(P
, C
->Dimension
);
422 C2
= DomainIntersection(C
, C1
, options
->MaxRays
);
425 POL_ENSURE_VERTICES(C2
);
426 if (C2
->NbBid
!= 0) {
428 Matrix
*M
, *Minv
, *M2
;
430 if (C2
->NbEq
|| P
->NbEq
) {
431 /* We remove all equalities to be sure all lines are unit vectors */
433 remove_all_equalities(&PT
, &CT
, &CP
, NULL
, C2
->Dimension
,
440 inv
= left_inverse(CP
, &eq
);
446 div
= Matrix_Alloc(inv
->NbRows
-1, inv
->NbColumns
+1);
447 for (int i
= 0; i
< inv
->NbRows
-1; ++i
) {
448 Vector_Gcd(inv
->p
[i
], inv
->NbColumns
, &tmp
);
449 if (mpz_divisible_p(tmp
,
450 inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]))
452 Vector_Copy(inv
->p
[i
], div
->p
[d
], inv
->NbColumns
);
453 value_assign(div
->p
[d
][inv
->NbColumns
],
454 inv
->p
[inv
->NbRows
-1][inv
->NbColumns
-1]);
466 POL_ENSURE_VERTICES(C2
);
468 /* Since we have "compressed" the parameters (in case there were
469 * any equalities), the result is independent of the coordinates in the
470 * coordinate subspace spanned by the lines. We can therefore assume
471 * these coordinates are zero and compute the inverse image of the map
472 * from a lower dimensional space that adds zeros in the appropriate
475 M
= Matrix_Alloc(C2
->Dimension
+1, C2
->Dimension
-C2
->NbBid
+1);
477 for (int i
= 0; i
< C2
->NbBid
; ++i
) {
478 int j
= First_Non_Zero(C2
->Ray
[i
]+1, C2
->Dimension
);
479 assert(First_Non_Zero(C2
->Ray
[i
]+1+j
+1, C2
->Dimension
-j
-1) == -1);
481 value_set_si(M
->p
[k
+i
][k
], 1);
483 for ( ; k
< C2
->Dimension
-C2
->NbBid
+1; k
++)
484 value_set_si(M
->p
[k
+C2
->NbBid
][k
], 1);
487 apply_transformation(&PT
, &C2
, PT
!= P
, C2
!= C
, M
, Minv
, &inv
, options
);
489 POL_ENSURE_VERTICES(C2
);
490 if (!Polyhedron_has_revlex_positive_rays(C2
, C2
->Dimension
)) {
494 Constraints
= Matrix_Alloc(C2
->NbConstraints
, C2
->Dimension
+1);
495 for (int i
= 0; i
< C2
->NbConstraints
; ++i
)
496 Vector_Copy(C2
->Constraint
[i
]+1, Constraints
->p
[i
], C2
->Dimension
);
497 left_hermite(Constraints
, &H
, &Q
, &U
);
499 for (int i
= 0; i
< C2
->Dimension
/2; ++i
)
500 Vector_Exchange(Q
->p
[i
], Q
->p
[C2
->Dimension
-1-i
], C2
->Dimension
);
503 Matrix
*M
= Matrix_Alloc(C2
->Dimension
+1, C2
->Dimension
+1);
505 int ok
= Matrix_Inverse(U
, M
);
509 apply_transformation(&PT
, &C2
, PT
!= P
, C2
!= C
, M
, Q
, &inv
, options
);
511 gf
= barvinok_series_with_options(PT
, C2
, options
);
515 return new skewed_gen_fun(gf
, inv
, eq
, div
);
518 int main(int argc
, char **argv
)
523 Enumeration
*en
= NULL
;
524 skewed_gen_fun
*gf
= NULL
;
526 int print_solution
= 1;
528 struct arguments options
;
529 static struct argp_child argp_children
[] = {
530 { &barvinok_argp
, 0, 0, 0 },
531 { &verify_argp
, 0, "verification", 1 },
534 static struct argp argp
= { argp_options
, parse_opt
, 0, 0, argp_children
};
535 struct barvinok_options
*bv_options
= barvinok_options_new_with_defaults();
537 options
.barvinok
= bv_options
;
538 argp_parse(&argp
, argc
, argv
, 0, 0, &options
);
541 A
= Constraints2Polyhedron(M
, bv_options
->MaxRays
);
544 C
= Constraints2Polyhedron(M
, bv_options
->MaxRays
);
546 param_name
= Read_ParamNames(stdin
, C
->Dimension
);
548 if (options
.verify
.verify
) {
549 verify_options_set_range(&options
.verify
, A
);
550 if (!options
.verbose
)
554 if (print_solution
) {
555 Polyhedron_Print(stdout
, P_VALUE_FMT
, A
);
556 Polyhedron_Print(stdout
, P_VALUE_FMT
, C
);
559 if (options
.series
) {
560 gf
= series(A
, C
, bv_options
);
561 if (print_solution
) {
562 gf
->print(stdout
, C
->Dimension
, param_name
);
565 if (options
.function
) {
568 print_evalue(stdout
, EP
, param_name
);
571 EP
= barvinok_enumerate_with_options(A
, C
, bv_options
);
573 print_evalue(stdout
, EP
, param_name
);
575 printf("\nSize: %d\n", evalue_size(EP
));
577 fprintf(stderr
, "WARNING: floor conversion not supported\n");
578 evalue_frac2floor2(EP
, 0);
579 print_evalue(stdout
, EP
, param_name
);
580 } else if (options
.convert
) {
581 evalue_mod2table(EP
, C
->Dimension
);
582 print_evalue(stdout
, EP
, param_name
);
584 printf("\nSize: %d\n", evalue_size(EP
));
588 if (options
.verify
.verify
) {
590 en
= partition2enumeration(EP
);
593 result
= verify(A
, &C
, en
, gf
, &options
);
597 Enumeration_Free(en
);
601 free_evalue_refs(EP
);
605 if (options
.print_stats
)
606 barvinok_stats_print(options
.barvinok
->stats
, stdout
);
608 Free_ParamNames(param_name
, C
->Dimension
);
611 barvinok_options_free(bv_options
);