diff --git a/integraal/src/parameters.rs b/integraal/src/parameters.rs index af6966c..35fa8a3 100644 --- a/integraal/src/parameters.rs +++ b/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, diff --git a/integraal/src/structure/implementations.rs b/integraal/src/structure/implementations.rs index 4fa8c58..123387e 100644 --- a/integraal/src/structure/implementations.rs +++ b/integraal/src/structure/implementations.rs @@ -145,6 +145,10 @@ fn values_explicit_arm( }) .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( }) .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( }) .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( 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( closure: impl Fn(X) -> X, domain: &DomainDescriptor, @@ -364,6 +400,31 @@ fn closure_uniform_arm( }) .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]); diff --git a/integraal/src/structure/tests.rs b/integraal/src/structure/tests.rs index 63e285d..3c15a0f 100644 --- a/integraal/src/structure/tests.rs +++ b/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 {