1 /*************************************************/
3 /* program to compare effective number of points */
4 /* in a polytope with the corresponding */
5 /* evaluation of the Ehrhart polynomial. */
6 /* Parameters vary in range -RANGE to RANGE */
7 /* (define below) by default. */
8 /* Can be overridden by specifying */
9 /* -r<RANGE>, or -m<min> and -M<max> */
11 /* written by Vincent Loechner (c) 2000. */
12 /* loechner@icps.u-strasbg.fr */
13 /*************************************************/
15 This file is part of PolyLib.
17 PolyLib is free software: you can redistribute it and/or modify
18 it under the terms of the GNU General Public License as published by
19 the Free Software Foundation, either version 3 of the License, or
20 (at your option) any later version.
22 PolyLib is distributed in the hope that it will be useful,
23 but WITHOUT ANY WARRANTY; without even the implied warranty of
24 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 GNU General Public License for more details.
27 You should have received a copy of the GNU General Public License
28 along with PolyLib. If not, see <http://www.gnu.org/licenses/>.
36 #include <polylib/polylib.h>
39 /* define this to print all the results */
40 /* else, only a progress bar is printed */
41 /* #define PRINT_ALL_RESULTS */
42 /* define this to continue the test after first error found */
43 /* #define DONT_BREAK_ON_ERROR */
45 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
48 /* SRANGE : small range for evalutations */
51 /* if dimension >= BIDDIM, use SRANGE */
54 /* VSRANGE : very small range for evalutations */
57 /* if dimension >= VBIDDIM, use VSRANGE */
62 #ifdef DONT_BREAK_ON_ERROR
63 #define PRINT_ALL_RESULTS
66 #ifndef PRINT_ALL_RESULTS
70 /****************************************************/
71 /* function check_poly : */
72 /* scans the parameter space from min to max (all */
73 /* directions). Computes the number of points in */
74 /* the polytope using both methods, and compare them*/
75 /* returns 1 on success */
76 /****************************************************/
78 int check_poly(Polyhedron
*S
,Polyhedron
*C
,Enumeration
*en
,
79 int nparam
,int pos
,Value
*z
) {
84 value_init(c
); value_init(tmp
);
88 /* Computes the ehrhart polynomial */
89 value_assign(c
,*(ctmp
=compute_poly(en
,&z
[S
->Dimension
-nparam
+1])));
91 /* if c=0 we may be out of context. */
92 /* scanning is useless in this case*/
93 if(!in_domain(C
,&z
[S
->Dimension
-nparam
+1])) {
99 #ifdef PRINT_ALL_RESULTS
101 value_print(stdout
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
102 for(k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
104 value_print(stdout
,VALUE_FMT
,z
[k
]);
107 value_print(stdout
,VALUE_FMT
,c
);
111 /* Count manually the number of points */
112 count_points(1,S
,z
,&tmp
);
113 #ifdef PRINT_ALL_RESULTS
114 printf(", count = ");
115 value_print(stdout
, P_VALUE_FMT
, tmp
);
119 if(value_ne(tmp
,c
)) {
122 fprintf(stderr
,"Error !\n");
123 fprintf(stderr
,"EP( ");
124 value_print(stderr
,VALUE_FMT
,z
[S
->Dimension
-nparam
+1]);
125 for(k
=S
->Dimension
-nparam
+2;k
<=S
->Dimension
;++k
) {
126 fprintf(stderr
,", ");
127 value_print(stderr
,VALUE_FMT
,z
[k
]);
129 fprintf(stderr
," ) should be ");
130 value_print(stderr
,VALUE_FMT
,tmp
);
131 fprintf(stderr
,", while EP eval gives ");
132 value_print(stderr
,VALUE_FMT
,c
);
133 fprintf(stderr
,".\n");
134 #ifndef DONT_BREAK_ON_ERROR
135 value_clear(c
); value_clear(tmp
);
140 #ifdef PRINT_ALL_RESULTS
147 for(value_assign(tmp
,min
); value_le(tmp
,max
); value_increment(tmp
,tmp
)) {
149 #ifndef PRINT_ALL_RESULTS
150 k
= VALUE_TO_INT(tmp
);
151 if(!pos
&& !(k
%st
)) {
157 value_assign(z
[pos
+S
->Dimension
-nparam
+1],tmp
);
158 if(!check_poly(S
,C
,en
,nparam
,pos
+1,z
)) {
159 value_clear(c
); value_clear(tmp
);
163 value_clear(c
); value_clear(tmp
);
167 int main(int argc
,char *argv
[]) {
170 Polyhedron
*C
, *P
, *S
;
177 /******* Read the input *********/
181 if(C1
->NbColumns
< 2) {
182 fprintf(stderr
,"Not enough parameters !\n");
186 P
= Constraints2Polyhedron(P1
,MAXRAYS
);
187 C
= Constraints2Polyhedron(C1
,MAXRAYS
);
191 /******* Read the options: initialize min and max ********/
192 if(P
->Dimension
>= VBIGDIM
)
194 else if(P
->Dimension
>= BIGDIM
)
200 for(i
=1;i
<argc
;i
++) {
201 if(!strncmp(argv
[i
],"-m",2)) {
204 m
= atoi(&argv
[i
][2]);
206 else if(!strncmp(argv
[i
],"-M",2)) {
209 M
= atoi(&argv
[i
][2]);
211 else if(!strncmp(argv
[i
], "-r", 2)) {
213 /* range specified */
214 M
= atoi(&argv
[i
][2]);
218 fprintf(stderr
,"Unknown option: %s\n",argv
[i
]);
219 fprintf(stderr
,"Usage: %s [-m<>][-M<>][-r<>]\n",argv
[0]);
225 fprintf(stderr
,"Nothing to do: min > max !\n");
234 /******* Compute true context *******/
235 CC
= align_context(C
,P
->Dimension
,MAXRAYS
);
236 PP
= DomainIntersection(P
,CC
,MAXRAYS
);
238 C1
= Matrix_Alloc(C
->Dimension
+1,P
->Dimension
+1);
240 for(i
=0;i
<C1
->NbRows
;i
++)
241 for(j
=0;j
<C1
->NbColumns
;j
++)
242 if(i
==j
-P
->Dimension
+C
->Dimension
)
243 value_set_si(C1
->p
[i
][j
],1);
245 value_set_si(C1
->p
[i
][j
],0);
246 CC
= Polyhedron_Image(PP
,C1
,MAXRAYS
);
250 /******* Compute EP *********/
251 en
= Polyhedron_Enumerate(P
,C
,MAXRAYS
,NULL
);
253 /******* Initializations for check *********/
254 p
= (Value
*)malloc(sizeof(Value
) * (P
->Dimension
+2));
255 for(i
=0;i
<=P
->Dimension
;i
++) {
257 value_set_si(p
[i
],0);
260 value_set_si(p
[i
],1);
262 /* S = scanning list of polyhedra */
263 S
= Polyhedron_Scan(P
,C
,MAXRAYS
);
265 #ifndef PRINT_ALL_RESULTS
266 if(C
->Dimension
> 0) {
267 value_subtract(tmp
,max
,min
);
268 if (VALUE_TO_INT(tmp
) > 80)
269 st
= 1+(VALUE_TO_INT(tmp
))/80;
272 for(i
=VALUE_TO_INT(min
);i
<=VALUE_TO_INT(max
);i
+=st
)
279 /******* CHECK NOW *********/
280 if(S
&& !check_poly(S
,C
,en
,C
->Dimension
,0,p
)) {
281 fprintf(stderr
,"Check failed !\n");
282 for(i
=0;i
<=(P
->Dimension
+1);i
++)
288 #ifndef PRINT_ALL_RESULTS
292 for(i
=0;i
<=(P
->Dimension
+1);i
++)