The R Project SVN R

/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1999, 2001 the R Development Core Team
 *
 *  This program is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program; if not, write to the Free Software
 *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/* from NETLIB c/brent.shar with max.iter, add'l info and convergence
   details hacked in by Peter Dalgaard */

/*************************************************************************
 *              C math library
 * function ZEROIN - obtain a function zero within the given range
 *
 * Input
 *  double zeroin(ax,bx,f,info,Tol,Maxit)
 *  double ax;          Root will be seeked for within
 *  double bx;          a range [ax,bx]
 *  double (*f)(double x, void *info); Name of the function whose zero
 *                  will be seeked for
 *  void *info;         Add'l info passed to f
 *  double *Tol;            Acceptable tolerance for the root
 *                  value.
 *                  May be specified as 0.0 to cause
 *                  the program to find the root as
 *                  accurate as possible
 *
 *  int *Maxit;         Max. iterations
 *
 *
 * Output
 *  Zeroin returns an estimate for the root with accuracy
 *  4*EPSILON*abs(x) + tol
 *  *Tol returns estimated precision
 *  *Maxit returns actual # of iterations
 *
 * Algorithm
 *  G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
 *  computations. M., Mir, 1980, p.180 of the Russian edition
 *
 *  The function makes use of the bisection procedure combined with
 *  the linear or quadric inverse interpolation.
 *  At every step program operates on three abscissae - a, b, and c.
 *  b - the last and the best approximation to the root
 *  a - the last but one approximation
 *  c - the last but one or even earlier approximation than a that
 *      1) |f(b)| <= |f(c)|
 *      2) f(b) and f(c) have opposite signs, i.e. b and c confine
 *         the root
 *  At every step Zeroin selects one of the two new approximations, the
 *  former being obtained by the bisection procedure and the latter
 *  resulting in the interpolation (if a,b, and c are all different
 *  the quadric interpolation is utilized, otherwise the linear one).
 *  If the latter (i.e. obtained by the interpolation) point is
 *  reasonable (i.e. lies within the current interval [b,c] not being
 *  too close to the boundaries) it is accepted. The bisection result
 *  is used in the other case. Therefore, the range of uncertainty is
 *  ensured to be reduced at least by the factor 1.6
 *
 ************************************************************************
 */

#include <float.h>
#include <math.h>

#include <R_ext/Applic.h>

#define EPSILON DBL_EPSILON

double R_zeroin(            /* An estimate of the root */
    double ax,              /* Left border | of the range   */
    double bx,              /* Right border| the root is seeked*/
    double (*f)(double x, void *info),  /* Function under investigation */
    void *info,             /* Add'l info passed on to f    */
    double *Tol,            /* Acceptable tolerance     */
    int *Maxit)             /* Max # of iterations */
{
    double a,b,c,           /* Abscissae, descr. see above  */
    fa, fb, fc;         /* f(a), f(b), f(c) */
    double tol;
    int maxit;

    a = ax;  b = bx;  fa = (*f)(a, info);  fb = (*f)(b, info);
    c = a;   fc = fa;
    maxit = *Maxit + 1; tol = * Tol;

    while(maxit--)      /* Main iteration loop  */
    {
    double prev_step = b-a;     /* Distance from the last but one
                       to the last approximation    */
    double tol_act;         /* Actual tolerance     */
    double p;           /* Interpolation step is calcu- */
    double q;           /* lated in the form p/q; divi-
                     * sion operations is delayed
                     * until the last moment    */
    double new_step;        /* Step at this iteration   */

    if( fabs(fc) < fabs(fb) )
    {               /* Swap data for b to be the    */
        a = b;  b = c;  c = a;  /* best approximation       */
        fa=fb;  fb=fc;  fc=fa;
    }
    tol_act = 2*EPSILON*fabs(b) + tol/2;
    new_step = (c-b)/2;

    if( fabs(new_step) <= tol_act || fb == (double)0 )
    {
        *Maxit -= maxit;
        *Tol = fabs(c-b);
        return b;           /* Acceptable approx. is found  */
    }

    /* Decide if the interpolation can be tried */
    if( fabs(prev_step) >= tol_act  /* If prev_step was large enough*/
        && fabs(fa) > fabs(fb) ) {  /* and was in true direction,
                     * Interpolation may be tried   */
        register double t1,cb,t2;
        cb = c-b;
        if( a==c ) {        /* If we have only two distinct */
                    /* points linear interpolation  */
        t1 = fb/fa;     /* can only be applied      */
        p = cb*t1;
        q = 1.0 - t1;
        }
        else {          /* Quadric inverse interpolation*/

        q = fa/fc;  t1 = fb/fc;  t2 = fb/fa;
        p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) );
        q = (q-1.0) * (t1-1.0) * (t2-1.0);
        }
        if( p>(double)0 )       /* p was calculated with the op-*/
        q = -q;         /* posite sign; make p positive */
        else            /* and assign possible minus to */
        p = -p;         /* q                */

        if( p < (0.75*cb*q-fabs(tol_act*q)/2) /* If b+p/q falls in [b,c]*/
        && p < fabs(prev_step*q/2) )    /* and isn't too large  */
        new_step = p/q;         /* it is accepted
                         * If p/q is too large then the
                         * bisection procedure can
                         * reduce [b,c] range to more
                         * extent */
    }

    if( fabs(new_step) < tol_act) { /* Adjust the step to be not less*/
        if( new_step > (double)0 )  /* than tolerance       */
        new_step = tol_act;
        else
        new_step = -tol_act;
    }
    a = b;  fa = fb;            /* Save the previous approx. */
    b += new_step;  fb = (*f)(b, info); /* Do step to a new approxim. */
    if( (fb > 0 && fc > 0) || (fb < 0 && fc < 0) ) {
        /* Adjust c for it to have a sign opposite to that of b */
        c = a;  fc = fa;
    }

    }
    /* failed! */
    *Tol = fabs(c-b);
    return b;
}