counter.cc

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