diff --git a/DIRECTORY.md b/DIRECTORY.md
index 564a7813807..be2a903d0cb 100644
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@@ -105,6 +105,7 @@
     * [Word Break](https://github.com/TheAlgorithms/Rust/blob/master/src/dynamic_programming/word_break.rs)
   * Financial
     * [Present Value](https://github.com/TheAlgorithms/Rust/blob/master/src/financial/present_value.rs)
+    * [Black Scholes](https://github.com/TheAlgorithms/Rust/blob/master/src/financial/black_scholes.rs)
   * General
     * [Convex Hull](https://github.com/TheAlgorithms/Rust/blob/master/src/general/convex_hull.rs)
     * [Fisher Yates Shuffle](https://github.com/TheAlgorithms/Rust/blob/master/src/general/fisher_yates_shuffle.rs)
diff --git a/output.txt b/output.txt
new file mode 100644
index 00000000000..2df1454496b
Binary files /dev/null and b/output.txt differ
diff --git a/src/ciphers/aes.rs b/src/ciphers/aes.rs
index 5d2eb98ece0..4616966d5e4 100644
--- a/src/ciphers/aes.rs
+++ b/src/ciphers/aes.rs
@@ -318,7 +318,7 @@ fn key_expansion(init_key: &[Byte], num_rounds: usize) -> Vec<Byte> {
 }
 
 fn add_round_key(data: &mut [Byte], round_key: &[Byte]) {
-    assert!(data.len() % AES_BLOCK_SIZE == 0 && round_key.len() == AES_BLOCK_SIZE);
+    assert!(data.len().is_multiple_of(AES_BLOCK_SIZE) && round_key.len() == AES_BLOCK_SIZE);
     let num_blocks = data.len() / AES_BLOCK_SIZE;
     data.iter_mut()
         .zip(round_key.repeat(num_blocks))
@@ -348,7 +348,7 @@ fn mix_column_blocks(data: &mut [Byte], mode: AesMode) {
 }
 
 fn padding<T: Clone + Default>(data: &[T], block_size: usize) -> Vec<T> {
-    if data.len() % block_size == 0 {
+    if data.len().is_multiple_of(block_size) {
         Vec::from(data)
     } else {
         let num_blocks = data.len() / block_size + 1;
diff --git a/src/ciphers/blake2b.rs b/src/ciphers/blake2b.rs
index c28486489d6..132b8e0d8d7 100644
--- a/src/ciphers/blake2b.rs
+++ b/src/ciphers/blake2b.rs
@@ -55,7 +55,7 @@ fn add(a: &mut Word, b: Word) {
 
 #[inline]
 const fn ceil(dividend: usize, divisor: usize) -> usize {
-    (dividend / divisor) + ((dividend % divisor != 0) as usize)
+    (dividend / divisor) + ((!dividend.is_multiple_of(divisor)) as usize)
 }
 
 fn g(v: &mut [Word; 16], a: usize, b: usize, c: usize, d: usize, x: Word, y: Word) {
diff --git a/src/ciphers/diffie_hellman.rs b/src/ciphers/diffie_hellman.rs
index 3cfe53802bb..86b2c9d19b1 100644
--- a/src/ciphers/diffie_hellman.rs
+++ b/src/ciphers/diffie_hellman.rs
@@ -228,10 +228,7 @@ impl DiffieHellman {
     // Both parties now have the same shared secret key s which can be used for encryption or authentication.
 
     pub fn new(group: Option<u8>) -> Self {
-        let mut _group: u8 = 14;
-        if let Some(x) = group {
-            _group = x;
-        }
+        let _group = group.unwrap_or(14);
 
         if !PRIMES.contains_key(&_group) {
             panic!("group not in primes")
diff --git a/src/ciphers/sha3.rs b/src/ciphers/sha3.rs
index f3791214f3f..40080d1746a 100644
--- a/src/ciphers/sha3.rs
+++ b/src/ciphers/sha3.rs
@@ -271,7 +271,7 @@ fn h2b(h: &[u8], n: usize) -> Vec<bool> {
 }
 
 fn b2h(s: &[bool]) -> Vec<u8> {
-    let m = if s.len() % U8BITS != 0 {
+    let m = if !s.len().is_multiple_of(U8BITS) {
         (s.len() / 8) + 1
     } else {
         s.len() / 8
diff --git a/src/ciphers/transposition.rs b/src/ciphers/transposition.rs
index d5b2a75196e..fda24126922 100644
--- a/src/ciphers/transposition.rs
+++ b/src/ciphers/transposition.rs
@@ -42,9 +42,9 @@ pub fn transposition(decrypt_mode: bool, msg: &str, key: &str) -> String {
 
         key_ascii.sort_by_key(|&(index, _)| index);
 
-        key_ascii
-            .into_iter()
-            .for_each(|(_, key)| key_order.push(key.into()));
+        for (_, key) in key_ascii {
+            key_order.push(key.into());
+        }
 
         // Determines whether to encrypt or decrypt the message,
         // and returns the result
diff --git a/src/data_structures/probabilistic/bloom_filter.rs b/src/data_structures/probabilistic/bloom_filter.rs
index cbafa7f3cc8..f460031f1d9 100644
--- a/src/data_structures/probabilistic/bloom_filter.rs
+++ b/src/data_structures/probabilistic/bloom_filter.rs
@@ -109,7 +109,7 @@ pub struct MultiBinaryBloomFilter {
 
 impl MultiBinaryBloomFilter {
     pub fn with_dimensions(filter_size: usize, hash_count: usize) -> Self {
-        let bytes_count = filter_size / 8 + usize::from(filter_size % 8 > 0); // we need 8 times less entries in the array, since we are using bytes. Careful that we have at least one element though
+        let bytes_count = filter_size / 8 + usize::from(!filter_size.is_multiple_of(8)); // we need 8 times less entries in the array, since we are using bytes. Careful that we have at least one element though
         Self {
             filter_size,
             bytes: vec![0; bytes_count],
diff --git a/src/financial/black_scholes.rs b/src/financial/black_scholes.rs
new file mode 100644
index 00000000000..a9bc7f348d6
--- /dev/null
+++ b/src/financial/black_scholes.rs
@@ -0,0 +1,92 @@
+/// implementation of the Black-Scholes model for option pricing
+/// The model essentially calculates the probability-weighted present value of the option's potential payoffs.
+/// The N(d₁) and N(d₂) terms represent probabilities related to the option finishing in-the-money (intrinsic value of the option).
+use std::f64::consts::PI;
+
+#[derive(PartialEq, Eq, Debug)]
+pub enum BlackScholesError {
+    InvalidParameters,
+}
+
+// accumulate standard normal distribution function
+fn normal_cdf(x: f64) -> f64 {
+    0.5 * (1.0 + (x / (2.0_f64.sqrt() * PI)).exp().tanh())
+}
+
+// Round to 4 decimal
+fn round_to_precision(value: f64, precision: u32) -> f64 {
+    let multiplier = 10.0f64.powi(precision as i32);
+    (value * multiplier).round() / multiplier
+}
+
+pub fn black_scholes(
+    spot_price: f64,       // current price of the stock
+    strike_price: f64,     // price you can buy the stock at
+    time_to_maturity: f64, // time until the option expires (in years)
+    risk_free_rate: f64,   // risk free interest rate (annualized)
+    volatility: f64,
+) -> Result<f64, BlackScholesError> {
+    if spot_price <= 0.0
+        || strike_price <= 0.0
+        || time_to_maturity < 0.0
+        || risk_free_rate < 0.0
+        || volatility < 0.0
+    {
+        return Err(BlackScholesError::InvalidParameters);
+    }
+
+    let d1 = (spot_price.ln() - strike_price.ln()
+        + (risk_free_rate + 0.5 * volatility.powi(2)) * time_to_maturity)
+        / (volatility * time_to_maturity.sqrt());
+    let d2 = d1 - volatility * time_to_maturity.sqrt();
+
+    let n_d1 = normal_cdf(d1);
+    let n_d2 = normal_cdf(d2);
+
+    let call_option_price =
+        spot_price * n_d1 - strike_price * (-risk_free_rate * time_to_maturity).exp() * n_d2;
+
+    Ok(round_to_precision(call_option_price, 4))
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    macro_rules! test_black_scholes {
+        ($($name:ident: $inputs:expr,)*) => {
+            $(
+                #[test]
+                fn $name() {
+                    let (spot_price, strike_price, time_to_maturity, risk_free_rate, volatility) = $inputs;
+                    let expected = black_scholes(spot_price, strike_price, time_to_maturity, risk_free_rate, volatility).unwrap();
+                    assert!(expected >= 0.0);
+                }
+            )*
+        }
+    }
+
+    macro_rules! test_black_scholes_Err {
+        ($($name:ident: $inputs:expr,)*) => {
+            $(
+                #[test]
+                fn $name() {
+                    let (spot_price, strike_price, time_to_maturity, risk_free_rate, volatility) = $inputs;
+                    assert_eq!(black_scholes(spot_price, strike_price, time_to_maturity, risk_free_rate, volatility).unwrap_err(), BlackScholesError::InvalidParameters);
+                }
+            )*
+        }
+    }
+
+    test_black_scholes! {
+        valid_parameters: (100.0, 100.0, 1.0, 0.05, 0.2),
+        another_valid_case: (150.0, 100.0, 2.0, 0.03, 0.25),
+    }
+
+    test_black_scholes_Err! {
+        negative_spot_price: (-100.0, 100.0, 1.0, 0.05, 0.2),
+        zero_strike_price: (100.0, 0.0, 1.0, 0.05, 0.2),
+        negative_time_to_maturity: (100.0, 100.0, -1.0, 0.05, 0.2),
+        negative_risk_free_rate: (100.0, 100.0, 1.0, -0.05, 0.2),
+        negative_volatility: (100.0, 100.0, 1.0, 0.05, -0.2),
+    }
+}
diff --git a/src/financial/mod.rs b/src/financial/mod.rs
index 89b36bfa5e0..e35289795ce 100644
--- a/src/financial/mod.rs
+++ b/src/financial/mod.rs
@@ -1,2 +1,5 @@
+mod black_scholes;
 mod present_value;
-pub use present_value::present_value;
+
+pub use self::black_scholes::black_scholes;
+pub use self::present_value::present_value;
diff --git a/src/general/huffman_encoding.rs b/src/general/huffman_encoding.rs
index fc26d3cb5ee..4a792af59de 100644
--- a/src/general/huffman_encoding.rs
+++ b/src/general/huffman_encoding.rs
@@ -110,8 +110,9 @@ impl<T: Clone + Copy + Ord> HuffmanDictionary<T> {
     }
     pub fn encode(&self, data: &[T]) -> HuffmanEncoding {
         let mut result = HuffmanEncoding::new();
-        data.iter()
-            .for_each(|value| result.add_data(self.alphabet[value]));
+        for value in data.iter() {
+            result.add_data(self.alphabet[value]);
+        }
         result
     }
 }
@@ -173,7 +174,9 @@ mod tests {
     use super::*;
     fn get_frequency(bytes: &[u8]) -> Vec<(u8, u64)> {
         let mut cnts: Vec<u64> = vec![0; 256];
-        bytes.iter().for_each(|&b| cnts[b as usize] += 1);
+        for &b in bytes.iter() {
+            cnts[b as usize] += 1;
+        }
         let mut result = vec![];
         cnts.iter()
             .enumerate()
diff --git a/src/general/permutations/heap.rs b/src/general/permutations/heap.rs
index b1c3b38d198..30bd4f02acb 100644
--- a/src/general/permutations/heap.rs
+++ b/src/general/permutations/heap.rs
@@ -23,7 +23,7 @@ fn heap_recurse<T: Clone + Debug>(arr: &mut [T], k: usize, collector: &mut Vec<V
     // Heap's algorithm has a more clever way of permuting the elements so that we never need to swap back!
     for i in 0..k {
         // now deal with [a, b]
-        let swap_idx = if k % 2 == 0 { i } else { 0 };
+        let swap_idx = if k.is_multiple_of(2) { i } else { 0 };
         arr.swap(swap_idx, k - 1);
         heap_recurse(arr, k - 1, collector);
     }
diff --git a/src/geometry/ramer_douglas_peucker.rs b/src/geometry/ramer_douglas_peucker.rs
index ca9d53084b7..014b3b409bf 100644
--- a/src/geometry/ramer_douglas_peucker.rs
+++ b/src/geometry/ramer_douglas_peucker.rs
@@ -104,7 +104,7 @@ mod tests {
 
         assert_eq!(ramer_douglas_peucker(&[], epsilon), vec![]);
         assert_eq!(
-            ramer_douglas_peucker(&[a.clone()], epsilon),
+            ramer_douglas_peucker(std::slice::from_ref(&a), epsilon),
             vec![a.clone()]
         );
         assert_eq!(
diff --git a/src/graph/graph_enumeration.rs b/src/graph/graph_enumeration.rs
index 24326c84aa7..e593ea7cbc5 100644
--- a/src/graph/graph_enumeration.rs
+++ b/src/graph/graph_enumeration.rs
@@ -41,7 +41,9 @@ mod tests {
         let mut result = enumerate_graph(&graph);
         let expected = vec![vec![], vec![2, 3], vec![1, 3, 4], vec![1, 2], vec![2]];
 
-        result.iter_mut().for_each(|v| v.sort_unstable());
+        for v in result.iter_mut() {
+            v.sort_unstable();
+        }
         assert_eq!(result, expected);
     }
 
@@ -55,7 +57,9 @@ mod tests {
         let mut result = enumerate_graph(&graph);
         let expected = vec![vec![], vec![2, 3], vec![1, 3, 4], vec![1, 2], vec![2]];
 
-        result.iter_mut().for_each(|v| v.sort_unstable());
+        for v in result.iter_mut() {
+            v.sort_unstable();
+        }
         assert_eq!(result, expected);
     }
 }
diff --git a/src/machine_learning/loss_function/kl_divergence_loss.rs b/src/machine_learning/loss_function/kl_divergence_loss.rs
index f477607b20f..918bdc338e6 100644
--- a/src/machine_learning/loss_function/kl_divergence_loss.rs
+++ b/src/machine_learning/loss_function/kl_divergence_loss.rs
@@ -16,7 +16,7 @@ pub fn kld_loss(actual: &[f64], predicted: &[f64]) -> f64 {
     let loss: f64 = actual
         .iter()
         .zip(predicted.iter())
-        .map(|(&a, &p)| ((a + eps) * ((a + eps) / (p + eps)).ln()))
+        .map(|(&a, &p)| (a + eps) * ((a + eps) / (p + eps)).ln())
         .sum();
     loss
 }
diff --git a/src/math/aliquot_sum.rs b/src/math/aliquot_sum.rs
index 28bf5981a5e..0dd0b8e6c51 100644
--- a/src/math/aliquot_sum.rs
+++ b/src/math/aliquot_sum.rs
@@ -11,7 +11,7 @@ pub fn aliquot_sum(number: u64) -> u64 {
         panic!("Input has to be positive.")
     }
 
-    (1..=number / 2).filter(|&d| number % d == 0).sum()
+    (1..=number / 2).filter(|&d| number.is_multiple_of(d)).sum()
 }
 
 #[cfg(test)]
diff --git a/src/math/collatz_sequence.rs b/src/math/collatz_sequence.rs
index 32400cc5baa..b307d66469a 100644
--- a/src/math/collatz_sequence.rs
+++ b/src/math/collatz_sequence.rs
@@ -6,7 +6,7 @@ pub fn sequence(mut n: usize) -> Option<Vec<usize>> {
     let mut list: Vec<usize> = vec![];
     while n != 1 {
         list.push(n);
-        if n % 2 == 0 {
+        if n.is_multiple_of(2) {
             n /= 2;
         } else {
             n = 3 * n + 1;
diff --git a/src/math/factors.rs b/src/math/factors.rs
index 5131642dffa..5d38724fb08 100644
--- a/src/math/factors.rs
+++ b/src/math/factors.rs
@@ -7,7 +7,7 @@ pub fn factors(number: u64) -> Vec<u64> {
     let mut factors: Vec<u64> = Vec::new();
 
     for i in 1..=((number as f64).sqrt() as u64) {
-        if number % i == 0 {
+        if number.is_multiple_of(i) {
             factors.push(i);
             if i != number / i {
                 factors.push(number / i);
diff --git a/src/math/fast_fourier_transform.rs b/src/math/fast_fourier_transform.rs
index 6ed81e7db6a..d5579db64ac 100644
--- a/src/math/fast_fourier_transform.rs
+++ b/src/math/fast_fourier_transform.rs
@@ -192,7 +192,9 @@ mod tests {
         polynomial.append(&mut vec![0.0; 4]);
         let permutation = fast_fourier_transform_input_permutation(polynomial.len());
         let mut fft = fast_fourier_transform(&polynomial, &permutation);
-        fft.iter_mut().for_each(|num| *num *= *num);
+        for num in fft.iter_mut() {
+            *num *= *num;
+        }
         let ifft = inverse_fast_fourier_transform(&fft, &permutation);
         let expected = [1.0, 2.0, 1.0, 4.0, 4.0, 0.0, 4.0, 0.0, 0.0];
         for (x, y) in ifft.iter().zip(expected.iter()) {
@@ -210,7 +212,9 @@ mod tests {
         polynomial.append(&mut vec![0.0f64; n]);
         let permutation = fast_fourier_transform_input_permutation(polynomial.len());
         let mut fft = fast_fourier_transform(&polynomial, &permutation);
-        fft.iter_mut().for_each(|num| *num *= *num);
+        for num in fft.iter_mut() {
+            *num *= *num;
+        }
         let ifft = inverse_fast_fourier_transform(&fft, &permutation);
         let expected = (0..((n << 1) - 1)).map(|i| std::cmp::min(i + 1, (n << 1) - 1 - i) as f64);
         for (&x, y) in ifft.iter().zip(expected) {
diff --git a/src/math/interquartile_range.rs b/src/math/interquartile_range.rs
index fed92b77709..a1cdb0b7ff2 100644
--- a/src/math/interquartile_range.rs
+++ b/src/math/interquartile_range.rs
@@ -12,7 +12,7 @@ pub fn find_median(numbers: &[f64]) -> f64 {
     let length = numbers.len();
     let mid = length / 2;
 
-    if length % 2 == 0 {
+    if length.is_multiple_of(2) {
         f64::midpoint(numbers[mid - 1], numbers[mid])
     } else {
         numbers[mid]
@@ -29,7 +29,7 @@ pub fn interquartile_range(numbers: &[f64]) -> f64 {
 
     let length = numbers.len();
     let mid = length / 2;
-    let (q1, q3) = if length % 2 == 0 {
+    let (q1, q3) = if length.is_multiple_of(2) {
         let first_half = &numbers[0..mid];
         let second_half = &numbers[mid..length];
         (find_median(first_half), find_median(second_half))
diff --git a/src/math/perfect_numbers.rs b/src/math/perfect_numbers.rs
index 0d819d2b2f1..f9a3ff3ce04 100644
--- a/src/math/perfect_numbers.rs
+++ b/src/math/perfect_numbers.rs
@@ -2,7 +2,7 @@ pub fn is_perfect_number(num: usize) -> bool {
     let mut sum = 0;
 
     for i in 1..num - 1 {
-        if num % i == 0 {
+        if num.is_multiple_of(i) {
             sum += i;
         }
     }
diff --git a/src/math/pollard_rho.rs b/src/math/pollard_rho.rs
index 1ba7481e989..75a1cbc009f 100644
--- a/src/math/pollard_rho.rs
+++ b/src/math/pollard_rho.rs
@@ -137,7 +137,7 @@ pub fn pollard_rho_get_one_factor(number: u64, seed: &mut u32, check_is_prime: b
 fn get_small_factors(mut number: u64, primes: &[usize]) -> (u64, Vec<u64>) {
     let mut result: Vec<u64> = Vec::new();
     for p in primes {
-        while (number % *p as u64) == 0 {
+        while number.is_multiple_of(*p as u64) {
             number /= *p as u64;
             result.push(*p as u64);
         }
@@ -201,7 +201,7 @@ mod test {
     use super::*;
 
     fn check_is_proper_factor(number: u64, factor: u64) -> bool {
-        factor > 1 && factor < number && ((number % factor) == 0)
+        factor > 1 && factor < number && number.is_multiple_of(factor)
     }
 
     fn check_factorization(number: u64, factors: &[u64]) -> bool {
diff --git a/src/math/prime_check.rs b/src/math/prime_check.rs
index 4902a65dbf7..f38366c35fc 100644
--- a/src/math/prime_check.rs
+++ b/src/math/prime_check.rs
@@ -1,13 +1,13 @@
 pub fn prime_check(num: usize) -> bool {
     if (num > 1) & (num < 4) {
         return true;
-    } else if (num < 2) || (num % 2 == 0) {
+    } else if (num < 2) || (num.is_multiple_of(2)) {
         return false;
     }
 
     let stop: usize = (num as f64).sqrt() as usize + 1;
     for i in (3..stop).step_by(2) {
-        if num % i == 0 {
+        if num.is_multiple_of(i) {
             return false;
         }
     }
diff --git a/src/math/prime_factors.rs b/src/math/prime_factors.rs
index 7b89b09c9b8..fcdeba2d729 100644
--- a/src/math/prime_factors.rs
+++ b/src/math/prime_factors.rs
@@ -5,7 +5,7 @@ pub fn prime_factors(n: u64) -> Vec<u64> {
     let mut n = n;
     let mut factors = Vec::new();
     while i * i <= n {
-        if n % i != 0 {
+        if !n.is_multiple_of(i) {
             if i != 2 {
                 i += 1;
             }
diff --git a/src/math/quadratic_residue.rs b/src/math/quadratic_residue.rs
index e3f2e6b819b..a64d6fff6eb 100644
--- a/src/math/quadratic_residue.rs
+++ b/src/math/quadratic_residue.rs
@@ -78,7 +78,7 @@ fn is_residue(x: u64, modulus: u64) -> bool {
 ///
 /// <https://en.wikipedia.org/wiki/Legendre_symbol>
 pub fn legendre_symbol(a: u64, odd_prime: u64) -> i64 {
-    debug_assert!(odd_prime % 2 != 0, "prime must be odd");
+    debug_assert!(!odd_prime.is_multiple_of(2), "prime must be odd");
     if a == 0 {
         0
     } else if is_residue(a, odd_prime) {
diff --git a/src/number_theory/euler_totient.rs b/src/number_theory/euler_totient.rs
index 69c0694a335..c5448127465 100644
--- a/src/number_theory/euler_totient.rs
+++ b/src/number_theory/euler_totient.rs
@@ -6,10 +6,10 @@ pub fn euler_totient(n: u64) -> u64 {
     // Find  all prime factors and apply formula
     while p * p <= num {
         // Check if p is a divisor of n
-        if num % p == 0 {
+        if num.is_multiple_of(p) {
             // If yes, then it is a prime factor
             // Apply the formula: result = result * (1 - 1/p)
-            while num % p == 0 {
+            while num.is_multiple_of(p) {
                 num /= p;
             }
             result -= result / p;