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(const signed_cone& sc, barvinok_options *options)
   9 {
  10     assert(sc.rays.NumRows() == dim);
  11     factor.n *= sc.sign;
  12     handle(sc.rays, current_vertex, factor, sc.closed, options);
  13     factor.n *= sc.sign;
  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, const 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(const mat_ZZ& den, Value *V, QQ c, int *closed,
 206                      barvinok_options *options)
 207 {
 208     lattice_point(V, den, vertex, closed);
 209
 210     reduce(c, vertex, den);
 211 }
 212
 213 void ireducer::split(vec_ZZ& num, ZZ& num_s, vec_ZZ& num_p,
 214                      const mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r)
 215 {
 216     unsigned len = den_f.NumRows();  // number of factors in den
 217     unsigned d = num.length() - 1;
 218
 219     den_s.SetLength(len);
 220     den_r.SetDims(len, d);
 221
 222     for (int r = 0; r < len; ++r) {
 223         den_s[r] = den_f[r][0];
 224         for (int k = 1; k <= d; ++k)
 225             den_r[r][k-1] = den_f[r][k];
 226     }
 227
 228     num_s = num[0];
 229     num_p.SetLength(d);
 230     for (int k = 1 ; k <= d; ++k)
 231         num_p[k-1] = num[k];
 232 }
 233
 234 void normalize(ZZ& sign, ZZ& num, vec_ZZ& den)
 235 {
 236     unsigned dim = den.length();
 237
 238     int change = 0;
 239
 240     for (int j = 0; j < den.length(); ++j) {
 241         if (den[j] > 0)
 242             change ^= 1;
 243         else {
 244             den[j] = abs(den[j]);
 245             num += den[j];
 246         }
 247     }
 248     if (change)
 249         sign = -sign;
 250 }
 251
 252 void icounter::base(QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
 253 {
 254     int r;
 255     unsigned len = den_f.NumRows();  // number of factors in den
 256     vec_ZZ den_s;
 257     den_s.SetLength(len);
 258     ZZ num_s = num[0];
 259     for (r = 0; r < len; ++r)
 260         den_s[r] = den_f[r][0];
 261     normalize(c.n, num_s, den_s);
 262
 263     dpoly n(len, num_s);
 264     dpoly D(len, den_s[0], 1);
 265     for (int k = 1; k < len; ++k) {
 266         dpoly fact(len, den_s[k], 1);
 267         D *= fact;
 268     }
 269     mpq_set_si(tcount, 0, 1);
 270     n.div(D, tcount, one);
 271     zz2value(c.n, tn);
 272     zz2value(c.d, td);
 273     mpz_mul(mpq_numref(tcount), mpq_numref(tcount), tn);
 274     mpz_mul(mpq_denref(tcount), mpq_denref(tcount), td);
 275     mpq_canonicalize(tcount);
 276     mpq_add(count, count, tcount);
 277 }
 278
 279 void infinite_icounter::base(QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
 280 {
 281     int r;
 282     unsigned len = den_f.NumRows();  // number of factors in den
 283     vec_ZZ den_s;
 284     den_s.SetLength(len);
 285     ZZ num_s = num[0];
 286
 287     for (r = 0; r < len; ++r)
 288         den_s[r] = den_f[r][0];
 289     normalize(c.n, num_s, den_s);
 290
 291     dpoly n(len, num_s);
 292     dpoly D(len, den_s[0], 1);
 293     for (int k = 1; k < len; ++k) {
 294         dpoly fact(len, den_s[k], 1);
 295         D *= fact;
 296     }
 297
 298     Value tmp;
 299     mpq_t factor;
 300     mpq_init(factor);
 301     value_init(tmp);
 302     zz2value(c.n, tmp);
 303     value_assign(mpq_numref(factor), tmp);
 304     zz2value(c.d, tmp);
 305     value_assign(mpq_denref(factor), tmp);
 306
 307     n.div(D, count, factor);
 308
 309     value_clear(tmp);
 310     mpq_clear(factor);
 311 }