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