- added ex::unitcontprim()
authorChristian Bauer <Christian.Bauer@uni-mainz.de>
Fri, 20 Aug 2004 20:14:10 +0000 (20:14 +0000)
committerChristian Bauer <Christian.Bauer@uni-mainz.de>
Fri, 20 Aug 2004 20:14:10 +0000 (20:14 +0000)
- fixed a problem where content() would return something that was not unit
  normal

doc/tutorial/ginac.texi
ginac/ex.h
ginac/normal.cpp

index b36b57fd2c3d6912190ac450db5896b2d894fc6e..35dbebc76b6c0eb0c553ddb1d373244dfa28e7ed 100644 (file)
@@ -4725,6 +4725,7 @@ in which case the value of @code{q} is undefined.
 @cindex @code{unit()}
 @cindex @code{content()}
 @cindex @code{primpart()}
+@cindex @code{unitcontprim()}
 
 The methods
 
@@ -4732,14 +4733,26 @@ The methods
 ex ex::unit(const ex & x);
 ex ex::content(const ex & x);
 ex ex::primpart(const ex & x);
+ex ex::primpart(const ex & x, const ex & c);
 @end example
 
 return the unit part, content part, and primitive polynomial of a multivariate
 polynomial with respect to the variable @samp{x} (the unit part being the sign
 of the leading coefficient, the content part being the GCD of the coefficients,
 and the primitive polynomial being the input polynomial divided by the unit and
-content parts). The product of unit, content, and primitive part is the
-original polynomial.
+content parts). The second variant of @code{primpart()} expects the previously
+calculated content part of the polynomial in @code{c}, which enables it to
+work faster in the case where the content part has already been computed. The
+product of unit, content, and primitive part is the original polynomial.
+
+Additionally, the method
+
+@example
+void ex::unitcontprim(const ex & x, ex & u, ex & c, ex & p);
+@end example
+
+computes the unit, content, and primitive parts in one go, returning them
+in @code{u}, @code{c}, and @code{p}, respectively.
 
 
 @subsection GCD, LCM and resultant
index 13177e28b9d6613fdeef29d6ebcdc51fc2193421..b274b49ee2c72fb05d5c3747a20a88fc3f33a248 100644 (file)
@@ -192,6 +192,7 @@ public:
        numeric integer_content() const;
        ex primpart(const ex &x) const;
        ex primpart(const ex &x, const ex &cont) const;
+       void unitcontprim(const ex &x, ex &u, ex &c, ex &p) const;
        ex smod(const numeric &xi) const { return bp->smod(xi); }
        numeric max_coefficient() const;
 
index 67fbf93612beaa629ed2df6a987152e56a546374..47fec5088ab59d62dfdefc8b4b7c851e40acd17e 100644 (file)
@@ -819,14 +819,14 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite
  *  polynomial in Q[x]. The product of unit part, content part, and primitive
  *  part is the polynomial itself.
  *
- *  @param x  variable in which to compute the unit part
+ *  @param x  main variable
  *  @return unit part
- *  @see ex::content, ex::primpart */
+ *  @see ex::content, ex::primpart, ex::unitcontprim */
 ex ex::unit(const ex &x) const
 {
        ex c = expand().lcoeff(x);
        if (is_exactly_a<numeric>(c))
-               return c < _ex0 ? _ex_1 : _ex1;
+               return c.info(info_flags::negative) ?_ex_1 : _ex1;
        else {
                ex y;
                if (get_first_symbol(c, y))
@@ -841,15 +841,14 @@ ex ex::unit(const ex &x) const
  *  multivariate polynomial in Q[x]. The product of unit part, content part,
  *  and primitive part is the polynomial itself.
  *
- *  @param x  variable in which to compute the content part
+ *  @param x  main variable
  *  @return content part
- *  @see ex::unit, ex::primpart */
+ *  @see ex::unit, ex::primpart, ex::unitcontprim */
 ex ex::content(const ex &x) const
 {
-       if (is_zero())
-               return _ex0;
        if (is_exactly_a<numeric>(*this))
                return info(info_flags::negative) ? -*this : *this;
+
        ex e = expand();
        if (e.is_zero())
                return _ex0;
@@ -858,15 +857,15 @@ ex ex::content(const ex &x) const
        // If the leading coefficient of the quotient is an integer, we are done.
        ex c = e.integer_content();
        ex r = e / c;
-       ex lcoeff = r.lcoeff(x);
+       int deg = r.degree(x);
+       ex lcoeff = r.coeff(x, deg);
        if (lcoeff.info(info_flags::integer))
                return c;
 
        // GCD of all coefficients
-       int deg = r.degree(x);
        int ldeg = r.ldegree(x);
        if (deg == ldeg)
-               return lcoeff * c;
+               return lcoeff * c / lcoeff.unit(x);
        ex cont = _ex0;
        for (int i=ldeg; i<=deg; i++)
                cont = gcd(r.coeff(x, i), cont, NULL, NULL, false);
@@ -874,28 +873,19 @@ ex ex::content(const ex &x) const
 }
 
 
-/** Compute primitive part of a multivariate polynomial in Q[x].
- *  The product of unit part, content part, and primitive part is the
- *  polynomial itself.
+/** Compute primitive part of a multivariate polynomial in Q[x]. The result
+ *  will be a unit-normal polynomial with a content part of 1. The product
+ *  of unit part, content part, and primitive part is the polynomial itself.
  *
- *  @param x  variable in which to compute the primitive part
+ *  @param x  main variable
  *  @return primitive part
- *  @see ex::unit, ex::content */
+ *  @see ex::unit, ex::content, ex::unitcontprim */
 ex ex::primpart(const ex &x) const
 {
-       if (is_zero())
-               return _ex0;
-       if (is_exactly_a<numeric>(*this))
-               return _ex1;
-
-       ex c = content(x);
-       if (c.is_zero())
-               return _ex0;
-       ex u = unit(x);
-       if (is_exactly_a<numeric>(c))
-               return *this / (c * u);
-       else
-               return quo(*this, c * u, x, false);
+       // We need to compute the unit and content anyway, so call unitcontprim()
+       ex u, c, p;
+       unitcontprim(x, u, c, p);
+       return p;
 }
 
 
@@ -903,18 +893,17 @@ ex ex::primpart(const ex &x) const
  *  content part is already known. This function is faster in computing the
  *  primitive part than the previous function.
  *
- *  @param x  variable in which to compute the primitive part
+ *  @param x  main variable
  *  @param c  previously computed content part
  *  @return primitive part */
 ex ex::primpart(const ex &x, const ex &c) const
 {
-       if (is_zero())
-               return _ex0;
-       if (c.is_zero())
+       if (is_zero() || c.is_zero())
                return _ex0;
        if (is_exactly_a<numeric>(*this))
                return _ex1;
 
+       // Divide by unit and content to get primitive part
        ex u = unit(x);
        if (is_exactly_a<numeric>(c))
                return *this / (c * u);
@@ -923,6 +912,61 @@ ex ex::primpart(const ex &x, const ex &c) const
 }
 
 
+/** Compute unit part, content part, and primitive part of a multivariate
+ *  polynomial in Q[x]. The product of the three parts is the polynomial
+ *  itself.
+ *
+ *  @param x  main variable
+ *  @param u  unit part (returned)
+ *  @param c  content part (returned)
+ *  @param p  primitive part (returned)
+ *  @see ex::unit, ex::content, ex::primpart */
+void ex::unitcontprim(const ex &x, ex &u, ex &c, ex &p) const
+{
+       // Quick check for zero (avoid expanding)
+       if (is_zero()) {
+               u = _ex1;
+               c = p = _ex0;
+               return;
+       }
+
+       // Special case: input is a number
+       if (is_exactly_a<numeric>(*this)) {
+               if (info(info_flags::negative)) {
+                       u = _ex_1;
+                       c = abs(ex_to<numeric>(*this));
+               } else {
+                       u = _ex1;
+                       c = *this;
+               }
+               p = _ex1;
+               return;
+       }
+
+       // Expand input polynomial
+       ex e = expand();
+       if (e.is_zero()) {
+               u = _ex1;
+               c = p = _ex0;
+               return;
+       }
+
+       // Compute unit and content
+       u = unit(x);
+       c = content(x);
+
+       // Divide by unit and content to get primitive part
+       if (c.is_zero()) {
+               p = _ex0;
+               return;
+       }
+       if (is_exactly_a<numeric>(c))
+               p = *this / (c * u);
+       else
+               p = quo(e, c * u, x, false);
+}
+
+
 /*
  *  GCD of multivariate polynomials
  */
@@ -1396,16 +1440,18 @@ factored_b:
 
        // Try to eliminate variables
        if (var->deg_a == 0) {
-               ex c = bex.content(x);
-               ex g = gcd(aex, c, ca, cb, false);
+               ex bex_u, bex_c, bex_p;
+               bex.unitcontprim(x, bex_u, bex_c, bex_p);
+               ex g = gcd(aex, bex_c, ca, cb, false);
                if (cb)
-                       *cb *= bex.unit(x) * bex.primpart(x, c);
+                       *cb *= bex_u * bex_p;
                return g;
        } else if (var->deg_b == 0) {
-               ex c = aex.content(x);
-               ex g = gcd(c, bex, ca, cb, false);
+               ex aex_u, aex_c, aex_p;
+               aex.unitcontprim(x, aex_u, aex_c, aex_p);
+               ex g = gcd(aex_c, bex, ca, cb, false);
                if (ca)
-                       *ca *= aex.unit(x) * aex.primpart(x, c);
+                       *ca *= aex_u * aex_p;
                return g;
        }