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
let eps = sqrt epsilon_float
let gold = 0.5 *. (3. -. sqrt 5.)
let gold' = 1. -. gold
exception Found of (float * float)
let default_tol = sqrt epsilon_float
let golden_search ?(tol=default_tol) f a b c =
let rec sort x1 x2 x3 f2 =
if abs_float(x1 -. x3) > tol *. (abs_float x1 +. abs_float x3)
then begin
if abs_float(x1 -. x2) < abs_float(x2 -. x3) then
let x = gold *. x1 +. gold' *. x3 in
next x1 x2 x x3 f2 (f x)
else
let x = gold' *. x1 +. gold *. x3 in
next x1 x x2 x3 (f x) f2
end else
(x2, f2)
and next x0 x1 x2 x3 f1 f2 =
if (f1: float) > f2 then sort x0 x1 x2 f1 else sort x1 x2 x3 f2 in
if not((a < b && b < c) || (a > b && b > c)) then
invalid_arg "Max1D.golden_search: b must be (strictly) between a and c";
let fb = f b in
if fb <= f a || fb <= f c then
invalid_arg "Max1D.golden_search: f(b) must be greater than f(a) and f(c)";
sort a b c fb
let do_brent f a b tol =
let v = ref(a +. gold *. (b -. a)) in
let x = ref !v
and w = ref !v in
let fv = ref(f !v) in
let fx = ref !fv
and fw = ref !fv
and a = ref a and b = ref b in
try
while true do
let m = 0.5 *. (!a +. !b)
and tol_act = eps *. abs_float !x +. tol /. 3. in
if abs_float(!x -. m) +. 0.5 *. (!b -. !a) <= 2. *. tol_act then
raise(Found(!x, !fx));
let gold_step = gold *. (if !x < m then !b -. !x else !a -. !x) in
let new_step =
if abs_float(!x -. !w) >= tol_act then
let r = (!x -. !w) *. (!fx -. !fv)
and q = (!x -. !v) *. (!fx -. !fw) in
let p = ref((!x -. !v) *. q -. (!x -. !w) *. r) in
let q = ref(2. *. (q -. r)) in
if !q > 0. then p := -. !p else q := -. !q;
if abs_float !p < abs_float(gold_step *. !q)
&& !p > !q *. (!a -. !x +. 2. *. tol_act)
&& !p < !q *. (!b -. !x -. 2. *. tol_act) then !p /. !q
else gold_step
else gold_step in
let new_step =
if new_step >= 0. then
if new_step < tol_act then tol_act else new_step
else
if -. new_step < tol_act then -. tol_act else new_step in
let t = !x +. new_step in
let ft = f t in
if ft >= !fx then (
if t < !x then b := !x else a := !x;
v := !w; w := !x; x := t;
fv := !fw; fw := !fx; fx := ft;
)
else (
if t < !x then a := t else b := t;
if ft >= !fw || !w = !x then (
v := !w; w := t;
fv := !fw; fw := ft;
)
else if ft >= !fv || !v = !x || !v = !w then (
v := t;
fv := ft;
)
)
done;
assert false
with Found m -> m
let brent ?(tol=default_tol) f a b =
if tol <= 0. then invalid_arg("Max1D.brent: tol <= 0.");
if (a:float) < b then do_brent f a b tol
else if a > b then do_brent f b a tol
else (a, f a)