00001 #ifndef NR_DBRENT_H
00002 #define NR_DBRENT_H
00003
00004
00005 #include <math.h>
00006 #include "utils.h"
00007
00008
00009 namespace nr {
00010
00011
00018 template < class Real, class Func, class DFunc >
00019 Real dbrent (
00020 Real ax, Real bx, Real cx,
00021 Func (*f),
00022 DFunc (*df),
00023 Real tol,
00024 Real *xmin,
00025 unsigned int ITMAX = 10000,
00026 Real ZEPS = 1.0e-10
00027 )
00028 {
00029 #define MOV3(a,b,c, d,e,f) (a)=(d);(b)=(e);(c)=(f);
00030 int ok1,ok2;
00031 Real a,b,d,d1,d2,du,dv,dw,dx,e=0.0;
00032 Real fu,fv,fw,fx,olde,tol1,tol2,u,u1,u2,v,w,x,xm;
00033
00034 a=(ax < cx ? ax : cx);
00035 b=(ax > cx ? ax : cx);
00036 x=w=v=bx;
00037 fw=fv=fx=(*f)(x);
00038 dw=dv=dx=(*df)(x);
00039 for (std::size_t iter=1;iter<=ITMAX;iter++) {
00040 xm=0.5*(a+b);
00041 tol1=tol*fabs(x)+ZEPS;
00042 tol2=2.0*tol1;
00043 if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
00044 *xmin=x;
00045 return fx;
00046 }
00047 if (fabs(e) > tol1) {
00048 d1=2.0*(b-a);
00049 d2=d1;
00050 if (dw != dx) d1=(w-x)*dx/(dx-dw);
00051 if (dv != dx) d2=(v-x)*dx/(dx-dv);
00052 u1=x+d1;
00053 u2=x+d2;
00054 ok1 = (a-u1)*(u1-b) > 0.0 && dx*d1 <= 0.0;
00055 ok2 = (a-u2)*(u2-b) > 0.0 && dx*d2 <= 0.0;
00056 olde=e;
00057 e=d;
00058 if (ok1 || ok2) {
00059 if (ok1 && ok2)
00060 d=(fabs(d1) < fabs(d2) ? d1 : d2);
00061 else if (ok1)
00062 d=d1;
00063 else
00064 d=d2;
00065 if (fabs(d) <= fabs(0.5*olde)) {
00066 u=x+d;
00067 if (u-a < tol2 || b-u < tol2)
00068 d=SIGN(tol1,xm-x);
00069 } else {
00070 d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
00071 }
00072 } else {
00073 d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
00074 }
00075 } else {
00076 d=0.5*(e=(dx >= 0.0 ? a-x : b-x));
00077 }
00078 if (fabs(d) >= tol1) {
00079 u=x+d;
00080 fu=(*f)(u);
00081 } else {
00082 u=x+SIGN(tol1,d);
00083 fu=(*f)(u);
00084 if (fu > fx) {
00085 *xmin=x;
00086 return fx;
00087 }
00088 }
00089 du=(*df)(u);
00090 if (fu <= fx) {
00091 if (u >= x) a=x; else b=x;
00092 MOV3(v,fv,dv, w,fw,dw)
00093 MOV3(w,fw,dw, x,fx,dx)
00094 MOV3(x,fx,dx, u,fu,du)
00095 } else {
00096 if (u < x) a=u; else b=u;
00097 if (fu <= fw || w == x) {
00098 MOV3(v,fv,dv, w,fw,dw)
00099 MOV3(w,fw,dw, u,fu,du)
00100 } else if (fu < fv || v == x || v == w) {
00101 MOV3(v,fv,dv, u,fu,du)
00102 }
00103 }
00104 }
00105 nrerror("Too many iterations in routine dbrent");
00106 return 0.0;
00107 #undef MOV3
00108 }
00109
00110 }
00111 #endif
00112
Documentation 
1.3.6