1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
use ::libc;
extern "C" {}
/* R_zeroin2() is faster for "expensive" f(), in those typical cases where
* f(ax) and f(bx) are available anyway : */
#[no_mangle]
pub unsafe extern "C" fn R_zeroin2(
mut ax: libc::c_double,
mut bx: libc::c_double,
mut fa: libc::c_double,
mut fb: libc::c_double,
mut f: Option<unsafe extern "C" fn(_: libc::c_double, _: *mut libc::c_void) -> libc::c_double>,
mut info: *mut libc::c_void,
mut Tol: *mut libc::c_double,
mut Maxit: *mut libc::c_int,
) -> libc::c_double
/* Max # of iterations */ {
let mut a: libc::c_double = 0.; /* Abscissae, descr. see above, f(c) */
let mut b: libc::c_double = 0.;
let mut c: libc::c_double = 0.;
let mut fc: libc::c_double = 0.;
let mut tol: libc::c_double = 0.;
let mut maxit: libc::c_int = 0;
a = ax;
b = bx;
c = a;
fc = fa;
maxit = *Maxit + 1 as libc::c_int;
tol = *Tol;
/* First test if we have found a root at an endpoint */
if fa == 0.0f64 {
*Tol = 0.0f64;
*Maxit = 0 as libc::c_int;
return a;
}
if fb == 0.0f64 {
*Tol = 0.0f64;
*Maxit = 0 as libc::c_int;
return b;
}
loop {
let fresh0 = maxit;
maxit = maxit - 1;
if !(fresh0 != 0) {
break;
}
/* Main iteration loop */
let mut prev_step: libc::c_double = b - a; /* Distance from the last but one
to the last approximation */
let mut tol_act: libc::c_double = 0.; /* Actual tolerance */
let mut p: libc::c_double = 0.; /* Interpolation step is calcu- */
let mut q: libc::c_double = 0.; /* lated in the form p/q; divi-
* sion operations is delayed
* until the last moment */
let mut new_step: libc::c_double = 0.; /* Step at this iteration */
if fc.abs() < fb.abs() {
/* Swap data for b to be the */
a = b; /* best approximation */
b = c;
c = a;
fa = fb;
fb = fc;
fc = fa
}
tol_act = 2 as libc::c_int as libc::c_double * 2.2204460492503131e-16f64 * b.abs()
+ tol / 2 as libc::c_int as libc::c_double;
new_step = (c - b) / 2 as libc::c_int as libc::c_double;
if new_step.abs() <= tol_act || fb == 0 as libc::c_int as libc::c_double {
*Maxit -= maxit;
*Tol = (c - b).abs();
return b;
/* Acceptable approx. is found */
}
if prev_step.abs() >= tol_act && fa.abs() > fb.abs() {
/* and was in true direction,
* Interpolation may be tried */
let mut t1: libc::c_double = 0.;
let mut cb: libc::c_double = 0.;
let mut t2: libc::c_double = 0.;
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.0f64 - t1
} else {
/* Quadric inverse interpolation*/
q = fa / fc; /* q */
t1 = fb / fc;
t2 = fb / fa;
p = t2 * (cb * q * (q - t1) - (b - a) * (t1 - 1.0f64));
q = (q - 1.0f64) * (t1 - 1.0f64) * (t2 - 1.0f64)
}
if p > 0 as libc::c_int as libc::c_double {
/* p was calculated with the */
q = -q
} else {
/* opposite sign; make p positive */
/* and assign possible minus to */
p = -p
}
if p < 0.75f64 * cb * q - (tol_act * q).abs() / 2 as libc::c_int as libc::c_double
&& p < (prev_step * q / 2 as libc::c_int as libc::c_double).abs()
{
/* and isn't too large */
new_step = p / q
}
}
if new_step.abs() < tol_act {
/* Decide if the interpolation can be tried */
/* Adjust the step to be not less*/
if new_step > 0 as libc::c_int as libc::c_double {
/* than tolerance */
new_step = tol_act
} else {
new_step = -tol_act
}
} /* Save the previous approx. */
a = b; /* Do step to a new approxim. */
fa = fb;
b += new_step;
fb = Some(f.expect("non-null function pointer")).expect("non-null function pointer")(
b, info,
);
if fb > 0 as libc::c_int as libc::c_double && fc > 0 as libc::c_int as libc::c_double
|| fb < 0 as libc::c_int as libc::c_double && fc < 0 as libc::c_int as libc::c_double
{
/* Adjust c for it to have a sign opposite to that of b */
c = a;
fc = fa
}
}
/* failed! */
*Tol = (c - b).abs();
*Maxit = -(1 as libc::c_int);
return b;
}