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 <barvinok/evalue.h>
  20 #include <barvinok/barvinok.h>
  21 #include "verif_ehrhart.h"
  22
  23 #undef CS   /* for Solaris 10 */
  24
  25 #include "config.h"
  26 #ifndef HAVE_COUNT_POINTS4
  27 #define count_points(a,b,c,d) {                         \
  28             int cc = count_points(a,b,c);               \
  29             value_set_si(*d,cc);                        \
  30         }
  31 #endif
  32
  33 char **params;
  34
  35 int st;
  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, int print_all)
  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     if (print_all) {
  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     }
  76
  77       /* Manually count the number of points */
  78       count_points(1,S,z,&tmp);
  79     if (print_all) {
  80         printf(", count = ");
  81         value_print(stdout, P_VALUE_FMT, tmp);
  82         printf(". ");
  83     }
  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       } else if (print_all)
 119         printf("OK.\n");
 120   }
 121   else {
 122     int ok =
 123         !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
 124     assert(ok);
 125     for(value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
 126
 127     if (!print_all) {
 128       k = VALUE_TO_INT(tmp);
 129       if(!pos && !(k%st)) {
 130         printf("o");
 131         fflush(stdout);
 132       }
 133     }
 134
 135       value_assign(z[pos+S->Dimension-nparam+1],tmp);
 136       if (!check_poly(S, CS->next, en, nparam, pos+1, z, print_all)) {
 137         value_clear(c); value_clear(tmp);
 138         value_clear(LB);
 139         value_clear(UB);
 140         return(0);
 141       }
 142     }
 143     value_set_si(z[pos+S->Dimension-nparam+1],0);
 144   }
 145   value_clear(c); value_clear(tmp);
 146   value_clear(LB);
 147   value_clear(UB);
 148   return(1);
 149 } /* check_poly */
 150