reducer.cc

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