-/*
- * Simplify the expression x by applying the rules
- *
- * 1/(z+n)/(z+m) = 1/(n-m)/(z+m) - 1/(n-m)/(z+n);
- * z^m/(z+n) = z^(m-1) - n*z^(m-1)/(z+n)
- *
- * as often as possible, where z is given as an argument to the function
- * and n and m are arbitrary positive numbers.
- */
-ex poles_simplify(const ex &x, const ex &z)
-{ if (is_a<mul>(x))
- { for (int i=0; i<x.nops(); ++i)
- { ex arg1;
- if (!is_z_pole(x.op(i), z, arg1))
- continue;
- for (int j=i+1; j<x.nops(); ++j)
- { ex arg2;
- if (!is_z_pole(x.op(j), z, arg2))
- continue;
- ex result = x/x.op(i)/(arg1-arg2)-x/x.op(j)/(arg1-arg2);
- result = poles_simplify(result, z);
- return result;
- }
- }
- ex result = expand(x);
- if (is_a<add>(result))
- return poles_simplify(result, z);
- for (int i=0; i<x.nops(); ++i)
- { ex arg1;
- if (!is_z_pole(x.op(i), z, arg1))
- continue;
- for (int j=0; j<x.nops(); ++j)
- { ex expon = 0;
- if (is_a<power>(x.op(j)) && x.op(j).op(0)==z
- && x.op(j).op(1).info(info_flags::posint))
- expon = x.op(j).op(1);
- if (x.op(j) == z)
- expon = 1;
- if (expon == 0)
- continue;
- ex result = x/x.op(i)/z - arg1*x/z;
- result = poles_simplify(result, z);
- return result;
- }
+ for (size_t i=2; i<n; ++i)
+ {
+ // next = previous*(1 - (2*i-1)*xi) =
+ // previous - (2*i-1)*xi - (2i-1)*(previous-1)*xi
+ // Such representation is convenient since previous contains only
+ // only x1...x[i-1]
+
+ fractions[i] = fractions[i-1];
+ fractions[i].push_back(-cl_I(2*i-1)); // - (2*i-1)*xi
+ // Now add -(2i+1)*(previous-1)*xi using formula
+ // 1/(z+j)*1/(z+i) = 1/(i-j)*(1/(z+j) - 1/(z+i)), that is
+ // xj*xi = 1/(i-j)(xj - xi)
+ for (size_t j=1; j < i; j++) {
+ cl_I curr = - exquo(cl_I(2*i-1)*fractions[i-1][j], cl_I(i-j));
+ fractions[i][j] = fractions[i][j] + curr;
+ fractions[i][i] = fractions[i][i] - curr;
+ // N.B. i-th polynomial has (i+1) non-zero coefficients