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 "gsl__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_math.h"
31 #include "gsl_errno.h"
32 #include "gsl_roots.h"
33
34 #include "gsl_roots__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
brent_init(void * vstate,gsl_function * f,double * root,double x_lower,double x_upper)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
brent_iterate(void * vstate,gsl_function * f,double * root,double * x_lower,double * x_upper)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;
245