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, unsigned MaxRays)
  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, MaxRays);
  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             vector< dpoly_r_term * >& final = rc->c[rc->len-1];
 179             int rows;
 180             for (int j = 0; j < final.size(); ++j) {
 181                 if (final[j]->coeff == 0)
 182                     continue;
 183                 rows = common;
 184                 pden.SetDims(rows, pden.NumCols());
 185                 for (int k = 0; k < rc->dim; ++k) {
 186                     int n = final[j]->powers[k];
 187                     if (n == 0)
 188                         continue;
 189                     pden.SetDims(rows+n, pden.NumCols());
 190                     for (int l = 0; l < n; ++l)
 191                         pden[rows+l] = den_r[k];
 192                     rows += n;
 193                 }
 194                 factor.n = final[j]->coeff *= c.n;
 195                 reduce(factor, num_p, pden);
 196             }
 197
 198             delete rc;
 199             delete r;
 200         }
 201     }
 202 }
 203
 204 void reducer::handle_polar(Polyhedron *C, Value *V, QQ c)
 205 {
 206     lattice_point(V, C, vertex);
 207
 208     mat_ZZ den;
 209     den.SetDims(dim, dim);
 210
 211     int r;
 212     for (r = 0; r < dim; ++r)
 213         values2zz(C->Ray[r]+1, den[r], dim);
 214
 215     reduce(c, vertex, den);
 216 }
 217
 218 void ireducer::split(vec_ZZ& num, ZZ& num_s, vec_ZZ& num_p,
 219                      mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r)
 220 {
 221     unsigned len = den_f.NumRows();  // number of factors in den
 222     unsigned d = num.length() - 1;
 223
 224     den_s.SetLength(len);
 225     den_r.SetDims(len, d);
 226
 227     for (int r = 0; r < len; ++r) {
 228         den_s[r] = den_f[r][0];
 229         for (int k = 1; k <= d; ++k)
 230             den_r[r][k-1] = den_f[r][k];
 231     }
 232
 233     num_s = num[0];
 234     num_p.SetLength(d);
 235     for (int k = 1 ; k <= d; ++k)
 236         num_p[k-1] = num[k];
 237 }
 238
 239 void normalize(ZZ& sign, ZZ& num, vec_ZZ& den)
 240 {
 241     unsigned dim = den.length();
 242
 243     int change = 0;
 244
 245     for (int j = 0; j < den.length(); ++j) {
 246         if (den[j] > 0)
 247             change ^= 1;
 248         else {
 249             den[j] = abs(den[j]);
 250             num += den[j];
 251         }
 252     }
 253     if (change)
 254         sign = -sign;
 255 }
 256
 257 void icounter::base(QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
 258 {
 259     int r;
 260     unsigned len = den_f.NumRows();  // number of factors in den
 261     vec_ZZ den_s;
 262     den_s.SetLength(len);
 263     ZZ num_s = num[0];
 264     for (r = 0; r < len; ++r)
 265         den_s[r] = den_f[r][0];
 266     normalize(c.n, num_s, den_s);
 267
 268     dpoly n(len, num_s);
 269     dpoly D(len, den_s[0], 1);
 270     for (int k = 1; k < len; ++k) {
 271         dpoly fact(len, den_s[k], 1);
 272         D *= fact;
 273     }
 274     mpq_set_si(tcount, 0, 1);
 275     n.div(D, tcount, one);
 276     zz2value(c.n, tn);
 277     zz2value(c.d, td);
 278     mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
 279     mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
 280     mpq_canonicalize(tcount);
 281     mpq_add(count, count, tcount);
 282 }
 283
 284 void infinite_icounter::base(QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
 285 {
 286     int r;
 287     unsigned len = den_f.NumRows();  // number of factors in den
 288     vec_ZZ den_s;
 289     den_s.SetLength(len);
 290     ZZ num_s = num[0];
 291
 292     for (r = 0; r < len; ++r)
 293         den_s[r] = den_f[r][0];
 294     normalize(c.n, num_s, den_s);
 295
 296     dpoly n(len, num_s);
 297     dpoly D(len, den_s[0], 1);
 298     for (int k = 1; k < len; ++k) {
 299         dpoly fact(len, den_s[k], 1);
 300         D *= fact;
 301     }
 302
 303     Value tmp;
 304     mpq_t factor;
 305     mpq_init(factor);
 306     value_init(tmp);
 307     zz2value(c.n, tmp);
 308     value_assign(mpq_numref(factor), tmp);
 309     zz2value(c.d, tmp);
 310     value_assign(mpq_denref(factor), tmp);
 311
 312     n.div(D, count, factor);
 313
 314     value_clear(tmp);
 315     mpq_clear(factor);
 316 }