update isl for change in space of isl_pw_* objects
[barvinok.git] / 4coins.cc
blob11c05297546c095bb92de699e95e00a5b2979096
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 int main(int argc, char **argv)
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();
70 coins = Vector_Read();
71 assert(coins);
72 assert(coins->Size == 4);
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);
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);
88 B = uncone(C, options->MaxRays);
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);
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);
104 Polyhedron_Print(stderr, P_VALUE_FMT, C);
105 P = Polyhedron_Image(C, T, options->MaxRays);
106 Matrix_Free(T);
107 Polyhedron_Free(C);
109 Polyhedron_Print(stderr, P_VALUE_FMT, P);
111 up.SetLength(2);
112 up[0] = 1;
113 up[1] = 0;
115 S = barvinok_enumerate_scarf_series(P, 2, 2, options);
116 S->print(std::cerr, 0, NULL);
117 cerr << endl;
119 S_shift = new gen_fun(S);
120 S_divide = new gen_fun(S);
121 S_divide->divide(up);
122 S_divide->shift(up);
124 value_init(c);
126 do {
127 S_shift->shift(up);
128 hp = S->Hadamard_product(S_shift, options);
129 S->add(mone, hp, options);
130 delete hp;
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);
140 value_clear(c);
142 frob = S->summate(1, options);
143 frob->print(std::cout, 0, NULL);
144 cout << endl;
145 delete frob;
147 delete S_shift;
148 delete S_divide;
149 delete S;
151 Polyhedron_Free(P);
152 barvinok_options_free(options);