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