evalue.c: export evalue_denom
[barvinok.git] / reducer.cc
blob646ad9f2829e10629040da4bb04f84f7d93c5bea
1 #include <vector>
2 #include <barvinok/util.h>
3 #include "reducer.h"
4 #include "lattice_point.h"
6 using std::vector;
8 void np_base::handle_polar(Polyhedron *C, int s)
10 assert(C->NbRays-1 == dim);
11 factor.n *= s;
12 handle_polar(C, current_vertex, factor);
13 factor.n *= s;
16 void np_base::start(Polyhedron *P, unsigned MaxRays)
18 QQ factor(1, 1);
19 init(P);
20 for (int i = 0; i < P->NbRays; ++i) {
21 if (!value_pos_p(P->Ray[i][dim+1]))
22 continue;
24 Polyhedron *C = supporting_cone(P, i);
25 do_vertex_cone(factor, C, P->Ray[i]+1, MaxRays);
29 /* input:
30 * f: the powers in the denominator for the remaining vars
31 * each row refers to a factor
32 * den_s: for each factor, the power of (s+1)
33 * sign
34 * num_s: powers in the numerator corresponding to the summed vars
35 * num_p: powers in the numerator corresponding to the remaining vars
36 * number of rays in cone: "dim" = "k"
37 * length of each ray: "dim" = "d"
38 * for now, it is assumed: k == d
39 * output:
40 * den_p: for each factor
41 * 0: independent of remaining vars
42 * 1: power corresponds to corresponding row in f
44 * all inputs are subject to change
46 void normalize(ZZ& sign, ZZ& num_s, vec_ZZ& num_p, vec_ZZ& den_s, vec_ZZ& den_p,
47 mat_ZZ& f)
49 unsigned dim = f.NumRows();
50 unsigned nparam = num_p.length();
51 unsigned nvar = dim - nparam;
53 int change = 0;
55 for (int j = 0; j < den_s.length(); ++j) {
56 if (den_s[j] == 0) {
57 den_p[j] = 1;
58 continue;
60 int k;
61 for (k = 0; k < nparam; ++k)
62 if (f[j][k] != 0)
63 break;
64 if (k < nparam) {
65 den_p[j] = 1;
66 if (den_s[j] > 0) {
67 f[j] = -f[j];
68 num_p += f[j];
70 } else
71 den_p[j] = 0;
72 if (den_s[j] > 0)
73 change ^= 1;
74 else {
75 den_s[j] = abs(den_s[j]);
76 num_s += den_s[j];
80 if (change)
81 sign = -sign;
84 void reducer::reduce(QQ c, vec_ZZ& num, mat_ZZ& den_f)
86 unsigned len = den_f.NumRows(); // number of factors in den
88 if (num.length() == lower) {
89 base(c, num, den_f);
90 return;
92 assert(num.length() > 1);
94 vec_ZZ den_s;
95 mat_ZZ den_r;
96 ZZ num_s;
97 vec_ZZ num_p;
99 split(num, num_s, num_p, den_f, den_s, den_r);
101 vec_ZZ den_p;
102 den_p.SetLength(len);
104 normalize(c.n, num_s, num_p, den_s, den_p, den_r);
106 int only_param = 0; // k-r-s from text
107 int no_param = 0; // r from text
108 for (int k = 0; k < len; ++k) {
109 if (den_p[k] == 0)
110 ++no_param;
111 else if (den_s[k] == 0)
112 ++only_param;
114 if (no_param == 0) {
115 reduce(c, num_p, den_r);
116 } else {
117 int k, l;
118 mat_ZZ pden;
119 pden.SetDims(only_param, den_r.NumCols());
121 for (k = 0, l = 0; k < len; ++k)
122 if (den_s[k] == 0)
123 pden[l++] = den_r[k];
125 for (k = 0; k < len; ++k)
126 if (den_p[k] == 0)
127 break;
129 dpoly n(no_param, num_s);
130 dpoly D(no_param, den_s[k], 1);
131 for ( ; ++k < len; )
132 if (den_p[k] == 0) {
133 dpoly fact(no_param, den_s[k], 1);
134 D *= fact;
137 if (no_param + only_param == len) {
138 mpq_set_si(tcount, 0, 1);
139 n.div(D, tcount, one);
141 QQ q;
142 value2zz(mpq_numref(tcount), q.n);
143 value2zz(mpq_denref(tcount), q.d);
145 q *= c;
147 if (q.n != 0)
148 reduce(q, num_p, pden);
149 } else {
150 dpoly_r * r = 0;
152 for (k = 0; k < len; ++k) {
153 if (den_s[k] == 0 || den_p[k] == 0)
154 continue;
156 dpoly pd(no_param-1, den_s[k], 1);
158 int l;
159 for (l = 0; l < k; ++l)
160 if (den_r[l] == den_r[k])
161 break;
163 if (r == 0)
164 r = new dpoly_r(n, pd, l, len);
165 else {
166 dpoly_r *nr = new dpoly_r(r, pd, l, len);
167 delete r;
168 r = nr;
172 dpoly_r *rc = r->div(D);
174 QQ factor;
175 factor.d = rc->denom * c.d;
177 int common = pden.NumRows();
178 vector< dpoly_r_term * >& final = rc->c[rc->len-1];
179 int rows;
180 for (int j = 0; j < final.size(); ++j) {
181 if (final[j]->coeff == 0)
182 continue;
183 rows = common;
184 pden.SetDims(rows, pden.NumCols());
185 for (int k = 0; k < rc->dim; ++k) {
186 int n = final[j]->powers[k];
187 if (n == 0)
188 continue;
189 pden.SetDims(rows+n, pden.NumCols());
190 for (int l = 0; l < n; ++l)
191 pden[rows+l] = den_r[k];
192 rows += n;
194 factor.n = final[j]->coeff *= c.n;
195 reduce(factor, num_p, pden);
198 delete rc;
199 delete r;
204 void reducer::handle_polar(Polyhedron *C, Value *V, QQ c)
206 lattice_point(V, C, vertex);
208 mat_ZZ den;
209 den.SetDims(dim, dim);
211 int r;
212 for (r = 0; r < dim; ++r)
213 values2zz(C->Ray[r]+1, den[r], dim);
215 reduce(c, vertex, den);
218 void ireducer::split(vec_ZZ& num, ZZ& num_s, vec_ZZ& num_p,
219 mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r)
221 unsigned len = den_f.NumRows(); // number of factors in den
222 unsigned d = num.length() - 1;
224 den_s.SetLength(len);
225 den_r.SetDims(len, d);
227 for (int r = 0; r < len; ++r) {
228 den_s[r] = den_f[r][0];
229 for (int k = 1; k <= d; ++k)
230 den_r[r][k-1] = den_f[r][k];
233 num_s = num[0];
234 num_p.SetLength(d);
235 for (int k = 1 ; k <= d; ++k)
236 num_p[k-1] = num[k];
239 void normalize(ZZ& sign, ZZ& num, vec_ZZ& den)
241 unsigned dim = den.length();
243 int change = 0;
245 for (int j = 0; j < den.length(); ++j) {
246 if (den[j] > 0)
247 change ^= 1;
248 else {
249 den[j] = abs(den[j]);
250 num += den[j];
253 if (change)
254 sign = -sign;
257 void icounter::base(QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
259 int r;
260 unsigned len = den_f.NumRows(); // number of factors in den
261 vec_ZZ den_s;
262 den_s.SetLength(len);
263 ZZ num_s = num[0];
264 for (r = 0; r < len; ++r)
265 den_s[r] = den_f[r][0];
266 normalize(c.n, num_s, den_s);
268 dpoly n(len, num_s);
269 dpoly D(len, den_s[0], 1);
270 for (int k = 1; k < len; ++k) {
271 dpoly fact(len, den_s[k], 1);
272 D *= fact;
274 mpq_set_si(tcount, 0, 1);
275 n.div(D, tcount, one);
276 zz2value(c.n, tn);
277 zz2value(c.d, td);
278 mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
279 mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
280 mpq_canonicalize(tcount);
281 mpq_add(count, count, tcount);
284 void infinite_icounter::base(QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
286 int r;
287 unsigned len = den_f.NumRows(); // number of factors in den
288 vec_ZZ den_s;
289 den_s.SetLength(len);
290 ZZ num_s = num[0];
292 for (r = 0; r < len; ++r)
293 den_s[r] = den_f[r][0];
294 normalize(c.n, num_s, den_s);
296 dpoly n(len, num_s);
297 dpoly D(len, den_s[0], 1);
298 for (int k = 1; k < len; ++k) {
299 dpoly fact(len, den_s[k], 1);
300 D *= fact;
303 Value tmp;
304 mpq_t factor;
305 mpq_init(factor);
306 value_init(tmp);
307 zz2value(c.n, tmp);
308 value_assign(mpq_numref(factor), tmp);
309 zz2value(c.d, tmp);
310 value_assign(mpq_denref(factor), tmp);
312 n.div(D, count, factor);
314 value_clear(tmp);
315 mpq_clear(factor);