blob 6f66cba65c989ec8739296065789dc307a2d7cae
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, options);
17 factor.n *= sc.sign;
20 void np_base::start(Polyhedron *P, barvinok_options *options)
22 int n_try = 0;
23 QQ factor(1, 1);
24 for (;;) {
25 try {
26 init(P, n_try);
27 for (int i = 0; i < P->NbRays; ++i) {
28 if (!value_pos_p(P->Ray[i][dim+1]))
29 continue;
31 Polyhedron *C = supporting_cone(P, i);
32 do_vertex_cone(factor, C, P->Ray[i]+1, options);
34 break;
35 } catch (OrthogonalException &e) {
36 n_try++;
37 reset();
42 /* input:
43 * f: the powers in the denominator for the remaining vars
44 * each row refers to a factor
45 * den_s: for each factor, the power of (s+1)
46 * sign
47 * num_s: powers in the numerator corresponding to the summed vars
48 * num_p: powers in the numerator corresponding to the remaining vars
49 * number of rays in cone: "dim" = "k"
50 * length of each ray: "dim" = "d"
51 * for now, it is assumed: k == d
52 * output:
53 * den_p: for each factor
54 * 0: independent of remaining vars
55 * 1: power corresponds to corresponding row in f
57 * all inputs are subject to change
59 void normalize(ZZ& sign, vec_ZZ& num_s, mat_ZZ& num_p, vec_ZZ& den_s, vec_ZZ& den_p,
60 mat_ZZ& f)
62 unsigned dim = f.NumRows();
63 unsigned nparam = num_p.NumCols();
64 unsigned nvar = dim - nparam;
66 int change = 0;
68 for (int j = 0; j < den_s.length(); ++j) {
69 if (den_s[j] == 0) {
70 den_p[j] = 1;
71 continue;
73 int k;
74 for (k = 0; k < nparam; ++k)
75 if (f[j][k] != 0)
76 break;
77 if (k < nparam) {
78 den_p[j] = 1;
79 if (den_s[j] > 0) {
80 f[j] = -f[j];
81 for (int i = 0; i < num_p.NumRows(); ++i)
82 num_p[i] += f[j];
84 } else
85 den_p[j] = 0;
86 if (den_s[j] > 0)
87 change ^= 1;
88 else {
89 den_s[j] = abs(den_s[j]);
90 for (int i = 0; i < num_p.NumRows(); ++i)
91 num_s[i] += den_s[j];
95 if (change)
96 sign = -sign;
99 void reducer::base(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den_f)
101 for (int i = 0; i < num.NumRows(); ++i)
102 base(c[i], num[i], den_f);
105 struct dpoly_r_scanner {
106 const dpoly_r *rc;
107 const dpoly * const *num;
108 int n;
109 int dim;
110 dpoly_r_term_list::iterator *iter;
111 vector<int> powers;
112 vec_ZZ coeff;
114 dpoly_r_scanner(const dpoly * const *num, int n, const dpoly_r *rc, int dim)
115 : num(num), rc(rc), n(n), dim(dim), powers(dim, 0) {
116 coeff.SetLength(n);
117 iter = new dpoly_r_term_list::iterator[rc->len];
118 for (int i = 0; i < rc->len; ++i) {
119 int k;
120 for (k = 0; k < n; ++k)
121 if (value_notzero_p(num[k]->coeff->p[rc->len-1-i]))
122 break;
123 if (k < n)
124 iter[i] = rc->c[i].begin();
125 else
126 iter[i] = rc->c[i].end();
129 bool next() {
130 int *pos;
131 int len = 0;
133 for (int i = 0; i < rc->len; ++i) {
134 if (iter[i] == rc->c[i].end())
135 continue;
136 if (!len) {
137 pos = new int[rc->len];
138 pos[len++] = i;
139 } else {
140 if ((*iter[i])->powers < (*iter[pos[0]])->powers) {
141 pos[0] = i;
142 len = 1;
143 } else if ((*iter[i])->powers == (*iter[pos[0]])->powers)
144 pos[len++] = i;
148 if (!len)
149 return false;
151 powers = (*iter[pos[0]])->powers;
152 for (int k = 0; k < n; ++k) {
153 value2zz(num[k]->coeff->p[rc->len-1-pos[0]], tmp);
154 mul(coeff[k], (*iter[pos[0]])->coeff, tmp);
156 ++iter[pos[0]];
157 for (int i = 1; i < len; ++i) {
158 for (int k = 0; k < n; ++k) {
159 value2zz(num[k]->coeff->p[rc->len-1-pos[i]], tmp);
160 mul(tmp, (*iter[pos[i]])->coeff, tmp);
163 ++iter[pos[i]];
166 delete [] pos;
167 return true;
169 ~dpoly_r_scanner() {
170 delete [] iter;
172 private:
173 ZZ tmp;
176 void reducer::reduce_canonical(const vec_QQ& c, const mat_ZZ& num,
177 const mat_ZZ& den_f)
179 vec_QQ c2 = c;
180 mat_ZZ num2 = num;
182 for (int i = 0; i < c2.length(); ++i) {
183 c2[i].canonicalize();
184 if (c2[i].n != 0)
185 continue;
187 if (i < c2.length()-1) {
188 num2[i] = num2[c2.length()-1];
189 c2[i] = c2[c2.length()-1];
191 num2.SetDims(num2.NumRows()-1, num2.NumCols());
192 c2.SetLength(c2.length()-1);
193 --i;
195 reduce(c2, num2, den_f);
198 void reducer::reduce(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den_f)
200 assert(c.length() == num.NumRows());
201 unsigned len = den_f.NumRows(); // number of factors in den
202 vec_QQ c2 = c;
204 if (num.NumCols() == lower) {
205 base(c, num, den_f);
206 return;
208 assert(num.NumCols() > 1);
209 assert(num.NumRows() > 0);
211 vec_ZZ den_s;
212 mat_ZZ den_r;
213 vec_ZZ num_s;
214 mat_ZZ num_p;
216 split(num, num_s, num_p, den_f, den_s, den_r);
218 vec_ZZ den_p;
219 den_p.SetLength(len);
221 ZZ sign(INIT_VAL, 1);
222 normalize(sign, num_s, num_p, den_s, den_p, den_r);
223 c2 *= sign;
225 int only_param = 0; // k-r-s from text
226 int no_param = 0; // r from text
227 for (int k = 0; k < len; ++k) {
228 if (den_p[k] == 0)
229 ++no_param;
230 else if (den_s[k] == 0)
231 ++only_param;
233 if (no_param == 0) {
234 reduce(c2, num_p, den_r);
235 } else {
236 int k, l;
237 mat_ZZ pden;
238 pden.SetDims(only_param, den_r.NumCols());
240 for (k = 0, l = 0; k < len; ++k)
241 if (den_s[k] == 0)
242 pden[l++] = den_r[k];
244 for (k = 0; k < len; ++k)
245 if (den_p[k] == 0)
246 break;
248 dpoly **n = new dpoly *[num_s.length()];
249 for (int i = 0; i < num_s.length(); ++i) {
250 zz2value(num_s[i], tz);
251 n[i] = new dpoly(no_param, tz);
252 /* Search for other numerator (j) with same num_p.
253 * If found, replace a[j]/b[j] * n[j] and a[i]/b[i] * n[i]
254 * by 1/(b[j]*b[i]/g) * (a[j]*b[i]/g * n[j] + a[i]*b[j]/g * n[i])
255 * where g = gcd(b[i], b[j].
257 for (int j = 0; j < i; ++j) {
258 if (num_p[i] != num_p[j])
259 continue;
260 ZZ g = GCD(c2[i].d, c2[j].d);
261 zz2value(c2[j].n * c2[i].d/g, tz);
262 *n[j] *= tz;
263 zz2value(c2[i].n * c2[j].d/g, tz);
264 *n[i] *= tz;
265 *n[j] += *n[i];
266 c2[j].n = 1;
267 c2[j].d *= c2[i].d/g;
268 delete n[i];
269 if (i < num_s.length()-1) {
270 num_s[i] = num_s[num_s.length()-1];
271 num_p[i] = num_p[num_s.length()-1];
272 c2[i] = c2[num_s.length()-1];
274 num_s.SetLength(num_s.length()-1);
275 c2.SetLength(c2.length()-1);
276 num_p.SetDims(num_p.NumRows()-1, num_p.NumCols());
277 --i;
278 break;
281 zz2value(den_s[k], tz);
282 dpoly D(no_param, tz, 1);
283 for ( ; ++k < len; )
284 if (den_p[k] == 0) {
285 zz2value(den_s[k], tz);
286 dpoly fact(no_param, tz, 1);
287 D *= fact;
290 if (no_param + only_param == len) {
291 vec_QQ q;
292 q.SetLength(num_s.length());
293 for (int i = 0; i < num_s.length(); ++i) {
294 mpq_set_si(tcount, 0, 1);
295 n[i]->div(D, tcount, 1);
297 value2zz(mpq_numref(tcount), q[i].n);
298 value2zz(mpq_denref(tcount), q[i].d);
299 q[i] *= c2[i];
301 for (int i = q.length()-1; i >= 0; --i) {
302 if (q[i].n == 0) {
303 q[i] = q[q.length()-1];
304 num_p[i] = num_p[q.length()-1];
305 q.SetLength(q.length()-1);
306 num_p.SetDims(num_p.NumRows()-1, num_p.NumCols());
310 if (q.length() != 0)
311 reduce(q, num_p, pden);
312 } else {
313 value_set_si(tz, 0);
314 dpoly one(no_param, tz);
315 dpoly_r *r = NULL;
317 for (k = 0; k < len; ++k) {
318 if (den_s[k] == 0 || den_p[k] == 0)
319 continue;
321 zz2value(den_s[k], tz);
322 dpoly pd(no_param-1, tz, 1);
324 int l;
325 for (l = 0; l < k; ++l)
326 if (den_r[l] == den_r[k])
327 break;
329 if (!r)
330 r = new dpoly_r(one, pd, l, len);
331 else {
332 dpoly_r *nr = new dpoly_r(r, pd, l, len);
333 delete r;
334 r = nr;
338 vec_QQ factor;
339 factor.SetLength(c2.length());
340 int common = pden.NumRows();
341 dpoly_r *rc = r->div(D);
342 for (int i = 0; i < num_s.length(); ++i) {
343 factor[i].d = c2[i].d;
344 factor[i].d *= rc->denom;
347 dpoly_r_scanner scanner(n, num_s.length(), rc, len);
348 int rows;
349 while (scanner.next()) {
350 int i;
351 for (i = 0; i < num_s.length(); ++i)
352 if (scanner.coeff[i] != 0)
353 break;
354 if (i == num_s.length())
355 continue;
356 rows = common;
357 pden.SetDims(rows, pden.NumCols());
358 for (int k = 0; k < rc->dim; ++k) {
359 int n = scanner.powers[k];
360 if (n == 0)
361 continue;
362 pden.SetDims(rows+n, pden.NumCols());
363 for (int l = 0; l < n; ++l)
364 pden[rows+l] = den_r[k];
365 rows += n;
367 /* The denominators in factor are kept constant
368 * over all iterations of the enclosing while loop.
369 * The rational numbers in factor may therefore not be
370 * canonicalized. Some may even be zero.
372 for (int i = 0; i < num_s.length(); ++i) {
373 factor[i].n = c2[i].n;
374 factor[i].n *= scanner.coeff[i];
376 reduce_canonical(factor, num_p, pden);
379 delete rc;
380 delete r;
382 for (int i = 0; i < num_s.length(); ++i)
383 delete n[i];
384 delete [] n;
388 void reducer::handle(const mat_ZZ& den, Value *V, const QQ& c,
389 unsigned long det, barvinok_options *options)
391 vec_QQ vc;
393 Matrix *points = Matrix_Alloc(det, dim);
394 Matrix* Rays = zz2matrix(den);
395 lattice_points_fixed(V, V, Rays, Rays, points, det);
396 Matrix_Free(Rays);
397 matrix2zz(points, vertex, points->NbRows, points->NbColumns);
398 Matrix_Free(points);
400 vc.SetLength(vertex.NumRows());
401 for (int i = 0; i < vc.length(); ++i)
402 vc[i] = c;
404 reduce(vc, vertex, den);
407 void split_one(const mat_ZZ& num, vec_ZZ& num_s, mat_ZZ& num_p,
408 const mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r)
410 unsigned len = den_f.NumRows(); // number of factors in den
411 unsigned d = num.NumCols() - 1;
413 den_s.SetLength(len);
414 den_r.SetDims(len, d);
416 for (int r = 0; r < len; ++r) {
417 den_s[r] = den_f[r][0];
418 for (int k = 1; k <= d; ++k)
419 den_r[r][k-1] = den_f[r][k];
422 num_s.SetLength(num.NumRows());
423 num_p.SetDims(num.NumRows(), d);
424 for (int i = 0; i < num.NumRows(); ++i) {
425 num_s[i] = num[i][0];
426 for (int k = 1 ; k <= d; ++k)
427 num_p[i][k-1] = num[i][k];
431 void icounter::base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
433 zz2value(c.n, tn);
434 zz2value(c.d, td);
436 unsigned len = den_f.NumRows(); // number of factors in den
438 if (len > 0) {
439 int r;
440 vec_ZZ den_s;
441 den_s.SetLength(len);
442 assert(num.length() == 1);
443 ZZ num_s = num[0];
444 for (r = 0; r < len; ++r)
445 den_s[r] = den_f[r][0];
446 int sign = (len % 2) ? -1 : 1;
448 zz2value(num_s, tz);
449 dpoly n(len, tz);
450 zz2value(den_s[0], tz);
451 dpoly D(len, tz, 1);
452 for (int k = 1; k < len; ++k) {
453 zz2value(den_s[k], tz);
454 dpoly fact(len, tz, 1);
455 D *= fact;
457 mpq_set_si(tcount, 0, 1);
458 n.div(D, tcount, 1);
459 if (sign == -1)
460 value_oppose(tn, tn);
462 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
463 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
464 mpq_canonicalize(tcount);
465 } else {
466 value_assign(mpq_numref(tcount), tn);
467 value_assign(mpq_denref(tcount), td);