genfun.cc

   1 #include <iostream>
   2 #include <assert.h>
   3 #include "config.h"
   4 #include <barvinok/genfun.h>
   5 #include <barvinok/barvinok.h>
   6 #include "conversion.h"
   7 #include "mat_util.h"
   8
   9 using std::cout;
  10
  11 static int lex_cmp(mat_ZZ& a, mat_ZZ& b)
  12 {
  13     assert(a.NumCols() == b.NumCols());
  14     int alen = a.NumRows();
  15     int blen = b.NumRows();
  16     int len = alen < blen ? alen : blen;
  17
  18     for (int i = 0; i < len; ++i) {
  19         int s = lex_cmp(a[i], b[i]);
  20         if (s)
  21             return s;
  22     }
  23     return alen-blen;
  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             tmp = rat->n.coeff[m];
  38             rat->n.coeff[m] = rat->n.coeff[i];
  39             rat->n.coeff[i] = tmp;
  40         }
  41     }
  42 }
  43
  44 void short_rat::add(short_rat *r)
  45 {
  46     for (int i = 0; i < r->n.power.NumRows(); ++i) {
  47         int len = n.coeff.NumRows();
  48         int j;
  49         for (j = 0; j < len; ++j)
  50             if (r->n.power[i] == n.power[j])
  51                 break;
  52         if (j < len) {
  53             ZZ g = GCD(r->n.coeff[i][1], n.coeff[j][1]);
  54             ZZ num = n.coeff[j][0] * (r->n.coeff[i][1] / g) +
  55                        (n.coeff[j][1] / g) * r->n.coeff[i][0];
  56             ZZ d = n.coeff[j][1] / g * r->n.coeff[i][1];
  57             if (num != 0) {
  58                 g = GCD(num,d);
  59                 n.coeff[j][0] = num/g;
  60                 n.coeff[j][1] = d/g;
  61             } else {
  62                 if (j < len-1) {
  63                     n.power[j] = n.power[len-1];
  64                     n.coeff[j] = n.coeff[len-1];
  65                 }
  66                 int dim = n.power.NumCols();
  67                 n.coeff.SetDims(len-1, 2);
  68                 n.power.SetDims(len-1, dim);
  69             }
  70         } else {
  71             int dim = n.power.NumCols();
  72             n.coeff.SetDims(len+1, 2);
  73             n.power.SetDims(len+1, dim);
  74             n.coeff[len] = r->n.coeff[i];
  75             n.power[len] = r->n.power[i];
  76         }
  77     }
  78 }
  79
  80 bool short_rat::reduced()
  81 {
  82     int dim = n.power.NumCols();
  83     lex_order_terms(this);
  84     if (n.power.NumRows() % 2 == 0) {
  85         if (n.coeff[0][0] == -n.coeff[1][0] &&
  86             n.coeff[0][1] == n.coeff[1][1]) {
  87             vec_ZZ step = n.power[1] - n.power[0];
  88             int k;
  89             for (k = 1; k < n.power.NumRows()/2; ++k) {
  90                 if (n.coeff[2*k][0] != -n.coeff[2*k+1][0] ||
  91                     n.coeff[2*k][1] != n.coeff[2*k+1][1])
  92                     break;
  93                 if (step != n.power[2*k+1] - n.power[2*k])
  94                     break;
  95             }
  96             if (k == n.power.NumRows()/2) {
  97                 for (k = 0; k < d.power.NumRows(); ++k)
  98                     if (d.power[k] == step)
  99                         break;
 100                 if (k < d.power.NumRows()) {
 101                     for (++k; k < d.power.NumRows(); ++k)
 102                         d.power[k-1] = d.power[k];
 103                     d.power.SetDims(k-1, dim);
 104                     for (k = 1; k < n.power.NumRows()/2; ++k) {
 105                         n.coeff[k] = n.coeff[2*k];
 106                         n.power[k] = n.power[2*k];
 107                     }
 108                     n.coeff.SetDims(k, 2);
 109                     n.power.SetDims(k, dim);
 110                     return true;
 111                 }
 112             }
 113         }
 114     }
 115     return false;
 116 }
 117
 118 void gen_fun::add(const ZZ& cn, const ZZ& cd, const vec_ZZ& num,
 119                   const mat_ZZ& den)
 120 {
 121     if (cn == 0)
 122         return;
 123
 124     short_rat * r = new short_rat;
 125     r->n.coeff.SetDims(1, 2);
 126     ZZ g = GCD(cn, cd);
 127     r->n.coeff[0][0] = cn/g;
 128     r->n.coeff[0][1] = cd/g;
 129     r->n.power.SetDims(1, num.length());
 130     r->n.power[0] = num;
 131     r->d.power = den;
 132
 133     /* Make all powers in denominator lexico-positive */
 134     for (int i = 0; i < r->d.power.NumRows(); ++i) {
 135         int j;
 136         for (j = 0; j < r->d.power.NumCols(); ++j)
 137             if (r->d.power[i][j] != 0)
 138                 break;
 139         if (r->d.power[i][j] < 0) {
 140             r->d.power[i] = -r->d.power[i];
 141             r->n.coeff[0][0] = -r->n.coeff[0][0];
 142             r->n.power[0] += r->d.power[i];
 143         }
 144     }
 145
 146     /* Order powers in denominator */
 147     lex_order_rows(r->d.power);
 148
 149     for (int i = 0; i < term.size(); ++i)
 150         if (lex_cmp(term[i]->d.power, r->d.power) == 0) {
 151             term[i]->add(r);
 152             if (term[i]->n.coeff.NumRows() == 0) {
 153                 delete term[i];
 154                 if (i != term.size()-1)
 155                     term[i] = term[term.size()-1];
 156                 term.pop_back();
 157             } else if (term[i]->reduced()) {
 158                 delete r;
 159                 /* we've modified term[i], so removed it
 160                  * and add it back again
 161                  */
 162                 r = term[i];
 163                 if (i != term.size()-1)
 164                     term[i] = term[term.size()-1];
 165                 term.pop_back();
 166                 i = -1;
 167                 continue;
 168             }
 169             delete r;
 170             return;
 171         }
 172
 173     term.push_back(r);
 174 }
 175
 176 void gen_fun::add(const ZZ& cn, const ZZ& cd, const gen_fun *gf)
 177 {
 178     ZZ n, d;
 179     for (int i = 0; i < gf->term.size(); ++i) {
 180         for (int j = 0; j < gf->term[i]->n.power.NumRows(); ++j) {
 181             n = cn * gf->term[i]->n.coeff[j][0];
 182             d = cd * gf->term[i]->n.coeff[j][1];
 183             add(n, d, gf->term[i]->n.power[j], gf->term[i]->d.power);
 184         }
 185     }
 186 }
 187
 188 /*
 189  * Perform the substitution specified by CP and (map, offset)
 190  *
 191  * CP is a homogeneous matrix that maps a set of "compressed parameters"
 192  * to the original set of parameters.
 193  *
 194  * This function is applied to a gen_fun computed with the compressed parameters
 195  * and adapts it to refer to the original parameters.
 196  *
 197  * That is, if y are the compressed parameters and x = A y + b are the original
 198  * parameters, then we want the coefficient of the monomial t^y in the original
 199  * generating function to be the coefficient of the monomial u^x in the resulting
 200  * generating function.
 201  * The original generating function has the form
 202  *
 203  *        a t^m/(1-t^n) = a t^m + a t^{m+n} + a t^{m+2n} + ...
 204  *
 205  * Since each term t^y should correspond to a term u^x, with x = A y + b, we want
 206  *
 207  *         a u^{A m + b} + a u^{A (m+n) + b} + a u^{A (m+2n) +b} + ... =
 208  *
 209  *         = a u^{A m + b}/(1-u^{A n})
 210  *
 211  * Therefore, we multiply the powers m and n in both numerator and denominator by A
 212  * and add b to the power in the numerator.
 213  * Since the above powers are stored as row vectors m^T and n^T,
 214  * we compute, say, m'^T = m^T A^T to obtain m' = A m.
 215  *
 216  * The pair (map, offset) contains the same information as CP.
 217  * map is the transpose of the linear part of CP, while offset is the constant part.
 218  */
 219 void gen_fun::substitute(Matrix *CP, const mat_ZZ& map, const vec_ZZ& offset)
 220 {
 221     Polyhedron *C = Polyhedron_Image(context, CP, 0);
 222     Polyhedron_Free(context);
 223     context = C;
 224     for (int i = 0; i < term.size(); ++i) {
 225         term[i]->d.power *= map;
 226         term[i]->n.power *= map;
 227         for (int j = 0; j < term[i]->n.power.NumRows(); ++j)
 228             term[i]->n.power[j] += offset;
 229     }
 230 }
 231
 232 gen_fun *gen_fun::Hadamard_product(gen_fun *gf, unsigned MaxRays)
 233 {
 234     Polyhedron *C = DomainIntersection(context, gf->context, MaxRays);
 235     Polyhedron *U = Universe_Polyhedron(C->Dimension);
 236     gen_fun *sum = new gen_fun(C);
 237     for (int i = 0; i < term.size(); ++i) {
 238         for (int i2 = 0; i2 < gf->term.size(); ++i2) {
 239             int d = term[i]->d.power.NumCols();
 240             int k1 = term[i]->d.power.NumRows();
 241             int k2 = gf->term[i2]->d.power.NumRows();
 242             assert(term[i]->d.power.NumCols() == gf->term[i2]->d.power.NumCols());
 243             for (int j = 0; j < term[i]->n.power.NumRows(); ++j) {
 244                 for (int j2 = 0; j2 < gf->term[i2]->n.power.NumRows(); ++j2) {
 245                     Matrix *M = Matrix_Alloc(k1+k2+d+d, 1+k1+k2+d+1);
 246                     for (int k = 0; k < k1+k2; ++k) {
 247                         value_set_si(M->p[k][0], 1);
 248                         value_set_si(M->p[k][1+k], 1);
 249                     }
 250                     for (int k = 0; k < d; ++k) {
 251                         value_set_si(M->p[k1+k2+k][1+k1+k2+k], -1);
 252                         zz2value(term[i]->n.power[j][k], M->p[k1+k2+k][1+k1+k2+d]);
 253                         for (int l = 0; l < k1; ++l)
 254                             zz2value(term[i]->d.power[l][k], M->p[k1+k2+k][1+l]);
 255                     }
 256                     for (int k = 0; k < d; ++k) {
 257                         value_set_si(M->p[k1+k2+d+k][1+k1+k2+k], -1);
 258                         zz2value(gf->term[i2]->n.power[j2][k],
 259                                  M->p[k1+k2+d+k][1+k1+k2+d]);
 260                         for (int l = 0; l < k2; ++l)
 261                             zz2value(gf->term[i2]->d.power[l][k],
 262                                      M->p[k1+k2+d+k][1+k1+l]);
 263                     }
 264                     Polyhedron *P = Constraints2Polyhedron(M, MaxRays);
 265                     Matrix_Free(M);
 266
 267                     gen_fun *t = barvinok_series(P, U, MaxRays);
 268
 269                     ZZ cn = term[i]->n.coeff[j][0] * gf->term[i2]->n.coeff[j2][0];
 270                     ZZ cd = term[i]->n.coeff[j][1] * gf->term[i2]->n.coeff[j2][1];
 271                     sum->add(cn, cd, t);
 272                     delete t;
 273
 274                     Polyhedron_Free(P);
 275                 }
 276             }
 277         }
 278     }
 279     Polyhedron_Free(U);
 280     return sum;
 281 }
 282
 283 void gen_fun::add_union(gen_fun *gf, unsigned MaxRays)
 284 {
 285     ZZ one, mone;
 286     one = 1;
 287     mone = -1;
 288
 289     gen_fun *hp = Hadamard_product(gf, MaxRays);
 290     add(one, one, gf);
 291     add(mone, one, hp);
 292     delete hp;
 293 }
 294
 295 static void Polyhedron_Shift(Polyhedron *P, Vector *offset)
 296 {
 297     Value tmp;
 298     value_init(tmp);
 299     for (int i = 0; i < P->NbConstraints; ++i) {
 300         Inner_Product(P->Constraint[i]+1, offset->p, P->Dimension, &tmp);
 301         value_subtract(P->Constraint[i][1+P->Dimension],
 302                        P->Constraint[i][1+P->Dimension], tmp);
 303     }
 304     for (int i = 0; i < P->NbRays; ++i) {
 305         if (value_notone_p(P->Ray[i][0]))
 306             continue;
 307         if (value_zero_p(P->Ray[i][1+P->Dimension]))
 308             continue;
 309         Vector_Combine(P->Ray[i]+1, offset->p, P->Ray[i]+1,
 310                        P->Ray[i][0], P->Ray[i][1+P->Dimension], P->Dimension);
 311     }
 312     value_clear(tmp);
 313 }
 314
 315 void gen_fun::shift(const vec_ZZ& offset)
 316 {
 317     for (int i = 0; i < term.size(); ++i)
 318         for (int j = 0; j < term[i]->n.power.NumRows(); ++j)
 319             term[i]->n.power[j] += offset;
 320
 321     Vector *v = Vector_Alloc(offset.length());
 322     zz2values(offset, v->p);
 323     Polyhedron_Shift(context, v);
 324     Vector_Free(v);
 325 }
 326
 327 static void print_power(vec_ZZ& c, vec_ZZ& p,
 328                         unsigned int nparam, char **param_name)
 329 {
 330     bool first = true;
 331
 332     for (int i = 0; i < p.length(); ++i) {
 333         if (p[i] == 0)
 334             continue;
 335         if (first) {
 336             if (c[0] == -1 && c[1] == 1)
 337                 cout << "-";
 338             else if (c[0] != 1 || c[1] != 1) {
 339                 cout << c[0];
 340                 if (c[1] != 1)
 341                     cout << " / " << c[1];
 342                 cout << "*";
 343             }
 344             first = false;
 345         } else
 346             cout << "*";
 347         if (i < nparam)
 348             cout << param_name[i];
 349         else
 350             cout << "x" << i;
 351         if (p[i] == 1)
 352             continue;
 353         if (p[i] < 0)
 354             cout << "^(" << p[i] << ")";
 355         else
 356             cout << "^" << p[i];
 357     }
 358     if (first) {
 359         cout << c[0];
 360         if (c[1] != 1)
 361             cout << " / " << c[1];
 362     }
 363 }
 364
 365 void gen_fun::print(unsigned int nparam, char **param_name) const
 366 {
 367     vec_ZZ mone;
 368     mone.SetLength(2);
 369     mone[0] = -1;
 370     mone[1] = 1;
 371     for (int i = 0; i < term.size(); ++i) {
 372         if (i != 0)
 373             cout << " + ";
 374         cout << "(";
 375         for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
 376             if (j != 0 && term[i]->n.coeff[j][0] > 0)
 377                 cout << "+";
 378             print_power(term[i]->n.coeff[j], term[i]->n.power[j],
 379                         nparam, param_name);
 380         }
 381         cout << ")/(";
 382         for (int j = 0; j < term[i]->d.power.NumRows(); ++j) {
 383             if (j != 0)
 384                 cout << " * ";
 385             cout << "(1";
 386             print_power(mone, term[i]->d.power[j],
 387                         nparam, param_name);
 388             cout << ")";
 389         }
 390         cout << ")";
 391     }
 392 }
 393
 394 gen_fun::operator evalue *() const
 395 {
 396     evalue *EP = NULL;
 397     evalue factor;
 398     value_init(factor.d);
 399     value_init(factor.x.n);
 400     for (int i = 0; i < term.size(); ++i) {
 401         unsigned nvar = term[i]->d.power.NumRows();
 402         unsigned nparam = term[i]->d.power.NumCols();
 403         Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + nparam + 1);
 404         mat_ZZ& d = term[i]->d.power;
 405         Polyhedron *U = context ? context : Universe_Polyhedron(nparam);
 406
 407         for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
 408             for (int r = 0; r < nparam; ++r) {
 409                 value_set_si(C->p[r][0], 0);
 410                 for (int c = 0; c < nvar; ++c) {
 411                     zz2value(d[c][r], C->p[r][1+c]);
 412                 }
 413                 Vector_Set(&C->p[r][1+nvar], 0, nparam);
 414                 value_set_si(C->p[r][1+nvar+r], -1);
 415                 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar+nparam]);
 416             }
 417             for (int r = 0; r < nvar; ++r) {
 418                 value_set_si(C->p[nparam+r][0], 1);
 419                 Vector_Set(&C->p[nparam+r][1], 0, nvar + nparam + 1);
 420                 value_set_si(C->p[nparam+r][1+r], 1);
 421             }
 422             Polyhedron *P = Constraints2Polyhedron(C, 0);
 423             evalue *E = barvinok_enumerate_ev(P, U, 0);
 424             Polyhedron_Free(P);
 425             if (EVALUE_IS_ZERO(*E)) {
 426                 free_evalue_refs(E);
 427                 free(E);
 428                 continue;
 429             }
 430             zz2value(term[i]->n.coeff[j][0], factor.x.n);
 431             zz2value(term[i]->n.coeff[j][1], factor.d);
 432             emul(&factor, E);
 433             /*
 434             Matrix_Print(stdout, P_VALUE_FMT, C);
 435             char *test[] = { "A", "B", "C", "D", "E", "F", "G" };
 436             print_evalue(stdout, E, test);
 437             */
 438             if (!EP)
 439                 EP = E;
 440             else {
 441                 eadd(E, EP);
 442                 free_evalue_refs(E);
 443                 free(E);
 444             }
 445         }
 446         Matrix_Free(C);
 447         if (!context)
 448             Polyhedron_Free(U);
 449     }
 450     value_clear(factor.d);
 451     value_clear(factor.x.n);
 452     return EP;
 453 }
 454
 455 void gen_fun::coefficient(Value* params, Value* c) const
 456 {
 457     if (context && !in_domain(context, params)) {
 458         value_set_si(*c, 0);
 459         return;
 460     }
 461
 462     evalue part;
 463     value_init(part.d);
 464     value_init(part.x.n);
 465     evalue sum;
 466     value_init(sum.d);
 467     evalue_set_si(&sum, 0, 1);
 468     Value tmp;
 469     value_init(tmp);
 470
 471     for (int i = 0; i < term.size(); ++i) {
 472         unsigned nvar = term[i]->d.power.NumRows();
 473         unsigned nparam = term[i]->d.power.NumCols();
 474         Matrix *C = Matrix_Alloc(nparam + nvar, 1 + nvar + 1);
 475         mat_ZZ& d = term[i]->d.power;
 476
 477         for (int j = 0; j < term[i]->n.coeff.NumRows(); ++j) {
 478             for (int r = 0; r < nparam; ++r) {
 479                 value_set_si(C->p[r][0], 0);
 480                 for (int c = 0; c < nvar; ++c) {
 481                     zz2value(d[c][r], C->p[r][1+c]);
 482                 }
 483                 zz2value(term[i]->n.power[j][r], C->p[r][1+nvar]);
 484                 value_subtract(C->p[r][1+nvar], C->p[r][1+nvar], params[r]);
 485             }
 486             for (int r = 0; r < nvar; ++r) {
 487                 value_set_si(C->p[nparam+r][0], 1);
 488                 Vector_Set(&C->p[nparam+r][1], 0, nvar + 1);
 489                 value_set_si(C->p[nparam+r][1+r], 1);
 490             }
 491             Polyhedron *P = Constraints2Polyhedron(C, 0);
 492             if (emptyQ(P)) {
 493                 Polyhedron_Free(P);
 494                 continue;
 495             }
 496             barvinok_count(P, &tmp, 0);
 497             Polyhedron_Free(P);
 498             if (value_zero_p(tmp))
 499                 continue;
 500             zz2value(term[i]->n.coeff[j][0], part.x.n);
 501             zz2value(term[i]->n.coeff[j][1], part.d);
 502             value_multiply(part.x.n, part.x.n, tmp);
 503             eadd(&part, &sum);
 504         }
 505         Matrix_Free(C);
 506     }
 507
 508     assert(value_one_p(sum.d));
 509     value_assign(*c, sum.x.n);
 510
 511     value_clear(tmp);
 512     value_clear(part.d);
 513     value_clear(part.x.n);
 514     value_clear(sum.d);
 515     value_clear(sum.x.n);
 516 }