30 #define GSL_SQRT_DBL_EPSILON 1.e-4
31 #define GSL_DBL_EPSILON 1.e-8
47 static int brent(
void *vstate,
double (*f)(
double),
double *x_minimum,
48 double *f_minimum,
double *x_lower,
double *f_lower,
49 double *x_upper,
double *f_upper)
51 brent_state_t *
state = (brent_state_t *)vstate;
53 const double x_left = *x_lower;
54 const double x_right = *x_upper;
56 const double z = *x_minimum;
60 const double v =
state->v;
61 const double w =
state->w;
62 const double f_v =
state->f_v;
63 const double f_w =
state->f_w;
64 const double f_z = *f_minimum;
66 const double golden = 0.3819660;
68 const double w_lower = (z - x_left);
69 const double w_upper = (x_right - z);
73 double p = 0, q = 0,
r = 0;
75 const double midpoint = 0.5 * (x_left + x_right);
77 if (fabs(e) > tolerance) {
80 r = (z - w) * (f_z - f_v);
81 q = (z - v) * (f_z - f_w);
82 p = (z - v) * q - (z - w) *
r;
96 if (fabs(p) < fabs(0.5 * q *
r) && p < q * w_lower && p < q * w_upper) {
97 double t2 = 2 * tolerance;
102 if ((u - x_left) < t2 || (x_right - z) < t2) {
103 d = (z < midpoint) ? tolerance : -tolerance;
107 e = (z < midpoint) ? x_right - z : -(z - x_left);
111 if (fabs(d) >= tolerance) {
115 u = z + ((d > 0) ? tolerance : -tolerance);
136 else if (f_u < f_z) {
154 else if (f_u <= f_w || w == z) {
161 else if (f_u <= f_v || v == z || v == w) {
177 double x_minimum = (x_upper + x_lower) / 2.;
183 brent_state_t itstate;
184 const double golden = 0.3819660;
185 double v = x_lower + golden * (x_upper - x_lower);
189 f_lower = (*f)(x_lower);
190 f_upper = (*f)(x_upper);
191 f_minimum = (*f)(x_minimum);
206 for (i = 0; i < maxiter; i++) {
207 brent(&itstate, f, &x_minimum, &f_minimum, &x_lower, &f_lower, &x_upper,
#define GSL_SQRT_DBL_EPSILON
double brent_iterate(double(*f)(double), double x_lower, double x_upper, int maxiter)