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