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 #include <math.h>
  19
  20 #include <barvinok/evalue.h>
  21 #include <barvinok/barvinok.h>
  22 #include "verif_ehrhart.h"
  23
  24 #undef CS   /* for Solaris 10 */
  25
  26 #include "config.h"
  27 #ifndef HAVE_COUNT_POINTS4
  28 #define count_points(a,b,c,d) {                         \
  29             int cc = count_points(a,b,c);               \
  30             value_set_si(*d,cc);                        \
  31         }
  32 #endif
  33
  34 struct check_poly_EP_data {
  35     struct check_poly_data   cp;
  36     Polyhedron              *S;
  37     const evalue            *EP;
  38     int                      exist;
  39 };
  40
  41 static int cp_EP(const struct check_poly_data *data, int nparam, Value *z,
  42                  const struct verify_options *options)
  43 {
  44     int k;
  45     int ok;
  46     Value c, tmp;
  47     int pa = options->barvinok->polynomial_approximation;
  48     struct check_poly_EP_data* EP_data = (struct check_poly_EP_data*) data;
  49     const evalue *EP = EP_data->EP;
  50     int exist = EP_data->exist;
  51     Polyhedron *S = EP_data->S;
  52
  53     value_init(c);
  54     value_init(tmp);
  55
  56     /* Computes the ehrhart polynomial */
  57     if (!options->exact) {
  58         double d = compute_evalue(EP, z);
  59         if (pa == BV_APPROX_SIGN_LOWER)
  60             d = ceil(d-0.1);
  61         else if (pa == BV_APPROX_SIGN_UPPER)
  62             d = floor(d+0.1);
  63         value_set_double(c, d+.25);
  64     } else {
  65         evalue *res = evalue_eval(EP, z);
  66         if (pa == BV_APPROX_SIGN_LOWER)
  67             mpz_cdiv_q(c, res->x.n, res->d);
  68         else if (pa == BV_APPROX_SIGN_UPPER)
  69             mpz_fdiv_q(c, res->x.n, res->d);
  70         else
  71             mpz_tdiv_q(c, res->x.n, res->d);
  72         free_evalue_refs(res);
  73         free(res);
  74     }
  75
  76     if (options->print_all) {
  77         printf("EP(");
  78         value_print(stdout, VALUE_FMT, z[0]);
  79         for (k = 1; k < nparam; ++k) {
  80             printf(", ");
  81             value_print(stdout, VALUE_FMT, z[k]);
  82         }
  83         printf(") = ");
  84         value_print(stdout, VALUE_FMT, c);
  85     }
  86
  87     /* Manually count the number of points */
  88     if (exist)
  89         count_points_e(1, S, exist, nparam, data->z, &tmp);
  90     else
  91         count_points(1, S, data->z, &tmp);
  92
  93     if (options->print_all) {
  94         printf(", count = ");
  95         value_print(stdout, VALUE_FMT, tmp);
  96         printf(". ");
  97     }
  98
  99     if (pa == BV_APPROX_SIGN_APPROX)
 100         /* just accept everything */
 101         ok = 1;
 102     else if (pa == BV_APPROX_SIGN_LOWER)
 103         ok = value_le(c, tmp);
 104     else if (pa == BV_APPROX_SIGN_UPPER)
 105         ok = value_ge(c, tmp);
 106     else
 107         ok = value_eq(c, tmp);
 108
 109     if (!ok) {
 110         printf("\n");
 111         fflush(stdout);
 112         fprintf(stderr, "Error !\n");
 113         fprintf(stderr, "EP(");
 114         value_print(stderr, VALUE_FMT, z[0]);
 115         for (k = 1; k < nparam; ++k) {
 116             fprintf(stderr,", ");
 117             value_print(stderr, VALUE_FMT, z[k]);
 118         }
 119         fprintf(stderr, ") should be ");
 120         value_print(stderr, VALUE_FMT, tmp);
 121         fprintf(stderr, ", while EP eval gives ");
 122         value_print(stderr, VALUE_FMT, c);
 123         fprintf(stderr, ".\n");
 124         print_evalue(stderr, EP, options->params);
 125         if (value_zero_p(EP->d) && EP->x.p->type == partition)
 126             for (k = 0; k < EP->x.p->size/2; ++k) {
 127                 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*k]);
 128                 if (in_domain(D, z)) {
 129                     Print_Domain(stderr, D, options->params);
 130                     print_evalue(stderr, &EP->x.p->arr[2*k+1], options->params);
 131                 }
 132         }
 133     } else if (options->print_all)
 134         printf("OK.\n");
 135
 136     value_clear(c);
 137     value_clear(tmp);
 138
 139     return ok;
 140 }
 141
 142 int check_poly_EP(Polyhedron *S, Polyhedron *CS, evalue *EP, int exist,
 143                int nparam, int pos, Value *z, const struct verify_options *options)
 144 {
 145     struct check_poly_EP_data data;
 146     data.cp.z = z;
 147     data.cp.check = cp_EP;
 148     data.S = S;
 149     data.EP = EP;
 150     data.exist = exist;
 151     return check_poly(CS, &data.cp, nparam, pos, z+S->Dimension-nparam+1, options);
 152 }