applications/verif_ehrhart.c

   1 /*************************************************/
   2 /*     verif_ehrhart.c                           */
   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>             */
  10 /*                                               */
  11 /* written by Vincent Loechner (c) 2000.         */
  12 /*  loechner@icps.u-strasbg.fr                   */
  13 /*************************************************/
  14
  15 #include <stdio.h>
  16 #include <string.h>
  17 #include <stdlib.h>
  18
  19 #include <polylib/polylib.h>
  20 #define MAXRAYS 1024
  21
  22 /* define this to print all the results */
  23 /* else, only a progress bar is printed */
  24 /* #define PRINT_ALL_RESULTS     */
  25 /* define this to continue the test after first error found */
  26 /* #define DONT_BREAK_ON_ERROR */
  27
  28 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
  29 #define RANGE 50
  30
  31 /* SRANGE : small range for evalutations */
  32 #define SRANGE 15
  33
  34 /* if dimension >= BIDDIM, use SRANGE */
  35 #define BIGDIM 5
  36
  37 /* VSRANGE : very small range for evalutations */
  38 #define VSRANGE 5
  39
  40 /* if dimension >= VBIDDIM, use VSRANGE */
  41 #define VBIGDIM 8
  42
  43 Value min, max;
  44
  45 #ifdef DONT_BREAK_ON_ERROR
  46 #define PRINT_ALL_RESULTS
  47 #endif
  48
  49 #ifndef PRINT_ALL_RESULTS
  50 int st;
  51 #endif
  52
  53 /****************************************************/
  54 /* function check_poly :                            */
  55 /* scans the parameter space from min to max (all   */
  56 /* directions). Computes the number of points in    */
  57 /* the polytope using both methods, and compare them*/
  58 /* returns 1 on success                             */
  59 /****************************************************/
  60
  61 int check_poly(Polyhedron *S,Polyhedron *C,Enumeration *en,
  62                int nparam,int pos,Value *z) {
  63
  64   int cc,k;
  65   Value c,tmp,*ctmp;
  66
  67   value_init(c); value_init(tmp);
  68
  69   if(pos == nparam) {
  70
  71     /* Computes the ehrhart polynomial */
  72     value_assign(c,*(ctmp=compute_poly(en,&z[S->Dimension-nparam+1])));
  73     free(ctmp);
  74     /* if c=0 we may be out of context. */
  75     /* scanning is useless in this case*/
  76     if(!in_domain(C,&z[S->Dimension-nparam+1])) {
  77
  78       /* ok */ ;
  79     }
  80     else {
  81
  82 #ifdef PRINT_ALL_RESULTS
  83       printf("EP( ");
  84       value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
  85       for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  86         printf(", ");
  87         value_print(stdout,VALUE_FMT,z[k]);
  88       }
  89       printf(" ) = ");
  90       value_print(stdout,VALUE_FMT,c);
  91       printf(" ");
  92 #endif
  93
  94       /* Count manually the number of points */
  95       count_points(1,S,z,&tmp);
  96 #ifdef PRINT_ALL_RESULTS
  97         printf(", count = ");
  98         value_print(stdout, P_VALUE_FMT, tmp);
  99         printf(". ");
 100 #endif
 101
 102       if(value_ne(tmp,c)) {
 103         printf("\n");
 104         fflush(stdout);
 105         fprintf(stderr,"Error !\n");
 106         fprintf(stderr,"EP( ");
 107         value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
 108         for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
 109           fprintf(stderr,", ");
 110           value_print(stderr,VALUE_FMT,z[k]);
 111         }
 112         fprintf(stderr," ) should be ");
 113         value_print(stderr,VALUE_FMT,tmp);
 114         fprintf(stderr,", while EP eval gives ");
 115         value_print(stderr,VALUE_FMT,c);
 116         fprintf(stderr,".\n");
 117 #ifndef DONT_BREAK_ON_ERROR
 118         value_clear(c); value_clear(tmp);
 119         return(0);
 120 #endif
 121       }
 122
 123 #ifdef PRINT_ALL_RESULTS
 124       else
 125         printf("OK.\n");
 126 #endif
 127     }
 128   }
 129   else
 130     for(value_assign(tmp,min); value_le(tmp,max); value_increment(tmp,tmp)) {
 131
 132 #ifndef PRINT_ALL_RESULTS
 133       k = VALUE_TO_INT(tmp);
 134       if(!pos && !(k%st)) {
 135         printf("o");
 136         fflush(stdout);
 137       }
 138 #endif
 139
 140       value_assign(z[pos+S->Dimension-nparam+1],tmp);
 141       if(!check_poly(S,C,en,nparam,pos+1,z)) {
 142         value_clear(c); value_clear(tmp);
 143         return(0);
 144       }
 145     }
 146   value_clear(c); value_clear(tmp);
 147   return(1);
 148 } /* check_poly */
 149
 150 int main(int argc,char *argv[]) {
 151
 152   Matrix *C1, *P1;
 153   Polyhedron *C, *P, *S;
 154   Polyhedron *CC, *PP;
 155   Enumeration *en;
 156   Value *p, tmp;
 157   int i,j;
 158   int m,M;
 159
 160 /******* Read the input *********/
 161   P1 = Matrix_Read();
 162   C1 = Matrix_Read();
 163
 164   if(C1->NbColumns < 2) {
 165     fprintf(stderr,"Not enough parameters !\n");
 166     exit(0);
 167   }
 168
 169   P = Constraints2Polyhedron(P1,MAXRAYS);
 170   C = Constraints2Polyhedron(C1,MAXRAYS);
 171   Matrix_Free(C1);
 172   Matrix_Free(P1);
 173
 174   /******* Read the options: initialize min and max ********/
 175   if(P->Dimension >= VBIGDIM)
 176     M = VSRANGE;
 177   else if(P->Dimension >= BIGDIM)
 178     M = SRANGE;
 179   else
 180     M = RANGE;
 181   m = -M;
 182   if(argc != 1 ) {
 183     for(i=1;i<argc;i++) {
 184       if(!strncmp(argv[i],"-m",2)) {
 185
 186         /* min specified */
 187         m = atoi(&argv[i][2]);
 188       }
 189       else if(!strncmp(argv[i],"-M",2)) {
 190
 191         /* max specified */
 192         M = atoi(&argv[i][2]);
 193       }
 194       else if(!strncmp(argv[i], "-r", 2)) {
 195
 196         /* range specified */
 197         M = atoi(&argv[i][2]);
 198         m = -M;
 199       }
 200       else {
 201         fprintf(stderr,"Unknown option: %s\n",argv[i]);
 202         fprintf(stderr,"Usage: %s [-m<>][-M<>][-r<>]\n",argv[0]);
 203         return(-1);
 204       }
 205     }
 206   }
 207   if(m > M) {
 208     fprintf(stderr,"Nothing to do: min > max !\n");
 209     return(0);
 210   }
 211   value_init(min);
 212   value_init(max);
 213   value_set_si(min,m);
 214   value_set_si(max,M);
 215   value_init(tmp);
 216
 217   /******* Compute true context *******/
 218   CC = align_context(C,P->Dimension,MAXRAYS);
 219   PP = DomainIntersection(P,CC,MAXRAYS);
 220   Domain_Free(CC);
 221   C1 = Matrix_Alloc(C->Dimension+1,P->Dimension+1);
 222
 223   for(i=0;i<C1->NbRows;i++)
 224     for(j=0;j<C1->NbColumns;j++)
 225       if(i==j-P->Dimension+C->Dimension)
 226         value_set_si(C1->p[i][j],1);
 227       else
 228         value_set_si(C1->p[i][j],0);
 229   CC = Polyhedron_Image(PP,C1,MAXRAYS);
 230   Domain_Free(C);
 231   C = CC;
 232
 233   /******* Compute EP *********/
 234   en = Polyhedron_Enumerate(P,C,MAXRAYS,NULL);
 235
 236   /******* Initializations for check *********/
 237   p = (Value *)malloc(sizeof(Value) * (P->Dimension+2));
 238   for(i=0;i<=P->Dimension;i++) {
 239     value_init(p[i]);
 240     value_set_si(p[i],0);
 241   }
 242   value_init(p[i]);
 243   value_set_si(p[i],1);
 244
 245   /* S = scanning list of polyhedra */
 246   S = Polyhedron_Scan(P,C,MAXRAYS);
 247
 248 #ifndef PRINT_ALL_RESULTS
 249   if(C->Dimension > 0) {
 250     value_subtract(tmp,max,min);
 251     if (VALUE_TO_INT(tmp) > 80)
 252       st = 1+(VALUE_TO_INT(tmp))/80;
 253     else
 254       st=1;
 255     for(i=VALUE_TO_INT(min);i<=VALUE_TO_INT(max);i+=st)
 256       printf(".");
 257     printf( "\r" );
 258     fflush(stdout);
 259   }
 260 #endif
 261
 262   /******* CHECK NOW *********/
 263   if(S && !check_poly(S,C,en,C->Dimension,0,p)) {
 264     fprintf(stderr,"Check failed !\n");
 265     for(i=0;i<=(P->Dimension+1);i++)
 266       value_clear(p[i]);
 267     value_clear(tmp);
 268     return(-1);
 269   }
 270
 271 #ifndef PRINT_ALL_RESULTS
 272   printf( "\n" );
 273 #endif
 274
 275   for(i=0;i<=(P->Dimension+1);i++)
 276     value_clear(p[i]);
 277   value_clear(tmp);
 278   return(0);
 279 } /* main */
 280
 281
 282
 283