Improve method of setting status_flags::dynallocated.

[ginac.git] / ginac / integral.cpp

diff --git a/ginac/integral.cpp b/ginac/integral.cpp
index 2f54b503bb9ef5b160b60e1f9193faf0c89b30b0..81e7adefcb3208f72af6e61746a1ec1ba239bd46 100644 (file)
--- a/ginac/integral.cpp
+++ b/ginac/integral.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's symbolic  integral. */
 
 /*
- *  GiNaC Copyright (C) 1999-2006 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2015 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<symbol>())
 {}
 
 //////////
@@ -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<symbol>(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
 //////////
@@ -160,11 +158,10 @@ ex integral::eval(int level) const
        if (ea==eb)
                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))
+       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 dynallocate<integral>(eintvar, ea, eb, ef).setflag(status_flags::evaluated);
 }
 
 ex integral::evalf(int level) const
@@ -186,20 +183,17 @@ ex integral::evalf(int level) const
        }
 
        // 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<numeric>(ea) && is_exactly_a<numeric>(eb) 
-                       && is_exactly_a<numeric>(ef.subs(x==12.34).evalf())) {
-               try {
+       // results after substituting a number for the integration variable.
+       if (is_exactly_a<numeric>(ea) && is_exactly_a<numeric>(eb) &&
+           is_exactly_a<numeric>(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<integral>(x, ea, eb, ef);
 }
 
 int integral::max_integration_level = 15;
@@ -223,7 +217,7 @@ struct error_and_integral
 
 struct error_and_integral_is_less
 {
-       bool operator()(const error_and_integral &e1,const error_and_integral &e2)
+       bool operator()(const error_and_integral &e1,const error_and_integral &e2) const
        {
                int c = e1.integral.compare(e2.integral);
                if(c < 0)
@@ -279,7 +273,7 @@ ex adaptivesimpson(const ex & x, const ex & a_in, const ex & b_in, const ex & f,
        fbvec[i] = subsvalue(x, b, f);
        svec[i] = hvec[i]*(favec[i]+4*fcvec[i]+fbvec[i])/3;
        lvec[i] = 1;
-       errorvec[i] = error*svec[i];
+       errorvec[i] = error*abs(svec[i]);
 
        while (i>0) {
                ex fd = subsvalue(x, avec[i]+hvec[i]/2, f);
@@ -292,7 +286,7 @@ ex adaptivesimpson(const ex & x, const ex & a_in, const ex & b_in, const ex & f,
                ex nu4 = fbvec[i];
                ex nu5 = hvec[i];
                // hopefully prevents a crash if the function is zero sometimes.
-               ex nu6 = max(errorvec[i], (s1+s2)*error);
+               ex nu6 = max(errorvec[i], abs(s1+s2)*error);
                ex nu7 = svec[i];
                int nu8 = lvec[i];
                --i;
@@ -410,15 +404,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<integral>(x, newa, newb, newf);
        if (options == 0)
                newint.setflag(status_flags::expanded);
        return newint;
@@ -436,7 +429,7 @@ unsigned integral::return_type() const
        return f.return_type();
 }
 
-const basic* integral::return_type_tinfo() const
+return_type_t integral::return_type_tinfo() const
 {
        return f.return_type_tinfo();
 }
@@ -447,12 +440,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<integral>(x, conja, conjb, conjf);
 }
 
 ex integral::eval_integ() const
@@ -474,4 +466,5 @@ ex integral::eval_integ() const
        return *this;
 }
 
+GINAC_BIND_UNARCHIVER(integral);
 } // namespace GiNaC