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
   8 static void rays(mat_ZZ& rays, Polyhedron *C)
   9 {
  10     unsigned dim = C->NbRays - 1; /* don't count zero vertex */
  11     assert(C->NbRays - 1 == C->Dimension);
  12     rays.SetDims(dim, dim);
  13
  14     int i, j;
  15     for (i = 0, j = 0; i < C->NbRays; ++i) {
  16         if (value_notzero_p(C->Ray[i][dim+1]))
  17             continue;
  18         values2zz(C->Ray[i]+1, rays[j], dim);
  19         ++j;
  20     }
  21 }
  22
  23 void np_base::handle(const signed_cone& sc)
  24 {
  25     assert(sc.C->NbRays-1 == dim);
  26     factor.n *= sc.sign;
  27     mat_ZZ r;
  28     rays(r, sc.C);
  29     handle(r, current_vertex, factor, sc.closed);
  30     factor.n *= sc.sign;
  31 }
  32
  33 void np_base::start(Polyhedron *P, barvinok_options *options)
  34 {
  35     QQ factor(1, 1);
  36     init(P);
  37     for (int i = 0; i < P->NbRays; ++i) {
  38         if (!value_pos_p(P->Ray[i][dim+1]))
  39             continue;
  40
  41         Polyhedron *C = supporting_cone(P, i);
  42         do_vertex_cone(factor, C, P->Ray[i]+1, options);
  43     }
  44 }
  45
  46 /* input:
  47  *      f: the powers in the denominator for the remaining vars
  48  *            each row refers to a factor
  49  *      den_s: for each factor, the power of  (s+1)
  50  *      sign
  51  *      num_s: powers in the numerator corresponding to the summed vars
  52  *      num_p: powers in the numerator corresponding to the remaining vars
  53  * number of rays in cone: "dim" = "k"
  54  * length of each ray: "dim" = "d"
  55  * for now, it is assumed: k == d
  56  * output:
  57  *      den_p: for each factor
  58  *              0: independent of remaining vars
  59  *              1: power corresponds to corresponding row in f
  60  *
  61  * all inputs are subject to change
  62  */
  63 void normalize(ZZ& sign, ZZ& num_s, vec_ZZ& num_p, vec_ZZ& den_s, vec_ZZ& den_p,
  64                mat_ZZ& f)
  65 {
  66     unsigned dim = f.NumRows();
  67     unsigned nparam = num_p.length();
  68     unsigned nvar = dim - nparam;
  69
  70     int change = 0;
  71
  72     for (int j = 0; j < den_s.length(); ++j) {
  73         if (den_s[j] == 0) {
  74             den_p[j] = 1;
  75             continue;
  76         }
  77         int k;
  78         for (k = 0; k < nparam; ++k)
  79             if (f[j][k] != 0)
  80                 break;
  81         if (k < nparam) {
  82             den_p[j] = 1;
  83             if (den_s[j] > 0) {
  84                 f[j] = -f[j];
  85                 num_p += 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             num_s += den_s[j];
  94         }
  95     }
  96
  97     if (change)
  98         sign = -sign;
  99 }
 100
 101 void reducer::reduce(QQ c, vec_ZZ& num, const mat_ZZ& den_f)
 102 {
 103     unsigned len = den_f.NumRows();  // number of factors in den
 104
 105     if (num.length() == lower) {
 106         base(c, num, den_f);
 107         return;
 108     }
 109     assert(num.length() > 1);
 110
 111     vec_ZZ den_s;
 112     mat_ZZ den_r;
 113     ZZ num_s;
 114     vec_ZZ num_p;
 115
 116     split(num, num_s, num_p, den_f, den_s, den_r);
 117
 118     vec_ZZ den_p;
 119     den_p.SetLength(len);
 120
 121     normalize(c.n, num_s, num_p, den_s, den_p, den_r);
 122
 123     int only_param = 0;     // k-r-s from text
 124     int no_param = 0;       // r from text
 125     for (int k = 0; k < len; ++k) {
 126         if (den_p[k] == 0)
 127             ++no_param;
 128         else if (den_s[k] == 0)
 129             ++only_param;
 130     }
 131     if (no_param == 0) {
 132         reduce(c, num_p, den_r);
 133     } else {
 134         int k, l;
 135         mat_ZZ pden;
 136         pden.SetDims(only_param, den_r.NumCols());
 137
 138         for (k = 0, l = 0; k < len; ++k)
 139             if (den_s[k] == 0)
 140                 pden[l++] = den_r[k];
 141
 142         for (k = 0; k < len; ++k)
 143             if (den_p[k] == 0)
 144                 break;
 145
 146         dpoly n(no_param, num_s);
 147         dpoly D(no_param, den_s[k], 1);
 148         for ( ; ++k < len; )
 149             if (den_p[k] == 0) {
 150                 dpoly fact(no_param, den_s[k], 1);
 151                 D *= fact;
 152             }
 153
 154         if (no_param + only_param == len) {
 155             mpq_set_si(tcount, 0, 1);
 156             n.div(D, tcount, one);
 157
 158             QQ q;
 159             value2zz(mpq_numref(tcount), q.n);
 160             value2zz(mpq_denref(tcount), q.d);
 161
 162             q *= c;
 163
 164             if (q.n != 0)
 165                 reduce(q, num_p, pden);
 166         } else {
 167             dpoly_r * r = 0;
 168
 169             for (k = 0; k < len; ++k) {
 170                 if (den_s[k] == 0 || den_p[k] == 0)
 171                     continue;
 172
 173                 dpoly pd(no_param-1, den_s[k], 1);
 174
 175                 int l;
 176                 for (l = 0; l < k; ++l)
 177                     if (den_r[l] == den_r[k])
 178                         break;
 179
 180                 if (r == 0)
 181                     r = new dpoly_r(n, pd, l, len);
 182                 else {
 183                     dpoly_r *nr = new dpoly_r(r, pd, l, len);
 184                     delete r;
 185                     r = nr;
 186                 }
 187             }
 188
 189             dpoly_r *rc = r->div(D);
 190
 191             QQ factor;
 192             factor.d = rc->denom * c.d;
 193
 194             int common = pden.NumRows();
 195             dpoly_r_term_list& final = rc->c[rc->len-1];
 196             int rows;
 197             dpoly_r_term_list::iterator j;
 198             for (j = final.begin(); j != final.end(); ++j) {
 199                 if ((*j)->coeff == 0)
 200                     continue;
 201                 rows = common;
 202                 pden.SetDims(rows, pden.NumCols());
 203                 for (int k = 0; k < rc->dim; ++k) {
 204                     int n = (*j)->powers[k];
 205                     if (n == 0)
 206                         continue;
 207                     pden.SetDims(rows+n, pden.NumCols());
 208                     for (int l = 0; l < n; ++l)
 209                         pden[rows+l] = den_r[k];
 210                     rows += n;
 211                 }
 212                 factor.n = (*j)->coeff *= c.n;
 213                 reduce(factor, num_p, pden);
 214             }
 215
 216             delete rc;
 217             delete r;
 218         }
 219     }
 220 }
 221
 222 void reducer::handle(const mat_ZZ& den, Value *V, QQ c, int *closed)
 223 {
 224     lattice_point(V, den, vertex, closed);
 225
 226     reduce(c, vertex, den);
 227 }
 228
 229 void ireducer::split(vec_ZZ& num, ZZ& num_s, vec_ZZ& num_p,
 230                      const mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r)
 231 {
 232     unsigned len = den_f.NumRows();  // number of factors in den
 233     unsigned d = num.length() - 1;
 234
 235     den_s.SetLength(len);
 236     den_r.SetDims(len, d);
 237
 238     for (int r = 0; r < len; ++r) {
 239         den_s[r] = den_f[r][0];
 240         for (int k = 1; k <= d; ++k)
 241             den_r[r][k-1] = den_f[r][k];
 242     }
 243
 244     num_s = num[0];
 245     num_p.SetLength(d);
 246     for (int k = 1 ; k <= d; ++k)
 247         num_p[k-1] = num[k];
 248 }
 249
 250 void normalize(ZZ& sign, ZZ& num, vec_ZZ& den)
 251 {
 252     unsigned dim = den.length();
 253
 254     int change = 0;
 255
 256     for (int j = 0; j < den.length(); ++j) {
 257         if (den[j] > 0)
 258             change ^= 1;
 259         else {
 260             den[j] = abs(den[j]);
 261             num += den[j];
 262         }
 263     }
 264     if (change)
 265         sign = -sign;
 266 }
 267
 268 void icounter::base(QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
 269 {
 270     int r;
 271     unsigned len = den_f.NumRows();  // number of factors in den
 272     vec_ZZ den_s;
 273     den_s.SetLength(len);
 274     ZZ num_s = num[0];
 275     for (r = 0; r < len; ++r)
 276         den_s[r] = den_f[r][0];
 277     normalize(c.n, num_s, den_s);
 278
 279     dpoly n(len, num_s);
 280     dpoly D(len, den_s[0], 1);
 281     for (int k = 1; k < len; ++k) {
 282         dpoly fact(len, den_s[k], 1);
 283         D *= fact;
 284     }
 285     mpq_set_si(tcount, 0, 1);
 286     n.div(D, tcount, one);
 287     zz2value(c.n, tn);
 288     zz2value(c.d, td);
 289     mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
 290     mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
 291     mpq_canonicalize(tcount);
 292     mpq_add(count, count, tcount);
 293 }
 294
 295 void infinite_icounter::base(QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
 296 {
 297     int r;
 298     unsigned len = den_f.NumRows();  // number of factors in den
 299     vec_ZZ den_s;
 300     den_s.SetLength(len);
 301     ZZ num_s = num[0];
 302
 303     for (r = 0; r < len; ++r)
 304         den_s[r] = den_f[r][0];
 305     normalize(c.n, num_s, den_s);
 306
 307     dpoly n(len, num_s);
 308     dpoly D(len, den_s[0], 1);
 309     for (int k = 1; k < len; ++k) {
 310         dpoly fact(len, den_s[k], 1);
 311         D *= fact;
 312     }
 313
 314     Value tmp;
 315     mpq_t factor;
 316     mpq_init(factor);
 317     value_init(tmp);
 318     zz2value(c.n, tmp);
 319     value_assign(mpq_numref(factor), tmp);
 320     zz2value(c.d, tmp);
 321     value_assign(mpq_denref(factor), tmp);
 322
 323     n.div(D, count, factor);
 324
 325     value_clear(tmp);
 326     mpq_clear(factor);
 327 }