From: Richard Kreckel Date: Sun, 11 Mar 2001 19:10:03 +0000 (+0000) Subject: normal.cpp, normal.h: Fix Yun's algorithm (there was a mistake in Geddes et X-Git-Tag: release_0-8-0~29 X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=commitdiff_plain;h=fae800e3a2dfe04b49bac90cfd7ab4e155e865e6;hp=e13c691c25abd739163d30de5e6d038977683690 normal.cpp, normal.h: Fix Yun's algorithm (there was a mistake in Geddes et al.) for square-free factorization. Also change the syntax to allow for square-free factoring in a recursive way in Q[X]. This is an example how it works in Z[x,y,z] in ginsh: > ?sqrfree sqrfree(expression [, symbol-list]) - square-free factorization of a polynomial > sqrfree(expand((x+y)*(x+y+z)^2*(y-z)^3)); (x+y)*(-x-y-z)^2*(y^3+3*y*z^2-3*y^2*z-z^3) > sqrfree(expand((x+y)*(x+y+z)^2*(y-z)^3),[y]); (x+y)*(y-z)^3*(x+y+z)^2 > sqrfree(expand((x+y)*(x+y+z)^2*(y-z)^3),[z]); (x+y)*(y-z)^3*(x+y+z)^2 > sqrfree(expand((x+y)*(x+y+z)^2*(y-z)^3),[x]); (x+y)*(-x-y-z)^2*(y^3+3*y*z^2-3*y^2*z-z^3) > sqrfree(expand((x+y)*(x+y+z)^2*(y-z)^3),[x,y]); (x+y)*(-x-y-z)^2*(y^3+3*y*z^2-3*y^2*z-z^3) > sqrfree(expand((x+y)*(x+y+z)^2*(y-z)^3),[x,y,z]); (x+y)*(-x-y-z)^2*(y^3+3*y*z^2-3*y^2*z-z^3) > sqrfree(expand((x+y)*(x+y+z)^2*(y-z)^3),[v]); 3*x*y^5+3*x^2*y^4+x^3*y^3+y^6-y*z^5-2*y^4*z^2+y^2*z^4-y^5*z+2*y^3*z^3-3*x^3*y^2*z+3*x^3*y*z^2-5*x*y^4*z-x*y*z^4+3*x^2*y*z^3+6*x*y^2*z^3+3*x^2*y^2*z^2-2*x*y^3*z^2-7*x^2*y^3*z-x^3*z^3-2*x^2*z^4-x*z^5 --- diff --git a/NEWS b/NEWS index 2b269c3e..2e80bc4f 100644 --- a/NEWS +++ b/NEWS @@ -11,6 +11,7 @@ This file records noteworthy changes. indexed objects like (a+b).i -> a.i + b.i * Renamed get_indices() to get_free_indices(), which no longer returns dummy indices and checks the consistency of indices in sums. +* sqrfree() factorization fixed and improved syntactically. * subs() works on matrices. 0.7.3 (28 February 2001) diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 669a59f8..e9d883f5 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1691,70 +1691,85 @@ ex lcm(const ex &a, const ex &b, bool check_args) * 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; i0) { + 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(); } diff --git a/ginac/normal.h b/ginac/normal.h index b461dcf0..fb6960f4 100644 --- a/ginac/normal.h +++ b/ginac/normal.h @@ -26,6 +26,8 @@ #ifndef __GINAC_NORMAL_H__ #define __GINAC_NORMAL_H__ +#include "lst.h" + namespace GiNaC { class ex; @@ -50,7 +52,7 @@ extern ex gcd(const ex &a, const ex &b, ex *ca = NULL, ex *cb = NULL, bool check extern ex lcm(const ex &a, const ex &b, bool check_args = true); // Square-free factorization of a polynomial a(x) -extern ex sqrfree(const ex &a, const symbol &x); +extern ex sqrfree(const ex &a, const lst &l = lst()); } // namespace GiNaC diff --git a/ginsh/ginsh.1.in b/ginsh/ginsh.1.in index 7857f662..9be13b43 100644 --- a/ginsh/ginsh.1.in +++ b/ginsh/ginsh.1.in @@ -319,7 +319,7 @@ detail here. Please refer to the GiNaC documentation. .BI series( expression ", " relation-or-symbol ", " order ) \- series expansion .br -.BI sqrfree( expression ", " symbol ) +.BI sqrfree( "expression [" ", " symbol-list] ) \- square-free factorization of a polynomial .br .BI sqrt( expression ) diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy index a83bbc32..784af1cb 100644 --- a/ginsh/ginsh_parser.yy +++ b/ginsh/ginsh_parser.yy @@ -442,10 +442,15 @@ static ex f_series(const exprseq &e) 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) @@ -533,7 +538,8 @@ static const fcn_init builtin_fcns[] = { {"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)},