X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=514423c97930f4eff87ad74205553735069f3b26;hp=9e24b99bc12e68bb2a5cdcbe263d4637bd9c64ea;hb=55d35dcf72dc411c8265628fcad2bd67d320a8c9;hpb=fc80c36140aaf126150b2dea811177d8af35dac3 diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 9e24b99b..514423c9 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -59,7 +59,8 @@ namespace GiNaC { #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) +// polynomial interpolation (always slower except for completely dense +// polynomials) #define USE_TRIAL_DIVISION 0 // Set this to enable some statistical output for the GCD routines @@ -93,7 +94,6 @@ static struct _stat_print { * @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)) { @@ -145,7 +145,7 @@ struct sym_desc { }; // Vector of sym_desc structures -typedef vector sym_desc_vec; +typedef std::vector sym_desc_vec; // Add symbol the sym_desc_vec (used internally by get_symbol_stats()) static void add_symbol(const symbol *s, sym_desc_vec &v) @@ -186,7 +186,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 @@ -204,11 +203,11 @@ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v) } sort(v.begin(), v.end()); #if 0 - clog << "Symbols:\n"; + std::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; + std::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; + std::clog << " lcoeff_a=" << a.lcoeff(*(it->sym)) << ", lcoeff_b=" << b.lcoeff(*(it->sym)) << endl; it++; } #endif @@ -247,7 +246,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()); @@ -258,7 +256,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)) { @@ -288,7 +285,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); @@ -349,7 +345,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()) @@ -400,7 +395,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()) @@ -452,7 +446,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()) @@ -497,6 +490,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) @@ -506,14 +550,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 (a.is_zero()) + return true; if (is_ex_exactly_of_type(b, numeric)) { q = a / b; return true; @@ -525,7 +568,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 @@ -564,8 +608,8 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) * Remembering */ -typedef pair ex2; -typedef pair exbool; +typedef std::pair ex2; +typedef std::pair exbool; struct ex2_less { bool operator() (const ex2 p, const ex2 q) const @@ -574,7 +618,7 @@ struct ex2_less { } }; -typedef map ex2_exbool_remember; +typedef std::map ex2_exbool_remember; #endif @@ -816,7 +860,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()) @@ -838,9 +881,9 @@ 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. +/** Compute GCD of polynomials in Q[X] using the Euclidean 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 @@ -850,7 +893,7 @@ ex ex::primpart(const symbol &x, const ex &c) const static ex eu_gcd(const ex &a, const ex &b, const symbol *x) { -//clog << "eu_gcd(" << a << "," << b << ")\n"; +//std::clog << "eu_gcd(" << a << "," << b << ")\n"; // Sort c and d so that c has higher degree ex c, d; @@ -863,13 +906,17 @@ static ex eu_gcd(const ex &a, const ex &b, const symbol *x) d = a; } + // Normalize in Q[x] + c = c / c.lcoeff(*x); + d = d / d.lcoeff(*x); + // Euclidean algorithm ex r; for (;;) { -//clog << " d = " << d << endl; +//std::clog << " d = " << d << endl; r = rem(c, d, *x, false); if (r.is_zero()) - return d.primpart(*x); + return d / d.lcoeff(*x); c = d; d = r; } @@ -888,7 +935,7 @@ static ex eu_gcd(const ex &a, const ex &b, const symbol *x) static ex euprem_gcd(const ex &a, const ex &b, const symbol *x) { -//clog << "euprem_gcd(" << a << "," << b << ")\n"; +//std::clog << "euprem_gcd(" << a << "," << b << ")\n"; // Sort c and d so that c has higher degree ex c, d; @@ -901,13 +948,16 @@ static ex euprem_gcd(const ex &a, const ex &b, const symbol *x) d = a; } + // Calculate GCD of contents + ex gamma = gcd(c.content(*x), d.content(*x), NULL, NULL, false); + // Euclidean algorithm with pseudo-remainders ex r; for (;;) { -//clog << " d = " << d << endl; +//std::clog << " d = " << d << endl; r = prem(c, d, *x, false); if (r.is_zero()) - return d.primpart(*x); + return d.primpart(*x) * gamma; c = d; d = r; } @@ -926,7 +976,7 @@ static ex euprem_gcd(const ex &a, const ex &b, const symbol *x) static ex peu_gcd(const ex &a, const ex &b, const symbol *x) { -//clog << "peu_gcd(" << a << "," << b << ")\n"; +//std::clog << "peu_gcd(" << a << "," << b << ")\n"; // Sort c and d so that c has higher degree ex c, d; @@ -954,7 +1004,7 @@ static ex peu_gcd(const ex &a, const ex &b, const symbol *x) // Euclidean algorithm with content removal ex r; for (;;) { -//clog << " d = " << d << endl; +//std::clog << " d = " << d << endl; r = prem(c, d, *x, false); if (r.is_zero()) return gamma * d; @@ -975,7 +1025,7 @@ static ex peu_gcd(const ex &a, const ex &b, const symbol *x) static ex red_gcd(const ex &a, const ex &b, const symbol *x) { -//clog << "red_gcd(" << a << "," << b << ")\n"; +//std::clog << "red_gcd(" << a << "," << b << ")\n"; // Sort c and d so that c has higher degree ex c, d; @@ -1002,13 +1052,13 @@ static ex red_gcd(const ex &a, const ex &b, const symbol *x) c = c.primpart(*x, cont_c); d = d.primpart(*x, cont_d); - // First element of subresultant sequence + // First element of divisor sequence ex r, ri = _ex1(); int delta = cdeg - ddeg; for (;;) { // Calculate polynomial pseudo-remainder -//clog << " d = " << d << endl; +//std::clog << " d = " << d << endl; r = prem(c, d, *x, false); if (r.is_zero()) return gamma * d.primpart(*x); @@ -1034,22 +1084,25 @@ static ex red_gcd(const ex &a, const ex &b, const symbol *x) /** 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 + * @param a first multivariate polynomial + * @param b second multivariate polynomial + * @param var iterator to first element of vector of sym_desc structs * @return the GCD as a new expression * @see gcd */ -static ex sr_gcd(const ex &a, const ex &b, const symbol *x) +static ex sr_gcd(const ex &a, const ex &b, sym_desc_vec::const_iterator var) { -//clog << "sr_gcd(" << a << "," << b << ")\n"; +//std::clog << "sr_gcd(" << a << "," << b << ")\n"; #if STATISTICS sr_gcd_called++; #endif + // The first symbol is our main variable + const symbol &x = *(var->sym); + // Sort c and d so that c has higher degree ex c, d; - int adeg = a.degree(*x), bdeg = b.degree(*x); + int adeg = a.degree(x), bdeg = b.degree(x); int cdeg, ddeg; if (adeg >= bdeg) { c = a; @@ -1064,14 +1117,14 @@ static ex sr_gcd(const ex &a, const ex &b, const symbol *x) } // 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 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); -//clog << " content " << gamma << " removed, continuing with sr_gcd(" << c << "," << d << ")\n"; + c = c.primpart(x, cont_c); + d = d.primpart(x, cont_d); +//std::clog << " content " << gamma << " removed, continuing with sr_gcd(" << c << "," << d << ")\n"; // First element of subresultant sequence ex r = _ex0(), ri = _ex1(), psi = _ex1(); @@ -1079,31 +1132,31 @@ 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); +//std::clog << " start of loop, psi = " << psi << ", calculating pseudo-remainder...\n"; +//std::clog << " d = " << d << endl; + r = prem(c, d, x, false); if (r.is_zero()) - return gamma * d.primpart(*x); + return gamma * d.primpart(x); c = d; cdeg = ddeg; -//clog << " dividing...\n"; - if (!divide(r, ri * pow(psi, delta), d, false)) +//std::clog << " dividing...\n"; + if (!divide_in_z(r, ri * pow(psi, delta), d, var)) throw(std::runtime_error("invalid expression in sr_gcd(), division failed")); - ddeg = d.degree(*x); + ddeg = d.degree(x); if (ddeg == 0) { if (is_ex_exactly_of_type(r, numeric)) return gamma; else - return gamma * r.primpart(*x); + return gamma * r.primpart(x); } // Next element of subresultant sequence -//clog << " calculating next subresultant...\n"; - ri = c.expand().lcoeff(*x); +//std::clog << " calculating next subresultant...\n"; + ri = c.expand().lcoeff(x); if (delta == 1) psi = ri; else if (delta) - divide(pow(ri, delta), pow(psi, delta-1), psi, false); + divide_in_z(pow(ri, delta), pow(psi, delta-1), psi, var+1); delta = cdeg - ddeg; } } @@ -1115,7 +1168,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); @@ -1171,7 +1223,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); @@ -1241,6 +1292,20 @@ ex mul::smod(const numeric &xi) const } +/** xi-adic polynomial interpolation */ +static ex interpolate(const ex &gamma, const numeric &xi, const symbol &x) +{ + ex g = _ex0(); + ex e = gamma; + numeric rxi = xi.inverse(); + for (int i=0; !e.is_zero(); i++) { + ex gi = e.smod(xi); + g += gi * power(x, i); + e = (e - gi) * rxi; + } + return g; +} + /** Exception thrown by heur_gcd() to signal failure. */ class gcdheu_failed {}; @@ -1259,36 +1324,36 @@ 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"; +//std::clog << "heur_gcd(" << a << "," << b << ")\n"; #if STATISTICS heur_gcd_called++; #endif + // Algorithms only works for non-vanishing input polynomials + if (a.is_zero() || b.is_zero()) + return *new ex(fail()); + // 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; - if (ca || cb) - rg = g.inverse(); if (ca) - *ca = ex_to_numeric(a).mul(rg); + *ca = ex_to_numeric(a) / g; if (cb) - *cb = ex_to_numeric(b).mul(rg); + *cb = ex_to_numeric(b) / g; return g; } // The first symbol is our main variable - const symbol *x = var->sym; + const symbol &x = *(var->sym); // Remove integer content numeric gc = gcd(a.integer_content(), b.integer_content()); numeric rgc = gc.inverse(); ex p = a * rgc; ex q = b * rgc; - int maxdeg = max(p.degree(*x), q.degree(*x)); + int maxdeg = max(p.degree(x), q.degree(x)); // Find evaluation point numeric mp = p.max_coefficient(), mq = q.max_coefficient(); @@ -1301,35 +1366,59 @@ 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 > 100000) { -//clog << "giving up heur_gcd, xi.int_length = " << xi.int_length() << ", maxdeg = " << maxdeg << endl; +//std::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(); + ex cp, cq; + ex gamma = heur_gcd(p.subs(x == xi), q.subs(x == xi), &cp, &cq, var+1).expand(); if (!is_ex_exactly_of_type(gamma, fail)) { // Reconstruct polynomial from GCD of mapped polynomials - ex g = _ex0(); - numeric rxi = xi.inverse(); - for (int i=0; !gamma.is_zero(); i++) { - ex gi = gamma.smod(xi); - g += gi * power(*x, i); - gamma = (gamma - gi) * rxi; - } + ex g = interpolate(gamma, xi, x); + // Remove integer content g /= g.integer_content(); - // If the calculated polynomial divides both a and b, this is the GCD + // If the calculated polynomial divides both p and q, this is the GCD ex dummy; if (divide_in_z(p, g, ca ? *ca : dummy, var) && divide_in_z(q, g, cb ? *cb : dummy, var)) { g *= gc; - ex lc = g.lcoeff(*x); + ex lc = g.lcoeff(x); if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) return -g; else return g; } +#if 0 + cp = interpolate(cp, xi, x); + if (divide_in_z(cp, p, g, var)) { + if (divide_in_z(g, q, cb ? *cb : dummy, var)) { + g *= gc; + if (ca) + *ca = cp; + ex lc = g.lcoeff(x); + if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) + return -g; + else + return g; + } + } + cq = interpolate(cq, xi, x); + if (divide_in_z(cq, q, g, var)) { + if (divide_in_z(g, p, ca ? *ca : dummy, var)) { + g *= gc; + if (cb) + *cb = cq; + ex lc = g.lcoeff(x); + if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) + return -g; + else + return g; + } + } +#endif } // Next evaluation point @@ -1347,10 +1436,9 @@ 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"; +//std::clog << "gcd(" << a << "," << b << ")\n"; #if STATISTICS gcd_called++; #endif @@ -1358,10 +1446,19 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) // 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; + if (ca || cb) { + if (g.is_zero()) { + if (ca) + *ca = _ex0(); + if (cb) + *cb = _ex0(); + } else { + if (ca) + *ca = ex_to_numeric(a) / g; + if (cb) + *cb = ex_to_numeric(b) / g; + } + } return g; } @@ -1493,32 +1590,32 @@ factored_b: // The symbol with least degree is our main variable sym_desc_vec::const_iterator var = sym_stats.begin(); - const symbol *x = var->sym; + const symbol &x = *(var->sym); // Cancel trivial common factor int ldeg_a = var->ldeg_a; int ldeg_b = var->ldeg_b; int min_ldeg = min(ldeg_a, ldeg_b); if (min_ldeg > 0) { - ex common = power(*x, min_ldeg); -//clog << "trivial common factor " << common << endl; + ex common = power(x, min_ldeg); +//std::clog << "trivial common factor " << common << endl; return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common; } // Try to eliminate variables if (var->deg_a == 0) { -//clog << "eliminating variable " << *x << " from b" << endl; - ex c = bex.content(*x); +//std::clog << "eliminating variable " << x << " from b" << endl; + ex c = bex.content(x); ex g = gcd(aex, c, ca, cb, false); if (cb) - *cb *= bex.unit(*x) * bex.primpart(*x, c); + *cb *= bex.unit(x) * bex.primpart(x, c); return g; } else if (var->deg_b == 0) { -//clog << "eliminating variable " << *x << " from a" << endl; - ex c = aex.content(*x); +//std::clog << "eliminating variable " << x << " from a" << endl; + ex c = aex.content(x); ex g = gcd(c, bex, ca, cb, false); if (ca) - *ca *= aex.unit(*x) * aex.primpart(*x, c); + *ca *= aex.unit(x) * aex.primpart(x, c); return g; } @@ -1531,17 +1628,17 @@ factored_b: g = *new ex(fail()); } if (is_ex_exactly_of_type(g, fail)) { -//clog << "heuristics failed" << endl; +//std::clog << "heuristics failed" << endl; #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); +// 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, var); if (g.is_equal(_ex1())) { // Keep cofactors factored if possible if (ca) @@ -1682,7 +1779,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 basis ex n = basis.bp->normal(sym_lst, repl_lst, level-1); @@ -2061,7 +2173,8 @@ ex basic::to_rational(lst &repl_lst) const } -/** Implementation of ex::to_rational() for symbols. This returns the unmodified symbol. +/** Implementation of ex::to_rational() for symbols. This returns the + * unmodified symbol. * @see ex::to_rational */ ex symbol::to_rational(lst &repl_lst) const { @@ -2069,17 +2182,18 @@ ex symbol::to_rational(lst &repl_lst) const } -/** 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. +/** 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()) + if (!is_rational()) return replace_with_symbol(*this, repl_lst); } else { // complex - numeric re = real(), im = imag(); + 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); @@ -2100,6 +2214,24 @@ ex power::to_rational(lst &repl_lst) const } +/** 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,