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

   1 /* roots/newton.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 /* newton.c -- newton root finding algorithm
  21
  22    This is the classical Newton-Raphson iteration.
  23
  24    x[i+1] = x[i] - f(x[i])/f'(x[i])
  25
  26 */
  27
  28 #include <config.h>
  29
  30 #include <stddef.h>
  31 #include <stdlib.h>
  32 #include <stdio.h>
  33 #include <math.h>
  34 #include <float.h>
  35
  36 #include <gsl/gsl_math.h>
  37 #include <gsl/gsl_errno.h>
  38 #include <gsl/gsl_roots.h>
  39
  40 #include "roots.h"
  41
  42 typedef struct
  43   {
  44     double f, df;
  45   }
  46 newton_state_t;
  47
  48 static int newton_init (void * vstate, gsl_function_fdf * fdf, double * root);
  49 static int newton_iterate (void * vstate, gsl_function_fdf * fdf, double * root);
  50
  51 static int
  52 newton_init (void * vstate, gsl_function_fdf * fdf, double * root)
  53 {
  54   newton_state_t * state = (newton_state_t *) vstate;
  55
  56   const double x = *root ;
  57
  58   state->f = GSL_FN_FDF_EVAL_F (fdf, x);
  59   state->df = GSL_FN_FDF_EVAL_DF (fdf, x) ;
  60
  61   return GSL_SUCCESS;
  62
  63 }
  64
  65 static int
  66 newton_iterate (void * vstate, gsl_function_fdf * fdf, double * root)
  67 {
  68   newton_state_t * state = (newton_state_t *) vstate;
  69
  70   double root_new, f_new, df_new;
  71
  72   if (state->df == 0.0)
  73     {
  74       GSL_ERROR("derivative is zero", GSL_EZERODIV);
  75     }
  76
  77   root_new = *root - (state->f / state->df);
  78
  79   *root = root_new ;
  80
  81   GSL_FN_FDF_EVAL_F_DF(fdf, root_new, &f_new, &df_new);
  82
  83   state->f = f_new ;
  84   state->df = df_new ;
  85
  86   if (!gsl_finite(f_new))
  87     {
  88       GSL_ERROR ("function value is not finite", GSL_EBADFUNC);
  89     }
  90
  91   if (!gsl_finite (df_new))
  92     {
  93       GSL_ERROR ("derivative value is not finite", GSL_EBADFUNC);
  94     }
  95
  96   return GSL_SUCCESS;
  97 }
  98
  99
 100 static const gsl_root_fdfsolver_type newton_type =
 101 {"newton",                              /* name */
 102  sizeof (newton_state_t),
 103  &newton_init,
 104  &newton_iterate};
 105
 106 const gsl_root_fdfsolver_type  * gsl_root_fdfsolver_newton = &newton_type;