- exvector yun = sqrfree_yun(denom, x);
-//clog << "yun factors: " << exprseq(yun) << endl;
- size_t num_yun = yun.size();
- exvector factor; factor.reserve(num_yun);
- exvector cofac; cofac.reserve(num_yun);
- for (size_t i=0; i<num_yun; i++) {
- if (!yun[i].is_equal(_ex1)) {
- for (size_t j=0; j<=i; j++) {
- factor.push_back(pow(yun[i], j+1));
- ex prod = _ex1;
- for (size_t k=0; k<num_yun; k++) {
- if (k == i)
- prod *= pow(yun[k], i-j);
- else
- prod *= pow(yun[k], k+1);
- }
- cofac.push_back(prod.expand());
+ epvector yun = sqrfree_yun(denom, x);
+ size_t yun_max_exponent = yun.empty() ? 0 : ex_to<numeric>(yun.back().coeff).to_int();
+ exvector factor, cofac;
+ for (size_t i=0; i<yun.size(); i++) {
+ numeric i_exponent = ex_to<numeric>(yun[i].coeff);
+ for (size_t j=0; j<i_exponent; j++) {
+ factor.push_back(pow(yun[i].rest, j+1));
+ ex prod = _ex1;
+ for (size_t k=0; k<yun.size(); k++) {
+ if (yun[k].coeff == i_exponent)
+ prod *= pow(yun[k].rest, i_exponent-1-j);
+ else
+ prod *= pow(yun[k].rest, yun[k].coeff);