genfun_constructor.cc

   1 #include <NTL/vec_ZZ.h>
   2 #include <NTL/mat_ZZ.h>
   3 #include <barvinok/NTL_QQ.h>
   4 #include <barvinok/polylib.h>
   5 #include "genfun_constructor.h"
   6
   7 gf_base *gf_base::create(Polyhedron *context, unsigned dim, unsigned nparam,
   8                          barvinok_options *options)
   9 {
  10     gf_base *red;
  11     if (options->incremental_specialization == 2)
  12         red = new partial_bfcounter(context, dim, nparam);
  13     else if (options->incremental_specialization == 1)
  14         red = new partial_ireducer(context, dim, nparam);
  15     else
  16         red = new partial_reducer(context, dim, nparam);
  17     return red;
  18 }
  19
  20 void partial_ireducer::base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
  21 {
  22     gf->add(c, num, den_f);
  23 }
  24
  25 void partial_reducer::split(const mat_ZZ& num, vec_ZZ& num_s, mat_ZZ& num_p,
  26                             const mat_ZZ& den_f, vec_ZZ& den_s, mat_ZZ& den_r)
  27 {
  28     unsigned len = den_f.NumRows();  // number of factors in den
  29     unsigned nvar = tmp.length();
  30
  31     den_s.SetLength(len);
  32     den_r.SetDims(len, lower);
  33
  34     for (int r = 0; r < len; ++r) {
  35         bool all_zero;
  36
  37         for (int k = 0; k < nvar; ++k)
  38             tmp[k] = den_f[r][k];
  39         den_s[r] = tmp * lambda;
  40         all_zero = den_s[r] == 0;
  41
  42         for (int k = nvar; k < dim; ++k) {
  43             den_r[r][k-nvar] = den_f[r][k];
  44             if (all_zero)
  45                 all_zero = den_r[r][k-nvar] == 0;
  46         }
  47         if (all_zero)
  48             throw Orthogonal;
  49     }
  50
  51     num_s.SetLength(num.NumRows());
  52     num_p.SetDims(num.NumRows(), lower);
  53     for (int i = 0; i < num.NumRows(); ++i) {
  54         for (int k = 0; k < nvar; ++k)
  55             tmp[k] = num[i][k];
  56         num_s[i] = tmp * lambda;
  57         for (int k = nvar ; k < dim; ++k)
  58             num_p[i][k-nvar] = num[i][k];
  59     }
  60 }
  61
  62 void partial_reducer::base(const QQ& c, const vec_ZZ& num, const mat_ZZ& den_f)
  63 {
  64     gf->add(c, num, den_f);
  65 }
  66
  67 void partial_reducer::base(const vec_QQ& c, const mat_ZZ& num, const mat_ZZ& den_f)
  68 {
  69     gf->add(new short_rat(c, num, den_f));
  70 }
  71
  72 void partial_bfcounter::base(mat_ZZ& factors, bfc_vec& v)
  73 {
  74     mat_ZZ den;
  75     unsigned nf = factors.NumRows();
  76
  77     for (int i = 0; i < v.size(); ++i) {
  78         bfc_term* bfct = static_cast<bfc_term *>(v[i]);
  79         den.SetDims(0, lower);
  80         int total_power = 0;
  81         int p = 0;
  82         for (int j = 0; j < nf; ++j) {
  83             total_power += v[i]->powers[j];
  84             den.SetDims(total_power, lower);
  85             for (int k = 0; k < v[i]->powers[j]; ++k)
  86                 den[p++] = factors[j];
  87         }
  88         for (int j = 0; j < v[i]->terms.NumRows(); ++j)
  89             gf->add(bfct->c[j], v[i]->terms[j], den);
  90         delete v[i];
  91     }
  92 }