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         mpq_canonicalize(count[i]);
 123     }
 124
 125     mpq_clear(tmp);
 126     Vector_Free(denom);
 127 }
 128
 129 void dpoly_r::add_term(int i, const vector<int>& powers, const ZZ& coeff)
 130 {
 131     if (coeff == 0)
 132         return;
 133
 134     dpoly_r_term tmp;
 135     tmp.powers = powers;
 136     dpoly_r_term_list::iterator k = c[i].find(&tmp);
 137     if (k != c[i].end()) {
 138         (*k)->coeff += coeff;
 139         return;
 140     }
 141     dpoly_r_term *t = new dpoly_r_term;
 142     t->powers = powers;
 143     t->coeff = coeff;
 144     c[i].insert(t);
 145 }
 146
 147 dpoly_r::dpoly_r(int len, int dim)
 148 {
 149     denom = 1;
 150     this->len = len;
 151     this->dim = dim;
 152     c = new dpoly_r_term_list[len];
 153 }
 154
 155 dpoly_r::dpoly_r(dpoly& num, int dim)
 156 {
 157     denom = 1;
 158     len = num.coeff->Size;
 159     c = new dpoly_r_term_list[len];
 160     this->dim = dim;
 161     vector<int> powers(dim, 0);
 162
 163     for (int i = 0; i < len; ++i) {
 164         ZZ coeff;
 165         value2zz(num.coeff->p[i], coeff);
 166         add_term(i, powers, coeff);
 167     }
 168 }
 169
 170 dpoly_r::dpoly_r(dpoly& num, dpoly& den, int pos, int dim)
 171 {
 172     denom = 1;
 173     len = num.coeff->Size;
 174     c = new dpoly_r_term_list[len];
 175     this->dim = dim;
 176     int powers[dim];
 177     ZZ coeff;
 178
 179     for (int i = 0; i < len; ++i) {
 180         vector<int> powers(dim, 0);
 181         powers[pos] = 1;
 182
 183         value2zz(num.coeff->p[i], coeff);
 184         add_term(i, powers, coeff);
 185
 186         for (int j = 1; j <= i; ++j) {
 187             dpoly_r_term_list::iterator k;
 188             for (k = c[i-j].begin(); k != c[i-j].end(); ++k) {
 189                 powers = (*k)->powers;
 190                 powers[pos]++;
 191                 value2zz(den.coeff->p[j-1], coeff);
 192                 negate(coeff, coeff);
 193                 coeff *= (*k)->coeff;
 194                 add_term(i, powers, coeff);
 195             }
 196         }
 197     }
 198     //dump();
 199 }
 200
 201 dpoly_r::dpoly_r(const dpoly_r* num, dpoly& den, int pos, int dim)
 202 {
 203     denom = num->denom;
 204     len = num->len;
 205     c = new dpoly_r_term_list[len];
 206     this->dim = dim;
 207     ZZ coeff;
 208
 209     for (int i = 0 ; i < len; ++i) {
 210         dpoly_r_term_list::iterator k;
 211         for (k = num->c[i].begin(); k != num->c[i].end(); ++k) {
 212             vector<int> powers = (*k)->powers;
 213             powers[pos]++;
 214             add_term(i, powers, (*k)->coeff);
 215         }
 216
 217         for (int j = 1; j <= i; ++j) {
 218             dpoly_r_term_list::iterator k;
 219             for (k = c[i-j].begin(); k != c[i-j].end(); ++k) {
 220                 vector<int> powers = (*k)->powers;
 221                 powers[pos]++;
 222                 value2zz(den.coeff->p[j-1], coeff);
 223                 negate(coeff, coeff);
 224                 coeff *= (*k)->coeff;
 225                 add_term(i, powers, coeff);
 226             }
 227         }
 228     }
 229 }
 230
 231 dpoly_r::~dpoly_r()
 232 {
 233     for (int i = 0 ; i < len; ++i)
 234         for (dpoly_r_term_list::iterator k = c[i].begin(); k != c[i].end(); ++k) {
 235             delete (*k);
 236         }
 237     delete [] c;
 238 }
 239
 240 dpoly_r *dpoly_r::div(const dpoly& d) const
 241 {
 242     dpoly_r *rc = new dpoly_r(len, dim);
 243     ZZ coeff;
 244     ZZ coeff0;
 245     value2zz(d.coeff->p[0], coeff0);
 246     rc->denom = power(coeff0, len);
 247     ZZ inv_d = rc->denom / coeff0;
 248
 249     for (int i = 0; i < len; ++i) {
 250         for (dpoly_r_term_list::iterator k = c[i].begin(); k != c[i].end(); ++k) {
 251             coeff = (*k)->coeff * inv_d;
 252             rc->add_term(i, (*k)->powers, coeff);
 253         }
 254
 255         for (int j = 1; j <= i; ++j) {
 256             dpoly_r_term_list::iterator k;
 257             for (k = rc->c[i-j].begin(); k != rc->c[i-j].end(); ++k) {
 258                 value2zz(d.coeff->p[j], coeff);
 259                 coeff = - coeff * (*k)->coeff / coeff0;
 260                 rc->add_term(i, (*k)->powers, coeff);
 261             }
 262         }
 263     }
 264     return rc;
 265 }
 266
 267 void dpoly_r::dump(void)
 268 {
 269     for (int i = 0; i < len; ++i) {
 270         cerr << endl;
 271         cerr << i << endl;
 272         cerr << c[i].size() << endl;
 273         for (dpoly_r_term_list::iterator j = c[i].begin(); j != c[i].end(); ++j) {
 274             for (int k = 0; k < dim; ++k) {
 275                 cerr << (*j)->powers[k] << " ";
 276             }
 277             cerr << ": " << (*j)->coeff << "/" << denom << endl;
 278         }
 279         cerr << endl;
 280     }
 281 }