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