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 <barvinok/evalue.h>
  21 #include <barvinok/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 char **params;
  32
  33 #ifndef PRINT_ALL_RESULTS
  34 int st;
  35 #endif
  36
  37 /****************************************************/
  38 /* function check_poly :                            */
  39 /* scans the parameter space from Min to Max (all   */
  40 /* directions). Computes the number of points in    */
  41 /* the polytope using both methods, and compare them*/
  42 /* returns 1 on success                             */
  43 /****************************************************/
  44
  45 int check_poly(Polyhedron *S,Polyhedron *CS,Enumeration *en,
  46                int nparam,int pos,Value *z) {
  47
  48   int k;
  49   Value c,tmp,*ctmp;
  50   Value LB, UB;
  51
  52   value_init(c); value_init(tmp);
  53   value_init(LB);
  54   value_init(UB);
  55
  56   if(pos == nparam) {
  57
  58     /* Computes the ehrhart polynomial */
  59     value_assign(c,*(ctmp=compute_poly(en,&z[S->Dimension-nparam+1])));
  60     value_clear(*ctmp);
  61     free(ctmp);
  62     /* if c=0 we may be out of context. */
  63     /* scanning is useless in this case*/
  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
  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         {
 101                  Enumeration *ee;
 102                  Enumeration_Print(stderr, en, params);
 103                  ee = en;
 104                  while (ee) {
 105                 if (in_domain(ee->ValidityDomain,&z[S->Dimension-nparam+1])) {
 106                          Print_Domain(stderr, ee->ValidityDomain, params);
 107                          print_evalue(stderr, &ee->EP, params);
 108                 }
 109                 ee = ee->next;
 110           }
 111         }
 112 #ifndef DONT_BREAK_ON_ERROR
 113         value_clear(c); value_clear(tmp);
 114         value_clear(LB);
 115         value_clear(UB);
 116         return(0);
 117 #endif
 118       }
 119
 120 #ifdef PRINT_ALL_RESULTS
 121       else
 122         printf("OK.\n");
 123 #endif
 124   }
 125   else {
 126     int ok =
 127         !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
 128     assert(ok);
 129     for(value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
 130
 131 #ifndef PRINT_ALL_RESULTS
 132       k = VALUE_TO_INT(tmp);
 133       if(!pos && !(k%st)) {
 134         printf("o");
 135         fflush(stdout);
 136       }
 137 #endif
 138
 139       value_assign(z[pos+S->Dimension-nparam+1],tmp);
 140       if(!check_poly(S, CS->next, en, nparam, pos+1, z)) {
 141         value_clear(c); value_clear(tmp);
 142         value_clear(LB);
 143         value_clear(UB);
 144         return(0);
 145       }
 146     }
 147     value_set_si(z[pos+S->Dimension-nparam+1],0);
 148   }
 149   value_clear(c); value_clear(tmp);
 150   value_clear(LB);
 151   value_clear(UB);
 152   return(1);
 153 } /* check_poly */
 154