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
cp_EP(const struct check_poly_data * data,int nparam,Value * z,const struct verify_options * options)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
check_poly_EP(Polyhedron * S,Polyhedron * CS,evalue * EP,int exist,int nparam,int pos,Value * z,const struct verify_options * options)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 }
116