barvinok_ehrhart.cc

   1 #include <unistd.h>
   2 #include <stdlib.h>
   3 #include <strings.h>
   4 #include <polylib/polylibgmp.h>
   5 #include <barvinok/util.h>
   6 #include <barvinok/barvinok.h>
   7 #include "config.h"
   8
   9 /* The input of this example program is a polytope in PolyLib notation,
  10  * i.e., an n by d+2 matrix of the n constraints A x + b >= 0 defining
  11  * the polytope * sitting in a d-dimensional space.  The first column
  12  * is 1 for an inequality and 0 for an equality.  b is placed in the
  13  * final column.
  14  * Alternatively, if the matrix is preceded by the word "vertices"
  15  * on a line by itself, it will be interpreted as a list of vertices
  16  * in PolyLib notation, i.e., an n by (d+2) matrix, where n is
  17  * the number of vertices/rays and d the dimension.  The first column is
  18  * 0 for lines and 1 for vertices/rays.  The final column is the denominator
  19  * or 0 for rays.  Note that for barvinok_ehrhart, the first column
  20  * should always be 1.
  21  */
  22
  23 #ifdef HAVE_GROWING_CHERNIKOVA
  24 #define MAXRAYS    POL_NO_DUAL
  25 #else
  26 #define MAXRAYS  600
  27 #endif
  28
  29 #ifndef HAVE_GETOPT_H
  30 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
  31 #else
  32 #include <getopt.h>
  33 struct option options[] = {
  34     { "convert",   no_argument,  0,  'c' },
  35     { "floor",     no_argument,  0,  'f' },
  36     { "series",    no_argument,  0,  's' },
  37     { "version",   no_argument,  0,  'V' },
  38     { 0, 0, 0, 0 }
  39 };
  40 #endif
  41
  42 static Polyhedron *Polyhedron_Read()
  43 {
  44     int vertices = 0;
  45     unsigned NbRows, NbColumns;
  46     Matrix *M;
  47     Polyhedron *P;
  48     char s[128];
  49
  50     while (fgets(s, sizeof(s), stdin)) {
  51         if (*s == '#')
  52             continue;
  53         if (strncasecmp(s, "vertices", sizeof("vertices")-1) == 0)
  54             vertices = 1;
  55         if (sscanf(s, "%u %u", &NbRows, &NbColumns) == 2)
  56             break;
  57     }
  58     if (feof(stdin))
  59         return NULL;
  60     M = Matrix_Alloc(NbRows,NbColumns);
  61     Matrix_Read_Input(M);
  62     if (vertices)
  63         P = Rays2Polyhedron(M, MAXRAYS);
  64     else
  65         P = Constraints2Polyhedron(M, MAXRAYS);
  66     Matrix_Free(M);
  67     return P;
  68 }
  69
  70 int main(int argc, char **argv)
  71 {
  72     Polyhedron *A, *C, *U;
  73     char **param_name;
  74     int c, ind = 0;
  75     int convert = 0;
  76     int floor = 0;
  77     int series = 0;
  78
  79     while ((c = getopt_long(argc, argv, "sfcV", options, &ind)) != -1) {
  80         switch (c) {
  81         case 's':
  82             series = 1;
  83             break;
  84         case 'c':
  85             convert = 1;
  86             break;
  87         case 'f':
  88             floor = 1;
  89             break;
  90         case 'V':
  91             printf(barvinok_version());
  92             exit(0);
  93             break;
  94         }
  95     }
  96
  97     A = Polyhedron_Read();
  98     param_name = Read_ParamNames(stdin, 1);
  99     Polyhedron_Print(stdout, P_VALUE_FMT, A);
 100     C = Cone_over_Polyhedron(A);
 101     U = Universe_Polyhedron(1);
 102     if (series) {
 103         gen_fun *gf;
 104         gf = barvinok_series(C, U, MAXRAYS);
 105         gf->print(U->Dimension, param_name);
 106         puts("");
 107         delete gf;
 108     } else {
 109         evalue *EP;
 110         /* A (conceptually) obvious optimization would be to pass in
 111          * the parametric vertices, which are just n times the original
 112          * vertices, rather than letting barvinok_enumerate_ev (re)compute
 113          * them through Polyhedron2Param_SimplifiedDomain.
 114          */
 115         EP = barvinok_enumerate_ev(C, U, MAXRAYS);
 116         print_evalue(stdout, EP, param_name);
 117         if (floor) {
 118             fprintf(stderr, "WARNING: floor conversion not supported\n");
 119             evalue_frac2floor(EP);
 120             print_evalue(stdout, EP, param_name);
 121         } else if (convert) {
 122             evalue_mod2table(EP, C->Dimension);
 123             print_evalue(stdout, EP, param_name);
 124         }
 125         free_evalue_refs(EP);
 126         free(EP);
 127     }
 128     Free_ParamNames(param_name, 1);
 129     Polyhedron_Free(A);
 130     Polyhedron_Free(C);
 131     Polyhedron_Free(U);
 132     return 0;
 133 }