X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fintegral.cpp;h=2a332083eeb7e97e26bc25451289bcf045dc901b;hp=041f823aaf43eb353e9cebb0c70c892c6c597e6c;hb=7d870583a6bf21a2ffb7b6f051b702064623892e;hpb=695f6ae955ec530cded8f21efd5569df39447f76

diff --git a/ginac/integral.cpp b/ginac/integral.cpp
index 041f823a..2a332083 100644
--- a/ginac/integral.cpp
+++ b/ginac/integral.cpp
@@ -3,7 +3,7 @@
  *  Implementation of GiNaC's symbolic  integral. */
 
 /*
- *  GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
+ *  GiNaC Copyright (C) 1999-2008 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
@@ -17,7 +17,7 @@
  *
  *  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
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  */
 
 #include "integral.h"
@@ -48,7 +48,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(integral, basic,
 //////////
 
 integral::integral()
-		: inherited(TINFO_integral),
+		: inherited(&integral::tinfo_static),
 		x((new symbol())->setflag(status_flags::dynallocated))
 {}
 
@@ -59,7 +59,7 @@ integral::integral()
 // public
 
 integral::integral(const ex & x_, const ex & a_, const ex & b_, const ex & f_)
-		: inherited(TINFO_integral), x(x_), a(a_), b(b_), f(f_)
+		: inherited(&integral::tinfo_static), 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"));
@@ -210,9 +210,32 @@ ex subsvalue(const ex & var, const ex & value, const ex & fun)
 	ex result = fun.subs(var==value).evalf();
 	if (is_a<numeric>(result))
 		return result;
-	throw logic_error("integrant does not evaluate to numeric");
+	throw logic_error("integrand does not evaluate to numeric");
 }
 
+struct error_and_integral
+{
+	error_and_integral(const ex &err, const ex &integ)
+		:error(err), integral(integ){}
+	ex error;
+	ex integral;
+};
+
+struct error_and_integral_is_less
+{
+	bool operator()(const error_and_integral &e1,const error_and_integral &e2) const
+	{
+		int c = e1.integral.compare(e2.integral);
+		if(c < 0)
+			return true;
+		if(c > 0)
+			return false;
+		return ex_is_less()(e1.error, e2.error);
+	}
+};
+
+typedef map<error_and_integral, ex, error_and_integral_is_less> lookup_map;
+
 /** Numeric integration routine based upon the "Adaptive Quadrature" one
   * in "Numerical Analysis" by Burden and Faires. Parameters are integration
   * variable, left boundary, right boundary, function to be integrated and
@@ -220,13 +243,21 @@ ex subsvalue(const ex & var, const ex & value, const ex & fun)
   * after substituting the integration variable by a number. Another thing
   * to note is that this implementation is no good at integrating functions
   * with discontinuities. */
-ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f, const ex & error)
+ex adaptivesimpson(const ex & x, const ex & a_in, const ex & b_in, const ex & f, const ex & error)
 {
-	// use lookup table to be potentially much faster.
-	static exmap lookup;
+	// Check whether boundaries and error are numbers.
+	ex a = is_exactly_a<numeric>(a_in) ? a_in : a_in.evalf();
+	ex b = is_exactly_a<numeric>(b_in) ? b_in : b_in.evalf();
+	if(!is_exactly_a<numeric>(a) || !is_exactly_a<numeric>(b))
+		throw std::runtime_error("For numerical integration the boundaries of the integral should evalf into numbers.");
+	if(!is_exactly_a<numeric>(error))
+		throw std::runtime_error("For numerical integration the error should be a number.");
+
+	// Use lookup table to be potentially much faster.
+	static lookup_map lookup;
 	static symbol ivar("ivar");
 	ex lookupex = integral(ivar,a,b,f.subs(x==ivar));
-	exmap::iterator emi = lookup.find(lookupex);
+	lookup_map::iterator emi = lookup.find(error_and_integral(error, lookupex));
 	if (emi!=lookup.end())
 		return emi->second;
 
@@ -248,7 +279,7 @@ ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f, const
 	fbvec[i] = subsvalue(x, b, f);
 	svec[i] = hvec[i]*(favec[i]+4*fcvec[i]+fbvec[i])/3;
 	lvec[i] = 1;
-	errorvec[i] = integral::relative_integration_error*svec[i];
+	errorvec[i] = error*abs(svec[i]);
 
 	while (i>0) {
 		ex fd = subsvalue(x, avec[i]+hvec[i]/2, f);
@@ -261,7 +292,7 @@ ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f, const
 		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)*integral::relative_integration_error);
+		ex nu6 = max(errorvec[i], abs(s1+s2)*error);
 		ex nu7 = svec[i];
 		int nu8 = lvec[i];
 		--i;
@@ -291,7 +322,7 @@ ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f, const
 		}
 	}
 
-	lookup[lookupex]=app;
+	lookup[error_and_integral(error, lookupex)]=app;
 	return app;
 }
 
@@ -405,7 +436,7 @@ unsigned integral::return_type() const
 	return f.return_type();
 }
 
-unsigned integral::return_type_tinfo() const
+tinfo_t integral::return_type_tinfo() const
 {
 	return f.return_type_tinfo();
 }