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 *CS, gen_fun *gf,
  49                  int nparam,int pos,Value *z)
  50 {
  51     int k;
  52     Value c, tmp;
  53     Value LB, UB;
  54
  55     value_init(c);
  56     value_init(tmp);
  57     value_init(LB);
  58     value_init(UB);
  59
  60     if(pos == nparam) {
  61         /* Computes the coefficient */
  62         gf->coefficient(&z[S->Dimension-nparam+1], &c);
  63
  64         /* if c=0 we may be out of context. */
  65         /* scanning is useless in this case*/
  66
  67 #ifdef PRINT_ALL_RESULTS
  68             printf("EP( ");
  69             value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
  70             for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  71                 printf(", ");
  72                 value_print(stdout,VALUE_FMT,z[k]);
  73             }
  74             printf(" ) = ");
  75             value_print(stdout,VALUE_FMT,c);
  76             printf(" ");
  77 #endif
  78             /* Manually count the number of points */
  79             count_points(1,S,z,&tmp);
  80 #ifdef PRINT_ALL_RESULTS
  81             printf(", count = ");
  82             value_print(stdout, P_VALUE_FMT, tmp);
  83             printf(". ");
  84 #endif
  85
  86             if(value_ne(tmp,c)) {
  87                 printf("\n");
  88                 fflush(stdout);
  89                 fprintf(stderr,"Error !\n");
  90                 fprintf(stderr,"EP( ");
  91                 value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
  92                 for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  93                   fprintf(stderr,", ");
  94                   value_print(stderr,VALUE_FMT,z[k]);
  95                 }
  96                 fprintf(stderr," ) should be ");
  97                 value_print(stderr,VALUE_FMT,tmp);
  98                 fprintf(stderr,", while EP eval gives ");
  99                 value_print(stderr,VALUE_FMT,c);
 100                 fprintf(stderr,".\n");
 101 #ifndef DONT_BREAK_ON_ERROR
 102                 value_clear(c); value_clear(tmp);
 103                 return 0;
 104 #endif
 105             }
 106 #ifdef PRINT_ALL_RESULTS
 107             else
 108                 printf("OK.\n");
 109 #endif
 110     } else {
 111         int ok =
 112           !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
 113         assert(ok);
 114         for(value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
 115 #ifndef PRINT_ALL_RESULTS
 116           k = VALUE_TO_INT(tmp);
 117           if(!pos && !(k%st)) {
 118             printf("o");
 119             fflush(stdout);
 120           }
 121 #endif
 122           value_assign(z[pos+S->Dimension-nparam+1],tmp);
 123           if(!check_series(S, CS, gf, nparam, pos+1, z)) {
 124             value_clear(c); value_clear(tmp);
 125             value_clear(LB);
 126             value_clear(UB);
 127             return(0);
 128           }
 129         }
 130         value_set_si(z[pos+S->Dimension-nparam+1],0);
 131      }
 132
 133     value_clear(c);
 134     value_clear(tmp);
 135     value_clear(LB);
 136     value_clear(UB);
 137     return 1;
 138 }
 139
 140 int main(int argc,char *argv[]) {
 141
 142   Matrix *C1, *P1, *MM;
 143   Polyhedron *C, *P, *S, *CS, *U;
 144   Polyhedron *CC, *PP;
 145   Enumeration *en;
 146   Value *p, tmp;
 147   Value Min, Max;
 148   int i,j;
 149     int m = INT_MAX, M = INT_MIN, r;
 150     int c, ind = 0;
 151     int series = 0;
 152     int verbose = 0;
 153     int function = 0;
 154     int result = 0;
 155
 156     while ((c = getopt_long(argc, argv, "m:M:r:sveV", options, &ind)) != -1) {
 157         switch (c) {
 158         case 'e':
 159             function = 1;
 160             break;
 161         case 's':
 162             series = 1;
 163             break;
 164         case 'v':
 165             verbose = 1;
 166             break;
 167         case 'm':
 168             m = atoi(optarg);
 169             break;
 170         case 'M':
 171             M = atoi(optarg);
 172             break;
 173         case 'r':
 174             M = atoi(optarg);
 175             m = -M;
 176             break;
 177         case 'V':
 178             printf(barvinok_version());
 179             exit(0);
 180             break;
 181         }
 182     }
 183
 184 /******* Read the input *********/
 185   P1 = Matrix_Read();
 186   C1 = Matrix_Read();
 187
 188   if(C1->NbColumns < 2) {
 189     fprintf(stderr,"Not enough parameters !\n");
 190     exit(0);
 191   }
 192
 193   P = Constraints2Polyhedron(P1,MAXRAYS);
 194   C = Constraints2Polyhedron(C1,MAXRAYS);
 195   params = Read_ParamNames(stdin, C->Dimension);
 196   Matrix_Free(C1);
 197   Matrix_Free(P1);
 198
 199   /******* Read the options: initialize Min and Max ********/
 200   if(P->Dimension >= VBIGDIM)
 201     r = VSRANGE;
 202   else if(P->Dimension >= BIGDIM)
 203     r = SRANGE;
 204   else
 205     r = RANGE;
 206   if (M == INT_MIN)
 207     M = r;
 208   if (m == INT_MAX)
 209     m = -r;
 210
 211
 212   if(m > M) {
 213     fprintf(stderr,"Nothing to do: Min > Max !\n");
 214     return(0);
 215   }
 216   value_init(Min);
 217   value_init(Max);
 218   value_set_si(Min,m);
 219   value_set_si(Max,M);
 220   value_init(tmp);
 221
 222   /******* Compute true context *******/
 223   CC = align_context(C,P->Dimension,MAXRAYS);
 224   PP = DomainIntersection(P,CC,MAXRAYS);
 225   Domain_Free(CC);
 226   C1 = Matrix_Alloc(C->Dimension+1,P->Dimension+1);
 227
 228   for(i=0;i<C1->NbRows;i++)
 229     for(j=0;j<C1->NbColumns;j++)
 230       if(i==j-P->Dimension+C->Dimension)
 231         value_set_si(C1->p[i][j],1);
 232       else
 233         value_set_si(C1->p[i][j],0);
 234   CC = Polyhedron_Image(PP,C1,MAXRAYS);
 235   Matrix_Free(C1);
 236   Domain_Free(PP);
 237   Domain_Free(C);
 238   C = CC;
 239
 240     /* Intersect context with range */
 241     if(C->Dimension > 0) {
 242         MM = Matrix_Alloc(2*C->Dimension, C->Dimension+2);
 243         for (i = 0; i < C->Dimension; ++i) {
 244             value_set_si(MM->p[2*i][0], 1);
 245             value_set_si(MM->p[2*i][1+i], 1);
 246             value_set_si(MM->p[2*i][1+C->Dimension], -m);
 247             value_set_si(MM->p[2*i+1][0], 1);
 248             value_set_si(MM->p[2*i+1][1+i], -1);
 249             value_set_si(MM->p[2*i+1][1+C->Dimension], M);
 250         }
 251         CC = AddConstraints(MM->p[0], 2*C->Dimension, C, MAXRAYS);
 252         U = Universe_Polyhedron(0);
 253         CS = Polyhedron_Scan(CC, U, MAXRAYS);
 254         Polyhedron_Free(U);
 255         Polyhedron_Free(CC);
 256         Matrix_Free(MM);
 257     } else
 258         CS = NULL;
 259
 260     gen_fun *gf = 0;
 261
 262     /******* Compute EP *********/
 263     if (!series)
 264         en = barvinok_enumerate(P,C,MAXRAYS);
 265     else {
 266         evalue *EP;
 267         gf = barvinok_series(P, C, MAXRAYS);
 268         if (verbose) {
 269             gf->print(C->Dimension, params);
 270             puts("");
 271         }
 272         if (function) {
 273             EP = *gf;
 274             en =  partition2enumeration(EP);
 275         }
 276     }
 277
 278   /******* Initializations for check *********/
 279   p = (Value *)malloc(sizeof(Value) * (P->Dimension+2));
 280   for(i=0;i<=P->Dimension;i++) {
 281     value_init(p[i]);
 282     value_set_si(p[i],0);
 283   }
 284   value_init(p[i]);
 285   value_set_si(p[i],1);
 286
 287   /* S = scanning list of polyhedra */
 288   S = Polyhedron_Scan(P,C,MAXRAYS);
 289
 290 #ifndef PRINT_ALL_RESULTS
 291   if(C->Dimension > 0) {
 292     value_subtract(tmp,Max,Min);
 293     if (VALUE_TO_INT(tmp) > 80)
 294       st = 1+(VALUE_TO_INT(tmp))/80;
 295     else
 296       st=1;
 297     for(i=VALUE_TO_INT(Min);i<=VALUE_TO_INT(Max);i+=st)
 298       printf(".");
 299     printf( "\r" );
 300     fflush(stdout);
 301   }
 302 #endif
 303
 304     /******* CHECK NOW *********/
 305     if(S) {
 306         if (!series || function) {
 307             if (!check_poly(S, CS,en,C->Dimension,0,p))
 308                 result = -1;
 309         } else {
 310             if (!check_series(S, CS,gf,C->Dimension,0,p))
 311                 result = -1;
 312         }
 313         Domain_Free(S);
 314     }
 315
 316     if (result == -1)
 317         fprintf(stderr,"Check failed !\n");
 318
 319 #ifndef PRINT_ALL_RESULTS
 320   printf( "\n" );
 321 #endif
 322
 323     if (gf)
 324         delete gf;
 325
 326   for(i=0;i<=(P->Dimension+1);i++)
 327     value_clear(p[i]);
 328   free(p);
 329   value_clear(tmp);
 330   Free_ParamNames(params, C->Dimension);
 331   Domain_Free(P);
 332   Domain_Free(C);
 333   if (CS)
 334     Domain_Free(CS);
 335   return result;
 336 } /* main */
 337
 338
 339
 340