dpoly.cc

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