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