reducer.cc

   1 #include <vector>
   2 #include <barvinok/util.h>
   3 #include "reducer.h"
   4 #include "lattice_point.h"
   5
   6 using std::vector;
   7 using std::cerr;
   8 using std::endl;
   9
  10 struct OrthogonalException Orthogonal;
  11
  12 void np_base::handle(const signed_cone& sc, barvinok_options *options)
  13 {
  14     assert(sc.rays.NumRows() == dim);
  15     factor.n *= sc.sign;
  16     handle(sc.rays, current_vertex, factor, sc.det, sc.closed, options);
  17     factor.n *= sc.sign;
  18 }
  19
  20 void np_base::start(Polyhedron *P, barvinok_options *options)
  21 {
  22     QQ factor(1, 1);
  23     for (;;) {
  24         try {
  25             init(P);
  26             for (int i = 0; i < P->NbRays; ++i) {
  27                 if (!value_pos_p(P->Ray[i][dim+1]))
  28                     continue;
  29
  30                 Polyhedron *C = supporting_cone(P, i);
  31                 do_vertex_cone(factor, C, P->Ray[i]+1, options);
  32             }
  33             break;
  34         } catch (OrthogonalException &e) {
  35             reset();
  36         }
  37     }
  38 }
  39
  40 /* input:
  41  *      f: the powers in the denominator for the remaining vars
  42  *            each row refers to a factor
  43  *      den_s: for each factor, the power of  (s+1)
  44  *      sign
  45  *      num_s: powers in the numerator corresponding to the summed vars
  46  *      num_p: powers in the numerator corresponding to the remaining vars
  47  * number of rays in cone: "dim" = "k"
  48  * length of each ray: "dim" = "d"
  49  * for now, it is assumed: k == d
  50  * output:
  51  *      den_p: for each factor
  52  *              0: independent of remaining vars
  53  *              1: power corresponds to corresponding row in f
  54  *
  55  * all inputs are subject to change
  56  */
  57 void normalize(ZZ& sign, vec_ZZ& num_s, mat_ZZ& num_p, vec_ZZ& den_s, vec_ZZ& den_p,
  58                mat_ZZ& f)
  59 {
  60     unsigned dim = f.NumRows();
  61     unsigned nparam = num_p.NumCols();
  62     unsigned nvar = dim - nparam;
  63
  64     int change = 0;
  65
  66     for (int j = 0; j < den_s.length(); ++j) {
  67         if (den_s[j] == 0) {
  68             den_p[j] = 1;
  69             continue;
  70         }
  71         int k;
  72         for (k = 0; k < nparam; ++k)
  73             if (f[j][k] != 0)
  74                 break;
  75         if (k < nparam) {
  76             den_p[j] = 1;
  77             if (den_s[j] > 0) {
  78                 f[j] = -f[j];
  79                 for (int i = 0; i < num_p.NumRows(); ++i)
  80                     num_p[i] += f[j];
  81             }
  82         } else
  83             den_p[j] = 0;
  84         if (den_s[j] > 0)
  85             change ^= 1;
  86         else {
  87             den_s[j] = abs(den_s[j]);
  88             for (int i = 0; i < num_p.NumRows(); ++i)
  89                 num_s[i] += den_s[j];
  90         }
  91     }
  92
  93     if (change)
  94         sign = -sign;
  95 }
  96
  97 void reducer::base(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den_f)
  98 {
  99     for (int i = 0; i < num.NumRows(); ++i)
 100         base(c[i], num[i], den_f);
 101 }
 102
 103 struct dpoly_r_scanner {
 104     const dpoly_r *rc;
 105     const dpoly * const *num;
 106     int n;
 107     int dim;
 108     dpoly_r_term_list::iterator *iter;
 109     vector<int> powers;
 110     vec_ZZ coeff;
 111
 112     dpoly_r_scanner(const dpoly * const *num, int n, const dpoly_r *rc, int dim)
 113                     : num(num), rc(rc), n(n), dim(dim), powers(dim, 0) {
 114         coeff.SetLength(n);
 115         iter = new dpoly_r_term_list::iterator[rc->len];
 116         for (int i = 0; i < rc->len; ++i) {
 117             int k;
 118             for (k = 0; k < n; ++k)
 119                 if (value_notzero_p(num[k]->coeff->p[rc->len-1-i]))
 120                     break;
 121             if (k < n)
 122                 iter[i] = rc->c[i].begin();
 123             else
 124                 iter[i] = rc->c[i].end();
 125         }
 126     }
 127     bool next() {
 128         int pos[rc->len];
 129         int len = 0;
 130
 131         for (int i = 0; i < rc->len; ++i) {
 132             if (iter[i] == rc->c[i].end())
 133                 continue;
 134             if (!len)
 135                 pos[len++] = i;
 136             else {
 137                 if ((*iter[i])->powers < (*iter[pos[0]])->powers) {
 138                     pos[0] = i;
 139                     len = 1;
 140                 } else if ((*iter[i])->powers == (*iter[pos[0]])->powers)
 141                     pos[len++] = i;
 142             }
 143         }
 144
 145         if (!len)
 146             return false;
 147
 148         powers = (*iter[pos[0]])->powers;
 149         for (int k = 0; k < n; ++k) {
 150             value2zz(num[k]->coeff->p[rc->len-1-pos[0]], tmp);
 151             mul(coeff[k], (*iter[pos[0]])->coeff, tmp);
 152         }
 153         ++iter[pos[0]];
 154         for (int i = 1; i < len; ++i) {
 155             for (int k = 0; k < n; ++k) {
 156                 value2zz(num[k]->coeff->p[rc->len-1-pos[i]], tmp);
 157                 mul(tmp, (*iter[pos[i]])->coeff, tmp);
 158                 add(coeff[k], coeff[k], tmp);
 159             }
 160             ++iter[pos[i]];
 161         }
 162
 163         return true;
 164     }
 165     ~dpoly_r_scanner() {
 166         delete [] iter;
 167     }
 168 private:
 169     ZZ tmp;
 170 };
 171
 172 void reducer::reduce(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den_f)
 173 {
 174     assert(c.length() == num.NumRows());
 175     unsigned len = den_f.NumRows();  // number of factors in den
 176     vec_QQ c2 = c;
 177
 178     if (num.NumCols() == lower) {
 179         base(c, num, den_f);
 180         return;
 181     }
 182     assert(num.NumCols() > 1);
 183     assert(num.NumRows() > 0);
 184
 185     vec_ZZ den_s;
 186     mat_ZZ den_r;
 187     vec_ZZ num_s;
 188     mat_ZZ num_p;
 189
 190     split(num, num_s, num_p, den_f, den_s, den_r);
 191
 192     vec_ZZ den_p;
 193     den_p.SetLength(len);
 194
 195     ZZ sign(INIT_VAL, 1);
 196     normalize(sign, num_s, num_p, den_s, den_p, den_r);
 197     c2 *= sign;
 198
 199     int only_param = 0;     // k-r-s from text
 200     int no_param = 0;       // r from text
 201     for (int k = 0; k < len; ++k) {
 202         if (den_p[k] == 0)
 203             ++no_param;
 204         else if (den_s[k] == 0)
 205             ++only_param;
 206     }
 207     if (no_param == 0) {
 208         reduce(c2, num_p, den_r);
 209     } else {
 210         int k, l;
 211         mat_ZZ pden;
 212         pden.SetDims(only_param, den_r.NumCols());
 213
 214         for (k = 0, l = 0; k < len; ++k)
 215             if (den_s[k] == 0)
 216                 pden[l++] = den_r[k];
 217
 218         for (k = 0; k < len; ++k)
 219             if (den_p[k] == 0)
 220                 break;
 221
 222         dpoly *n[num_s.length()];
 223         for (int i = 0; i < num_s.length(); ++i) {
 224             zz2value(num_s[i], tz);
 225             n[i] = new dpoly(no_param, tz);
 226             /* Search for other numerator (j) with same num_p.
 227              * If found, replace a[j]/b[j] * n[j] and a[i]/b[i] * n[i]
 228              * by 1/(b[j]*b[i]/g) * (a[j]*b[i]/g * n[j] + a[i]*b[j]/g * n[i])
 229              * where g = gcd(b[i], b[j].
 230              */
 231             for (int j = 0; j < i; ++j) {
 232                 if (num_p[i] != num_p[j])
 233                     continue;
 234                 ZZ g = GCD(c2[i].d, c2[j].d);
 235                 zz2value(c2[j].n * c2[i].d/g, tz);
 236                 *n[j] *= tz;
 237                 zz2value(c2[i].n * c2[j].d/g, tz);
 238                 *n[i] *= tz;
 239                 *n[j] += *n[i];
 240                 c2[j].n = 1;
 241                 c2[j].d *= c2[i].d/g;
 242                 delete n[i];
 243                 if (i < num_s.length()-1) {
 244                     num_s[i] = num_s[num_s.length()-1];
 245                     num_p[i] = num_p[num_s.length()-1];
 246                     c2[i] = c2[num_s.length()-1];
 247                 }
 248                 num_s.SetLength(num_s.length()-1);
 249                 c2.SetLength(c2.length()-1);
 250                 num_p.SetDims(num_p.NumRows()-1, num_p.NumCols());
 251                 --i;
 252                 break;
 253             }
 254         }
 255         zz2value(den_s[k], tz);
 256         dpoly D(no_param, tz, 1);
 257         for ( ; ++k < len; )
 258             if (den_p[k] == 0) {
 259                 zz2value(den_s[k], tz);
 260                 dpoly fact(no_param, tz, 1);
 261                 D *= fact;
 262             }
 263
 264         if (no_param + only_param == len) {
 265             vec_QQ q;
 266             q.SetLength(num_s.length());
 267             for (int i = 0; i < num_s.length(); ++i) {
 268                 mpq_set_si(tcount, 0, 1);
 269                 n[i]->div(D, tcount, 1);
 270
 271                 value2zz(mpq_numref(tcount), q[i].n);
 272                 value2zz(mpq_denref(tcount), q[i].d);
 273                 q[i] *= c2[i];
 274             }
 275             for (int i = q.length()-1; i >= 0; --i) {
 276                 if (q[i].n == 0) {
 277                     q[i] = q[q.length()-1];
 278                     num_p[i] = num_p[q.length()-1];
 279                     q.SetLength(q.length()-1);
 280                     num_p.SetDims(num_p.NumRows()-1, num_p.NumCols());
 281                 }
 282             }
 283
 284             if (q.length() != 0)
 285                 reduce(q, num_p, pden);
 286         } else {
 287             value_set_si(tz, 0);
 288             dpoly one(no_param, tz);
 289             dpoly_r *r = NULL;
 290
 291             for (k = 0; k < len; ++k) {
 292                 if (den_s[k] == 0 || den_p[k] == 0)
 293                     continue;
 294
 295                 zz2value(den_s[k], tz);
 296                 dpoly pd(no_param-1, tz, 1);
 297
 298                 int l;
 299                 for (l = 0; l < k; ++l)
 300                     if (den_r[l] == den_r[k])
 301                         break;
 302
 303                 if (!r)
 304                     r = new dpoly_r(one, pd, l, len);
 305                 else {
 306                     dpoly_r *nr = new dpoly_r(r, pd, l, len);
 307                     delete r;
 308                     r = nr;
 309                 }
 310             }
 311
 312             vec_QQ factor;
 313             factor.SetLength(c2.length());
 314             int common = pden.NumRows();
 315             dpoly_r *rc = r->div(D);
 316             for (int i = 0; i < num_s.length(); ++i) {
 317                 factor[i].d = c2[i].d;
 318                 factor[i].d *= rc->denom;
 319             }
 320
 321             dpoly_r_scanner scanner(n, num_s.length(), rc, len);
 322             int rows;
 323             while (scanner.next()) {
 324                 int i;
 325                 for (i = 0; i < num_s.length(); ++i)
 326                     if (scanner.coeff[i] != 0)
 327                         break;
 328                 if (i == num_s.length())
 329                     continue;
 330                 rows = common;
 331                 pden.SetDims(rows, pden.NumCols());
 332                 for (int k = 0; k < rc->dim; ++k) {
 333                     int n = scanner.powers[k];
 334                     if (n == 0)
 335                         continue;
 336                     pden.SetDims(rows+n, pden.NumCols());
 337                     for (int l = 0; l < n; ++l)
 338                         pden[rows+l] = den_r[k];
 339                     rows += n;
 340                 }
 341                 for (int i = 0; i < num_s.length(); ++i) {
 342                     factor[i].n = c2[i].n;
 343                     factor[i].n *= scanner.coeff[i];
 344                 }
 345                 reduce(factor, num_p, pden);
 346             }
 347
 348             delete rc;
 349             delete r;
 350         }
 351         for (int i = 0; i < num_s.length(); ++i)
 352             delete n[i];
 353     }
 354 }
 355
 356 void reducer::handle(const mat_ZZ& den, Value *V, const QQ& c, unsigned long det,
 357                      int *closed, barvinok_options *options)
 358 {
 359     vec_QQ vc;
 360
 361     Matrix *points = Matrix_Alloc(det, dim);
 362     Matrix* Rays = zz2matrix(den);
 363     lattice_points_fixed(V, V, Rays, Rays, points, det, closed);
 364     Matrix_Free(Rays);
 365     matrix2zz(points, vertex, points->NbRows, points->NbColumns);
 366     Matrix_Free(points);
 367
 368     vc.SetLength(vertex.NumRows());
 369     for (int i = 0; i < vc.length(); ++i)
 370         vc[i] = c;
 371
 372     reduce(vc, vertex, den);
 373 }
 374
 375 void split_one(const mat_ZZ& num, vec_ZZ& num_s, mat_ZZ& num_p,
 376                const mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r)
 377 {
 378     unsigned len = den_f.NumRows();  // number of factors in den
 379     unsigned d = num.NumCols() - 1;
 380
 381     den_s.SetLength(len);
 382     den_r.SetDims(len, d);
 383
 384     for (int r = 0; r < len; ++r) {
 385         den_s[r] = den_f[r][0];
 386         for (int k = 1; k <= d; ++k)
 387             den_r[r][k-1] = den_f[r][k];
 388     }
 389
 390     num_s.SetLength(num.NumRows());
 391     num_p.SetDims(num.NumRows(), d);
 392     for (int i = 0; i < num.NumRows(); ++i) {
 393         num_s[i] = num[i][0];
 394         for (int k = 1 ; k <= d; ++k)
 395             num_p[i][k-1] = num[i][k];
 396     }
 397 }
 398
 399 void normalize(ZZ& sign, ZZ& num, vec_ZZ& den)
 400 {
 401     unsigned dim = den.length();
 402
 403     int change = 0;
 404
 405     for (int j = 0; j < den.length(); ++j) {
 406         if (den[j] > 0)
 407             change ^= 1;
 408         else {
 409             den[j] = abs(den[j]);
 410             num += den[j];
 411         }
 412     }
 413     if (change)
 414         sign = -sign;
 415 }
 416
 417 void icounter::base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
 418 {
 419     int r;
 420     unsigned len = den_f.NumRows();  // number of factors in den
 421     vec_ZZ den_s;
 422     den_s.SetLength(len);
 423     assert(num.length() == 1);
 424     ZZ num_s = num[0];
 425     for (r = 0; r < len; ++r)
 426         den_s[r] = den_f[r][0];
 427     ZZ sign = ZZ(INIT_VAL, 1);
 428     normalize(sign, num_s, den_s);
 429
 430     zz2value(num_s, tz);
 431     dpoly n(len, tz);
 432     zz2value(den_s[0], tz);
 433     dpoly D(len, tz, 1);
 434     for (int k = 1; k < len; ++k) {
 435         zz2value(den_s[k], tz);
 436         dpoly fact(len, tz, 1);
 437         D *= fact;
 438     }
 439     mpq_set_si(tcount, 0, 1);
 440     n.div(D, tcount, 1);
 441     zz2value(c.n, tn);
 442     if (sign == -1)
 443         value_oppose(tn, tn);
 444     zz2value(c.d, td);
 445     mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
 446     mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
 447     mpq_canonicalize(tcount);
 448     mpq_add(count, count, tcount);
 449 }
 450
 451 void infinite_icounter::base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
 452 {
 453     int r;
 454     unsigned len = den_f.NumRows();  // number of factors in den
 455     vec_ZZ den_s;
 456     den_s.SetLength(len);
 457     assert(num.length() == 1);
 458     ZZ num_s = num[0];
 459
 460     for (r = 0; r < len; ++r)
 461         den_s[r] = den_f[r][0];
 462     ZZ sign = ZZ(INIT_VAL, 1);
 463     normalize(sign, num_s, den_s);
 464
 465     zz2value(num_s, tz);
 466     dpoly n(len, tz);
 467     zz2value(den_s[0], tz);
 468     dpoly D(len, tz, 1);
 469     for (int k = 1; k < len; ++k) {
 470         zz2value(den_s[k], tz);
 471         dpoly fact(len, tz, 1);
 472         D *= fact;
 473     }
 474
 475     Value tmp;
 476     mpq_t factor;
 477     mpq_init(factor);
 478     value_init(tmp);
 479     zz2value(c.n, tmp);
 480     if (sign == -1)
 481         value_oppose(tmp, tmp);
 482     value_assign(mpq_numref(factor), tmp);
 483     zz2value(c.d, tmp);
 484     value_assign(mpq_denref(factor), tmp);
 485
 486     n.div(D, count, factor);
 487
 488     value_clear(tmp);
 489     mpq_clear(factor);
 490 }