* Make lookup in adaptivesimpson precision-aware. (Chris Dams)

diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi
index 0d337cdb74275d10cb7961edb3d2c0a87f9e830a..1da5cdaa49c04b03e7793141977f5937fd18dd6d 100644 (file)
--- a/doc/tutorial/ginac.texi
+++ b/doc/tutorial/ginac.texi
@@ -1910,7 +1910,7 @@ much work if an expression contains the same integral multiple times,
 a lookup table is used.
 
 If you know that an expression holds an integral, you can get the
-integration variable, the left boundary, right boundary and integrant by
+integration variable, the left boundary, right boundary and integrand by
 respectively calling @code{.op(0)}, @code{.op(1)}, @code{.op(2)}, and
 @code{.op(3)}. Differentiating integrals with respect to variables works
 as expected. Note that it makes no sense to differentiate an integral

diff --git a/ginac/integral.cpp b/ginac/integral.cpp
index bc28ce99aef33b0d36d5f63e19bb382e9e719fac..3daad3efcf8d0419dbb792443df353e668555c13 100644 (file)
--- a/ginac/integral.cpp
+++ b/ginac/integral.cpp
@@ -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)
+       {
+               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;
 
@@ -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;
 }