genfun.cc

   1 #include <iostream>
   2 #include <vector>
   3 #include <assert.h>
   4 #include "config.h"
   5 #include <barvinok/genfun.h>
   6 #include <barvinok/barvinok.h>
   7 #include "conversion.h"
   8 #include "genfun_constructor.h"
   9 #include "mat_util.h"
  10
  11 using std::cout;
  12 using std::cerr;
  13 using std::endl;
  14 using std::pair;
  15 using std::vector;
  16
  17 bool short_rat_lex_smaller_denominator::operator()(const short_rat* r1,
  18                                                    const short_rat* r2) const
  19 {
  20     return lex_cmp(r1->d.power, r2->d.power) < 0;
  21 }
  22
  23 static void lex_order_terms(struct short_rat* rat)
  24 {
  25     for (int i = 0; i < rat->n.power.NumRows(); ++i) {
  26         int m = i;
  27         for (int j = i+1; j < rat->n.power.NumRows(); ++j)
  28             if (lex_cmp(rat->n.power[j], rat->n.power[m]) < 0)
  29                 m = j;
  30         if (m != i) {
  31             vec_ZZ tmp = rat->n.power[m];
  32             rat->n.power[m] = rat->n.power[i];
  33             rat->n.power[i] = tmp;
  34             QQ tmp_coeff = rat->n.coeff[m];
  35             rat->n.coeff[m] = rat->n.coeff[i];
  36             rat->n.coeff[i] = tmp_coeff;
  37         }
  38     }
  39 }
  40
  41 short_rat::short_rat(Value c)
  42 {
  43     n.coeff.SetLength(1);
  44     value2zz(c, n.coeff[0].n);
  45     n.coeff[0].d = 1;
  46     n.power.SetDims(1, 0);
  47     d.power.SetDims(0, 0);
  48 }
  49
  50 short_rat::short_rat(const QQ& c, const vec_ZZ& num, const mat_ZZ& den)
  51 {
  52     n.coeff.SetLength(1);
  53     ZZ g = GCD(c.n, c.d);
  54     n.coeff[0].n = c.n/g;
  55     n.coeff[0].d = c.d/g;
  56     n.power.SetDims(1, num.length());
  57     n.power[0] = num;
  58     d.power = den;
  59     normalize();
  60 }
  61
  62 void short_rat::normalize()
  63 {
  64     /* Make all powers in denominator reverse-lexico-positive */
  65     for (int i = 0; i < d.power.NumRows(); ++i) {
  66         int j;
  67         for (j = d.power.NumCols()-1; j >= 0; --j)
  68             if (d.power[i][j] != 0)
  69                 break;
  70         assert(j >= 0);
  71         if (d.power[i][j] < 0) {
  72             d.power[i] = -d.power[i];
  73             for (int k = 0; k < n.coeff.length(); ++k) {
  74                 n.coeff[k].n = -n.coeff[k].n;
  75                 n.power[k] += d.power[i];
  76             }
  77         }
  78     }
  79
  80     /* Order powers in denominator */
  81     lex_order_rows(d.power);
  82 }
  83
  84 void short_rat::add(short_rat *r)
  85 {
  86     for (int i = 0; i < r->n.power.NumRows(); ++i) {
  87         int len = n.coeff.length();
  88         int j;
  89         for (j = 0; j < len; ++j)
  90             if (r->n.power[i] == n.power[j])
  91                 break;
  92         if (j < len) {
  93             n.coeff[j] += r->n.coeff[i];
  94             if (n.coeff[j].n == 0) {
  95                 if (j < len-1) {
  96                     n.power[j] = n.power[len-1];
  97                     n.coeff[j] = n.coeff[len-1];
  98                 }
  99                 int dim = n.power.NumCols();
 100                 n.coeff.SetLength(len-1);
 101                 n.power.SetDims(len-1, dim);
 102             }
 103         } else {
 104             int dim = n.power.NumCols();
 105             n.coeff.SetLength(len+1);
 106             n.power.SetDims(len+1, dim);
 107             n.coeff[len] = r->n.coeff[i];
 108             n.power[len] = r->n.power[i];
 109         }
 110     }
 111 }
 112
 113 bool short_rat::reduced()
 114 {
 115     int dim = n.power.NumCols();
 116     lex_order_terms(this);
 117     if (n.power.NumRows() % 2 == 0) {
 118         if (n.coeff[0].n == -n.coeff[1].n &&
 119             n.coeff[0].d == n.coeff[1].d) {
 120             vec_ZZ step = n.power[1] - n.power[0];
 121             int k;
 122             for (k = 1; k < n.power.NumRows()/2; ++k) {
 123                 if (n.coeff[2*k].n != -n.coeff[2*k+1].n ||
 124                     n.coeff[2*k].d != n.coeff[2*k+1].d)
 125                     break;
 126                 if (step != n.power[2*k+1] - n.power[2*k])
 127                     break;
 128             }
 129             if (k == n.power.NumRows()/2) {
 130                 for (k = 0; k < d.power.NumRows(); ++k)
 131                     if (d.power[k] == step)
 132                         break;
 133                 if (k < d.power.NumRows()) {
 134                     for (++k; k < d.power.NumRows(); ++k)
 135                         d.power[k-1] = d.power[k];
 136                     d.power.SetDims(k-1, dim);
 137                     for (k = 1; k < n.power.NumRows()/2; ++k) {
 138                         n.coeff[k] = n.coeff[2*k];
 139                         n.power[k] = n.power[2*k];
 140                     }
 141                     n.coeff.SetLength(k);
 142                     n.power.SetDims(k, dim);
 143                     return true;
 144                 }
 145             }
 146         }
 147     }
 148     return false;
 149 }
 150
 151 gen_fun::gen_fun(Value c)
 152 {
 153     short_rat *r = new short_rat(c);
 154     context = Universe_Polyhedron(0);
 155     term.insert(r);
 156 }
 157
 158 void gen_fun::add(const QQ& c, const vec_ZZ& num, const mat_ZZ& den)
 159 {
 160     if (c.n == 0)
 161         return;
 162
 163     short_rat * r = new short_rat(c, num, den);
 164
 165     short_rat_list::iterator i = term.find(r);
 166     while (i != term.end()) {
 167         (*i)->add(r);
 168         if ((*i)->n.coeff.length() == 0) {
 169             delete *i;
 170             term.erase(i);
 171         } else if ((*i)->reduced()) {
 172             delete r;
 173             /* we've modified term[i], so remove it
 174              * and add it back again
 175              */
 176             r = *i;
 177             term.erase(i);
 178             i = term.find(r);
 179             continue;
 180         }
 181         delete r;
 182         return;
 183     }
 184
 185     term.insert(r);
 186 }
 187
 188 void gen_fun::add(const QQ& c, const gen_fun *gf)
 189 {
 190     QQ p;
 191     for (short_rat_list::iterator i = gf->term.begin(); i != gf->term.end(); ++i) {
 192         for (int j = 0; j < (*i)->n.power.NumRows(); ++j) {
 193             p = c;
 194             p *= (*i)->n.coeff[j];
 195             add(p, (*i)->n.power[j], (*i)->d.power);
 196         }
 197     }
 198 }
 199
 200 static void split_param_compression(Matrix *CP, mat_ZZ& map, vec_ZZ& offset)
 201 {
 202     Matrix *T = Transpose(CP);
 203     matrix2zz(T, map, T->NbRows-1, T->NbColumns-1);
 204     values2zz(T->p[T->NbRows-1], offset, T->NbColumns-1);
 205     Matrix_Free(T);
 206 }
 207
 208 /*
 209  * Perform the substitution specified by CP
 210  *
 211  * CP is a homogeneous matrix that maps a set of "compressed parameters"
 212  * to the original set of parameters.
 213  *
 214  * This function is applied to a gen_fun computed with the compressed parameters
 215  * and adapts it to refer to the original parameters.
 216  *
 217  * That is, if y are the compressed parameters and x = A y + b are the original
 218  * parameters, then we want the coefficient of the monomial t^y in the original
 219  * generating function to be the coefficient of the monomial u^x in the resulting
 220  * generating function.
 221  * The original generating function has the form
 222  *
 223  *        a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
 224  *
 225  * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
 226  *
 227  *         a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
 228  *
 229  *         = a u^{A m + b}/(1-u^{A n})
 230  *
 231  * Therefore, we multiply the powers m and n in both numerator and denominator by A
 232  * and add b to the power in the numerator.
 233  * Since the above powers are stored as row vectors m^T and n^T,
 234  * we compute, say, m'^T = m^T A^T to obtain m' = A m.
 235  *
 236  * The pair (map, offset) contains the same information as CP.
 237  * map is the transpose of the linear part of CP, while offset is the constant part.
 238  */
 239 void gen_fun::substitute(Matrix *CP)
 240 {
 241     mat_ZZ map;
 242     vec_ZZ offset;
 243     split_param_compression(CP, map, offset);
 244     Polyhedron *C = Polyhedron_Image(context, CP, 0);
 245     Polyhedron_Free(context);
 246     context = C;
 247
 248     short_rat_list new_term;
 249     for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
 250         short_rat *r = (*i);
 251         r->d.power *= map;
 252         r->n.power *= map;
 253         for (int j = 0; j < r->n.power.NumRows(); ++j)
 254             r->n.power[j] += offset;
 255         r->normalize();
 256         new_term.insert(r);
 257     }
 258     term.swap(new_term);
 259 }
 260
 261 struct cone {
 262     int     *pos;
 263     vector<pair<Vector *, QQ> >     vertices;
 264     cone(int *pos) : pos(pos) {}
 265 };
 266
 267 #ifndef HAVE_COMPRESS_PARMS
 268 static Matrix *compress_parms(Matrix *M, unsigned nparam)
 269 {
 270     assert(0);
 271 }
 272 #endif
 273
 274 struct parallel_polytopes {
 275     gf_base *red;
 276     Polyhedron *context;
 277     Matrix *Constraints;
 278     Matrix *CP, *T;
 279     int dim;
 280     int nparam;
 281     vector<cone>    cones;
 282     barvinok_options    *options;
 283
 284     parallel_polytopes(int n, Polyhedron *context, int nparam,
 285                        barvinok_options *options) :
 286                         context(context), dim(-1), nparam(nparam),
 287                         options(options) {
 288         red = NULL;
 289         Constraints = NULL;
 290         CP = NULL;
 291         T = NULL;
 292     }
 293     bool add(const QQ& c, Polyhedron *P) {
 294         int i;
 295
 296         for (i = 0; i < P->NbEq; ++i)
 297             if (First_Non_Zero(P->Constraint[i]+1,
 298                                P->Dimension-nparam) == -1)
 299                 break;
 300         if (i < P->NbEq)
 301             return false;
 302
 303         Polyhedron *Q = remove_equalities_p(Polyhedron_Copy(P), P->Dimension-nparam,
 304                                             NULL);
 305         POL_ENSURE_VERTICES(Q);
 306         if (emptyQ(Q)) {
 307             Polyhedron_Free(Q);
 308             return true;
 309         }
 310
 311         if (Q->NbEq != 0) {
 312             Polyhedron *R;
 313             if (!CP) {
 314                 Matrix *M;
 315                 M = Matrix_Alloc(Q->NbEq, Q->Dimension+2);
 316                 Vector_Copy(Q->Constraint[0], M->p[0], Q->NbEq * (Q->Dimension+2));
 317                 CP = compress_parms(M, nparam);
 318                 T = align_matrix(CP, Q->Dimension+1);
 319                 Matrix_Free(M);
 320             }
 321             R = Polyhedron_Preimage(Q, T, options->MaxRays);
 322             Polyhedron_Free(Q);
 323             Q = remove_equalities_p(R, R->Dimension-nparam, NULL);
 324         }
 325         assert(Q->NbEq == 0);
 326
 327         if (First_Non_Zero(Q->Constraint[Q->NbConstraints-1]+1, Q->Dimension) == -1)
 328             Q->NbConstraints--;
 329
 330         if (!Constraints) {
 331             dim = Q->Dimension;
 332             red = gf_base::create(Polyhedron_Copy(context), dim, nparam, options);
 333             red->base->init(Q);
 334             Constraints = Matrix_Alloc(Q->NbConstraints, Q->Dimension);
 335             for (int i = 0; i < Q->NbConstraints; ++i) {
 336                 Vector_Copy(Q->Constraint[i]+1, Constraints->p[i], Q->Dimension);
 337             }
 338         } else {
 339             assert(Q->Dimension == dim);
 340             for (int i = 0; i < Q->NbConstraints; ++i) {
 341                 int j;
 342                 for (j = 0; j < Constraints->NbRows; ++j)
 343                     if (Vector_Equal(Q->Constraint[i]+1, Constraints->p[j],
 344                                         Q->Dimension))
 345                         break;
 346                 assert(j < Constraints->NbRows);
 347             }
 348         }
 349
 350         for (int i = 0; i < Q->NbRays; ++i) {
 351             if (!value_pos_p(Q->Ray[i][dim+1]))
 352                 continue;
 353
 354             Polyhedron *C = supporting_cone(Q, i);
 355
 356             if (First_Non_Zero(C->Constraint[C->NbConstraints-1]+1,
 357                                C->Dimension) == -1)
 358                 C->NbConstraints--;
 359
 360             int *pos = new int[1+C->NbConstraints];
 361             pos[0] = C->NbConstraints;
 362             int l = 0;
 363             for (int k = 0; k < Constraints->NbRows; ++k) {
 364                 for (int j = 0; j < C->NbConstraints; ++j) {
 365                     if (Vector_Equal(C->Constraint[j]+1, Constraints->p[k],
 366                                      C->Dimension)) {
 367                         pos[1+l++] = k;
 368                         break;
 369                     }
 370                 }
 371             }
 372             assert(l == C->NbConstraints);
 373
 374             int j;
 375             for (j = 0; j < cones.size(); ++j)
 376                 if (!memcmp(pos, cones[j].pos, (1+C->NbConstraints)*sizeof(int)))
 377                     break;
 378             if (j == cones.size())
 379                 cones.push_back(cone(pos));
 380             else
 381                 delete [] pos;
 382
 383             Polyhedron_Free(C);
 384
 385             int k;
 386             for (k = 0; k < cones[j].vertices.size(); ++k)
 387                 if (Vector_Equal(Q->Ray[i]+1, cones[j].vertices[k].first->p,
 388                                  Q->Dimension+1))
 389                     break;
 390
 391             if (k == cones[j].vertices.size()) {
 392                 Vector *vertex = Vector_Alloc(Q->Dimension+1);
 393                 Vector_Copy(Q->Ray[i]+1, vertex->p, Q->Dimension+1);
 394                 cones[j].vertices.push_back(pair<Vector*,QQ>(vertex, c));
 395             } else {
 396                 cones[j].vertices[k].second += c;
 397                 if (cones[j].vertices[k].second.n == 0) {
 398                     int size = cones[j].vertices.size();
 399                     Vector_Free(cones[j].vertices[k].first);
 400                     if (k < size-1)
 401                         cones[j].vertices[k] = cones[j].vertices[size-1];
 402                     cones[j].vertices.pop_back();
 403                 }
 404             }
 405         }
 406
 407         Polyhedron_Free(Q);
 408         return true;
 409     }
 410     gen_fun *compute() {
 411         if (!red)
 412             return NULL;
 413         for (int i = 0; i < cones.size(); ++i) {
 414             Matrix *M = Matrix_Alloc(cones[i].pos[0], 1+Constraints->NbColumns+1);
 415             Polyhedron *Cone;
 416             for (int j = 0; j <cones[i].pos[0]; ++j) {
 417                 value_set_si(M->p[j][0], 1);
 418                 Vector_Copy(Constraints->p[cones[i].pos[1+j]], M->p[j]+1,
 419                             Constraints->NbColumns);
 420             }
 421             Cone = Constraints2Polyhedron(M, options->MaxRays);
 422             Matrix_Free(M);
 423             for (int j = 0; j < cones[i].vertices.size(); ++j) {
 424                 red->base->do_vertex_cone(cones[i].vertices[j].second,
 425                                           Polyhedron_Copy(Cone),
 426                                           cones[i].vertices[j].first->p, options);
 427             }
 428             Polyhedron_Free(Cone);
 429         }
 430         if (CP)
 431             red->gf->substitute(CP);
 432         return red->gf;
 433     }
 434     void print(std::ostream& os) const {
 435         for (int i = 0; i < cones.size(); ++i) {
 436             os << "[";
 437             for (int j = 0; j < cones[i].pos[0]; ++j) {
 438                 if (j)
 439                     os << ", ";
 440                 os << cones[i].pos[1+j];
 441             }
 442             os << "]" << endl;
 443             for (int j = 0; j < cones[i].vertices.size(); ++j) {
 444                 Vector_Print(stderr, P_VALUE_FMT, cones[i].vertices[j].first);
 445                 os << cones[i].vertices[j].second << endl;
 446             }
 447         }
 448     }
 449     ~parallel_polytopes() {
 450         for (int i = 0; i < cones.size(); ++i) {
 451             delete [] cones[i].pos;
 452             for (int j = 0; j < cones[i].vertices.size(); ++j)
 453                 Vector_Free(cones[i].vertices[j].first);
 454         }
 455         if (Constraints)
 456             Matrix_Free(Constraints);
 457         if (CP)
 458             Matrix_Free(CP);
 459         if (T)
 460             Matrix_Free(T);
 461         delete red;
 462     }
 463 };
 464
 465 gen_fun *gen_fun::Hadamard_product(const gen_fun *gf, barvinok_options *options)
 466 {
 467     QQ one(1, 1);
 468     Polyhedron *C = DomainIntersection(context, gf->context, options->MaxRays);
 469     Polyhedron *U = Universe_Polyhedron(C->Dimension);
 470     gen_fun *sum = new gen_fun(C);
 471     for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
 472         for (short_rat_list::iterator i2 = gf->term.begin(); i2 != gf->term.end();
 473                                                              ++i2) {
 474             int d = (*i)->d.power.NumCols();
 475             int k1 = (*i)->d.power.NumRows();
 476             int k2 = (*i2)->d.power.NumRows();
 477             assert((*i)->d.power.NumCols() == (*i2)->d.power.NumCols());
 478
 479             parallel_polytopes pp((*i)->n.power.NumRows() *
 480                                   (*i2)->n.power.NumRows(),
 481                                   sum->context, d, options);
 482
 483             for (int j = 0; j < (*i)->n.power.NumRows(); ++j) {
 484                 for (int j2 = 0; j2 < (*i2)->n.power.NumRows(); ++j2) {
 485                     Matrix *M = Matrix_Alloc(k1+k2+d+d, 1+k1+k2+d+1);
 486                     for (int k = 0; k < k1+k2; ++k) {
 487                         value_set_si(M->p[k][0], 1);
 488                         value_set_si(M->p[k][1+k], 1);
 489                     }
 490                     for (int k = 0; k < d; ++k) {
 491                         value_set_si(M->p[k1+k2+k][1+k1+k2+k], -1);
 492                         zz2value((*i)->n.power[j][k], M->p[k1+k2+k][1+k1+k2+d]);
 493                         for (int l = 0; l < k1; ++l)
 494                             zz2value((*i)->d.power[l][k], M->p[k1+k2+k][1+l]);
 495                     }
 496                     for (int k = 0; k < d; ++k) {
 497                         value_set_si(M->p[k1+k2+d+k][1+k1+k2+k], -1);
 498                         zz2value((*i2)->n.power[j2][k],
 499                                  M->p[k1+k2+d+k][1+k1+k2+d]);
 500                         for (int l = 0; l < k2; ++l)
 501                             zz2value((*i2)->d.power[l][k],
 502                                      M->p[k1+k2+d+k][1+k1+l]);
 503                     }
 504                     Polyhedron *P = Constraints2Polyhedron(M, options->MaxRays);
 505                     Matrix_Free(M);
 506
 507                     QQ c = (*i)->n.coeff[j];
 508                     c *= (*i2)->n.coeff[j2];
 509                     if (!pp.add(c, P)) {
 510                         gen_fun *t = barvinok_series(P, U, options->MaxRays);
 511                         sum->add(c, t);
 512                         delete t;
 513                     }
 514
 515                     Polyhedron_Free(P);
 516                 }
 517             }
 518
 519             gen_fun *t = pp.compute();
 520             if (t) {
 521                 sum->add(one, t);
 522                 delete t;
 523             }
 524         }
 525     }
 526     Polyhedron_Free(U);
 527     return sum;
 528 }
 529
 530 void gen_fun::add_union(gen_fun *gf, barvinok_options *options)
 531 {
 532     QQ one(1, 1), mone(-1, 1);
 533
 534     gen_fun *hp = Hadamard_product(gf, options);
 535     add(one, gf);
 536     add(mone, hp);
 537     delete hp;
 538 }
 539
 540 static void Polyhedron_Shift(Polyhedron *P, Vector *offset)
 541 {
 542     Value tmp;
 543     value_init(tmp);
 544     for (int i = 0; i < P->NbConstraints; ++i) {
 545         Inner_Product(P->Constraint[i]+1, offset->p, P->Dimension, &tmp);
 546         value_subtract(P->Constraint[i][1+P->Dimension],
 547                        P->Constraint[i][1+P->Dimension], tmp);
 548     }
 549     for (int i = 0; i < P->NbRays; ++i) {
 550         if (value_notone_p(P->Ray[i][0]))
 551             continue;
 552         if (value_zero_p(P->Ray[i][1+P->Dimension]))
 553             continue;
 554         Vector_Combine(P->Ray[i]+1, offset->p, P->Ray[i]+1,
 555                        P->Ray[i][0], P->Ray[i][1+P->Dimension], P->Dimension);
 556     }
 557     value_clear(tmp);
 558 }
 559
 560 void gen_fun::shift(const vec_ZZ& offset)
 561 {
 562     for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
 563         for (int j = 0; j < (*i)->n.power.NumRows(); ++j)
 564             (*i)->n.power[j] += offset;
 565
 566     Vector *v = Vector_Alloc(offset.length());
 567     zz2values(offset, v->p);
 568     Polyhedron_Shift(context, v);
 569     Vector_Free(v);
 570 }
 571
 572 /* Divide the generating functin by 1/(1-z^power).
 573  * The effect on the corresponding explicit function f(x) is
 574  * f'(x) = \sum_{i=0}^\infty f(x - i * power)
 575  */
 576 void gen_fun::divide(const vec_ZZ& power)
 577 {
 578     for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
 579         int r = (*i)->d.power.NumRows();
 580         int c = (*i)->d.power.NumCols();
 581         (*i)->d.power.SetDims(r+1, c);
 582         (*i)->d.power[r] = power;
 583     }
 584
 585     Vector *v = Vector_Alloc(1+power.length()+1);
 586     value_set_si(v->p[0], 1);
 587     zz2values(power, v->p+1);
 588     Polyhedron *C = AddRays(v->p, 1, context, context->NbConstraints+1);
 589     Vector_Free(v);
 590     Polyhedron_Free(context);
 591     context = C;
 592 }
 593
 594 static void print_power(std::ostream& os, QQ& c, vec_ZZ& p,
 595                         unsigned int nparam, char **param_name)
 596 {
 597     bool first = true;
 598
 599     for (int i = 0; i < p.length(); ++i) {
 600         if (p[i] == 0)
 601             continue;
 602         if (first) {
 603             if (c.n == -1 && c.d == 1)
 604                 os << "-";
 605             else if (c.n != 1 || c.d != 1) {
 606                 os << c.n;
 607                 if (c.d != 1)
 608                     os << " / " << c.d;
 609                 os << "*";
 610             }
 611             first = false;
 612         } else
 613             os << "*";
 614         if (i < nparam)
 615             os << param_name[i];
 616         else
 617             os << "x" << i;
 618         if (p[i] == 1)
 619             continue;
 620         if (p[i] < 0)
 621             os << "^(" << p[i] << ")";
 622         else
 623             os << "^" << p[i];
 624     }
 625     if (first) {
 626         os << c.n;
 627         if (c.d != 1)
 628             os << " / " << c.d;
 629     }
 630 }
 631
 632 void gen_fun::print(std::ostream& os, unsigned int nparam, char **param_name) const
 633 {
 634     QQ mone(-1, 1);
 635     for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
 636         if (i != term.begin())
 637             os << " + ";
 638         os << "(";
 639         for (int j = 0; j < (*i)->n.coeff.length(); ++j) {
 640             if (j != 0 && (*i)->n.coeff[j].n > 0)
 641                 os << "+";
 642             print_power(os, (*i)->n.coeff[j], (*i)->n.power[j],
 643                         nparam, param_name);
 644         }
 645         os << ")/(";
 646         for (int j = 0; j < (*i)->d.power.NumRows(); ++j) {
 647             if (j != 0)
 648                 os << " * ";
 649             os << "(1";
 650             print_power(os, mone, (*i)->d.power[j], nparam, param_name);
 651             os << ")";
 652         }
 653         os << ")";
 654     }
 655 }
 656
 657 gen_fun::operator evalue *() const
 658 {
 659     evalue *EP = NULL;
 660     evalue factor;
 661     value_init(factor.d);
 662     value_init(factor.x.n);
 663     for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
 664         unsigned nvar = (*i)->d.power.NumRows();
 665         unsigned nparam = (*i)->d.power.NumCols();
 666         Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + nparam + 1);
 667         mat_ZZ& d = (*i)->d.power;
 668         Polyhedron *U = context ? context : Universe_Polyhedron(nparam);
 669
 670         for (int j = 0; j < (*i)->n.coeff.length(); ++j) {
 671             for (int r = 0; r < nparam; ++r) {
 672                 value_set_si(C->p[r][0], 0);
 673                 for (int c = 0; c < nvar; ++c) {
 674                     zz2value(d[c][r], C->p[r][1+c]);
 675                 }
 676                 Vector_Set(&C->p[r][1+nvar], 0, nparam);
 677                 value_set_si(C->p[r][1+nvar+r], -1);
 678                 zz2value((*i)->n.power[j][r], C->p[r][1+nvar+nparam]);
 679             }
 680             for (int r = 0; r < nvar; ++r) {
 681                 value_set_si(C->p[nparam+r][0], 1);
 682                 Vector_Set(&C->p[nparam+r][1], 0, nvar + nparam + 1);
 683                 value_set_si(C->p[nparam+r][1+r], 1);
 684             }
 685             Polyhedron *P = Constraints2Polyhedron(C, 0);
 686             evalue *E = barvinok_enumerate_ev(P, U, 0);
 687             Polyhedron_Free(P);
 688             if (EVALUE_IS_ZERO(*E)) {
 689                 free_evalue_refs(E);
 690                 free(E);
 691                 continue;
 692             }
 693             zz2value((*i)->n.coeff[j].n, factor.x.n);
 694             zz2value((*i)->n.coeff[j].d, factor.d);
 695             emul(&factor, E);
 696             /*
 697             Matrix_Print(stdout, P_VALUE_FMT, C);
 698             char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
 699             print_evalue(stdout, E, test);
 700             */
 701             if (!EP)
 702                 EP = E;
 703             else {
 704                 eadd(E, EP);
 705                 free_evalue_refs(E);
 706                 free(E);
 707             }
 708         }
 709         Matrix_Free(C);
 710         if (!context)
 711             Polyhedron_Free(U);
 712     }
 713     value_clear(factor.d);
 714     value_clear(factor.x.n);
 715     return EP;
 716 }
 717
 718 void gen_fun::coefficient(Value* params, Value* c) const
 719 {
 720     if (context && !in_domain(context, params)) {
 721         value_set_si(*c, 0);
 722         return;
 723     }
 724
 725     evalue part;
 726     value_init(part.d);
 727     value_init(part.x.n);
 728     evalue sum;
 729     value_init(sum.d);
 730     evalue_set_si(&sum, 0, 1);
 731     Value tmp;
 732     value_init(tmp);
 733
 734     for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i) {
 735         unsigned nvar = (*i)->d.power.NumRows();
 736         unsigned nparam = (*i)->d.power.NumCols();
 737         Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + 1);
 738         mat_ZZ& d = (*i)->d.power;
 739
 740         for (int j = 0; j < (*i)->n.coeff.length(); ++j) {
 741             C->NbRows = nparam+nvar;
 742             for (int r = 0; r < nparam; ++r) {
 743                 value_set_si(C->p[r][0], 0);
 744                 for (int c = 0; c < nvar; ++c) {
 745                     zz2value(d[c][r], C->p[r][1+c]);
 746                 }
 747                 zz2value((*i)->n.power[j][r], C->p[r][1+nvar]);
 748                 value_subtract(C->p[r][1+nvar], C->p[r][1+nvar], params[r]);
 749             }
 750             for (int r = 0; r < nvar; ++r) {
 751                 value_set_si(C->p[nparam+r][0], 1);
 752                 Vector_Set(&C->p[nparam+r][1], 0, nvar + 1);
 753                 value_set_si(C->p[nparam+r][1+r], 1);
 754             }
 755             Polyhedron *P = Constraints2Polyhedron(C, 0);
 756             if (emptyQ(P)) {
 757                 Polyhedron_Free(P);
 758                 continue;
 759             }
 760             barvinok_count(P, &tmp, 0);
 761             Polyhedron_Free(P);
 762             if (value_zero_p(tmp))
 763                 continue;
 764             zz2value((*i)->n.coeff[j].n, part.x.n);
 765             zz2value((*i)->n.coeff[j].d, part.d);
 766             value_multiply(part.x.n, part.x.n, tmp);
 767             eadd(&part, &sum);
 768         }
 769         Matrix_Free(C);
 770     }
 771
 772     assert(value_one_p(sum.d));
 773     value_assign(*c, sum.x.n);
 774
 775     value_clear(tmp);
 776     value_clear(part.d);
 777     value_clear(part.x.n);
 778     value_clear(sum.d);
 779     value_clear(sum.x.n);
 780 }
 781
 782 gen_fun *gen_fun::summate(int nvar, barvinok_options *options) const
 783 {
 784     int dim = context->Dimension;
 785     int nparam = dim - nvar;
 786     reducer *red;
 787     gen_fun *gf;
 788
 789     if (options->incremental_specialization == 1) {
 790         red = new partial_ireducer(Polyhedron_Project(context, nparam), dim, nparam);
 791     } else
 792         red = new partial_reducer(Polyhedron_Project(context, nparam), dim, nparam);
 793     red->init(context);
 794     for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
 795         for (int j = 0; j < (*i)->n.power.NumRows(); ++j)
 796             red->reduce((*i)->n.coeff[j], (*i)->n.power[j], (*i)->d.power);
 797     gf = red->get_gf();
 798     delete red;
 799     return gf;
 800 }
 801
 802 /* returns true if the set was finite and false otherwise */
 803 bool gen_fun::summate(Value *sum) const
 804 {
 805     if (term.size() == 0) {
 806         value_set_si(*sum, 0);
 807         return true;
 808     }
 809
 810     int maxlen = 0;
 811     for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
 812         if ((*i)->d.power.NumRows() > maxlen)
 813             maxlen = (*i)->d.power.NumRows();
 814
 815     infinite_icounter cnt((*term.begin())->d.power.NumCols(), maxlen);
 816     for (short_rat_list::iterator i = term.begin(); i != term.end(); ++i)
 817         for (int j = 0; j < (*i)->n.power.NumRows(); ++j)
 818             cnt.reduce((*i)->n.coeff[j], (*i)->n.power[j], (*i)->d.power);
 819
 820     for (int i = 1; i <= maxlen; ++i)
 821         if (value_notzero_p(mpq_numref(cnt.count[i]))) {
 822             value_set_si(*sum, -1);
 823             return false;
 824         }
 825
 826     assert(value_one_p(mpq_denref(cnt.count[0])));
 827     value_assign(*sum, mpq_numref(cnt.count[0]));
 828     return true;
 829 }