genfun.cc

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