* Square-free factorization
*/
-// Univariate GCD of polynomials in Q[x] (used internally by sqrfree()).
-// a and b can be multivariate polynomials but they are treated as univariate polynomials in x.
-static ex univariate_gcd(const ex &a, const ex &b, const symbol &x)
+/** Compute square-free factorization of multivariate polynomial a(x) using
+ * Yun´s algorithm. Used internally by sqrfree().
+ *
+ * @param a multivariate polynomial over Z[X], treated here as univariate
+ * polynomial in x.
+ * @param x variable to factor in
+ * @return vector of factors sorted in ascending degree */
+exvector sqrfree_yun(const ex &a, const symbol &x)
{
- if (a.is_zero())
- return b;
- if (b.is_zero())
- return a;
- if (a.is_equal(_ex1()) || b.is_equal(_ex1()))
- return _ex1();
- if (is_ex_of_type(a, numeric) && is_ex_of_type(b, numeric))
- return gcd(ex_to_numeric(a), ex_to_numeric(b));
- if (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))
- throw(std::invalid_argument("univariate_gcd: arguments must be polynomials over the rationals"));
-
- // Euclidean algorithm
- ex c, d, r;
- if (a.degree(x) >= b.degree(x)) {
- c = a;
- d = b;
- } else {
- c = b;
- d = a;
- }
- for (;;) {
- r = rem(c, d, x, false);
- if (r.is_zero())
- break;
- c = d;
- d = r;
- }
- return d / d.lcoeff(x);
+ int i = 0;
+ exvector res;
+ ex w = a;
+ ex z = w.diff(x);
+ ex g = gcd(w, z);
+ if (g.is_equal(_ex1())) {
+ res.push_back(a);
+ return res;
+ }
+ ex y;
+ do {
+ w = quo(w, g, x);
+ y = quo(z, g, x);
+ z = y - w.diff(x);
+ g = gcd(w, z);
+ res.push_back(g);
+ ++i;
+ } while (!z.is_zero());
+ return res;
}
-
-
-/** Compute square-free factorization of multivariate polynomial a(x) using
- * Yun´s algorithm.
+/** Compute square-free factorization of multivariate polynomial in Q[X].
*
- * @param a multivariate polynomial
- * @param x variable to factor in
- * @return factored polynomial */
-ex sqrfree(const ex &a, const symbol &x)
+ * @param a multivariate polynomial over Q[X]
+ * @param x lst of variables to factor in, may be left empty for autodetection
+ * @return polynomail a in square-free factored form. */
+ex sqrfree(const ex &a, const lst &l)
{
- int i = 1;
- ex res = _ex1();
- ex b = a.diff(x);
- ex c = univariate_gcd(a, b, x);
- ex w;
- if (c.is_equal(_ex1())) {
- w = a;
+ if (is_ex_of_type(a,numeric) || // algorithm does not trap a==0
+ is_ex_of_type(a,symbol)) // 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);
+ for (sym_desc_vec::iterator it=sdv.begin(); it!=sdv.end(); ++it)
+ args.append(*it->sym);
} else {
- w = quo(a, c, x);
- ex y = quo(b, c, x);
- ex z = y - w.diff(x);
- while (!z.is_zero()) {
- ex g = univariate_gcd(w, z, x);
- res *= power(g, i);
- w = quo(w, g, x);
- y = quo(z, g, x);
- z = y - w.diff(x);
- i++;
- }
- }
- return res * power(w, i);
+ 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]
+ numeric lcm = lcm_of_coefficients_denominators(a);
+ 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;
+ for (int i=1; i<args.nops(); ++i)
+ newargs.append(args.op(i));
+ // recurse down the factors in remaining vars
+ if (newargs.nops()>0) {
+ for (exvector::iterator i=factors.begin(); i!=factors.end(); ++i)
+ *i = sqrfree(*i, newargs);
+ }
+ // Done with recursion, now construct the final result
+ ex result = _ex1();
+ exvector::iterator it = factors.begin();
+ for (int p = 1; it!=factors.end(); ++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:
+ result = result * quo(tmp, result, x);
+ return result * lcm.inverse();
}
return e[0].series(e[1], ex_to_numeric(e[2]).to_int());
}
-static ex f_sqrfree(const exprseq &e)
+static ex f_sqrfree1(const exprseq &e)
{
- CHECK_ARG(1, symbol, sqrfree);
- return sqrfree(e[0], ex_to_symbol(e[1]));
+ return sqrfree(e[0]);
+}
+
+static ex f_sqrfree2(const exprseq &e)
+{
+ CHECK_ARG(1, lst, sqrfree);
+ return sqrfree(e[0], ex_to_lst(e[1]));
}
static ex f_subs3(const exprseq &e)
{"quo", fcn_desc(f_quo, 3)},
{"rem", fcn_desc(f_rem, 3)},
{"series", fcn_desc(f_series, 3)},
- {"sqrfree", fcn_desc(f_sqrfree, 2)},
+ {"sqrfree", fcn_desc(f_sqrfree1, 1)},
+ {"sqrfree", fcn_desc(f_sqrfree2, 2)},
{"sqrt", fcn_desc(f_sqrt, 1)},
{"subs", fcn_desc(f_subs2, 2)},
{"subs", fcn_desc(f_subs3, 3)},