From d7d0bcda91b647db9588f3aa1a465f1570d088c4 Mon Sep 17 00:00:00 2001 From: Christian Bauer Date: Fri, 20 Aug 2004 20:14:10 +0000 Subject: [PATCH] - added ex::unitcontprim() - fixed a problem where content() would return something that was not unit normal --- doc/tutorial/ginac.texi | 17 +++++- ginac/ex.h | 1 + ginac/normal.cpp | 122 +++++++++++++++++++++++++++------------- 3 files changed, 100 insertions(+), 40 deletions(-) diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index b36b57fd..35dbebc7 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -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 diff --git a/ginac/ex.h b/ginac/ex.h index 13177e28..b274b49e 100644 --- a/ginac/ex.h +++ b/ginac/ex.h @@ -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; diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 67fbf936..47fec508 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -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(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(*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(*this)) - return _ex1; - - ex c = content(x); - if (c.is_zero()) - return _ex0; - ex u = unit(x); - if (is_exactly_a(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(*this)) return _ex1; + // Divide by unit and content to get primitive part ex u = unit(x); if (is_exactly_a(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(*this)) { + if (info(info_flags::negative)) { + u = _ex_1; + c = abs(ex_to(*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(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; } -- 2.38.1