counter.cc

   1 #include <assert.h>
   2 #include <barvinok/util.h>
   3 #include "counter.h"
   4 #include "lattice_point.h"
   5
   6 using std::cerr;
   7 using std::endl;
   8
   9 /* Computes the integer points in the fundamental parallelepiped and
  10  * passes them along (in num) to the counter specific (i.e., specialization
  11  * specific) add_lattice_points.
  12  */
  13 void counter_base::handle(const mat_ZZ& rays, Value *V, const QQ& c,
  14                           unsigned long det, barvinok_options *options)
  15 {
  16     Matrix* Rays = zz2matrix(rays);
  17
  18     assert(c.d == 1);
  19     assert(c.n == 1 || c.n == -1);
  20     int sign = to_int(c.n);
  21
  22     Matrix_Vector_Product(Rays, lambda->p, den->p_Init);
  23     for (int k = 0; k < dim; ++k)
  24         if (value_zero_p(den->p_Init[k])) {
  25             Matrix_Free(Rays);
  26             throw Orthogonal;
  27         }
  28     Inner_Product(lambda->p, V, dim, &tmp);
  29     lattice_points_fixed(V, &tmp, Rays, den, num, det);
  30     num->NbRows = det;
  31     Matrix_Free(Rays);
  32
  33     if (dim % 2)
  34         sign = -sign;
  35
  36     add_lattice_points(sign);
  37 }
  38
  39 static void add_falling_powers(dpoly& n, Value degree)
  40 {
  41     value_increment(n.coeff->p[0], n.coeff->p[0]);
  42     if (n.coeff->Size == 1)
  43         return;
  44
  45     int min = n.coeff->Size-1;
  46     if (value_posz_p(degree) && value_cmp_si(degree, min) < 0)
  47         min = VALUE_TO_INT(degree);
  48
  49     Value tmp;
  50     value_init(tmp);
  51     value_assign(tmp, degree);
  52     value_addto(n.coeff->p[1], n.coeff->p[1], tmp);
  53     for (int i = 2; i <= min; ++i) {
  54         value_decrement(degree, degree);
  55         value_multiply(tmp, tmp, degree);
  56         mpz_divexact_ui(tmp, tmp, i);
  57         value_addto(n.coeff->p[i], n.coeff->p[i], tmp);
  58     }
  59     value_clear(tmp);
  60 }
  61
  62 void counter::add_lattice_points(int sign)
  63 {
  64     dpoly d(dim);
  65     for (int k = 0; k < num->NbRows; ++k)
  66         add_falling_powers(d, num->p_Init[k]);
  67     dpoly n(dim, den->p_Init[0], 1);
  68     for (int k = 1; k < dim; ++k) {
  69         dpoly fact(dim, den->p_Init[k], 1);
  70         n *= fact;
  71     }
  72     d.div(n, count, sign);
  73 }
  74
  75
  76
  77
  78 void tcounter::setup_todd(unsigned dim)
  79 {
  80     value_set_si(todd.coeff->p[0], 1);
  81
  82     dpoly d(dim);
  83     value_set_si(d.coeff->p[dim], 1);
  84     for (int i = dim-1; i >= 0; --i)
  85         mpz_mul_ui(d.coeff->p[i], d.coeff->p[i+1], i+2);
  86
  87     todd_denom = todd.div(d);
  88     /* shift denominators up -> divide by (dim+1)! */
  89     for (int i = todd.coeff->Size-1; i >= 1; --i)
  90         value_assign(todd_denom->p[i], todd_denom->p[i-1]);
  91     value_set_si(todd_denom->p[0], 1);
  92 }
  93
  94 void tcounter::adapt_todd(dpoly& t, const Value c)
  95 {
  96     if (t.coeff->Size <= 1)
  97         return;
  98     value_assign(tmp, c);
  99     value_multiply(t.coeff->p[1], t.coeff->p[1], tmp);
 100     for (int i = 2; i < t.coeff->Size; ++i) {
 101         value_multiply(tmp, tmp, c);
 102         value_multiply(t.coeff->p[i], t.coeff->p[i], tmp);
 103     }
 104 }
 105
 106 void tcounter::add_powers(dpoly& n, const Value c)
 107 {
 108     value_increment(n.coeff->p[0], n.coeff->p[0]);
 109     if (n.coeff->Size == 1)
 110         return;
 111     value_assign(tmp, c);
 112     value_addto(n.coeff->p[1], n.coeff->p[1], tmp);
 113     for (int i = 2; i < n.coeff->Size; ++i) {
 114         value_multiply(tmp, tmp, c);
 115         value_addto(n.coeff->p[i], n.coeff->p[i], tmp);
 116     }
 117 }
 118
 119 void tcounter::add_lattice_points(int sign)
 120 {
 121     dpoly t(todd);
 122     value_assign(denom, den->p_Init[0]);
 123     adapt_todd(t, den->p_Init[0]);
 124     for (int k = 1; k < dim; ++k) {
 125         dpoly fact(todd);
 126         value_multiply(denom, denom, den->p_Init[k]);
 127         adapt_todd(fact, den->p_Init[k]);
 128         t *= fact;
 129     }
 130
 131     dpoly n(dim);
 132     for (int k = 0; k < num->NbRows; ++k)
 133         add_powers(n, num->p_Init[k]);
 134
 135     for (int i = 0; i < n.coeff->Size; ++i)
 136         value_multiply(n.coeff->p[i], n.coeff->p[i], todd_denom->p[i]);
 137     value_multiply(denom, denom, todd_denom->p[todd_denom->Size-1]);
 138
 139     value_set_si(tmp, 1);
 140     for (int i = 2; i < n.coeff->Size; ++i) {
 141         mpz_mul_ui(tmp, tmp, i);
 142         mpz_divexact(n.coeff->p[i], n.coeff->p[i], tmp);
 143     }
 144
 145     value_multiply(tmp, t.coeff->p[0], n.coeff->p[n.coeff->Size-1]);
 146     for (int i = 1; i < n.coeff->Size; ++i)
 147         value_addmul(tmp, t.coeff->p[i], n.coeff->p[n.coeff->Size-1-i]);
 148
 149     value_assign(mpq_numref(tcount), tmp);
 150     value_assign(mpq_denref(tcount), denom);
 151     mpq_canonicalize(tcount);
 152     if (sign == -1)
 153         mpq_sub(count, count, tcount);
 154     else
 155         mpq_add(count, count, tcount);
 156 }