nekolib/algo/
extremum_float.rs

1//! 三分探索（実数）。
2
3use std::ops::RangeInclusive;
4
5/// 三分探索で極値を探す。
6///
7/// 関数 $f$ の $[x\_l, x\_r]$ における極大値を $x\^\\ast$ として、
8/// $|x-x\^\\ast| \\le \\varepsilon$ なる $x$ を求め、$(x, f(x))$ を返す。
9///
10/// # Requirements
11/// $f$ は凸である。
12///
13/// # Notes
14/// 黄金比を用いて分割する実装のため、関数値を使い回すことができる。
15///
16/// # Complexity
17/// `f` の呼び出しを $\\log\_{\\varphi}(\\frac{x\_r-x\_l}{\\varepsilon}) + 1$ 回行う。
18///
19/// # Suggestions
20/// `f64` に限らずジェネリックにするべき？ 関数の返り値も `T: PartialOrd` にする？
21///
22/// # Examples
23/// $f(x) = x\^x$ の最小値を求める。
24///
25/// $x = 1/e$ のとき、最小値 $e\^{-1/e}$ をとる。
26/// cf. [Wolfram|Alpha](https://www.wolframalpha.com/input/?i=y+%3D+x**x)
27/// ```
28/// use nekolib::algo::extremum_float;
29///
30/// let p = 3.0_f64;
31/// let f = |x: f64| -x.powf(x);
32///
33/// let xl = 0.0;
34/// let xr = 140.0;
35/// let eps = 1.0e-8;
36/// let (x, y) = extremum_float(xl..=xr, eps, f);
37/// let y = -y;
38///
39/// let e = std::f64::consts::E;
40/// assert!(((1.0 / e) - x).abs() < eps);
41/// assert!((e.powf(-1.0 / e) - y).abs() < eps);
42/// ```
43pub fn extremum_float(
44    range: RangeInclusive<f64>,
45    eps: f64,
46    mut f: impl FnMut(f64) -> f64,
47) -> (f64, f64) {
48    let phi: f64 = (1.0 + 5.0_f64.sqrt()) / 2.0;
49    let phi_p1 = phi + 1.0;
50
51    let &(mut xl) = range.start();
52    let &(mut xr) = range.end();
53
54    let iter = ((xr - xl) / eps).log(phi) as u32 + 1;
55
56    let mut xml = (phi * xl + xr) / phi_p1;
57    let mut xmr = (xl + phi * xr) / phi_p1;
58    let mut yml = f(xml);
59    let mut ymr = f(xmr);
60
61    for _ in 0..iter {
62        if yml > ymr {
63            xr = std::mem::replace(&mut xmr, xml);
64            ymr = yml;
65            xml = (phi * xl + xr) / phi_p1;
66            yml = f(xml);
67        } else {
68            xl = std::mem::replace(&mut xml, xmr);
69            yml = ymr;
70            xmr = (xl + phi * xr) / phi_p1;
71            ymr = f(xmr);
72        }
73    }
74    (xml, yml)
75}
nekolib/algo/extremum_float.rs

nekolib/algo/
extremum_float.rs
