counter.cc

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