feat: implement Boole's method (#31)

* add `Boole` enum variant

* implement boole's method

* add tests

integraal/src/parameters.rs

@@ -65,19 +65,21 @@ pub enum ComputeMethod {
     /// Rectangle method, using the right rule --
     /// [reference](https://en.wikipedia.org/wiki/Riemann_sum#Right_rule)
     RectangleRight,
-    /// Trapezoid method [reference](https://en.wikipedia.org/wiki/Trapezoidal_rule)
+    /// Trapezoid method -- [reference](https://en.wikipedia.org/wiki/Trapezoidal_rule)
     Trapezoid,
-    /// Simpson's third rule [reference](https://en.wikipedia.org/wiki/Simpson%27s_rule),
+    /// Simpson's third rule -- [reference](https://en.wikipedia.org/wiki/Simpson%27s_rule),
     /// [issue](https://github.com/imrn99/integraal/issues/23)
     SimpsonsThird,
-    /// Romberg's method [reference](https://en.wikipedia.org/wiki/Romberg%27s_method#Implementation)
+    /// Boole's method -- [reference](https://en.wikipedia.org/wiki/Boole%27s_rule#Composite_Boole's_Rule)
+    Boole,
+    /// Romberg's method -- [reference](https://en.wikipedia.org/wiki/Romberg%27s_method#Implementation)
     #[cfg(feature = "romberg")]
     Romberg {
         /// Maximum number of iteration done by the algorithm
         max_steps: usize,
     },
     #[cfg(feature = "montecarlo")]
-    /// Monte-Carlo method [reference](https://en.wikipedia.org/wiki/Monte_Carlo_integration)
+    /// Monte-Carlo method -- [reference](https://en.wikipedia.org/wiki/Monte_Carlo_integration)
     MonteCarlo {
         /// Number of random number sample per step computation.
         n_sample: usize,

integraal/src/structure/implementations.rs

@@ -145,6 +145,10 @@ fn values_explicit_arm<X: Scalar>(
                 })
                 .sum()
         }
+        ComputeMethod::Boole => {
+            // FIXME: replace by an error
+            unimplemented!("E: Romberg's method isn't implemented for non-uniform domains");
+        }
         #[cfg(feature = "romberg")] // FIXME: replace by an error
         ComputeMethod::Romberg { .. } => {
             unimplemented!("E: Romberg's method isn't implemented for non-uniform domains");
@@ -209,6 +213,30 @@ fn values_uniform_arm<X: Scalar>(
                     })
                     .sum()
         }
+        ComputeMethod::Boole => {
+            // FIXME: replace with an error
+            assert_eq!(
+                *n_step % 4,
+                0,
+                "E: domain should be divided into a multiple of 4 segments for Boole's method"
+            );
+            let c = X::from(2.0 / 45.0).unwrap() * *step;
+
+            let m1 = X::from(7.0).unwrap() * (vals[0] + vals[*n_step - 1]);
+
+            let c1 = X::from(14.0).unwrap();
+            let c2 = X::from(12.0).unwrap();
+            let c3 = X::from(32.0).unwrap();
+            let m2: X = (1..*n_step - 1)
+                .map(|id| match id % 4 {
+                    0 => c1 * vals[id],     // multiple of 4
+                    2 => c2 * vals[id],     // pair, non-multiple of 4
+                    1 | 3 => c3 * vals[id], // odd
+                    _ => unreachable!(),
+                })
+                .sum();
+            c * (m1 + m2)
+        }
         #[cfg(feature = "romberg")]
         ComputeMethod::Romberg { max_steps } => {
             let (mut r1, mut r2) = (vec![X::zero(); *max_steps], vec![X::zero(); *max_steps]);
@@ -298,6 +326,10 @@ fn closure_explicit_arm<X: Scalar>(
                 })
                 .sum()
         }
+        ComputeMethod::Boole => {
+            // FIXME: replace by an error
+            unimplemented!("E: Romberg's method isn't implemented for non-uniform domains");
+        }
         #[cfg(feature = "romberg")] // FIXME: replace by an error
         ComputeMethod::Romberg { .. } => {
             unimplemented!("E: Romberg's method isn't implemented for non-uniform domains");
@@ -310,7 +342,11 @@ fn closure_explicit_arm<X: Scalar>(
     Ok(res)
 }
 
-#[allow(clippy::unnecessary_wraps, clippy::cast_possible_truncation)]
+#[allow(
+    clippy::unnecessary_wraps,
+    clippy::cast_possible_truncation,
+    clippy::too_many_lines
+)]
 fn closure_uniform_arm<X: Scalar>(
     closure: impl Fn(X) -> X,
     domain: &DomainDescriptor<X>,
@@ -364,6 +400,31 @@ fn closure_uniform_arm<X: Scalar>(
                     })
                     .sum()
         }
+        ComputeMethod::Boole => {
+            // FIXME: replace with an error
+            assert_eq!(
+                *n_step % 4,
+                0,
+                "E: domain should be divided into a multiple of 4 segments for Boole's method"
+            );
+            let c = X::from(2.0 / 45.0).unwrap() * *step;
+
+            let m1 = X::from(7.0).unwrap()
+                * (closure(*start) + closure(*start + X::from(*n_step - 1).unwrap() * *step));
+
+            let c1 = X::from(14.0).unwrap();
+            let c2 = X::from(12.0).unwrap();
+            let c3 = X::from(32.0).unwrap();
+            let m2: X = (1..*n_step - 1)
+                .map(|id| match id % 4 {
+                    0 => c1 * closure(*start + X::from(id).unwrap() * *step), // multiple of 4
+                    2 => c2 * closure(*start + X::from(id).unwrap() * *step), // pair, non-multiple of 4
+                    1 | 3 => c3 * closure(*start + X::from(id).unwrap() * *step), // odd
+                    _ => unreachable!(),
+                })
+                .sum();
+            c * (m1 + m2)
+        }
         #[cfg(feature = "romberg")]
         ComputeMethod::Romberg { max_steps } => {
             let (mut r1, mut r2) = (vec![X::zero(); *max_steps], vec![X::zero(); *max_steps]);

integraal/src/structure/tests.rs

@@ -381,6 +381,39 @@ mod a_simpsons3rd {
     );
 }
 
+#[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
+mod a_boole {
+    use super::*;
+
+    generate_test!(
+        ClosureUniform,
+        FunctionDescriptor::Closure(Box::new(f64::sin)),
+        DomainDescriptor::Uniform {
+            start: 0.,
+            step: STEP,
+            n_step: (1000. * std::f64::consts::PI) as usize - 1,
+        },
+        ComputeMethod::Boole,
+        TRAPEZOID_TOLERANCE // FIXME: update tol
+    );
+
+    generate_test!(
+        ValuesUniform,
+        FunctionDescriptor::Values(
+            (0..(1000. * std::f64::consts::PI) as usize - 1)
+                .map(|step_id| (step_id as f64 * STEP).sin())
+                .collect()
+        ),
+        DomainDescriptor::Uniform {
+            start: 0.,
+            step: STEP,
+            n_step: (1000. * std::f64::consts::PI) as usize - 1,
+        },
+        ComputeMethod::Boole,
+        TRAPEZOID_TOLERANCE // FIXME: update tol
+    );
+}
+
 #[cfg(feature = "romberg")]
 #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
 mod a_romberg {