4coins.cc

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