1 | #include <math.h> |
2 | #define NRANSI |
3 | #include "nrutil.h" |
4 | #define ITMAX 100 |
5 | #define EPS 3.0e-8 |
6 | |
7 | double zbrent(double (*func)(double), double x1, double x2, double tol) |
8 | { |
9 | int iter; |
10 | double a=x1,b=x2,c=x2,d,e,min1,min2; |
11 | double fa=(*func)(a),fb=(*func)(b),fc,p,q,r,s,tol1,xm; |
12 | |
13 | if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) |
14 | nrerror("Root must be bracketed in zbrent"); |
15 | fc=fb; |
16 | for (iter=1;iter<=ITMAX;iter++) { |
| 1 | Loop condition is true. Entering loop body | |
|
17 | if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) { |
18 | c=a; |
19 | fc=fa; |
20 | e=d=b-a; |
21 | } |
22 | if (fabs(fc) < fabs(fb)) { |
| |
23 | a=b; |
24 | b=c; |
25 | c=a; |
26 | fa=fb; |
27 | fb=fc; |
28 | fc=fa; |
29 | } |
30 | tol1=2.0*EPS*fabs(b)+0.5*tol; |
31 | xm=0.5*(c-b); |
32 | if (fabs(xm) <= tol1 || fb == 0.0) return b; |
| |
33 | if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) { |
| 4 | | Function call argument is an uninitialized value |
|
34 | s=fb/fa; |
35 | if (a == c) { |
36 | p=2.0*xm*s; |
37 | q=1.0-s; |
38 | } else { |
39 | q=fa/fc; |
40 | r=fb/fc; |
41 | p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0)); |
42 | q=(q-1.0)*(r-1.0)*(s-1.0); |
43 | } |
44 | if (p > 0.0) q = -q; |
45 | p=fabs(p); |
46 | min1=3.0*xm*q-fabs(tol1*q); |
47 | min2=fabs(e*q); |
48 | if (2.0*p < (min1 < min2 ? min1 : min2)) { |
49 | e=d; |
50 | d=p/q; |
51 | } else { |
52 | d=xm; |
53 | e=d; |
54 | } |
55 | } else { |
56 | d=xm; |
57 | e=d; |
58 | } |
59 | a=b; |
60 | fa=fb; |
61 | if (fabs(d) > tol1) |
62 | b += d; |
63 | else |
64 | b += SIGN(tol1,xm)((xm) >= 0.0 ? fabs(tol1) : -fabs(tol1)); |
65 | fb=(*func)(b); |
66 | } |
67 | nrerror("Maximum number of iterations exceeded in zbrent"); |
68 | return 0.0; |
69 | } |
70 | #undef ITMAX |
71 | #undef EPS |
72 | #undef NRANSI |