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