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