update AUTHORS
[barvinok.git] / laurent.cc
blob48021e3a4ab322113673a2b92c3aa64b07e13a4f
1 #include <assert.h>
2 #include <iostream>
3 #include <vector>
4 #include <barvinok/options.h>
5 #include <barvinok/util.h>
6 #include "binomial.h"
7 #include "decomposer.h"
8 #include "laurent.h"
9 #include "param_polynomial.h"
10 #include "param_util.h"
11 #include "reduce_domain.h"
12 #include "vertex_cone.h"
14 #define ALLOC(type) (type*)malloc(sizeof(type))
16 void multinomial(const std::vector<int> &k, Value *m);
18 using std::cerr;
19 using std::endl;
20 using std::ostream;
21 using std::vector;
23 template <typename T>
24 static ostream & operator<< (ostream & os, const vector<T> & v)
26 os << "[";
27 for (int i = 0; i < v.size(); ++i) {
28 if (i)
29 os << ", ";
30 os << v[i];
32 os << "]";
33 return os;
36 struct laurent_summator : public signed_cone_consumer,
37 public vertex_decomposer {
38 const evalue *polynomial;
39 unsigned dim;
40 vertex_cone vc;
41 param_polynomial poly;
42 evalue *result;
43 unsigned max_power;
45 /* For each row, first column with non-zero element */
46 vector<int> first;
47 /* For each row, last column with non-zero element */
48 vector<int> last;
49 /* For each column, number of rows that have this column as first */
50 vector<int> n_first;
51 /* For each column, last row that has this column as first */
52 vector<int> last_first;
53 /* For each column i,
54 * sum_{j >= i} power[j] + n_first[j] - [ sum_k selected[k][j] ]_+
56 vector<int> acc_power;
57 /* The exponents that we still need to match by subsequent factors */
58 vector<int> left;
59 /* For each factor,variable, the minimum allowed exponent */
60 vector< vector<int> > min;
61 /* For each factor,variable, the maximum allowed exponent */
62 vector< vector<int> > max;
63 /* For each factor,variable, the currently selected exponent */
64 vector< vector<int> > selected;
66 laurent_summator(const evalue *e, unsigned dim,
67 Param_Polyhedron *PP) :
68 vertex_decomposer(PP, *this), polynomial(e), dim(dim),
69 vc(dim), poly(e, dim) {
70 max_power = dim + poly.degree();
71 result = NULL;
72 first = vector<int>(dim);
73 last = vector<int>(dim);
74 n_first = vector<int>(dim);
75 last_first = vector<int>(dim);
76 acc_power = vector<int>(dim);
77 for (int i = 0; i < dim; ++i) {
78 min.push_back(vector<int>(dim));
79 max.push_back(vector<int>(dim));
80 selected.push_back(vector<int>(dim));
83 ~laurent_summator() {
84 if (result)
85 evalue_free(result);
87 virtual void handle(const signed_cone& sc, barvinok_options *options);
89 evalue *handle_term(vector<int> &power);
90 void set_min_max(int row, int col);
91 evalue *selection_contribution();
94 /* Change (row,col) to the previous entry in a dim x dim matrix */
95 static void prev(int &row, int &col, int dim, int &init)
97 init = 0;
98 if (--col >= 0)
99 return;
100 --row;
101 col = dim - 1;
104 /* Change (row,col) to the next entry in a dim x dim matrix */
105 static void next(int &row, int &col, int dim, int &init)
107 init = 1;
108 if (++col < dim)
109 return;
110 ++row;
111 col = 0;
114 /* Set min[row][col] and max[row][col] to the minimum and maximum
115 * possible exponents for variable col in factor row, given a choice
116 * of exponents for the previous factors and the previous variables in
117 * the same factor.
119 * There are two kinds of terms in each factor i.
120 * In the first kind, the exponent a_ij of first[i] is negative
121 * and equal to -1 - sum_{k > j} a_ik
122 * In the second kind, all exponents are non-negative.
124 * If the coefficient of the corresponding ray is zero, in particular
125 * if col < first[row], then we can only construct terms that are
126 * constant in the given variable, i.e., min = max = 0.
128 * If col == first[row], then we can assign both negative and
129 * positive exponents to this variable.
130 * If we assign a negative exponent a_ij, then we need to make
131 * sure that we can distribute -a_ij - 1 over the remaining variables.
132 * The total sum of exponents we can distribute is equal to
133 * the total sum in the target exponents for the subsequent variables
134 * + the total number of times any of the subsequent variables is
135 * the first variable with non-zero coefficient in a ray
136 * (since we can increase the total exponent by one,
137 * by making the exponent negative)
138 * - the total sum of (positive) exponents we have already selected
139 * for these variables in earlier factors.
140 * This total sum is maintained in acc_power[col + 1].
141 * However, if this column is not only the first, but also the last (and only)
142 * non-zero column, then the only negative exponent we can allocate is -1.
144 * If col > first[row], then we can only assign non-negative exponents,
145 * so we set min[row][col] to 0.
147 * The maximal exponent is the target exponent minus the exponents we
148 * have already selected in previous factors.
149 * However, if the given column appears first in at least one row,
150 * then we can move exponents of later columns to this column by selecting
151 * a negative exponent in later factors, and then the upper bound is
152 * given by acc_power[col]. If the current factor corresponds to one
153 * of these rows, then we can subtract 1 from acc_power[col], because
154 * if the current exponent is non-negative, then it cannot be used
155 * to increase the total exponent by 1.
157 * If the current column is not the first non-zero column and
158 * a negative exponent was chosen for the first non-zero column,
159 * then we need to make sure that we don't exceed
160 * -selected[row][first[row]] - 1 in this row. The maximal exponent
161 * is adjusted accordingly.
162 * Furthermore, if this is the last column in such a row, then we
163 * have to select what is left and adjust the minimal exponent accordingly.
165 * Finally, if this is the last row, then we have to select what is left.
166 * The minimum and maximum exponent are adjusted accordingly.
168 void laurent_summator::set_min_max(int row, int col)
170 if (col == first[row]) {
171 if (col < last[row])
172 min[row][col] = -acc_power[col + 1] - 1;
173 else
174 min[row][col] = -1;
175 max[row][col] = acc_power[col] - 1;
177 if (col < first[row]) {
178 min[row][col] = 0;
179 max[row][col] = 0;
181 if (col > first[row]) {
182 min[row][col] = 0;
183 max[row][col] = n_first[col] > 0 ? acc_power[col] : left[col];
185 if (value_zero_p(vc.Rays.p[row][col])) {
186 if (max[row][col] > 0)
187 max[row][col] = 0;
188 if (min[row][col] < 0)
189 min[row][col] = 0;
191 if (col > first[row] && selected[row][first[row]] < 0) {
192 int l = -selected[row][first[row]] - 1;
193 for (int i = first[row] + 1; i < col; ++i)
194 l -= selected[row][i];
195 if (l < max[row][col])
196 max[row][col] = l;
197 if (col >= last[row] && l > min[row][col])
198 min[row][col] = l;
200 if (row == dim - 1 || (col == first[row] && last_first[col] == row)) {
201 if (left[col] < max[row][col])
202 max[row][col] = left[col];
203 if (left[col] > min[row][col])
204 min[row][col] = left[col];
208 /* Compute and return the contribution of the selected distribution of
209 * exponents over the factors. The total contribution is the product
210 * of the contribution of each factor.
211 * The contributions of a factor depends on whether the first exponent
212 * is negative or not.
214 * If the first exponent is negative, then it is equal to
215 * -m = - 1 - \sum_k n_k, with n_k the remaining exponents.
216 * The contribution of this factor is then
218 * (m-1 \choose n) (-1)^m b_j^n b_{jf}^{-m}
220 * Otherwise, i.e., when all exponents are non-negative, let
221 * m = 1 + \sum_k n_k, with n_k all exponents.
222 * Then the contribution of this factor is
224 * bernoulli(m, s_j)/(m * \prod_k n_k!) b_j^n
226 * In these formulae, b_j represents the coefficients of ray j
227 * and s_j is such that the exponent of the numerator of the
228 * vertex cone is w = <s_j, b_j>.
230 evalue *laurent_summator::selection_contribution()
232 evalue *c = NULL;
234 for (int i = 0; i < dim; ++i) {
235 evalue *f = ALLOC(evalue);
236 value_init(f->d);
237 evalue_set_si(f, 1, 1);
238 if (selected[i][first[i]] < 0) {
239 int d_exp = -selected[i][first[i]];
240 selected[i][first[i]] = 0;
241 multinomial(selected[i], &f->x.n);
242 selected[i][first[i]] = -d_exp;
243 if (d_exp % 2)
244 value_oppose(f->x.n, f->x.n);
245 for (int j = first[i] + 1; j < dim; ++j) {
246 if (selected[i][j] == 0)
247 continue;
248 value_multiply(f->x.n, f->x.n,
249 *(*vc.coeff_power[i][j])[selected[i][j]]);
251 value_multiply(f->d, f->d,
252 *(*vc.coeff_power[i][first[i]])[d_exp]);
253 if (value_neg_p(f->d)) {
254 value_oppose(f->d, f->d);
255 value_oppose(f->x.n, f->x.n);
257 } else {
258 int sum = 0;
259 for (int j = 0; j < dim; ++j)
260 sum += selected[i][j];
261 value_set_si(f->x.n, -1);
262 value_set_si(f->d, 1 + sum);
263 for (int j = 0; j < dim; ++j) {
264 if (selected[i][j] == 0)
265 continue;
266 value_multiply(f->x.n, f->x.n,
267 *(*vc.coeff_power[i][j])[selected[i][j]]);
268 value_multiply(f->d, f->d,
269 *factorial(selected[i][j]));
271 emul(vc.get_bernoulli(1 + sum, i), f);
273 if (!c)
274 c = f;
275 else {
276 emul(f, c);
277 evalue_free(f);
280 return c;
283 /* Compute and return the coefficient of the term with exponents "power"
284 * in the Laurent expansion.
286 * We perform a depth-first search over all distributions of "power"
287 * over the factors and collect (sum) all the corresponding coefficients.
288 * Every time we enter a node for the first time, we compute the minimum
289 * and maximum possible exponent and select the minimum.
290 * Every time we return to the node during backtracking, we increase
291 * the exponent by 1 until we have reached the maximum.
293 evalue *laurent_summator::handle_term(vector<int> &power)
295 evalue *s = NULL;
297 left = power;
298 for (int i = 0; i < dim; ++i)
299 n_first[i] = 0;
300 for (int i = 0; i < dim; ++i) {
301 first[i] = First_Non_Zero(vc.Rays.p[i], dim);
302 last[i] = Last_Non_Zero(vc.Rays.p[i], dim);
303 n_first[first[i]]++;
304 last_first[first[i]] = i;
306 acc_power[dim - 1] = power[dim - 1] + n_first[dim - 1];
307 for (int i = dim - 2; i >= 0; --i)
308 acc_power[i] = power[i] + n_first[i] + acc_power[i + 1];
310 int row = 0, col = 0, init = 1;
312 while (row >= 0) {
313 if (row >= dim) {
314 evalue *t;
316 t = selection_contribution();
317 if (!s)
318 s = t;
319 else {
320 eadd(t, s);
321 evalue_free(t);
323 prev(row, col, dim, init);
325 if (init) {
326 set_min_max(row, col);
327 if (max[row][col] < min[row][col]) {
328 prev(row, col, dim, init);
329 continue;
331 selected[row][col] = min[row][col];
332 left[col] -= selected[row][col];
333 if (selected[row][col] >= 0)
334 acc_power[col] -= selected[row][col];
335 } else {
336 if (selected[row][col] >= max[row][col]) {
337 left[col] += selected[row][col];
338 if (selected[row][col] >= 0)
339 acc_power[col] += selected[row][col];
340 prev(row, col, dim, init);
341 continue;
343 if (selected[row][col] >= 0)
344 acc_power[col]--;
345 selected[row][col]++;
346 left[col]--;
348 next(row, col, dim, init);
351 return s;
354 /* Compute and accumulate the contribution of the given vertex cone
355 * to the total sum.
357 * We compute the corresponding coefficient in the Laurent expansion
358 * and divide it by n!, with n the vector of exponents.
360 void laurent_summator::handle(const signed_cone& sc, barvinok_options *options)
362 assert(sc.det == 1);
364 vc.init(sc.rays, V, max_power);
365 for (int i = 0; i < poly.terms.size(); ++i) {
366 evalue *t, *f;
368 t = handle_term(poly.terms[i].powers);
370 f = evalue_dup(poly.terms[i].coeff);
371 if (sc.sign < 0)
372 evalue_negate(f);
373 for (int j = 0; j < dim; ++j)
374 evalue_mul(f, *factorial(poly.terms[i].powers[j]));
375 evalue_shift_variables(f, 0, -dim);
376 emul(f, t);
377 evalue_free(f);
378 if (!result)
379 result = t;
380 else {
381 eadd(t, result);
382 evalue_free(t);
385 vc.clear();
388 evalue *laurent_summate(Polyhedron *P, evalue *e, unsigned nvar,
389 struct barvinok_options *options)
391 Polyhedron *U, *TC;
392 Param_Polyhedron *PP;
393 struct evalue_section *s;
394 int nd = -1;
395 Param_Domain *PD;
396 evalue *sum;
398 U = Universe_Polyhedron(P->Dimension - nvar);
399 PP = Polyhedron2Param_Polyhedron(P, U, options);
401 for (nd = 0, PD = PP->D; PD; ++nd, PD = PD->next);
402 s = ALLOCN(struct evalue_section, nd);
404 TC = true_context(P, U, options->MaxRays);
405 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, PD, rVD)
406 Param_Vertices *V;
407 laurent_summator ls(e, nvar, PP);
409 FORALL_PVertex_in_ParamPolyhedron(V, PD, PP) // _i internal
410 ls.decompose_at_vertex(V, _i, options);
411 END_FORALL_PVertex_in_ParamPolyhedron;
413 s[i].D = rVD;
414 s[i].E = ls.result;
415 ls.result = NULL;
416 END_FORALL_REDUCED_DOMAIN
417 Polyhedron_Free(TC);
418 Polyhedron_Free(U);
419 Param_Polyhedron_Free(PP);
421 sum = evalue_from_section_array(s, nd);
422 free(s);
424 return sum;