tcounter.cc

   1 #include <barvinok/util.h>
   2 #include "tcounter.h"
   3 #include "lattice_point.h"
   4
   5 using std::cerr;
   6 using std::endl;
   7
   8 void tcounter::setup_todd(unsigned dim)
   9 {
  10     value_set_si(todd.coeff->p[0], 1);
  11
  12     dpoly d(dim);
  13     value_set_si(d.coeff->p[dim], 1);
  14     for (int i = dim-1; i >= 0; --i)
  15         mpz_mul_ui(d.coeff->p[i], d.coeff->p[i+1], i+2);
  16
  17     todd_denom = todd.div(d);
  18     /* shift denominators up -> divide by (dim+1)! */
  19     for (int i = todd.coeff->Size-1; i >= 1; --i)
  20         value_assign(todd_denom->p[i], todd_denom->p[i-1]);
  21     value_set_si(todd_denom->p[0], 1);
  22 }
  23
  24 void tcounter::adapt_todd(dpoly& t, const Value c)
  25 {
  26     if (t.coeff->Size <= 1)
  27         return;
  28     value_assign(tmp, c);
  29     value_multiply(t.coeff->p[1], t.coeff->p[1], tmp);
  30     for (int i = 2; i < t.coeff->Size; ++i) {
  31         value_multiply(tmp, tmp, c);
  32         value_multiply(t.coeff->p[i], t.coeff->p[i], tmp);
  33     }
  34 }
  35
  36 void tcounter::add_powers(dpoly& n, const Value c)
  37 {
  38     value_increment(n.coeff->p[0], n.coeff->p[0]);
  39     if (n.coeff->Size == 1)
  40         return;
  41     value_assign(tmp, c);
  42     value_addto(n.coeff->p[1], n.coeff->p[1], tmp);
  43     for (int i = 2; i < n.coeff->Size; ++i) {
  44         value_multiply(tmp, tmp, c);
  45         value_addto(n.coeff->p[i], n.coeff->p[i], tmp);
  46     }
  47 }
  48
  49 void tcounter::handle(const mat_ZZ& rays, Value *V, const QQ& c,
  50                       unsigned long det, barvinok_options *options)
  51 {
  52     Matrix* Rays = zz2matrix(rays);
  53
  54     assert(c.d == 1);
  55     assert(c.n == 1 || c.n == -1);
  56     sign = c.n;
  57
  58     if (dim % 2)
  59         sign = -sign;
  60
  61     Matrix_Vector_Product(Rays, lambda->p, den->p_Init);
  62     for (int k = 0; k < dim; ++k)
  63         if (value_zero_p(den->p_Init[k])) {
  64             Matrix_Free(Rays);
  65             throw Orthogonal;
  66         }
  67     Inner_Product(lambda->p, V, dim, &tmp);
  68     lattice_points_fixed(V, &tmp, Rays, den, num, det);
  69     num->NbRows = det;
  70     Matrix_Free(Rays);
  71
  72     dpoly t(todd);
  73     value_assign(denom, den->p_Init[0]);
  74     adapt_todd(t, den->p_Init[0]);
  75     for (int k = 1; k < dim; ++k) {
  76         dpoly fact(todd);
  77         value_multiply(denom, denom, den->p_Init[k]);
  78         adapt_todd(fact, den->p_Init[k]);
  79         t *= fact;
  80     }
  81
  82     dpoly n(dim);
  83     for (int k = 0; k < num->NbRows; ++k)
  84         add_powers(n, num->p_Init[k]);
  85
  86     for (int i = 0; i < n.coeff->Size; ++i)
  87         value_multiply(n.coeff->p[i], n.coeff->p[i], todd_denom->p[i]);
  88     value_multiply(denom, denom, todd_denom->p[todd_denom->Size-1]);
  89
  90     value_set_si(tmp, 1);
  91     for (int i = 2; i < n.coeff->Size; ++i) {
  92         mpz_mul_ui(tmp, tmp, i);
  93         mpz_divexact(n.coeff->p[i], n.coeff->p[i], tmp);
  94     }
  95
  96     value_multiply(tmp, t.coeff->p[0], n.coeff->p[n.coeff->Size-1]);
  97     for (int i = 1; i < n.coeff->Size; ++i)
  98         value_addmul(tmp, t.coeff->p[i], n.coeff->p[n.coeff->Size-1-i]);
  99
 100     value_assign(mpq_numref(tcount), tmp);
 101     value_assign(mpq_denref(tcount), denom);
 102     mpq_canonicalize(tcount);
 103     if (sign == -1)
 104         mpq_sub(count, count, tcount);
 105     else
 106         mpq_add(count, count, tcount);
 107 }