/*
	Combined Newton-Raphson and Bisection root-finder

	Original Function Name:  Inv_Xnrbis()
	Author:  Dr. John Spouge
	Location:  NCBI
	Received:  July 16, 1991
*/
#include <ncbi.h>
#include <gish.h>
#include <gishlib.h>

#define F(x)  ((*f)(x)-y)
#define DF(x) ((*df)(x))
#define ITMAX 100

double fct_bisnr(y,f,df,p,x,q,tol) /* Bisection & Newton-Raphson combination */
/*
Routine locates x : f(x)=y to within absolute error +-tol.
It always calls df(x) after calling f(x) for the same x.
Assumptions:
	(Error) f(p)<=y<=f(q) or f(q)<=y<=f(p)
*/
	double	y;			/* f(x)=y */
	double	(*f) PROTO((double));		/* function */
	double	(*df) PROTO((double));	/* function derivative */
	double	p;			/* end-point */
	double	x;			/* initial root estimate */
	double	q;			/* end-point */
	double	tol;		/* tolerance */
{
	double	temp;	/* for swapping end-points if necessary */
	double	dx;		/* present interval length */
	double	dxold;	/* old interval length */
	double	fx;		/* f(x)-y */
	double	dfx;	/* Df(x) */
	int	j;                       /* loop index */
/* Preliminary checks for improper bracketing and end-point root. */
	if(F(p)*F(q)>0.)
		fatal(GISH_ERR_INVAL, "fct_bisnr: root not bracketed");
	if(F(p)==0.)
		return p;
	else if(F(q)==0.)
		return q;
/* Swaps end-points if necessary to make F(p)<0.<F(q) */
	if(F(p)>0.) {
		temp=p;
		p=q;
		q=temp;
	}
/* Set up the Bisection & Newton-Raphson iteration. */
	if((x-p)*(x-q)>0.)
		x=0.5*(p+q);
	dxold=dx=p-q;
	for(j=1;j<=ITMAX;j++) {
		fx=F(x);    
		if(fx==0.) {               /* Check for termination. */
			return x;
		} else if(fx<0.) {
			p=x;
		} else {
			q=x;
		}
		dfx=DF(x);
/* Check: root out of bounds or bisection faster than Newton-Raphson? */
		if((dfx*(x-p)-fx)*(dfx*(x-q)-fx)>=0. 
			|| 2.*ABS(fx)>ABS(dfx*dx)) {
			dx=dxold;               /* Bisect */
			dxold=0.5*(p-q);
			x=0.5*(p+q);
			if(ABS(dxold)<=tol) 
				return x;
		} else {     
			dx=dxold;               /* Newton-Raphson */
			dxold=fx/dfx;
			x-=dxold;
			if(ABS(dxold)<tol) {
				if(F(x-SIGN(dxold)*tol)*fx<0.)	
					return x;
			}
		}
	}
	fatal(GISH_ERR_INVAL, "fct_bisnr: iterations > ITMAX");
	/*NOTREACHED*/
}