3 #include <polylib/polylibgmp.h>
4 #include <barvinok/util.h>
5 #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 { "convert", no_argument
, 0, 'c' },
38 { "floor", no_argument
, 0, 'f' },
39 { "range-reduction", no_argument
, 0, 'R' },
40 { "verify", no_argument
, 0, 'T' },
41 { "min", required_argument
, 0, 'm' },
42 { "max", required_argument
, 0, 'M' },
43 { "range", required_argument
, 0, 'r' },
44 { "version", no_argument
, 0, 'V' },
50 #define PIPLIB_OPT "p"
58 Polyhedron
*Omega_simplify(Polyhedron
*P
,
59 unsigned exist
, unsigned nparam
, char **parms
)
63 Relation r
= Polyhedron2relation(P
, exist
, nparam
, parms
);
65 return relation2Domain(r
, varv
, paramv
);
69 Polyhedron
*Omega_simplify(Polyhedron
*P
,
70 unsigned exist
, unsigned nparam
, char **parms
)
76 /* define this to print all the results */
77 /* else, only a progress bar is printed */
78 /* #define PRINT_ALL_RESULTS
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
;
101 #ifdef DONT_BREAK_ON_ERROR
102 #define PRINT_ALL_RESULTS
105 #ifndef PRINT_ALL_RESULTS
109 static int check_poly(Polyhedron
*S
, Polyhedron
*C
, evalue
*EP
,
110 int exist
, int nparam
, int pos
, Value
*z
);
111 static void verify_results(Polyhedron
*P
, evalue
*EP
, int exist
, int nparam
,
114 int main(int argc
, char **argv
)
119 int exist
, nparam
, nvar
;
129 int m
= INT_MAX
, M
= INT_MIN
, r
;
130 int print_solution
= 1;
132 while ((c
= getopt_long(argc
, argv
,
133 OMEGA_OPT PIPLIB_OPT
"fcRTm:M:r:V", options
, &ind
)) != -1) {
167 printf(barvinok_version());
174 A
= Constraints2Polyhedron(MA
, MAXRAYS
);
177 fgets(s
, 128, stdin
);
178 while ((*s
=='#') || (sscanf(s
, "E %d", &exist
)<1))
179 fgets(s
, 128, stdin
);
181 fgets(s
, 128, stdin
);
182 while ((*s
=='#') || (sscanf(s
, "P %d", &nparam
)<1))
183 fgets(s
, 128, stdin
);
185 /******* Read the options: initialize Min and Max ********/
186 if (A
->Dimension
>= VBIGDIM
)
188 else if (A
->Dimension
>= BIGDIM
)
197 if (verify
&& m
> M
) {
198 fprintf(stderr
,"Nothing to do: min > max !\n");
204 if (print_solution
) {
205 Polyhedron_Print(stdout
, P_VALUE_FMT
, A
);
206 printf("exist: %d, nparam: %d\n", exist
, nparam
);
208 param_name
= Read_ParamNames(stdin
, nparam
);
209 nvar
= A
->Dimension
- exist
- nparam
;
211 A
= Omega_simplify(A
, exist
, nparam
, param_name
);
213 exist
= A
->Dimension
- nvar
- nparam
;
215 if (pip
&& exist
> 0)
216 EP
= barvinok_enumerate_pip(A
, exist
, nparam
, MAXRAYS
);
218 EP
= barvinok_enumerate_e(A
, exist
, nparam
, MAXRAYS
);
222 evalue_range_reduction(EP
);
224 print_evalue(stdout
, EP
, param_name
);
226 fprintf(stderr
, "WARNING: floor conversion not supported\n");
227 evalue_frac2floor(EP
);
229 print_evalue(stdout
, EP
, param_name
);
230 } else if (convert
) {
231 evalue_mod2table(EP
, nparam
);
233 print_evalue(stdout
, EP
, param_name
);
236 verify_results(A
, EP
, exist
, nparam
, m
, M
);
237 free_evalue_refs(EP
);
239 Free_ParamNames(param_name
, nparam
);
244 void verify_results(Polyhedron
*P
, evalue
*EP
, int exist
, int nparam
, int m
, int M
)
250 Polyhedron
*C
= Polyhedron_Project(P
, nparam
);
253 value_set_si(min_val
,m
);
254 value_set_si(max_val
,M
);
257 p
= (Value
*)malloc(sizeof(Value
) * (P
->Dimension
+2));
258 for(i
=0;i
<=P
->Dimension
;i
++) {
260 value_set_si(p
[i
],0);
263 value_set_si(p
[i
],1);
265 /* S = scanning list of polyhedra */
266 S
= Polyhedron_Scan(P
, C
, MAXRAYS
& POL_NO_DUAL
? 0 : MAXRAYS
);
268 #ifndef PRINT_ALL_RESULTS
269 if (C
->Dimension
> 0) {
270 value_subtract(tmp
,max_val
,min_val
);
271 if (VALUE_TO_INT(tmp
) > 80)
272 st
= 1+(VALUE_TO_INT(tmp
))/80;
275 for(i
=VALUE_TO_INT(min_val
);i
<=VALUE_TO_INT(max_val
);i
+=st
)
282 /******* CHECK NOW *********/
284 if(S
&& !check_poly(S
, C
, EP
, exist
, nparam
, 0, p
)) {
285 fprintf(stderr
,"Check failed !\n");
289 #ifndef PRINT_ALL_RESULTS
293 for(i
=0;i
<=(P
->Dimension
+1);i
++)
301 /****************************************************/
302 /* function check_poly : */
303 /* scans the parameter space from min to max (all */
304 /* directions). Computes the number of points in */
305 /* the polytope using both methods, and compare them*/
306 /* returns 1 on success */
307 /****************************************************/
309 int check_poly(Polyhedron
*S
, Polyhedron
*C
, evalue
*EP
,
310 int exist
, int nparam
, int pos
, Value
*z
)
315 value_init(c
); value_init(tmp
);
319 /* Computes the ehrhart polynomial */
320 value_set_double(c
, compute_evalue(EP
,&z
[S
->Dimension
-nparam
+1])+.25);
321 /* if c=0 we may be out of context. */
322 /* scanning is useless in this case*/
323 if(!in_domain(C
,&z
[S
->Dimension
-nparam
+1])) {
329 #ifdef PRINT_ALL_RESULTS
331 value_print(stdout
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
332 for(k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
334 value_print(stdout
,VALUE_FMT
,z
[k
]);
337 value_print(stdout
,VALUE_FMT
,c
);
341 /* Manually count the number of points */
342 count_points_e(1, S
, exist
, nparam
, z
, &tmp
);
343 #ifdef PRINT_ALL_RESULTS
344 printf(", count = ");
345 value_print(stdout
, P_VALUE_FMT
, tmp
);
349 if(value_ne(tmp
,c
)) {
352 fprintf(stderr
,"Error !\n");
353 fprintf(stderr
,"EP( ");
354 value_print(stderr
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
355 for(k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
356 fprintf(stderr
,", ");
357 value_print(stderr
,VALUE_FMT
,z
[k
]);
359 fprintf(stderr
," ) should be ");
360 value_print(stderr
,VALUE_FMT
,tmp
);
361 fprintf(stderr
,", while EP eval gives ");
362 value_print(stderr
,VALUE_FMT
,c
);
363 fprintf(stderr
,".\n");
364 print_evalue(stderr
, EP
, params
);
365 #ifndef DONT_BREAK_ON_ERROR
366 value_clear(c
); value_clear(tmp
);
371 #ifdef PRINT_ALL_RESULTS
378 for(value_assign(tmp
,min_val
); value_le(tmp
,max_val
); value_increment(tmp
,tmp
)) {
380 #ifndef PRINT_ALL_RESULTS
381 k
= VALUE_TO_INT(tmp
);
382 if(!pos
&& !(k
%st
)) {
388 value_assign(z
[pos
+S
->Dimension
-nparam
+1],tmp
);
389 if(!check_poly(S
, C
, EP
, exist
, nparam
, pos
+1, z
)) {
390 value_clear(c
); value_clear(tmp
);
394 value_clear(c
); value_clear(tmp
);