4 #include <barvinok/util.h>
5 #include <barvinok/barvinok.h>
9 #include "omega/convert.h"
12 /* The input of this example program is a polytope in combined
13 * data and parameter space followed by two lines indicating
14 * the number of existential variables and parameters respectively.
15 * The first lines starts with "E ", followed by a number.
16 * The second lines starts with "P ", followed by a number.
17 * These two lines are (optionally) followed by the names of the parameters.
18 * The polytope is in PolyLib notation.
21 struct argp_option argp_options
[] = {
23 { "omega", 'o', 0, 0 },
28 { "series", 's', 0, 0 },
29 { "scarf", 'S', 0, 0 },
30 { "convert", 'c', 0, 0 },
31 { "floor", 'f', 0, 0 },
32 { "range-reduction", 'R', 0, 0 },
33 { "verify", 'T', 0, 0 },
34 { "print-all", 'A', 0, 0 },
35 { "min", 'm', "int", 0 },
36 { "max", 'M', "int", 0 },
37 { "range", 'r', "int", 0 },
42 struct barvinok_options
*options
;
56 error_t
parse_opt(int key
, char *arg
, struct argp_state
*state
)
58 struct arguments
*arguments
= (struct arguments
*)(state
->input
);
62 state
->child_inputs
[0] = arguments
->options
;
65 arguments
->series
= 1;
80 arguments
->convert
= 1;
86 arguments
->verify
= 1;
89 arguments
->print_all
= 1;
92 arguments
->m
= atoi(arg
);
93 arguments
->verify
= 1;
96 arguments
->M
= atoi(arg
);
97 arguments
->verify
= 1;
100 arguments
->M
= atoi(arg
);
101 arguments
->m
= -arguments
->M
;
102 arguments
->verify
= 1;
105 return ARGP_ERR_UNKNOWN
;
112 Polyhedron
*Omega_simplify(Polyhedron
*P
,
113 unsigned exist
, unsigned nparam
, char **parms
)
117 Relation r
= Polyhedron2relation(P
, exist
, nparam
, parms
);
119 return relation2Domain(r
, varv
, paramv
);
122 Polyhedron
*Omega_simplify(Polyhedron
*P
,
123 unsigned exist
, unsigned nparam
, char **parms
)
129 /* define this to continue the test after first error found */
130 /* #define DONT_BREAK_ON_ERROR */
132 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
135 /* SRANGE : small range for evalutations */
138 /* if dimension >= BIDDIM, use SRANGE */
141 /* VSRANGE : very small range for evalutations */
144 /* if dimension >= VBIDDIM, use VSRANGE */
147 static Value min_val
, max_val
;
149 static char **params
;
153 static int check_poly(Polyhedron
*S
, Polyhedron
*C
, evalue
*EP
,
154 int exist
, int nparam
, int pos
, Value
*z
, int print_all
);
155 static void verify_results(Polyhedron
*P
, evalue
*EP
, int exist
, int nparam
,
156 int m
, int M
, int print_all
, unsigned MaxRays
);
158 int main(int argc
, char **argv
)
163 int exist
, nparam
, nvar
;
166 int print_solution
= 1;
168 struct arguments arguments
;
169 static struct argp_child argp_children
[] = {
170 { &barvinok_argp
, 0, 0, 0 },
173 static struct argp argp
= { argp_options
, parse_opt
, 0, 0, argp_children
};
174 struct barvinok_options
*options
= barvinok_options_new_with_defaults();
176 arguments
.options
= options
;
178 arguments
.convert
= 0;
182 arguments
.series
= 0;
184 arguments
.verify
= 0;
185 arguments
.print_all
= 0;
186 arguments
.m
= INT_MAX
;
187 arguments
.M
= INT_MIN
;
189 argp_parse(&argp
, argc
, argv
, ARGP_NO_EXIT
, 0, &arguments
);
191 if (arguments
.series
&& !arguments
.scarf
) {
193 "--series currently only available if --scarf is specified\n");
198 A
= Constraints2Polyhedron(MA
, options
->MaxRays
);
201 fgets(s
, 128, stdin
);
202 while ((*s
=='#') || (sscanf(s
, "E %d", &exist
)<1))
203 fgets(s
, 128, stdin
);
205 fgets(s
, 128, stdin
);
206 while ((*s
=='#') || (sscanf(s
, "P %d", &nparam
)<1))
207 fgets(s
, 128, stdin
);
209 /******* Read the options: initialize Min and Max ********/
210 if (A
->Dimension
>= VBIGDIM
)
212 else if (A
->Dimension
>= BIGDIM
)
216 if (arguments
.M
== INT_MIN
)
218 if (arguments
.m
== INT_MAX
)
221 if (arguments
.verify
&& arguments
.m
> arguments
.M
) {
222 fprintf(stderr
,"Nothing to do: min > max !\n");
225 if (arguments
.verify
)
228 if (print_solution
) {
229 Polyhedron_Print(stdout
, P_VALUE_FMT
, A
);
230 printf("exist: %d, nparam: %d\n", exist
, nparam
);
232 param_name
= Read_ParamNames(stdin
, nparam
);
233 nvar
= A
->Dimension
- exist
- nparam
;
234 if (arguments
.omega
) {
235 A
= Omega_simplify(A
, exist
, nparam
, param_name
);
237 exist
= A
->Dimension
- nvar
- nparam
;
239 if (arguments
.series
) {
241 barvinok_options
*options
= barvinok_options_new_with_defaults();
242 assert(arguments
.scarf
);
243 gf
= barvinok_enumerate_scarf_series(A
, exist
, nparam
, options
);
244 if (print_solution
) {
245 gf
->print(std::cout
, nparam
, param_name
);
251 if (arguments
.scarf
) {
252 barvinok_options
*options
= barvinok_options_new_with_defaults();
253 EP
= barvinok_enumerate_scarf(A
, exist
, nparam
, options
);
255 } else if (arguments
.pip
&& exist
> 0)
256 EP
= barvinok_enumerate_pip(A
, exist
, nparam
, options
->MaxRays
);
258 EP
= barvinok_enumerate_e(A
, exist
, nparam
, options
->MaxRays
);
262 evalue_range_reduction(EP
);
264 print_evalue(stdout
, EP
, param_name
);
265 if (arguments
.floor
) {
266 fprintf(stderr
, "WARNING: floor conversion not supported\n");
267 evalue_frac2floor2(EP
, 0);
269 print_evalue(stdout
, EP
, param_name
);
270 } else if (arguments
.convert
) {
271 evalue_mod2table(EP
, nparam
);
273 print_evalue(stdout
, EP
, param_name
);
275 if (arguments
.verify
)
276 verify_results(A
, EP
, exist
, nparam
, arguments
.m
, arguments
.M
,
277 arguments
.print_all
, options
->MaxRays
);
278 free_evalue_refs(EP
);
281 Free_ParamNames(param_name
, nparam
);
286 void verify_results(Polyhedron
*P
, evalue
*EP
, int exist
, int nparam
, int m
, int M
,
287 int print_all
, unsigned MaxRays
)
293 Polyhedron
*C
= Polyhedron_Project(P
, nparam
);
296 value_set_si(min_val
,m
);
297 value_set_si(max_val
,M
);
300 p
= (Value
*)malloc(sizeof(Value
) * (P
->Dimension
+2));
301 for(i
=0;i
<=P
->Dimension
;i
++) {
303 value_set_si(p
[i
],0);
306 value_set_si(p
[i
],1);
308 /* S = scanning list of polyhedra */
309 S
= Polyhedron_Scan(P
, C
, MaxRays
& POL_NO_DUAL
? 0 : MaxRays
);
312 if (C
->Dimension
> 0) {
313 value_subtract(tmp
,max_val
,min_val
);
314 if (VALUE_TO_INT(tmp
) > 80)
315 st
= 1+(VALUE_TO_INT(tmp
))/80;
318 for(i
=VALUE_TO_INT(min_val
);i
<=VALUE_TO_INT(max_val
);i
+=st
)
325 /******* CHECK NOW *********/
327 if(S
&& !check_poly(S
, C
, EP
, exist
, nparam
, 0, p
, print_all
)) {
328 fprintf(stderr
,"Check failed !\n");
335 for(i
=0;i
<=(P
->Dimension
+1);i
++)
343 /****************************************************/
344 /* function check_poly : */
345 /* scans the parameter space from min to max (all */
346 /* directions). Computes the number of points in */
347 /* the polytope using both methods, and compare them*/
348 /* returns 1 on success */
349 /****************************************************/
351 int check_poly(Polyhedron
*S
, Polyhedron
*C
, evalue
*EP
,
352 int exist
, int nparam
, int pos
, Value
*z
, int print_all
)
357 value_init(c
); value_init(tmp
);
361 /* Computes the ehrhart polynomial */
362 value_set_double(c
, compute_evalue(EP
,&z
[S
->Dimension
-nparam
+1])+.25);
363 /* if c=0 we may be out of context. */
364 /* scanning is useless in this case*/
365 if(!in_domain(C
,&z
[S
->Dimension
-nparam
+1])) {
373 value_print(stdout
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
374 for(k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
376 value_print(stdout
,VALUE_FMT
,z
[k
]);
379 value_print(stdout
,VALUE_FMT
,c
);
383 /* Manually count the number of points */
384 count_points_e(1, S
, exist
, nparam
, z
, &tmp
);
386 printf(", count = ");
387 value_print(stdout
, P_VALUE_FMT
, tmp
);
391 if(value_ne(tmp
,c
)) {
394 fprintf(stderr
,"Error !\n");
395 fprintf(stderr
,"EP( ");
396 value_print(stderr
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
397 for(k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
398 fprintf(stderr
,", ");
399 value_print(stderr
,VALUE_FMT
,z
[k
]);
401 fprintf(stderr
," ) should be ");
402 value_print(stderr
,VALUE_FMT
,tmp
);
403 fprintf(stderr
,", while EP eval gives ");
404 value_print(stderr
,VALUE_FMT
,c
);
405 fprintf(stderr
,".\n");
406 print_evalue(stderr
, EP
, params
);
407 #ifndef DONT_BREAK_ON_ERROR
408 value_clear(c
); value_clear(tmp
);
417 for(value_assign(tmp
,min_val
); value_le(tmp
,max_val
); value_increment(tmp
,tmp
)) {
419 k
= VALUE_TO_INT(tmp
);
420 if(!pos
&& !(k
%st
)) {
426 value_assign(z
[pos
+S
->Dimension
-nparam
+1],tmp
);
427 if(!check_poly(S
, C
, EP
, exist
, nparam
, pos
+1, z
, print_all
)) {
428 value_clear(c
); value_clear(tmp
);
432 value_clear(c
); value_clear(tmp
);