X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=9a3cd159a2bd72858cecfaab982d450c14005213;hp=614b1db1bfc4a00d5630ecfeeadcbee69fdc5956;hb=721cff7b2f8d555e085202f6f3f47801b92d0413;hpb=6b3768e8c544739ae53321539cb4d1e3112ded1b diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 614b1db1..9a3cd159 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -6,9 +6,49 @@ * computation, square-free factorization and rational function normalization. */ +/* + * GiNaC Copyright (C) 1999 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 + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA + */ + #include +#include +#include -#include "ginac.h" +#include "normal.h" +#include "basic.h" +#include "ex.h" +#include "add.h" +#include "constant.h" +#include "expairseq.h" +#include "fail.h" +#include "indexed.h" +#include "inifcns.h" +#include "lst.h" +#include "mul.h" +#include "ncmul.h" +#include "numeric.h" +#include "power.h" +#include "relational.h" +#include "series.h" +#include "symbol.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 @@ -48,8 +88,6 @@ static bool get_first_symbol(const ex &e, const symbol *&x) * Statistical information about symbols in polynomials */ -#include - /** This structure holds information about the highest and lowest degrees * in which a symbol appears in two multivariate polynomials "a" and "b". * A vector of these structures with information about all symbols in @@ -184,7 +222,7 @@ static numeric lcm_of_coefficients_denominators(const ex &e) numeric ex::integer_content(void) const { - ASSERT(bp!=0); + GINAC_ASSERT(bp!=0); return bp->integer_content(); } @@ -204,27 +242,27 @@ numeric add::integer_content(void) const epvector::const_iterator itend = seq.end(); numeric c = numZERO(); while (it != itend) { - ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); - ASSERT(is_ex_exactly_of_type(it->coeff,numeric)); + GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); + GINAC_ASSERT(is_ex_exactly_of_type(it->coeff,numeric)); c = gcd(ex_to_numeric(it->coeff), c); it++; } - ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); + GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); c = gcd(ex_to_numeric(overall_coeff),c); return c; } numeric mul::integer_content(void) const { -#ifdef DOASSERT +#ifdef DO_GINAC_ASSERT epvector::const_iterator it = seq.begin(); epvector::const_iterator itend = seq.end(); while (it != itend) { - ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); + GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); ++it; } -#endif // def DOASSERT - ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); +#endif // def DO_GINAC_ASSERT + GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); return abs(ex_to_numeric(overall_coeff)); } @@ -455,8 +493,6 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) * Remembering */ -#include - typedef pair ex2; typedef pair exbool; @@ -807,7 +843,7 @@ static ex sr_gcd(const ex &a, const ex &b, const symbol *x) numeric ex::max_coefficient(void) const { - ASSERT(bp!=0); + GINAC_ASSERT(bp!=0); return bp->max_coefficient(); } @@ -825,11 +861,11 @@ numeric add::max_coefficient(void) const { epvector::const_iterator it = seq.begin(); epvector::const_iterator itend = seq.end(); - ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); + GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); numeric cur_max = abs(ex_to_numeric(overall_coeff)); while (it != itend) { numeric a; - ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); + GINAC_ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); a = abs(ex_to_numeric(it->coeff)); if (a > cur_max) cur_max = a; @@ -840,15 +876,15 @@ numeric add::max_coefficient(void) const numeric mul::max_coefficient(void) const { -#ifdef DOASSERT +#ifdef DO_GINAC_ASSERT epvector::const_iterator it = seq.begin(); epvector::const_iterator itend = seq.end(); while (it != itend) { - ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); + GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); it++; } -#endif // def DOASSERT - ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); +#endif // def DO_GINAC_ASSERT + GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); return abs(ex_to_numeric(overall_coeff)); } @@ -863,7 +899,7 @@ numeric mul::max_coefficient(void) const ex ex::smod(const numeric &xi) const { - ASSERT(bp!=0); + GINAC_ASSERT(bp!=0); return bp->smod(xi); } @@ -874,7 +910,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 @@ -884,37 +924,49 @@ ex add::smod(const numeric &xi) const epvector::const_iterator it = seq.begin(); epvector::const_iterator itend = seq.end(); while (it != itend) { - ASSERT(!is_ex_exactly_of_type(it->rest,numeric)); + 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++; } - ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); + 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); } ex mul::smod(const numeric &xi) const { -#ifdef DOASSERT +#ifdef DO_GINAC_ASSERT epvector::const_iterator it = seq.begin(); epvector::const_iterator itend = seq.end(); while (it != itend) { - ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); + GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*it),numeric)); it++; } -#endif // def DOASSERT +#endif // def DO_GINAC_ASSERT mul * mulcopyp=new mul(*this); - ASSERT(is_ex_exactly_of_type(overall_coeff,numeric)); - mulcopyp->overall_coeff=::smod(ex_to_numeric(overall_coeff),xi); + 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. @@ -1016,21 +1068,22 @@ 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) { // Some trivial cases - if (a.is_zero()) { + ex aex = a.expand(), bex = b.expand(); + if (aex.is_zero()) { if (ca) *ca = exZERO(); if (cb) *cb = exONE(); return b; } - if (b.is_zero()) { + if (bex.is_zero()) { if (ca) *ca = exONE(); if (cb) *cb = exZERO(); return a; } - if (a.is_equal(exONE()) || b.is_equal(exONE())) { + if (aex.is_equal(exONE()) || bex.is_equal(exONE())) { if (ca) *ca = a; if (cb) @@ -1046,17 +1099,15 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) return a; } #endif - 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 (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(a) / g; + *ca = ex_to_numeric(aex) / g; if (cb) - *cb = ex_to_numeric(b) / g; + *cb = ex_to_numeric(bex) / g; return g; } if (check_args && !a.info(info_flags::rational_polynomial) || !b.info(info_flags::rational_polynomial)) { - cerr << "a=" << a << endl; - cerr << "b=" << b << endl; throw(std::invalid_argument("gcd: arguments must be polynomials over the rationals")); } @@ -1075,40 +1126,40 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args) if (min_ldeg > 0) { ex common = power(*x, min_ldeg); //clog << "trivial common factor " << common << endl; - return gcd((a / common).expand(), (b / common).expand(), ca, cb, false) * common; + 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 = b.content(*x); - ex g = gcd(a, c, ca, cb, false); + ex c = bex.content(*x); + ex g = gcd(aex, c, ca, cb, false); if (cb) - *cb *= b.unit(*x) * b.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 = a.content(*x); - ex g = gcd(c, b, ca, cb, false); + ex c = aex.content(*x); + ex g = gcd(c, bex, ca, cb, false); if (ca) - *ca *= a.unit(*x) * a.primpart(*x, c); + *ca *= aex.unit(*x) * aex.primpart(*x, c); return g; } // Try heuristic algorithm first, fall back to PRS if that failed ex g; try { - g = heur_gcd(a.expand(), b.expand(), ca, cb, var); + g = heur_gcd(aex, bex, ca, cb, var); } catch (gcdheu_failed) { g = *new ex(fail()); } if (is_ex_exactly_of_type(g, fail)) { //clog << "heuristics failed\n"; - g = sr_gcd(a, b, x); + g = sr_gcd(aex, bex, x); if (ca) - divide(a, g, *ca, false); + divide(aex, g, *ca, false); if (cb) - divide(b, g, *cb, false); + divide(bex, g, *cb, false); } return g; } @@ -1257,21 +1308,20 @@ 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 + return replace_with_symbol(*this, sym_lst, repl_lst); else { // complex numeric re = real(), im = 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); - } + 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); + } } /* * Helper function for fraction cancellation (returns cancelled fraction n/d) */ - static ex frac_cancel(const ex &n, const ex &d) { ex num = n; @@ -1450,3 +1500,7 @@ ex ex::normal(int level) const else return e; } + +#ifndef NO_GINAC_NAMESPACE +} // namespace GiNaC +#endif // ndef NO_GINAC_NAMESPACE