verify_main.cc

   1 #include <gmp.h>
   2 #include <NTL/mat_ZZ.h>
   3 extern "C" {
   4 #include <polylib/polylibgmp.h>
   5 #include "ev_operations.h"
   6 }
   7 #include <util.h>
   8 #include <barvinok.h>
   9 #include <barvinok2.h>
  10
  11 #include "verif_ehrhart.h"
  12
  13 #ifdef HAVE_GROWING_CHERNIKOVA
  14 #define MAXRAYS    0
  15 #else
  16 #define MAXRAYS  600
  17 #endif
  18
  19 #include "config.h"
  20 #ifndef HAVE_GETOPT_H
  21 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
  22 #else
  23 #include <getopt.h>
  24 struct option options[] = {
  25     { "explicit",  no_argument,  0,  'e' },
  26     { "series",  no_argument,  0,  's' },
  27     { "verbose",  no_argument,  0,  'v' },
  28     { "version",   no_argument,  0,  'V' },
  29     { 0, 0, 0, 0 }
  30 };
  31 #endif
  32
  33 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
  34 #define RANGE 50
  35
  36 /* SRANGE : small range for evalutations */
  37 #define SRANGE 15
  38
  39 /* if dimension >= BIDDIM, use SRANGE */
  40 #define BIGDIM 5
  41
  42 /* VSRANGE : very small range for evalutations */
  43 #define VSRANGE 5
  44
  45 /* if dimension >= VBIDDIM, use VSRANGE */
  46 #define VBIGDIM 8
  47
  48 int check_series(Polyhedron *S,Polyhedron *C, gen_fun *gf,
  49                  int nparam,int pos,Value *z)
  50 {
  51     int k;
  52     Value c, tmp;
  53     value_init(c);
  54     value_init(tmp);
  55
  56     if(pos == nparam) {
  57         /* Computes the coefficient */
  58         gf->coefficient(&z[S->Dimension-nparam+1], &c);
  59
  60         /* if c=0 we may be out of context. */
  61         /* scanning is useless in this case*/
  62         if(!in_domain(C,&z[S->Dimension-nparam+1])) {
  63             /* ok */ ;
  64         } else {
  65
  66 #ifdef PRINT_ALL_RESULTS
  67             printf("EP( ");
  68             value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
  69             for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  70                 printf(", ");
  71                 value_print(stdout,VALUE_FMT,z[k]);
  72             }
  73             printf(" ) = ");
  74             value_print(stdout,VALUE_FMT,c);
  75             printf(" ");
  76 #endif
  77             /* Manually count the number of points */
  78             count_points(1,S,z,&tmp);
  79 #ifdef PRINT_ALL_RESULTS
  80             printf(", count = ");
  81             value_print(stdout, P_VALUE_FMT, tmp);
  82             printf(". ");
  83 #endif
  84
  85             if(value_ne(tmp,c)) {
  86                 printf("\n");
  87                 fflush(stdout);
  88                 fprintf(stderr,"Error !\n");
  89                 fprintf(stderr,"EP( ");
  90                 value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
  91                 for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  92                   fprintf(stderr,", ");
  93                   value_print(stderr,VALUE_FMT,z[k]);
  94                 }
  95                 fprintf(stderr," ) should be ");
  96                 value_print(stderr,VALUE_FMT,tmp);
  97                 fprintf(stderr,", while EP eval gives ");
  98                 value_print(stderr,VALUE_FMT,c);
  99                 fprintf(stderr,".\n");
 100 #ifndef DONT_BREAK_ON_ERROR
 101                 value_clear(c); value_clear(tmp);
 102                 return 0;
 103 #endif
 104             }
 105 #ifdef PRINT_ALL_RESULTS
 106             else
 107                 printf("OK.\n");
 108 #endif
 109         }
 110     } else
 111         for(value_assign(tmp,Min); value_le(tmp,Max);
 112                                    value_increment(tmp,tmp)) {
 113 #ifndef PRINT_ALL_RESULTS
 114           k = VALUE_TO_INT(tmp);
 115           if(!pos && !(k%st)) {
 116             printf("o");
 117             fflush(stdout);
 118           }
 119 #endif
 120           value_assign(z[pos+S->Dimension-nparam+1],tmp);
 121           if(!check_series(S,C,gf,nparam,pos+1,z)) {
 122             value_clear(c); value_clear(tmp);
 123             return(0);
 124           }
 125         }
 126
 127     value_clear(c);
 128     value_clear(tmp);
 129     return 1;
 130 }
 131
 132 int main(int argc,char *argv[]) {
 133
 134   Matrix *C1, *P1;
 135   Polyhedron *C, *P, *S;
 136   Polyhedron *CC, *PP;
 137   Enumeration *en;
 138   Value *p, tmp;
 139   int i,j;
 140     int m = INT_MAX, M = INT_MIN, r;
 141     int c, ind = 0;
 142     int series = 0;
 143     int verbose = 0;
 144     int function = 0;
 145     int result = 0;
 146
 147     while ((c = getopt_long(argc, argv, "m:M:r:sveV", options, &ind)) != -1) {
 148         switch (c) {
 149         case 'e':
 150             function = 1;
 151             break;
 152         case 's':
 153             series = 1;
 154             break;
 155         case 'v':
 156             verbose = 1;
 157             break;
 158         case 'm':
 159             m = atoi(optarg);
 160             break;
 161         case 'M':
 162             M = atoi(optarg);
 163             break;
 164         case 'r':
 165             M = atoi(optarg);
 166             m = -M;
 167             break;
 168         case 'V':
 169             printf(barvinok_version());
 170             exit(0);
 171             break;
 172         }
 173     }
 174
 175 /******* Read the input *********/
 176   P1 = Matrix_Read();
 177   C1 = Matrix_Read();
 178
 179   if(C1->NbColumns < 2) {
 180     fprintf(stderr,"Not enough parameters !\n");
 181     exit(0);
 182   }
 183
 184   P = Constraints2Polyhedron(P1,MAXRAYS);
 185   C = Constraints2Polyhedron(C1,MAXRAYS);
 186   params = Read_ParamNames(stdin, C->Dimension);
 187   Matrix_Free(C1);
 188   Matrix_Free(P1);
 189
 190   /******* Read the options: initialize Min and Max ********/
 191   if(P->Dimension >= VBIGDIM)
 192     r = VSRANGE;
 193   else if(P->Dimension >= BIGDIM)
 194     r = SRANGE;
 195   else
 196     r = RANGE;
 197   if (M == INT_MIN)
 198     M = r;
 199   if (m == INT_MAX)
 200     m = -r;
 201
 202
 203   if(m > M) {
 204     fprintf(stderr,"Nothing to do: Min > Max !\n");
 205     return(0);
 206   }
 207   value_init(Min);
 208   value_init(Max);
 209   value_set_si(Min,m);
 210   value_set_si(Max,M);
 211   value_init(tmp);
 212
 213   /******* Compute true context *******/
 214   CC = align_context(C,P->Dimension,MAXRAYS);
 215   PP = DomainIntersection(P,CC,MAXRAYS);
 216   Domain_Free(CC);
 217   C1 = Matrix_Alloc(C->Dimension+1,P->Dimension+1);
 218
 219   for(i=0;i<C1->NbRows;i++)
 220     for(j=0;j<C1->NbColumns;j++)
 221       if(i==j-P->Dimension+C->Dimension)
 222         value_set_si(C1->p[i][j],1);
 223       else
 224         value_set_si(C1->p[i][j],0);
 225   CC = Polyhedron_Image(PP,C1,MAXRAYS);
 226   Domain_Free(C);
 227   C = CC;
 228
 229     gen_fun *gf = 0;
 230
 231     /******* Compute EP *********/
 232     if (!series)
 233         en = barvinok_enumerate(P,C,MAXRAYS);
 234     else {
 235         evalue *EP;
 236         gf = barvinok_series(P, C, MAXRAYS);
 237         if (verbose) {
 238             gf->print(C->Dimension, params);
 239             puts("");
 240         }
 241         if (function) {
 242             EP = *gf;
 243             en =  partition2enumeration(EP);
 244         }
 245     }
 246
 247   /******* Initializations for check *********/
 248   p = (Value *)malloc(sizeof(Value) * (P->Dimension+2));
 249   for(i=0;i<=P->Dimension;i++) {
 250     value_init(p[i]);
 251     value_set_si(p[i],0);
 252   }
 253   value_init(p[i]);
 254   value_set_si(p[i],1);
 255
 256   /* S = scanning list of polyhedra */
 257   S = Polyhedron_Scan(P,C,MAXRAYS);
 258
 259 #ifndef PRINT_ALL_RESULTS
 260   if(C->Dimension > 0) {
 261     value_substract(tmp,Max,Min);
 262     if (VALUE_TO_INT(tmp) > 80)
 263       st = 1+(VALUE_TO_INT(tmp))/80;
 264     else
 265       st=1;
 266     for(i=VALUE_TO_INT(Min);i<=VALUE_TO_INT(Max);i+=st)
 267       printf(".");
 268     printf( "\r" );
 269     fflush(stdout);
 270   }
 271 #endif
 272
 273     /******* CHECK NOW *********/
 274     if(S) {
 275         if (!series || function) {
 276             if (!check_poly(S,C,en,C->Dimension,0,p))
 277                 result = -1;
 278         } else {
 279             if (!check_series(S,C,gf,C->Dimension,0,p))
 280                 result = -1;
 281         }
 282     }
 283
 284     if (result == -1)
 285         fprintf(stderr,"Check failed !\n");
 286
 287 #ifndef PRINT_ALL_RESULTS
 288   printf( "\n" );
 289 #endif
 290
 291     if (gf)
 292         delete gf;
 293
 294   for(i=0;i<=(P->Dimension+1);i++)
 295     value_clear(p[i]);
 296   value_clear(tmp);
 297   Free_ParamNames(params, C->Dimension);
 298   return result;
 299 } /* main */
 300
 301
 302
 303