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