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