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 */
19 if (value_posz_p(degree
) && value_cmp_si(degree
, min
) < 0)
20 min
= VALUE_TO_INT(degree
);
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
);
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
]);
64 Vector
*dpoly::div(const dpoly
& d
)
66 int len
= coeff
->Size
;
67 Vector
*denom
= Vector_Alloc(len
);
70 value_assign(denom
->p
[0], d
.coeff
->p
[0]);
71 for (int i
= 1; i
< len
; ++i
) {
72 value_multiply(denom
->p
[i
], denom
->p
[i
-1], denom
->p
[0]);
73 value_multiply(coeff
->p
[i
], coeff
->p
[i
], denom
->p
[i
-1]);
75 mpz_submul(coeff
->p
[i
], d
.coeff
->p
[1], coeff
->p
[i
-1]);
76 for (int j
= 2; j
<= i
; ++j
) {
77 value_multiply(tmp
, denom
->p
[j
-2], coeff
->p
[i
-j
]);
78 mpz_submul(coeff
->p
[i
], d
.coeff
->p
[j
], tmp
);
86 void dpoly::div(const dpoly
& d
, mpq_t count
, int sign
)
88 int len
= coeff
->Size
;
89 Vector
*denom
= div(d
);
92 value_assign(mpq_numref(tmp
), coeff
->p
[len
-1]);
93 value_assign(mpq_denref(tmp
), denom
->p
[len
-1]);
94 mpq_canonicalize(tmp
);
97 mpq_sub(count
, count
, tmp
);
99 mpq_add(count
, count
, tmp
);
105 void dpoly::div(const dpoly
& d
, mpq_t
*count
, const mpq_t
& factor
)
107 int len
= coeff
->Size
;
108 Vector
*denom
= div(d
);
112 for (int i
= 0; i
< len
; ++i
) {
113 value_multiply(mpq_numref(tmp
), coeff
->p
[len
-1 - i
], mpq_numref(factor
));
114 value_multiply(mpq_denref(tmp
), denom
->p
[len
-1 - i
], mpq_denref(factor
));
115 mpq_add(count
[i
], count
[i
], tmp
);
122 void dpoly_r::add_term(int i
, const vector
<int>& powers
, const ZZ
& coeff
)
129 dpoly_r_term_list::iterator k
= c
[i
].find(&tmp
);
130 if (k
!= c
[i
].end()) {
131 (*k
)->coeff
+= coeff
;
134 dpoly_r_term
*t
= new dpoly_r_term
;
140 dpoly_r::dpoly_r(int len
, int dim
)
145 c
= new dpoly_r_term_list
[len
];
148 dpoly_r::dpoly_r(dpoly
& num
, int dim
)
151 len
= num
.coeff
->Size
;
152 c
= new dpoly_r_term_list
[len
];
154 vector
<int> powers(dim
, 0);
156 for (int i
= 0; i
< len
; ++i
) {
158 value2zz(num
.coeff
->p
[i
], coeff
);
159 add_term(i
, powers
, coeff
);
163 dpoly_r::dpoly_r(dpoly
& num
, dpoly
& den
, int pos
, int dim
)
166 len
= num
.coeff
->Size
;
167 c
= new dpoly_r_term_list
[len
];
172 for (int i
= 0; i
< len
; ++i
) {
173 vector
<int> powers(dim
, 0);
176 value2zz(num
.coeff
->p
[i
], coeff
);
177 add_term(i
, powers
, coeff
);
179 for (int j
= 1; j
<= i
; ++j
) {
180 dpoly_r_term_list::iterator k
;
181 for (k
= c
[i
-j
].begin(); k
!= c
[i
-j
].end(); ++k
) {
182 powers
= (*k
)->powers
;
184 value2zz(den
.coeff
->p
[j
-1], coeff
);
185 negate(coeff
, coeff
);
186 coeff
*= (*k
)->coeff
;
187 add_term(i
, powers
, coeff
);
194 dpoly_r::dpoly_r(const dpoly_r
* num
, dpoly
& den
, int pos
, int dim
)
198 c
= new dpoly_r_term_list
[len
];
202 for (int i
= 0 ; i
< len
; ++i
) {
203 dpoly_r_term_list::iterator k
;
204 for (k
= num
->c
[i
].begin(); k
!= num
->c
[i
].end(); ++k
) {
205 vector
<int> powers
= (*k
)->powers
;
207 add_term(i
, powers
, (*k
)->coeff
);
210 for (int j
= 1; j
<= i
; ++j
) {
211 dpoly_r_term_list::iterator k
;
212 for (k
= c
[i
-j
].begin(); k
!= c
[i
-j
].end(); ++k
) {
213 vector
<int> powers
= (*k
)->powers
;
215 value2zz(den
.coeff
->p
[j
-1], coeff
);
216 negate(coeff
, coeff
);
217 coeff
*= (*k
)->coeff
;
218 add_term(i
, powers
, coeff
);
226 for (int i
= 0 ; i
< len
; ++i
)
227 for (dpoly_r_term_list::iterator k
= c
[i
].begin(); k
!= c
[i
].end(); ++k
) {
233 dpoly_r
*dpoly_r::div(const dpoly
& d
) const
235 dpoly_r
*rc
= new dpoly_r(len
, dim
);
238 value2zz(d
.coeff
->p
[0], coeff0
);
239 rc
->denom
= power(coeff0
, len
);
240 ZZ inv_d
= rc
->denom
/ coeff0
;
242 for (int i
= 0; i
< len
; ++i
) {
243 for (dpoly_r_term_list::iterator k
= c
[i
].begin(); k
!= c
[i
].end(); ++k
) {
244 coeff
= (*k
)->coeff
* inv_d
;
245 rc
->add_term(i
, (*k
)->powers
, coeff
);
248 for (int j
= 1; j
<= i
; ++j
) {
249 dpoly_r_term_list::iterator k
;
250 for (k
= rc
->c
[i
-j
].begin(); k
!= rc
->c
[i
-j
].end(); ++k
) {
251 value2zz(d
.coeff
->p
[j
], coeff
);
252 coeff
= - coeff
* (*k
)->coeff
/ coeff0
;
253 rc
->add_term(i
, (*k
)->powers
, coeff
);
260 void dpoly_r::dump(void)
262 for (int i
= 0; i
< len
; ++i
) {
265 cerr
<< c
[i
].size() << endl
;
266 for (dpoly_r_term_list::iterator j
= c
[i
].begin(); j
!= c
[i
].end(); ++j
) {
267 for (int k
= 0; k
< dim
; ++k
) {
268 cerr
<< (*j
)->powers
[k
] << " ";
270 cerr
<< ": " << (*j
)->coeff
<< "/" << denom
<< endl
;