genfun.cc

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