X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=9e24b99bc12e68bb2a5cdcbe263d4637bd9c64ea;hp=7d6ddca7cf3f4da9f0d212aa3dc3ebcda1621d3f;hb=fc80c36140aaf126150b2dea811177d8af35dac3;hpb=79eebbd4ebb7b5f7c6517298c728a16a68282ace diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 7d6ddca7..9e24b99b 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -3,8 +3,7 @@ * This file implements several functions that work on univariate and * multivariate polynomials and rational functions. * These functions include polynomial quotient and remainder, GCD and LCM - * computation, square-free factorization and rational function normalization. - */ + * computation, square-free factorization and rational function normalization. */ /* * GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany @@ -138,11 +137,11 @@ struct sym_desc { /** Lowest degree of symbol in polynomial "b" */ int ldeg_b; - /** Minimum of ldeg_a and ldeg_b (Used for sorting) */ - int min_deg; + /** Maximum of deg_a and deg_b (Used for sorting) */ + int max_deg; /** Commparison operator for sorting */ - bool operator<(const sym_desc &x) const {return min_deg < x.min_deg;} + bool operator<(const sym_desc &x) const {return max_deg < x.max_deg;} }; // Vector of sym_desc structures @@ -198,12 +197,21 @@ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v) int deg_b = b.degree(*(it->sym)); it->deg_a = deg_a; it->deg_b = deg_b; - it->min_deg = min(deg_a, deg_b); + it->max_deg = max(deg_a, deg_b); it->ldeg_a = a.ldegree(*(it->sym)); it->ldeg_b = b.ldegree(*(it->sym)); it++; } sort(v.begin(), v.end()); +#if 0 + clog << "Symbols:\n"; + it = v.begin(); itend = v.end(); + while (it != itend) { + clog << " " << *it->sym << ": deg_a=" << it->deg_a << ", deg_b=" << it->deg_b << ", ldeg_a=" << it->ldeg_a << ", ldeg_b=" << it->ldeg_b << ", max_deg=" << it->max_deg << endl; + clog << " lcoeff_a=" << a.lcoeff(*(it->sym)) << ", lcoeff_b=" << b.lcoeff(*(it->sym)) << endl; + it++; + } +#endif } @@ -504,6 +512,8 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) q = _ex0(); if (b.is_zero()) throw(std::overflow_error("divide: division by zero")); + if (a.is_zero()) + return true; if (is_ex_exactly_of_type(b, numeric)) { q = a / b; return true; @@ -828,8 +838,201 @@ ex ex::primpart(const symbol &x, const ex &c) const * GCD of multivariate polynomials */ +/** Compute GCD of multivariate polynomials using the Euclidean PRS algorithm + * (not really suited for multivariate GCDs). This function is only provided + * for testing purposes. + * + * @param a first multivariate polynomial + * @param b second multivariate polynomial + * @param x pointer to symbol (main variable) in which to compute the GCD in + * @return the GCD as a new expression + * @see gcd */ + +static ex eu_gcd(const ex &a, const ex &b, const symbol *x) +{ +//clog << "eu_gcd(" << a << "," << b << ")\n"; + + // Sort c and d so that c has higher degree + ex c, d; + int adeg = a.degree(*x), bdeg = b.degree(*x); + if (adeg >= bdeg) { + c = a; + d = b; + } else { + c = b; + d = a; + } + + // Euclidean algorithm + ex r; + for (;;) { +//clog << " d = " << d << endl; + r = rem(c, d, *x, false); + if (r.is_zero()) + return d.primpart(*x); + c = d; + d = r; + } +} + + +/** Compute GCD of multivariate polynomials using the Euclidean PRS algorithm + * with pseudo-remainders ("World's Worst GCD Algorithm", staying in Z[X]). + * This function is only provided for testing purposes. + * + * @param a first multivariate polynomial + * @param b second multivariate polynomial + * @param x pointer to symbol (main variable) in which to compute the GCD in + * @return the GCD as a new expression + * @see gcd */ + +static ex euprem_gcd(const ex &a, const ex &b, const symbol *x) +{ +//clog << "euprem_gcd(" << a << "," << b << ")\n"; + + // Sort c and d so that c has higher degree + ex c, d; + int adeg = a.degree(*x), bdeg = b.degree(*x); + if (adeg >= bdeg) { + c = a; + d = b; + } else { + c = b; + d = a; + } + + // Euclidean algorithm with pseudo-remainders + ex r; + for (;;) { +//clog << " d = " << d << endl; + r = prem(c, d, *x, false); + if (r.is_zero()) + return d.primpart(*x); + c = d; + d = r; + } +} + + +/** Compute GCD of multivariate polynomials using the primitive Euclidean + * PRS algorithm (complete content removal at each step). This function is + * only provided for testing purposes. + * + * @param a first multivariate polynomial + * @param b second multivariate polynomial + * @param x pointer to symbol (main variable) in which to compute the GCD in + * @return the GCD as a new expression + * @see gcd */ + +static ex peu_gcd(const ex &a, const ex &b, const symbol *x) +{ +//clog << "peu_gcd(" << a << "," << b << ")\n"; + + // Sort c and d so that c has higher degree + ex c, d; + int adeg = a.degree(*x), bdeg = b.degree(*x); + int ddeg; + if (adeg >= bdeg) { + c = a; + d = b; + ddeg = bdeg; + } else { + c = b; + d = a; + ddeg = adeg; + } + + // Remove content from c and d, to be attached to GCD later + ex cont_c = c.content(*x); + ex cont_d = d.content(*x); + ex gamma = gcd(cont_c, cont_d, NULL, NULL, false); + if (ddeg == 0) + return gamma; + c = c.primpart(*x, cont_c); + d = d.primpart(*x, cont_d); + + // Euclidean algorithm with content removal + ex r; + for (;;) { +//clog << " d = " << d << endl; + r = prem(c, d, *x, false); + if (r.is_zero()) + return gamma * d; + c = d; + d = r.primpart(*x); + } +} + + +/** Compute GCD of multivariate polynomials using the reduced PRS algorithm. + * This function is only provided for testing purposes. + * + * @param a first multivariate polynomial + * @param b second multivariate polynomial + * @param x pointer to symbol (main variable) in which to compute the GCD in + * @return the GCD as a new expression + * @see gcd */ + +static ex red_gcd(const ex &a, const ex &b, const symbol *x) +{ +//clog << "red_gcd(" << a << "," << b << ")\n"; + + // Sort c and d so that c has higher degree + ex c, d; + int adeg = a.degree(*x), bdeg = b.degree(*x); + int cdeg, ddeg; + if (adeg >= bdeg) { + c = a; + d = b; + cdeg = adeg; + ddeg = bdeg; + } else { + c = b; + d = a; + cdeg = bdeg; + ddeg = adeg; + } + + // Remove content from c and d, to be attached to GCD later + ex cont_c = c.content(*x); + ex cont_d = d.content(*x); + ex gamma = gcd(cont_c, cont_d, NULL, NULL, false); + if (ddeg == 0) + return gamma; + c = c.primpart(*x, cont_c); + d = d.primpart(*x, cont_d); + + // First element of subresultant sequence + ex r, ri = _ex1(); + int delta = cdeg - ddeg; + + for (;;) { + // Calculate polynomial pseudo-remainder +//clog << " d = " << d << endl; + r = prem(c, d, *x, false); + if (r.is_zero()) + return gamma * d.primpart(*x); + c = d; + cdeg = ddeg; + + if (!divide(r, pow(ri, delta), d, false)) + throw(std::runtime_error("invalid expression in red_gcd(), division failed")); + ddeg = d.degree(*x); + if (ddeg == 0) { + if (is_ex_exactly_of_type(r, numeric)) + return gamma; + else + return gamma * r.primpart(*x); + } + + ri = c.expand().lcoeff(*x); + delta = cdeg - ddeg; + } +} + + /** Compute GCD of multivariate polynomials using the subresultant PRS - * algorithm. This function is used internally gy gcd(). + * algorithm. This function is used internally by gcd(). * * @param a first multivariate polynomial * @param b second multivariate polynomial @@ -877,13 +1080,14 @@ static ex sr_gcd(const ex &a, const ex &b, const symbol *x) for (;;) { // Calculate polynomial pseudo-remainder //clog << " start of loop, psi = " << psi << ", calculating pseudo-remainder...\n"; +//clog << " d = " << d << endl; r = prem(c, d, *x, false); if (r.is_zero()) return gamma * d.primpart(*x); c = d; cdeg = ddeg; //clog << " dividing...\n"; - if (!divide(r, ri * power(psi, delta), d, false)) + if (!divide(r, ri * pow(psi, delta), d, false)) throw(std::runtime_error("invalid expression in sr_gcd(), division failed")); ddeg = d.degree(*x); if (ddeg == 0) { @@ -899,7 +1103,7 @@ static ex sr_gcd(const ex &a, const ex &b, const symbol *x) if (delta == 1) psi = ri; else if (delta) - divide(power(ri, delta), power(psi, delta-1), psi, false); + divide(pow(ri, delta), pow(psi, delta-1), psi, false); delta = cdeg - ddeg; } } @@ -1318,8 +1522,9 @@ factored_b: return g; } - // Try heuristic algorithm first, fall back to PRS if that failed ex g; +#if 1 + // Try heuristic algorithm first, fall back to PRS if that failed try { g = heur_gcd(aex, bex, ca, cb, var); } catch (gcdheu_failed) { @@ -1330,7 +1535,13 @@ factored_b: #if STATISTICS heur_gcd_failed++; #endif - g = sr_gcd(aex, bex, x); +#endif +// g = heur_gcd(aex, bex, ca, cb, var); +// g = eu_gcd(aex, bex, x); +// g = euprem_gcd(aex, bex, x); +// g = peu_gcd(aex, bex, x); +// g = red_gcd(aex, bex, x); + g = sr_gcd(aex, bex, x); if (g.is_equal(_ex1())) { // Keep cofactors factored if possible if (ca) @@ -1343,6 +1554,7 @@ factored_b: if (cb) divide(bex, g, *cb, false); } +#if 1 } else { if (g.is_equal(_ex1())) { // Keep cofactors factored if possible @@ -1351,8 +1563,9 @@ factored_b: if (cb) *cb = b; } - return g; } +#endif + return g; } @@ -1460,8 +1673,8 @@ ex sqrfree(const ex &a, const symbol &x) */ /** Create a symbol for replacing the expression "e" (or return a previously - * assigned symbol). The symbol is appended to sym_list and returned, the - * expression is appended to repl_list. + * assigned symbol). The symbol is appended to sym_lst and returned, the + * expression is appended to repl_lst. * @see ex::normal */ static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst) { @@ -1481,6 +1694,26 @@ static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst) return es; } +/** Create a symbol for replacing the expression "e" (or return a previously + * assigned symbol). An expression of the form "symbol == expression" is added + * to repl_lst and the symbol is returned. + * @see ex::to_rational */ +static ex replace_with_symbol(const ex &e, lst &repl_lst) +{ + // Expression already in repl_lst? Then return the assigned symbol + for (unsigned i=0; isetflag(status_flags::dynallocated); +} + + /** Normalization of rational functions. * This function converts an expression to its normal form * "numerator/denominator", where numerator and denominator are (relatively * prime) polynomials. Any subexpressions which are not rational functions - * (like non-rational numbers, non-integer powers or functions like Sin(), - * Cos() etc.) are replaced by temporary symbols which are re-substituted by + * (like non-rational numbers, non-integer powers or functions like sin(), + * cos() etc.) are replaced by temporary symbols which are re-substituted by * the (normalized) subexpressions before normal() returns (this way, any * expression can be treated as a rational function). normal() is applied * recursively to arguments of functions etc. @@ -1810,6 +2051,75 @@ ex ex::denom(void) const return e.op(1); } + +/** Default implementation of ex::to_rational(). It replaces the object with a + * temporary symbol. + * @see ex::to_rational */ +ex basic::to_rational(lst &repl_lst) const +{ + return replace_with_symbol(*this, repl_lst); +} + + +/** Implementation of ex::to_rational() for symbols. This returns the unmodified symbol. + * @see ex::to_rational */ +ex symbol::to_rational(lst &repl_lst) const +{ + return *this; +} + + +/** Implementation of ex::to_rational() for a numeric. It splits complex numbers + * into re+I*im and replaces I and non-rational real numbers with a temporary + * symbol. + * @see ex::to_rational */ +ex numeric::to_rational(lst &repl_lst) const +{ + if (is_real()) { + if (!is_integer()) + return replace_with_symbol(*this, repl_lst); + } else { // complex + numeric re = real(), im = imag(); + ex re_ex = re.is_rational() ? re : replace_with_symbol(re, repl_lst); + ex im_ex = im.is_rational() ? im : replace_with_symbol(im, repl_lst); + return re_ex + im_ex * replace_with_symbol(I, repl_lst); + } + return *this; +} + + +/** Implementation of ex::to_rational() for powers. It replaces non-integer + * powers by temporary symbols. + * @see ex::to_rational */ +ex power::to_rational(lst &repl_lst) const +{ + if (exponent.info(info_flags::integer)) + return power(basis.to_rational(repl_lst), exponent); + else + return replace_with_symbol(*this, repl_lst); +} + + +/** Rationalization of non-rational functions. + * This function converts a general expression to a rational polynomial + * by replacing all non-rational subexpressions (like non-rational numbers, + * non-integer powers or functions like sin(), cos() etc.) to temporary + * symbols. This makes it possible to use functions like gcd() and divide() + * on non-rational functions by applying to_rational() on the arguments, + * calling the desired function and re-substituting the temporary symbols + * in the result. To make the last step possible, all temporary symbols and + * their associated expressions are collected in the list specified by the + * repl_lst parameter in the form {symbol == expression}, ready to be passed + * as an argument to ex::subs(). + * + * @param repl_lst collects a list of all temporary symbols and their replacements + * @return rationalized expression */ +ex ex::to_rational(lst &repl_lst) const +{ + return bp->to_rational(repl_lst); +} + + #ifndef NO_NAMESPACE_GINAC } // namespace GiNaC #endif // ndef NO_NAMESPACE_GINAC