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 struct check_poly_EP_data {
  34     struct check_poly_data   cp;
  35     Polyhedron              *S;
  36     const evalue            *EP;
  37     int                      exist;
  38 };
  39
  40 static int cp_EP(const struct check_poly_data *data, int nparam, Value *z,
  41                  const struct verify_options *options)
  42 {
  43     int k;
  44     int ok;
  45     Value c, tmp;
  46     int pa = options->barvinok->polynomial_approximation;
  47     struct check_poly_EP_data* EP_data = (struct check_poly_EP_data*) data;
  48     const evalue *EP = EP_data->EP;
  49     int exist = EP_data->exist;
  50     Polyhedron *S = EP_data->S;
  51
  52     value_init(c);
  53     value_init(tmp);
  54
  55     /* Computes the ehrhart polynomial */
  56     value_set_double(c, compute_evalue(EP, z)+.25);
  57
  58     if (options->print_all) {
  59         printf("EP(");
  60         value_print(stdout, VALUE_FMT, z[0]);
  61         for (k = 1; k < nparam; ++k) {
  62             printf(", ");
  63             value_print(stdout, VALUE_FMT, z[k]);
  64         }
  65         printf(") = ");
  66         value_print(stdout, VALUE_FMT, c);
  67     }
  68
  69     /* Manually count the number of points */
  70     if (exist)
  71         count_points_e(1, S, exist, nparam, data->z, &tmp);
  72     else
  73         count_points(1, S, data->z, &tmp);
  74
  75     if (options->print_all) {
  76         printf(", count = ");
  77         value_print(stdout, VALUE_FMT, tmp);
  78         printf(". ");
  79     }
  80
  81     if (pa == BV_POLAPPROX_PRE_APPROX)
  82         /* just accept everything */
  83         ok = 1;
  84     else if (pa == BV_POLAPPROX_PRE_LOWER || pa == BV_POLAPPROX_LOWER)
  85         ok = value_le(c, tmp);
  86     else if (pa == BV_POLAPPROX_PRE_UPPER || pa == BV_POLAPPROX_UPPER)
  87         ok = value_ge(c, tmp);
  88     else
  89         ok = value_eq(c, tmp);
  90
  91     if (!ok) {
  92         printf("\n");
  93         fflush(stdout);
  94         fprintf(stderr, "Error !\n");
  95         fprintf(stderr, "EP(");
  96         value_print(stderr, VALUE_FMT, z[0]);
  97         for (k = 1; k < nparam; ++k) {
  98             fprintf(stderr,", ");
  99             value_print(stderr, VALUE_FMT, z[k]);
 100         }
 101         fprintf(stderr, ") should be ");
 102         value_print(stderr, VALUE_FMT, tmp);
 103         fprintf(stderr, ", while EP eval gives ");
 104         value_print(stderr, VALUE_FMT, c);
 105         fprintf(stderr, ".\n");
 106         print_evalue(stderr, EP, options->params);
 107         if (value_zero_p(EP->d) && EP->x.p->type == partition)
 108             for (k = 0; k < EP->x.p->size/2; ++k) {
 109                 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*k]);
 110                 if (in_domain(D, z)) {
 111                     Print_Domain(stderr, D, options->params);
 112                     print_evalue(stderr, &EP->x.p->arr[2*k+1], options->params);
 113                 }
 114         }
 115     } else if (options->print_all)
 116         printf("OK.\n");
 117
 118     value_clear(c);
 119     value_clear(tmp);
 120
 121     return ok;
 122 }
 123
 124 int check_poly_EP(Polyhedron *S, Polyhedron *CS, evalue *EP, int exist,
 125                int nparam, int pos, Value *z, const struct verify_options *options)
 126 {
 127     struct check_poly_EP_data data;
 128     data.cp.z = z;
 129     data.cp.check = cp_EP;
 130     data.S = S;
 131     data.EP = EP;
 132     data.exist = exist;
 133     return check_poly(CS, &data.cp, nparam, pos, z+S->Dimension-nparam+1, options);
 134 }