bisection_search/
lib.rs

1// Copyright (C) The Strand-Braid Authors
2// SPDX-License-Identifier: MIT OR Apache-2.0
3
4use simba::scalar::RealField;
5
6#[derive(Debug, Clone)]
7pub struct Interval<T> {
8    a: T,
9    b: T,
10}
11
12impl<T: RealField + Copy> Interval<T> {
13    pub fn new(a: T, b: T) -> Option<Self> {
14        if b >= a { Some(Self { a, b }) } else { None }
15    }
16}
17
18impl<T> Interval<T> {
19    pub fn a(&self) -> &T {
20        &self.a
21    }
22    pub fn b(&self) -> &T {
23        &self.b
24    }
25}
26
27impl<T: RealField + Copy> Interval<T> {
28    pub fn new_from_range<RB: std::ops::RangeBounds<T>>(range: RB) -> Option<Self> {
29        if let std::ops::Bound::Included(start) = range.start_bound()
30            && let std::ops::Bound::Included(end) = range.end_bound()
31        {
32            return Interval::new(*start, *end);
33        }
34        None
35    }
36}
37
38impl<T: RealField + Copy> Interval<T> {
39    pub fn size(&self) -> T {
40        self.b - self.a
41    }
42}
43
44#[derive(Clone)]
45pub struct BisectionSearch<T, F>
46where
47    F: Fn(&T) -> T,
48{
49    pub interval: Interval<T>,
50    fa: T,
51    fb: T,
52    f: F,
53}
54
55impl<T, F> BisectionSearch<T, F>
56where
57    F: Fn(&T) -> T,
58{
59    pub fn new(interval: Interval<T>, f: F) -> Self {
60        let fa = f(&interval.a);
61        let fb = f(&interval.b);
62        Self {
63            interval,
64            fa,
65            fb,
66            f,
67        }
68    }
69}
70
71impl<T, F> BisectionSearch<T, F>
72where
73    T: RealField + Copy,
74    F: Fn(&T) -> T,
75{
76    pub fn step(mut self) -> Self {
77        let two = T::one() + T::one();
78        let c = (self.interval.a + self.interval.b) / two;
79        let fc = (self.f)(&c);
80        if fc == T::zero() {
81            return BisectionSearch {
82                interval: Interval { a: c, b: c },
83                fa: fc,
84                fb: fc,
85                f: self.f,
86            };
87        }
88
89        if self.fa.is_sign_positive() != fc.is_sign_positive() {
90            self.interval.b = c;
91            self.fb = fc;
92            return self;
93        }
94
95        self.interval.a = c;
96        self.fa = fc;
97        self
98    }
99}
100
101#[cfg(test)]
102mod tests {
103    use crate::*;
104
105    #[test]
106    fn wikipedia_example() {
107        // example at https://en.wikipedia.org/wiki/Bisection_method
108        let mut bisect =
109            BisectionSearch::new(Interval::new(1.0, 2.0).unwrap(), |x| x * x * x - x - 2.0);
110
111        for _ in 0..15 {
112            dbg!((&bisect.interval.a(), &bisect.interval.b()));
113            bisect = bisect.step();
114        }
115
116        assert!(f64::abs(bisect.interval.a - 1.521) < 0.001);
117    }
118}
bisection_search/lib.rs

bisection_search/
lib.rs
