verif_ehrhart.c

   1 /*************************************************/
   2 /*     verif_ehrhart.c                           */
   3 /* program to compare effective number of points */
   4 /* in a polytope with the corresponding          */
   5 /* evaluation of the Ehrhart polynomial.         */
   6 /* Parameters vary in range -RANGE to RANGE      */
   7 /* (define below) by default.                    */
   8 /* Can be overridden by specifying               */
   9 /* -r<RANGE>, or -m<min> and -M<max>             */
  10 /*                                               */
  11 /* written by Vincent Loechner (c) 2000.         */
  12 /*  loechner@icps.u-strasbg.fr                   */
  13 /*************************************************/
  14
  15 #include <stdio.h>
  16 #include <string.h>
  17 #include <stdlib.h>
  18
  19 #include <polylib/polylibgmp.h>
  20 #include "ev_operations.h"
  21 #include <barvinok.h>
  22
  23 #include "config.h"
  24 #ifndef HAVE_COUNT_POINTS4
  25 #define count_points(a,b,c,d) {                         \
  26             int cc = count_points(a,b,c);               \
  27             value_set_si(*d,cc);                        \
  28         }
  29 #endif
  30
  31 Value Min, Max;
  32
  33 char **params;
  34
  35 #ifndef PRINT_ALL_RESULTS
  36 int st;
  37 #endif
  38
  39 /****************************************************/
  40 /* function check_poly :                            */
  41 /* scans the parameter space from Min to Max (all   */
  42 /* directions). Computes the number of points in    */
  43 /* the polytope using both methods, and compare them*/
  44 /* returns 1 on success                             */
  45 /****************************************************/
  46
  47 int check_poly(Polyhedron *S,Polyhedron *C,Enumeration *en,
  48                int nparam,int pos,Value *z) {
  49
  50   int k;
  51   Value c,tmp,*ctmp;
  52
  53   value_init(c); value_init(tmp);
  54
  55   if(pos == nparam) {
  56
  57     /* Computes the ehrhart polynomial */
  58     value_assign(c,*(ctmp=compute_poly(en,&z[S->Dimension-nparam+1])));
  59     free(ctmp);
  60     /* if c=0 we may be out of context. */
  61     /* scanning is useless in this case*/
  62     if(!in_domain(C,&z[S->Dimension-nparam+1])) {
  63
  64       /* ok */ ;
  65     }
  66     else {
  67
  68 #ifdef PRINT_ALL_RESULTS
  69       printf("EP( ");
  70       value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
  71       for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  72         printf(", ");
  73         value_print(stdout,VALUE_FMT,z[k]);
  74       }
  75       printf(" ) = ");
  76       value_print(stdout,VALUE_FMT,c);
  77       printf(" ");
  78 #endif
  79
  80       /* Manually count the number of points */
  81       count_points(1,S,z,&tmp);
  82 #ifdef PRINT_ALL_RESULTS
  83         printf(", count = ");
  84         value_print(stdout, P_VALUE_FMT, tmp);
  85         printf(". ");
  86 #endif
  87
  88       if(value_ne(tmp,c)) {
  89         printf("\n");
  90         fflush(stdout);
  91         fprintf(stderr,"Error !\n");
  92         fprintf(stderr,"EP( ");
  93         value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
  94         for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  95           fprintf(stderr,", ");
  96           value_print(stderr,VALUE_FMT,z[k]);
  97         }
  98         fprintf(stderr," ) should be ");
  99         value_print(stderr,VALUE_FMT,tmp);
 100         fprintf(stderr,", while EP eval gives ");
 101         value_print(stderr,VALUE_FMT,c);
 102         fprintf(stderr,".\n");
 103         {
 104                  Enumeration *ee;
 105                  Enumeration_Print(stderr, en, params);
 106                  ee = en;
 107                  while (ee) {
 108                 if (in_domain(ee->ValidityDomain,&z[S->Dimension-nparam+1])) {
 109                          Print_Domain(stderr, ee->ValidityDomain, params);
 110                          print_evalue(stderr, &ee->EP, params);
 111                 }
 112                 ee = ee->next;
 113           }
 114         }
 115 #ifndef DONT_BREAK_ON_ERROR
 116         value_clear(c); value_clear(tmp);
 117         return(0);
 118 #endif
 119       }
 120
 121 #ifdef PRINT_ALL_RESULTS
 122       else
 123         printf("OK.\n");
 124 #endif
 125     }
 126   }
 127   else
 128     for(value_assign(tmp,Min); value_le(tmp,Max); value_increment(tmp,tmp)) {
 129
 130 #ifndef PRINT_ALL_RESULTS
 131       k = VALUE_TO_INT(tmp);
 132       if(!pos && !(k%st)) {
 133         printf("o");
 134         fflush(stdout);
 135       }
 136 #endif
 137
 138       value_assign(z[pos+S->Dimension-nparam+1],tmp);
 139       if(!check_poly(S,C,en,nparam,pos+1,z)) {
 140         value_clear(c); value_clear(tmp);
 141         return(0);
 142       }
 143     }
 144   value_clear(c); value_clear(tmp);
 145   return(1);
 146 } /* check_poly */
 147