X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=514423c97930f4eff87ad74205553735069f3b26;hp=4e814361ed0b99d98ec3617dafa9746babf5eea1;hb=55d35dcf72dc411c8265628fcad2bd67d320a8c9;hpb=61e847cdfd8455a7ba67a0038add8dd35b4e0bff diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 4e814361..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 @@ -144,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) @@ -202,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 @@ -489,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) @@ -516,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 @@ -555,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 @@ -565,7 +618,7 @@ struct ex2_less { } }; -typedef map ex2_exbool_remember; +typedef std::map ex2_exbool_remember; #endif @@ -828,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 @@ -840,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; @@ -853,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; } @@ -878,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; @@ -891,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; } @@ -916,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; @@ -944,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; @@ -965,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; @@ -992,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); @@ -1024,21 +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; @@ -1053,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(); @@ -1068,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; } } @@ -1228,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 {}; @@ -1248,33 +1326,34 @@ class gcdheu_failed {}; * @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(); @@ -1287,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 @@ -1335,7 +1438,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const * @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 @@ -1343,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; } @@ -1478,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; } @@ -1516,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) @@ -1751,7 +1863,7 @@ static ex frac_cancel(const ex &n, const ex &d) ex den = d; numeric pre_factor = _num1(); -//clog << "frac_cancel num = " << num << ", den = " << den << endl; +//std::clog << "frac_cancel num = " << num << ", den = " << den << endl; // Handle special cases where numerator or denominator is 0 if (num.is_zero()) @@ -1786,7 +1898,7 @@ static ex frac_cancel(const ex &n, const ex &d) } // Return result as list -//clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << endl; +//std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << endl; return (new lst(num * pre_factor.numer(), den * pre_factor.denom()))->setflag(status_flags::dynallocated); } @@ -1796,6 +1908,11 @@ static ex frac_cancel(const ex &n, const ex &d) * @see ex::normal */ ex add::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 and expand children, chop into summands exvector o; o.reserve(seq.size()+1); @@ -1829,21 +1946,21 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const // Determine common denominator ex den = _ex1(); exvector::const_iterator ait = o.begin(), aitend = o.end(); -//clog << "add::normal uses the following summands:\n"; +//std::clog << "add::normal uses the following summands:\n"; while (ait != aitend) { -//clog << " num = " << ait->op(0) << ", den = " << ait->op(1) << endl; +//std::clog << " num = " << ait->op(0) << ", den = " << ait->op(1) << endl; den = lcm(ait->op(1), den, false); ait++; } -//clog << " common denominator = " << den << endl; +//std::clog << " common denominator = " << den << endl; // Add fractions if (den.is_equal(_ex1())) { - // Common denominator is 1, simply add all numerators + // Common denominator is 1, simply add all fractions exvector num_seq; for (ait=o.begin(); ait!=aitend; ait++) { - num_seq.push_back(ait->op(0)); + num_seq.push_back(ait->op(0) / ait->op(1)); } return (new lst((new add(num_seq))->setflag(status_flags::dynallocated), den))->setflag(status_flags::dynallocated); @@ -1872,6 +1989,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(); @@ -1898,6 +2020,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);