add some iscc tests
[barvinok.git] / dpoly.cc
bloba3c90b1114b9d7ee36c4d7d846d16ed0a3e85fc0
1 #include <assert.h>
2 #include <iostream>
3 #include <vector>
4 #include "dpoly.h"
6 using std::cerr;
7 using std::endl;
8 using std::vector;
10 /* Construct truncated expansion of (1+t)^(degree),
11 * computing the first 1+d coefficients
13 dpoly::dpoly(int d, const Value degree, int offset)
15 coeff = Vector_Alloc(d+1);
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);
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);
35 value_clear(c);
36 value_clear(tmp);
39 void dpoly::operator += (const dpoly& t)
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]);
46 void dpoly::operator *= (const Value f)
48 for (int i = 0; i < coeff->Size; ++i)
49 value_multiply(coeff->p[i], coeff->p[i], f);
52 void dpoly::operator *= (const dpoly& f)
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);
64 Vector *dpoly::div(const dpoly& d)
66 int len = coeff->Size;
67 Vector *denom = Vector_Alloc(len);
68 Value tmp;
69 value_init(tmp);
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);
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]);
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);
87 value_clear(tmp);
89 return denom;
92 void dpoly::div(const dpoly& d, mpq_t count, int sign)
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);
102 if (sign == -1)
103 mpq_sub(count, count, tmp);
104 else
105 mpq_add(count, count, tmp);
107 mpq_clear(tmp);
108 Vector_Free(denom);
111 void dpoly::div(const dpoly& d, mpq_t *count, const mpq_t& factor)
113 int len = coeff->Size;
114 Vector *denom = div(d);
115 mpq_t tmp;
116 mpq_init(tmp);
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]);
125 mpq_clear(tmp);
126 Vector_Free(denom);
129 void dpoly_r::add_term(int i, const vector<int>& powers, const ZZ& coeff)
131 if (coeff == 0)
132 return;
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;
141 dpoly_r_term *t = new dpoly_r_term;
142 t->powers = powers;
143 t->coeff = coeff;
144 c[i].insert(t);
147 dpoly_r::dpoly_r(int len, int dim)
149 denom = 1;
150 this->len = len;
151 this->dim = dim;
152 c = new dpoly_r_term_list[len];
155 dpoly_r::dpoly_r(dpoly& num, int dim)
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);
163 for (int i = 0; i < len; ++i) {
164 ZZ coeff;
165 value2zz(num.coeff->p[i], coeff);
166 add_term(i, powers, coeff);
170 dpoly_r::dpoly_r(dpoly& num, dpoly& den, int pos, int dim)
172 denom = 1;
173 len = num.coeff->Size;
174 c = new dpoly_r_term_list[len];
175 this->dim = dim;
176 int *powers = new int[dim];
177 ZZ coeff;
179 for (int i = 0; i < len; ++i) {
180 vector<int> powers(dim, 0);
181 powers[pos] = 1;
183 value2zz(num.coeff->p[i], coeff);
184 add_term(i, powers, coeff);
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);
198 delete [] powers;
199 //dump();
202 dpoly_r::dpoly_r(const dpoly_r* num, dpoly& den, int pos, int dim)
204 denom = num->denom;
205 len = num->len;
206 c = new dpoly_r_term_list[len];
207 this->dim = dim;
208 ZZ coeff;
210 for (int i = 0 ; i < len; ++i) {
211 dpoly_r_term_list::iterator k;
212 for (k = num->c[i].begin(); k != num->c[i].end(); ++k) {
213 vector<int> powers = (*k)->powers;
214 powers[pos]++;
215 add_term(i, powers, (*k)->coeff);
218 for (int j = 1; j <= i; ++j) {
219 dpoly_r_term_list::iterator k;
220 for (k = c[i-j].begin(); k != c[i-j].end(); ++k) {
221 vector<int> powers = (*k)->powers;
222 powers[pos]++;
223 value2zz(den.coeff->p[j-1], coeff);
224 negate(coeff, coeff);
225 coeff *= (*k)->coeff;
226 add_term(i, powers, coeff);
232 dpoly_r::~dpoly_r()
234 for (int i = 0 ; i < len; ++i)
235 for (dpoly_r_term_list::iterator k = c[i].begin(); k != c[i].end(); ++k) {
236 delete (*k);
238 delete [] c;
241 dpoly_r *dpoly_r::div(const dpoly& d) const
243 dpoly_r *rc = new dpoly_r(len, dim);
244 ZZ coeff;
245 ZZ coeff0;
246 value2zz(d.coeff->p[0], coeff0);
247 rc->denom = power(coeff0, len);
248 ZZ inv_d = rc->denom / coeff0;
250 for (int i = 0; i < len; ++i) {
251 for (dpoly_r_term_list::iterator k = c[i].begin(); k != c[i].end(); ++k) {
252 coeff = (*k)->coeff * inv_d;
253 rc->add_term(i, (*k)->powers, coeff);
256 for (int j = 1; j <= i; ++j) {
257 dpoly_r_term_list::iterator k;
258 for (k = rc->c[i-j].begin(); k != rc->c[i-j].end(); ++k) {
259 value2zz(d.coeff->p[j], coeff);
260 coeff = - coeff * (*k)->coeff / coeff0;
261 rc->add_term(i, (*k)->powers, coeff);
265 return rc;
268 void dpoly_r::dump(void)
270 for (int i = 0; i < len; ++i) {
271 cerr << endl;
272 cerr << i << endl;
273 cerr << c[i].size() << endl;
274 for (dpoly_r_term_list::iterator j = c[i].begin(); j != c[i].end(); ++j) {
275 for (int k = 0; k < dim; ++k) {
276 cerr << (*j)->powers[k] << " ";
278 cerr << ": " << (*j)->coeff << "/" << denom << endl;
280 cerr << endl;