X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=f57e4f0dbb60b31c3fc7bb053b905881ffee24f2;hp=391b80b6dd88c6c77615bda512427b718d175cc1;hb=14405559f3674c95cae22c5238730c375bc965e4;hpb=24fe247f9ed16114a765a01c593fec5c4a2f591c diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 391b80b6..f57e4f0d 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -7,7 +7,7 @@ */ /* - * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -43,10 +43,13 @@ #include "numeric.h" #include "power.h" #include "relational.h" -#include "series.h" +#include "pseries.h" #include "symbol.h" +#include "utils.h" +#ifndef NO_GINAC_NAMESPACE namespace GiNaC { +#endif // ndef NO_GINAC_NAMESPACE // If comparing expressions (ex::compare()) is fast, you can set this to 1. // Some routines like quo(), rem() and gcd() will then return a quick answer @@ -71,7 +74,7 @@ static bool get_first_symbol(const ex &e, const symbol *&x) x = static_cast(e.bp); return true; } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) { - for (int i=0; i(e.bp), v); } else if (is_ex_exactly_of_type(e, add) || is_ex_exactly_of_type(e, mul)) { - for (int i=0; irest,numeric)); GINAC_ASSERT(is_ex_exactly_of_type(it->coeff,numeric)); @@ -287,13 +318,13 @@ ex quo(const ex &a, const ex &b, const symbol &x, bool check_args) return a / b; #if FAST_COMPARE if (a.is_equal(b)) - return exONE(); + return _ex1(); #endif if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) throw(std::invalid_argument("quo: arguments must be polynomials over the rationals")); // Polynomial long division - ex q = exZERO(); + ex q = _ex0(); ex r = a.expand(); if (r.is_zero()) return r; @@ -336,13 +367,13 @@ ex rem(const ex &a, const ex &b, const symbol &x, bool check_args) throw(std::overflow_error("rem: division by zero")); if (is_ex_exactly_of_type(a, numeric)) { if (is_ex_exactly_of_type(b, numeric)) - return exZERO(); + return _ex0(); else return b; } #if FAST_COMPARE if (a.is_equal(b)) - return exZERO(); + return _ex0(); #endif if (check_args && (!a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial))) throw(std::invalid_argument("rem: arguments must be polynomials over the rationals")); @@ -388,7 +419,7 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) 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 exZERO(); + return _ex0(); else return b; } @@ -404,18 +435,18 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) if (bdeg <= rdeg) { blcoeff = eb.coeff(x, bdeg); if (bdeg == 0) - eb = exZERO(); + eb = _ex0(); else eb -= blcoeff * power(x, bdeg); } else - blcoeff = exONE(); + blcoeff = _ex1(); int delta = rdeg - bdeg + 1, i = 0; 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 = exZERO(); + r = _ex0(); else r -= rlcoeff * power(x, rdeg); r = (blcoeff * r).expand() - term; @@ -438,7 +469,7 @@ ex prem(const ex &a, const ex &b, const symbol &x, bool check_args) bool divide(const ex &a, const ex &b, ex &q, bool check_args) { - q = exZERO(); + q = _ex0(); if (b.is_zero()) throw(std::overflow_error("divide: division by zero")); if (is_ex_exactly_of_type(b, numeric)) { @@ -448,7 +479,7 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) return false; #if FAST_COMPARE if (a.is_equal(b)) { - q = exONE(); + q = _ex1(); return true; } #endif @@ -523,10 +554,10 @@ typedef map ex2_exbool_remember; * @see get_symbol_stats, heur_gcd */ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_iterator var) { - q = exZERO(); + q = _ex0(); if (b.is_zero()) throw(std::overflow_error("divide_in_z: division by zero")); - if (b.is_equal(exONE())) { + if (b.is_equal(_ex1())) { q = a; return true; } @@ -539,7 +570,7 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite } #if FAST_COMPARE if (a.is_equal(b)) { - q = exONE(); + q = _ex1(); return true; } #endif @@ -599,19 +630,19 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite // Compute values at evaluation points 0..adeg vector alpha; alpha.reserve(adeg + 1); exvector u; u.reserve(adeg + 1); - numeric point = numZERO(); + numeric point = _num0(); ex c; for (i=0; i<=adeg; i++) { ex bs = b.subs(*x == point); while (bs.is_zero()) { - point += numONE(); + point += _num1(); bs = b.subs(*x == point); } if (!divide_in_z(a.subs(*x == point), bs, c, var+1)) return false; alpha.push_back(point); u.push_back(c); - point += numONE(); + point += _num1(); } // Compute inverses @@ -663,7 +694,7 @@ ex ex::unit(const symbol &x) const { ex c = expand().lcoeff(x); if (is_ex_exactly_of_type(c, numeric)) - return c < exZERO() ? exMINUSONE() : exONE(); + return c < _ex0() ? _ex_1() : _ex1(); else { const symbol *y; if (get_first_symbol(c, y)) @@ -684,12 +715,12 @@ ex ex::unit(const symbol &x) const ex ex::content(const symbol &x) const { if (is_zero()) - return exZERO(); + return _ex0(); if (is_ex_exactly_of_type(*this, numeric)) return info(info_flags::negative) ? -*this : *this; ex e = expand(); if (e.is_zero()) - return exZERO(); + return _ex0(); // First, try the integer content ex c = e.integer_content(); @@ -703,7 +734,7 @@ ex ex::content(const symbol &x) const int ldeg = e.ldegree(x); if (deg == ldeg) return e.lcoeff(x) / e.unit(x); - c = exZERO(); + c = _ex0(); for (int i=ldeg; i<=deg; i++) c = gcd(e.coeff(x, i), c, NULL, NULL, false); return c; @@ -720,13 +751,13 @@ ex ex::content(const symbol &x) const ex ex::primpart(const symbol &x) const { if (is_zero()) - return exZERO(); + return _ex0(); if (is_ex_exactly_of_type(*this, numeric)) - return exONE(); + return _ex1(); ex c = content(x); if (c.is_zero()) - return exZERO(); + return _ex0(); ex u = unit(x); if (is_ex_exactly_of_type(c, numeric)) return *this / (c * u); @@ -746,11 +777,11 @@ ex ex::primpart(const symbol &x) const ex ex::primpart(const symbol &x, const ex &c) const { if (is_zero()) - return exZERO(); + return _ex0(); if (c.is_zero()) - return exZERO(); + return _ex0(); if (is_ex_exactly_of_type(*this, numeric)) - return exONE(); + return _ex1(); ex u = unit(x); if (is_ex_exactly_of_type(c, numeric)) @@ -801,7 +832,7 @@ static ex sr_gcd(const ex &a, const ex &b, const symbol *x) d = d.primpart(*x, cont_d); // First element of subresultant sequence - ex r = exZERO(), ri = exONE(), psi = exONE(); + ex r = _ex0(), ri = _ex1(), psi = _ex1(); int delta = cdeg - ddeg; for (;;) { @@ -847,7 +878,7 @@ numeric ex::max_coefficient(void) const numeric basic::max_coefficient(void) const { - return numONE(); + return _num1(); } numeric numeric::max_coefficient(void) const @@ -908,7 +939,11 @@ ex basic::smod(const numeric &xi) const ex numeric::smod(const numeric &xi) const { +#ifndef NO_GINAC_NAMESPACE return GiNaC::smod(*this, xi); +#else // ndef NO_GINAC_NAMESPACE + return ::smod(*this, xi); +#endif // ndef NO_GINAC_NAMESPACE } ex add::smod(const numeric &xi) const @@ -919,13 +954,21 @@ ex add::smod(const numeric &xi) const epvector::const_iterator itend = seq.end(); while (it != itend) { GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); +#ifndef NO_GINAC_NAMESPACE numeric coeff = GiNaC::smod(ex_to_numeric(it->coeff), xi); +#else // ndef NO_GINAC_NAMESPACE + numeric coeff = ::smod(ex_to_numeric(it->coeff), xi); +#endif // ndef NO_GINAC_NAMESPACE if (!coeff.is_zero()) newseq.push_back(expair(it->rest, coeff)); it++; } GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); +#ifndef NO_GINAC_NAMESPACE numeric coeff = GiNaC::smod(ex_to_numeric(overall_coeff), xi); +#else // ndef NO_GINAC_NAMESPACE + numeric coeff = ::smod(ex_to_numeric(overall_coeff), xi); +#endif // ndef NO_GINAC_NAMESPACE return (new add(newseq,coeff))->setflag(status_flags::dynallocated); } @@ -941,14 +984,18 @@ ex mul::smod(const numeric &xi) const #endif // def DO_GINAC_ASSERT mul * mulcopyp=new mul(*this); GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); +#ifndef NO_GINAC_NAMESPACE mulcopyp->overall_coeff = GiNaC::smod(ex_to_numeric(overall_coeff),xi); +#else // ndef NO_GINAC_NAMESPACE + mulcopyp->overall_coeff = ::smod(ex_to_numeric(overall_coeff),xi); +#endif // ndef NO_GINAC_NAMESPACE mulcopyp->clearflag(status_flags::evaluated); mulcopyp->clearflag(status_flags::hash_calculated); return mulcopyp->setflag(status_flags::dynallocated); } -/** Exception thrown by heur_gcd() to signal failure */ +/** Exception thrown by heur_gcd() to signal failure. */ class gcdheu_failed {}; /** Compute GCD of multivariate polynomials using the heuristic GCD algorithm. @@ -995,9 +1042,9 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const numeric mp = p.max_coefficient(), mq = q.max_coefficient(); numeric xi; if (mp > mq) - xi = mq * numTWO() + numTWO(); + xi = mq * _num2() + _num2(); else - xi = mp * numTWO() + numTWO(); + xi = mp * _num2() + _num2(); // 6 tries maximum for (int t=0; t<6; t++) { @@ -1009,7 +1056,7 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const if (!is_ex_exactly_of_type(gamma, fail)) { // Reconstruct polynomial from GCD of mapped polynomials - ex g = exZERO(); + ex g = _ex0(); numeric rxi = xi.inverse(); for (int i=0; !gamma.is_zero(); i++) { ex gi = gamma.smod(xi); @@ -1024,7 +1071,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(exZERO()) < 0) + if (is_ex_exactly_of_type(lc, numeric) && lc.compare(_ex0()) < 0) return -g; else return g; @@ -1049,35 +1096,74 @@ static ex heur_gcd(const ex &a, const ex &b, ex *ca, ex *cb, sym_desc_vec::const ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) { + // 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()) + goto factored_b; +factored_a: + ex g = _ex1(); + ex acc_ca = _ex1(); + ex part_b = b; + for (int i=0; i b.nops()) + goto factored_a; +factored_b: + ex g = _ex1(); + ex acc_cb = _ex1(); + ex part_a = a; + for (int i=0; inormal(sym_lst, repl_lst, level-1)); // Determine common denominator - ex den = exONE(); + ex den = _ex1(); exvector::const_iterator ait = o.begin(), aitend = o.end(); while (ait != aitend) { den = lcm((*ait).denom(false), den, false); @@ -1384,7 +1468,7 @@ ex add::normal(lst &sym_lst, lst &repl_lst, int level) const } // Add fractions - if (den.is_equal(exONE())) + if (den.is_equal(_ex1())) return (new add(o))->setflag(status_flags::dynallocated); else { exvector num_seq; @@ -1443,10 +1527,10 @@ ex power::normal(lst &sym_lst, lst &repl_lst, int level) const } -/** Implementation of ex::normal() for series. It normalizes each coefficient and +/** Implementation of ex::normal() for pseries. It normalizes each coefficient and * replaces the series by a temporary symbol. * @see ex::normal */ -ex series::normal(lst &sym_lst, lst &repl_lst, int level) const +ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) const { epvector new_seq; new_seq.reserve(seq.size()); @@ -1457,7 +1541,7 @@ ex series::normal(lst &sym_lst, lst &repl_lst, int level) const it++; } - ex n = series(var, point, new_seq); + ex n = pseries(var, point, new_seq); return replace_with_symbol(n, sym_lst, repl_lst); } @@ -1484,4 +1568,6 @@ ex ex::normal(int level) const return e; } +#ifndef NO_GINAC_NAMESPACE } // namespace GiNaC +#endif // ndef NO_GINAC_NAMESPACE