4
\$\begingroup\$

In an attempt to learn Rust, I've written up implementations of the bisection method and Newton's method for finding roots of an equation. Both methods come in two variants: the first one searches for a single root on a given interval (with assumption that there is at most one root there), and the second one tries to find all roots in the interval by splitting it into a number of smaller intervals, and calling the first variant on each.

Everything seems to work as intended, but my main concern is whether or not it's idiomatic Rust.

All of the code can be seen here, in the frozen branch review-2018-09-01 (in particular, src/roots.rs).

First, here's the Epsilon trait for approximately comparing things (because comparing floats exactly is almost always a bad idea), with a blanket implementation for everything that implements Float and Signed from the num package:

/// A trait for things that can be approximately equal.
pub trait Epsilon {
    type RHS;
    type Precision;

    /// Return true if self and `other` differ no more than by a given amount.
    fn close(&self, other: Self::RHS, precision: Self::Precision) -> bool;

    /// Return true if self is close to zero.
    fn near_zero(&self, precision: Self::Precision) -> bool;
}

Bisection method for a single root is implemented as follows:

/// Configuration structure for the bisection method (one root version).
#[derive(Debug, Clone, Copy)]
pub struct OneRootBisectCfg<T> {
    /// The real root, if any, will be no further than this from the reported
    /// root.
    pub precision: T,
    /// A limit on the number of iterations to perform. Pass `None` if you
    /// don't want a limit.
    pub max_iters: Option<u32>
}

/// Find a root for a given function in a given interval, assuming there is
/// only one root there.
pub fn bisect_one<T, F>(config: OneRootBisectCfg<T>,
                        left: T,
                        right: T,
                        target: &F)
    -> Option<T>
    where T: Float + FromPrimitive + Signed,
          F: Fn(T) -> T
{
    let mut iter = 0;
    let mut left = left;
    let mut right = right;
    let mut left_val = target(left);
    let mut right_val = target(right);

    if left_val * right_val > T::zero() {
        return None;
    }

    let mut mid = (left + right) / T::from_i32(2).unwrap();
    let mut mid_val = target(mid);
    let max = config.max_iters;
    while right - left > config.precision && max.map_or(true, |m| iter < m) {
        if left_val * mid_val <= T::zero() {
            right = mid;
            right_val = mid_val;
        } else if mid_val * right_val <= T::zero() {
            left = mid;
            left_val = mid_val;
        } else {
            return None;
        }
        iter += 1;
        mid = (left + right) / T::from_i32(2).unwrap();
        mid_val = target(mid);
    }

    if abs(left_val) < abs(mid_val) {
        Some(left)
    } else if abs(right_val) < abs(mid_val) {
        Some(right)
    } else {
        Some(mid)
    }
}

And here's the version that tries to find all roots in the given interval. It's written as an iterator. I have an additional question regarding this: is it better to export just bisect_multi, or export both bisect_multi and MultiRootBisectState (perhaps with a name change)?

/// Configuration structure for the bisection method (multiple roots version).
#[derive(Debug, Clone, Copy)]
pub struct MultiRootBisectCfg<T> {
    /// Real roots will be no further than this from the reported roots.
    pub precision: T,
    /// A limit on the number of iterations to perform. Pass `None` if you
    /// don't want a limit.
    pub max_iters: Option<u32>,
    /// The requested interval will be split into this many chunks, and each
    /// chunk will be tried for a root.
    pub num_intervals: usize
}

#[derive(Debug, Clone, Copy)]
struct MultiRootBisectState<'a, T, F: 'a> {
    cfg: MultiRootBisectCfg<T>,
    left: T,
    right: T,
    target: &'a F,
    cur_interval: usize,
    last_root: Option<T>
}

pub fn bisect_multi<'a, T, F>(config: MultiRootBisectCfg<T>,
                              left: T,
                              right: T,
                              target: &'a F)
    -> impl Iterator<Item=T> + 'a
    where T: 'a + Float + FromPrimitive + Epsilon<RHS=T, Precision=T> + Signed,
          F: Fn(T) -> T
{
    MultiRootBisectState {
        cfg: config,
        left,
        right,
        target,
        cur_interval: 0,
        last_root: None
    }
}

impl<'a, T, F> Iterator for MultiRootBisectState<'a, T, F>
    where T: Float + FromPrimitive + Signed + Epsilon<RHS=T, Precision=T>,
          F: 'a + Fn(T) -> T
{
    type Item = T;

    fn next(&mut self) -> Option<T> {
        if self.cur_interval > self.cfg.num_intervals {
            return None
        }
        let num_ints = T::from_usize(self.cfg.num_intervals)
            .expect("Failed to convert the number of intervals into a float");
        let interval_width = (self.right - self.left) / num_ints;
        while self.cur_interval < self.cfg.num_intervals {
            let int = T::from_usize(self.cur_interval)
                .expect("Failed to convert an index into a float");
            let left = self.left + interval_width * int;
            let right = left + interval_width;
            let one_cfg = OneRootBisectCfg {
                precision: self.cfg.precision,
                max_iters: self.cfg.max_iters
            };
            let res = bisect_one(one_cfg, left, right, &self.target);
            self.cur_interval += 1;
            if let Some(root) = res {
                let two = T::from_i32(2).unwrap();
                let double_prec = two * self.cfg.precision;
                let mapper = |old: T| old.close(root, double_prec);
                let duplicate = self.last_root.map_or(false, mapper);
                if duplicate {
                    continue
                }
                self.last_root = Some(root);
                return Some(root)
            }
        }
        None
    }
}

Here's the Newton's method, single root version. It falls back to linearly interpolating the target function on the (possibly shrunk) interval and finding the root of that if the derivative of the target function gets too small, and fails altogether if it doesn't help either.

/// Configuration structure for the Newton's method (one root version).
#[derive(Debug, Clone, Copy)]
pub struct OneRootNewtonCfg<T> {
    /// The real root, if any, is most likely to be within this distance from
    /// the reported root, but this is not guaranteed.
    pub precision: T,
    /// A limit on the number of iterations to perform. Pass `None` if you
    /// don't want a limit.
    pub max_iters: Option<u32>
}

pub fn newton_one<T, F, D>(config: OneRootNewtonCfg<T>,
                           left: T,
                           right: T,
                           first_approx: T,
                           target: &F,
                           derivative: &D)
    -> Option<T>
    where T: Float + Epsilon<RHS=T, Precision=T>,
          F: Fn(T) -> T,
          D: Fn(T) -> T
{
    let mut left = left;
    let mut right = right;
    let mut left_val = target(left);
    let mut right_val = target(right);
    let mut root = first_approx;
    let mut prev_root = None;
    let mut iter = 0;
    while prev_root.map_or(true, |old| !root.close(old, config.precision))
        && config.max_iters.map_or(true, |max| iter < max) {
            iter += 1;
            if let Some(next) = next_newton_iter(config.precision,
                                                 left, 
                                                 right, 
                                                 root, 
                                                 target, 
                                                 derivative) {
                prev_root = Some(root);
                root = next;
            } else if let Some(fallback_root) 
                = linear_fallback(left, right, left_val, right_val) {
                    prev_root = Some(root);
                    root = fallback_root;
            } else {
                return None
            }
            let val_at_root = target(root);
            if left_val * val_at_root <= T::zero() {
                right = root;
                right_val = val_at_root;
            } else {
                left = root;
                left_val = val_at_root;
            }
    }
    Some(root)
}

fn next_newton_iter<T, F, D>(prec: T,
                             left: T,
                             right: T,
                             old: T,
                             target: &F,
                             derivative: &D)
    -> Option<T>
    where T: Float + Epsilon<RHS=T, Precision=T>,
          F: Fn(T) -> T,
          D: Fn(T) -> T
{
    let d = derivative(old);
    if d.near_zero(prec) {
        return None
    }
    let res = old - target(old) / d;
    if res < left {
        None
    } else if res > right {
        None
    } else {
        Some(res)
    }
}

fn linear_fallback<T: Float>(x1: T , x2: T, y1: T, y2: T) -> Option<T>
{
    let res = ((y2 - y1) * x1 - (x2 - x1) * y1) / (y2 - y1);
    if res < x1 {
        None
    } else if res > x2 {
        None
    } else {
        Some(res)
    }
}

Multiple roots version of Newton's method more or less the same as the multiple roots version of the bisection algorithm - split, find roots, weed out duplicates. Again, should I export the iterator struct, or is exporting newton_multi sufficient?

/// Configuration structure for the Newton's method (multiple roots version).
#[derive(Debug, Clone, Copy)]
pub struct MultiRootNewtonCfg<T> {
    /// Real root will most likely be no further that this from the reported
    /// roots, but it's not guaranteed.
    pub precision: T,
    /// A limit on the number of iterations to perform. Pass `None` if you 
    /// don't want to limit it. Note that this option governs maximum number of
    /// iterations on *each* chunk of the requested interval, not the number
    /// of iterations total.
    pub max_iters: Option<u32>,
    /// The requested interval will be split into this many chunks, and each
    /// chunk will be tested for a root separately.
    pub num_intervals: usize
}

#[derive(Debug, Clone, Copy)]
struct MultiRootNewtonState<'a, T, F: 'a, D: 'a> {
    cfg: MultiRootNewtonCfg<T>,
    left: T,
    right: T,
    target: &'a F,
    derivative: &'a D,
    last_root: Option<T>,
    cur_interval: usize
}

pub fn newton_multi<'a, T, F, D>(cfg: MultiRootNewtonCfg<T>,
                                 left: T,
                                 right: T,
                                 target: &'a F,
                                 derivative: &'a D)
    -> impl 'a + Iterator<Item=T> 
    where T: 'a + Float + FromPrimitive + Epsilon<RHS=T, Precision=T>,
          F: 'a + Fn(T) -> T,
          D: 'a + Fn(T) -> T
{
    MultiRootNewtonState {
        cfg,
        left,
        right,
        target,
        derivative,
        last_root: None,
        cur_interval: 0
    }
}

impl<'a, T, F, D> Iterator for MultiRootNewtonState<'a, T, F, D>
    where T: Float + FromPrimitive + Epsilon<RHS=T, Precision=T>,
          F: 'a + Fn(T) -> T,
          D: 'a + Fn(T) -> T
{
    type Item = T;

    fn next(&mut self) -> Option<T> {
        if self.cur_interval > self.cfg.num_intervals {
            return None
        }
        let intervals = T::from_usize(self.cfg.num_intervals)
            .expect("Failed to convert the number of intervals into a float");
        let interval_width = (self.right - self.left) / intervals;
        while self.cur_interval < self.cfg.num_intervals {
            let int = T::from_usize(self.cur_interval)
                .expect("Failed to convert an index into a float");
            let left = self.left + interval_width * int;
            let right = left + interval_width;
            let one_cfg = OneRootNewtonCfg {
                precision: self.cfg.precision,
                max_iters: self.cfg.max_iters
            };
            let two = T::from_i32(2).unwrap();
            let res = newton_one(one_cfg, left, right, (right - left) / two,
                self.target,
                self.derivative);
            self.cur_interval += 1;
            if let Some(root) = res {
                let double_prec = self.cfg.precision * two;
                let is_close = |prev: T| prev.close(root, double_prec);
                let duplicate = self.last_root.map_or(false, is_close);
                if duplicate {
                    continue
                }
                self.last_root = Some(root);
                return Some(root);
            }
        }
        None
    }
}
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.