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