00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00036 #ifndef MAT_BISECTION
00037 #define MAT_BISECTION
00038 #include <cmath>
00039 namespace mat {
00047 template<typename Treal>
00048 inline int sign(Treal value) {
00049 if (value > 0)
00050 return 1;
00051 else if (value < 0)
00052 return -1;
00053 else
00054 return 0;
00055 }
00056
00057
00069 template<typename Treal, typename Tfun>
00070 Treal bisection(Tfun const & fun, Treal min, Treal max, Treal const tol) {
00071 int sign_min = sign(fun.eval(min));
00072 int sign_max = sign(fun.eval(max));
00073 if (sign_min == sign_max)
00074 throw Failure("bisection(Tfun&, Treal, Treal, Treal): interval "
00075 "incorrect");
00076 Treal middle = (max + min) / 2;
00077 int sign_middle = sign(fun.eval(middle));
00078 while (template_blas_fabs(max - min) > tol * 2 && sign_middle != 0) {
00079 if (sign_middle == sign_min) {
00080 min = middle;
00081 sign_min = sign_middle;
00082 }
00083 else {
00084 max = middle;
00085 sign_max = sign_middle;
00086 }
00087 middle = (max + min) / 2;
00088 sign_middle = sign(fun.eval(middle));
00089 }
00090 return middle;
00091 }
00092
00093 }
00094 #endif