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