struct Mnbrak { inline void shft3(Doub &a, Doub &b, Doub &c, const Doub d) { a=b; b=c; c=d; } template 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); } } };