1 /* roots/bisection.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 21 /* bisection.c -- bisection root finding algorithm */ 22 23 #include "gsl__config.h" 24 25 #include <stddef.h> 26 #include <stdlib.h> 27 #include <stdio.h> 28 #include <math.h> 29 #include <float.h> 30 31 #include "gsl_math.h" 32 #include "gsl_errno.h" 33 #include "gsl_roots.h" 34 35 #include "gsl_roots__roots.h" 36 37 typedef struct 38 { 39 double f_lower, f_upper; 40 } 41 bisection_state_t; 42 43 static int bisection_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper); 44 static int bisection_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper); 45 46 static int 47 bisection_init (void * vstate, gsl_function * f, double * root, double x_lower, double x_upper) 48 { 49 bisection_state_t * state = (bisection_state_t *) vstate; 50 51 double f_lower, f_upper ; 52 53 *root = 0.5 * (x_lower + x_upper) ; 54 55 SAFE_FUNC_CALL (f, x_lower, &f_lower); 56 SAFE_FUNC_CALL (f, x_upper, &f_upper); 57 58 state->f_lower = f_lower; 59 state->f_upper = f_upper; 60 61 if ((f_lower < 0.0 && f_upper < 0.0) || (f_lower > 0.0 && f_upper > 0.0)) 62 { 63 GSL_ERROR ("endpoints do not straddle y=0", GSL_EINVAL); 64 } 65 66 return GSL_SUCCESS; 67 68 } 69 70 static int 71 bisection_iterate (void * vstate, gsl_function * f, double * root, double * x_lower, double * x_upper) 72 { 73 bisection_state_t * state = (bisection_state_t *) vstate; 74 75 double x_bisect, f_bisect; 76 77 const double x_left = *x_lower ; nied2_state_size(unsigned int dimension)78 const double x_right = *x_upper ; 79 80 const double f_lower = state->f_lower; 81 const double f_upper = state->f_upper; 82 83 if (f_lower == 0.0) 84 { 85 *root = x_left ; 86 *x_upper = x_left; 87 return GSL_SUCCESS; 88 } 89 poly_multiply(const int pa[],int pa_degree,const int pb[],int pb_degree,int pc[],int * pc_degree)90 if (f_upper == 0.0) 91 { 92 *root = x_right ; 93 *x_lower = x_right; 94 return GSL_SUCCESS; 95 } 96 97 x_bisect = (x_left + x_right) / 2.0; 98 99 SAFE_FUNC_CALL (f, x_bisect, &f_bisect); 100 101 if (f_bisect == 0.0) 102 { 103 *root = x_bisect; 104 *x_lower = x_bisect; 105 *x_upper = x_bisect; 106 return GSL_SUCCESS; 107 } 108 109 /* Discard the half of the interval which doesn't contain the root. */ 110 111 if ((f_lower > 0.0 && f_bisect < 0.0) || (f_lower < 0.0 && f_bisect > 0.0)) 112 { 113 *root = 0.5 * (x_left + x_bisect) ; 114 *x_upper = x_bisect; 115 state->f_upper = f_bisect; 116 } 117 else 118 { 119 *root = 0.5 * (x_bisect + x_right) ; 120 *x_lower = x_bisect; 121 state->f_lower = f_bisect; 122 } 123 124 return GSL_SUCCESS; 125 } 126 calculate_v(const int px[],int px_degree,int pb[],int * pb_degree,int v[],int maxv)127 128 static const gsl_root_fsolver_type bisection_type = 129 {"bisection", /* name */ 130 sizeof (bisection_state_t), 131 &bisection_init, 132 &bisection_iterate}; 133 134 const gsl_root_fsolver_type * gsl_root_fsolver_bisection = &bisection_type; 135