X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=074759df4db92413e902b5ec6182af426dde48b6;hp=852400608fe263832822931758c949725c80fe06;hb=e7a89b57b362bf02a13c0e72292ff490924a550c;hpb=f293ecba8b6026a7754795256b2f23910bf70507 diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 85240060..074759df 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 @@ -57,7 +56,34 @@ namespace GiNaC { #define FAST_COMPARE 1 // Set this if you want divide_in_z() to use remembering -#define USE_REMEMBER 1 +#define USE_REMEMBER 0 + +// Set this if you want divide_in_z() to use trial division followed by +// polynomial interpolation (usually slower except for very large problems) +#define USE_TRIAL_DIVISION 0 + +// Set this to enable some statistical output for the GCD routines +#define STATISTICS 0 + + +#if STATISTICS +// Statistics variables +static int gcd_called = 0; +static int sr_gcd_called = 0; +static int heur_gcd_called = 0; +static int heur_gcd_failed = 0; + +// Print statistics at end of program +static struct _stat_print { + _stat_print() {} + ~_stat_print() { + cout << "gcd() called " << gcd_called << " times\n"; + cout << "sr_gcd() called " << sr_gcd_called << " times\n"; + cout << "heur_gcd() called " << heur_gcd_called << " times\n"; + cout << "heur_gcd() failed " << heur_gcd_failed << " times\n"; + } +} stat_print; +#endif /** Return pointer to first symbol found in expression. Due to GiNaC´s @@ -67,7 +93,6 @@ namespace GiNaC { * @param e expression to search * @param x pointer to first symbol found (returned) * @return "false" if no symbol was found, "true" otherwise */ - static bool get_first_symbol(const ex &e, const symbol *&x) { if (is_ex_exactly_of_type(e, symbol)) { @@ -111,11 +136,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 @@ -160,7 +185,6 @@ static void collect_symbols(const ex &e, sym_desc_vec &v) * @param a first multivariate polynomial * @param b second multivariate polynomial * @param v vector of sym_desc structs (filled in) */ - static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v) { collect_symbols(a.eval(), v); // eval() to expand assigned symbols @@ -171,12 +195,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 } @@ -212,7 +245,6 @@ static numeric lcmcoeff(const ex &e, const numeric &l) * * @param e multivariate polynomial (need not be expanded) * @return LCM of denominators of coefficients */ - static numeric lcm_of_coefficients_denominators(const ex &e) { return lcmcoeff(e, _num1()); @@ -223,7 +255,6 @@ static numeric lcm_of_coefficients_denominators(const ex &e) * * @param e multivariate polynomial (need not be expanded) * @param lcm LCM to multiply in */ - static ex multiply_lcm(const ex &e, const numeric &lcm) { if (is_ex_exactly_of_type(e, mul)) { @@ -253,7 +284,6 @@ static ex multiply_lcm(const ex &e, const numeric &lcm) * * @param e expanded polynomial * @return integer content */ - numeric ex::integer_content(void) const { GINAC_ASSERT(bp!=0); @@ -314,7 +344,6 @@ numeric mul::integer_content(void) const * @param check_args check whether a and b are polynomials with rational * coefficients (defaults to "true") * @return quotient of a and b in Q[x] */ - ex quo(const ex &a, const ex &b, const symbol &x, bool check_args) { if (b.is_zero()) @@ -365,7 +394,6 @@ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args) * @param check_args check whether a and b are polynomials with rational * coefficients (defaults to "true") * @return remainder of a(x) and b(x) in Q[x] */ - ex rem(const ex &a, const ex &b, const symbol &x, bool check_args) { if (b.is_zero()) @@ -417,7 +445,6 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args) * @param check_args check whether a and b are polynomials with rational * coefficients (defaults to "true") * @return pseudo-remainder of a(x) and b(x) in Z[x] */ - ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) { if (b.is_zero()) @@ -462,6 +489,57 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) } +/** Sparse pseudo-remainder of polynomials a(x) and b(x) in Z[x]. + * + * @param a first polynomial in x (dividend) + * @param b second polynomial in x (divisor) + * @param x a and b are polynomials in x + * @param check_args check whether a and b are polynomials with rational + * coefficients (defaults to "true") + * @return sparse pseudo-remainder of a(x) and b(x) in Z[x] */ + +ex sprem(const ex &a, const ex &b, const symbol &x, bool check_args) +{ + if (b.is_zero()) + throw(std::overflow_error("prem: division by zero")); + if (is_ex_exactly_of_type(a, numeric)) { + if (is_ex_exactly_of_type(b, numeric)) + return _ex0(); + else + return b; + } + if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) + throw(std::invalid_argument("prem: arguments must be polynomials over the rationals")); + + // Polynomial long division + ex r = a.expand(); + ex eb = b.expand(); + int rdeg = r.degree(x); + int bdeg = eb.degree(x); + ex blcoeff; + if (bdeg <= rdeg) { + blcoeff = eb.coeff(x, bdeg); + if (bdeg == 0) + eb = _ex0(); + else + eb -= blcoeff * power(x, bdeg); + } else + blcoeff = _ex1(); + + while (rdeg >= bdeg && !r.is_zero()) { + ex rlcoeff = r.coeff(x, rdeg); + ex term = (power(x, rdeg - bdeg) * eb * rlcoeff).expand(); + if (rdeg == 0) + r = _ex0(); + else + r -= rlcoeff * power(x, rdeg); + r = (blcoeff * r).expand() - term; + rdeg = r.degree(x); + } + return r; +} + + /** Exact polynomial division of a(X) by b(X) in Q[X]. * * @param a first multivariate polynomial (dividend) @@ -471,12 +549,13 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) * coefficients (defaults to "true") * @return "true" when exact division succeeds (quotient returned in q), * "false" otherwise */ - 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; @@ -488,7 +567,8 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) return true; } #endif - if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) + if (check_args && (!a.info(info_flags::rational_polynomial) || + !b.info(info_flags::rational_polynomial))) throw(std::invalid_argument("divide: arguments must be polynomials over the rationals")); // Find first symbol @@ -598,38 +678,9 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite if (bdeg > adeg) return false; -#if 1 - - // Polynomial long division (recursive) - ex r = a.expand(); - if (r.is_zero()) - return true; - int rdeg = adeg; - ex eb = b.expand(); - ex blcoeff = eb.coeff(*x, bdeg); - while (rdeg >= bdeg) { - ex term, rcoeff = r.coeff(*x, rdeg); - if (!divide_in_z(rcoeff, blcoeff, term, var+1)) - break; - term = (term * power(*x, rdeg - bdeg)).expand(); - q += term; - r -= (term * eb).expand(); - if (r.is_zero()) { -#if USE_REMEMBER - dr_remember[ex2(a, b)] = exbool(q, true); -#endif - return true; - } - rdeg = r.degree(*x); - } -#if USE_REMEMBER - dr_remember[ex2(a, b)] = exbool(q, false); -#endif - return false; - -#else +#if USE_TRIAL_DIVISION - // Trial division using polynomial interpolation + // Trial division with polynomial interpolation int i, k; // Compute values at evaluation points 0..adeg @@ -652,7 +703,7 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite // Compute inverses vector rcp; rcp.reserve(adeg + 1); - rcp.push_back(0); + rcp.push_back(_num0()); for (k=1; k<=adeg; k++) { numeric product = alpha[k] - alpha[0]; for (i=1; i= bdeg) { + ex term, rcoeff = r.coeff(*x, rdeg); + if (!divide_in_z(rcoeff, blcoeff, term, var+1)) + break; + term = (term * power(*x, rdeg - bdeg)).expand(); + q += term; + r -= (term * eb).expand(); + if (r.is_zero()) { +#if USE_REMEMBER + dr_remember[ex2(a, b)] = exbool(q, true); +#endif + return true; + } + rdeg = r.degree(*x); + } +#if USE_REMEMBER + dr_remember[ex2(a, b)] = exbool(q, false); +#endif + return false; + #endif } @@ -778,7 +859,6 @@ ex ex::primpart(const symbol &x) const * @param x variable in which to compute the primitive part * @param c previously computed content part * @return primitive part */ - ex ex::primpart(const symbol &x, const ex &c) const { if (is_zero()) @@ -800,8 +880,85 @@ ex ex::primpart(const symbol &x, const ex &c) const * GCD of multivariate polynomials */ -/** Compute GCD of multivariate polynomials using the subresultant PRS - * algorithm. This function is used internally gy gcd(). +/** 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 @@ -809,8 +966,128 @@ ex ex::primpart(const symbol &x, const ex &c) const * @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 by gcd(). + * + * @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 sr_gcd(const ex &a, const ex &b, const symbol *x) { +//clog << "sr_gcd(" << a << "," << b << ")\n"; +#if STATISTICS + sr_gcd_called++; +#endif + // Sort c and d so that c has higher degree ex c, d; int adeg = a.degree(*x), bdeg = b.degree(*x); @@ -835,6 +1112,7 @@ static ex sr_gcd(const ex &a, const ex &b, const symbol *x) return gamma; c = c.primpart(*x, cont_c); d = d.primpart(*x, cont_d); +//clog << " content " << gamma << " removed, continuing with sr_gcd(" << c << "," << d << ")\n"; // First element of subresultant sequence ex r = _ex0(), ri = _ex1(), psi = _ex1(); @@ -842,12 +1120,15 @@ 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; - if (!divide(r, ri * power(psi, delta), d, false)) +//clog << " dividing...\n"; + 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) { @@ -858,11 +1139,12 @@ static ex sr_gcd(const ex &a, const ex &b, const symbol *x) } // Next element of subresultant sequence +//clog << " calculating next subresultant...\n"; ri = c.expand().lcoeff(*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; } } @@ -874,7 +1156,6 @@ static ex sr_gcd(const ex &a, const ex &b, const symbol *x) * @param e expanded multivariate polynomial * @return maximum coefficient * @see heur_gcd */ - numeric ex::max_coefficient(void) const { GINAC_ASSERT(bp!=0); @@ -930,7 +1211,6 @@ numeric mul::max_coefficient(void) const * @param xi modulus * @return mapped polynomial * @see heur_gcd */ - ex ex::smod(const numeric &xi) const { GINAC_ASSERT(bp!=0); @@ -1018,9 +1298,14 @@ class gcdheu_failed {}; * @return the GCD as a new expression * @see gcd * @exception gcdheu_failed() */ - static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const_iterator var) { +//clog << "heur_gcd(" << a << "," << b << ")\n"; +#if STATISTICS + heur_gcd_called++; +#endif + + // GCD of two numeric values -> CLN if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) { numeric g = gcd(ex_to_numeric(a), ex_to_numeric(b)); numeric rg; @@ -1053,8 +1338,10 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const // 6 tries maximum for (int t=0; t<6; t++) { - if (xi.int_length() * maxdeg > 50000) + if (xi.int_length() * maxdeg > 100000) { +//clog << "giving up heur_gcd, xi.int_length = " << xi.int_length() << ", maxdeg = " << maxdeg << endl; throw gcdheu_failed(); + } // Apply evaluation homomorphism and calculate GCD ex gamma = heur_gcd(p.subs(*x == xi), q.subs(*x == xi), NULL, NULL, var+1).expand(); @@ -1098,10 +1385,27 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const * @param check_args check whether a and b are polynomials with rational * coefficients (defaults to "true") * @return the GCD as a new expression */ - ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) { //clog << "gcd(" << a << "," << b << ")\n"; +#if STATISTICS + gcd_called++; +#endif + + // GCD of numerics -> CLN + if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) { + numeric g = gcd(ex_to_numeric(a), ex_to_numeric(b)); + if (ca) + *ca = ex_to_numeric(a) / g; + if (cb) + *cb = ex_to_numeric(b) / g; + return g; + } + + // Check arguments + if (check_args && !a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) { + throw(std::invalid_argument("gcd: arguments must be polynomials over the rationals")); + } // Partially factored cases (to avoid expanding large expressions) if (is_ex_exactly_of_type(a, mul)) { @@ -1142,6 +1446,51 @@ factored_b: return g; } +#if FAST_COMPARE + // Input polynomials of the form poly^n are sometimes also trivial + if (is_ex_exactly_of_type(a, power)) { + ex p = a.op(0); + if (is_ex_exactly_of_type(b, power)) { + if (p.is_equal(b.op(0))) { + // a = p^n, b = p^m, gcd = p^min(n, m) + ex exp_a = a.op(1), exp_b = b.op(1); + if (exp_a < exp_b) { + if (ca) + *ca = _ex1(); + if (cb) + *cb = power(p, exp_b - exp_a); + return power(p, exp_a); + } else { + if (ca) + *ca = power(p, exp_a - exp_b); + if (cb) + *cb = _ex1(); + return power(p, exp_b); + } + } + } else { + if (p.is_equal(b)) { + // a = p^n, b = p, gcd = p + if (ca) + *ca = power(p, a.op(1) - 1); + if (cb) + *cb = _ex1(); + return p; + } + } + } else if (is_ex_exactly_of_type(b, power)) { + ex p = b.op(0); + if (p.is_equal(a)) { + // a = p, b = p^n, gcd = p + if (ca) + *ca = _ex1(); + if (cb) + *cb = power(p, b.op(1) - 1); + return p; + } + } +#endif + // Some trivial cases ex aex = a.expand(), bex = b.expand(); if (aex.is_zero()) { @@ -1174,17 +1523,6 @@ factored_b: return a; } #endif - if (is_ex_exactly_of_type(aex, numeric) && is_ex_exactly_of_type(bex, numeric)) { - numeric g = gcd(ex_to_numeric(aex), ex_to_numeric(bex)); - if (ca) - *ca = ex_to_numeric(aex) / g; - if (cb) - *cb = ex_to_numeric(bex) / g; - return g; - } - if (check_args && !a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) { - throw(std::invalid_argument("gcd: arguments must be polynomials over the rationals")); - } // Gather symbol statistics sym_desc_vec sym_stats; @@ -1221,8 +1559,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) { @@ -1230,12 +1569,39 @@ factored_b: } if (is_ex_exactly_of_type(g, fail)) { //clog << "heuristics failed" << endl; - g = sr_gcd(aex, bex, x); - if (ca) - divide(aex, g, *ca, false); - if (cb) - divide(bex, g, *cb, false); - } +#if STATISTICS + heur_gcd_failed++; +#endif +#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) + *ca = a; + if (cb) + *cb = b; + } else { + if (ca) + divide(aex, g, *ca, false); + if (cb) + divide(bex, g, *cb, false); + } +#if 1 + } else { + if (g.is_equal(_ex1())) { + // Keep cofactors factored if possible + if (ca) + *ca = a; + if (cb) + *cb = b; + } + } +#endif return g; } @@ -1344,8 +1710,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) { @@ -1353,7 +1719,7 @@ static ex replace_with_symbol(const ex &e, lst &sym_lst, lst &repl_lst) for (unsigned i=0; isetflag(status_flags::dynallocated); + else if (level == -max_recursion_level) + throw(std::runtime_error("max recursion level reached")); + // Normalize and expand children, chop into summands exvector o; o.reserve(seq.size()+1); @@ -1538,6 +1929,11 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const * @see ex::normal() */ ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const { + if (level == 1) + return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated); + else if (level == -max_recursion_level) + throw(std::runtime_error("max recursion level reached")); + // Normalize children, separate into numerator and denominator ex num = _ex1(); ex den = _ex1(); @@ -1564,6 +1960,11 @@ ex mul::normal(lst &sym_lst, lst &repl_lst, int level) const * @see ex::normal */ ex power::normal(lst &sym_lst, lst &repl_lst, int level) const { + if (level == 1) + return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated); + else if (level == -max_recursion_level) + throw(std::runtime_error("max recursion level reached")); + // Normalize basis ex n = basis.bp->normal(sym_lst, repl_lst, level-1); @@ -1627,12 +2028,20 @@ ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const } +/** Implementation of ex::normal() for relationals. It normalizes both sides. + * @see ex::normal */ +ex relational::normal(lst &sym_lst, lst &repl_lst, int level) const +{ + return (new lst(relational(lh.normal(), rh.normal(), o), _ex1()))->setflag(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. @@ -1694,6 +2103,95 @@ 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_rational()) + return replace_with_symbol(*this, repl_lst); + } else { // complex + numeric re = real(); + numeric 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); +} + + +/** Implementation of ex::to_rational() for expairseqs. + * @see ex::to_rational */ +ex expairseq::to_rational(lst &repl_lst) const +{ + epvector s; + s.reserve(seq.size()); + for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) { + s.push_back(split_ex_to_pair(recombine_pair_to_ex(*it).to_rational(repl_lst))); + // s.push_back(combine_ex_with_coeff_to_pair((*it).rest.to_rational(repl_lst), + } + ex oc = overall_coeff.to_rational(repl_lst); + if (oc.info(info_flags::numeric)) + return thisexpairseq(s, overall_coeff); + else s.push_back(combine_ex_with_coeff_to_pair(oc,_ex1())); + return thisexpairseq(s, default_overall_coeff()); +} + + +/** 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