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