1 #include <iostream>
2 #include <ostream>
3 #include <vector>
4 #include <NTL/ZZ.h>
5 #include <NTL/vec_ZZ.h>
6 #include <NTL/mat_ZZ.h>
7 #include <barvinok/polylib.h>
8 #include <barvinok/util.h>
9 #include "reducer.h"
10 #include "lattice_point.h"
11
12 using std::vector;
13 using std::cerr;
14 using std::endl;
15
16 struct OrthogonalException Orthogonal;
17
handle(const signed_cone & sc,barvinok_options * options)18 void np_base::handle(const signed_cone& sc, barvinok_options *options)
19 {
20 assert(sc.rays.NumRows() == dim);
21 factor.n *= sc.sign;
22 handle(sc.rays, current_vertex, factor, sc.det, options);
23 factor.n *= sc.sign;
24 }
25
start(Polyhedron * P,barvinok_options * options)26 void np_base::start(Polyhedron *P, barvinok_options *options)
27 {
28 int n_try = 0;
29 QQ factor(1, 1);
30 for (;;) {
31 try {
32 init(P, n_try);
33 for (int i = 0; i < P->NbRays; ++i) {
34 if (!value_pos_p(P->Ray[i][dim+1]))
35 continue;
36
37 Polyhedron *C = supporting_cone(P, i);
38 do_vertex_cone(factor, C, P->Ray[i]+1, options);
39 }
40 break;
41 } catch (OrthogonalException &e) {
42 n_try++;
43 reset();
44 }
45 }
46 }
47
48 /* input:
49 * f: the powers in the denominator for the remaining vars
50 * each row refers to a factor
51 * den_s: for each factor, the power of (s+1)
52 * sign
53 * num_s: powers in the numerator corresponding to the summed vars
54 * num_p: powers in the numerator corresponding to the remaining vars
55 * number of rays in cone: "dim" = "k"
56 * length of each ray: "dim" = "d"
57 * for now, it is assumed: k == d
58 * output:
59 * den_p: for each factor
60 * 0: independent of remaining vars
61 * 1: power corresponds to corresponding row in f
62 *
63 * all inputs are subject to change
64 */
normalize(ZZ & sign,vec_ZZ & num_s,mat_ZZ & num_p,vec_ZZ & den_s,vec_ZZ & den_p,mat_ZZ & f)65 void normalize(ZZ& sign, vec_ZZ& num_s, mat_ZZ& num_p, vec_ZZ& den_s, vec_ZZ& den_p,
66 mat_ZZ& f)
67 {
68 unsigned nparam = num_p.NumCols();
69 int change = 0;
70
71 for (int j = 0; j < den_s.length(); ++j) {
72 if (den_s[j] == 0) {
73 den_p[j] = 1;
74 continue;
75 }
76 int k;
77 for (k = 0; k < nparam; ++k)
78 if (f[j][k] != 0)
79 break;
80 if (k < nparam) {
81 den_p[j] = 1;
82 if (den_s[j] > 0) {
83 f[j] = -f[j];
84 for (int i = 0; i < num_p.NumRows(); ++i)
85 num_p[i] += f[j];
86 }
87 } else
88 den_p[j] = 0;
89 if (den_s[j] > 0)
90 change ^= 1;
91 else {
92 den_s[j] = abs(den_s[j]);
93 for (int i = 0; i < num_p.NumRows(); ++i)
94 num_s[i] += den_s[j];
95 }
96 }
97
98 if (change)
99 sign = -sign;
100 }
101
base(const vec_QQ & c,const mat_ZZ & num,const mat_ZZ & den_f)102 void reducer::base(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den_f)
103 {
104 for (int i = 0; i < num.NumRows(); ++i)
105 base(c[i], num[i], den_f);
106 }
107
108 struct dpoly_r_scanner {
109 const dpoly_r *rc;
110 const dpoly * const *num;
111 int n;
112 int dim;
113 dpoly_r_term_list::iterator *iter;
114 vector<int> powers;
115 vec_ZZ coeff;
116
dpoly_r_scannerdpoly_r_scanner117 dpoly_r_scanner(const dpoly * const *num, int n, const dpoly_r *rc, int dim)
118 : rc(rc), num(num), n(n), dim(dim), powers(dim, 0) {
119 coeff.SetLength(n);
120 iter = new dpoly_r_term_list::iterator[rc->len];
121 for (int i = 0; i < rc->len; ++i) {
122 int k;
123 for (k = 0; k < n; ++k)
124 if (value_notzero_p(num[k]->coeff->p[rc->len-1-i]))
125 break;
126 if (k < n)
127 iter[i] = rc->c[i].begin();
128 else
129 iter[i] = rc->c[i].end();
130 }
131 }
nextdpoly_r_scanner132 bool next() {
133 int *pos;
134 int len = 0;
135
136 for (int i = 0; i < rc->len; ++i) {
137 if (iter[i] == rc->c[i].end())
138 continue;
139 if (!len) {
140 pos = new int[rc->len];
141 pos[len++] = i;
142 } else {
143 if ((*iter[i])->powers < (*iter[pos[0]])->powers) {
144 pos[0] = i;
145 len = 1;
146 } else if ((*iter[i])->powers == (*iter[pos[0]])->powers)
147 pos[len++] = i;
148 }
149 }
150
151 if (!len)
152 return false;
153
154 powers = (*iter[pos[0]])->powers;
155 for (int k = 0; k < n; ++k) {
156 value2zz(num[k]->coeff->p[rc->len-1-pos[0]], tmp);
157 mul(coeff[k], (*iter[pos[0]])->coeff, tmp);
158 }
159 ++iter[pos[0]];
160 for (int i = 1; i < len; ++i) {
161 for (int k = 0; k < n; ++k) {
162 value2zz(num[k]->coeff->p[rc->len-1-pos[i]], tmp);
163 mul(tmp, (*iter[pos[i]])->coeff, tmp);
164 add(coeff[k], coeff[k], tmp);
165 }
166 ++iter[pos[i]];
167 }
168
169 delete [] pos;
170 return true;
171 }
~dpoly_r_scannerdpoly_r_scanner172 ~dpoly_r_scanner() {
173 delete [] iter;
174 }
175 private:
176 ZZ tmp;
177 };
178
reduce_canonical(const vec_QQ & c,const mat_ZZ & num,const mat_ZZ & den_f)179 void reducer::reduce_canonical(const vec_QQ& c, const mat_ZZ& num,
180 const mat_ZZ& den_f)
181 {
182 vec_QQ c2 = c;
183 mat_ZZ num2 = num;
184
185 for (int i = 0; i < c2.length(); ++i) {
186 c2[i].canonicalize();
187 if (c2[i].n != 0)
188 continue;
189
190 if (i < c2.length()-1) {
191 num2[i] = num2[c2.length()-1];
192 c2[i] = c2[c2.length()-1];
193 }
194 num2.SetDims(num2.NumRows()-1, num2.NumCols());
195 c2.SetLength(c2.length()-1);
196 --i;
197 }
198 reduce(c2, num2, den_f);
199 }
200
reduce(const vec_QQ & c,const mat_ZZ & num,const mat_ZZ & den_f)201 void reducer::reduce(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den_f)
202 {
203 assert(c.length() == num.NumRows());
204 unsigned len = den_f.NumRows(); // number of factors in den
205 vec_QQ c2 = c;
206
207 if (num.NumCols() == lower) {
208 base(c, num, den_f);
209 return;
210 }
211 assert(num.NumCols() > 1);
212 assert(num.NumRows() > 0);
213
214 vec_ZZ den_s;
215 mat_ZZ den_r;
216 vec_ZZ num_s;
217 mat_ZZ num_p;
218
219 split(num, num_s, num_p, den_f, den_s, den_r);
220
221 vec_ZZ den_p;
222 den_p.SetLength(len);
223
224 ZZ sign(INIT_VAL, 1);
225 normalize(sign, num_s, num_p, den_s, den_p, den_r);
226 c2 *= sign;
227
228 int only_param = 0; // k-r-s from text
229 int no_param = 0; // r from text
230 for (int k = 0; k < len; ++k) {
231 if (den_p[k] == 0)
232 ++no_param;
233 else if (den_s[k] == 0)
234 ++only_param;
235 }
236 if (no_param == 0) {
237 reduce(c2, num_p, den_r);
238 } else {
239 int k, l;
240 mat_ZZ pden;
241 pden.SetDims(only_param, den_r.NumCols());
242
243 for (k = 0, l = 0; k < len; ++k)
244 if (den_s[k] == 0)
245 pden[l++] = den_r[k];
246
247 for (k = 0; k < len; ++k)
248 if (den_p[k] == 0)
249 break;
250
251 dpoly **n = new dpoly *[num_s.length()];
252 for (int i = 0; i < num_s.length(); ++i) {
253 zz2value(num_s[i], tz);
254 n[i] = new dpoly(no_param, tz);
255 /* Search for other numerator (j) with same num_p.
256 * If found, replace a[j]/b[j] * n[j] and a[i]/b[i] * n[i]
257 * by 1/(b[j]*b[i]/g) * (a[j]*b[i]/g * n[j] + a[i]*b[j]/g * n[i])
258 * where g = gcd(b[i], b[j].
259 */
260 for (int j = 0; j < i; ++j) {
261 if (num_p[i] != num_p[j])
262 continue;
263 ZZ g = GCD(c2[i].d, c2[j].d);
264 zz2value(c2[j].n * c2[i].d/g, tz);
265 *n[j] *= tz;
266 zz2value(c2[i].n * c2[j].d/g, tz);
267 *n[i] *= tz;
268 *n[j] += *n[i];
269 c2[j].n = 1;
270 c2[j].d *= c2[i].d/g;
271 delete n[i];
272 if (i < num_s.length()-1) {
273 num_s[i] = num_s[num_s.length()-1];
274 num_p[i] = num_p[num_s.length()-1];
275 c2[i] = c2[num_s.length()-1];
276 }
277 num_s.SetLength(num_s.length()-1);
278 c2.SetLength(c2.length()-1);
279 num_p.SetDims(num_p.NumRows()-1, num_p.NumCols());
280 --i;
281 break;
282 }
283 }
284 zz2value(den_s[k], tz);
285 dpoly D(no_param, tz, 1);
286 for ( ; ++k < len; )
287 if (den_p[k] == 0) {
288 zz2value(den_s[k], tz);
289 dpoly fact(no_param, tz, 1);
290 D *= fact;
291 }
292
293 if (no_param + only_param == len) {
294 vec_QQ q;
295 q.SetLength(num_s.length());
296 for (int i = 0; i < num_s.length(); ++i) {
297 mpq_set_si(tcount, 0, 1);
298 n[i]->div(D, tcount, 1);
299
300 value2zz(mpq_numref(tcount), q[i].n);
301 value2zz(mpq_denref(tcount), q[i].d);
302 q[i] *= c2[i];
303 }
304 for (int i = q.length()-1; i >= 0; --i) {
305 if (q[i].n == 0) {
306 q[i] = q[q.length()-1];
307 num_p[i] = num_p[q.length()-1];
308 q.SetLength(q.length()-1);
309 num_p.SetDims(num_p.NumRows()-1, num_p.NumCols());
310 }
311 }
312
313 if (q.length() != 0)
314 reduce(q, num_p, pden);
315 } else {
316 value_set_si(tz, 0);
317 dpoly one(no_param, tz);
318 dpoly_r *r = NULL;
319
320 for (k = 0; k < len; ++k) {
321 if (den_s[k] == 0 || den_p[k] == 0)
322 continue;
323
324 zz2value(den_s[k], tz);
325 dpoly pd(no_param-1, tz, 1);
326
327 int l;
328 for (l = 0; l < k; ++l)
329 if (den_r[l] == den_r[k])
330 break;
331
332 if (!r)
333 r = new dpoly_r(one, pd, l, len);
334 else {
335 dpoly_r *nr = new dpoly_r(r, pd, l, len);
336 delete r;
337 r = nr;
338 }
339 }
340
341 vec_QQ factor;
342 factor.SetLength(c2.length());
343 int common = pden.NumRows();
344 dpoly_r *rc = r->div(D);
345 for (int i = 0; i < num_s.length(); ++i) {
346 factor[i].d = c2[i].d;
347 factor[i].d *= rc->denom;
348 }
349
350 dpoly_r_scanner scanner(n, num_s.length(), rc, len);
351 int rows;
352 while (scanner.next()) {
353 int i;
354 for (i = 0; i < num_s.length(); ++i)
355 if (scanner.coeff[i] != 0)
356 break;
357 if (i == num_s.length())
358 continue;
359 rows = common;
360 pden.SetDims(rows, pden.NumCols());
361 for (int k = 0; k < rc->dim; ++k) {
362 int n = scanner.powers[k];
363 if (n == 0)
364 continue;
365 pden.SetDims(rows+n, pden.NumCols());
366 for (int l = 0; l < n; ++l)
367 pden[rows+l] = den_r[k];
368 rows += n;
369 }
370 /* The denominators in factor are kept constant
371 * over all iterations of the enclosing while loop.
372 * The rational numbers in factor may therefore not be
373 * canonicalized. Some may even be zero.
374 */
375 for (int i = 0; i < num_s.length(); ++i) {
376 factor[i].n = c2[i].n;
377 factor[i].n *= scanner.coeff[i];
378 }
379 reduce_canonical(factor, num_p, pden);
380 }
381
382 delete rc;
383 delete r;
384 }
385 for (int i = 0; i < num_s.length(); ++i)
386 delete n[i];
387 delete [] n;
388 }
389 }
390
handle(const mat_ZZ & den,Value * V,const QQ & c,unsigned long det,barvinok_options * options)391 void reducer::handle(const mat_ZZ& den, Value *V, const QQ& c,
392 unsigned long det, barvinok_options *options)
393 {
394 vec_QQ vc;
395
396 Matrix *points = Matrix_Alloc(det, dim);
397 Matrix* Rays = zz2matrix(den);
398 lattice_points_fixed(V, V, Rays, Rays, points, det);
399 Matrix_Free(Rays);
400 matrix2zz(points, vertex, points->NbRows, points->NbColumns);
401 Matrix_Free(points);
402
403 vc.SetLength(vertex.NumRows());
404 for (int i = 0; i < vc.length(); ++i)
405 vc[i] = c;
406
407 reduce(vc, vertex, den);
408 }
409
split_one(const mat_ZZ & num,vec_ZZ & num_s,mat_ZZ & num_p,const mat_ZZ & den_f,vec_ZZ & den_s,mat_ZZ & den_r)410 void split_one(const mat_ZZ& num, vec_ZZ& num_s, mat_ZZ& num_p,
411 const mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r)
412 {
413 unsigned len = den_f.NumRows(); // number of factors in den
414 unsigned d = num.NumCols() - 1;
415
416 den_s.SetLength(len);
417 den_r.SetDims(len, d);
418
419 for (int r = 0; r < len; ++r) {
420 den_s[r] = den_f[r][0];
421 for (int k = 1; k <= d; ++k)
422 den_r[r][k-1] = den_f[r][k];
423 }
424
425 num_s.SetLength(num.NumRows());
426 num_p.SetDims(num.NumRows(), d);
427 for (int i = 0; i < num.NumRows(); ++i) {
428 num_s[i] = num[i][0];
429 for (int k = 1 ; k <= d; ++k)
430 num_p[i][k-1] = num[i][k];
431 }
432 }
433
base(const QQ & c,const vec_ZZ & num,const mat_ZZ & den_f)434 void icounter::base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
435 {
436 zz2value(c.n, tn);
437 zz2value(c.d, td);
438
439 unsigned len = den_f.NumRows(); // number of factors in den
440
441 if (len > 0) {
442 int r;
443 vec_ZZ den_s;
444 den_s.SetLength(len);
445 assert(num.length() == 1);
446 ZZ num_s = num[0];
447 for (r = 0; r < len; ++r)
448 den_s[r] = den_f[r][0];
449 int sign = (len % 2) ? -1 : 1;
450
451 zz2value(num_s, tz);
452 dpoly n(len, tz);
453 zz2value(den_s[0], tz);
454 dpoly D(len, tz, 1);
455 for (int k = 1; k < len; ++k) {
456 zz2value(den_s[k], tz);
457 dpoly fact(len, tz, 1);
458 D *= fact;
459 }
460 mpq_set_si(tcount, 0, 1);
461 n.div(D, tcount, 1);
462 if (sign == -1)
463 value_oppose(tn, tn);
464
465 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
466 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
467 mpq_canonicalize(tcount);
468 } else {
469 value_assign(mpq_numref(tcount), tn);
470 value_assign(mpq_denref(tcount), td);
471 }
472 mpq_add(count, count, tcount);
473 }
474