polytope_scan.c

   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/options.h>
   7 #include <barvinok/basis_reduction.h>
   8 #include "config.h"
   9
  10 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
  11
  12 #ifndef HAVE_GETOPT_H
  13 #define getopt_long(a,b,c,d,e) getopt(a,b,c)
  14 #else
  15 #include <getopt.h>
  16 struct option options[] = {
  17     { "direct",    no_argument,  0,  'd' },
  18     { "version",   no_argument,  0,  'V' },
  19     { 0, 0, 0, 0 }
  20 };
  21 #endif
  22
  23 static void scan_poly(Polyhedron *S, int pos, Value *z, Matrix *T)
  24 {
  25     if (!S) {
  26         int k;
  27         Vector *v;
  28
  29         v = Vector_Alloc(T->NbRows);
  30         Matrix_Vector_Product(T, z+1, v->p);
  31         value_print(stdout, VALUE_FMT, v->p[0]);
  32         for (k=1; k < pos; ++k) {
  33             printf(", ");
  34             value_print(stdout,VALUE_FMT, v->p[k]);
  35         }
  36         Vector_Free(v);
  37         printf("\n");
  38     } else {
  39         int ok;
  40         Value LB, UB, tmp;
  41         value_init(LB);
  42         value_init(UB);
  43         value_init(tmp);
  44         ok = !(lower_upper_bounds(1+pos, S, z, &LB, &UB));
  45         assert(ok);
  46         for (value_assign(tmp,LB); value_le(tmp,UB); value_increment(tmp,tmp)) {
  47             value_assign(z[pos+1], tmp);
  48             scan_poly(S->next, pos+1, z, T);
  49         }
  50         value_set_si(z[pos+1], 0);
  51         value_clear(LB);
  52         value_clear(UB);
  53         value_clear(tmp);
  54     }
  55 }
  56
  57 int main(int argc, char **argv)
  58 {
  59     Polyhedron *A, *P, *U, *S;
  60     Value *p;
  61     int i, j, ok;
  62     Matrix *basis, *T, *inv;
  63     int c, ind = 0;
  64     int direct = 0;
  65     struct barvinok_options *bv_options = barvinok_options_new_with_defaults();
  66
  67     while ((c = getopt_long(argc, argv, "dV", options, &ind)) != -1) {
  68         switch (c) {
  69         case 'd':
  70             direct = 1;
  71             break;
  72         case 'V':
  73             printf(barvinok_version());
  74             exit(0);
  75             break;
  76         }
  77     }
  78
  79     A = Polyhedron_Read(bv_options->MaxRays);
  80
  81     if (direct) {
  82         inv = Identity(A->Dimension+1);
  83         P = A;
  84     } else {
  85         basis = Polyhedron_Reduced_Basis(A);
  86
  87         T = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
  88         inv = Matrix_Alloc(A->Dimension+1, A->Dimension+1);
  89         for (i = 0; i < A->Dimension; ++i)
  90             for (j = 0; j < A->Dimension; ++j)
  91                 value_assign(T->p[i][j], basis->p[i][j]);
  92         value_set_si(T->p[A->Dimension][A->Dimension], 1);
  93         Matrix_Free(basis);
  94
  95         ok = Matrix_Inverse(T, inv);
  96         assert(ok);
  97         Matrix_Free(T);
  98
  99         P = Polyhedron_Preimage(A, inv, bv_options->MaxRays);
 100         Polyhedron_Free(A);
 101     }
 102
 103     U = Universe_Polyhedron(0);
 104     S = Polyhedron_Scan(P, U, bv_options->MaxRays);
 105
 106     p = ALLOCN(Value, P->Dimension+2);
 107     for(i=0;i<=P->Dimension;i++) {
 108         value_init(p[i]);
 109         value_set_si(p[i],0);
 110     }
 111     value_init(p[i]);
 112     value_set_si(p[i],1);
 113
 114     scan_poly(S, 0, p, inv);
 115
 116     Matrix_Free(inv);
 117     for(i=0;i<=(P->Dimension+1);i++)
 118         value_clear(p[i]);
 119     free(p);
 120     Domain_Free(S);
 121     Polyhedron_Free(P);
 122     Polyhedron_Free(U);
 123     free(bv_options);
 124     return 0;
 125 }