dpoly.cc

   1 #include <assert.h>
   2 #include <iostream>
   3 #include <vector>
   4 #include <NTL/ZZ.h>
   5 #include <barvinok/polylib.h>
   6 #include "dpoly.h"
   7
   8 using std::cerr;
   9 using std::endl;
  10 using std::vector;
  11
  12 /* Construct truncated expansion of (1+t)^(degree),
  13  * computing the first 1+d coefficients
  14  */
  15 dpoly::dpoly(int d, const Value degree, int offset)
  16 {
  17     coeff = Vector_Alloc(d+1);
  18
  19     /* For small degrees, we only need to compute some coefficients */
  20     int min = d + offset;
  21     if (value_posz_p(degree) && value_cmp_si(degree, min) < 0)
  22         min = VALUE_TO_INT(degree);
  23
  24     Value c, tmp;
  25     value_init(c);
  26     value_init(tmp);
  27     value_set_si(c, 1);
  28     if (!offset)
  29         value_assign(coeff->p[0], c);
  30     value_assign(tmp, degree);
  31     for (int i = 1; i <= min; ++i) {
  32         value_multiply(c, c, tmp);
  33         value_decrement(tmp, tmp);
  34         mpz_divexact_ui(c, c, i);
  35         value_assign(coeff->p[i-offset], c);
  36     }
  37     value_clear(c);
  38     value_clear(tmp);
  39 }
  40
  41 void dpoly::operator += (const dpoly& t)
  42 {
  43     assert(coeff->Size == t.coeff->Size);
  44     for (int i = 0; i < coeff->Size; ++i)
  45         value_addto(coeff->p[i], coeff->p[i], t.coeff->p[i]);
  46 }
  47
  48 void dpoly::operator *= (const Value f)
  49 {
  50     for (int i = 0; i < coeff->Size; ++i)
  51         value_multiply(coeff->p[i], coeff->p[i], f);
  52 }
  53
  54 void dpoly::operator *= (const dpoly& f)
  55 {
  56     assert(coeff->Size == f.coeff->Size);
  57     Vector *old = Vector_Alloc(coeff->Size);
  58     Vector_Copy(coeff->p, old->p, coeff->Size);
  59     Vector_Scale(coeff->p, coeff->p, f.coeff->p[0], coeff->Size);
  60     for (int i = 1; i < coeff->Size; ++i)
  61         for (int j = 0; i+j < coeff->Size; ++j)
  62             value_addmul(coeff->p[i+j], f.coeff->p[i], old->p[j]);
  63     Vector_Free(old);
  64 }
  65
  66 Vector *dpoly::div(const dpoly& d)
  67 {
  68     int len = coeff->Size;
  69     Vector *denom = Vector_Alloc(len);
  70     Value tmp;
  71     value_init(tmp);
  72
  73     /* Make sure denominators are positive */
  74     if (value_neg_p(d.coeff->p[0])) {
  75         Vector_Oppose(d.coeff->p, d.coeff->p, d.coeff->Size);
  76         Vector_Oppose(coeff->p, coeff->p, coeff->Size);
  77     }
  78     value_assign(denom->p[0], d.coeff->p[0]);
  79     for (int i = 1; i < len; ++i) {
  80         value_multiply(denom->p[i], denom->p[i-1], denom->p[0]);
  81         value_multiply(coeff->p[i], coeff->p[i], denom->p[i-1]);
  82
  83         mpz_submul(coeff->p[i], d.coeff->p[1], coeff->p[i-1]);
  84         for (int j = 2; j <= i; ++j) {
  85             value_multiply(tmp, denom->p[j-2], coeff->p[i-j]);
  86             mpz_submul(coeff->p[i], d.coeff->p[j], tmp);
  87         }
  88     }
  89     value_clear(tmp);
  90
  91     return denom;
  92 }
  93
  94 void dpoly::div(const dpoly& d, mpq_t count, int sign)
  95 {
  96     int len = coeff->Size;
  97     Vector *denom = div(d);
  98     mpq_t tmp;
  99     mpq_init(tmp);
 100     value_assign(mpq_numref(tmp), coeff->p[len-1]);
 101     value_assign(mpq_denref(tmp), denom->p[len-1]);
 102     mpq_canonicalize(tmp);
 103
 104     if (sign == -1)
 105         mpq_sub(count, count, tmp);
 106     else
 107         mpq_add(count, count, tmp);
 108
 109     mpq_clear(tmp);
 110     Vector_Free(denom);
 111 }
 112
 113 void dpoly::div(const dpoly& d, mpq_t *count, const mpq_t& factor)
 114 {
 115     int len = coeff->Size;
 116     Vector *denom = div(d);
 117     mpq_t tmp;
 118     mpq_init(tmp);
 119
 120     for (int i = 0; i < len; ++i) {
 121         value_multiply(mpq_numref(tmp), coeff->p[len-1 - i], mpq_numref(factor));
 122         value_multiply(mpq_denref(tmp), denom->p[len-1 - i], mpq_denref(factor));
 123         mpq_add(count[i], count[i], tmp);
 124         mpq_canonicalize(count[i]);
 125     }
 126
 127     mpq_clear(tmp);
 128     Vector_Free(denom);
 129 }
 130
 131 void dpoly_r::add_term(int i, const vector<int>& powers, const ZZ& coeff)
 132 {
 133     if (coeff == 0)
 134         return;
 135
 136     dpoly_r_term tmp;
 137     tmp.powers = powers;
 138     dpoly_r_term_list::iterator k = c[i].find(&tmp);
 139     if (k != c[i].end()) {
 140         (*k)->coeff += coeff;
 141         return;
 142     }
 143     dpoly_r_term *t = new dpoly_r_term;
 144     t->powers = powers;
 145     t->coeff = coeff;
 146     c[i].insert(t);
 147 }
 148
 149 dpoly_r::dpoly_r(int len, int dim)
 150 {
 151     denom = 1;
 152     this->len = len;
 153     this->dim = dim;
 154     c = new dpoly_r_term_list[len];
 155 }
 156
 157 dpoly_r::dpoly_r(dpoly& num, int dim)
 158 {
 159     denom = 1;
 160     len = num.coeff->Size;
 161     c = new dpoly_r_term_list[len];
 162     this->dim = dim;
 163     vector<int> powers(dim, 0);
 164
 165     for (int i = 0; i < len; ++i) {
 166         ZZ coeff;
 167         value2zz(num.coeff->p[i], coeff);
 168         add_term(i, powers, coeff);
 169     }
 170 }
 171
 172 dpoly_r::dpoly_r(dpoly& num, dpoly& den, int pos, int dim)
 173 {
 174     denom = 1;
 175     len = num.coeff->Size;
 176     c = new dpoly_r_term_list[len];
 177     this->dim = dim;
 178     int *powers = new int[dim];
 179     ZZ coeff;
 180
 181     for (int i = 0; i < len; ++i) {
 182         vector<int> powers(dim, 0);
 183         powers[pos] = 1;
 184
 185         value2zz(num.coeff->p[i], coeff);
 186         add_term(i, powers, coeff);
 187
 188         for (int j = 1; j <= i; ++j) {
 189             dpoly_r_term_list::iterator k;
 190             for (k = c[i-j].begin(); k != c[i-j].end(); ++k) {
 191                 powers = (*k)->powers;
 192                 powers[pos]++;
 193                 value2zz(den.coeff->p[j-1], coeff);
 194                 negate(coeff, coeff);
 195                 coeff *= (*k)->coeff;
 196                 add_term(i, powers, coeff);
 197             }
 198         }
 199     }
 200     delete [] powers;
 201     //dump();
 202 }
 203
 204 dpoly_r::dpoly_r(const dpoly_r* num, dpoly& den, int pos, int dim)
 205 {
 206     denom = num->denom;
 207     len = num->len;
 208     c = new dpoly_r_term_list[len];
 209     this->dim = dim;
 210     ZZ coeff;
 211
 212     for (int i = 0 ; i < len; ++i) {
 213         dpoly_r_term_list::iterator k;
 214         for (k = num->c[i].begin(); k != num->c[i].end(); ++k) {
 215             vector<int> powers = (*k)->powers;
 216             powers[pos]++;
 217             add_term(i, powers, (*k)->coeff);
 218         }
 219
 220         for (int j = 1; j <= i; ++j) {
 221             dpoly_r_term_list::iterator k;
 222             for (k = c[i-j].begin(); k != c[i-j].end(); ++k) {
 223                 vector<int> powers = (*k)->powers;
 224                 powers[pos]++;
 225                 value2zz(den.coeff->p[j-1], coeff);
 226                 negate(coeff, coeff);
 227                 coeff *= (*k)->coeff;
 228                 add_term(i, powers, coeff);
 229             }
 230         }
 231     }
 232 }
 233
 234 dpoly_r::~dpoly_r()
 235 {
 236     for (int i = 0 ; i < len; ++i)
 237         for (dpoly_r_term_list::iterator k = c[i].begin(); k != c[i].end(); ++k) {
 238             delete (*k);
 239         }
 240     delete [] c;
 241 }
 242
 243 dpoly_r *dpoly_r::div(const dpoly& d) const
 244 {
 245     dpoly_r *rc = new dpoly_r(len, dim);
 246     ZZ coeff;
 247     ZZ coeff0;
 248     value2zz(d.coeff->p[0], coeff0);
 249     rc->denom = power(coeff0, len);
 250     ZZ inv_d = rc->denom / coeff0;
 251
 252     for (int i = 0; i < len; ++i) {
 253         for (dpoly_r_term_list::iterator k = c[i].begin(); k != c[i].end(); ++k) {
 254             coeff = (*k)->coeff * inv_d;
 255             rc->add_term(i, (*k)->powers, coeff);
 256         }
 257
 258         for (int j = 1; j <= i; ++j) {
 259             dpoly_r_term_list::iterator k;
 260             for (k = rc->c[i-j].begin(); k != rc->c[i-j].end(); ++k) {
 261                 value2zz(d.coeff->p[j], coeff);
 262                 coeff = - coeff * (*k)->coeff / coeff0;
 263                 rc->add_term(i, (*k)->powers, coeff);
 264             }
 265         }
 266     }
 267     return rc;
 268 }
 269
 270 void dpoly_r::dump(void)
 271 {
 272     for (int i = 0; i < len; ++i) {
 273         cerr << endl;
 274         cerr << i << endl;
 275         cerr << c[i].size() << endl;
 276         for (dpoly_r_term_list::iterator j = c[i].begin(); j != c[i].end(); ++j) {
 277             for (int k = 0; k < dim; ++k) {
 278                 cerr << (*j)->powers[k] << " ";
 279             }
 280             cerr << ": " << (*j)->coeff << "/" << denom << endl;
 281         }
 282         cerr << endl;
 283     }
 284 }