verify_main.cc

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