X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fintegral.cpp;h=0b3800488bbad870dec9a30371cdadfb4fa7927d;hp=220047568f160f36a9b522bb34dd6bb1f94f1040;hb=9f3e78a13c72ca73ade4d0829186edcdc1c1b2a0;hpb=12fefbca9b424cb8e9ae05d83883b96e17c7b96e diff --git a/ginac/integral.cpp b/ginac/integral.cpp index 22004756..0b380048 100644 --- a/ginac/integral.cpp +++ b/ginac/integral.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's symbolic integral. */ /* - * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2019 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 @@ -48,8 +48,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integral, basic, ////////// integral::integral() - : inherited(&integral::tinfo_static), - x((new symbol())->setflag(status_flags::dynallocated)) + : x(dynallocate()) {} ////////// @@ -59,7 +58,7 @@ integral::integral() // public integral::integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_) - : inherited(&integral::tinfo_static), x(x_), a(a_), b(b_), f(f_) + : x(x_), a(a_), b(b_), f(f_) { if (!is_a(x)) { throw(std::invalid_argument("first argument of integral must be of type symbol")); @@ -70,8 +69,9 @@ integral::integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_) // archiving ////////// -integral::integral(const archive_node & n, lst & sym_lst) : inherited(n, sym_lst) +void integral::read_archive(const archive_node& n, lst& sym_lst) { + inherited::read_archive(n, sym_lst); n.find_ex("x", x, sym_lst); n.find_ex("a", a, sym_lst); n.find_ex("b", b, sym_lst); @@ -87,8 +87,6 @@ void integral::archive(archive_node & n) const n.add_ex("f", f); } -DEFAULT_UNARCHIVE(integral) - ////////// // functions overriding virtual functions from base classes ////////// @@ -142,64 +140,38 @@ int integral::compare_same_type(const basic & other) const return f.compare(o.f); } -ex integral::eval(int level) const +ex integral::eval() const { - if ((level==1) && (flags & status_flags::evaluated)) + if (flags & status_flags::evaluated) return *this; - if (level == -max_recursion_level) - throw(std::runtime_error("max recursion level reached")); - - ex eintvar = (level==1) ? x : x.eval(level-1); - ex ea = (level==1) ? a : a.eval(level-1); - ex eb = (level==1) ? b : b.eval(level-1); - ex ef = (level==1) ? f : f.eval(level-1); - if (!ef.has(eintvar) && !haswild(ef)) - return eb*ef-ea*ef; + if (!f.has(x) && !haswild(f)) + return b*f-a*f; - if (ea==eb) + if (a==b) return _ex0; - if (are_ex_trivially_equal(eintvar,x) && are_ex_trivially_equal(ea,a) - && are_ex_trivially_equal(eb,b) && are_ex_trivially_equal(ef,f)) - return this->hold(); - return (new integral(eintvar, ea, eb, ef)) - ->setflag(status_flags::dynallocated | status_flags::evaluated); + return this->hold(); } -ex integral::evalf(int level) const +ex integral::evalf() const { - ex ea; - ex eb; - ex ef; - - if (level==1) { - ea = a; - eb = b; - ef = f; - } else if (level == -max_recursion_level) { - throw(runtime_error("max recursion level reached")); - } else { - ea = a.evalf(level-1); - eb = b.evalf(level-1); - ef = f.evalf(level-1); - } + ex ea = a.evalf(); + ex eb = b.evalf(); + ex ef = f.evalf(); // 12.34 is just an arbitrary number used to check whether a number - // results after subsituting a number for the integration variable. - if (is_exactly_a(ea) && is_exactly_a(eb) - && is_exactly_a(ef.subs(x==12.34).evalf())) { - try { + // results after substituting a number for the integration variable. + if (is_exactly_a(ea) && is_exactly_a(eb) && + is_exactly_a(ef.subs(x==12.34).evalf())) { return adaptivesimpson(x, ea, eb, ef); - } catch (runtime_error &rte) {} } - if (are_ex_trivially_equal(a, ea) && are_ex_trivially_equal(b, eb) - && are_ex_trivially_equal(f, ef)) + if (are_ex_trivially_equal(a, ea) && are_ex_trivially_equal(b, eb) && + are_ex_trivially_equal(f, ef)) return *this; else - return (new integral(x, ea, eb, ef)) - ->setflag(status_flags::dynallocated); + return dynallocate(x, ea, eb, ef); } int integral::max_integration_level = 15; @@ -257,7 +229,7 @@ ex adaptivesimpson(const ex & x, const ex & a_in, const ex & b_in, const ex & f, static lookup_map lookup; static symbol ivar("ivar"); ex lookupex = integral(ivar,a,b,f.subs(x==ivar)); - lookup_map::iterator emi = lookup.find(error_and_integral(error, lookupex)); + auto emi = lookup.find(error_and_integral(error, lookupex)); if (emi!=lookup.end()) return emi->second; @@ -410,15 +382,14 @@ ex integral::expand(unsigned options) const return (prefactor*integral(x, newa, newb, rest)).expand(options); } - if (are_ex_trivially_equal(a, newa) && are_ex_trivially_equal(b, newb) - && are_ex_trivially_equal(f, newf)) { + if (are_ex_trivially_equal(a, newa) && are_ex_trivially_equal(b, newb) && + are_ex_trivially_equal(f, newf)) { if (options==0) this->setflag(status_flags::expanded); return *this; } - const basic & newint = (new integral(x, newa, newb, newf)) - ->setflag(status_flags::dynallocated); + const integral & newint = dynallocate(x, newa, newb, newf); if (options == 0) newint.setflag(status_flags::expanded); return newint; @@ -447,12 +418,11 @@ ex integral::conjugate() const ex conjb = b.conjugate(); ex conjf = f.conjugate().subs(x.conjugate()==x); - if (are_ex_trivially_equal(a, conja) && are_ex_trivially_equal(b, conjb) - && are_ex_trivially_equal(f, conjf)) + if (are_ex_trivially_equal(a, conja) && are_ex_trivially_equal(b, conjb) && + are_ex_trivially_equal(f, conjf)) return *this; - return (new integral(x, conja, conjb, conjf)) - ->setflag(status_flags::dynallocated); + return dynallocate(x, conja, conjb, conjf); } ex integral::eval_integ() const @@ -474,4 +444,5 @@ ex integral::eval_integ() const return *this; } +GINAC_BIND_UNARCHIVER(integral); } // namespace GiNaC