| /* integrate.hpp header file |
| * |
| * Copyright Jens Maurer 2000 |
| * Distributed under the Boost Software License, Version 1.0. (See |
| * accompanying file LICENSE_1_0.txt or copy at |
| * http://www.boost.org/LICENSE_1_0.txt) |
| * |
| * $Id: integrate.hpp 60755 2010-03-22 00:45:06Z steven_watanabe $ |
| * |
| * Revision history |
| * 01 April 2001: Modified to use new <boost/limits.hpp> header. (JMaddock) |
| */ |
| |
| #ifndef INTEGRATE_HPP |
| #define INTEGRATE_HPP |
| |
| #include <boost/limits.hpp> |
| |
| template<class UnaryFunction> |
| inline typename UnaryFunction::result_type |
| trapezoid(UnaryFunction f, typename UnaryFunction::argument_type a, |
| typename UnaryFunction::argument_type b, int n) |
| { |
| typename UnaryFunction::result_type tmp = 0; |
| for(int i = 1; i <= n-1; ++i) |
| tmp += f(a+(b-a)/n*i); |
| return (b-a)/2/n * (f(a) + f(b) + 2*tmp); |
| } |
| |
| template<class UnaryFunction> |
| inline typename UnaryFunction::result_type |
| simpson(UnaryFunction f, typename UnaryFunction::argument_type a, |
| typename UnaryFunction::argument_type b, int n) |
| { |
| typename UnaryFunction::result_type tmp1 = 0; |
| for(int i = 1; i <= n-1; ++i) |
| tmp1 += f(a+(b-a)/n*i); |
| typename UnaryFunction::result_type tmp2 = 0; |
| for(int i = 1; i <= n ; ++i) |
| tmp2 += f(a+(b-a)/2/n*(2*i-1)); |
| |
| return (b-a)/6/n * (f(a) + f(b) + 2*tmp1 + 4*tmp2); |
| } |
| |
| // compute b so that f(b) = y; assume f is monotone increasing |
| template<class UnaryFunction, class T> |
| inline T |
| invert_monotone_inc(UnaryFunction f, typename UnaryFunction::result_type y, |
| T lower = -1, |
| T upper = 1) |
| { |
| while(upper-lower > 1e-6) { |
| double middle = (upper+lower)/2; |
| if(f(middle) > y) |
| upper = middle; |
| else |
| lower = middle; |
| } |
| return (upper+lower)/2; |
| } |
| |
| // compute b so that I(f(x), a, b) == y |
| template<class UnaryFunction> |
| inline typename UnaryFunction::argument_type |
| quantil(UnaryFunction f, typename UnaryFunction::argument_type a, |
| typename UnaryFunction::result_type y, |
| typename UnaryFunction::argument_type step) |
| { |
| typedef typename UnaryFunction::result_type result_type; |
| if(y >= 1.0) |
| return std::numeric_limits<result_type>::infinity(); |
| typename UnaryFunction::argument_type b = a; |
| for(result_type result = 0; result < y; b += step) |
| result += step*f(b); |
| return b; |
| } |
| |
| |
| #endif /* INTEGRATE_HPP */ |