counter.cc

   1 #include <assert.h>
   2 #include "counter.h"
   3 #include "lattice_point.h"
   4
   5 void counter::add_falling_powers(dpoly& n, Value degree)
   6 {
   7     value_increment(n.coeff->p[0], n.coeff->p[0]);
   8     if (n.coeff->Size == 1)
   9         return;
  10
  11     int min = n.coeff->Size-1;
  12     if (value_posz_p(degree) && value_cmp_si(degree, min) < 0)
  13         min = VALUE_TO_INT(degree);
  14
  15     Value tmp;
  16     value_init(tmp);
  17     value_assign(tmp, degree);
  18     value_addto(n.coeff->p[1], n.coeff->p[1], tmp);
  19     for (int i = 2; i <= min; ++i) {
  20         value_decrement(degree, degree);
  21         value_multiply(tmp, tmp, degree);
  22         mpz_divexact_ui(tmp, tmp, i);
  23         value_addto(n.coeff->p[i], n.coeff->p[i], tmp);
  24     }
  25     value_clear(tmp);
  26 }
  27
  28 void counter::handle(const mat_ZZ& rays, Value *V, const QQ& c, unsigned long det,
  29                      barvinok_options *options)
  30 {
  31     Matrix* Rays = zz2matrix(rays);
  32
  33     assert(c.d == 1);
  34     assert(c.n == 1 || c.n == -1);
  35     int sign = to_int(c.n);
  36
  37     Matrix_Vector_Product(Rays, lambda->p, den->p_Init);
  38     for (int k = 0; k < dim; ++k)
  39         if (value_zero_p(den->p_Init[k])) {
  40             Matrix_Free(Rays);
  41             throw Orthogonal;
  42         }
  43     Inner_Product(lambda->p, V, dim, &tmp);
  44     lattice_points_fixed(V, &tmp, Rays, den, num, det);
  45     num->NbRows = det;
  46     Matrix_Free(Rays);
  47
  48     if (dim % 2)
  49         sign = -sign;
  50
  51     dpoly d(dim);
  52     for (int k = 0; k < num->NbRows; ++k)
  53         add_falling_powers(d, num->p_Init[k]);
  54     dpoly n(dim, den->p_Init[0], 1);
  55     for (int k = 1; k < dim; ++k) {
  56         dpoly fact(dim, den->p_Init[k], 1);
  57         n *= fact;
  58     }
  59     d.div(n, count, sign);
  60 }