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