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