htswanalysis/MACS/lib/gsl/gsl-1.11/roots/brent.c

   1 /* roots/brent.c
   2  *
   3  * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Reid Priedhorsky, Brian Gough
   4  *
   5  * This program is free software; you can redistribute it and/or modify
   6  * it under the terms of the GNU General Public License as published by
   7  * the Free Software Foundation; either version 3 of the License, or (at
   8  * your option) any later version.
   9  *
  10  * This program is distributed in the hope that it will be useful, but
  11  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13  * General Public License for more details.
  14  *
  15  * You should have received a copy of the GNU General Public License
  16  * along with this program; if not, write to the Free Software
  17  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18  */
  19
  20 /* brent.c -- brent root finding algorithm */
  21
  22 #include <config.h>
  23
  24 #include <stddef.h>
  25 #include <stdlib.h>
  26 #include <stdio.h>
  27 #include <math.h>
  28 #include <float.h>
  29
  30 #include <gsl/gsl_math.h>
  31 #include <gsl/gsl_errno.h>
  32 #include <gsl/gsl_roots.h>
  33
  34 #include "roots.h"
  35
  36
  37 typedef struct
  38   {
  39     double a, b, c, d, e;
  40     double fa, fb, fc;
  41   }
  42 brent_state_t;
  43
  44 static int brent_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper);
  45 static int brent_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper);
  46
  47
  48 static int
  49 brent_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper)
  50 {
  51   brent_state_t * state = (brent_state_t *) vstate;
  52
  53   double f_lower, f_upper ;
  54
  55   *root = 0.5 * (x_lower + x_upper) ;
  56
  57   SAFE_FUNC_CALL (f, x_lower, &f_lower);
  58   SAFE_FUNC_CALL (f, x_upper, &f_upper);
  59
  60   state->a = x_lower;
  61   state->fa = f_lower;
  62
  63   state->b = x_upper;
  64   state->fb = f_upper;
  65
  66   state->c = x_upper;
  67   state->fc = f_upper;
  68
  69   state->d = x_upper - x_lower ;
  70   state->e = x_upper - x_lower ;
  71
  72   if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0))
  73     {
  74       GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL);
  75     }
  76
  77   return GSL_SUCCESS;
  78
  79 }
  80
  81 static int
  82 brent_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper)
  83 {
  84   brent_state_t * state = (brent_state_t *) vstate;
  85
  86   double tol, m;
  87
  88   int ac_equal = 0;
  89
  90   double a = state->a, b = state->b, c = state->c;
  91   double fa = state->fa, fb = state->fb, fc = state->fc;
  92   double d = state->d, e = state->e;
  93
  94   if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0))
  95     {
  96       ac_equal = 1;
  97       c = a;
  98       fc = fa;
  99       d = b - a;
 100       e = b - a;
 101     }
 102
 103   if (fabs (fc) < fabs (fb))
 104     {
 105       ac_equal = 1;
 106       a = b;
 107       b = c;
 108       c = a;
 109       fa = fb;
 110       fb = fc;
 111       fc = fa;
 112     }
 113
 114   tol = 0.5 * GSL_DBL_EPSILON * fabs (b);
 115   m = 0.5 * (c - b);
 116
 117   if (fb == 0)
 118     {
 119       *root = b;
 120       *x_lower = b;
 121       *x_upper = b;
 122
 123       return GSL_SUCCESS;
 124     }
 125
 126   if (fabs (m) <= tol)
 127     {
 128       *root = b;
 129
 130       if (b < c)
 131         {
 132           *x_lower = b;
 133           *x_upper = c;
 134         }
 135       else
 136         {
 137           *x_lower = c;
 138           *x_upper = b;
 139         }
 140
 141       return GSL_SUCCESS;
 142     }
 143
 144   if (fabs (e) < tol || fabs (fa) <= fabs (fb))
 145     {
 146       d = m;            /* use bisection */
 147       e = m;
 148     }
 149   else
 150     {
 151       double p, q, r;   /* use inverse cubic interpolation */
 152       double s = fb / fa;
 153
 154       if (ac_equal)
 155         {
 156           p = 2 * m * s;
 157           q = 1 - s;
 158         }
 159       else
 160         {
 161           q = fa / fc;
 162           r = fb / fc;
 163           p = s * (2 * m * q * (q - r) - (b - a) * (r - 1));
 164           q = (q - 1) * (r - 1) * (s - 1);
 165         }
 166
 167       if (p > 0)
 168         {
 169           q = -q;
 170         }
 171       else
 172         {
 173           p = -p;
 174         }
 175
 176       if (2 * p < GSL_MIN (3 * m * q - fabs (tol * q), fabs (e * q)))
 177         {
 178           e = d;
 179           d = p / q;
 180         }
 181       else
 182         {
 183           /* interpolation failed, fall back to bisection */
 184
 185           d = m;
 186           e = m;
 187         }
 188     }
 189
 190   a = b;
 191   fa = fb;
 192
 193   if (fabs (d) > tol)
 194     {
 195       b += d;
 196     }
 197   else
 198     {
 199       b += (m > 0 ? +tol : -tol);
 200     }
 201
 202   SAFE_FUNC_CALL (f, b, &fb);
 203
 204   state->a = a ;
 205   state->b = b ;
 206   state->c = c ;
 207   state->d = d ;
 208   state->e = e ;
 209   state->fa = fa ;
 210   state->fb = fb ;
 211   state->fc = fc ;
 212
 213   /* Update the best estimate of the root and bounds on each
 214      iteration */
 215
 216   *root = b;
 217
 218   if ((fb < 0 && fc < 0) || (fb > 0 && fc > 0))
 219     {
 220       c = a;
 221     }
 222
 223   if (b < c)
 224     {
 225       *x_lower = b;
 226       *x_upper = c;
 227     }
 228   else
 229     {
 230       *x_lower = c;
 231       *x_upper = b;
 232     }
 233
 234   return GSL_SUCCESS ;
 235 }
 236
 237
 238 static const gsl_root_fsolver_type brent_type =
 239 {"brent",                               /* name */
 240  sizeof (brent_state_t),
 241  &brent_init,
 242  &brent_iterate};
 243
 244 const gsl_root_fsolver_type  * gsl_root_fsolver_brent = &brent_type;