nekolib/algo/
extremum.rs

1//! 三分探索。
2
3use std::ops::Range;
4
5/// 三分探索で極値を探す。
6///
7/// # Notes
8/// [`extremum`] を参照されたい。
9///
10/// # Examples
11/// ```
12/// use nekolib::algo::extremum_slice;
13///
14/// let buf = [1, 3, 4, 6, 5, 2, 0, 1, 3];
15/// //         <------ f ------>
16/// //                  <------ g ------>
17///
18/// let f = |&x: &i32| x * x;
19/// assert_eq!(extremum_slice(&buf[..6], f), (3_usize, 36));
20/// let g = |&x: &i32| -x;
21/// assert_eq!(extremum_slice(&buf[3..], g), (3_usize, 0));
22/// ```
23pub fn extremum_slice<T: Ord>(
24    buf: &[T],
25    mut f: impl FnMut(&T) -> T,
26) -> (usize, T) {
27    extremum(0..buf.len(), |i| f(&buf[i]))
28}
29
30/// 三分探索で極値を探す。
31///
32/// 離散値の区間 $[l, r)$ において、以下を満たす $i$ および $f(i)$ を返す。
33/// $$ f(i-1) < f(i)\\text{ and } f(i) > f(i+1). $$
34///
35/// # Requirements
36/// 凸である。すなわち、ある $i$ が存在して以下の二つが成り立つ。
37/// - ${}^\\forall j \\in [l, i)$ に対して $f(j) < f(j+1)$。
38/// - ${}^\\forall j \\in [i, r-1)$ に対して $f(j) > f(j+1)$。
39///
40/// # Implementation notes
41/// 連続値の場合における黄金比分割のように、Fibonacci
42/// 数列に基づいて区間を分割していくため、関数値を使い回すことができる。
43///
44/// 関数 $f$ の呼び出し回数を、区間を三等分する素直な実装と比較する。
45/// 三等分の実装では
46/// $2\\cdot\\log_{3/2}(r-l)+O(1)$ 回（係数の $2$ に注意）であり、こちらは
47/// $\\log_{\\varphi}(r-l)+O(1)$ 回である。
48/// ただし $\\varphi$ は黄金比 $(1+\\sqrt{5})/2 = 1.618\\dots$ である。
49/// $$ \\sqrt{3/2} < 1.225 < 1.618 < \\varphi $$
50/// であり、$\\log$ の底は大きい方がうれしいので、こちらの実装の方が定数倍が軽い。
51///
52/// コード長はやや長くなるものの、単純な例での実測では三等分のものよりわずかに高速であった。
53///
54/// # Complexity
55/// $F\_0 = 1$, $F\_1 = 2$, $F\_i = F\_{i-1} + F\_{i-2}$ ($i \\ge 2$) で定義される数列 $\\{F\_k\\}$ を考える。
56/// 区間幅 $n$ がある $k$ に対して $n \\le F\_k$ と抑えられるとき、$f$ の呼び出しを高々
57/// $k$ 回、関数値同士の比較を高々 $k-1$ 回行う。
58/// なお、この $k$ は $\\lceil\\log\_{\\varphi}(n)\\rceil$ で抑えられる。
59///
60/// # Suggestions
61/// 引数は `Range<usize>` を渡すことにしているものの、実際には
62/// `RangeBounds<I: {integer}>` を渡せるようにする方がよさそう？
63/// ただし、両端とも `Unbounded` であっては困りそう（特に多倍長を視野に入れる場合？）。
64/// 多倍長だと `Copy` がないから、計算結果自体を使い回せても `.clone()` でつらい？
65///
66/// # Examples
67/// ```
68/// use nekolib::algo::extremum;
69///
70/// let buf = [1, 3, 4, 6, 5, 2, 0, 1, 3];
71/// //         <------ f ------>
72/// //                  <------ g ------>
73///
74/// let f = |i: usize| buf[i] * buf[i];
75/// assert_eq!(extremum(0..6, f), (3_usize, 36));
76/// let g = |i: usize| -buf[i];
77/// assert_eq!(extremum(3..8, g), (6_usize, 0));
78/// ```
79///
80/// ```
81/// use nekolib::algo::extremum;
82///
83/// let n = 1500;
84/// for k in 0..n {
85///     let mut count = 0;
86///     let f = |i| { count += 1; -(i as i32 - k as i32).abs() };
87///     assert_eq!(extremum(0..n, f), (k, 0));
88///     assert!(count <= 15);
89/// }
90/// ```
91pub fn extremum<T: Ord>(
92    Range { start, end }: Range<usize>,
93    mut f: impl FnMut(usize) -> T,
94) -> (usize, T) {
95    let n = end - start;
96    if n == 0 {
97        panic!("range must be non-empty");
98    }
99    if n == 1 {
100        return (start, f(start));
101    }
102
103    let (mut i0, mut i1) =
104        std::iter::successors(Some((1, 2)), |&(i, j)| Some((j, i + j)))
105            .find(|&(i, j)| i + j > n)
106            .unwrap();
107    let mut g = |i| if i <= n { Some(f(start + i - 1)) } else { None }; // None means -inf
108    let mut d = i0;
109    let mut g0 = g(i0);
110    let mut g1 = g(i1);
111    while d > 1 {
112        match (g0, g1) {
113            (Some(f0), Some(f1)) if f0 < f1 => {
114                // |lo  i0 < i1     hi|
115                //     |lo   i0 i1  hi|
116                let tmp = i0 + d;
117                i0 = std::mem::replace(&mut i1, tmp);
118                g0 = Some(f1);
119                g1 = g(i1);
120            }
121            (Some(f0), _) => {
122                // |lo     i0 > i1  hi|
123                // |lo  i0 i1   hi|
124                let tmp = i1 - d;
125                i1 = std::mem::replace(&mut i0, tmp);
126                g1 = Some(f0);
127                g0 = g(i0);
128            }
129            (None, _) => unreachable!(),
130        }
131        d -= i1 - i0;
132    }
133
134    match (g0, g1) {
135        (Some(f0), Some(f1)) if f0 < f1 => (start + i1 - 1, f1),
136        (Some(f0), _) => (start + i0 - 1, f0),
137        (None, _) => unreachable!(),
138    }
139}
140
141#[test]
142fn extremum_count() {
143    let mut fl = 1;
144    let mut fr = 2;
145    for i in 1..=15 {
146        for n in fl..fr {
147            for k in 0..n {
148                let mut count = 0;
149                let f = |i| {
150                    count += 1;
151                    -(i as i32 - k as i32).abs()
152                };
153                let res = extremum(0..n, f);
154                assert_eq!(res, (k, 0));
155                assert!(count <= i);
156            }
157        }
158        let tmp = fl + fr;
159        fl = std::mem::replace(&mut fr, tmp);
160    }
161}
nekolib/algo/extremum.rs

nekolib/algo/
extremum.rs
