The R Project SVN – R – Diff – /trunk/src/appl/zeroin.c


/*
 *  R : A Computer Language for Statistical Data Analysis
 *  Copyright (C) 1999, 2001 the R Development Core Team
 *
 *  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., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */
 
/* from NETLIB c/brent.shar with max.iter, add'l info and convergence
   details hacked in by Peter Dalgaard */
 
/*************************************************************************
 *			    C math library
 * function ZEROIN - obtain a function zero within the given range
 *
 * Input
 *	double zeroin(ax,bx,f,info,Tol,Maxit)
 *	double ax;			Root will be seeked for within
 *	double bx;			a range [ax,bx]
 *	double (*f)(double x, void *info); Name of the function whose zero
 *					will be seeked for
 *	void *info;			Add'l info passed to f
 *	double *Tol;			Acceptable tolerance for the root
 *					value.
 *					May be specified as 0.0 to cause
 *					the program to find the root as
 *					accurate as possible
 *
 *	int *Maxit;			Max. iterations
 *
 *
 * Output
 *	Zeroin returns an estimate for the root with accuracy
 *	4*EPSILON*abs(x) + tol
 *	*Tol returns estimated precision
 *	*Maxit returns actual # of iterations
 *
 * Algorithm
 *	G.Forsythe, M.Malcolm, C.Moler, Computer methods for mathematical
 *	computations. M., Mir, 1980, p.180 of the Russian edition
 *
 *	The function makes use of the bisection procedure combined with
 *	the linear or quadric inverse interpolation.
 *	At every step program operates on three abscissae - a, b, and c.
 *	b - the last and the best approximation to the root
 *	a - the last but one approximation
 *	c - the last but one or even earlier approximation than a that
 *		1) |f(b)| <= |f(c)|
 *		2) f(b) and f(c) have opposite signs, i.e. b and c confine
 *		   the root
 *	At every step Zeroin selects one of the two new approximations, the
 *	former being obtained by the bisection procedure and the latter
 *	resulting in the interpolation (if a,b, and c are all different
 *	the quadric interpolation is utilized, otherwise the linear one).
 *	If the latter (i.e. obtained by the interpolation) point is
 *	reasonable (i.e. lies within the current interval [b,c] not being
 *	too close to the boundaries) it is accepted. The bisection result
 *	is used in the other case. Therefore, the range of uncertainty is
 *	ensured to be reduced at least by the factor 1.6
 *
 ************************************************************************
 */
 
#include <float.h>
 
#include <math.h>
 
 
 
 
#include <R_ext/Applic.h>
 
#define EPSILON DBL_EPSILON
 
double R_zeroin(			/* An estimate of the root */
    double ax,				/* Left border | of the range	*/
    double bx,				/* Right border| the root is seeked*/
    double (*f)(double x, void *info),	/* Function under investigation	*/
    void *info,				/* Add'l info passed on to f	*/
    double *Tol,			/* Acceptable tolerance		*/
    int *Maxit)				/* Max # of iterations */
{
    double a,b,c,			/* Abscissae, descr. see above	*/
	fa, fb, fc;			/* f(a), f(b), f(c) */
    double tol;
    int maxit;
 
    a = ax;  b = bx;  fa = (*f)(a, info);  fb = (*f)(b, info);
    c = a;   fc = fa;
    maxit = *Maxit + 1; tol = * Tol;
 
    while(maxit--)		/* Main iteration loop	*/
    {
	double prev_step = b-a;		/* Distance from the last but one
					   to the last approximation	*/
	double tol_act;			/* Actual tolerance		*/
	double p;			/* Interpolation step is calcu- */
	double q;			/* lated in the form p/q; divi-
					 * sion operations is delayed
					 * until the last moment	*/
	double new_step;		/* Step at this iteration	*/
 
	if( fabs(fc) < fabs(fb) )
	{				/* Swap data for b to be the	*/
	    a = b;  b = c;  c = a;	/* best approximation		*/
	    fa=fb;  fb=fc;  fc=fa;
	}
	tol_act = 2*EPSILON*fabs(b) + tol/2;
	new_step = (c-b)/2;
 
	if( fabs(new_step) <= tol_act || fb == (double)0 )
	{
	    *Maxit -= maxit;
	    *Tol = fabs(c-b);
	    return b;			/* Acceptable approx. is found	*/
	}
 
	/* Decide if the interpolation can be tried	*/
	if( fabs(prev_step) >= tol_act	/* If prev_step was large enough*/
	    && fabs(fa) > fabs(fb) ) {	/* and was in true direction,
					 * Interpolation may be tried	*/
	    register double t1,cb,t2;
	    cb = c-b;
	    if( a==c ) {		/* If we have only two distinct	*/
					/* points linear interpolation	*/
		t1 = fb/fa;		/* can only be applied		*/
		p = cb*t1;
		q = 1.0 - t1;
	    }
	    else {			/* Quadric inverse interpolation*/
 
		q = fa/fc;  t1 = fb/fc;	 t2 = fb/fa;
		p = t2 * ( cb*q*(q-t1) - (b-a)*(t1-1.0) );
		q = (q-1.0) * (t1-1.0) * (t2-1.0);
	    }
	    if( p>(double)0 )		/* p was calculated with the op-*/
		q = -q;			/* posite sign; make p positive	*/
	    else			/* and assign possible minus to	*/
		p = -p;			/* q				*/
 
	    if( p < (0.75*cb*q-fabs(tol_act*q)/2) /* If b+p/q falls in [b,c]*/
		&& p < fabs(prev_step*q/2) )	/* and isn't too large	*/
		new_step = p/q;			/* it is accepted
						 * If p/q is too large then the
						 * bisection procedure can
						 * reduce [b,c] range to more
						 * extent */
	}
 
	if( fabs(new_step) < tol_act) {	/* Adjust the step to be not less*/
	    if( new_step > (double)0 )	/* than tolerance		*/
		new_step = tol_act;
	    else
		new_step = -tol_act;
	}
	a = b;	fa = fb;			/* Save the previous approx. */
	b += new_step;	fb = (*f)(b, info);	/* Do step to a new approxim. */
	if( (fb > 0 && fc > 0) || (fb < 0 && fc < 0) ) {
	    /* Adjust c for it to have a sign opposite to that of b */
	    c = a;  fc = fa;
	}
 
    }
    /* failed! */
    *Tol = fabs(c-b);
    return b;
}

The R Project SVN R

(root)/trunk/src/appl/zeroin.c @ 28362 – Rev 19500 → 24068

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