keep track of number of Bernoulli sums
[barvinok.git] / 4coins.cc
blob88d00a91d2ebcda20fbcf22f168268c8a3d6b347
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"
9 using std::cout;
10 using std::cerr;
11 using std::endl;
13 static Polyhedron *uncone(Polyhedron *C, unsigned MaxRays)
15 Matrix *Constraints;
16 Polyhedron *P;
17 int i;
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;
27 assert(value_zero_p(C->Constraint[i][1+C->Dimension]));
28 Vector_Copy(C->Constraint[i], Constraints->p[i], 1+C->Dimension);
30 P = Constraints2Polyhedron(Constraints, MaxRays);
31 Matrix_Free(Constraints);
33 return P;
36 void static scan(gen_fun *S)
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");
52 value_clear(c);
53 Vector_Free(params);
56 bool is_indicator(gen_fun *gf, barvinok_options *options)
58 gen_fun *test;
59 QQ mone(-1, 1);
61 test = gf->Hadamard_product(gf, options);
62 test->add(mone, gf, options);
64 test->print(std::cerr, 0, NULL);
65 cerr << endl;
68 int main(int argc, char **argv)
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();
82 coins = Vector_Read();
83 assert(coins);
84 assert(coins->Size == 4);
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);
94 C = Constraints2Polyhedron(M, options->MaxRays);
95 Matrix_Free(M);
96 C = remove_equalities_p(C, 4, NULL, options->MaxRays);
97 assert(C->NbEq == 0);
98 Polyhedron_Print(stderr, P_VALUE_FMT, C);
100 B = uncone(C, options->MaxRays);
102 basis = Polyhedron_Reduced_Basis(B, options);
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);
108 T = Matrix_Alloc(small->Size, small->Size);
109 Vector_Copy(small->p, T->p[0], small->Size);
110 ok = unimodular_complete(T, 1);
111 assert(ok);
112 Vector_Free(small);
113 Vector_Exchange(T->p[0], T->p[2], T->NbColumns);
114 Matrix_Print(stderr, P_VALUE_FMT, T);
116 Polyhedron_Print(stderr, P_VALUE_FMT, C);
117 P = Polyhedron_Image(C, T, options->MaxRays);
118 Matrix_Free(T);
119 Polyhedron_Free(C);
121 Polyhedron_Print(stderr, P_VALUE_FMT, P);
123 up.SetLength(2);
124 up[0] = 1;
125 up[1] = 0;
127 S = barvinok_enumerate_scarf_series(P, 2, 2, options);
128 S->print(std::cerr, 0, NULL);
129 cerr << endl;
131 S_shift = new gen_fun(S);
132 S_divide = new gen_fun(S);
133 S_divide->divide(up);
134 S_divide->shift(up);
136 value_init(c);
138 do {
139 S_shift->shift(up);
140 hp = S->Hadamard_product(S_shift, options);
141 S->add(mone, hp, options);
142 delete hp;
144 S_divide->shift(up);
145 hp = S->Hadamard_product(S_divide, options);
146 ok = hp->summate(&c);
147 if (ok)
148 assert(value_zero_p(c));
149 delete hp;
150 } while(!ok);
152 value_clear(c);
154 frob = S->summate(1, options);
155 frob->print(std::cout, 0, NULL);
156 cout << endl;
157 delete frob;
159 delete S_shift;
160 delete S_divide;
161 delete S;
163 Polyhedron_Free(P);
164 barvinok_options_free(options);