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