From 7e97059a3b604126a1545c65f3ad62fb201ca3d6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=A9r=C3=B4me?= <furious.jet@gmail.com>
Date: Sat, 25 Mar 2023 23:21:10 +0100
Subject: [PATCH] Added bisection

---
 src/main.rs               |  2 ++
 src/univariate_solvers.rs | 27 +++++++++++++++++++++++++++
 2 files changed, 29 insertions(+)

diff --git a/src/main.rs b/src/main.rs
index a4e4b3c..55ccc3f 100644
--- a/src/main.rs
+++ b/src/main.rs
@@ -28,7 +28,9 @@ fn main() {
     let x_mathematica = -3.26650043678562449167148755288;
     let x_newton = univariate_solvers::newton_solve(&(fct as fn(f64) -> f64), &(dfct as fn(f64) -> f64), x0, tol, max_iter);
     let x_newton_num: f64 = univariate_solvers::newton_solve_num(&(fct as fn(f64) -> f64), x0, tol, dx_num, max_iter);
+    let x_bisection : f64 = univariate_solvers::bisection_solve(&(fct as fn(f64) -> f64), -5.0, 1.0, tol).unwrap();
     println!("Mathematica           : x = {}", x_mathematica);
     println!("Newton's method       : x = {}", x_newton);
     println!("Newton's method (num) : x = {}", x_newton_num);
+    println!("Bisection             : x = {}", x_bisection);
 }
diff --git a/src/univariate_solvers.rs b/src/univariate_solvers.rs
index a6b03ce..189eca1 100644
--- a/src/univariate_solvers.rs
+++ b/src/univariate_solvers.rs
@@ -43,3 +43,30 @@ where F : Fn(f64) -> f64
     }, x0, tol, max_iter);
 }
 
+// --------------------------------------------------------------------
+// ------------------------ Bracketing methods ------------------------
+// --------------------------------------------------------------------
+
+pub fn bisection_solve<F>(f : F, mut a : f64, mut b : f64, tol : f64) -> Result<f64, &'static str>
+where F : Fn(f64) -> f64
+{
+    let mut c: f64;
+    let mut fa: f64 = f(a);
+    let mut fb: f64 = f(b);
+    let mut fc: f64;
+    let max_iter: u32 = (f64::log2((b-a)/tol)).ceil() as u32;
+    for _iter in 0..max_iter {
+        c = (a + b)/2.0;
+        fc = f(c);
+        if fa*fc < 0.0 {
+            b = c;
+            fb = fc;
+        } else if fb*fc < 0.0 {
+            a = c;
+            fa = fc;
+        } else {
+            return Result::Err("The interval does not backet the root.");
+        }
+    }
+    return Result::Ok((a + b)/2.0);
+}
\ No newline at end of file