nlbit.h

/********************************************************************/
/* Copyright (c) 2015 DG Technologies Inc. and Yuzi Mizuno          */
/* All rights reserved.                                             */
/********************************************************************/
#ifndef _MGNLBIT_HH_
#define _MGNLBIT_HH_

#include <math.h>

template<class func>
double mgNlbit(
    func& fn,   
    double xl,  
    double xr,  
    double eps, 
    int itr,    
    int& ier    
){
    ier = 0;
    double epsHalf = eps*.5;
    double x = xl;
    double fl = fn(x);
    if (fabs(fl)<=eps) return x;
    x = xr;
    double fr = fn(x);
    if(fabs(fr)<=eps) return x;
    if(fr * fl > 0.) {
        x = (xl + xr) * .5;
        ier = 2;
        return x;
    }

    for(int n=1; n<=itr; ++n){
        x=(xl+xr)*.5;
        double f = fn(x);
        if(fabs(f)<=eps)
            return x;

        if(f*fl < 0.){
            xr = x;
            fr = f;
        }else{
            xl = x;
            fl = f;
        }
        double df=fl-fr;
        if(fabs(df)<=epsHalf){
            x=(xr+xl)*.5;
            return x;
        }
        double dx = (xr-xl)*fl/df;
        x = xl+dx;
        f = fn(x);
        if (fabs(f)<=eps)
            return x;

        if(fl*f<0.){
            xr = x;
            fr = f;
        }else{
            xl = x;
            fl = f;
        }
    }
    ier = 1;
    return x;
}
 // end of ALGORITHM group
#endif