bfcounter.cc: move some more code from barvinok.cc
[barvinok.git] / bfcounter.cc
blob820ed87981f62f8a2359e79037fc849db759b843
1 #include <vector>
2 #include "bfcounter.h"
3 #include "lattice_point.h"
5 using std::vector;
7 static int lex_cmp(vec_ZZ& a, vec_ZZ& b)
9 assert(a.length() == b.length());
11 for (int j = 0; j < a.length(); ++j)
12 if (a[j] != b[j])
13 return a[j] < b[j] ? -1 : 1;
14 return 0;
17 void bf_base::add_term(bfc_term_base *t, vec_ZZ& num_orig, vec_ZZ& extra_num)
19 vec_ZZ num;
20 int d = num_orig.length();
21 num.SetLength(d-1);
22 for (int l = 0; l < d-1; ++l)
23 num[l] = num_orig[l+1] + extra_num[l];
25 add_term(t, num);
28 void bf_base::add_term(bfc_term_base *t, vec_ZZ& num)
30 int len = t->terms.NumRows();
31 int i, r;
32 for (i = 0; i < len; ++i) {
33 r = lex_cmp(t->terms[i], num);
34 if (r >= 0)
35 break;
37 if (i == len || r > 0) {
38 t->terms.SetDims(len+1, num.length());
39 insert_term(t, i);
40 t->terms[i] = num;
41 } else {
42 // i < len && r == 0
43 update_term(t, i);
47 bfc_term_base* bf_base::find_bfc_term(bfc_vec& v, int *powers, int len)
49 bfc_vec::iterator i;
50 for (i = v.begin(); i != v.end(); ++i) {
51 int j;
52 for (j = 0; j < len; ++j)
53 if ((*i)->powers[j] != powers[j])
54 break;
55 if (j == len)
56 return (*i);
57 if ((*i)->powers[j] > powers[j])
58 break;
61 bfc_term_base* t = new_bf_term(len);
62 v.insert(i, t);
63 memcpy(t->powers, powers, len * sizeof(int));
65 return t;
68 void bf_base::reduce(mat_ZZ& factors, bfc_vec& v)
70 assert(v.size() > 0);
71 unsigned nf = factors.NumRows();
72 unsigned d = factors.NumCols();
74 if (d == lower)
75 return base(factors, v);
77 bf_reducer bfr(factors, v, this);
79 bfr.reduce();
81 if (bfr.vn.size() > 0)
82 reduce(bfr.nfactors, bfr.vn);
85 int bf_base::setup_factors(Polyhedron *C, mat_ZZ& factors,
86 bfc_term_base* t, int s)
88 factors.SetDims(dim, dim);
90 int r;
92 for (r = 0; r < dim; ++r)
93 t->powers[r] = 1;
95 for (r = 0; r < dim; ++r) {
96 values2zz(C->Ray[r]+1, factors[r], dim);
97 int k;
98 for (k = 0; k < dim; ++k)
99 if (factors[r][k] != 0)
100 break;
101 if (factors[r][k] < 0) {
102 factors[r] = -factors[r];
103 t->terms[0] += factors[r];
104 s = -s;
108 return s;
111 void bf_base::handle_polar(Polyhedron *C, Value *vertex, QQ c)
113 bfc_term* t = new bfc_term(dim);
114 vector< bfc_term_base * > v;
115 v.push_back(t);
117 t->c.SetLength(1);
119 t->terms.SetDims(1, dim);
120 lattice_point(vertex, C, t->terms[0]);
122 // the elements of factors are always lexpositive
123 mat_ZZ factors;
124 int s = setup_factors(C, factors, t, 1);
126 t->c[0].n = s * c.n;
127 t->c[0].d = c.d;
129 reduce(factors, v);
132 bfc_term_base* bfcounter_base::new_bf_term(int len)
134 bfc_term* t = new bfc_term(len);
135 t->c.SetLength(0);
136 return t;
139 void bfcounter_base::set_factor(bfc_term_base *t, int k, int change)
141 bfc_term* bfct = static_cast<bfc_term *>(t);
142 c = bfct->c[k];
143 if (change)
144 c.n = -c.n;
147 void bfcounter_base::set_factor(bfc_term_base *t, int k, mpq_t &f, int change)
149 bfc_term* bfct = static_cast<bfc_term *>(t);
150 value2zz(mpq_numref(f), c.n);
151 value2zz(mpq_denref(f), c.d);
152 c *= bfct->c[k];
153 if (change)
154 c.n = -c.n;
157 void bfcounter_base::set_factor(bfc_term_base *t, int k, const QQ& c_factor,
158 int change)
160 bfc_term* bfct = static_cast<bfc_term *>(t);
161 c = bfct->c[k];
162 c *= c_factor;
163 if (change)
164 c.n = -c.n;
167 void bfcounter_base::insert_term(bfc_term_base *t, int i)
169 bfc_term* bfct = static_cast<bfc_term *>(t);
170 int len = t->terms.NumRows()-1; // already increased by one
172 bfct->c.SetLength(len+1);
173 for (int j = len; j > i; --j) {
174 bfct->c[j] = bfct->c[j-1];
175 t->terms[j] = t->terms[j-1];
177 bfct->c[i] = c;
180 void bfcounter_base::update_term(bfc_term_base *t, int i)
182 bfc_term* bfct = static_cast<bfc_term *>(t);
184 bfct->c[i] += c;
187 void bf_reducer::compute_extra_num(int i)
189 clear(extra_num);
190 changes = 0;
191 no_param = 0; // r from text
192 only_param = 0; // k-r-s from text
193 total_power = 0; // k from text
195 for (int j = 0; j < nf; ++j) {
196 if (v[i]->powers[j] == 0)
197 continue;
199 total_power += v[i]->powers[j];
200 if (factors[j][0] == 0) {
201 only_param += v[i]->powers[j];
202 continue;
205 if (old2new[j] == -1)
206 no_param += v[i]->powers[j];
207 else
208 extra_num += -sign[j] * v[i]->powers[j] * nfactors[old2new[j]];
209 changes += v[i]->powers[j];
213 void bf_reducer::update_powers(int *powers, int len)
215 for (int l = 0; l < nnf; ++l)
216 npowers[l] = bpowers[l];
218 l_extra_num = extra_num;
219 l_changes = changes;
221 for (int l = 0; l < len; ++l) {
222 int n = powers[l];
223 if (n == 0)
224 continue;
225 assert(old2new[l] != -1);
227 npowers[old2new[l]] += n;
228 // interpretation of sign has been inverted
229 // since we inverted the power for specialization
230 if (sign[l] == 1) {
231 l_extra_num += n * nfactors[old2new[l]];
232 l_changes += n;
238 void bf_reducer::compute_reduced_factors()
240 unsigned nf = factors.NumRows();
241 unsigned d = factors.NumCols();
242 nnf = 0;
243 nfactors.SetDims(nnf, d-1);
245 for (int i = 0; i < nf; ++i) {
246 int j;
247 int s = 1;
248 for (j = 0; j < nnf; ++j) {
249 int k;
250 for (k = 1; k < d; ++k)
251 if (factors[i][k] != 0 || nfactors[j][k-1] != 0)
252 break;
253 if (k < d && factors[i][k] == -nfactors[j][k-1])
254 s = -1;
255 for (; k < d; ++k)
256 if (factors[i][k] != s * nfactors[j][k-1])
257 break;
258 if (k == d)
259 break;
261 old2new[i] = j;
262 if (j == nnf) {
263 int k;
264 for (k = 1; k < d; ++k)
265 if (factors[i][k] != 0)
266 break;
267 if (k < d) {
268 if (factors[i][k] < 0)
269 s = -1;
270 nfactors.SetDims(++nnf, d-1);
271 for (int k = 1; k < d; ++k)
272 nfactors[j][k-1] = s * factors[i][k];
273 } else
274 old2new[i] = -1;
276 sign[i] = s;
278 npowers = new int[nnf];
279 bpowers = new int[nnf];
282 void bf_reducer::reduce()
284 compute_reduced_factors();
286 for (int i = 0; i < v.size(); ++i) {
287 compute_extra_num(i);
289 if (no_param == 0) {
290 vec_ZZ extra_num;
291 extra_num.SetLength(d-1);
292 int changes = 0;
293 int npowers[nnf];
294 for (int k = 0; k < nnf; ++k)
295 npowers[k] = 0;
296 for (int k = 0; k < nf; ++k) {
297 assert(old2new[k] != -1);
298 npowers[old2new[k]] += v[i]->powers[k];
299 if (sign[k] == -1) {
300 extra_num += v[i]->powers[k] * nfactors[old2new[k]];
301 changes += v[i]->powers[k];
305 bfc_term_base * t = bf->find_bfc_term(vn, npowers, nnf);
306 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
307 bf->set_factor(v[i], k, changes % 2);
308 bf->add_term(t, v[i]->terms[k], extra_num);
310 } else {
311 // powers of "constant" part
312 for (int k = 0; k < nnf; ++k)
313 bpowers[k] = 0;
314 for (int k = 0; k < nf; ++k) {
315 if (factors[k][0] != 0)
316 continue;
317 assert(old2new[k] != -1);
318 bpowers[old2new[k]] += v[i]->powers[k];
319 if (sign[k] == -1) {
320 extra_num += v[i]->powers[k] * nfactors[old2new[k]];
321 changes += v[i]->powers[k];
325 int j;
326 for (j = 0; j < nf; ++j)
327 if (old2new[j] == -1 && v[i]->powers[j] > 0)
328 break;
330 dpoly D(no_param, factors[j][0], 1);
331 for (int k = 1; k < v[i]->powers[j]; ++k) {
332 dpoly fact(no_param, factors[j][0], 1);
333 D *= fact;
335 for ( ; ++j < nf; )
336 if (old2new[j] == -1)
337 for (int k = 0; k < v[i]->powers[j]; ++k) {
338 dpoly fact(no_param, factors[j][0], 1);
339 D *= fact;
342 if (no_param + only_param == total_power &&
343 bf->constant_vertex(d)) {
344 bfc_term_base * t = NULL;
345 vec_ZZ num;
346 num.SetLength(d-1);
347 ZZ cn;
348 ZZ cd;
349 for (int k = 0; k < v[i]->terms.NumRows(); ++k) {
350 dpoly n(no_param, v[i]->terms[k][0]);
351 mpq_set_si(bf->tcount, 0, 1);
352 n.div(D, bf->tcount, bf->one);
354 if (value_zero_p(mpq_numref(bf->tcount)))
355 continue;
357 if (!t)
358 t = bf->find_bfc_term(vn, bpowers, nnf);
359 bf->set_factor(v[i], k, bf->tcount, changes % 2);
360 bf->add_term(t, v[i]->terms[k], extra_num);
362 } else {
363 for (int j = 0; j < v[i]->terms.NumRows(); ++j) {
364 dpoly n(no_param, v[i]->terms[j][0]);
366 dpoly_r * r = 0;
367 if (no_param + only_param == total_power)
368 r = new dpoly_r(n, nf);
369 else
370 for (int k = 0; k < nf; ++k) {
371 if (v[i]->powers[k] == 0)
372 continue;
373 if (factors[k][0] == 0 || old2new[k] == -1)
374 continue;
376 dpoly pd(no_param-1, factors[k][0], 1);
378 for (int l = 0; l < v[i]->powers[k]; ++l) {
379 int q;
380 for (q = 0; q < k; ++q)
381 if (old2new[q] == old2new[k] &&
382 sign[q] == sign[k])
383 break;
385 if (r == 0)
386 r = new dpoly_r(n, pd, q, nf);
387 else {
388 dpoly_r *nr = new dpoly_r(r, pd, q, nf);
389 delete r;
390 r = nr;
395 dpoly_r *rc = r->div(D);
396 delete r;
397 QQ factor;
398 factor.d = rc->denom;
400 if (bf->constant_vertex(d)) {
401 vector< dpoly_r_term * >& final = rc->c[rc->len-1];
403 for (int k = 0; k < final.size(); ++k) {
404 if (final[k]->coeff == 0)
405 continue;
407 update_powers(final[k]->powers, rc->dim);
409 bfc_term_base * t = bf->find_bfc_term(vn, npowers, nnf);
410 factor.n = final[k]->coeff;
411 bf->set_factor(v[i], j, factor, l_changes % 2);
412 bf->add_term(t, v[i]->terms[j], l_extra_num);
414 } else
415 bf->cum(this, v[i], j, rc);
417 delete rc;
421 delete v[i];