bracketing.c

该文件为c++的数学函数库!是一个非常有用的编程工具.它含有各种数学函数,为科学计算、工程应用等程序编写提供方便！

/* min/bracketing.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Fabrice Rossi
 * 
 * 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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* bracketing.c -- find an initial bracketing interval for a function to minimize */

#include <config.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_min.h>
#include <gsl/gsl_machine.h>

#include "min.h"

int 
gsl_min_find_bracket(gsl_function *f,double *x_minimum,double * f_minimum,
                     double * x_lower, double * f_lower, 
                     double * x_upper, double * f_upper,
                     size_t eval_max)
{
  /* The three following variables must be declared volatile to avoid storage
     in extended precision registers available on some architecture. The code
     relies on the ability to compare double values. As the values will be
     store in regular memory, the extended precision will then be lost and
     values that are different in extended precision might have equal
     representation in double precision. This behavior might break the
     algorithm. 
   */
  volatile double f_left = *f_lower;
  volatile double f_right = *f_upper;
  volatile double f_center;
  double x_left = *x_lower;
  double x_right= *x_upper; 
  double x_center;
  const double golden = 0.3819660;      /* golden = (3 - sqrt(5))/2 */
  size_t nb_eval = 0;
  
  
  if (f_right >= f_left) 
    {
      x_center = (x_right - x_left) * golden + x_left;
      nb_eval++;
      SAFE_FUNC_CALL (f, x_center, &f_center);
    }
  else
    {
      x_center = x_right ;
      f_center = f_right ;
      x_right = (x_center - x_left) / golden + x_left;
      nb_eval++;
      SAFE_FUNC_CALL (f, x_right, &f_right);
    }
  
  do
    {
      if (f_center < f_left )
        {
          if (f_center < f_right)
            {
              *x_lower = x_left;
              *x_upper = x_right;
              *x_minimum = x_center;
              *f_lower = f_left;
              *f_upper = f_right;
              *f_minimum = f_center;
/*            gsl_ieee_printf_double (&f_left);
              printf(" ");
              gsl_ieee_printf_double (&f_center);
              printf("\n");*/
              return GSL_SUCCESS;
            }
          else if (f_center > f_right)
            {
              x_left = x_center;
              f_left = f_center;
              x_center = x_right;
              f_center = f_right;
              x_right = (x_center - x_left) / golden + x_left;
              nb_eval++;
              SAFE_FUNC_CALL (f, x_right, &f_right);
            }
          else /* f_center == f_right */
            {
              x_right = x_center;
              f_right = f_center;
              x_center = (x_right - x_left) * golden + x_left;
              nb_eval++;
              SAFE_FUNC_CALL (f, x_center, &f_center);
            }
        }
      else /* f_center >= f_left */
        {
          x_right = x_center;
          f_right = f_center;
          x_center = (x_right - x_left) * golden + x_left;
          nb_eval++;
          SAFE_FUNC_CALL (f, x_center, &f_center);
        }
    }
  while (nb_eval < eval_max 
         && (x_right - x_left) > GSL_SQRT_DBL_EPSILON * ( (x_right + x_left) * 0.5 ) + GSL_SQRT_DBL_EPSILON);
  *x_lower = x_left;
  *x_upper = x_right;
  *x_minimum = x_center;
  *f_lower = f_left;
  *f_upper = f_right;
  *f_minimum = f_center;
  return GSL_FAILURE;
}