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