3 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
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.
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.
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.
20 /* bracketing.c -- find an initial bracketing interval for a function to minimize */
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_errno.h>
25 #include <gsl/gsl_min.h>
26 #include <gsl/gsl_machine.h>
31 gsl_min_find_bracket(gsl_function *f,double *x_minimum,double * f_minimum,
32 double * x_lower, double * f_lower,
33 double * x_upper, double * f_upper,
36 /* The three following variables must be declared volatile to avoid storage
37 in extended precision registers available on some architecture. The code
38 relies on the ability to compare double values. As the values will be
39 store in regular memory, the extended precision will then be lost and
40 values that are different in extended precision might have equal
41 representation in double precision. This behavior might break the
44 volatile double f_left = *f_lower;
45 volatile double f_right = *f_upper;
46 volatile double f_center;
47 double x_left = *x_lower;
48 double x_right= *x_upper;
50 const double golden = 0.3819660; /* golden = (3 - sqrt(5))/2 */
54 if (f_right >= f_left)
56 x_center = (x_right - x_left) * golden + x_left;
58 SAFE_FUNC_CALL (f, x_center, &f_center);
64 x_right = (x_center - x_left) / golden + x_left;
66 SAFE_FUNC_CALL (f, x_right, &f_right);
71 if (f_center < f_left )
73 if (f_center < f_right)
77 *x_minimum = x_center;
80 *f_minimum = f_center;
81 /* gsl_ieee_printf_double (&f_left);
83 gsl_ieee_printf_double (&f_center);
87 else if (f_center > f_right)
93 x_right = (x_center - x_left) / golden + x_left;
95 SAFE_FUNC_CALL (f, x_right, &f_right);
97 else /* f_center == f_right */
101 x_center = (x_right - x_left) * golden + x_left;
103 SAFE_FUNC_CALL (f, x_center, &f_center);
106 else /* f_center >= f_left */
110 x_center = (x_right - x_left) * golden + x_left;
112 SAFE_FUNC_CALL (f, x_center, &f_center);
115 while (nb_eval < eval_max
116 && (x_right - x_left) > GSL_SQRT_DBL_EPSILON * ( (x_right + x_left) * 0.5 ) + GSL_SQRT_DBL_EPSILON);
119 *x_minimum = x_center;
122 *f_minimum = f_center;