X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=1126d27c2fffb738e3c6cf57decf92af6dc430f0;hp=28ea9e6c1ff0d18f573094f58d38e6da00fe80d1;hb=a74473453218570d22f8932cc39ab48c7f0021ae;hpb=7a473d58cbf2ec744ef230e0503a071f9ab7aeec diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 28ea9e6c..1126d27c 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()) @@ -471,12 +498,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; @@ -598,38 +626,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 +651,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 +807,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 +828,9 @@ 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 @@ -809,8 +838,204 @@ ex ex::primpart(const symbol &x, const ex &c) const * @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 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 +1060,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 +1068,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 +1087,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 +1104,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 +1159,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 +1246,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 +1286,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(); @@ -1076,7 +1311,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const 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); - if (is_ex_exactly_of_type(lc, numeric) && lc.compare(_ex0()) < 0) + if (is_ex_exactly_of_type(lc, numeric) && ex_to_numeric(lc).is_negative()) return -g; else return g; @@ -1098,9 +1333,28 @@ 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)) { if (is_ex_exactly_of_type(b, mul) && b.nops() > a.nops()) @@ -1140,6 +1394,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()) { @@ -1172,17 +1471,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; @@ -1219,8 +1507,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) { @@ -1228,12 +1517,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; } @@ -1248,7 +1564,7 @@ factored_b: ex lcm(const ex &a, const ex &b, bool check_args) { if (is_ex_exactly_of_type(a, numeric) && is_ex_exactly_of_type(b, numeric)) - return gcd(ex_to_numeric(a), ex_to_numeric(b)); + return lcm(ex_to_numeric(a), ex_to_numeric(b)); if (check_args && !a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) throw(std::invalid_argument("lcm: arguments must be polynomials over the rationals")); @@ -1333,16 +1649,25 @@ ex sqrfree(const ex &a, const symbol &x) * Normal form of rational functions */ -// 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. +/* + * Note: The internal normal() functions (= basic::normal() and overloaded + * functions) all return lists of the form {numerator, denominator}. This + * is to get around mul::eval()'s automatic expansion of numeric coefficients. + * E.g. (a+b)/3 is automatically converted to a/3+b/3 but we want to keep + * the information that (a+b) is the numerator and 3 is the denominator. + */ + +/** Create a symbol for replacing the expression "e" (or return a previously + * 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) { // Expression already in repl_lst? Then return the assigned symbol for (unsigned i=0; isetflag(status_flags::dynallocated); } -/** Implementation of ex::normal() for symbols. This returns the unmodifies symbol. +/** Implementation of ex::normal() for symbols. This returns the unmodified symbol. * @see ex::normal */ ex symbol::normal(lst &sym_lst, lst &repl_lst, int level) const { - return *this; + return (new lst(*this, _ex1()))->setflag(status_flags::dynallocated); } @@ -1378,40 +1723,42 @@ ex symbol::normal(lst &sym_lst, lst &repl_lst, int level) const * @see ex::normal */ ex numeric::normal(lst &sym_lst, lst &repl_lst, int level) const { - if (is_real()) - if (is_rational()) - return *this; - else - return replace_with_symbol(*this, sym_lst, repl_lst); - else { // complex - numeric re = real(), im = imag(); + numeric num = numer(); + ex numex = num; + + if (num.is_real()) { + if (!num.is_integer()) + numex = replace_with_symbol(numex, sym_lst, repl_lst); + } else { // complex + numeric re = num.real(), im = num.imag(); ex re_ex = re.is_rational() ? re : replace_with_symbol(re, sym_lst, repl_lst); ex im_ex = im.is_rational() ? im : replace_with_symbol(im, sym_lst, repl_lst); - return re_ex + im_ex * replace_with_symbol(I, sym_lst, repl_lst); + numex = re_ex + im_ex * replace_with_symbol(I, sym_lst, repl_lst); } + + // Denominator is always a real integer (see numeric::denom()) + return (new lst(numex, denom()))->setflag(status_flags::dynallocated); } /** Fraction cancellation. * @param n numerator * @param d denominator - * @return cancelled fraction n/d */ + * @return cancelled fraction {n, d} as a list */ static ex frac_cancel(const ex &n, const ex &d) { ex num = n; ex den = d; numeric pre_factor = _num1(); +//clog << "frac_cancel num = " << num << ", den = " << den << endl; + // Handle special cases where numerator or denominator is 0 if (num.is_zero()) - return _ex0(); + return (new lst(_ex0(), _ex1()))->setflag(status_flags::dynallocated); if (den.expand().is_zero()) throw(std::overflow_error("frac_cancel: division by zero in frac_cancel")); - // More special cases - if (is_ex_exactly_of_type(den, numeric)) - return num / den; - // Bring numerator and denominator to Z[X] by multiplying with // LCM of all coefficients' denominators numeric num_lcm = lcm_of_coefficients_denominators(num); @@ -1431,12 +1778,16 @@ static ex frac_cancel(const ex &n, const ex &d) // as defined by get_first_symbol() is made positive) const symbol *x; if (get_first_symbol(den, x)) { - if (den.unit(*x).compare(_ex0()) < 0) { + GINAC_ASSERT(is_ex_exactly_of_type(den.unit(*x),numeric)); + if (ex_to_numeric(den.unit(*x)).is_negative()) { num *= _ex_1(); den *= _ex_1(); } } - return pre_factor * num / den; + + // Return result as list +//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); } @@ -1445,47 +1796,70 @@ 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 { - // Normalize and expand children + // Normalize and expand children, chop into summands exvector o; o.reserve(seq.size()+1); epvector::const_iterator it = seq.begin(), itend = seq.end(); while (it != itend) { + + // Normalize and expand child ex n = recombine_pair_to_ex(*it).bp->normal(sym_lst, repl_lst, level-1).expand(); - if (is_ex_exactly_of_type(n, add)) { - epvector::const_iterator bit = (static_cast(n.bp))->seq.begin(), bitend = (static_cast(n.bp))->seq.end(); + + // If numerator is a sum, chop into summands + if (is_ex_exactly_of_type(n.op(0), add)) { + epvector::const_iterator bit = ex_to_add(n.op(0)).seq.begin(), bitend = ex_to_add(n.op(0)).seq.end(); while (bit != bitend) { - o.push_back(recombine_pair_to_ex(*bit)); + o.push_back((new lst(recombine_pair_to_ex(*bit), n.op(1)))->setflag(status_flags::dynallocated)); bit++; } - o.push_back((static_cast(n.bp))->overall_coeff); + + // The overall_coeff is already normalized (== rational), we just + // split it into numerator and denominator + GINAC_ASSERT(ex_to_numeric(ex_to_add(n.op(0)).overall_coeff).is_rational()); + numeric overall = ex_to_numeric(ex_to_add(n.op(0)).overall_coeff); + o.push_back((new lst(overall.numer(), overall.denom() * n.op(1)))->setflag(status_flags::dynallocated)); } else o.push_back(n); it++; } o.push_back(overall_coeff.bp->normal(sym_lst, repl_lst, level-1)); + // o is now a vector of {numerator, denominator} lists + // Determine common denominator ex den = _ex1(); exvector::const_iterator ait = o.begin(), aitend = o.end(); +//clog << "add::normal uses the following summands:\n"; while (ait != aitend) { - den = lcm((*ait).denom(false), den, false); +//clog << " num = " << ait->op(0) << ", den = " << ait->op(1) << endl; + den = lcm(ait->op(1), den, false); ait++; } +//clog << " common denominator = " << den << endl; // Add fractions - if (den.is_equal(_ex1())) - return (new add(o))->setflag(status_flags::dynallocated); - else { + if (den.is_equal(_ex1())) { + + // Common denominator is 1, simply add all numerators + exvector num_seq; + for (ait=o.begin(); ait!=aitend; ait++) { + num_seq.push_back(ait->op(0)); + } + return (new lst((new add(num_seq))->setflag(status_flags::dynallocated), den))->setflag(status_flags::dynallocated); + + } else { + + // Perform fractional addition exvector num_seq; for (ait=o.begin(); ait!=aitend; ait++) { ex q; - if (!divide(den, (*ait).denom(false), q, false)) { + if (!divide(den, ait->op(1), q, false)) { // should not happen throw(std::runtime_error("invalid expression in add::normal, division failed")); } - num_seq.push_back((*ait).numer(false) * q); + num_seq.push_back((ait->op(0) * q).expand()); } - ex num = add(num_seq); + ex num = (new add(num_seq))->setflag(status_flags::dynallocated); // Cancel common factors from num/den return frac_cancel(num, den); @@ -1498,17 +1872,23 @@ 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 { - // Normalize children - exvector o; - o.reserve(seq.size()+1); + // Normalize children, separate into numerator and denominator + ex num = _ex1(); + ex den = _ex1(); + ex n; epvector::const_iterator it = seq.begin(), itend = seq.end(); while (it != itend) { - o.push_back(recombine_pair_to_ex(*it).bp->normal(sym_lst, repl_lst, level-1)); + n = recombine_pair_to_ex(*it).bp->normal(sym_lst, repl_lst, level-1); + num *= n.op(0); + den *= n.op(1); it++; } - o.push_back(overall_coeff.bp->normal(sym_lst, repl_lst, level-1)); - ex n = (new mul(o))->setflag(status_flags::dynallocated); - return frac_cancel(n.numer(false), n.denom(false)); + n = overall_coeff.bp->normal(sym_lst, repl_lst, level-1); + num *= n.op(0); + den *= n.op(1); + + // Perform fraction cancellation + return frac_cancel(num, den); } @@ -1518,16 +1898,47 @@ 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 (exponent.info(info_flags::integer)) { - // Integer powers are distributed - ex n = basis.bp->normal(sym_lst, repl_lst, level-1); - ex num = n.numer(false); - ex den = n.denom(false); - return power(num, exponent) / power(den, exponent); - } else { - // Non-integer powers are replaced by temporary symbol (after normalizing basis) - ex n = power(basis.bp->normal(sym_lst, repl_lst, level-1), exponent); - return replace_with_symbol(n, sym_lst, repl_lst); + // Normalize basis + ex n = basis.bp->normal(sym_lst, repl_lst, level-1); + + if (exponent.info(info_flags::integer)) { + + if (exponent.info(info_flags::positive)) { + + // (a/b)^n -> {a^n, b^n} + return (new lst(power(n.op(0), exponent), power(n.op(1), exponent)))->setflag(status_flags::dynallocated); + + } else if (exponent.info(info_flags::negative)) { + + // (a/b)^-n -> {b^n, a^n} + return (new lst(power(n.op(1), -exponent), power(n.op(0), -exponent)))->setflag(status_flags::dynallocated); + } + + } else { + + if (exponent.info(info_flags::positive)) { + + // (a/b)^x -> {sym((a/b)^x), 1} + return (new lst(replace_with_symbol(power(n.op(0) / n.op(1), exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); + + } else if (exponent.info(info_flags::negative)) { + + if (n.op(1).is_equal(_ex1())) { + + // a^-x -> {1, sym(a^x)} + return (new lst(_ex1(), replace_with_symbol(power(n.op(0), -exponent), sym_lst, repl_lst)))->setflag(status_flags::dynallocated); + + } else { + + // (a/b)^-x -> {sym((b/a)^x), 1} + return (new lst(replace_with_symbol(power(n.op(1) / n.op(0), -exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); + } + + } else { // exponent not numeric + + // (a/b)^x -> {sym((a/b)^x, 1} + return (new lst(replace_with_symbol(power(n.op(0) / n.op(1), exponent), sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); + } } } @@ -1545,9 +1956,16 @@ ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const new_seq.push_back(expair(it->rest.normal(), it->coeff)); it++; } + ex n = pseries(relational(var,point), new_seq); + return (new lst(replace_with_symbol(n, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated); +} - ex n = pseries(var, point, new_seq); - return replace_with_symbol(n, sym_lst, repl_lst); + +/** 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); } @@ -1555,8 +1973,8 @@ ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const * 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. @@ -1566,13 +1984,146 @@ ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const ex ex::normal(int level) const { lst sym_lst, repl_lst; + ex e = bp->normal(sym_lst, repl_lst, level); + GINAC_ASSERT(is_ex_of_type(e, lst)); + + // Re-insert replaced symbols if (sym_lst.nops() > 0) - return e.subs(sym_lst, repl_lst); - else - return e; + e = e.subs(sym_lst, repl_lst); + + // Convert {numerator, denominator} form back to fraction + return e.op(0) / e.op(1); } +/** Numerator of an expression. If the expression is not of the normal form + * "numerator/denominator", it is first converted to this form and then the + * numerator is returned. + * + * @see ex::normal + * @return numerator */ +ex ex::numer(void) const +{ + lst sym_lst, repl_lst; + + ex e = bp->normal(sym_lst, repl_lst, 0); + GINAC_ASSERT(is_ex_of_type(e, lst)); + + // Re-insert replaced symbols + if (sym_lst.nops() > 0) + return e.op(0).subs(sym_lst, repl_lst); + else + return e.op(0); +} + +/** Denominator of an expression. If the expression is not of the normal form + * "numerator/denominator", it is first converted to this form and then the + * denominator is returned. + * + * @see ex::normal + * @return denominator */ +ex ex::denom(void) const +{ + lst sym_lst, repl_lst; + + ex e = bp->normal(sym_lst, repl_lst, 0); + GINAC_ASSERT(is_ex_of_type(e, lst)); + + // Re-insert replaced symbols + if (sym_lst.nops() > 0) + return e.op(1).subs(sym_lst, repl_lst); + else + 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(); + 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