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 2b269c3e4e5df26e8ac596b492c3358989fc04f0..2e80bc4f36f462d43a0fd9a2a4f337a50b65d94b 100644 (file)
--- 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.
   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)
 * subs() works on matrices.
 
 0.7.3 (28 February 2001)

diff --git a/ginac/normal.cpp b/ginac/normal.cpp
index 669a59f8de12d0a7cfb3d3eba694eaa9703d275d..e9d883f53b12a74af830f2be921c7c8a6e3296ab 100644 (file)
--- 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
  */
 
  *  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 {
        } 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();

 }
 
 
 }
 
 

diff --git a/ginac/normal.h b/ginac/normal.h
index b461dcf0723407930f2b9ca9e9d1ade6f50e9bc5..fb6960f41f7fd13a813a30c507f202ce75fe371a 100644 (file)
--- a/ginac/normal.h
+++ b/ginac/normal.h
@@ -26,6 +26,8 @@
 #ifndef __GINAC_NORMAL_H__
 #define __GINAC_NORMAL_H__
 
 #ifndef __GINAC_NORMAL_H__
 #define __GINAC_NORMAL_H__
 

+#include "lst.h"
+

 namespace GiNaC {
 
 class ex;
 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 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
 
 
 } // namespace GiNaC
 

diff --git a/ginsh/ginsh.1.in b/ginsh/ginsh.1.in
index 7857f662788093f3b57f5fd68316eaf786ef7a62..9be13b4350461f0e3d1a75e8ec2706e30946e61c 100644 (file)
--- 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 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 )
 \- square-free factorization of a polynomial
 .br
 .BI sqrt( expression )

diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy
index a83bbc32a444ecbd7ee0cede115558c16876b95f..784af1cb22902a1724aaf10dad7a3c853d6ccd15 100644 (file)
--- 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());
 }
 
        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)
 }
 
 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)},
        {"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)},
        {"sqrt", fcn_desc(f_sqrt, 1)},
        {"subs", fcn_desc(f_subs2, 2)},
        {"subs", fcn_desc(f_subs3, 3)},