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