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
 167     for (int i = 0; i < len; ++i) {
 168         ZZ coeff = num.coeff[i];
 169         vector<int> powers(dim, 0);
 170         powers[pos] = 1;
 171
 172         add_term(i, powers, coeff);
 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                 coeff = -den.coeff[j-1] * (*k)->coeff;
 180                 add_term(i, powers, coeff);
 181             }
 182         }
 183     }
 184     //dump();
 185 }
 186
 187 dpoly_r::dpoly_r(dpoly_r* num, dpoly& den, int pos, int dim)
 188 {
 189     denom = num->denom;
 190     len = num->len;
 191     c = new dpoly_r_term_list[len];
 192     this->dim = dim;
 193     ZZ coeff;
 194
 195     for (int i = 0 ; i < len; ++i) {
 196         dpoly_r_term_list::iterator k;
 197         for (k = num->c[i].begin(); k != num->c[i].end(); ++k) {
 198             vector<int> powers = (*k)->powers;
 199             powers[pos]++;
 200             add_term(i, powers, (*k)->coeff);
 201         }
 202
 203         for (int j = 1; j <= i; ++j) {
 204             dpoly_r_term_list::iterator k;
 205             for (k = c[i-j].begin(); k != c[i-j].end(); ++k) {
 206                 vector<int> powers = (*k)->powers;
 207                 powers[pos]++;
 208                 coeff = -den.coeff[j-1] * (*k)->coeff;
 209                 add_term(i, powers, coeff);
 210             }
 211         }
 212     }
 213 }
 214
 215 dpoly_r::~dpoly_r()
 216 {
 217     for (int i = 0 ; i < len; ++i)
 218         for (dpoly_r_term_list::iterator k = c[i].begin(); k != c[i].end(); ++k) {
 219             delete (*k);
 220         }
 221     delete [] c;
 222 }
 223
 224 dpoly_r *dpoly_r::div(dpoly& d)
 225 {
 226     dpoly_r *rc = new dpoly_r(len, dim);
 227     rc->denom = power(d.coeff[0], len);
 228     ZZ inv_d = rc->denom / d.coeff[0];
 229     ZZ coeff;
 230
 231     for (int i = 0; i < len; ++i) {
 232         for (dpoly_r_term_list::iterator k = c[i].begin(); k != c[i].end(); ++k) {
 233             coeff = (*k)->coeff * inv_d;
 234             rc->add_term(i, (*k)->powers, coeff);
 235         }
 236
 237         for (int j = 1; j <= i; ++j) {
 238             dpoly_r_term_list::iterator k;
 239             for (k = rc->c[i-j].begin(); k != rc->c[i-j].end(); ++k) {
 240                 coeff = - d.coeff[j] * (*k)->coeff / d.coeff[0];
 241                 rc->add_term(i, (*k)->powers, coeff);
 242             }
 243         }
 244     }
 245     return rc;
 246 }
 247
 248 void dpoly_r::dump(void)
 249 {
 250     for (int i = 0; i < len; ++i) {
 251         cerr << endl;
 252         cerr << i << endl;
 253         cerr << c[i].size() << endl;
 254         for (dpoly_r_term_list::iterator j = c[i].begin(); j != c[i].end(); ++j) {
 255             for (int k = 0; k < dim; ++k) {
 256                 cerr << (*j)->powers[k] << " ";
 257             }
 258             cerr << ": " << (*j)->coeff << "/" << denom << endl;
 259         }
 260         cerr << endl;
 261     }
 262 }