4coins.cc

   1 #include <assert.h>
   2 #include <iostream>
   3 #include <barvinok/genfun.h>
   4 #include <barvinok/util.h>
   5 #include <barvinok/barvinok.h>
   6 #include <barvinok/basis_reduction.h>
   7 #include "config.h"
   8
   9 using std::cout;
  10 using std::cerr;
  11 using std::endl;
  12
  13 static Polyhedron *uncone(Polyhedron *C, unsigned MaxRays)
  14 {
  15     Matrix *Constraints;
  16     Polyhedron *P;
  17     int i;
  18
  19     Constraints = Matrix_Alloc(C->NbConstraints, 1+(C->Dimension-1)+1);
  20     for (i = 0; i < C->NbConstraints; ++i) {
  21         /* positivity constraints */
  22         if (First_Non_Zero(C->Constraint[i]+1, C->Dimension) == -1) {
  23             assert(i == C->NbConstraints-1);
  24             Constraints->NbRows = i;
  25             break;
  26         }
  27         assert(value_zero_p(C->Constraint[i][1+C->Dimension]));
  28         Vector_Copy(C->Constraint[i], Constraints->p[i], 1+C->Dimension);
  29     }
  30     P = Constraints2Polyhedron(Constraints, MaxRays);
  31     Matrix_Free(Constraints);
  32
  33     return P;
  34 }
  35
  36 int main(int argc, char **argv)
  37 {
  38     Matrix *M;
  39     Polyhedron *P, *C, *B;
  40     Vector *small;
  41     Matrix *T, *basis;
  42     gen_fun *S, *S_shift, *hp, *S_divide, *frob;
  43     vec_ZZ up;
  44     Value c;
  45     bool ok;
  46     QQ mone(-1, 1);
  47     Vector *coins;
  48     barvinok_options *options = barvinok_options_new_with_defaults();
  49
  50     coins = Vector_Read();
  51     assert(coins);
  52     assert(coins->Size == 4);
  53
  54     M = Matrix_Alloc(5, 7);
  55     Vector_Copy(coins->p, M->p[0]+1, 4);
  56     Vector_Free(coins);
  57     value_set_si(M->p[0][1+4], -1);
  58     for (int i = 0; i < 4; ++i) {
  59         value_set_si(M->p[1+i][0], 1);
  60         value_set_si(M->p[1+i][1+i], 1);
  61     }
  62     C = Constraints2Polyhedron(M, options->MaxRays);
  63     Matrix_Free(M);
  64     C = remove_equalities_p(C, 4, NULL, options->MaxRays);
  65     assert(C->NbEq == 0);
  66     Polyhedron_Print(stderr, P_VALUE_FMT, C);
  67
  68     B = uncone(C, options->MaxRays);
  69
  70     basis = Polyhedron_Reduced_Basis(B, options);
  71     small = Vector_Alloc(B->Dimension + 2);
  72     Vector_Copy(basis->p[0], small->p, B->Dimension);
  73     Matrix_Free(basis);
  74     Polyhedron_Free(B);
  75
  76     T = Matrix_Alloc(small->Size, small->Size);
  77     Vector_Copy(small->p, T->p[0], small->Size);
  78     ok = unimodular_complete(T, 1);
  79     assert(ok);
  80     Vector_Free(small);
  81     Vector_Exchange(T->p[0], T->p[2], T->NbColumns);
  82     Matrix_Print(stderr, P_VALUE_FMT, T);
  83
  84     Polyhedron_Print(stderr, P_VALUE_FMT, C);
  85     P = Polyhedron_Image(C, T, options->MaxRays);
  86     Matrix_Free(T);
  87     Polyhedron_Free(C);
  88
  89     Polyhedron_Print(stderr, P_VALUE_FMT, P);
  90
  91     up.SetLength(2);
  92     up[0] = 1;
  93     up[1] = 0;
  94
  95     S = barvinok_enumerate_scarf_series(P, 2, 2, options);
  96     S->print(std::cerr, 0, NULL);
  97     cerr << endl;
  98
  99     S_shift = new gen_fun(S);
 100     S_divide = new gen_fun(S);
 101     S_divide->divide(up);
 102     S_divide->shift(up);
 103
 104     value_init(c);
 105
 106     do {
 107         S_shift->shift(up);
 108         hp = S->Hadamard_product(S_shift, options);
 109         S->add(mone, hp, options);
 110         delete hp;
 111
 112         S_divide->shift(up);
 113         hp = S->Hadamard_product(S_divide, options);
 114         ok = hp->summate(&c);
 115         if (ok)
 116             assert(value_zero_p(c));
 117         delete hp;
 118     } while(!ok);
 119
 120     value_clear(c);
 121
 122     frob = S->summate(1, options);
 123     frob->print(std::cout, 0, NULL);
 124     cout << endl;
 125     delete frob;
 126
 127     delete S_shift;
 128     delete S_divide;
 129     delete S;
 130
 131     Polyhedron_Free(P);
 132     barvinok_options_free(options);
 133 }