[GiNaC-devel] Patch for adaptivesimpson.

[Chris Dams]

Dear developers,

Could you apply the attached patch to integral.cpp? It adds some checks to
adaptivesimpson to see whether the boundaries of integration and the
integration error really are numbers and, as discussed with Richy, it
makes the lookup table take into account that adaptivesimpson may be
called with the same integral but different precision.

Best wishes,
Chris
-------------- next part --------------
Index: ginac/integral.cpp
===================================================================
RCS file: /home/cvs/GiNaC/ginac/integral.cpp,v
retrieving revision 1.6
diff -c -r1.6 integral.cpp
*** ginac/integral.cpp	13 Sep 2005 20:41:25 -0000	1.6
--- ginac/integral.cpp	15 Sep 2005 16:39:04 -0000
***************
*** 213,218 ****
--- 213,241 ----
  	throw logic_error("integrant 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,232 ****
    * 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)
  {
! 	// use lookup table to be potentially much faster.
! 	static exmap lookup;
  	static symbol ivar("ivar");
  	ex lookupex = integral(ivar,a,b,f.subs(x==ivar));
! 	exmap::iterator emi = lookup.find(lookupex);
  	if (emi!=lookup.end())
  		return emi->second;
  
--- 243,263 ----
    * 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_in, const ex & b_in, const ex & f, const ex & error)
  {
! 	// 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));
! 	lookup_map::iterator emi = lookup.find(error_and_integral(error, lookupex));
  	if (emi!=lookup.end())
  		return emi->second;
  
***************
*** 291,297 ****
  		}
  	}
  
! 	lookup[lookupex]=app;
  	return app;
  }
  
--- 322,328 ----
  		}
  	}
  
! 	lookup[error_and_integral(error, lookupex)]=app;
  	return app;
  }