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