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 struct check_poly_EP_data {
  27     struct check_poly_data   cp;
  28     Polyhedron              *S;
  29     const evalue            *EP;
  30     int                      exist;
  31 };
  32
  33 static int cp_EP(const struct check_poly_data *data, int nparam, Value *z,
  34                  const struct verify_options *options)
  35 {
  36     int k;
  37     int ok;
  38     Value c, tmp, one;
  39     int pa = options->barvinok->polynomial_approximation;
  40     struct check_poly_EP_data* EP_data = (struct check_poly_EP_data*) data;
  41     const evalue *EP = EP_data->EP;
  42     int exist = EP_data->exist;
  43     Polyhedron *S = EP_data->S;
  44
  45     value_init(c);
  46     value_init(tmp);
  47     value_init(one);
  48     value_set_si(one, 1);
  49
  50     /* Computes the ehrhart polynomial */
  51     if (!options->exact) {
  52         double d = compute_evalue(EP, z);
  53         if (pa == BV_APPROX_SIGN_LOWER)
  54             d = ceil(d-0.1);
  55         else if (pa == BV_APPROX_SIGN_UPPER)
  56             d = floor(d+0.1);
  57         value_set_double(c, d+.25);
  58     } else {
  59         evalue *res = evalue_eval(EP, z);
  60         if (pa == BV_APPROX_SIGN_LOWER)
  61             mpz_cdiv_q(c, res->x.n, res->d);
  62         else if (pa == BV_APPROX_SIGN_UPPER)
  63             mpz_fdiv_q(c, res->x.n, res->d);
  64         else
  65             mpz_tdiv_q(c, res->x.n, res->d);
  66         evalue_free(res);
  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 (pa == BV_APPROX_SIGN_APPROX)
  76         /* just accept everything */
  77         ok = 1;
  78     else if (pa == BV_APPROX_SIGN_LOWER)
  79         ok = value_le(c, tmp);
  80     else if (pa == BV_APPROX_SIGN_UPPER)
  81         ok = value_ge(c, tmp);
  82     else
  83         ok = value_eq(c, tmp);
  84
  85     check_poly_print(ok, nparam, z, tmp, one, c, one,
  86                      "EP", "count", "EP eval", options);
  87
  88     if (!ok) {
  89         print_evalue(stderr, EP, options->params);
  90         if (value_zero_p(EP->d) && EP->x.p->type == partition)
  91             for (k = 0; k < EP->x.p->size/2; ++k) {
  92                 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*k]);
  93                 if (in_domain(D, z)) {
  94                     Print_Domain(stderr, D, options->params);
  95                     print_evalue(stderr, &EP->x.p->arr[2*k+1], options->params);
  96                 }
  97         }
  98     }
  99
 100     value_clear(c);
 101     value_clear(tmp);
 102     value_clear(one);
 103
 104     return ok;
 105 }
 106
 107 int check_poly_EP(Polyhedron *S, Polyhedron *CS, evalue *EP, int exist,
 108                int nparam, int pos, Value *z, const struct verify_options *options)
 109 {
 110     struct check_poly_EP_data data;
 111     data.cp.z = z;
 112     data.cp.check = cp_EP;
 113     data.S = S;
 114     data.EP = EP;
 115     data.exist = exist;
 116     return check_poly(CS, &data.cp, nparam, pos, z+S->Dimension-nparam+1, options);
 117 }