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;
  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
  48     /* Computes the ehrhart polynomial */
  49     if (!options->exact) {
  50         double d = compute_evalue(EP, z);
  51         if (pa == BV_APPROX_SIGN_LOWER)
  52             d = ceil(d-0.1);
  53         else if (pa == BV_APPROX_SIGN_UPPER)
  54             d = floor(d+0.1);
  55         value_set_double(c, d+.25);
  56     } else {
  57         evalue *res = evalue_eval(EP, z);
  58         if (pa == BV_APPROX_SIGN_LOWER)
  59             mpz_cdiv_q(c, res->x.n, res->d);
  60         else if (pa == BV_APPROX_SIGN_UPPER)
  61             mpz_fdiv_q(c, res->x.n, res->d);
  62         else
  63             mpz_tdiv_q(c, res->x.n, res->d);
  64         free_evalue_refs(res);
  65         free(res);
  66     }
  67
  68     if (options->print_all) {
  69         printf("EP(");
  70         value_print(stdout, VALUE_FMT, z[0]);
  71         for (k = 1; k < nparam; ++k) {
  72             printf(", ");
  73             value_print(stdout, VALUE_FMT, z[k]);
  74         }
  75         printf(") = ");
  76         value_print(stdout, VALUE_FMT, c);
  77     }
  78
  79     /* Manually count the number of points */
  80     if (exist)
  81         count_points_e(1, S, exist, nparam, data->z, &tmp);
  82     else
  83         count_points(1, S, data->z, &tmp);
  84
  85     if (options->print_all) {
  86         printf(", count = ");
  87         value_print(stdout, VALUE_FMT, tmp);
  88         printf(". ");
  89     }
  90
  91     if (pa == BV_APPROX_SIGN_APPROX)
  92         /* just accept everything */
  93         ok = 1;
  94     else if (pa == BV_APPROX_SIGN_LOWER)
  95         ok = value_le(c, tmp);
  96     else if (pa == BV_APPROX_SIGN_UPPER)
  97         ok = value_ge(c, tmp);
  98     else
  99         ok = value_eq(c, tmp);
 100
 101     if (!ok) {
 102         printf("\n");
 103         fflush(stdout);
 104         fprintf(stderr, "Error !\n");
 105         fprintf(stderr, "EP(");
 106         value_print(stderr, VALUE_FMT, z[0]);
 107         for (k = 1; k < nparam; ++k) {
 108             fprintf(stderr,", ");
 109             value_print(stderr, VALUE_FMT, z[k]);
 110         }
 111         fprintf(stderr, ") should be ");
 112         value_print(stderr, VALUE_FMT, tmp);
 113         fprintf(stderr, ", while EP eval gives ");
 114         value_print(stderr, VALUE_FMT, c);
 115         fprintf(stderr, ".\n");
 116         print_evalue(stderr, EP, options->params);
 117         if (value_zero_p(EP->d) && EP->x.p->type == partition)
 118             for (k = 0; k < EP->x.p->size/2; ++k) {
 119                 Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*k]);
 120                 if (in_domain(D, z)) {
 121                     Print_Domain(stderr, D, options->params);
 122                     print_evalue(stderr, &EP->x.p->arr[2*k+1], options->params);
 123                 }
 124         }
 125     } else if (options->print_all)
 126         printf("OK.\n");
 127
 128     value_clear(c);
 129     value_clear(tmp);
 130
 131     return ok;
 132 }
 133
 134 int check_poly_EP(Polyhedron *S, Polyhedron *CS, evalue *EP, int exist,
 135                int nparam, int pos, Value *z, const struct verify_options *options)
 136 {
 137     struct check_poly_EP_data data;
 138     data.cp.z = z;
 139     data.cp.check = cp_EP;
 140     data.S = S;
 141     data.EP = EP;
 142     data.exist = exist;
 143     return check_poly(CS, &data.cp, nparam, pos, z+S->Dimension-nparam+1, options);
 144 }