mnbrak.h

经典numerical receip 配套代码

struct Mnbrak {
	inline void shft3(Doub &a, Doub &b, Doub &c, const Doub d)
	{
		a=b;
		b=c;
		c=d;
	}
	template <class T>
	void bracket(Doub &ax, Doub &bx, Doub &cx, Doub &fa, Doub &fb, Doub &fc,
		T &func)
	{
		const Doub GOLD=1.618034,GLIMIT=100.0,TINY=1.0e-20;
		Doub fu;
		fa=func(ax);
		fb=func(bx);
		if (fb > fa) {
			SWAP(ax,bx);
			SWAP(fb,fa);
		}
		cx=bx+GOLD*(bx-ax);
		fc=func(cx);
		while (fb > fc) {
			Doub r=(bx-ax)*(fb-fc);
			Doub q=(bx-cx)*(fb-fa);
			Doub u=bx-((bx-cx)*q-(bx-ax)*r)/
				(2.0*SIGN(MAX(abs(q-r),TINY),q-r));
			Doub ulim=bx+GLIMIT*(cx-bx);
			if ((bx-u)*(u-cx) > 0.0) {
				fu=func(u);
				if (fu < fc) {
					ax=bx;
					bx=u;
					fa=fb;
					fb=fu;
					return;
				} else if (fu > fb) {
					cx=u;
					fc=fu;
					return;
				}
				u=cx+GOLD*(cx-bx);
				fu=func(u);
			} else if ((cx-u)*(u-ulim) > 0.0) {
				fu=func(u);
				if (fu < fc) {
					shft3(bx,cx,u,u+GOLD*(u-cx));
					shft3(fb,fc,fu,func(u));
				}
			} else if ((u-ulim)*(ulim-cx) >= 0.0) {
				u=ulim;
				fu=func(u);
			} else {
				u=cx+GOLD*(cx-bx);
				fu=func(u);
			}
			shft3(ax,bx,cx,u);
			shft3(fa,fb,fc,fu);
		}
	}
};