genfun.cc

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