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