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 #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 int st;
  34
  35 /****************************************************/
  36 /* function check_poly :                            */
  37 /* scans the parameter space from Min to Max (all   */
  38 /* directions). Computes the number of points in    */
  39 /* the polytope using both methods, and compare them*/
  40 /* returns 1 on success                             */
  41 /****************************************************/
  42
  43 int check_poly(Polyhedron *S,Polyhedron *CS,Enumeration *en,
  44                int nparam, int pos, Value *z, int print_all)
  45 {
  46   int k;
  47   Value c,tmp,*ctmp;
  48   Value LB, UB;
  49
  50   value_init(c); value_init(tmp);
  51   value_init(LB);
  52   value_init(UB);
  53
  54   if(pos == nparam) {
  55
  56     /* Computes the ehrhart polynomial */
  57     value_assign(c,*(ctmp=compute_poly(en,&z[S->Dimension-nparam+1])));
  58     value_clear(*ctmp);
  59     free(ctmp);
  60     /* if c=0 we may be out of context. */
  61     /* scanning is useless in this case*/
  62
  63     if (print_all) {
  64       printf("EP( ");
  65       value_print(stdout,VALUE_FMT,z[S->Dimension-nparam+1]);
  66       for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  67         printf(", ");
  68         value_print(stdout,VALUE_FMT,z[k]);
  69       }
  70       printf(" ) = ");
  71       value_print(stdout,VALUE_FMT,c);
  72       printf(" ");
  73     }
  74
  75       /* Manually count the number of points */
  76       count_points(1,S,z,&tmp);
  77     if (print_all) {
  78         printf(", count = ");
  79         value_print(stdout, P_VALUE_FMT, tmp);
  80         printf(". ");
  81     }
  82
  83       if(value_ne(tmp,c)) {
  84         printf("\n");
  85         fflush(stdout);
  86         fprintf(stderr,"Error !\n");
  87         fprintf(stderr,"EP( ");
  88         value_print(stderr,VALUE_FMT,z[S->Dimension-nparam+1]);
  89         for(k=S->Dimension-nparam+2;k<=S->Dimension;++k) {
  90           fprintf(stderr,", ");
  91           value_print(stderr,VALUE_FMT,z[k]);
  92         }
  93         fprintf(stderr," ) should be ");
  94         value_print(stderr,VALUE_FMT,tmp);
  95         fprintf(stderr,", while EP eval gives ");
  96         value_print(stderr,VALUE_FMT,c);
  97         fprintf(stderr,".\n");
  98         {
  99                  Enumeration *ee;
 100                  Enumeration_Print(stderr, en, params);
 101                  ee = en;
 102                  while (ee) {
 103                 if (in_domain(ee->ValidityDomain,&z[S->Dimension-nparam+1])) {
 104                          Print_Domain(stderr, ee->ValidityDomain, params);
 105                          print_evalue(stderr, &ee->EP, params);
 106                 }
 107                 ee = ee->next;
 108           }
 109         }
 110 #ifndef DONT_BREAK_ON_ERROR
 111         value_clear(c); value_clear(tmp);
 112         value_clear(LB);
 113         value_clear(UB);
 114         return(0);
 115 #endif
 116       } else if (print_all)
 117         printf("OK.\n");
 118   }
 119   else {
 120     int ok =
 121         !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
 122     assert(ok);
 123     for(value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
 124
 125     if (!print_all) {
 126       k = VALUE_TO_INT(tmp);
 127       if(!pos && !(k%st)) {
 128         printf("o");
 129         fflush(stdout);
 130       }
 131     }
 132
 133       value_assign(z[pos+S->Dimension-nparam+1],tmp);
 134       if (!check_poly(S, CS->next, en, nparam, pos+1, z, print_all)) {
 135         value_clear(c); value_clear(tmp);
 136         value_clear(LB);
 137         value_clear(UB);
 138         return(0);
 139       }
 140     }
 141     value_set_si(z[pos+S->Dimension-nparam+1],0);
 142   }
 143   value_clear(c); value_clear(tmp);
 144   value_clear(LB);
 145   value_clear(UB);
 146   return(1);
 147 } /* check_poly */
 148