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
00034 #ifndef MAT_BISECTION
00035 #define MAT_BISECTION
00036 #include <cmath>
00037 namespace mat {
00045 template<typename Treal>
00046 inline int sign(Treal value) {
00047 if (value > 0)
00048 return 1;
00049 else if (value < 0)
00050 return -1;
00051 else
00052 return 0;
00053 }
00054
00055
00067 template<typename Treal, typename Tfun>
00068 Treal bisection(Tfun const & fun, Treal min, Treal max, Treal const tol) {
00069 int sign_min = sign(fun.eval(min));
00070 int sign_max = sign(fun.eval(max));
00071 if (sign_min == sign_max)
00072 throw Failure("bisection(Tfun&, Treal, Treal, Treal): interval "
00073 "incorrect");
00074 Treal middle = (max + min) / 2;
00075 int sign_middle = sign(fun.eval(middle));
00076 while (template_blas_fabs(max - min) > tol * 2 && sign_middle != 0) {
00077 if (sign_middle == sign_min) {
00078 min = middle;
00079 sign_min = sign_middle;
00080 }
00081 else {
00082 max = middle;
00083 sign_max = sign_middle;
00084 }
00085 middle = (max + min) / 2;
00086 sign_middle = sign(fun.eval(middle));
00087 }
00088 return middle;
00089 }
00090
00091 }
00092 #endif