3 #include <barvinok/util.h>
4 #include <barvinok/barvinok.h>
8 #include "omega/convert.h"
11 /* The input of this example program is a polytope in combined
12 * data and parameter space followed by two lines indicating
13 * the number of existential variables and parameters respectively.
14 * The first lines starts with "E ", followed by a number.
15 * The second lines starts with "P ", followed by a number.
16 * These two lines are (optionally) followed by the names of the parameters.
17 * The polytope is in PolyLib notation.
20 #ifdef HAVE_GROWING_CHERNIKOVA
21 #define MAXRAYS POL_NO_DUAL
27 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
30 struct option options
[] = {
32 { "omega", no_argument
, 0, 'o' },
35 { "pip", no_argument
, 0, 'p' },
37 { "series", no_argument
, 0, 's' },
38 { "scarf", no_argument
, 0, 'S' },
39 { "convert", no_argument
, 0, 'c' },
40 { "floor", no_argument
, 0, 'f' },
41 { "range-reduction", no_argument
, 0, 'R' },
42 { "verify", no_argument
, 0, 'T' },
43 { "print-all", no_argument
, 0, 'A' },
44 { "min", required_argument
, 0, 'm' },
45 { "max", required_argument
, 0, 'M' },
46 { "range", required_argument
, 0, 'r' },
47 { "version", no_argument
, 0, 'V' },
53 #define PIPLIB_OPT "p"
61 Polyhedron
*Omega_simplify(Polyhedron
*P
,
62 unsigned exist
, unsigned nparam
, char **parms
)
66 Relation r
= Polyhedron2relation(P
, exist
, nparam
, parms
);
68 return relation2Domain(r
, varv
, paramv
);
72 Polyhedron
*Omega_simplify(Polyhedron
*P
,
73 unsigned exist
, unsigned nparam
, char **parms
)
79 /* define this to continue the test after first error found */
80 /* #define DONT_BREAK_ON_ERROR */
82 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
85 /* SRANGE : small range for evalutations */
88 /* if dimension >= BIDDIM, use SRANGE */
91 /* VSRANGE : very small range for evalutations */
94 /* if dimension >= VBIDDIM, use VSRANGE */
97 static Value min_val
, max_val
;
103 static int check_poly(Polyhedron
*S
, Polyhedron
*C
, evalue
*EP
,
104 int exist
, int nparam
, int pos
, Value
*z
, int print_all
);
105 static void verify_results(Polyhedron
*P
, evalue
*EP
, int exist
, int nparam
,
106 int m
, int M
, int print_all
);
108 int main(int argc
, char **argv
)
113 int exist
, nparam
, nvar
;
126 int m
= INT_MAX
, M
= INT_MIN
, r
;
127 int print_solution
= 1;
129 while ((c
= getopt_long(argc
, argv
,
130 OMEGA_OPT PIPLIB_OPT
"sSfcRTAm:M:r:V", options
, &ind
)) != -1) {
173 printf(barvinok_version());
179 if (series
&& !scarf
) {
181 "--series currently only available if --scarf is specified\n");
186 A
= Constraints2Polyhedron(MA
, MAXRAYS
);
189 fgets(s
, 128, stdin
);
190 while ((*s
=='#') || (sscanf(s
, "E %d", &exist
)<1))
191 fgets(s
, 128, stdin
);
193 fgets(s
, 128, stdin
);
194 while ((*s
=='#') || (sscanf(s
, "P %d", &nparam
)<1))
195 fgets(s
, 128, stdin
);
197 /******* Read the options: initialize Min and Max ********/
198 if (A
->Dimension
>= VBIGDIM
)
200 else if (A
->Dimension
>= BIGDIM
)
209 if (verify
&& m
> M
) {
210 fprintf(stderr
,"Nothing to do: min > max !\n");
216 if (print_solution
) {
217 Polyhedron_Print(stdout
, P_VALUE_FMT
, A
);
218 printf("exist: %d, nparam: %d\n", exist
, nparam
);
220 param_name
= Read_ParamNames(stdin
, nparam
);
221 nvar
= A
->Dimension
- exist
- nparam
;
223 A
= Omega_simplify(A
, exist
, nparam
, param_name
);
225 exist
= A
->Dimension
- nvar
- nparam
;
230 gf
= barvinok_enumerate_scarf_series(A
, exist
, nparam
, MAXRAYS
);
231 if (print_solution
) {
232 gf
->print(std::cout
, nparam
, param_name
);
238 EP
= barvinok_enumerate_scarf(A
, exist
, nparam
, MAXRAYS
);
239 } else if (pip
&& exist
> 0)
240 EP
= barvinok_enumerate_pip(A
, exist
, nparam
, MAXRAYS
);
242 EP
= barvinok_enumerate_e(A
, exist
, nparam
, MAXRAYS
);
246 evalue_range_reduction(EP
);
248 print_evalue(stdout
, EP
, param_name
);
250 fprintf(stderr
, "WARNING: floor conversion not supported\n");
251 evalue_frac2floor(EP
);
253 print_evalue(stdout
, EP
, param_name
);
254 } else if (convert
) {
255 evalue_mod2table(EP
, nparam
);
257 print_evalue(stdout
, EP
, param_name
);
260 verify_results(A
, EP
, exist
, nparam
, m
, M
, print_all
);
261 free_evalue_refs(EP
);
264 Free_ParamNames(param_name
, nparam
);
269 void verify_results(Polyhedron
*P
, evalue
*EP
, int exist
, int nparam
, int m
, int M
,
276 Polyhedron
*C
= Polyhedron_Project(P
, nparam
);
279 value_set_si(min_val
,m
);
280 value_set_si(max_val
,M
);
283 p
= (Value
*)malloc(sizeof(Value
) * (P
->Dimension
+2));
284 for(i
=0;i
<=P
->Dimension
;i
++) {
286 value_set_si(p
[i
],0);
289 value_set_si(p
[i
],1);
291 /* S = scanning list of polyhedra */
292 S
= Polyhedron_Scan(P
, C
, MAXRAYS
& POL_NO_DUAL
? 0 : MAXRAYS
);
295 if (C
->Dimension
> 0) {
296 value_subtract(tmp
,max_val
,min_val
);
297 if (VALUE_TO_INT(tmp
) > 80)
298 st
= 1+(VALUE_TO_INT(tmp
))/80;
301 for(i
=VALUE_TO_INT(min_val
);i
<=VALUE_TO_INT(max_val
);i
+=st
)
308 /******* CHECK NOW *********/
310 if(S
&& !check_poly(S
, C
, EP
, exist
, nparam
, 0, p
, print_all
)) {
311 fprintf(stderr
,"Check failed !\n");
318 for(i
=0;i
<=(P
->Dimension
+1);i
++)
326 /****************************************************/
327 /* function check_poly : */
328 /* scans the parameter space from min to max (all */
329 /* directions). Computes the number of points in */
330 /* the polytope using both methods, and compare them*/
331 /* returns 1 on success */
332 /****************************************************/
334 int check_poly(Polyhedron
*S
, Polyhedron
*C
, evalue
*EP
,
335 int exist
, int nparam
, int pos
, Value
*z
, int print_all
)
340 value_init(c
); value_init(tmp
);
344 /* Computes the ehrhart polynomial */
345 value_set_double(c
, compute_evalue(EP
,&z
[S
->Dimension
-nparam
+1])+.25);
346 /* if c=0 we may be out of context. */
347 /* scanning is useless in this case*/
348 if(!in_domain(C
,&z
[S
->Dimension
-nparam
+1])) {
356 value_print(stdout
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
357 for(k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
359 value_print(stdout
,VALUE_FMT
,z
[k
]);
362 value_print(stdout
,VALUE_FMT
,c
);
366 /* Manually count the number of points */
367 count_points_e(1, S
, exist
, nparam
, z
, &tmp
);
369 printf(", count = ");
370 value_print(stdout
, P_VALUE_FMT
, tmp
);
374 if(value_ne(tmp
,c
)) {
377 fprintf(stderr
,"Error !\n");
378 fprintf(stderr
,"EP( ");
379 value_print(stderr
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
380 for(k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
381 fprintf(stderr
,", ");
382 value_print(stderr
,VALUE_FMT
,z
[k
]);
384 fprintf(stderr
," ) should be ");
385 value_print(stderr
,VALUE_FMT
,tmp
);
386 fprintf(stderr
,", while EP eval gives ");
387 value_print(stderr
,VALUE_FMT
,c
);
388 fprintf(stderr
,".\n");
389 print_evalue(stderr
, EP
, params
);
390 #ifndef DONT_BREAK_ON_ERROR
391 value_clear(c
); value_clear(tmp
);
400 for(value_assign(tmp
,min_val
); value_le(tmp
,max_val
); value_increment(tmp
,tmp
)) {
402 k
= VALUE_TO_INT(tmp
);
403 if(!pos
&& !(k
%st
)) {
409 value_assign(z
[pos
+S
->Dimension
-nparam
+1],tmp
);
410 if(!check_poly(S
, C
, EP
, exist
, nparam
, pos
+1, z
, print_all
)) {
411 value_clear(c
); value_clear(tmp
);
415 value_clear(c
); value_clear(tmp
);