verify_main.cc

   1 #include <gmp.h>
   2 #include <NTL/mat_ZZ.h>
   3 extern "C" {
   4 #include <polylib/polylibgmp.h>
   5 #include <barvinok/evalue.h>
   6 }
   7 #include <barvinok/util.h>
   8 #include <barvinok/barvinok.h>
   9
  10 #include "verif_ehrhart.h"
  11
  12 #ifdef HAVE_GROWING_CHERNIKOVA
  13 #define MAXRAYS    0
  14 #else
  15 #define MAXRAYS  600
  16 #endif
  17
  18 #include "config.h"
  19 #ifndef HAVE_GETOPT_H
  20 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
  21 #else
  22 #include <getopt.h>
  23 struct option options[] = {
  24     { "explicit",  no_argument,  0,  'e' },
  25     { "series",  no_argument,  0,  's' },
  26     { "verbose",  no_argument,  0,  'v' },
  27     { "version",   no_argument,  0,  'V' },
  28     { 0, 0, 0, 0 }
  29 };
  30 #endif
  31
  32 /* RANGE : normal range for evalutations (-RANGE -> RANGE) */
  33 #define RANGE 50
  34
  35 /* SRANGE : small range for evalutations */
  36 #define SRANGE 15
  37
  38 /* if dimension >= BIDDIM, use SRANGE */
  39 #define BIGDIM 5
  40
  41 /* VSRANGE : very small range for evalutations */
  42 #define VSRANGE 5
  43
  44 /* if dimension >= VBIDDIM, use VSRANGE */
  45 #define VBIGDIM 8
  46
  47 int check_series(Polyhedron *S, Polyhedron *CS, gen_fun *gf,
  48                  int nparam,int pos,Value *z)
  49 {
  50     int k;
  51     Value c, tmp;
  52     Value LB, UB;
  53
  54     value_init(c);
  55     value_init(tmp);
  56     value_init(LB);
  57     value_init(UB);
  58
  59     if(pos == nparam) {
  60         /* Computes the coefficient */
  61         gf->coefficient(&z[S->Dimension-nparam+1], &c);
  62
  63         /* if c=0 we may be out of context. */
  64         /* scanning is useless in this case*/
  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     } else {
 110         int ok =
 111           !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
 112         assert(ok);
 113         for(value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
 114 #ifndef PRINT_ALL_RESULTS
 115           k = VALUE_TO_INT(tmp);
 116           if(!pos && !(k%st)) {
 117             printf("o");
 118             fflush(stdout);
 119           }
 120 #endif
 121           value_assign(z[pos+S->Dimension-nparam+1],tmp);
 122           if(!check_series(S, CS->next, gf, nparam, pos+1, z)) {
 123             value_clear(c); value_clear(tmp);
 124             value_clear(LB);
 125             value_clear(UB);
 126             return(0);
 127           }
 128         }
 129         value_set_si(z[pos+S->Dimension-nparam+1],0);
 130      }
 131
 132     value_clear(c);
 133     value_clear(tmp);
 134     value_clear(LB);
 135     value_clear(UB);
 136     return 1;
 137 }
 138
 139 int main(int argc,char *argv[]) {
 140
 141   Matrix *C1, *P1, *MM;
 142   Polyhedron *C, *P, *S, *CS, *U;
 143   Polyhedron *CC, *PP;
 144   Enumeration *en = NULL;
 145   Value *p, tmp;
 146   Value Min, Max;
 147   int i,j;
 148     int m = INT_MAX, M = INT_MIN, r;
 149     int c, ind = 0;
 150     int series = 0;
 151     int verbose = 0;
 152     int function = 0;
 153     int result = 0;
 154
 155     while ((c = getopt_long(argc, argv, "m:M:r:sveV", options, &ind)) != -1) {
 156         switch (c) {
 157         case 'e':
 158             function = 1;
 159             break;
 160         case 's':
 161             series = 1;
 162             break;
 163         case 'v':
 164             verbose = 1;
 165             break;
 166         case 'm':
 167             m = atoi(optarg);
 168             break;
 169         case 'M':
 170             M = atoi(optarg);
 171             break;
 172         case 'r':
 173             M = atoi(optarg);
 174             m = -M;
 175             break;
 176         case 'V':
 177             printf(barvinok_version());
 178             exit(0);
 179             break;
 180         }
 181     }
 182
 183 /******* Read the input *********/
 184   P1 = Matrix_Read();
 185   C1 = Matrix_Read();
 186
 187   if(C1->NbColumns < 2) {
 188     fprintf(stderr,"Not enough parameters !\n");
 189     exit(0);
 190   }
 191
 192   P = Constraints2Polyhedron(P1,MAXRAYS);
 193   C = Constraints2Polyhedron(C1,MAXRAYS);
 194   params = Read_ParamNames(stdin, C->Dimension);
 195   Matrix_Free(C1);
 196   Matrix_Free(P1);
 197
 198   /******* Read the options: initialize Min and Max ********/
 199   if(P->Dimension >= VBIGDIM)
 200     r = VSRANGE;
 201   else if(P->Dimension >= BIGDIM)
 202     r = SRANGE;
 203   else
 204     r = RANGE;
 205   if (M == INT_MIN)
 206     M = r;
 207   if (m == INT_MAX)
 208     m = -r;
 209
 210
 211   if(m > M) {
 212     fprintf(stderr,"Nothing to do: Min > Max !\n");
 213     return(0);
 214   }
 215   value_init(Min);
 216   value_init(Max);
 217   value_set_si(Min,m);
 218   value_set_si(Max,M);
 219   value_init(tmp);
 220
 221   /******* Compute true context *******/
 222   CC = align_context(C,P->Dimension,MAXRAYS);
 223   PP = DomainIntersection(P,CC,MAXRAYS);
 224   Domain_Free(CC);
 225   C1 = Matrix_Alloc(C->Dimension+1,P->Dimension+1);
 226
 227   for(i=0;i<C1->NbRows;i++)
 228     for(j=0;j<C1->NbColumns;j++)
 229       if(i==j-P->Dimension+C->Dimension)
 230         value_set_si(C1->p[i][j],1);
 231       else
 232         value_set_si(C1->p[i][j],0);
 233   CC = Polyhedron_Image(PP,C1,MAXRAYS);
 234   Matrix_Free(C1);
 235   Domain_Free(PP);
 236   Domain_Free(C);
 237   C = CC;
 238
 239     /* Intersect context with range */
 240     if(C->Dimension > 0) {
 241         MM = Matrix_Alloc(2*C->Dimension, C->Dimension+2);
 242         for (i = 0; i < C->Dimension; ++i) {
 243             value_set_si(MM->p[2*i][0], 1);
 244             value_set_si(MM->p[2*i][1+i], 1);
 245             value_set_si(MM->p[2*i][1+C->Dimension], -m);
 246             value_set_si(MM->p[2*i+1][0], 1);
 247             value_set_si(MM->p[2*i+1][1+i], -1);
 248             value_set_si(MM->p[2*i+1][1+C->Dimension], M);
 249         }
 250         CC = AddConstraints(MM->p[0], 2*C->Dimension, C, MAXRAYS);
 251         U = Universe_Polyhedron(0);
 252         CS = Polyhedron_Scan(CC, U, MAXRAYS);
 253         Polyhedron_Free(U);
 254         Polyhedron_Free(CC);
 255         Matrix_Free(MM);
 256     } else
 257         CS = NULL;
 258
 259     gen_fun *gf = 0;
 260
 261     /******* Compute EP *********/
 262     if (!series) {
 263         en = barvinok_enumerate(P,C,MAXRAYS);
 264         if (verbose)
 265             Enumeration_Print(stdout, en, params);
 266     } else {
 267         evalue *EP;
 268         gf = barvinok_series(P, C, MAXRAYS);
 269         if (verbose) {
 270             gf->print(std::cout, C->Dimension, params);
 271             puts("");
 272         }
 273         if (function) {
 274             EP = *gf;
 275             en =  partition2enumeration(EP);
 276             if (verbose)
 277                 Enumeration_Print(stdout, en, params);
 278         }
 279     }
 280
 281   /******* Initializations for check *********/
 282   p = (Value *)malloc(sizeof(Value) * (P->Dimension+2));
 283   for(i=0;i<=P->Dimension;i++) {
 284     value_init(p[i]);
 285     value_set_si(p[i],0);
 286   }
 287   value_init(p[i]);
 288   value_set_si(p[i],1);
 289
 290   /* S = scanning list of polyhedra */
 291   S = Polyhedron_Scan(P,C,MAXRAYS);
 292
 293 #ifndef PRINT_ALL_RESULTS
 294   if(C->Dimension > 0) {
 295     value_subtract(tmp,Max,Min);
 296     if (VALUE_TO_INT(tmp) > 80)
 297       st = 1+(VALUE_TO_INT(tmp))/80;
 298     else
 299       st=1;
 300     for(i=VALUE_TO_INT(Min);i<=VALUE_TO_INT(Max);i+=st)
 301       printf(".");
 302     printf( "\r" );
 303     fflush(stdout);
 304   }
 305 #endif
 306
 307     /******* CHECK NOW *********/
 308     if(S) {
 309         if (!series || function) {
 310             if (!check_poly(S, CS,en,C->Dimension,0,p))
 311                 result = -1;
 312         } else {
 313             if (!check_series(S, CS,gf,C->Dimension,0,p))
 314                 result = -1;
 315         }
 316         Domain_Free(S);
 317     }
 318
 319     if (result == -1)
 320         fprintf(stderr,"Check failed !\n");
 321
 322 #ifndef PRINT_ALL_RESULTS
 323   printf( "\n" );
 324 #endif
 325
 326     if (en)
 327         Enumeration_Free(en);
 328     if (gf)
 329         delete gf;
 330
 331   for(i=0;i<=(P->Dimension+1);i++)
 332     value_clear(p[i]);
 333   free(p);
 334   value_clear(Min);
 335   value_clear(Max);
 336   value_clear(tmp);
 337   Free_ParamNames(params, C->Dimension);
 338   Domain_Free(P);
 339   Domain_Free(C);
 340   if (CS)
 341     Domain_Free(CS);
 342   return result;
 343 } /* main */
 344
 345
 346
 347