00001 00002 /* 00003 * THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE 00004 * CHECKED IT OUT USING THE COMMAND CHECKOUT. 00005 * 00006 * $Id: brent_8c-source.html 2161 2006-05-19 16:55:03Z paulf $ 00007 * 00008 * Revision history: 00009 * $Log$ 00009 * Revision 1.1 2006/05/19 16:55:01 paulf 00009 * first inclusion 00009 * 00010 * Revision 1.1 2000/02/14 18:51:48 lucky 00011 * Initial revision 00012 * 00013 * 00014 */ 00015 00016 #include <math.h> 00017 00018 #define ITMAX 100 00019 #define CGOLD 0.3819660 00020 #define ZEPS 1.0e-10 00021 #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a)) 00022 #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d); 00023 00024 float brent(float ax, float bx, float cx, float (*f)(float), float tol, float *xmin) 00025 { 00026 int iter; 00027 float a,b,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm; 00028 float d=0.0; 00029 float e=0.0; 00030 00031 a=((ax < cx) ? ax : cx); 00032 b=((ax > cx) ? ax : cx); 00033 x=w=v=bx; 00034 fw=fv=fx=(*f)(x); 00035 00036 for (iter=1;iter<=ITMAX;iter++) 00037 { 00038 xm=(float)0.5*(a+b); 00039 00040 00041 tol2=(float)2.0*(tol1=(float)tol*(float)fabs(x)+(float)ZEPS); 00042 00043 if (fabs(x-xm) <= (tol2-0.5*(b-a))) 00044 { 00045 *xmin=x; 00046 return fx; 00047 } 00048 if (fabs(e) > tol1) 00049 { 00050 r=(x-w)*(fx-fv); 00051 q=(x-v)*(fx-fw); 00052 p=(x-v)*q-(x-w)*r; 00053 q=(float)2.0*(q-r); 00054 if (q > 0.0) p = -p; 00055 q=(float)fabs(q); 00056 etemp=e; 00057 e=d; 00058 if (fabs(p) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x)) 00059 d=(float)CGOLD*(e=(x >= xm ? a-x : b-x)); 00060 else 00061 { 00062 d=p/q; 00063 u=x+d; 00064 if (u-a < tol2 || b-u < tol2) 00065 d=(float)SIGN(tol1,xm-x); 00066 } 00067 } 00068 else 00069 d=(float)CGOLD*(e=(x >= xm ? a-x : b-x)); 00070 00071 u=(float)(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d)); 00072 fu=(*f)(u); 00073 if (fu <= fx) 00074 { 00075 if (u >= x) a=x; else b=x; 00076 SHFT(v,w,x,u) 00077 SHFT(fv,fw,fx,fu) 00078 } 00079 else 00080 { 00081 if (u < x) a=u; else b=u; 00082 if (fu <= fw || w == x) 00083 { 00084 v=w; 00085 w=u; 00086 fv=fw; 00087 fw=fu; 00088 } 00089 else if (fu <= fv || v == x || v == w) 00090 { 00091 v=u; 00092 fv=fu; 00093 } 00094 } 00095 } 00096 *xmin=x; 00097 return fx; 00098 } 00099 00100 #undef ITMAX 00101 #undef CGOLD 00102 #undef ZEPS 00103 #undef SIGN