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
  19 #include <barvinok/evalue.h>
  20 #include <barvinok/barvinok.h>
  21 #include "verif_ehrhart.h"
  22
  23 #undef CS   /* for Solaris 10 */
  24
  25 #include "config.h"
  26 #ifndef HAVE_COUNT_POINTS4
  27 #define count_points(a,b,c,d) {                         \
  28             int cc = count_points(a,b,c);               \
  29             value_set_si(*d,cc);                        \
  30         }
  31 #endif
  32
  33 Polyhedron *check_poly_context_scan(Polyhedron *C, struct verify_options *options)
  34 {
  35     int i;
  36     Matrix *MM;
  37     Polyhedron *CC, *CS, *U;
  38
  39     if (C->Dimension <= 0)
  40         return NULL;
  41
  42     /* Intersect context with range */
  43     MM = Matrix_Alloc(2*C->Dimension, C->Dimension+2);
  44     for (i = 0; i < C->Dimension; ++i) {
  45         value_set_si(MM->p[2*i][0], 1);
  46         value_set_si(MM->p[2*i][1+i], 1);
  47         value_set_si(MM->p[2*i][1+C->Dimension], -options->m);
  48         value_set_si(MM->p[2*i+1][0], 1);
  49         value_set_si(MM->p[2*i+1][1+i], -1);
  50         value_set_si(MM->p[2*i+1][1+C->Dimension], options->M);
  51     }
  52     CC = AddConstraints(MM->p[0], 2*C->Dimension, C, options->barvinok->MaxRays);
  53     U = Universe_Polyhedron(0);
  54     CS = Polyhedron_Scan(CC, U, options->barvinok->MaxRays);
  55     Polyhedron_Free(U);
  56     Polyhedron_Free(CC);
  57     Matrix_Free(MM);
  58     return CS;
  59 }
  60
  61 void check_poly_init(Polyhedron *C, struct verify_options *options)
  62 {
  63     int d, i;
  64
  65     if (options->print_all)
  66         return;
  67     if (C->Dimension <= 0)
  68         return;
  69
  70     d = options->M - options->m;
  71     if (d > 80)
  72         options->st = 1+d/80;
  73     else
  74         options->st = 1;
  75     for (i = options->m; i <= options->M; i += options->st)
  76         printf(".");
  77     printf( "\r" );
  78     fflush(stdout);
  79 }
  80
  81 /****************************************************/
  82 /* function check_poly :                            */
  83 /* scans the parameter space from Min to Max (all   */
  84 /* directions). Computes the number of points in    */
  85 /* the polytope using both methods, and compare them*/
  86 /* returns 1 on success                             */
  87 /****************************************************/
  88
  89 int check_poly(Polyhedron *S, Polyhedron *CS, evalue *EP, int exist,
  90                int nparam, int pos, Value *z, const struct verify_options *options)
  91 {
  92     int k;
  93     Value c, tmp;
  94     int ok;
  95     int pa = options->barvinok->polynomial_approximation;
  96
  97     value_init(c);
  98     value_init(tmp);
  99
 100     if (pos == nparam) {
 101         /* Computes the ehrhart polynomial */
 102         value_set_double(c, compute_evalue(EP, &z[S->Dimension-nparam+1])+.25);
 103
 104         if (options->print_all) {
 105             printf("EP(");
 106             value_print(stdout, VALUE_FMT, z[S->Dimension-nparam+1]);
 107             for(k=S->Dimension-nparam+2; k<=S->Dimension; ++k) {
 108                 printf(", ");
 109                 value_print(stdout,VALUE_FMT,z[k]);
 110             }
 111             printf(") = ");
 112             value_print(stdout, VALUE_FMT, c);
 113         }
 114
 115         /* Manually count the number of points */
 116         if (exist)
 117             count_points_e(1, S, exist, nparam, z, &tmp);
 118         else
 119             count_points(1, S, z, &tmp);
 120
 121         if (options->print_all) {
 122             printf(", count = ");
 123             value_print(stdout, VALUE_FMT, tmp);
 124             printf(". ");
 125         }
 126
 127         if (pa == BV_POLAPPROX_PRE_APPROX)
 128             /* just accept everything */
 129             ok = 1;
 130         else if (pa == BV_POLAPPROX_PRE_LOWER || pa == BV_POLAPPROX_LOWER)
 131             ok = value_le(c, tmp);
 132         else if (pa == BV_POLAPPROX_PRE_UPPER || pa == BV_POLAPPROX_UPPER)
 133             ok = value_ge(c, tmp);
 134         else
 135             ok = value_eq(c, tmp);
 136
 137         if (!ok) {
 138             printf("\n");
 139             fflush(stdout);
 140             fprintf(stderr, "Error !\n");
 141             fprintf(stderr, "EP(");
 142             value_print(stderr, VALUE_FMT, z[S->Dimension-nparam+1]);
 143             for(k=S->Dimension-nparam+2; k<=S->Dimension; ++k) {
 144                 fprintf(stderr,", ");
 145                 value_print(stderr,VALUE_FMT,z[k]);
 146             }
 147             fprintf(stderr, ") should be ");
 148             value_print(stderr, VALUE_FMT, tmp);
 149             fprintf(stderr, ", while EP eval gives ");
 150             value_print(stderr, VALUE_FMT, c);
 151             fprintf(stderr, ".\n");
 152             print_evalue(stderr, EP, options->params);
 153             if (value_zero_p(EP->d) && EP->x.p->type == partition)
 154                 for (k = 0; k < EP->x.p->size/2; ++k) {
 155                     Polyhedron *D = EVALUE_DOMAIN(EP->x.p->arr[2*k]);
 156                     if (in_domain(D, &z[S->Dimension-nparam+1])) {
 157                         Print_Domain(stderr, D, options->params);
 158                         print_evalue(stderr, &EP->x.p->arr[2*k+1], options->params);
 159                     }
 160             }
 161             if (!options->continue_on_error) {
 162                 value_clear(c); value_clear(tmp);
 163                 return 0;
 164             }
 165         } else if (options->print_all)
 166             printf("OK.\n");
 167     } else {
 168         Value LB, UB;
 169         int ok;
 170         value_init(LB);
 171         value_init(UB);
 172         ok = !(lower_upper_bounds(1+pos, CS, &z[S->Dimension-nparam], &LB, &UB));
 173         assert(ok);
 174         for(value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
 175             if (!options->print_all) {
 176                 k = VALUE_TO_INT(tmp);
 177                 if (!pos && !(k % options->st)) {
 178                     printf("o");
 179                     fflush(stdout);
 180                 }
 181             }
 182
 183             value_assign(z[pos+S->Dimension-nparam+1],tmp);
 184             if (!check_poly(S, CS->next, EP, exist, nparam, pos+1, z, options)) {
 185                 value_clear(c);
 186                 value_clear(tmp);
 187                 value_clear(LB);
 188                 value_clear(UB);
 189                 return 0;
 190             }
 191         }
 192         value_set_si(z[pos+S->Dimension-nparam+1], 0);
 193         value_clear(LB);
 194         value_clear(UB);
 195     }
 196     value_clear(c);
 197     value_clear(tmp);
 198     return 1;
 199 } /* check_poly */