+ * \note
+ * A polynomial \f$p(X) \in C[X]\f$ is said <EM>square-free</EM>
+ * if, whenever any two polynomials \f$q(X)\f$ and \f$r(X)\f$
+ * are such that
+ * \f[
+ * p(X) = q(X)^2 r(X),
+ * \f]
+ * we have \f$q(X) \in C\f$.
+ * This means that \f$p(X)\f$ has no repeated factors, apart
+ * eventually from constants.
+ * Given a polynomial \f$p(X) \in C[X]\f$, we say that the
+ * decomposition
+ * \f[
+ * p(X) = b \cdot p_1(X)^{a_1} \cdot p_2(X)^{a_2} \cdots p_r(X)^{a_r}
+ * \f]
+ * is a <EM>square-free factorization</EM> of \f$p(X)\f$ if the
+ * following conditions hold:
+ * -# \f$b \in C\f$ and \f$b \neq 0\f$;
+ * -# \f$a_i\f$ is a positive integer for \f$i = 1, \ldots, r\f$;
+ * -# the degree of the polynomial \f$p_i\f$ is strictly positive
+ * for \f$i = 1, \ldots, r\f$;
+ * -# the polynomial \f$\Pi_{i=1}^r p_i(X)\f$ is square-free.
+ *
+ * Square-free factorizations need not be unique. For example, if
+ * \f$a_i\f$ is even, we could change the polynomial \f$p_i(X)\f$
+ * into \f$-p_i(X)\f$.
+ * Observe also that the factors \f$p_i(X)\f$ need not be irreducible
+ * polynomials.
+ */
+ex sqrfree(const ex &a, const lst &l)
+{
+ if (is_exactly_a<numeric>(a) || // algorithm does not trap a==0
+ is_a<symbol>(a)) // shortcut
+ return a;
+
+ // If no lst of variables to factorize in was specified we have to
+ // invent one now. Maybe one can optimize here by reversing the order
+ // or so, I don't know.
+ lst args;
+ if (l.nops()==0) {
+ sym_desc_vec sdv;
+ get_symbol_stats(a, _ex0, sdv);
+ sym_desc_vec::const_iterator it = sdv.begin(), itend = sdv.end();
+ while (it != itend) {
+ args.append(*it->sym);
+ ++it;
+ }
+ } else {
+ args = l;
+ }
+
+ // Find the symbol to factor in at this stage
+ if (!is_ex_of_type(args.op(0), symbol))
+ throw (std::runtime_error("sqrfree(): invalid factorization variable"));
+ const symbol &x = ex_to<symbol>(args.op(0));
+
+ // convert the argument from something in Q[X] to something in Z[X]
+ const numeric lcm = lcm_of_coefficients_denominators(a);
+ const ex tmp = multiply_lcm(a,lcm);
+
+ // find the factors
+ exvector factors = sqrfree_yun(tmp,x);
+
+ // construct the next list of symbols with the first element popped
+ lst newargs = args;
+ newargs.remove_first();
+
+ // recurse down the factors in remaining variables
+ if (newargs.nops()>0) {
+ exvector::iterator i = factors.begin();
+ while (i != factors.end()) {
+ *i = sqrfree(*i, newargs);
+ ++i;
+ }
+ }
+
+ // Done with recursion, now construct the final result
+ ex result = _ex1;
+ exvector::const_iterator it = factors.begin(), itend = factors.end();
+ for (int p = 1; it!=itend; ++it, ++p)
+ result *= power(*it, p);
+
+ // Yun's algorithm does not account for constant factors. (For univariate
+ // polynomials it works only in the monic case.) We can correct this by
+ // inserting what has been lost back into the result. For completeness
+ // we'll also have to recurse down that factor in the remaining variables.
+ if (newargs.nops()>0)
+ result *= sqrfree(quo(tmp, result, x), newargs);
+ else
+ result *= quo(tmp, result, x);
+
+ // Put in the reational overall factor again and return
+ return result * lcm.inverse();
+}
+
+
+/** Compute square-free partial fraction decomposition of rational function
+ * a(x).
+ *
+ * @param a rational function over Z[x], treated as univariate polynomial
+ * in x
+ * @param x variable to factor in
+ * @return decomposed rational function */
+ex sqrfree_parfrac(const ex & a, const symbol & x)
+{
+ // Find numerator and denominator
+ ex nd = numer_denom(a);
+ ex numer = nd.op(0), denom = nd.op(1);
+//clog << "numer = " << numer << ", denom = " << denom << endl;
+
+ // Convert N(x)/D(x) -> Q(x) + R(x)/D(x), so degree(R) < degree(D)
+ ex red_poly = quo(numer, denom, x), red_numer = rem(numer, denom, x).expand();
+//clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << endl;
+
+ // Factorize denominator and compute cofactors
+ exvector yun = sqrfree_yun(denom, x);
+//clog << "yun factors: " << exprseq(yun) << endl;
+ unsigned num_yun = yun.size();
+ exvector factor; factor.reserve(num_yun);
+ exvector cofac; cofac.reserve(num_yun);
+ for (unsigned i=0; i<num_yun; i++) {
+ if (!yun[i].is_equal(_ex1)) {
+ for (unsigned j=0; j<=i; j++) {
+ factor.push_back(pow(yun[i], j+1));
+ ex prod = _ex1;
+ for (unsigned 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());
+ }
+ }
+ }
+ unsigned num_factors = factor.size();
+//clog << "factors : " << exprseq(factor) << endl;
+//clog << "cofactors: " << exprseq(cofac) << endl;
+
+ // Construct coefficient matrix for decomposition
+ int max_denom_deg = denom.degree(x);
+ matrix sys(max_denom_deg + 1, num_factors);
+ matrix rhs(max_denom_deg + 1, 1);
+ for (int i=0; i<=max_denom_deg; i++) {
+ for (unsigned j=0; j<num_factors; j++)
+ sys(i, j) = cofac[j].coeff(x, i);
+ rhs(i, 0) = red_numer.coeff(x, i);
+ }
+//clog << "coeffs: " << sys << endl;
+//clog << "rhs : " << rhs << endl;
+
+ // Solve resulting linear system
+ matrix vars(num_factors, 1);
+ for (unsigned i=0; i<num_factors; i++)
+ vars(i, 0) = symbol();
+ matrix sol = sys.solve(vars, rhs);
+
+ // Sum up decomposed fractions
+ ex sum = 0;
+ for (unsigned i=0; i<num_factors; i++)
+ sum += sol(i, 0) / factor[i];
+
+ return red_poly + sum;