4coins.cc

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