1 /* -*- mode: C -*- */
2 /*
3 IGraph library.
4 Copyright (C) 2007-2012 Gabor Csardi <csardi.gabor@gmail.com>
5 334 Harvard street, Cambridge, MA 02139 USA
6
7 This program is free software; you can redistribute it and/or modify
8 it under the terms of the GNU General Public License as published by
9 the Free Software Foundation; either version 2 of the License, or
10 (at your option) any later version.
11
12 This program is distributed in the hope that it will be useful,
13 but WITHOUT ANY WARRANTY; without even the implied warranty of
14 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 GNU General Public License for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software
19 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
20 02110-1301 USA
21
22 */
23
24 /* from GNU R's zeroin.c, minor modifications by Gabor Csardi */
25
26 /* from NETLIB c/brent.shar with max.iter, add'l info and convergence
27 details hacked in by Peter Dalgaard */
28
29 /*************************************************************************
30 * C math library
31 * function ZEROIN - obtain a function zero within the given range
32 *
33 * Input
34 * double zeroin(ax,bx,f,info,Tol,Maxit)
35 * double ax; Root will be seeked for within
36 * double bx; a range [ax,bx]
37 * double (*f)(double x, void *info); Name of the function whose zero
38 * will be seeked for
39 * void *info; Add'l info passed to f
40 * double *Tol; Acceptable tolerance for the root
41 * value.
42 * May be specified as 0.0 to cause
43 * the program to find the root as
44 * accurate as possible
45 *
46 * int *Maxit; Max. iterations
47 *
48 *
49 * Output
50 * Zeroin returns an estimate for the root with accuracy
51 * 4*EPSILON*abs(x) + tol
52 * *Tol returns estimated precision
53 * *Maxit returns actual # of iterations, or -1 if maxit was
54 * reached without convergence.
55 *
56 * Algorithm
57 * G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
58 * computations. M., Mir, 1980, p.180 of the Russian edition
59 *
60 * The function makes use of the bisection procedure combined with
61 * the linear or quadric inverse interpolation.
62 * At every step program operates on three abscissae - a, b, and c.
63 * b - the last and the best approximation to the root
64 * a - the last but one approximation
65 * c - the last but one or even earlier approximation than a that
66 * 1) |f(b)| <= |f(c)|
67 * 2) f(b) and f(c) have opposite signs, i.e. b and c confine
68 * the root
69 * At every step Zeroin selects one of the two new approximations, the
70 * former being obtained by the bisection procedure and the latter
71 * resulting in the interpolation (if a,b, and c are all different
72 * the quadric interpolation is utilized, otherwise the linear one).
73 * If the latter (i.e. obtained by the interpolation) point is
74 * reasonable (i.e. lies within the current interval [b,c] not being
75 * too close to the boundaries) it is accepted. The bisection result
76 * is used in the other case. Therefore, the range of uncertainty is
77 * ensured to be reduced at least by the factor 1.6
78 *
79 ************************************************************************
80 */
81
82 #include "igraph_nongraph.h"
83 #include "igraph_types.h"
84 #include "igraph_interrupt_internal.h"
85
86 #include <float.h>
87 #include <math.h>
88
89 #define EPSILON DBL_EPSILON
90
igraph_zeroin(igraph_real_t * ax,igraph_real_t * bx,igraph_real_t (* f)(igraph_real_t x,void * info),void * info,igraph_real_t * Tol,int * Maxit,igraph_real_t * res)91 int igraph_zeroin( /* An estimate of the root */
92 igraph_real_t *ax, /* Left border | of the range */
93 igraph_real_t *bx, /* Right border| the root is seeked*/
94 igraph_real_t (*f)(igraph_real_t x, void *info), /* Function under investigation */
95 void *info, /* Add'l info passed on to f */
96 igraph_real_t *Tol, /* Acceptable tolerance */
97 int *Maxit, /* Max # of iterations */
98 igraph_real_t *res) { /* Result is stored here */
99 igraph_real_t a, b, c, /* Abscissae, descr. see above */
100 fa, fb, fc; /* f(a), f(b), f(c) */
101 igraph_real_t tol;
102 int maxit;
103
104 a = *ax; b = *bx; fa = (*f)(a, info); fb = (*f)(b, info);
105 c = a; fc = fa;
106 maxit = *Maxit + 1; tol = * Tol;
107
108 /* First test if we have found a root at an endpoint */
109 if (fa == 0.0) {
110 *Tol = 0.0;
111 *Maxit = 0;
112 *res = a;
113 return 0;
114 }
115 if (fb == 0.0) {
116 *Tol = 0.0;
117 *Maxit = 0;
118 *res = b;
119 return 0;
120 }
121
122 while (maxit--) { /* Main iteration loop */
123 igraph_real_t prev_step = b - a; /* Distance from the last but one
124 to the last approximation */
125 igraph_real_t tol_act; /* Actual tolerance */
126 igraph_real_t p; /* Interpolation step is calcu- */
127 igraph_real_t q; /* lated in the form p/q; divi-
128 * sion operations is delayed
129 * until the last moment */
130 igraph_real_t new_step; /* Step at this iteration */
131
132 IGRAPH_ALLOW_INTERRUPTION();
133
134 if ( fabs(fc) < fabs(fb) ) {
135 /* Swap data for b to be the */
136 a = b; b = c; c = a; /* best approximation */
137 fa = fb; fb = fc; fc = fa;
138 }
139 tol_act = 2 * EPSILON * fabs(b) + tol / 2;
140 new_step = (c - b) / 2;
141
142 if ( fabs(new_step) <= tol_act || fb == (igraph_real_t)0 ) {
143 *Maxit -= maxit;
144 *Tol = fabs(c - b);
145 *res = b;
146 return 0; /* Acceptable approx. is found */
147 }
148
149 /* Decide if the interpolation can be tried */
150 if ( fabs(prev_step) >= tol_act /* If prev_step was large enough*/
151 && fabs(fa) > fabs(fb) ) {
152 /* and was in true direction,
153 * Interpolation may be tried */
154 register igraph_real_t t1, cb, t2;
155 cb = c - b;
156 if ( a == c ) { /* If we have only two distinct */
157 /* points linear interpolation */
158 t1 = fb / fa; /* can only be applied */
159 p = cb * t1;
160 q = 1.0 - t1;
161 } else { /* Quadric inverse interpolation*/
162
163 q = fa / fc; t1 = fb / fc; t2 = fb / fa;
164 p = t2 * ( cb * q * (q - t1) - (b - a) * (t1 - 1.0) );
165 q = (q - 1.0) * (t1 - 1.0) * (t2 - 1.0);
166 }
167 if ( p > (igraph_real_t)0 ) { /* p was calculated with the */
168 q = -q; /* opposite sign; make p positive */
169 } else { /* and assign possible minus to */
170 p = -p; /* q */
171 }
172
173 if ( p < (0.75 * cb * q - fabs(tol_act * q) / 2) /* If b+p/q falls in [b,c]*/
174 && p < fabs(prev_step * q / 2) ) { /* and isn't too large */
175 new_step = p / q;
176 } /* it is accepted
177 * If p/q is too large then the
178 * bisection procedure can
179 * reduce [b,c] range to more
180 * extent */
181 }
182
183 if ( fabs(new_step) < tol_act) { /* Adjust the step to be not less*/
184 if ( new_step > (igraph_real_t)0 ) { /* than tolerance */
185 new_step = tol_act;
186 } else {
187 new_step = -tol_act;
188 }
189 }
190 a = b; fa = fb; /* Save the previous approx. */
191 b += new_step; fb = (*f)(b, info); /* Do step to a new approxim. */
192 if ( (fb > 0 && fc > 0) || (fb < 0 && fc < 0) ) {
193 /* Adjust c for it to have a sign opposite to that of b */
194 c = a; fc = fa;
195 }
196
197 }
198 /* failed! */
199 *Tol = fabs(c - b);
200 *Maxit = -1;
201 *res = b;
202 return IGRAPH_DIVERGED;
203 }
204
205