counter.cc

   1 #include <assert.h>
   2 #include <barvinok/util.h>
   3 #include "counter.h"
   4 #include "lattice_point.h"
   5
   6 /* Computes the integer points in the fundamental parallelepiped and
   7  * passes them along (in num) to the counter specific (i.e., specialization
   8  * specific) add_lattice_points.
   9  */
  10 void counter_base::handle(const mat_ZZ& rays, Value *V, const QQ& c,
  11                           unsigned long det, barvinok_options *options)
  12 {
  13     Matrix* Rays = zz2matrix(rays);
  14
  15     assert(c.d == 1);
  16     assert(c.n == 1 || c.n == -1);
  17     int sign = to_int(c.n);
  18
  19     Matrix_Vector_Product(Rays, lambda->p, den->p_Init);
  20     for (int k = 0; k < dim; ++k)
  21         if (value_zero_p(den->p_Init[k])) {
  22             Matrix_Free(Rays);
  23             throw Orthogonal;
  24         }
  25     Inner_Product(lambda->p, V, dim, &tmp);
  26     lattice_points_fixed(V, &tmp, Rays, den, num, det);
  27     num->NbRows = det;
  28     Matrix_Free(Rays);
  29
  30     if (dim % 2)
  31         sign = -sign;
  32
  33     add_lattice_points(sign);
  34 }
  35
  36 static void add_falling_powers(dpoly& n, Value degree)
  37 {
  38     value_increment(n.coeff->p[0], n.coeff->p[0]);
  39     if (n.coeff->Size == 1)
  40         return;
  41
  42     int min = n.coeff->Size-1;
  43     if (value_posz_p(degree) && value_cmp_si(degree, min) < 0)
  44         min = VALUE_TO_INT(degree);
  45
  46     Value tmp;
  47     value_init(tmp);
  48     value_assign(tmp, degree);
  49     value_addto(n.coeff->p[1], n.coeff->p[1], tmp);
  50     for (int i = 2; i <= min; ++i) {
  51         value_decrement(degree, degree);
  52         value_multiply(tmp, tmp, degree);
  53         mpz_divexact_ui(tmp, tmp, i);
  54         value_addto(n.coeff->p[i], n.coeff->p[i], tmp);
  55     }
  56     value_clear(tmp);
  57 }
  58
  59 void counter::add_lattice_points(int sign)
  60 {
  61     dpoly d(dim);
  62     for (int k = 0; k < num->NbRows; ++k)
  63         add_falling_powers(d, num->p_Init[k]);
  64     dpoly n(dim, den->p_Init[0], 1);
  65     for (int k = 1; k < dim; ++k) {
  66         dpoly fact(dim, den->p_Init[k], 1);
  67         n *= fact;
  68     }
  69     d.div(n, count, sign);
  70 }
  71
  72
  73
  74
  75 void tcounter::setup_todd(unsigned dim)
  76 {
  77     value_set_si(todd.coeff->p[0], 1);
  78
  79     dpoly d(dim);
  80     value_set_si(d.coeff->p[dim], 1);
  81     for (int i = dim-1; i >= 0; --i)
  82         mpz_mul_ui(d.coeff->p[i], d.coeff->p[i+1], i+2);
  83
  84     todd_denom = todd.div(d);
  85     /* shift denominators up -> divide by (dim+1)! */
  86     for (int i = todd.coeff->Size-1; i >= 1; --i)
  87         value_assign(todd_denom->p[i], todd_denom->p[i-1]);
  88     value_set_si(todd_denom->p[0], 1);
  89 }
  90
  91 void tcounter::adapt_todd(dpoly& t, const Value c)
  92 {
  93     if (t.coeff->Size <= 1)
  94         return;
  95     value_assign(tmp, c);
  96     value_multiply(t.coeff->p[1], t.coeff->p[1], tmp);
  97     for (int i = 2; i < t.coeff->Size; ++i) {
  98         value_multiply(tmp, tmp, c);
  99         value_multiply(t.coeff->p[i], t.coeff->p[i], tmp);
 100     }
 101 }
 102
 103 void tcounter::add_powers(dpoly& n, const Value c)
 104 {
 105     value_increment(n.coeff->p[0], n.coeff->p[0]);
 106     if (n.coeff->Size == 1)
 107         return;
 108     value_assign(tmp, c);
 109     value_addto(n.coeff->p[1], n.coeff->p[1], tmp);
 110     for (int i = 2; i < n.coeff->Size; ++i) {
 111         value_multiply(tmp, tmp, c);
 112         value_addto(n.coeff->p[i], n.coeff->p[i], tmp);
 113     }
 114 }
 115
 116 void tcounter::add_lattice_points(int sign)
 117 {
 118     dpoly t(todd);
 119     value_assign(denom, den->p_Init[0]);
 120     adapt_todd(t, den->p_Init[0]);
 121     for (int k = 1; k < dim; ++k) {
 122         dpoly fact(todd);
 123         value_multiply(denom, denom, den->p_Init[k]);
 124         adapt_todd(fact, den->p_Init[k]);
 125         t *= fact;
 126     }
 127
 128     dpoly n(dim);
 129     for (int k = 0; k < num->NbRows; ++k)
 130         add_powers(n, num->p_Init[k]);
 131
 132     for (int i = 0; i < n.coeff->Size; ++i)
 133         value_multiply(n.coeff->p[i], n.coeff->p[i], todd_denom->p[i]);
 134     value_multiply(denom, denom, todd_denom->p[todd_denom->Size-1]);
 135
 136     value_set_si(tmp, 1);
 137     for (int i = 2; i < n.coeff->Size; ++i) {
 138         mpz_mul_ui(tmp, tmp, i);
 139         mpz_divexact(n.coeff->p[i], n.coeff->p[i], tmp);
 140     }
 141
 142     value_multiply(tmp, t.coeff->p[0], n.coeff->p[n.coeff->Size-1]);
 143     for (int i = 1; i < n.coeff->Size; ++i)
 144         value_addmul(tmp, t.coeff->p[i], n.coeff->p[n.coeff->Size-1-i]);
 145
 146     value_assign(mpq_numref(tcount), tmp);
 147     value_assign(mpq_denref(tcount), denom);
 148     mpq_canonicalize(tcount);
 149     if (sign == -1)
 150         mpq_sub(count, count, tcount);
 151     else
 152         mpq_add(count, count, tcount);
 153 }
 154
 155
 156 /*
 157  * Set lambda to a random vector that has a positive inner product
 158  * with all the rays of the context { x | A x + b >= 0 }.
 159  *
 160  * To do so, we take d rows A' from the constraint matrix A.
 161  * For every ray, we have
 162  *              A' r >= 0
 163  * We compute a random positive row vector x' and set x = x' A'.
 164  * We then have, for each ray r,
 165  *              x r = x' A' r >= 0
 166  * Although we can take any d rows from A, we choose linearly
 167  * independent rows from A to avoid the elements of the transformed
 168  * random vector to have simple linear relations, which would
 169  * increase the risk of the vector being orthogonal to one of
 170  * powers in the denominator of one of the terms in the generating
 171  * function.
 172  */
 173 void infinite_counter::init(Polyhedron *context, int n_try)
 174 {
 175     Matrix *M, *H, *Q, *U;
 176     mat_ZZ A;
 177
 178     randomvector(context, lambda, context->Dimension, n_try);
 179
 180     M = Matrix_Alloc(context->NbConstraints, context->Dimension);
 181     for (int i = 0; i < context->NbConstraints; ++i)
 182         Vector_Copy(context->Constraint[i]+1, M->p[i], context->Dimension);
 183     left_hermite(M, &H, &Q, &U);
 184     Matrix_Free(Q);
 185     Matrix_Free(U);
 186
 187     for (int col = 0, row = 0; col < H->NbColumns; ++col, ++row) {
 188         for (; row < H->NbRows; ++row)
 189             if (value_notzero_p(H->p[row][col]))
 190                 break;
 191         assert(row < H->NbRows);
 192         Vector_Copy(M->p[row], M->p[col], M->NbColumns);
 193     }
 194     matrix2zz(M, A, context->Dimension, context->Dimension);
 195     Matrix_Free(H);
 196     Matrix_Free(M);
 197
 198     for (int i = 0; i < lambda.length(); ++i)
 199         lambda[i] = abs(lambda[i]);
 200     lambda = lambda * A;
 201 }
 202
 203 static ZZ LCM(const ZZ& a, const ZZ& b)
 204 {
 205     return a * b / GCD(a, b);
 206 }
 207
 208 /* Normalize the powers in the denominator to be positive
 209  * and return -1 is the sign has to be changed.
 210  */
 211 static int normalized_sign(vec_ZZ& num, vec_ZZ& den)
 212 {
 213     int change = 0;
 214
 215     for (int j = 0; j < den.length(); ++j) {
 216         if (den[j] > 0)
 217             change ^= 1;
 218         else {
 219             den[j] = abs(den[j]);
 220             for (int k = 0; k < num.length(); ++k)
 221                 num[k] += den[j];
 222         }
 223     }
 224     return change ? -1 : 1;
 225 }
 226
 227 void infinite_counter::reduce(const vec_QQ& c, const mat_ZZ& num,
 228                               const mat_ZZ& den_f)
 229 {
 230     mpq_t factor;
 231     mpq_init(factor);
 232     unsigned len = den_f.NumRows();
 233
 234     ZZ lcm = c[0].d;
 235     for (int i = 1; i < c.length(); ++i)
 236         lcm = LCM(lcm, c[i].d);
 237
 238     vec_ZZ coeff;
 239     coeff.SetLength(c.length());
 240     for (int i = 0; i < c.length(); ++i)
 241         coeff[i] = c[i].n * lcm/c[i].d;
 242
 243     if (len == 0) {
 244         for (int i = 0; i < c.length(); ++i) {
 245             zz2value(coeff[i], tz);
 246             value_addto(mpq_numref(factor), mpq_numref(factor), tz);
 247         }
 248         zz2value(lcm, tz);
 249         value_assign(mpq_denref(factor), tz);
 250         mpq_add(count[0], count[0], factor);
 251         mpq_clear(factor);
 252         return;
 253     }
 254
 255     vec_ZZ num_s = num * lambda;
 256     vec_ZZ den_s = den_f * lambda;
 257     for (int i = 0; i < den_s.length(); ++i)
 258         assert(den_s[i] != 0);
 259     int sign = normalized_sign(num_s, den_s);
 260     if (sign < 0)
 261         coeff = -coeff;
 262
 263     dpoly n(len);
 264     zz2value(num_s[0], tz);
 265     add_falling_powers(n, tz);
 266     zz2value(coeff[0], tz);
 267     n *= tz;
 268     for (int i = 1; i < c.length(); ++i) {
 269         dpoly t(len);
 270         zz2value(num_s[i], tz);
 271         add_falling_powers(t, tz);
 272         zz2value(coeff[i], tz);
 273         t *= tz;
 274         n += t;
 275     }
 276     zz2value(den_s[0], tz);
 277     dpoly d(len, tz, 1);
 278     for (int i = 1; i < len; ++i) {
 279         zz2value(den_s[i], tz);
 280         dpoly fact(len, tz, 1);
 281         d *= fact;
 282     }
 283     value_set_si(mpq_numref(factor), 1);
 284     zz2value(lcm, tz);
 285     value_assign(mpq_denref(factor), tz);
 286     n.div(d, count, factor);
 287
 288     mpq_clear(factor);
 289 }