1 /* integrate.hpp header file
2 *
3 * Copyright Jens Maurer 2000
4 * Distributed under the Boost Software License, Version 1.0. (See
5 * accompanying file LICENSE_1_0.txt or copy at
6 * http://www.boost.org/LICENSE_1_0.txt)
7 *
8 * $Id$
9 *
10 * Revision history
11 * 01 April 2001: Modified to use new <boost/limits.hpp> header. (JMaddock)
12 */
13
14 #ifndef INTEGRATE_HPP
15 #define INTEGRATE_HPP
16
17 #include <boost/limits.hpp>
18
19 template<class UnaryFunction>
20 inline typename UnaryFunction::result_type
trapezoid(UnaryFunction f,typename UnaryFunction::argument_type a,typename UnaryFunction::argument_type b,int n)21 trapezoid(UnaryFunction f, typename UnaryFunction::argument_type a,
22 typename UnaryFunction::argument_type b, int n)
23 {
24 typename UnaryFunction::result_type tmp = 0;
25 for(int i = 1; i <= n-1; ++i)
26 tmp += f(a+(b-a)/n*i);
27 return (b-a)/2/n * (f(a) + f(b) + 2*tmp);
28 }
29
30 template<class UnaryFunction>
31 inline typename UnaryFunction::result_type
simpson(UnaryFunction f,typename UnaryFunction::argument_type a,typename UnaryFunction::argument_type b,int n)32 simpson(UnaryFunction f, typename UnaryFunction::argument_type a,
33 typename UnaryFunction::argument_type b, int n)
34 {
35 typename UnaryFunction::result_type tmp1 = 0;
36 for(int i = 1; i <= n-1; ++i)
37 tmp1 += f(a+(b-a)/n*i);
38 typename UnaryFunction::result_type tmp2 = 0;
39 for(int i = 1; i <= n ; ++i)
40 tmp2 += f(a+(b-a)/2/n*(2*i-1));
41
42 return (b-a)/6/n * (f(a) + f(b) + 2*tmp1 + 4*tmp2);
43 }
44
45 // compute b so that f(b) = y; assume f is monotone increasing
46 template<class UnaryFunction, class T>
47 inline T
invert_monotone_inc(UnaryFunction f,typename UnaryFunction::result_type y,T lower=-1,T upper=1)48 invert_monotone_inc(UnaryFunction f, typename UnaryFunction::result_type y,
49 T lower = -1,
50 T upper = 1)
51 {
52 while(upper-lower > 1e-6) {
53 double middle = (upper+lower)/2;
54 if(f(middle) > y)
55 upper = middle;
56 else
57 lower = middle;
58 }
59 return (upper+lower)/2;
60 }
61
62 // compute b so that I(f(x), a, b) == y
63 template<class UnaryFunction>
64 inline typename UnaryFunction::argument_type
quantil(UnaryFunction f,typename UnaryFunction::argument_type a,typename UnaryFunction::result_type y,typename UnaryFunction::argument_type step)65 quantil(UnaryFunction f, typename UnaryFunction::argument_type a,
66 typename UnaryFunction::result_type y,
67 typename UnaryFunction::argument_type step)
68 {
69 typedef typename UnaryFunction::result_type result_type;
70 if(y >= 1.0)
71 return std::numeric_limits<result_type>::infinity();
72 typename UnaryFunction::argument_type b = a;
73 for(result_type result = 0; result < y; b += step)
74 result += step*f(b);
75 return b;
76 }
77
78
79 #endif /* INTEGRATE_HPP */
80