9 /* Construct truncated expansion of (1+t)^(degree),
10 * computing the first 1+d coefficients
12 dpoly::dpoly(int d
, ZZ
& degree
, int offset
)
16 /* For small degrees, we only need to compute some coefficients */
18 if (degree
>= 0 && degree
< ZZ(INIT_VAL
, min
))
21 ZZ c
= ZZ(INIT_VAL
, 1);
24 for (int i
= 1; i
<= min
; ++i
) {
31 void dpoly::operator += (const dpoly
& t
)
33 assert(coeff
.length() == t
.coeff
.length());
34 for (int i
= 0; i
< coeff
.length(); ++i
)
35 coeff
[i
] += t
.coeff
[i
];
38 void dpoly::operator *= (const ZZ
& f
)
40 for (int i
= 0; i
< coeff
.length(); ++i
)
44 void dpoly::operator *= (dpoly
& f
)
46 assert(coeff
.length() == f
.coeff
.length());
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
];
54 mpq_t
*dpoly::div(dpoly
& d
) const
56 int len
= coeff
.length();
59 mpq_t
* c
= new mpq_t
[coeff
.length()];
62 for (int i
= 0; i
< len
; ++i
) {
64 zz2value(coeff
[i
], tmp
);
67 for (int j
= 1; j
<= i
; ++j
) {
68 zz2value(d
.coeff
[j
], tmp
);
70 mpq_mul(qtmp
, qtmp
, c
[i
-j
]);
71 mpq_sub(c
[i
], c
[i
], qtmp
);
74 zz2value(d
.coeff
[0], tmp
);
76 mpq_div(c
[i
], c
[i
], qtmp
);
84 void dpoly::clear_div(mpq_t
*c
) const
86 int len
= coeff
.length();
88 for (int i
= 0; i
< len
; ++i
)
93 void dpoly::div(dpoly
& d
, mpq_t count
, ZZ
& sign
)
95 int len
= coeff
.length();
99 mpq_sub(count
, count
, c
[len
-1]);
101 mpq_add(count
, count
, c
[len
-1]);
106 void dpoly::div(dpoly
& d
, mpq_t
*count
, const mpq_t
& factor
)
108 int len
= coeff
.length();
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
]);
119 void dpoly_r::add_term(int i
, const vector
<int>& powers
, const ZZ
& coeff
)
126 dpoly_r_term_list::iterator k
= c
[i
].find(&tmp
);
127 if (k
!= c
[i
].end()) {
128 (*k
)->coeff
+= coeff
;
131 dpoly_r_term
*t
= new dpoly_r_term
;
137 dpoly_r::dpoly_r(int len
, int dim
)
142 c
= new dpoly_r_term_list
[len
];
145 dpoly_r::dpoly_r(dpoly
& num
, int dim
)
148 len
= num
.coeff
.length();
149 c
= new dpoly_r_term_list
[len
];
151 vector
<int> powers(dim
, 0);
153 for (int i
= 0; i
< len
; ++i
) {
154 ZZ coeff
= num
.coeff
[i
];
155 add_term(i
, powers
, coeff
);
159 dpoly_r::dpoly_r(dpoly
& num
, dpoly
& den
, int pos
, int dim
)
162 len
= num
.coeff
.length();
163 c
= new dpoly_r_term_list
[len
];
167 for (int i
= 0; i
< len
; ++i
) {
168 ZZ coeff
= num
.coeff
[i
];
169 vector
<int> powers(dim
, 0);
172 add_term(i
, powers
, coeff
);
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
;
179 coeff
= -den
.coeff
[j
-1] * (*k
)->coeff
;
180 add_term(i
, powers
, coeff
);
187 dpoly_r::dpoly_r(const dpoly_r
* num
, dpoly
& den
, int pos
, int dim
)
191 c
= new dpoly_r_term_list
[len
];
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
;
200 add_term(i
, powers
, (*k
)->coeff
);
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
;
208 coeff
= -den
.coeff
[j
-1] * (*k
)->coeff
;
209 add_term(i
, powers
, coeff
);
217 for (int i
= 0 ; i
< len
; ++i
)
218 for (dpoly_r_term_list::iterator k
= c
[i
].begin(); k
!= c
[i
].end(); ++k
) {
224 dpoly_r
*dpoly_r::div(const dpoly
& d
) const
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];
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
);
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
);
248 void dpoly_r::dump(void)
250 for (int i
= 0; i
< len
; ++i
) {
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
] << " ";
258 cerr
<< ": " << (*j
)->coeff
<< "/" << denom
<< endl
;