iscc: only print valid results
[barvinok.git] / laurent_old.cc
blob51585c81bc101b43dc91c503d42e5c95664a0291
1 #include <assert.h>
2 #include <iostream>
3 #include <vector>
4 #include <barvinok/options.h>
5 #include "binomial.h"
6 #include "conversion.h"
7 #include "decomposer.h"
8 #include "laurent_old.h"
9 #include "param_polynomial.h"
10 #include "param_util.h"
11 #include "reduce_domain.h"
12 #include "vertex_cone.h"
13 #include "config.h"
15 using std::cerr;
16 using std::ostream;
17 using std::endl;
18 using std::vector;
20 #ifdef HAVE_GNUCXX_HASHMAP
22 #include <ext/hash_map>
24 #define HASH_MAP __gnu_cxx::hash_map
26 namespace __gnu_cxx
28 template<> struct hash< std::vector<int> >
30 size_t operator()( const std::vector<int>& x ) const
32 unsigned long __h = 0;
33 for (int i = 0; i < x.size(); ++i)
34 __h = 5 * __h + x[i];
35 return size_t(__h);
40 #else
42 #warning "no hash_map available"
43 #include <map>
44 #define HASH_MAP std::map
46 #endif
48 #define ALLOC(type) (type*)malloc(sizeof(type))
49 #define ALLOCN(type,n) (type*)malloc((n) * sizeof(type))
51 static ostream & operator<< (ostream & os, const vector<int> & v)
53 os << "[";
54 for (int i = 0; i < v.size(); ++i) {
55 if (i)
56 os << ", ";
57 os << v[i];
59 os << "]";
60 return os;
63 struct todd_product {
64 vertex_cone &vc;
65 unsigned dim;
66 /* The non-zero coefficients in the rays of the vertex cone */
67 vector< vector<int> > mask;
68 /* For each ray, the power of each variable it contributes */
69 vector< vector<int> > selected;
70 /* The powers of each variable that still need to be selected */
71 vector<int> left;
72 /* For each variable, the last ray that has a non-zero coefficient */
73 vector<int> last_level;
75 HASH_MAP<vector<int>, const evalue *> cache;
77 todd_product(vertex_cone &vc);
78 evalue *add(evalue *sum);
79 const evalue *get_coefficient(const vector<int> &powers);
81 ~todd_product() {
82 HASH_MAP<vector<int>, const evalue *>::iterator j;
83 for (j = cache.begin(); j != cache.end(); ++j)
84 if ((*j).second)
85 evalue_free(const_cast<evalue *>((*j).second));
89 todd_product::todd_product(vertex_cone &vc) : vc(vc)
91 dim = vc.dim;
92 for (int i = 0; i < dim; ++i) {
93 mask.push_back(vector<int>(dim));
94 selected.push_back(vector<int>(dim));
96 last_level = vector<int>(dim);
97 for (int i = 0; i < dim; ++i) {
98 for (int j = 0; j < dim; ++j) {
99 if (vc.coeff_power[i][j]) {
100 mask[i][j] = 1;
101 last_level[j] = i;
107 void multinomial(const vector<int> &k, Value *m)
109 int s = 0;
110 value_set_si(*m, 1);
111 for (int i = 0; i < k.size(); ++i) {
112 if (k[i] == 0)
113 continue;
114 s += k[i];
115 value_multiply(*m, *m, *binomial(s, k[i]));
119 /* Add coefficient of selected powers of variables to sum
120 * and return the result.
121 * The contribution of each ray j of the vertex cone is
123 * ( \sum k )
124 * b(\sum k, \ceil{v}) ( k_1, \ldots, k_n ) c_1^{k_1} \cdots c_n^{k_n}
126 * with k_i the selected powers, c_i the coefficients of the ray
127 * and \ceil{v} the coordinate E_vertex[j] corresponding to the ray.
129 evalue *todd_product::add(evalue *sum)
131 evalue *c = NULL;
132 for (int i = 0; i < dim; ++i) {
133 int s = 0;
134 evalue *f = ALLOC(evalue);
135 value_init(f->d);
136 evalue_set_si(f, 1, 1);
137 for (int j = 0; j < dim; ++j) {
138 if (!selected[i][j])
139 continue;
140 value_multiply(f->x.n, f->x.n,
141 *(*vc.coeff_power[i][j])[selected[i][j]]);
142 value_multiply(f->d, f->d, *factorial(selected[i][j]));
143 s += selected[i][j];
145 if (s > 0)
146 emul(vc.get_bernoulli(s, i), f);
147 if (!c)
148 c = f;
149 else {
150 emul(f, c);
151 evalue_free(f);
154 if (!sum)
155 sum = c;
156 else {
157 eadd(c, sum);
158 evalue_free(c);
160 return sum;
163 static int first_non_zero(const vector<int>& row)
165 for (int i = 0; i < row.size(); ++i)
166 if (row[i] != 0)
167 return i;
168 return -1;
171 /* Return coefficient of the term with exponents powers by
172 * iterating over all combinations of exponents for each ray
173 * that sum up to the given exponents.
175 const evalue *todd_product::get_coefficient(const vector<int> &powers)
177 evalue *c = NULL;
179 HASH_MAP<vector<int>, const evalue *>::iterator i;
180 i = cache.find(powers);
181 if (i != cache.end())
182 return (*i).second;
184 for (int i = 0; i < vc.dim; ++i)
185 for (int j = 0; j < vc.dim; ++j)
186 selected[i][j] = 0;
188 left = powers;
189 int nz = first_non_zero(left);
190 int level = last_level[nz];
191 int p = nz;
192 while (level >= 0) {
193 if (mask[level][p] && left[p] > 0) {
194 selected[level][p]++;
195 left[p]--;
196 /* Fill out remaining powers and make sure we backtrack from
197 * the right position.
199 for (int i = 0; i < vc.dim; ++i) {
200 if (left[i] == 0)
201 continue;
202 selected[last_level[i]][i] += left[i];
203 left[i] = 0;
204 if (last_level[i] >= level) {
205 level = last_level[i];
206 p = i;
210 c = add(c);
212 if (selected[level][p]) {
213 left[p] += selected[level][p];
214 selected[level][p] = 0;
216 if (--p < 0) {
217 --level;
218 p = dim-1;
221 cache[powers] = c;
222 return c;
225 /* Maintains the coefficients of the reciprocals of the
226 * (negated) rays of the vertex cone vc.
228 struct reciprocal {
229 vertex_cone &vc;
231 /* can_borrow_from[i][j] = 1 if there is a ray
232 * with first non-zero coefficient i and a subsequent
233 * non-zero coefficient j.
235 vector< vector<int> > can_borrow_from;
236 /* Total exponent that a variable can borrow from subsequent vars */
237 vector<int> can_borrow;
238 /* has_borrowed_from[i][j] is the exponent borrowed by i from j */
239 vector< vector<int> > has_borrowed_from;
240 /* Total exponent borrowed from subsequent vars */
241 vector<int> has_borrowed;
242 /* The last variable that can borrow from subsequent vars */
243 int last;
245 /* Position of the first non-zero coefficient in each ray. */
246 vector<int> neg;
248 /* Power without any "borrowing" */
249 vector<int> base_power;
250 /* Power after "borrowing" */
251 vector<int> power;
253 /* The non-zero coefficients in the rays of the vertex cone,
254 * except the first.
256 vector< vector<int> > mask;
257 /* For each ray, the (positive) power of each variable it contributes */
258 vector< vector<int> > selected;
259 /* The powers of each variable that still need to be selected */
260 vector<int> left;
262 HASH_MAP<vector<int>, const evalue *> cache;
264 reciprocal(vertex_cone &vc);
265 void start(vector<int> &power);
266 int next();
268 evalue *add(evalue *sum);
269 const evalue *get_coefficient();
270 ~reciprocal() {
271 HASH_MAP<vector<int>, const evalue *>::iterator j;
272 for (j = cache.begin(); j != cache.end(); ++j)
273 if ((*j).second)
274 evalue_free(const_cast<evalue *>((*j).second));
278 reciprocal::reciprocal(vertex_cone &vc) : vc(vc)
280 for (int i = 0; i < vc.dim; ++i) {
281 can_borrow_from.push_back(vector<int>(vc.dim));
282 has_borrowed_from.push_back(vector<int>(vc.dim));
283 mask.push_back(vector<int>(vc.dim));
284 selected.push_back(vector<int>(vc.dim));
286 can_borrow = vector<int>(vc.dim);
287 has_borrowed = vector<int>(vc.dim);
288 neg = vector<int>(vc.dim);
289 left = vector<int>(vc.dim);
291 for (int i = 0; i < vc.dim; ++i) {
292 int pos = First_Non_Zero(vc.coeff[i]->p, vc.coeff[i]->Size);
293 neg[i] = pos;
294 for (int j = pos+1; j < vc.dim; ++j)
295 if (value_notzero_p(vc.coeff[i]->p[j])) {
296 mask[i][j] = 1;
297 can_borrow_from[neg[i]][j] = 1;
302 /* Initialize power to the exponent of the todd product
303 * required to compute the coefficient in the full product
304 * of the term with exponent req_power, without any
305 * "borrowing".
307 void reciprocal::start(vector<int> &req_power)
309 power = req_power;
310 for (int j = 0; j < vc.dim; ++j)
311 power[neg[j]]++;
313 base_power = power;
315 last = -1;
316 for (int i = 0; i < vc.dim; ++i) {
317 can_borrow[i] = 0;
318 has_borrowed[i] = 0;
319 for (int j = i+1; j < vc.dim; ++j) {
320 has_borrowed_from[i][j] = 0;
321 if (can_borrow_from[i][j])
322 can_borrow[i] += power[j];
324 if (can_borrow[i])
325 last = i;
329 /* Set power to the next exponent of the todd product required
330 * and return 1 as long as there is any such exponent left.
332 int reciprocal::next()
334 int p = last;
335 while (p >= 0) {
336 if (has_borrowed[p] < can_borrow[p]) {
337 int j;
338 for (j = p+1; j < vc.dim; ++j)
339 if (can_borrow_from[p][j]) {
340 if (power[j] > 0)
341 break;
342 else if (has_borrowed_from[p][j]) {
343 power[j] += has_borrowed_from[p][j];
344 power[p] -= has_borrowed_from[p][j];
345 has_borrowed[p] -= has_borrowed_from[p][j];
346 has_borrowed_from[p][j] = 0;
349 if (j < vc.dim) {
350 has_borrowed_from[p][j]++;
351 has_borrowed[p]++;
352 power[p]++;
353 power[j]--;
354 return 1;
357 if (has_borrowed[p]) {
358 for (int j = p+1; j < vc.dim; ++j)
359 if (has_borrowed_from[p][j]) {
360 power[j] += has_borrowed_from[p][j];
361 has_borrowed_from[p][j] = 0;
363 power[p] -= has_borrowed[p];
364 has_borrowed[p] = 0;
366 --p;
368 return 0;
371 /* Add coefficient of selected powers of variables to sum
372 * and return the result.
373 * The contribution of each ray j of the vertex cone is
375 * ( K )
376 * ( k_{f+1}, \ldots, k_n ) (-1)^{K+1}
377 * c_f^{-K-1} c_{f+1}^{k_{f+1}} \cdots c_n^{k_n}
379 * K = \sum_{i=f+1}^n k_i
381 evalue *reciprocal::add(evalue *sum)
383 evalue *t = NULL;
384 for (int i = 0; i < vc.dim; ++i) {
385 evalue *c = ALLOC(evalue);
386 value_init(c->d);
387 value_init(c->x.n);
388 value_set_si(c->d, 1);
389 multinomial(selected[i], &c->x.n);
390 int s = 0;
391 for (int j = 0; j < vc.dim; ++j) {
392 if (selected[i][j] == 0)
393 continue;
394 value_multiply(c->x.n, c->x.n,
395 *(*vc.coeff_power[i][j])[selected[i][j]]);
396 s += selected[i][j];
398 value_multiply(c->d, c->d, *(*vc.coeff_power[i][neg[i]])[s+1]);
399 if (!(s % 2))
400 value_oppose(c->x.n, c->x.n);
401 if (value_neg_p(c->d)) {
402 value_oppose(c->d, c->d);
403 value_oppose(c->x.n, c->x.n);
405 if (!t)
406 t = c;
407 else {
408 emul(c, t);
409 evalue_free(c);
412 if (!sum)
413 sum = t;
414 else {
415 eadd(t, sum);
416 evalue_free(t);
418 return sum;
421 /* Return coefficient of the term with exponents powers by
422 * iterating over all combinations of exponents for each ray
423 * that sum up to the given exponents.
425 const evalue *reciprocal::get_coefficient()
427 evalue *c = NULL;
429 for (int j = 0; j < vc.dim; ++j)
430 left[j] = base_power[j] - power[j];
432 HASH_MAP<vector<int>, const evalue *>::iterator i;
433 i = cache.find(left);
434 if (i != cache.end())
435 return (*i).second;
437 int nz = first_non_zero(left);
438 if (nz == -1)
439 return cache[power] = add(c);
440 if (left[nz] > 0)
441 return NULL;
443 for (int i = 0; i < vc.dim; ++i)
444 for (int j = 0; j < vc.dim; ++j)
445 selected[i][j] = 0;
447 int level = vc.dim-1;
448 int p = vc.dim-1;
449 while (level >= 0) {
450 int nz = first_non_zero(left);
451 if (nz < neg[level] || (nz == neg[level] && left[nz] > 0)) {
452 assert(p == vc.dim-1);
453 --level;
454 continue;
456 if (nz == neg[level] && mask[level][p]) {
457 selected[level][p]++;
458 left[p]--;
459 left[neg[level]]++;
460 int nz = first_non_zero(left);
461 if (nz == -1)
462 c = add(c);
463 else if (left[nz] < 0) {
464 level = vc.dim-1;
465 p = vc.dim-1;
466 continue;
469 if (selected[level][p]) {
470 left[p] += selected[level][p];
471 left[neg[level]] -= selected[level][p];
472 selected[level][p] = 0;
474 if (--p < 0) {
475 --level;
476 p = vc.dim-1;
479 cache[left] = c;
481 return c;
484 struct laurent_summator_old : public signed_cone_consumer,
485 public vertex_decomposer {
486 const evalue *polynomial;
487 unsigned dim;
488 vertex_cone vc;
489 param_polynomial poly;
490 evalue *result;
491 unsigned max_power;
493 laurent_summator_old(const evalue *e, unsigned dim, Param_Polyhedron *PP) :
494 polynomial(e), dim(dim), vertex_decomposer(PP, *this),
495 vc(dim), poly(e, dim) {
496 max_power = dim + poly.degree();
497 result = NULL;
499 ~laurent_summator_old() {
500 if (result)
501 evalue_free(result);
503 virtual void handle(const signed_cone& sc, barvinok_options *options);
506 static int first_non_zero(const vec_ZZ& row)
508 for (int i = 0; i < row.length(); ++i)
509 if (row[i] != 0)
510 return i;
511 return -1;
514 void laurent_summator_old::handle(const signed_cone& sc, barvinok_options *options)
516 assert(sc.det == 1);
518 vc.init(sc.rays, V, max_power);
519 reciprocal recip(vc);
520 todd_product tp(vc);
521 for (int i = 0; i < poly.terms.size(); ++i) {
522 recip.start(poly.terms[i].powers);
523 do {
524 const evalue *c = recip.get_coefficient();
525 if (!c)
526 continue;
528 const evalue *t = tp.get_coefficient(recip.power);
530 evalue *f = evalue_dup(poly.terms[i].coeff);
531 if (sc.sign < 0)
532 evalue_negate(f);
533 for (int j = 0; j < dim; ++j)
534 evalue_mul(f, *factorial(poly.terms[i].powers[j]));
535 evalue_shift_variables(f, 0, -dim);
536 emul(c, f);
537 emul(t, f);
538 if (!result)
539 result = f;
540 else {
541 eadd(f, result);
542 evalue_free(f);
544 } while (recip.next());
546 vc.clear();
549 evalue *laurent_summate_old(Polyhedron *P, evalue *e, unsigned nvar,
550 struct barvinok_options *options)
552 Polyhedron *U, *TC;
553 Param_Polyhedron *PP;
554 struct evalue_section *s;
555 int nd = -1;
556 Param_Domain *PD;
557 evalue *sum;
559 U = Universe_Polyhedron(P->Dimension - nvar);
560 PP = Polyhedron2Param_Polyhedron(P, U, options);
562 for (nd = 0, PD = PP->D; PD; ++nd, PD = PD->next);
563 s = ALLOCN(struct evalue_section, nd);
565 TC = true_context(P, U, options->MaxRays);
566 FORALL_REDUCED_DOMAIN(PP, TC, nd, options, i, PD, rVD)
567 Param_Vertices *V;
568 laurent_summator_old ls(e, nvar, PP);
570 FORALL_PVertex_in_ParamPolyhedron(V, PD, PP) // _i internal counter
571 ls.decompose_at_vertex(V, _i, options);
572 END_FORALL_PVertex_in_ParamPolyhedron;
574 s[i].D = rVD;
575 s[i].E = ls.result;
576 ls.result = NULL;
577 END_FORALL_REDUCED_DOMAIN
578 Polyhedron_Free(TC);
579 Polyhedron_Free(U);
580 Param_Polyhedron_Free(PP);
582 sum = evalue_from_section_array(s, nd);
583 free(s);
585 return sum;