diff --git a/Cargo.toml b/Cargo.toml
index ba76fed..10be95c 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -19,7 +19,7 @@ path = "src/bin/main.rs"
 required-features = ["bin"]
 
 [dependencies]
-rand = { version = "0.8", optional = true }
+rand = { version = "0.9", optional = true }
 
 num-traits = { version = "0.2", default-features = false, features = ["std"], optional = true }
 
@@ -35,7 +35,7 @@ clap = { version = "4.0", optional = true, features = ["cargo"] }
 clap_autocomplete = { version = "0.4", optional = true }
 poloto = { version = "19", optional = true, default-features = false }
 tagu = "0.1.6"
-rand_xorshift = { version = "0.3.0", optional = true }
+rand_xorshift = { version = "0.4.0", optional = true }
 
 [features]
 default = ["bin", "pretty", "completion", "regression", "ols", "percentile-rand", "generic-impls", "binary_search_rng", "random_subset_regression"]
@@ -87,7 +87,7 @@ pretty = ["bin", "atty", "colored"]
 completion = ["clap_autocomplete"]
 
 [target.'cfg(all(target_arch = "wasm32", target_os = "unknown"))'.dependencies]
-getrandom = { version = "0.2", features = ["js"] }
+getrandom = { version = "0.3", features = ["wasm_js"] }
 
 # Build with `--profile production`
 [profile.production]
diff --git a/src/percentile.rs b/src/percentile.rs
index 7143017..65bf9f7 100644
--- a/src/percentile.rs
+++ b/src/percentile.rs
@@ -541,9 +541,9 @@ pub mod pivot_fn {
     #[cfg(feature = "percentile-rand")]
     #[inline]
     pub fn rand<T: Clone, S: SliceSubset<T> + ?Sized>() -> impl FnMut(&mut S) -> Cow<'_, T> {
-        let mut rng = rand::thread_rng();
+        let mut rng = rand::rng();
         move |slice| {
-            let idx = rng.sample(rand::distributions::Uniform::new(0_usize, slice.len()));
+            let idx = rng.random_range(0..slice.len());
             // UNWRAP: it's less than `slice.len`.
             // We assume `!slice.is_empty()`.
             Cow::Borrowed(slice.get(idx).unwrap())
@@ -614,9 +614,9 @@ pub mod cluster {
         #[cfg(feature = "percentile-rand")]
         #[inline]
         pub fn rand() -> impl FnMut(&ClusterList) -> f64 {
-            let mut rng = rand::thread_rng();
+            let mut rng = rand::rng();
             move |slice| {
-                let idx = rng.sample(rand::distributions::Uniform::new(0_usize, slice.len()));
+                let idx = rng.random_range(0..slice.len());
                 // Panic (index call): it's less than `slice.len`.
                 // We assume `!slice.is_empty()`.
                 *slice.index(idx)
diff --git a/src/regression.rs b/src/regression.rs
index 00f8ec3..6ba2ad3 100644
--- a/src/regression.rs
+++ b/src/regression.rs
@@ -2313,7 +2313,7 @@ pub mod theil_sen {
             let mut s1 = [0.0; 20];
             let mut s2 = [0.0; 20];
 
-            let mut rng = rand::thread_rng();
+            let mut rng = rand::rng();
             rng.fill(&mut s1);
             rng.fill(&mut s2);
 
@@ -3495,8 +3495,8 @@ pub mod binary_search {
                 // decrease randomness at the end
                 let progress = 1.0 - iter as f64 / self.iterations as f64;
                 // gen f32 since that takes less bytes
-                let rng_factor =
-                    1. + (2.0 * rng.gen::<f32>() as f64 - 1.) * self.randomness_factor * progress;
+                let rng_factor = 1.
+                    + (2.0 * rng.random::<f32>() as f64 - 1.) * self.randomness_factor * progress;
 
                 // for each variable to optimize
                 for i in 0..n {
@@ -3573,7 +3573,7 @@ pub mod binary_search {
                 impl $name for Options {
                     fn $method(&self, predictors: &[f64], outcomes: &[f64]) -> $ret {
                         use rand::SeedableRng;
-                        let mut rng = rand_xorshift::XorShiftRng::from_rng(rand::thread_rng()).unwrap();
+                        let mut rng = rand_xorshift::XorShiftRng::from_rng(&mut rand::rng());
 
                         #[cfg(feature = "random_subset_regression")]
                         if let Some(random_config) = &self.random_subset_regression {
@@ -3619,7 +3619,7 @@ pub mod binary_search {
                 impl $name for Options {
                     fn $method(&self, predictors: &[f64], outcomes: &[f64], max_frequency: f64) -> $ret {
                         use rand::SeedableRng;
-                        let mut rng = rand_xorshift::XorShiftRng::from_rng(rand::thread_rng()).unwrap();
+                        let mut rng = rand_xorshift::XorShiftRng::from_rng(&mut rand::rng());
 
                         #[cfg(feature = "random_subset_regression")]
                         if let Some(random_config) = &self.random_subset_regression {
@@ -3730,7 +3730,7 @@ pub mod binary_search {
             degree: usize,
         ) -> PolynomialCoefficients {
             use rand::SeedableRng;
-            let mut rng = rand_xorshift::XorShiftRng::from_rng(rand::thread_rng()).unwrap();
+            let mut rng = rand_xorshift::XorShiftRng::from_rng(&mut rand::rng());
 
             #[cfg(feature = "random_subset_regression")]
             if let Some(random_config) = &self.random_subset_regression {
@@ -3793,7 +3793,7 @@ pub mod binary_search {
         #[test]
         #[cfg(feature = "binary_search_rng")]
         fn two_variable_regression() {
-            let mut rng = rand::thread_rng();
+            let mut rng = rand::rng();
             let now = std::time::Instant::now();
             let x = [1.3, 4.7, 9.4];
             let y = [4., 5.3, 6.7];
@@ -3822,7 +3822,7 @@ pub mod binary_search {
         #[cfg(feature = "binary_search_rng")]
         fn second_degree_regression() {
             // init thread rng
-            let _rng = rand::thread_rng();
+            let _rng = rand::rng();
             let now = std::time::Instant::now();
             let x = [1.3, 4.7, 9.4];
             let y = [4., 5.3, 6.7];
@@ -3839,7 +3839,7 @@ pub mod binary_search {
             use rand::SeedableRng;
             // init thread rng
 
-            let mut rng = rand_xorshift::XorShiftRng::from_rng(rand::thread_rng()).unwrap();
+            let mut rng = rand_xorshift::XorShiftRng::from_rng(&mut rand::rng());
             let now = std::time::Instant::now();
             let coeffs = Options::default()
                 .max_precision()
@@ -3912,7 +3912,7 @@ pub mod random_subset_regression {
                 );
                 return None;
             }
-            let distribution = rand::distributions::Uniform::new(0, x.len());
+            let distribution = rand::distr::Uniform::new(0, x.len()).unwrap();
             let subsets = (0..config.subsets_count)
                 .map(|_| {
                     let mut new_x = Vec::with_capacity(config.subset_length);