diff --git a/integraal/src/structure/implementations.rs b/integraal/src/structure/implementations.rs
index d390b23..524cd89 100644
--- a/integraal/src/structure/implementations.rs
+++ b/integraal/src/structure/implementations.rs
@@ -5,6 +5,8 @@
 use crate::{
     ComputeMethod, DomainDescriptor, FunctionDescriptor, Integraal, IntegraalError, Scalar,
 };
+#[cfg(feature = "montecarlo")]
+use rand::Rng;
 
 // ------ CONTENT
 
@@ -158,15 +160,41 @@ fn values_explicit_arm<X: Scalar>(
             ));
         }
         #[cfg(feature = "montecarlo")]
-        ComputeMethod::MonteCarlo { n_sample: _ } => {
-            todo!()
+        ComputeMethod::MonteCarlo { n_sample } => {
+            let (mut min, mut max) = (vals[0], vals[0]);
+            for v in vals {
+                min = min.min(*v);
+                max = max.max(*v);
+            }
+            let height: X = max - min;
+            let widths = args.windows(2).enumerate().map(|(i, slice)| {
+                let [a, b] = slice else { unreachable!() };
+                (vals[i].min(X::zero())..vals[i].max(X::zero()), (*b - *a))
+            });
+            let mut rng = rand::thread_rng();
+            let random_numbers: Vec<X> = (&mut rng)
+                .sample_iter(
+                    rand::distr::Uniform::new(min.to_f64().unwrap(), max.to_f64().unwrap())
+                        .unwrap(),
+                )
+                .take(*n_sample * (args.len() - 1))
+                .map(|s| X::from(s).unwrap())
+                .collect();
+            let total_in: X = random_numbers
+                .chunks_exact(*n_sample)
+                .zip(widths)
+                .map(|(samples, (range, width))| {
+                    X::from(samples.iter().filter(|s| range.contains(s)).count()).unwrap() * width
+                })
+                .sum();
+            height * total_in / X::from(*n_sample * (args.len() - 1)).unwrap()
         }
     };
 
     Ok(res)
 }
 
-#[allow(clippy::cast_possible_truncation)]
+#[allow(clippy::cast_possible_truncation, clippy::too_many_lines)]
 fn values_uniform_arm<X: Scalar>(
     vals: &[X],
     domain: &DomainDescriptor<X>,
@@ -278,8 +306,29 @@ fn values_uniform_arm<X: Scalar>(
             }
         }
         #[cfg(feature = "montecarlo")]
-        ComputeMethod::MonteCarlo { n_sample: _ } => {
-            todo!()
+        ComputeMethod::MonteCarlo { n_sample } => {
+            let (mut min, mut max) = (vals[0], vals[0]);
+            for v in vals {
+                min = min.min(*v);
+                max = max.max(*v);
+            }
+            let intervals = vals.iter().map(|v| v.min(X::zero())..v.max(X::zero()));
+            let volume: X = (max - min) * (X::from(*n_step).unwrap() * *step);
+            let mut rng = rand::thread_rng();
+            let random_numbers: Vec<X> = (&mut rng)
+                .sample_iter(
+                    rand::distr::Uniform::new(min.to_f64().unwrap(), max.to_f64().unwrap())
+                        .unwrap(),
+                )
+                .take(*n_sample * *n_step)
+                .map(|s| X::from(s).unwrap())
+                .collect();
+            let total_in: usize = random_numbers
+                .chunks_exact(*n_sample)
+                .zip(intervals)
+                .map(|(samples, range)| samples.iter().filter(|s| range.contains(s)).count())
+                .sum();
+            volume * X::from(total_in as f64 / (*n_sample * *n_step) as f64).unwrap()
         }
     };
 
@@ -345,8 +394,37 @@ fn closure_explicit_arm<X: Scalar>(
             ));
         }
         #[cfg(feature = "montecarlo")]
-        ComputeMethod::MonteCarlo { n_sample: _ } => {
-            todo!()
+        ComputeMethod::MonteCarlo { n_sample } => {
+            let (mut min, mut max) = (closure(args[0]), closure(args[0]));
+            for a in args {
+                min = min.min(closure(*a));
+                max = max.max(closure(*a));
+            }
+            let height: X = max - min;
+            let widths = args.windows(2).map(|slice| {
+                let [a, b] = slice else { unreachable!() };
+                (
+                    closure(*a).min(X::zero())..closure(*a).max(X::zero()),
+                    (*b - *a),
+                )
+            });
+            let mut rng = rand::thread_rng();
+            let random_numbers: Vec<X> = (&mut rng)
+                .sample_iter(
+                    rand::distr::Uniform::new(min.to_f64().unwrap(), max.to_f64().unwrap())
+                        .unwrap(),
+                )
+                .take(*n_sample * (args.len() - 1))
+                .map(|s| X::from(s).unwrap())
+                .collect();
+            let total_in: X = random_numbers
+                .chunks_exact(*n_sample)
+                .zip(widths)
+                .map(|(samples, (range, width))| {
+                    X::from(samples.iter().filter(|s| range.contains(s)).count()).unwrap() * width
+                })
+                .sum();
+            height * total_in / X::from(*n_sample * (args.len() - 1)).unwrap()
         }
     };
     Ok(res)
@@ -476,8 +554,33 @@ fn closure_uniform_arm<X: Scalar>(
             }
         }
         #[cfg(feature = "montecarlo")]
-        ComputeMethod::MonteCarlo { n_sample: _ } => {
-            todo!()
+        ComputeMethod::MonteCarlo { n_sample } => {
+            // FIXME: nuke this temp allocation
+            let vals: Vec<_> = (0..*n_step)
+                .map(|i| closure(*start + X::from(i).unwrap() * *step))
+                .collect();
+            let (mut min, mut max) = (vals[0], vals[0]);
+            for v in &vals {
+                min = min.min(*v);
+                max = max.max(*v);
+            }
+            let intervals = vals.iter().map(|v| v.min(X::zero())..v.max(X::zero()));
+            let volume: X = (max - min) * (X::from(*n_step).unwrap() * *step);
+            let mut rng = rand::thread_rng();
+            let random_numbers: Vec<X> = (&mut rng)
+                .sample_iter(
+                    rand::distr::Uniform::new(min.to_f64().unwrap(), max.to_f64().unwrap())
+                        .unwrap(),
+                )
+                .take(*n_sample * *n_step)
+                .map(|s| X::from(s).unwrap())
+                .collect();
+            let total_in: usize = random_numbers
+                .chunks_exact(*n_sample)
+                .zip(intervals)
+                .map(|(samples, range)| samples.iter().filter(|s| range.contains(s)).count())
+                .sum();
+            volume * X::from(total_in as f64 / (*n_sample * *n_step) as f64).unwrap()
         }
     };
 
diff --git a/integraal/src/tests/function_a.rs b/integraal/src/tests/function_a.rs
index 9ee9b34..46b4410 100644
--- a/integraal/src/tests/function_a.rs
+++ b/integraal/src/tests/function_a.rs
@@ -18,6 +18,7 @@ mod double {
     const SIMPSON_TOLERANCE: f64 = 1e-5;
     const BOOLE_TOLERANCE: f64 = 1e-5;
     const ROMBERG_TOLERANCE: f64 = 1e-5;
+    const MONTECARLO_TOLERANCE: f64 = 1e-5;
 
     all_tests!(
         f64,
@@ -50,6 +51,7 @@ mod simple {
     const SIMPSON_TOLERANCE: f32 = 1e-5;
     const BOOLE_TOLERANCE: f32 = 1e-5;
     const ROMBERG_TOLERANCE: f32 = 1e-5;
+    const MONTECARLO_TOLERANCE: f32 = 1e-5;
 
     all_tests!(
         f32,
diff --git a/integraal/src/tests/function_b.rs b/integraal/src/tests/function_b.rs
index 100b41d..f024b7d 100644
--- a/integraal/src/tests/function_b.rs
+++ b/integraal/src/tests/function_b.rs
@@ -20,6 +20,7 @@ mod double {
     const SIMPSON_TOLERANCE: f64 = 1e-5;
     const BOOLE_TOLERANCE: f64 = 1e-5;
     const ROMBERG_TOLERANCE: f64 = 1e-5;
+    const MONTECARLO_TOLERANCE: f64 = 1e-5;
 
     all_tests!(
         f64,
@@ -64,6 +65,7 @@ mod simple {
     const SIMPSON_TOLERANCE: f32 = 1e-5;
     const BOOLE_TOLERANCE: f32 = 1e-5;
     const ROMBERG_TOLERANCE: f32 = 1e-5;
+    const MONTECARLO_TOLERANCE: f32 = 1e-5;
 
     all_tests!(
         f32,
diff --git a/integraal/src/tests/function_c.rs b/integraal/src/tests/function_c.rs
index ec49720..1b4e017 100644
--- a/integraal/src/tests/function_c.rs
+++ b/integraal/src/tests/function_c.rs
@@ -20,6 +20,7 @@ mod double {
     const SIMPSON_TOLERANCE: f64 = 1e-5;
     const BOOLE_TOLERANCE: f64 = 1e-5;
     const ROMBERG_TOLERANCE: f64 = 1e-5;
+    const MONTECARLO_TOLERANCE: f64 = 1e-5;
 
     all_tests!(
         f64,
@@ -62,6 +63,7 @@ mod simple {
     const SIMPSON_TOLERANCE: f32 = 1e-5;
     const BOOLE_TOLERANCE: f32 = 1e-5;
     const ROMBERG_TOLERANCE: f32 = 1e-5;
+    const MONTECARLO_TOLERANCE: f32 = 1e-5;
 
     all_tests!(
         f32,
diff --git a/integraal/src/tests/function_d.rs b/integraal/src/tests/function_d.rs
index 53b44b3..26bbac4 100644
--- a/integraal/src/tests/function_d.rs
+++ b/integraal/src/tests/function_d.rs
@@ -19,6 +19,7 @@ mod double {
     const SIMPSON_TOLERANCE: f64 = 1e-5;
     const BOOLE_TOLERANCE: f64 = 1e-5;
     const ROMBERG_TOLERANCE: f64 = 1e-5;
+    const MONTECARLO_TOLERANCE: f64 = 1e-5;
 
     all_tests!(
         f64,
@@ -51,6 +52,7 @@ mod simple {
     const SIMPSON_TOLERANCE: f32 = 1e-5;
     const BOOLE_TOLERANCE: f32 = 1e-5;
     const ROMBERG_TOLERANCE: f32 = 1e-5;
+    const MONTECARLO_TOLERANCE: f32 = 1e-5;
 
     all_tests!(
         f32,
diff --git a/integraal/src/tests/function_e.rs b/integraal/src/tests/function_e.rs
index 285c897..3f2f878 100644
--- a/integraal/src/tests/function_e.rs
+++ b/integraal/src/tests/function_e.rs
@@ -19,6 +19,7 @@ mod double {
     const SIMPSON_TOLERANCE: f64 = 1e-5;
     const BOOLE_TOLERANCE: f64 = 1e-5;
     const ROMBERG_TOLERANCE: f64 = 1e-5;
+    const MONTECARLO_TOLERANCE: f64 = 1e-5;
 
     all_tests!(
         f64,
@@ -51,6 +52,7 @@ mod simple {
     const SIMPSON_TOLERANCE: f32 = 1e-5;
     const BOOLE_TOLERANCE: f32 = 1e-5;
     const ROMBERG_TOLERANCE: f32 = 1e-5;
+    const MONTECARLO_TOLERANCE: f32 = 1e-5;
 
     all_tests!(
         f32,
diff --git a/integraal/src/tests/mod.rs b/integraal/src/tests/mod.rs
index 3dc871c..3385581 100644
--- a/integraal/src/tests/mod.rs
+++ b/integraal/src/tests/mod.rs
@@ -360,6 +360,54 @@ macro_rules! all_tests {
                 ROMBERG_TOLERANCE
             );
         }
+
+        #[cfg(feature = "montecarlo")]
+        #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
+        mod montecarlo {
+            use super::*;
+
+            generate_test!(
+                $ft,
+                ClosureExplicit,
+                $dm,
+                $fnd_cls,
+                $dmd_xpl,
+                ComputeMethod::MonteCarlo { n_sample: 100 },
+                RES,
+                MONTECARLO_TOLERANCE
+            );
+
+            generate_test!(
+                $ft,
+                ClosureUniform,
+                $fnd_cls,
+                $dmd_uni,
+                ComputeMethod::MonteCarlo { n_sample: 100 },
+                RES,
+                MONTECARLO_TOLERANCE
+            );
+
+            generate_test!(
+                $ft,
+                ValuesExplicit,
+                $dm,
+                $fnd_val,
+                $dmd_xpl,
+                ComputeMethod::MonteCarlo { n_sample: 100 },
+                RES,
+                MONTECARLO_TOLERANCE
+            );
+
+            generate_test!(
+                $ft,
+                ValuesUniform,
+                $fnd_val,
+                $dmd_uni,
+                ComputeMethod::MonteCarlo { n_sample: 100 },
+                RES,
+                MONTECARLO_TOLERANCE
+            );
+        }
     };
 }