dpoly.cc

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