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 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 int main(int argc, char **argv)
  57 {
  58     Matrix *M;
  59     Polyhedron *P, *C, *B, *Pold;
  60     Vector *small;
  61     Matrix *T, *basis;
  62     gen_fun *S, *S_shift, *hp, *S_divide, *frob;
  63     vec_ZZ up;
  64     Value c;
  65     bool ok;
  66     QQ mone(-1, 1);
  67     Vector *coins;
  68     barvinok_options *options = barvinok_options_new_with_defaults();
  69
  70     coins = Vector_Read();
  71     assert(coins);
  72     assert(coins->Size == 4);
  73
  74     M = Matrix_Alloc(5, 7);
  75     Vector_Copy(coins->p, M->p[0]+1, 4);
  76     Vector_Free(coins);
  77     value_set_si(M->p[0][1+4], -1);
  78     for (int i = 0; i < 4; ++i) {
  79         value_set_si(M->p[1+i][0], 1);
  80         value_set_si(M->p[1+i][1+i], 1);
  81     }
  82     C = Constraints2Polyhedron(M, options->MaxRays);
  83     Matrix_Free(M);
  84     C = remove_equalities_p(C, 4, NULL, options->MaxRays);
  85     assert(C->NbEq == 0);
  86     Polyhedron_Print(stderr, P_VALUE_FMT, C);
  87
  88     B = uncone(C, options->MaxRays);
  89
  90     basis = Polyhedron_Reduced_Basis(B, options);
  91     small = Vector_Alloc(B->Dimension + 2);
  92     Vector_Copy(basis->p[0], small->p, B->Dimension);
  93     Matrix_Free(basis);
  94     Polyhedron_Free(B);
  95
  96     T = Matrix_Alloc(small->Size, small->Size);
  97     Vector_Copy(small->p, T->p[0], small->Size);
  98     ok = unimodular_complete(T, 1);
  99     assert(ok);
 100     Vector_Free(small);
 101     Vector_Exchange(T->p[0], T->p[2], T->NbColumns);
 102     Matrix_Print(stderr, P_VALUE_FMT, T);
 103
 104     Polyhedron_Print(stderr, P_VALUE_FMT, C);
 105     P = Polyhedron_Image(C, T, options->MaxRays);
 106     Matrix_Free(T);
 107     Polyhedron_Free(C);
 108
 109     Polyhedron_Print(stderr, P_VALUE_FMT, P);
 110
 111     up.SetLength(2);
 112     up[0] = 1;
 113     up[1] = 0;
 114
 115     S = barvinok_enumerate_scarf_series(P, 2, 2, options);
 116     S->print(std::cerr, 0, NULL);
 117     cerr << endl;
 118
 119     S_shift = new gen_fun(S);
 120     S_divide = new gen_fun(S);
 121     S_divide->divide(up);
 122     S_divide->shift(up);
 123
 124     value_init(c);
 125
 126     do {
 127         S_shift->shift(up);
 128         hp = S->Hadamard_product(S_shift, options);
 129         S->add(mone, hp, options);
 130         delete hp;
 131
 132         S_divide->shift(up);
 133         hp = S->Hadamard_product(S_divide, options);
 134         ok = hp->summate(&c);
 135         if (ok)
 136             assert(value_zero_p(c));
 137         delete hp;
 138     } while(!ok);
 139
 140     value_clear(c);
 141
 142     frob = S->summate(1, options);
 143     frob->print(std::cout, 0, NULL);
 144     cout << endl;
 145     delete frob;
 146
 147     delete S_shift;
 148     delete S_divide;
 149     delete S;
 150
 151     Polyhedron_Free(P);
 152     barvinok_options_free(options);
 153 }