add an ehrhart example with interesting chambers
[barvinok.git] / reducer.cc
blob5fc454ae3943bd2eec02fed264cf6fc8a49672e4
1 #include <vector>
2 #include <barvinok/util.h>
3 #include "reducer.h"
4 #include "lattice_point.h"
6 using std::vector;
7 using std::cerr;
8 using std::endl;
10 struct OrthogonalException Orthogonal;
12 void np_base::handle(const signed_cone& sc, barvinok_options *options)
14 assert(sc.rays.NumRows() == dim);
15 factor.n *= sc.sign;
16 handle(sc.rays, current_vertex, factor, sc.det, sc.closed, options);
17 factor.n *= sc.sign;
20 void np_base::start(Polyhedron *P, barvinok_options *options)
22 QQ factor(1, 1);
23 for (;;) {
24 try {
25 init(P);
26 for (int i = 0; i < P->NbRays; ++i) {
27 if (!value_pos_p(P->Ray[i][dim+1]))
28 continue;
30 Polyhedron *C = supporting_cone(P, i);
31 do_vertex_cone(factor, C, P->Ray[i]+1, options);
33 break;
34 } catch (OrthogonalException &e) {
35 reset();
40 /* input:
41 * f: the powers in the denominator for the remaining vars
42 * each row refers to a factor
43 * den_s: for each factor, the power of (s+1)
44 * sign
45 * num_s: powers in the numerator corresponding to the summed vars
46 * num_p: powers in the numerator corresponding to the remaining vars
47 * number of rays in cone: "dim" = "k"
48 * length of each ray: "dim" = "d"
49 * for now, it is assumed: k == d
50 * output:
51 * den_p: for each factor
52 * 0: independent of remaining vars
53 * 1: power corresponds to corresponding row in f
55 * all inputs are subject to change
57 void normalize(ZZ& sign, vec_ZZ& num_s, mat_ZZ& num_p, vec_ZZ& den_s, vec_ZZ& den_p,
58 mat_ZZ& f)
60 unsigned dim = f.NumRows();
61 unsigned nparam = num_p.NumCols();
62 unsigned nvar = dim - nparam;
64 int change = 0;
66 for (int j = 0; j < den_s.length(); ++j) {
67 if (den_s[j] == 0) {
68 den_p[j] = 1;
69 continue;
71 int k;
72 for (k = 0; k < nparam; ++k)
73 if (f[j][k] != 0)
74 break;
75 if (k < nparam) {
76 den_p[j] = 1;
77 if (den_s[j] > 0) {
78 f[j] = -f[j];
79 for (int i = 0; i < num_p.NumRows(); ++i)
80 num_p[i] += f[j];
82 } else
83 den_p[j] = 0;
84 if (den_s[j] > 0)
85 change ^= 1;
86 else {
87 den_s[j] = abs(den_s[j]);
88 for (int i = 0; i < num_p.NumRows(); ++i)
89 num_s[i] += den_s[j];
93 if (change)
94 sign = -sign;
97 void reducer::base(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den_f)
99 for (int i = 0; i < num.NumRows(); ++i)
100 base(c[i], num[i], den_f);
103 struct dpoly_r_scanner {
104 const dpoly_r *rc;
105 const dpoly * const *num;
106 int n;
107 int dim;
108 dpoly_r_term_list::iterator *iter;
109 vector<int> powers;
110 vec_ZZ coeff;
112 dpoly_r_scanner(const dpoly * const *num, int n, const dpoly_r *rc, int dim)
113 : num(num), rc(rc), n(n), dim(dim), powers(dim, 0) {
114 coeff.SetLength(n);
115 iter = new dpoly_r_term_list::iterator[rc->len];
116 for (int i = 0; i < rc->len; ++i) {
117 int k;
118 for (k = 0; k < n; ++k)
119 if (num[k]->coeff[rc->len-1-i] != 0)
120 break;
121 if (k < n)
122 iter[i] = rc->c[i].begin();
123 else
124 iter[i] = rc->c[i].end();
127 bool next() {
128 int pos[rc->len];
129 int len = 0;
131 for (int i = 0; i < rc->len; ++i) {
132 if (iter[i] == rc->c[i].end())
133 continue;
134 if (!len)
135 pos[len++] = i;
136 else {
137 if ((*iter[i])->powers < (*iter[pos[0]])->powers) {
138 pos[0] = i;
139 len = 1;
140 } else if ((*iter[i])->powers == (*iter[pos[0]])->powers)
141 pos[len++] = i;
145 if (!len)
146 return false;
148 powers = (*iter[pos[0]])->powers;
149 for (int k = 0; k < n; ++k)
150 mul(coeff[k], (*iter[pos[0]])->coeff, num[k]->coeff[rc->len-1-pos[0]]);
151 ++iter[pos[0]];
152 for (int i = 1; i < len; ++i) {
153 for (int k = 0; k < n; ++k) {
154 mul(tmp, (*iter[pos[i]])->coeff, num[k]->coeff[rc->len-1-pos[i]]);
155 add(coeff[k], coeff[k], tmp);
157 ++iter[pos[i]];
160 return true;
162 ~dpoly_r_scanner() {
163 delete [] iter;
165 private:
166 ZZ tmp;
169 void reducer::reduce(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den_f)
171 assert(c.length() == num.NumRows());
172 unsigned len = den_f.NumRows(); // number of factors in den
173 vec_QQ c2 = c;
175 if (num.NumCols() == lower) {
176 base(c, num, den_f);
177 return;
179 assert(num.NumCols() > 1);
180 assert(num.NumRows() > 0);
182 vec_ZZ den_s;
183 mat_ZZ den_r;
184 vec_ZZ num_s;
185 mat_ZZ num_p;
187 split(num, num_s, num_p, den_f, den_s, den_r);
189 vec_ZZ den_p;
190 den_p.SetLength(len);
192 ZZ sign(INIT_VAL, 1);
193 normalize(sign, num_s, num_p, den_s, den_p, den_r);
194 c2 *= sign;
196 int only_param = 0; // k-r-s from text
197 int no_param = 0; // r from text
198 for (int k = 0; k < len; ++k) {
199 if (den_p[k] == 0)
200 ++no_param;
201 else if (den_s[k] == 0)
202 ++only_param;
204 if (no_param == 0) {
205 reduce(c2, num_p, den_r);
206 } else {
207 int k, l;
208 mat_ZZ pden;
209 pden.SetDims(only_param, den_r.NumCols());
211 for (k = 0, l = 0; k < len; ++k)
212 if (den_s[k] == 0)
213 pden[l++] = den_r[k];
215 for (k = 0; k < len; ++k)
216 if (den_p[k] == 0)
217 break;
219 dpoly *n[num_s.length()];
220 for (int i = 0; i < num_s.length(); ++i) {
221 n[i] = new dpoly(no_param, num_s[i]);
222 /* Search for other numerator (j) with same num_p.
223 * If found, replace a[j]/b[j] * n[j] and a[i]/b[i] * n[i]
224 * by 1/(b[j]*b[i]/g) * (a[j]*b[i]/g * n[j] + a[i]*b[j]/g * n[i])
225 * where g = gcd(b[i], b[j].
227 for (int j = 0; j < i; ++j) {
228 if (num_p[i] != num_p[j])
229 continue;
230 ZZ g = GCD(c2[i].d, c2[j].d);
231 *n[j] *= c2[j].n * c2[i].d/g;
232 *n[i] *= c2[i].n * c2[j].d/g;
233 *n[j] += *n[i];
234 c2[j].n = 1;
235 c2[j].d *= c2[i].d/g;
236 delete n[i];
237 if (i < num_s.length()-1) {
238 num_s[i] = num_s[num_s.length()-1];
239 num_p[i] = num_p[num_s.length()-1];
240 c2[i] = c2[num_s.length()-1];
242 num_s.SetLength(num_s.length()-1);
243 c2.SetLength(c2.length()-1);
244 num_p.SetDims(num_p.NumRows()-1, num_p.NumCols());
245 --i;
246 break;
249 dpoly D(no_param, den_s[k], 1);
250 for ( ; ++k < len; )
251 if (den_p[k] == 0) {
252 dpoly fact(no_param, den_s[k], 1);
253 D *= fact;
256 if (no_param + only_param == len) {
257 vec_QQ q;
258 q.SetLength(num_s.length());
259 for (int i = 0; i < num_s.length(); ++i) {
260 mpq_set_si(tcount, 0, 1);
261 n[i]->div(D, tcount, one);
263 value2zz(mpq_numref(tcount), q[i].n);
264 value2zz(mpq_denref(tcount), q[i].d);
265 q[i] *= c2[i];
267 for (int i = q.length()-1; i >= 0; --i) {
268 if (q[i].n == 0) {
269 q[i] = q[q.length()-1];
270 num_p[i] = num_p[q.length()-1];
271 q.SetLength(q.length()-1);
272 num_p.SetDims(num_p.NumRows()-1, num_p.NumCols());
276 if (q.length() != 0)
277 reduce(q, num_p, pden);
278 } else {
279 ZZ zz_zero(INIT_VAL, 0);
280 dpoly one(no_param, zz_zero);
281 dpoly_r *r = NULL;
283 for (k = 0; k < len; ++k) {
284 if (den_s[k] == 0 || den_p[k] == 0)
285 continue;
287 dpoly pd(no_param-1, den_s[k], 1);
289 int l;
290 for (l = 0; l < k; ++l)
291 if (den_r[l] == den_r[k])
292 break;
294 if (!r)
295 r = new dpoly_r(one, pd, l, len);
296 else {
297 dpoly_r *nr = new dpoly_r(r, pd, l, len);
298 delete r;
299 r = nr;
303 vec_QQ factor;
304 factor.SetLength(c2.length());
305 int common = pden.NumRows();
306 dpoly_r *rc = r->div(D);
307 for (int i = 0; i < num_s.length(); ++i) {
308 factor[i].d = c2[i].d;
309 factor[i].d *= rc->denom;
312 dpoly_r_scanner scanner(n, num_s.length(), rc, len);
313 int rows;
314 while (scanner.next()) {
315 int i;
316 for (i = 0; i < num_s.length(); ++i)
317 if (scanner.coeff[i] != 0)
318 break;
319 if (i == num_s.length())
320 continue;
321 rows = common;
322 pden.SetDims(rows, pden.NumCols());
323 for (int k = 0; k < rc->dim; ++k) {
324 int n = scanner.powers[k];
325 if (n == 0)
326 continue;
327 pden.SetDims(rows+n, pden.NumCols());
328 for (int l = 0; l < n; ++l)
329 pden[rows+l] = den_r[k];
330 rows += n;
332 for (int i = 0; i < num_s.length(); ++i) {
333 factor[i].n = c2[i].n;
334 factor[i].n *= scanner.coeff[i];
336 reduce(factor, num_p, pden);
339 delete rc;
340 delete r;
342 for (int i = 0; i < num_s.length(); ++i)
343 delete n[i];
347 void reducer::handle(const mat_ZZ& den, Value *V, const QQ& c, unsigned long det,
348 int *closed, barvinok_options *options)
350 vec_QQ vc;
352 lattice_point(V, den, vertex, det, closed);
354 vc.SetLength(vertex.NumRows());
355 for (int i = 0; i < vc.length(); ++i)
356 vc[i] = c;
358 reduce(vc, vertex, den);
361 void split_one(const mat_ZZ& num, vec_ZZ& num_s, mat_ZZ& num_p,
362 const mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r)
364 unsigned len = den_f.NumRows(); // number of factors in den
365 unsigned d = num.NumCols() - 1;
367 den_s.SetLength(len);
368 den_r.SetDims(len, d);
370 for (int r = 0; r < len; ++r) {
371 den_s[r] = den_f[r][0];
372 for (int k = 1; k <= d; ++k)
373 den_r[r][k-1] = den_f[r][k];
376 num_s.SetLength(num.NumRows());
377 num_p.SetDims(num.NumRows(), d);
378 for (int i = 0; i < num.NumRows(); ++i) {
379 num_s[i] = num[i][0];
380 for (int k = 1 ; k <= d; ++k)
381 num_p[i][k-1] = num[i][k];
385 void normalize(ZZ& sign, ZZ& num, vec_ZZ& den)
387 unsigned dim = den.length();
389 int change = 0;
391 for (int j = 0; j < den.length(); ++j) {
392 if (den[j] > 0)
393 change ^= 1;
394 else {
395 den[j] = abs(den[j]);
396 num += den[j];
399 if (change)
400 sign = -sign;
403 void icounter::base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
405 int r;
406 unsigned len = den_f.NumRows(); // number of factors in den
407 vec_ZZ den_s;
408 den_s.SetLength(len);
409 assert(num.length() == 1);
410 ZZ num_s = num[0];
411 for (r = 0; r < len; ++r)
412 den_s[r] = den_f[r][0];
413 ZZ sign = ZZ(INIT_VAL, 1);
414 normalize(sign, num_s, den_s);
416 dpoly n(len, num_s);
417 dpoly D(len, den_s[0], 1);
418 for (int k = 1; k < len; ++k) {
419 dpoly fact(len, den_s[k], 1);
420 D *= fact;
422 mpq_set_si(tcount, 0, 1);
423 n.div(D, tcount, one);
424 zz2value(c.n, tn);
425 if (sign == -1)
426 value_oppose(tn, tn);
427 zz2value(c.d, td);
428 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
429 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
430 mpq_canonicalize(tcount);
431 mpq_add(count, count, tcount);
434 void infinite_icounter::base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
436 int r;
437 unsigned len = den_f.NumRows(); // number of factors in den
438 vec_ZZ den_s;
439 den_s.SetLength(len);
440 assert(num.length() == 1);
441 ZZ num_s = num[0];
443 for (r = 0; r < len; ++r)
444 den_s[r] = den_f[r][0];
445 ZZ sign = ZZ(INIT_VAL, 1);
446 normalize(sign, num_s, den_s);
448 dpoly n(len, num_s);
449 dpoly D(len, den_s[0], 1);
450 for (int k = 1; k < len; ++k) {
451 dpoly fact(len, den_s[k], 1);
452 D *= fact;
455 Value tmp;
456 mpq_t factor;
457 mpq_init(factor);
458 value_init(tmp);
459 zz2value(c.n, tmp);
460 if (sign == -1)
461 value_oppose(tmp, tmp);
462 value_assign(mpq_numref(factor), tmp);
463 zz2value(c.d, tmp);
464 value_assign(mpq_denref(factor), tmp);
466 n.div(D, count, factor);
468 value_clear(tmp);
469 mpq_clear(factor);