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     value_assign(denom->p[0], d.coeff->p[0]);
  71     for (int i = 1; i < len; ++i) {
  72         value_multiply(denom->p[i], denom->p[i-1], denom->p[0]);
  73         value_multiply(coeff->p[i], coeff->p[i], denom->p[i-1]);
  74
  75         mpz_submul(coeff->p[i], d.coeff->p[1], coeff->p[i-1]);
  76         for (int j = 2; j <= i; ++j) {
  77             value_multiply(tmp, denom->p[j-2], coeff->p[i-j]);
  78             mpz_submul(coeff->p[i], d.coeff->p[j], tmp);
  79         }
  80     }
  81     value_clear(tmp);
  82
  83     return denom;
  84 }
  85
  86 void dpoly::div(const dpoly& d, mpq_t count, int sign)
  87 {
  88     int len = coeff->Size;
  89     Vector *denom = div(d);
  90     mpq_t tmp;
  91     mpq_init(tmp);
  92     value_assign(mpq_numref(tmp), coeff->p[len-1]);
  93     value_assign(mpq_denref(tmp), denom->p[len-1]);
  94     mpq_canonicalize(tmp);
  95
  96     if (sign == -1)
  97         mpq_sub(count, count, tmp);
  98     else
  99         mpq_add(count, count, tmp);
 100
 101     mpq_clear(tmp);
 102     Vector_Free(denom);
 103 }
 104
 105 void dpoly::div(const dpoly& d, mpq_t *count, const mpq_t& factor)
 106 {
 107     int len = coeff->Size;
 108     Vector *denom = div(d);
 109     mpq_t tmp;
 110     mpq_init(tmp);
 111
 112     for (int i = 0; i < len; ++i) {
 113         value_multiply(mpq_numref(tmp), coeff->p[len-1 - i], mpq_numref(factor));
 114         value_multiply(mpq_denref(tmp), denom->p[len-1 - i], mpq_denref(factor));
 115         mpq_add(count[i], count[i], tmp);
 116     }
 117
 118     mpq_clear(tmp);
 119     Vector_Free(denom);
 120 }
 121
 122 void dpoly_r::add_term(int i, const vector<int>& powers, const ZZ& coeff)
 123 {
 124     if (coeff == 0)
 125         return;
 126
 127     dpoly_r_term tmp;
 128     tmp.powers = powers;
 129     dpoly_r_term_list::iterator k = c[i].find(&tmp);
 130     if (k != c[i].end()) {
 131         (*k)->coeff += coeff;
 132         return;
 133     }
 134     dpoly_r_term *t = new dpoly_r_term;
 135     t->powers = powers;
 136     t->coeff = coeff;
 137     c[i].insert(t);
 138 }
 139
 140 dpoly_r::dpoly_r(int len, int dim)
 141 {
 142     denom = 1;
 143     this->len = len;
 144     this->dim = dim;
 145     c = new dpoly_r_term_list[len];
 146 }
 147
 148 dpoly_r::dpoly_r(dpoly& num, int dim)
 149 {
 150     denom = 1;
 151     len = num.coeff->Size;
 152     c = new dpoly_r_term_list[len];
 153     this->dim = dim;
 154     vector<int> powers(dim, 0);
 155
 156     for (int i = 0; i < len; ++i) {
 157         ZZ coeff;
 158         value2zz(num.coeff->p[i], coeff);
 159         add_term(i, powers, coeff);
 160     }
 161 }
 162
 163 dpoly_r::dpoly_r(dpoly& num, dpoly& den, int pos, int dim)
 164 {
 165     denom = 1;
 166     len = num.coeff->Size;
 167     c = new dpoly_r_term_list[len];
 168     this->dim = dim;
 169     int powers[dim];
 170     ZZ coeff;
 171
 172     for (int i = 0; i < len; ++i) {
 173         vector<int> powers(dim, 0);
 174         powers[pos] = 1;
 175
 176         value2zz(num.coeff->p[i], coeff);
 177         add_term(i, powers, coeff);
 178
 179         for (int j = 1; j <= i; ++j) {
 180             dpoly_r_term_list::iterator k;
 181             for (k = c[i-j].begin(); k != c[i-j].end(); ++k) {
 182                 powers = (*k)->powers;
 183                 powers[pos]++;
 184                 value2zz(den.coeff->p[j-1], coeff);
 185                 negate(coeff, coeff);
 186                 coeff *= (*k)->coeff;
 187                 add_term(i, powers, coeff);
 188             }
 189         }
 190     }
 191     //dump();
 192 }
 193
 194 dpoly_r::dpoly_r(const dpoly_r* num, dpoly& den, int pos, int dim)
 195 {
 196     denom = num->denom;
 197     len = num->len;
 198     c = new dpoly_r_term_list[len];
 199     this->dim = dim;
 200     ZZ coeff;
 201
 202     for (int i = 0 ; i < len; ++i) {
 203         dpoly_r_term_list::iterator k;
 204         for (k = num->c[i].begin(); k != num->c[i].end(); ++k) {
 205             vector<int> powers = (*k)->powers;
 206             powers[pos]++;
 207             add_term(i, powers, (*k)->coeff);
 208         }
 209
 210         for (int j = 1; j <= i; ++j) {
 211             dpoly_r_term_list::iterator k;
 212             for (k = c[i-j].begin(); k != c[i-j].end(); ++k) {
 213                 vector<int> powers = (*k)->powers;
 214                 powers[pos]++;
 215                 value2zz(den.coeff->p[j-1], coeff);
 216                 negate(coeff, coeff);
 217                 coeff *= (*k)->coeff;
 218                 add_term(i, powers, coeff);
 219             }
 220         }
 221     }
 222 }
 223
 224 dpoly_r::~dpoly_r()
 225 {
 226     for (int i = 0 ; i < len; ++i)
 227         for (dpoly_r_term_list::iterator k = c[i].begin(); k != c[i].end(); ++k) {
 228             delete (*k);
 229         }
 230     delete [] c;
 231 }
 232
 233 dpoly_r *dpoly_r::div(const dpoly& d) const
 234 {
 235     dpoly_r *rc = new dpoly_r(len, dim);
 236     ZZ coeff;
 237     ZZ coeff0;
 238     value2zz(d.coeff->p[0], coeff0);
 239     rc->denom = power(coeff0, len);
 240     ZZ inv_d = rc->denom / coeff0;
 241
 242     for (int i = 0; i < len; ++i) {
 243         for (dpoly_r_term_list::iterator k = c[i].begin(); k != c[i].end(); ++k) {
 244             coeff = (*k)->coeff * inv_d;
 245             rc->add_term(i, (*k)->powers, coeff);
 246         }
 247
 248         for (int j = 1; j <= i; ++j) {
 249             dpoly_r_term_list::iterator k;
 250             for (k = rc->c[i-j].begin(); k != rc->c[i-j].end(); ++k) {
 251                 value2zz(d.coeff->p[j], coeff);
 252                 coeff = - coeff * (*k)->coeff / coeff0;
 253                 rc->add_term(i, (*k)->powers, coeff);
 254             }
 255         }
 256     }
 257     return rc;
 258 }
 259
 260 void dpoly_r::dump(void)
 261 {
 262     for (int i = 0; i < len; ++i) {
 263         cerr << endl;
 264         cerr << i << endl;
 265         cerr << c[i].size() << endl;
 266         for (dpoly_r_term_list::iterator j = c[i].begin(); j != c[i].end(); ++j) {
 267             for (int k = 0; k < dim; ++k) {
 268                 cerr << (*j)->powers[k] << " ";
 269             }
 270             cerr << ": " << (*j)->coeff << "/" << denom << endl;
 271         }
 272         cerr << endl;
 273     }
 274 }