diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 8a88dcc07..7340a0fd0 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,6 +1,6 @@
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.4.0
+    rev: v4.5.0
     hooks:
       - id: check-executables-have-shebangs
       - id: check-toml
diff --git a/DIRECTORY.md b/DIRECTORY.md
index a975b9264..015efb3c7 100644
--- a/DIRECTORY.md
+++ b/DIRECTORY.md
@@ -54,13 +54,12 @@
   * [Largest Pow Of Two Le Num](bit_manipulation/largest_pow_of_two_le_num.py)
   * [Missing Number](bit_manipulation/missing_number.py)
   * [Numbers Different Signs](bit_manipulation/numbers_different_signs.py)
+  * [Power Of 4](bit_manipulation/power_of_4.py)
   * [Reverse Bits](bit_manipulation/reverse_bits.py)
   * [Single Bit Manipulation Operations](bit_manipulation/single_bit_manipulation_operations.py)
 
 ## Blockchain
-  * [Chinese Remainder Theorem](blockchain/chinese_remainder_theorem.py)
   * [Diophantine Equation](blockchain/diophantine_equation.py)
-  * [Modular Division](blockchain/modular_division.py)
 
 ## Boolean Algebra
   * [And Gate](boolean_algebra/and_gate.py)
@@ -101,11 +100,13 @@
   * [Diffie Hellman](ciphers/diffie_hellman.py)
   * [Elgamal Key Generator](ciphers/elgamal_key_generator.py)
   * [Enigma Machine2](ciphers/enigma_machine2.py)
+  * [Fractionated Morse Cipher](ciphers/fractionated_morse_cipher.py)
   * [Hill Cipher](ciphers/hill_cipher.py)
   * [Mixed Keyword Cypher](ciphers/mixed_keyword_cypher.py)
   * [Mono Alphabetic Ciphers](ciphers/mono_alphabetic_ciphers.py)
   * [Morse Code](ciphers/morse_code.py)
   * [Onepad Cipher](ciphers/onepad_cipher.py)
+  * [Permutation Cipher](ciphers/permutation_cipher.py)
   * [Playfair Cipher](ciphers/playfair_cipher.py)
   * [Polybius](ciphers/polybius.py)
   * [Porta Cipher](ciphers/porta_cipher.py)
@@ -172,6 +173,7 @@
 
 ## Data Structures
   * Arrays
+    * [Equilibrium Index In Array](data_structures/arrays/equilibrium_index_in_array.py)
     * [Median Two Array](data_structures/arrays/median_two_array.py)
     * [Permutations](data_structures/arrays/permutations.py)
     * [Prefix Sum](data_structures/arrays/prefix_sum.py)
@@ -352,6 +354,7 @@
   * [Smith Waterman](dynamic_programming/smith_waterman.py)
   * [Subset Generation](dynamic_programming/subset_generation.py)
   * [Sum Of Subset](dynamic_programming/sum_of_subset.py)
+  * [Trapped Water](dynamic_programming/trapped_water.py)
   * [Tribonacci](dynamic_programming/tribonacci.py)
   * [Viterbi](dynamic_programming/viterbi.py)
   * [Word Break](dynamic_programming/word_break.py)
@@ -360,6 +363,7 @@
   * [Apparent Power](electronics/apparent_power.py)
   * [Builtin Voltage](electronics/builtin_voltage.py)
   * [Carrier Concentration](electronics/carrier_concentration.py)
+  * [Charging Capacitor](electronics/charging_capacitor.py)
   * [Circular Convolution](electronics/circular_convolution.py)
   * [Coulombs Law](electronics/coulombs_law.py)
   * [Electric Conductivity](electronics/electric_conductivity.py)
@@ -466,6 +470,8 @@
     * [Test Min Spanning Tree Prim](graphs/tests/test_min_spanning_tree_prim.py)
 
 ## Greedy Methods
+  * [Best Time To Buy And Sell Stock](greedy_methods/best_time_to_buy_and_sell_stock.py)
+  * [Fractional Cover Problem](greedy_methods/fractional_cover_problem.py)
   * [Fractional Knapsack](greedy_methods/fractional_knapsack.py)
   * [Fractional Knapsack 2](greedy_methods/fractional_knapsack_2.py)
   * [Gas Station](greedy_methods/gas_station.py)
@@ -524,6 +530,10 @@
   * Local Weighted Learning
     * [Local Weighted Learning](machine_learning/local_weighted_learning/local_weighted_learning.py)
   * [Logistic Regression](machine_learning/logistic_regression.py)
+  * Loss Functions
+    * [Binary Cross Entropy](machine_learning/loss_functions/binary_cross_entropy.py)
+    * [Huber Loss](machine_learning/loss_functions/huber_loss.py)
+    * [Mean Squared Error](machine_learning/loss_functions/mean_squared_error.py)
   * [Mfcc](machine_learning/mfcc.py)
   * [Multilayer Perceptron Classifier](machine_learning/multilayer_perceptron_classifier.py)
   * [Polynomial Regression](machine_learning/polynomial_regression.py)
@@ -556,7 +566,7 @@
   * [Bell Numbers](maths/bell_numbers.py)
   * [Binary Exp Mod](maths/binary_exp_mod.py)
   * [Binary Exponentiation](maths/binary_exponentiation.py)
-  * [Binary Exponentiation 3](maths/binary_exponentiation_3.py)
+  * [Binary Exponentiation 2](maths/binary_exponentiation_2.py)
   * [Binary Multiplication](maths/binary_multiplication.py)
   * [Binomial Coefficient](maths/binomial_coefficient.py)
   * [Binomial Distribution](maths/binomial_distribution.py)
@@ -564,7 +574,9 @@
   * [Carmichael Number](maths/carmichael_number.py)
   * [Catalan Number](maths/catalan_number.py)
   * [Ceil](maths/ceil.py)
+  * [Chebyshev Distance](maths/chebyshev_distance.py)
   * [Check Polygon](maths/check_polygon.py)
+  * [Chinese Remainder Theorem](maths/chinese_remainder_theorem.py)
   * [Chudnovsky Algorithm](maths/chudnovsky_algorithm.py)
   * [Collatz Sequence](maths/collatz_sequence.py)
   * [Combinations](maths/combinations.py)
@@ -588,10 +600,10 @@
   * [Find Min](maths/find_min.py)
   * [Floor](maths/floor.py)
   * [Gamma](maths/gamma.py)
-  * [Gamma Recursive](maths/gamma_recursive.py)
   * [Gaussian](maths/gaussian.py)
   * [Gaussian Error Linear Unit](maths/gaussian_error_linear_unit.py)
   * [Gcd Of N Numbers](maths/gcd_of_n_numbers.py)
+  * [Germain Primes](maths/germain_primes.py)
   * [Greatest Common Divisor](maths/greatest_common_divisor.py)
   * [Greedy Coin Change](maths/greedy_coin_change.py)
   * [Hamming Numbers](maths/hamming_numbers.py)
@@ -619,7 +631,9 @@
   * [Matrix Exponentiation](maths/matrix_exponentiation.py)
   * [Max Sum Sliding Window](maths/max_sum_sliding_window.py)
   * [Median Of Two Arrays](maths/median_of_two_arrays.py)
+  * [Minkowski Distance](maths/minkowski_distance.py)
   * [Mobius Function](maths/mobius_function.py)
+  * [Modular Division](maths/modular_division.py)
   * [Modular Exponential](maths/modular_exponential.py)
   * [Monte Carlo](maths/monte_carlo.py)
   * [Monte Carlo Dice](maths/monte_carlo_dice.py)
@@ -721,11 +735,16 @@
 ## Neural Network
   * [2 Hidden Layers Neural Network](neural_network/2_hidden_layers_neural_network.py)
   * Activation Functions
+    * [Binary Step](neural_network/activation_functions/binary_step.py)
     * [Exponential Linear Unit](neural_network/activation_functions/exponential_linear_unit.py)
     * [Leaky Rectified Linear Unit](neural_network/activation_functions/leaky_rectified_linear_unit.py)
+    * [Mish](neural_network/activation_functions/mish.py)
     * [Rectified Linear Unit](neural_network/activation_functions/rectified_linear_unit.py)
     * [Scaled Exponential Linear Unit](neural_network/activation_functions/scaled_exponential_linear_unit.py)
     * [Sigmoid Linear Unit](neural_network/activation_functions/sigmoid_linear_unit.py)
+    * [Soboleva Modified Hyperbolic Tangent](neural_network/activation_functions/soboleva_modified_hyperbolic_tangent.py)
+    * [Softplus](neural_network/activation_functions/softplus.py)
+    * [Squareplus](neural_network/activation_functions/squareplus.py)
   * [Back Propagation Neural Network](neural_network/back_propagation_neural_network.py)
   * [Convolution Neural Network](neural_network/convolution_neural_network.py)
   * [Perceptron](neural_network/perceptron.py)
@@ -748,6 +767,7 @@
   * [Linear Congruential Generator](other/linear_congruential_generator.py)
   * [Lru Cache](other/lru_cache.py)
   * [Magicdiamondpattern](other/magicdiamondpattern.py)
+  * [Majority Vote Algorithm](other/majority_vote_algorithm.py)
   * [Maximum Subsequence](other/maximum_subsequence.py)
   * [Nested Brackets](other/nested_brackets.py)
   * [Number Container System](other/number_container_system.py)
@@ -778,6 +798,7 @@
   * [Newtons Second Law Of Motion](physics/newtons_second_law_of_motion.py)
   * [Photoelectric Effect](physics/photoelectric_effect.py)
   * [Potential Energy](physics/potential_energy.py)
+  * [Reynolds Number](physics/reynolds_number.py)
   * [Rms Speed Of Molecule](physics/rms_speed_of_molecule.py)
   * [Shear Stress](physics/shear_stress.py)
   * [Speed Of Sound](physics/speed_of_sound.py)
@@ -1100,6 +1121,7 @@
   * [Interpolation Search](searches/interpolation_search.py)
   * [Jump Search](searches/jump_search.py)
   * [Linear Search](searches/linear_search.py)
+  * [Median Of Medians](searches/median_of_medians.py)
   * [Quick Select](searches/quick_select.py)
   * [Sentinel Linear Search](searches/sentinel_linear_search.py)
   * [Simple Binary Search](searches/simple_binary_search.py)
@@ -1196,11 +1218,11 @@
   * [Rabin Karp](strings/rabin_karp.py)
   * [Remove Duplicate](strings/remove_duplicate.py)
   * [Reverse Letters](strings/reverse_letters.py)
-  * [Reverse Long Words](strings/reverse_long_words.py)
   * [Reverse Words](strings/reverse_words.py)
   * [Snake Case To Camel Pascal Case](strings/snake_case_to_camel_pascal_case.py)
   * [Split](strings/split.py)
   * [String Switch Case](strings/string_switch_case.py)
+  * [Strip](strings/strip.py)
   * [Text Justification](strings/text_justification.py)
   * [Top K Frequent Words](strings/top_k_frequent_words.py)
   * [Upper](strings/upper.py)
diff --git a/arithmetic_analysis/gaussian_elimination.py b/arithmetic_analysis/gaussian_elimination.py
index 13f509a4f..a1a35131b 100644
--- a/arithmetic_analysis/gaussian_elimination.py
+++ b/arithmetic_analysis/gaussian_elimination.py
@@ -34,7 +34,7 @@ def retroactive_resolution(
     x: NDArray[float64] = np.zeros((rows, 1), dtype=float)
     for row in reversed(range(rows)):
         total = np.dot(coefficients[row, row + 1 :], x[row + 1 :])
-        x[row, 0] = (vector[row] - total) / coefficients[row, row]
+        x[row, 0] = (vector[row][0] - total[0]) / coefficients[row, row]
 
     return x
 
diff --git a/backtracking/word_search.py b/backtracking/word_search.py
index c9d52012b..8a9b2f1b5 100644
--- a/backtracking/word_search.py
+++ b/backtracking/word_search.py
@@ -98,13 +98,7 @@ def word_exists(board: list[list[str]], word: str) -> bool:
     False
     >>> word_exists([["A"]], "A")
     True
-    >>> word_exists([["A","A","A","A","A","A"],
-    ...              ["A","A","A","A","A","A"],
-    ...              ["A","A","A","A","A","A"],
-    ...              ["A","A","A","A","A","A"],
-    ...              ["A","A","A","A","A","B"],
-    ...              ["A","A","A","A","B","A"]],
-    ...             "AAAAAAAAAAAAABB")
+    >>> word_exists([["B", "A", "A"], ["A", "A", "A"], ["A", "B", "A"]], "ABB")
     False
     >>> word_exists([["A"]], 123)
     Traceback (most recent call last):
diff --git a/bit_manipulation/power_of_4.py b/bit_manipulation/power_of_4.py
new file mode 100644
index 000000000..09e6e2862
--- /dev/null
+++ b/bit_manipulation/power_of_4.py
@@ -0,0 +1,67 @@
+"""
+
+Task:
+Given a positive int number. Return True if this number is power of 4
+or False otherwise.
+
+Implementation notes: Use bit manipulation.
+For example if the number is the power of 2 it's bits representation:
+n     = 0..100..00
+n - 1 = 0..011..11
+
+n & (n - 1) - no intersections = 0
+If the number is a power of 4 then it should be a power of 2
+and the set bit should be at an odd position.
+"""
+
+
+def power_of_4(number: int) -> bool:
+    """
+    Return True if this number is power of 4 or False otherwise.
+
+    >>> power_of_4(0)
+    Traceback (most recent call last):
+        ...
+    ValueError: number must be positive
+    >>> power_of_4(1)
+    True
+    >>> power_of_4(2)
+    False
+    >>> power_of_4(4)
+    True
+    >>> power_of_4(6)
+    False
+    >>> power_of_4(8)
+    False
+    >>> power_of_4(17)
+    False
+    >>> power_of_4(64)
+    True
+    >>> power_of_4(-1)
+    Traceback (most recent call last):
+        ...
+    ValueError: number must be positive
+    >>> power_of_4(1.2)
+    Traceback (most recent call last):
+        ...
+    TypeError: number must be an integer
+
+    """
+    if not isinstance(number, int):
+        raise TypeError("number must be an integer")
+    if number <= 0:
+        raise ValueError("number must be positive")
+    if number & (number - 1) == 0:
+        c = 0
+        while number:
+            c += 1
+            number >>= 1
+        return c % 2 == 1
+    else:
+        return False
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/blockchain/diophantine_equation.py b/blockchain/diophantine_equation.py
index 22b0cad75..7110d9023 100644
--- a/blockchain/diophantine_equation.py
+++ b/blockchain/diophantine_equation.py
@@ -1,11 +1,13 @@
 from __future__ import annotations
 
+from maths.greatest_common_divisor import greatest_common_divisor
+
 
 def diophantine(a: int, b: int, c: int) -> tuple[float, float]:
     """
     Diophantine Equation : Given integers a,b,c ( at least one of a and b != 0), the
     diophantine equation a*x + b*y = c has a solution (where x and y are integers)
-    iff gcd(a,b) divides c.
+    iff greatest_common_divisor(a,b) divides c.
 
     GCD ( Greatest Common Divisor ) or HCF ( Highest Common Factor )
 
@@ -22,7 +24,7 @@ def diophantine(a: int, b: int, c: int) -> tuple[float, float]:
 
     assert (
         c % greatest_common_divisor(a, b) == 0
-    )  # greatest_common_divisor(a,b) function implemented below
+    )  # greatest_common_divisor(a,b) is in maths directory
     (d, x, y) = extended_gcd(a, b)  # extended_gcd(a,b) function implemented below
     r = c / d
     return (r * x, r * y)
@@ -69,32 +71,6 @@ def diophantine_all_soln(a: int, b: int, c: int, n: int = 2) -> None:
         print(x, y)
 
 
-def greatest_common_divisor(a: int, b: int) -> int:
-    """
-    Euclid's Lemma :  d divides a and b, if and only if d divides a-b and b
-
-    Euclid's Algorithm
-
-    >>> greatest_common_divisor(7,5)
-    1
-
-    Note : In number theory, two integers a and b are said to be relatively prime,
-           mutually prime, or co-prime if the only positive integer (factor) that
-           divides both of them is 1  i.e., gcd(a,b) = 1.
-
-    >>> greatest_common_divisor(121, 11)
-    11
-
-    """
-    if a < b:
-        a, b = b, a
-
-    while a % b != 0:
-        a, b = b, a % b
-
-    return b
-
-
 def extended_gcd(a: int, b: int) -> tuple[int, int, int]:
     """
     Extended Euclid's Algorithm : If d divides a and b and d = a*x + b*y for integers
diff --git a/ciphers/affine_cipher.py b/ciphers/affine_cipher.py
index cd1e33b88..10d16367c 100644
--- a/ciphers/affine_cipher.py
+++ b/ciphers/affine_cipher.py
@@ -1,6 +1,8 @@
 import random
 import sys
 
+from maths.greatest_common_divisor import gcd_by_iterative
+
 from . import cryptomath_module as cryptomath
 
 SYMBOLS = (
@@ -26,7 +28,7 @@ def check_keys(key_a: int, key_b: int, mode: str) -> None:
             "Key A must be greater than 0 and key B must "
             f"be between 0 and {len(SYMBOLS) - 1}."
         )
-    if cryptomath.gcd(key_a, len(SYMBOLS)) != 1:
+    if gcd_by_iterative(key_a, len(SYMBOLS)) != 1:
         sys.exit(
             f"Key A {key_a} and the symbol set size {len(SYMBOLS)} "
             "are not relatively prime. Choose a different key."
@@ -76,7 +78,7 @@ def get_random_key() -> int:
     while True:
         key_b = random.randint(2, len(SYMBOLS))
         key_b = random.randint(2, len(SYMBOLS))
-        if cryptomath.gcd(key_b, len(SYMBOLS)) == 1 and key_b % len(SYMBOLS) != 0:
+        if gcd_by_iterative(key_b, len(SYMBOLS)) == 1 and key_b % len(SYMBOLS) != 0:
             return key_b * len(SYMBOLS) + key_b
 
 
diff --git a/ciphers/cryptomath_module.py b/ciphers/cryptomath_module.py
index 6f15f7b73..02e94e4b9 100644
--- a/ciphers/cryptomath_module.py
+++ b/ciphers/cryptomath_module.py
@@ -1,11 +1,8 @@
-def gcd(a: int, b: int) -> int:
-    while a != 0:
-        a, b = b % a, a
-    return b
+from maths.greatest_common_divisor import gcd_by_iterative
 
 
 def find_mod_inverse(a: int, m: int) -> int:
-    if gcd(a, m) != 1:
+    if gcd_by_iterative(a, m) != 1:
         msg = f"mod inverse of {a!r} and {m!r} does not exist"
         raise ValueError(msg)
     u1, u2, u3 = 1, 0, a
diff --git a/ciphers/diffie.py b/ciphers/diffie.py
index 4ff90be00..1e1e86899 100644
--- a/ciphers/diffie.py
+++ b/ciphers/diffie.py
@@ -1,11 +1,28 @@
 from __future__ import annotations
 
 
-def find_primitive(n: int) -> int | None:
-    for r in range(1, n):
+def find_primitive(modulus: int) -> int | None:
+    """
+    Find a primitive root modulo modulus, if one exists.
+
+    Args:
+        modulus : The modulus for which to find a primitive root.
+
+    Returns:
+        The primitive root if one exists, or None if there is none.
+
+    Examples:
+    >>> find_primitive(7)  # Modulo 7 has primitive root 3
+    3
+    >>> find_primitive(11)  # Modulo 11 has primitive root 2
+    2
+    >>> find_primitive(8) == None # Modulo 8 has no primitive root
+    True
+    """
+    for r in range(1, modulus):
         li = []
-        for x in range(n - 1):
-            val = pow(r, x, n)
+        for x in range(modulus - 1):
+            val = pow(r, x, modulus)
             if val in li:
                 break
             li.append(val)
@@ -15,18 +32,22 @@ def find_primitive(n: int) -> int | None:
 
 
 if __name__ == "__main__":
-    q = int(input("Enter a prime number q: "))
-    a = find_primitive(q)
-    if a is None:
-        print(f"Cannot find the primitive for the value: {a!r}")
+    import doctest
+
+    doctest.testmod()
+
+    prime = int(input("Enter a prime number q: "))
+    primitive_root = find_primitive(prime)
+    if primitive_root is None:
+        print(f"Cannot find the primitive for the value: {primitive_root!r}")
     else:
         a_private = int(input("Enter private key of A: "))
-        a_public = pow(a, a_private, q)
+        a_public = pow(primitive_root, a_private, prime)
         b_private = int(input("Enter private key of B: "))
-        b_public = pow(a, b_private, q)
+        b_public = pow(primitive_root, b_private, prime)
 
-        a_secret = pow(b_public, a_private, q)
-        b_secret = pow(a_public, b_private, q)
+        a_secret = pow(b_public, a_private, prime)
+        b_secret = pow(a_public, b_private, prime)
 
         print("The key value generated by A is: ", a_secret)
         print("The key value generated by B is: ", b_secret)
diff --git a/ciphers/hill_cipher.py b/ciphers/hill_cipher.py
index b4424e822..1201fda90 100644
--- a/ciphers/hill_cipher.py
+++ b/ciphers/hill_cipher.py
@@ -39,19 +39,7 @@ import string
 
 import numpy
 
-
-def greatest_common_divisor(a: int, b: int) -> int:
-    """
-    >>> greatest_common_divisor(4, 8)
-    4
-    >>> greatest_common_divisor(8, 4)
-    4
-    >>> greatest_common_divisor(4, 7)
-    1
-    >>> greatest_common_divisor(0, 10)
-    10
-    """
-    return b if a == 0 else greatest_common_divisor(b % a, a)
+from maths.greatest_common_divisor import greatest_common_divisor
 
 
 class HillCipher:
diff --git a/ciphers/permutation_cipher.py b/ciphers/permutation_cipher.py
new file mode 100644
index 000000000..c3f3fd1f7
--- /dev/null
+++ b/ciphers/permutation_cipher.py
@@ -0,0 +1,142 @@
+"""
+The permutation cipher, also called the transposition cipher, is a simple encryption
+technique that rearranges the characters in a message based on a secret key. It
+divides the message into blocks and applies a permutation to the characters within
+each block according to the key. The key is a sequence of unique integers that
+determine the order of character rearrangement.
+
+For more info: https://www.nku.edu/~christensen/1402%20permutation%20ciphers.pdf
+"""
+import random
+
+
+def generate_valid_block_size(message_length: int) -> int:
+    """
+    Generate a valid block size that is a factor of the message length.
+
+    Args:
+        message_length (int): The length of the message.
+
+    Returns:
+        int: A valid block size.
+
+    Example:
+        >>> random.seed(1)
+        >>> generate_valid_block_size(12)
+        3
+    """
+    block_sizes = [
+        block_size
+        for block_size in range(2, message_length + 1)
+        if message_length % block_size == 0
+    ]
+    return random.choice(block_sizes)
+
+
+def generate_permutation_key(block_size: int) -> list[int]:
+    """
+    Generate a random permutation key of a specified block size.
+
+    Args:
+        block_size (int): The size of each permutation block.
+
+    Returns:
+        list[int]: A list containing a random permutation of digits.
+
+    Example:
+        >>> random.seed(0)
+        >>> generate_permutation_key(4)
+        [2, 0, 1, 3]
+    """
+    digits = list(range(block_size))
+    random.shuffle(digits)
+    return digits
+
+
+def encrypt(
+    message: str, key: list[int] | None = None, block_size: int | None = None
+) -> tuple[str, list[int]]:
+    """
+    Encrypt a message using a permutation cipher with block rearrangement using a key.
+
+    Args:
+        message (str): The plaintext message to be encrypted.
+        key (list[int]): The permutation key for decryption.
+        block_size (int): The size of each permutation block.
+
+    Returns:
+        tuple: A tuple containing the encrypted message and the encryption key.
+
+    Example:
+        >>> encrypted_message, key = encrypt("HELLO WORLD")
+        >>> decrypted_message = decrypt(encrypted_message, key)
+        >>> decrypted_message
+        'HELLO WORLD'
+    """
+    message = message.upper()
+    message_length = len(message)
+
+    if key is None or block_size is None:
+        block_size = generate_valid_block_size(message_length)
+        key = generate_permutation_key(block_size)
+
+    encrypted_message = ""
+
+    for i in range(0, message_length, block_size):
+        block = message[i : i + block_size]
+        rearranged_block = [block[digit] for digit in key]
+        encrypted_message += "".join(rearranged_block)
+
+    return encrypted_message, key
+
+
+def decrypt(encrypted_message: str, key: list[int]) -> str:
+    """
+    Decrypt an encrypted message using a permutation cipher with block rearrangement.
+
+    Args:
+        encrypted_message (str): The encrypted message.
+        key (list[int]): The permutation key for decryption.
+
+    Returns:
+        str: The decrypted plaintext message.
+
+    Example:
+        >>> encrypted_message, key = encrypt("HELLO WORLD")
+        >>> decrypted_message = decrypt(encrypted_message, key)
+        >>> decrypted_message
+        'HELLO WORLD'
+    """
+    key_length = len(key)
+    decrypted_message = ""
+
+    for i in range(0, len(encrypted_message), key_length):
+        block = encrypted_message[i : i + key_length]
+        original_block = [""] * key_length
+        for j, digit in enumerate(key):
+            original_block[digit] = block[j]
+        decrypted_message += "".join(original_block)
+
+    return decrypted_message
+
+
+def main() -> None:
+    """
+    Driver function to pass message to get encrypted, then decrypted.
+
+    Example:
+    >>> main()
+    Decrypted message: HELLO WORLD
+    """
+    message = "HELLO WORLD"
+    encrypted_message, key = encrypt(message)
+
+    decrypted_message = decrypt(encrypted_message, key)
+    print(f"Decrypted message: {decrypted_message}")
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
+    main()
diff --git a/ciphers/rsa_key_generator.py b/ciphers/rsa_key_generator.py
index eedc73368..44970e8cb 100644
--- a/ciphers/rsa_key_generator.py
+++ b/ciphers/rsa_key_generator.py
@@ -2,6 +2,8 @@ import os
 import random
 import sys
 
+from maths.greatest_common_divisor import gcd_by_iterative
+
 from . import cryptomath_module, rabin_miller
 
 
@@ -27,7 +29,7 @@ def generate_key(key_size: int) -> tuple[tuple[int, int], tuple[int, int]]:
     # Generate e that is relatively prime to (p - 1) * (q - 1)
     while True:
         e = random.randrange(2 ** (key_size - 1), 2 ** (key_size))
-        if cryptomath_module.gcd(e, (p - 1) * (q - 1)) == 1:
+        if gcd_by_iterative(e, (p - 1) * (q - 1)) == 1:
             break
 
     # Calculate d that is mod inverse of e
diff --git a/data_structures/arrays/equilibrium_index_in_array.py b/data_structures/arrays/equilibrium_index_in_array.py
new file mode 100644
index 000000000..4099896d2
--- /dev/null
+++ b/data_structures/arrays/equilibrium_index_in_array.py
@@ -0,0 +1,59 @@
+"""
+Find the Equilibrium Index of an Array.
+Reference: https://www.geeksforgeeks.org/equilibrium-index-of-an-array/
+
+Python doctests can be run with the following command:
+python -m doctest -v equilibrium_index.py
+
+Given a sequence arr[] of size n, this function returns
+an equilibrium index (if any) or -1 if no equilibrium index exists.
+
+The equilibrium index of an array is an index such that the sum of
+elements at lower indexes is equal to the sum of elements at higher indexes.
+
+
+
+Example Input:
+arr = [-7, 1, 5, 2, -4, 3, 0]
+Output: 3
+
+"""
+
+
+def equilibrium_index(arr: list[int], size: int) -> int:
+    """
+    Find the equilibrium index of an array.
+
+    Args:
+        arr : The input array of integers.
+        size : The size of the array.
+
+    Returns:
+        int: The equilibrium index or -1 if no equilibrium index exists.
+
+    Examples:
+        >>> equilibrium_index([-7, 1, 5, 2, -4, 3, 0], 7)
+        3
+        >>> equilibrium_index([1, 2, 3, 4, 5], 5)
+        -1
+        >>> equilibrium_index([1, 1, 1, 1, 1], 5)
+        2
+        >>> equilibrium_index([2, 4, 6, 8, 10, 3], 6)
+        -1
+    """
+    total_sum = sum(arr)
+    left_sum = 0
+
+    for i in range(size):
+        total_sum -= arr[i]
+        if left_sum == total_sum:
+            return i
+        left_sum += arr[i]
+
+    return -1
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/data_structures/linked_list/is_palindrome.py b/data_structures/linked_list/is_palindrome.py
index d540fb69f..f949d9a2f 100644
--- a/data_structures/linked_list/is_palindrome.py
+++ b/data_structures/linked_list/is_palindrome.py
@@ -1,65 +1,169 @@
-def is_palindrome(head):
+from __future__ import annotations
+
+from dataclasses import dataclass
+
+
+@dataclass
+class ListNode:
+    val: int = 0
+    next_node: ListNode | None = None
+
+
+def is_palindrome(head: ListNode | None) -> bool:
+    """
+    Check if a linked list is a palindrome.
+
+    Args:
+        head: The head of the linked list.
+
+    Returns:
+        bool: True if the linked list is a palindrome, False otherwise.
+
+    Examples:
+        >>> is_palindrome(None)
+        True
+
+        >>> is_palindrome(ListNode(1))
+        True
+
+        >>> is_palindrome(ListNode(1, ListNode(2)))
+        False
+
+        >>> is_palindrome(ListNode(1, ListNode(2, ListNode(1))))
+        True
+
+        >>> is_palindrome(ListNode(1, ListNode(2, ListNode(2, ListNode(1)))))
+        True
+    """
     if not head:
         return True
     # split the list to two parts
-    fast, slow = head.next, head
-    while fast and fast.next:
-        fast = fast.next.next
-        slow = slow.next
-    second = slow.next
-    slow.next = None  # Don't forget here! But forget still works!
+    fast: ListNode | None = head.next_node
+    slow: ListNode | None = head
+    while fast and fast.next_node:
+        fast = fast.next_node.next_node
+        slow = slow.next_node if slow else None
+    if slow:
+        # slow will always be defined,
+        # adding this check to resolve mypy static check
+        second = slow.next_node
+        slow.next_node = None  # Don't forget here! But forget still works!
     # reverse the second part
-    node = None
+    node: ListNode | None = None
     while second:
-        nxt = second.next
-        second.next = node
+        nxt = second.next_node
+        second.next_node = node
         node = second
         second = nxt
     # compare two parts
     # second part has the same or one less node
-    while node:
+    while node and head:
         if node.val != head.val:
             return False
-        node = node.next
-        head = head.next
+        node = node.next_node
+        head = head.next_node
     return True
 
 
-def is_palindrome_stack(head):
-    if not head or not head.next:
+def is_palindrome_stack(head: ListNode | None) -> bool:
+    """
+    Check if a linked list is a palindrome using a stack.
+
+    Args:
+        head (ListNode): The head of the linked list.
+
+    Returns:
+        bool: True if the linked list is a palindrome, False otherwise.
+
+    Examples:
+        >>> is_palindrome_stack(None)
+        True
+
+        >>> is_palindrome_stack(ListNode(1))
+        True
+
+        >>> is_palindrome_stack(ListNode(1, ListNode(2)))
+        False
+
+        >>> is_palindrome_stack(ListNode(1, ListNode(2, ListNode(1))))
+        True
+
+        >>> is_palindrome_stack(ListNode(1, ListNode(2, ListNode(2, ListNode(1)))))
+        True
+    """
+    if not head or not head.next_node:
         return True
 
     # 1. Get the midpoint (slow)
-    slow = fast = cur = head
-    while fast and fast.next:
-        fast, slow = fast.next.next, slow.next
+    slow: ListNode | None = head
+    fast: ListNode | None = head
+    while fast and fast.next_node:
+        fast = fast.next_node.next_node
+        slow = slow.next_node if slow else None
 
-    # 2. Push the second half into the stack
-    stack = [slow.val]
-    while slow.next:
-        slow = slow.next
-        stack.append(slow.val)
+    # slow will always be defined,
+    # adding this check to resolve mypy static check
+    if slow:
+        stack = [slow.val]
 
-    # 3. Comparison
-    while stack:
-        if stack.pop() != cur.val:
-            return False
-        cur = cur.next
+        # 2. Push the second half into the stack
+        while slow.next_node:
+            slow = slow.next_node
+            stack.append(slow.val)
+
+        # 3. Comparison
+        cur: ListNode | None = head
+        while stack and cur:
+            if stack.pop() != cur.val:
+                return False
+            cur = cur.next_node
 
     return True
 
 
-def is_palindrome_dict(head):
-    if not head or not head.next:
+def is_palindrome_dict(head: ListNode | None) -> bool:
+    """
+    Check if a linked list is a palindrome using a dictionary.
+
+    Args:
+        head (ListNode): The head of the linked list.
+
+    Returns:
+        bool: True if the linked list is a palindrome, False otherwise.
+
+    Examples:
+        >>> is_palindrome_dict(None)
+        True
+
+        >>> is_palindrome_dict(ListNode(1))
+        True
+
+        >>> is_palindrome_dict(ListNode(1, ListNode(2)))
+        False
+
+        >>> is_palindrome_dict(ListNode(1, ListNode(2, ListNode(1))))
+        True
+
+        >>> is_palindrome_dict(ListNode(1, ListNode(2, ListNode(2, ListNode(1)))))
+        True
+
+        >>> is_palindrome_dict(
+        ...     ListNode(
+        ...         1, ListNode(2, ListNode(1, ListNode(3, ListNode(2, ListNode(1)))))
+        ...     )
+        ... )
+        False
+    """
+    if not head or not head.next_node:
         return True
-    d = {}
+    d: dict[int, list[int]] = {}
     pos = 0
     while head:
         if head.val in d:
             d[head.val].append(pos)
         else:
             d[head.val] = [pos]
-        head = head.next
+        head = head.next_node
         pos += 1
     checksum = pos - 1
     middle = 0
@@ -75,3 +179,9 @@ def is_palindrome_dict(head):
         if middle > 1:
             return False
     return True
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/digital_image_processing/test_digital_image_processing.py b/digital_image_processing/test_digital_image_processing.py
index 528b4bc3b..2e5630458 100644
--- a/digital_image_processing/test_digital_image_processing.py
+++ b/digital_image_processing/test_digital_image_processing.py
@@ -96,9 +96,16 @@ def test_nearest_neighbour(
 
 
 def test_local_binary_pattern():
-    file_path = "digital_image_processing/image_data/lena.jpg"
+    # pull request 10161 before:
+    # "digital_image_processing/image_data/lena.jpg"
+    # after: "digital_image_processing/image_data/lena_small.jpg"
 
-    # Reading the image and converting it to grayscale.
+    from os import getenv  # Speed up our Continuous Integration tests
+
+    file_name = "lena_small.jpg" if getenv("CI") else "lena.jpg"
+    file_path = f"digital_image_processing/image_data/{file_name}"
+
+    # Reading the image and converting it to grayscale
     image = imread(file_path, 0)
 
     # Test for get_neighbors_pixel function() return not None
diff --git a/dynamic_programming/minimum_steps_to_one.py b/dynamic_programming/minimum_steps_to_one.py
index 8785027fb..68eaf56e2 100644
--- a/dynamic_programming/minimum_steps_to_one.py
+++ b/dynamic_programming/minimum_steps_to_one.py
@@ -1,7 +1,7 @@
 """
 YouTube Explanation: https://www.youtube.com/watch?v=f2xi3c1S95M
 
-Given an integer n, return the minimum steps to 1
+Given an integer n, return the minimum steps from n to 1
 
 AVAILABLE STEPS:
     * Decrement by 1
diff --git a/greedy_methods/best_time_to_buy_and_sell_stock.py b/greedy_methods/best_time_to_buy_and_sell_stock.py
new file mode 100644
index 000000000..4aea19172
--- /dev/null
+++ b/greedy_methods/best_time_to_buy_and_sell_stock.py
@@ -0,0 +1,42 @@
+"""
+Given a list of stock prices calculate the maximum profit that can be made from a
+single buy and sell of one share of stock.  We only allowed to complete one buy
+transaction and one sell transaction but must buy before we sell.
+
+Example : prices = [7, 1, 5, 3, 6, 4]
+max_profit will return 5 - which is by buying at price 1 and selling at price 6.
+
+This problem can be solved using the concept of "GREEDY ALGORITHM".
+
+We iterate over the price array once, keeping track of the lowest price point
+(buy) and the maximum profit we can get at each point.  The greedy choice at each point
+is to either buy at the current price if it's less than our current buying price, or
+sell at the current price if the profit is more than our current maximum profit.
+"""
+
+
+def max_profit(prices: list[int]) -> int:
+    """
+    >>> max_profit([7, 1, 5, 3, 6, 4])
+    5
+    >>> max_profit([7, 6, 4, 3, 1])
+    0
+    """
+    if not prices:
+        return 0
+
+    min_price = prices[0]
+    max_profit: int = 0
+
+    for price in prices:
+        min_price = min(price, min_price)
+        max_profit = max(price - min_price, max_profit)
+
+    return max_profit
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
+    print(max_profit([7, 1, 5, 3, 6, 4]))
diff --git a/machine_learning/loss_functions/binary_cross_entropy.py b/machine_learning/loss_functions/binary_cross_entropy.py
new file mode 100644
index 000000000..4ebca7f21
--- /dev/null
+++ b/machine_learning/loss_functions/binary_cross_entropy.py
@@ -0,0 +1,59 @@
+"""
+Binary Cross-Entropy (BCE) Loss Function
+
+Description:
+Quantifies dissimilarity between true labels (0 or 1) and predicted probabilities.
+It's widely used in binary classification tasks.
+
+Formula:
+BCE = -Σ(y_true * log(y_pred) + (1 - y_true) * log(1 - y_pred))
+
+Source:
+[Wikipedia - Cross entropy](https://en.wikipedia.org/wiki/Cross_entropy)
+"""
+
+import numpy as np
+
+
+def binary_cross_entropy(
+    y_true: np.ndarray, y_pred: np.ndarray, epsilon: float = 1e-15
+) -> float:
+    """
+    Calculate the BCE Loss between true labels and predicted probabilities.
+
+    Parameters:
+    - y_true: True binary labels (0 or 1).
+    - y_pred: Predicted probabilities for class 1.
+    - epsilon: Small constant to avoid numerical instability.
+
+    Returns:
+    - bce_loss: Binary Cross-Entropy Loss.
+
+    Example Usage:
+    >>> true_labels = np.array([0, 1, 1, 0, 1])
+    >>> predicted_probs = np.array([0.2, 0.7, 0.9, 0.3, 0.8])
+    >>> binary_cross_entropy(true_labels, predicted_probs)
+    0.2529995012327421
+    >>> true_labels = np.array([0, 1, 1, 0, 1])
+    >>> predicted_probs = np.array([0.3, 0.8, 0.9, 0.2])
+    >>> binary_cross_entropy(true_labels, predicted_probs)
+    Traceback (most recent call last):
+        ...
+    ValueError: Input arrays must have the same length.
+    """
+    if len(y_true) != len(y_pred):
+        raise ValueError("Input arrays must have the same length.")
+    # Clip predicted probabilities to avoid log(0) and log(1)
+    y_pred = np.clip(y_pred, epsilon, 1 - epsilon)
+
+    # Calculate binary cross-entropy loss
+    bce_loss = -(y_true * np.log(y_pred) + (1 - y_true) * np.log(1 - y_pred))
+
+    # Take the mean over all samples
+    return np.mean(bce_loss)
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/machine_learning/loss_functions/categorical_cross_entropy.py b/machine_learning/loss_functions/categorical_cross_entropy.py
new file mode 100644
index 000000000..68f98902b
--- /dev/null
+++ b/machine_learning/loss_functions/categorical_cross_entropy.py
@@ -0,0 +1,85 @@
+"""
+Categorical Cross-Entropy Loss
+
+This function calculates the Categorical Cross-Entropy Loss between true class
+labels and predicted class probabilities.
+
+Formula:
+Categorical Cross-Entropy Loss = -Σ(y_true * ln(y_pred))
+
+Resources:
+- [Wikipedia - Cross entropy](https://en.wikipedia.org/wiki/Cross_entropy)
+"""
+
+import numpy as np
+
+
+def categorical_cross_entropy(
+    y_true: np.ndarray, y_pred: np.ndarray, epsilon: float = 1e-15
+) -> float:
+    """
+    Calculate Categorical Cross-Entropy Loss between true class labels and
+    predicted class probabilities.
+
+    Parameters:
+    - y_true: True class labels (one-hot encoded) as a NumPy array.
+    - y_pred: Predicted class probabilities as a NumPy array.
+    - epsilon: Small constant to avoid numerical instability.
+
+    Returns:
+    - ce_loss: Categorical Cross-Entropy Loss as a floating-point number.
+
+    Example:
+    >>> true_labels = np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]])
+    >>> pred_probs = np.array([[0.9, 0.1, 0.0], [0.2, 0.7, 0.1], [0.0, 0.1, 0.9]])
+    >>> categorical_cross_entropy(true_labels, pred_probs)
+    0.567395975254385
+
+    >>> y_true = np.array([[1, 0], [0, 1]])
+    >>> y_pred = np.array([[0.9, 0.1, 0.0], [0.2, 0.7, 0.1]])
+    >>> categorical_cross_entropy(y_true, y_pred)
+    Traceback (most recent call last):
+        ...
+    ValueError: Input arrays must have the same shape.
+
+    >>> y_true = np.array([[2, 0, 1], [1, 0, 0]])
+    >>> y_pred = np.array([[0.9, 0.1, 0.0], [0.2, 0.7, 0.1]])
+    >>> categorical_cross_entropy(y_true, y_pred)
+    Traceback (most recent call last):
+        ...
+    ValueError: y_true must be one-hot encoded.
+
+    >>> y_true = np.array([[1, 0, 1], [1, 0, 0]])
+    >>> y_pred = np.array([[0.9, 0.1, 0.0], [0.2, 0.7, 0.1]])
+    >>> categorical_cross_entropy(y_true, y_pred)
+    Traceback (most recent call last):
+        ...
+    ValueError: y_true must be one-hot encoded.
+
+    >>> y_true = np.array([[1, 0, 0], [0, 1, 0]])
+    >>> y_pred = np.array([[0.9, 0.1, 0.1], [0.2, 0.7, 0.1]])
+    >>> categorical_cross_entropy(y_true, y_pred)
+    Traceback (most recent call last):
+        ...
+    ValueError: Predicted probabilities must sum to approximately 1.
+    """
+    if y_true.shape != y_pred.shape:
+        raise ValueError("Input arrays must have the same shape.")
+
+    if np.any((y_true != 0) & (y_true != 1)) or np.any(y_true.sum(axis=1) != 1):
+        raise ValueError("y_true must be one-hot encoded.")
+
+    if not np.all(np.isclose(np.sum(y_pred, axis=1), 1, rtol=epsilon, atol=epsilon)):
+        raise ValueError("Predicted probabilities must sum to approximately 1.")
+
+    # Clip predicted probabilities to avoid log(0)
+    y_pred = np.clip(y_pred, epsilon, 1)
+
+    # Calculate categorical cross-entropy loss
+    return -np.sum(y_true * np.log(y_pred))
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/machine_learning/loss_functions/huber_loss.py b/machine_learning/loss_functions/huber_loss.py
new file mode 100644
index 000000000..202e013f2
--- /dev/null
+++ b/machine_learning/loss_functions/huber_loss.py
@@ -0,0 +1,52 @@
+"""
+Huber Loss Function
+
+Description:
+Huber loss function describes the penalty incurred by an estimation procedure.
+It serves as a measure of the model's accuracy in regression tasks.
+
+Formula:
+Huber Loss = if |y_true - y_pred| <= delta then 0.5 * (y_true - y_pred)^2
+             else delta * |y_true - y_pred| - 0.5 * delta^2
+
+Source:
+[Wikipedia - Huber Loss](https://en.wikipedia.org/wiki/Huber_loss)
+"""
+
+import numpy as np
+
+
+def huber_loss(y_true: np.ndarray, y_pred: np.ndarray, delta: float) -> float:
+    """
+    Calculate the mean of Huber Loss.
+
+    Parameters:
+    - y_true: The true values (ground truth).
+    - y_pred: The predicted values.
+
+    Returns:
+    - huber_loss: The mean of Huber Loss between y_true and y_pred.
+
+    Example usage:
+    >>> true_values = np.array([0.9, 10.0, 2.0, 1.0, 5.2])
+    >>> predicted_values = np.array([0.8, 2.1, 2.9, 4.2, 5.2])
+    >>> np.isclose(huber_loss(true_values, predicted_values, 1.0), 2.102)
+    True
+    >>> true_labels = np.array([11.0, 21.0, 3.32, 4.0, 5.0])
+    >>> predicted_probs = np.array([8.3, 20.8, 2.9, 11.2, 5.0])
+    >>> np.isclose(huber_loss(true_labels, predicted_probs, 1.0), 1.80164)
+    True
+    """
+
+    if len(y_true) != len(y_pred):
+        raise ValueError("Input arrays must have the same length.")
+
+    huber_mse = 0.5 * (y_true - y_pred) ** 2
+    huber_mae = delta * (np.abs(y_true - y_pred) - 0.5 * delta)
+    return np.where(np.abs(y_true - y_pred) <= delta, huber_mse, huber_mae).mean()
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/machine_learning/loss_functions/mean_squared_error.py b/machine_learning/loss_functions/mean_squared_error.py
new file mode 100644
index 000000000..d2b0e1e15
--- /dev/null
+++ b/machine_learning/loss_functions/mean_squared_error.py
@@ -0,0 +1,51 @@
+"""
+Mean Squared Error (MSE) Loss Function
+
+Description:
+MSE measures the mean squared difference between true values and predicted values.
+It serves as a measure of the model's accuracy in regression tasks.
+
+Formula:
+MSE = (1/n) * Σ(y_true - y_pred)^2
+
+Source:
+[Wikipedia - Mean squared error](https://en.wikipedia.org/wiki/Mean_squared_error)
+"""
+
+import numpy as np
+
+
+def mean_squared_error(y_true: np.ndarray, y_pred: np.ndarray) -> float:
+    """
+    Calculate the Mean Squared Error (MSE) between two arrays.
+
+    Parameters:
+    - y_true: The true values (ground truth).
+    - y_pred: The predicted values.
+
+    Returns:
+    - mse: The Mean Squared Error between y_true and y_pred.
+
+    Example usage:
+    >>> true_values = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
+    >>> predicted_values = np.array([0.8, 2.1, 2.9, 4.2, 5.2])
+    >>> mean_squared_error(true_values, predicted_values)
+    0.028000000000000032
+    >>> true_labels = np.array([1.0, 2.0, 3.0, 4.0, 5.0])
+    >>> predicted_probs = np.array([0.3, 0.8, 0.9, 0.2])
+    >>> mean_squared_error(true_labels, predicted_probs)
+    Traceback (most recent call last):
+        ...
+    ValueError: Input arrays must have the same length.
+    """
+    if len(y_true) != len(y_pred):
+        raise ValueError("Input arrays must have the same length.")
+
+    squared_errors = (y_true - y_pred) ** 2
+    return np.mean(squared_errors)
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/maths/binary_exponentiation_2.py b/maths/binary_exponentiation_2.py
new file mode 100644
index 000000000..edb6b66b2
--- /dev/null
+++ b/maths/binary_exponentiation_2.py
@@ -0,0 +1,61 @@
+"""
+Binary Exponentiation
+This is a method to find a^b in O(log b) time complexity
+This is one of the most commonly used methods of exponentiation
+It's also useful when the solution to (a^b) % c is required because a, b, c may be
+over the computer's calculation limits
+
+Let's say you need to calculate a ^ b
+- RULE 1 : a ^ b = (a*a) ^ (b/2) ---- example : 4 ^ 4 = (4*4) ^ (4/2) = 16 ^ 2
+- RULE 2 : IF b is odd, then a ^ b = a * (a ^ (b - 1)), where b - 1 is even
+Once b is even, repeat the process until b = 1 or b = 0, because a^1 = a and a^0 = 1
+
+For modular exponentiation, we use the fact that (a*b) % c = ((a%c) * (b%c)) % c
+Now apply RULE 1 or 2 as required
+
+@author chinmoy159
+"""
+
+
+def b_expo(a: int, b: int) -> int:
+    """
+    >>> b_expo(2, 10)
+    1024
+    >>> b_expo(9, 0)
+    1
+    >>> b_expo(0, 12)
+    0
+    >>> b_expo(4, 12)
+    16777216
+    """
+    res = 1
+    while b > 0:
+        if b & 1:
+            res *= a
+
+        a *= a
+        b >>= 1
+
+    return res
+
+
+def b_expo_mod(a: int, b: int, c: int) -> int:
+    """
+    >>> b_expo_mod(2, 10, 1000000007)
+    1024
+    >>> b_expo_mod(11, 13, 19)
+    11
+    >>> b_expo_mod(0, 19, 20)
+    0
+    >>> b_expo_mod(15, 5, 4)
+    3
+    """
+    res = 1
+    while b > 0:
+        if b & 1:
+            res = ((res % c) * (a % c)) % c
+
+        a *= a
+        b >>= 1
+
+    return res
diff --git a/maths/binary_exponentiation_3.py b/maths/binary_exponentiation_3.py
deleted file mode 100644
index 9cd143e09..000000000
--- a/maths/binary_exponentiation_3.py
+++ /dev/null
@@ -1,50 +0,0 @@
-"""
-* Binary Exponentiation for Powers
-* This is a method to find a^b in a time complexity of O(log b)
-* This is one of the most commonly used methods of finding powers.
-* Also useful in cases where solution to (a^b)%c is required,
-* where a,b,c can be numbers over the computers calculation limits.
-* Done using iteration, can also be done using recursion
-
-* @author chinmoy159
-* @version 1.0 dated 10/08/2017
-"""
-
-
-def b_expo(a: int, b: int) -> int:
-    res = 1
-    while b > 0:
-        if b & 1:
-            res *= a
-
-        a *= a
-        b >>= 1
-
-    return res
-
-
-def b_expo_mod(a: int, b: int, c: int) -> int:
-    res = 1
-    while b > 0:
-        if b & 1:
-            res = ((res % c) * (a % c)) % c
-
-        a *= a
-        b >>= 1
-
-    return res
-
-
-"""
-* Wondering how this method works !
-* It's pretty simple.
-* Let's say you need to calculate a ^ b
-* RULE 1 : a ^ b = (a*a) ^ (b/2) ---- example : 4 ^ 4 = (4*4) ^ (4/2) = 16 ^ 2
-* RULE 2 : IF b is ODD, then ---- a ^ b = a * (a ^ (b - 1)) :: where (b - 1) is even.
-* Once b is even, repeat the process to get a ^ b
-* Repeat the process till b = 1 OR b = 0, because a^1 = a AND a^0 = 1
-*
-* As far as the modulo is concerned,
-* the fact : (a*b) % c = ((a%c) * (b%c)) % c
-* Now apply RULE 1 OR 2 whichever is required.
-"""
diff --git a/maths/carmichael_number.py b/maths/carmichael_number.py
index c9c144759..81712520f 100644
--- a/maths/carmichael_number.py
+++ b/maths/carmichael_number.py
@@ -10,14 +10,7 @@ satisfies the following modular arithmetic condition:
 Examples of Carmichael Numbers: 561, 1105, ...
 https://en.wikipedia.org/wiki/Carmichael_number
 """
-
-
-def gcd(a: int, b: int) -> int:
-    if a < b:
-        return gcd(b, a)
-    if a % b == 0:
-        return b
-    return gcd(b, a % b)
+from maths.greatest_common_divisor import greatest_common_divisor
 
 
 def power(x: int, y: int, mod: int) -> int:
@@ -33,7 +26,7 @@ def power(x: int, y: int, mod: int) -> int:
 def is_carmichael_number(n: int) -> bool:
     b = 2
     while b < n:
-        if gcd(b, n) == 1 and power(b, n - 1, n) != 1:
+        if greatest_common_divisor(b, n) == 1 and power(b, n - 1, n) != 1:
             return False
         b += 1
     return True
diff --git a/maths/chebyshev_distance.py b/maths/chebyshev_distance.py
new file mode 100644
index 000000000..4801d3916
--- /dev/null
+++ b/maths/chebyshev_distance.py
@@ -0,0 +1,20 @@
+def chebyshev_distance(point_a: list[float], point_b: list[float]) -> float:
+    """
+    This function calculates the Chebyshev distance (also known as the
+    Chessboard distance) between two n-dimensional points represented as lists.
+
+    https://en.wikipedia.org/wiki/Chebyshev_distance
+
+    >>> chebyshev_distance([1.0, 1.0], [2.0, 2.0])
+    1.0
+    >>> chebyshev_distance([1.0, 1.0, 9.0], [2.0, 2.0, -5.2])
+    14.2
+    >>> chebyshev_distance([1.0], [2.0, 2.0])
+    Traceback (most recent call last):
+        ...
+    ValueError: Both points must have the same dimension.
+    """
+    if len(point_a) != len(point_b):
+        raise ValueError("Both points must have the same dimension.")
+
+    return max(abs(a - b) for a, b in zip(point_a, point_b))
diff --git a/blockchain/chinese_remainder_theorem.py b/maths/chinese_remainder_theorem.py
similarity index 100%
rename from blockchain/chinese_remainder_theorem.py
rename to maths/chinese_remainder_theorem.py
diff --git a/maths/combinations.py b/maths/combinations.py
index a2324012c..6e9e1a807 100644
--- a/maths/combinations.py
+++ b/maths/combinations.py
@@ -1,7 +1,6 @@
 """
 https://en.wikipedia.org/wiki/Combination
 """
-from math import factorial
 
 
 def combinations(n: int, k: int) -> int:
@@ -35,7 +34,11 @@ def combinations(n: int, k: int) -> int:
     # to calculate a factorial of a negative number, which is not possible
     if n < k or k < 0:
         raise ValueError("Please enter positive integers for n and k where n >= k")
-    return factorial(n) // (factorial(k) * factorial(n - k))
+    res = 1
+    for i in range(k):
+        res *= n - i
+        res //= i + 1
+    return res
 
 
 if __name__ == "__main__":
diff --git a/maths/factorial.py b/maths/factorial.py
index 18cacdef9..aaf90f384 100644
--- a/maths/factorial.py
+++ b/maths/factorial.py
@@ -1,4 +1,5 @@
-"""Factorial of a positive integer -- https://en.wikipedia.org/wiki/Factorial
+"""
+Factorial of a positive integer -- https://en.wikipedia.org/wiki/Factorial
 """
 
 
diff --git a/maths/germain_primes.py b/maths/germain_primes.py
new file mode 100644
index 000000000..078d1967f
--- /dev/null
+++ b/maths/germain_primes.py
@@ -0,0 +1,72 @@
+"""
+A Sophie Germain prime is any prime p, where 2p + 1 is also prime.
+The second number, 2p + 1 is called a safe prime.
+
+Examples of Germain primes include: 2, 3, 5, 11, 23
+
+Their corresponding safe primes: 5, 7, 11, 23, 47
+https://en.wikipedia.org/wiki/Safe_and_Sophie_Germain_primes
+"""
+
+from maths.prime_check import is_prime
+
+
+def is_germain_prime(number: int) -> bool:
+    """Checks if input number and 2*number + 1 are prime.
+
+    >>> is_germain_prime(3)
+    True
+    >>> is_germain_prime(11)
+    True
+    >>> is_germain_prime(4)
+    False
+    >>> is_germain_prime(23)
+    True
+    >>> is_germain_prime(13)
+    False
+    >>> is_germain_prime(20)
+    False
+    >>> is_germain_prime('abc')
+    Traceback (most recent call last):
+        ...
+    TypeError: Input value must be a positive integer. Input value: abc
+    """
+    if not isinstance(number, int) or number < 1:
+        msg = f"Input value must be a positive integer. Input value: {number}"
+        raise TypeError(msg)
+
+    return is_prime(number) and is_prime(2 * number + 1)
+
+
+def is_safe_prime(number: int) -> bool:
+    """Checks if input number and (number - 1)/2 are prime.
+    The smallest safe prime is 5, with the Germain prime is 2.
+
+    >>> is_safe_prime(5)
+    True
+    >>> is_safe_prime(11)
+    True
+    >>> is_safe_prime(1)
+    False
+    >>> is_safe_prime(2)
+    False
+    >>> is_safe_prime(3)
+    False
+    >>> is_safe_prime(47)
+    True
+    >>> is_safe_prime('abc')
+    Traceback (most recent call last):
+        ...
+    TypeError: Input value must be a positive integer. Input value: abc
+    """
+    if not isinstance(number, int) or number < 1:
+        msg = f"Input value must be a positive integer. Input value: {number}"
+        raise TypeError(msg)
+
+    return (number - 1) % 2 == 0 and is_prime(number) and is_prime((number - 1) // 2)
+
+
+if __name__ == "__main__":
+    from doctest import testmod
+
+    testmod()
diff --git a/maths/least_common_multiple.py b/maths/least_common_multiple.py
index 10cc63ac7..4f28da8ab 100644
--- a/maths/least_common_multiple.py
+++ b/maths/least_common_multiple.py
@@ -1,6 +1,8 @@
 import unittest
 from timeit import timeit
 
+from maths.greatest_common_divisor import greatest_common_divisor
+
 
 def least_common_multiple_slow(first_num: int, second_num: int) -> int:
     """
@@ -20,26 +22,6 @@ def least_common_multiple_slow(first_num: int, second_num: int) -> int:
     return common_mult
 
 
-def greatest_common_divisor(a: int, b: int) -> int:
-    """
-    Calculate Greatest Common Divisor (GCD).
-    see greatest_common_divisor.py
-    >>> greatest_common_divisor(24, 40)
-    8
-    >>> greatest_common_divisor(1, 1)
-    1
-    >>> greatest_common_divisor(1, 800)
-    1
-    >>> greatest_common_divisor(11, 37)
-    1
-    >>> greatest_common_divisor(3, 5)
-    1
-    >>> greatest_common_divisor(16, 4)
-    4
-    """
-    return b if a == 0 else greatest_common_divisor(b % a, a)
-
-
 def least_common_multiple_fast(first_num: int, second_num: int) -> int:
     """
     Find the least common multiple of two numbers.
diff --git a/maths/minkowski_distance.py b/maths/minkowski_distance.py
new file mode 100644
index 000000000..3237124e8
--- /dev/null
+++ b/maths/minkowski_distance.py
@@ -0,0 +1,45 @@
+def minkowski_distance(
+    point_a: list[float],
+    point_b: list[float],
+    order: int,
+) -> float:
+    """
+    This function calculates the Minkowski distance for a given order between
+    two n-dimensional points represented as lists. For the case of order = 1,
+    the Minkowski distance degenerates to the Manhattan distance. For
+    order = 2, the usual Euclidean distance is obtained.
+
+    https://en.wikipedia.org/wiki/Minkowski_distance
+
+    Note: due to floating point calculation errors the output of this
+    function may be inaccurate.
+
+    >>> minkowski_distance([1.0, 1.0], [2.0, 2.0], 1)
+    2.0
+    >>> minkowski_distance([1.0, 2.0, 3.0, 4.0], [5.0, 6.0, 7.0, 8.0], 2)
+    8.0
+    >>> import numpy as np
+    >>> np.isclose(5.0, minkowski_distance([5.0], [0.0], 3))
+    True
+    >>> minkowski_distance([1.0], [2.0], -1)
+    Traceback (most recent call last):
+        ...
+    ValueError: The order must be greater than or equal to 1.
+    >>> minkowski_distance([1.0], [1.0, 2.0], 1)
+    Traceback (most recent call last):
+        ...
+    ValueError: Both points must have the same dimension.
+    """
+    if order < 1:
+        raise ValueError("The order must be greater than or equal to 1.")
+
+    if len(point_a) != len(point_b):
+        raise ValueError("Both points must have the same dimension.")
+
+    return sum(abs(a - b) ** order for a, b in zip(point_a, point_b)) ** (1 / order)
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/blockchain/modular_division.py b/maths/modular_division.py
similarity index 100%
rename from blockchain/modular_division.py
rename to maths/modular_division.py
diff --git a/maths/number_of_digits.py b/maths/number_of_digits.py
index 86bc67f72..bb9c0d248 100644
--- a/maths/number_of_digits.py
+++ b/maths/number_of_digits.py
@@ -16,7 +16,15 @@ def num_digits(n: int) -> int:
     1
     >>> num_digits(-123456)
     6
+    >>> num_digits('123')  # Raises a TypeError for non-integer input
+    Traceback (most recent call last):
+        ...
+    TypeError: Input must be an integer
     """
+
+    if not isinstance(n, int):
+        raise TypeError("Input must be an integer")
+
     digits = 0
     n = abs(n)
     while True:
@@ -42,7 +50,15 @@ def num_digits_fast(n: int) -> int:
     1
     >>> num_digits_fast(-123456)
     6
+    >>> num_digits('123')  # Raises a TypeError for non-integer input
+    Traceback (most recent call last):
+        ...
+    TypeError: Input must be an integer
     """
+
+    if not isinstance(n, int):
+        raise TypeError("Input must be an integer")
+
     return 1 if n == 0 else math.floor(math.log(abs(n), 10) + 1)
 
 
@@ -61,7 +77,15 @@ def num_digits_faster(n: int) -> int:
     1
     >>> num_digits_faster(-123456)
     6
+    >>> num_digits('123')  # Raises a TypeError for non-integer input
+    Traceback (most recent call last):
+        ...
+    TypeError: Input must be an integer
     """
+
+    if not isinstance(n, int):
+        raise TypeError("Input must be an integer")
+
     return len(str(abs(n)))
 
 
diff --git a/maths/primelib.py b/maths/primelib.py
index 28b5aee9d..cf01750cf 100644
--- a/maths/primelib.py
+++ b/maths/primelib.py
@@ -21,7 +21,6 @@ get_primes_between(pNumber1, pNumber2)
 
 is_even(number)
 is_odd(number)
-gcd(number1, number2)  // greatest common divisor
 kg_v(number1, number2)  // least common multiple
 get_divisors(number)    // all divisors of 'number' inclusive 1, number
 is_perfect_number(number)
@@ -40,6 +39,8 @@ goldbach(number)  // Goldbach's assumption
 
 from math import sqrt
 
+from maths.greatest_common_divisor import gcd_by_iterative
+
 
 def is_prime(number: int) -> bool:
     """
@@ -317,39 +318,6 @@ def goldbach(number):
 # ----------------------------------------------
 
 
-def gcd(number1, number2):
-    """
-    Greatest common divisor
-    input: two positive integer 'number1' and 'number2'
-    returns the greatest common divisor of 'number1' and 'number2'
-    """
-
-    # precondition
-    assert (
-        isinstance(number1, int)
-        and isinstance(number2, int)
-        and (number1 >= 0)
-        and (number2 >= 0)
-    ), "'number1' and 'number2' must been positive integer."
-
-    rest = 0
-
-    while number2 != 0:
-        rest = number1 % number2
-        number1 = number2
-        number2 = rest
-
-    # precondition
-    assert isinstance(number1, int) and (
-        number1 >= 0
-    ), "'number' must been from type int and positive"
-
-    return number1
-
-
-# ----------------------------------------------------
-
-
 def kg_v(number1, number2):
     """
     Least common multiple
@@ -567,14 +535,14 @@ def simplify_fraction(numerator, denominator):
     ), "The arguments must been from type int and 'denominator' != 0"
 
     # build the greatest common divisor of numerator and denominator.
-    gcd_of_fraction = gcd(abs(numerator), abs(denominator))
+    gcd_of_fraction = gcd_by_iterative(abs(numerator), abs(denominator))
 
     # precondition
     assert (
         isinstance(gcd_of_fraction, int)
         and (numerator % gcd_of_fraction == 0)
         and (denominator % gcd_of_fraction == 0)
-    ), "Error in function gcd(...,...)"
+    ), "Error in function gcd_by_iterative(...,...)"
 
     return (numerator // gcd_of_fraction, denominator // gcd_of_fraction)
 
diff --git a/neural_network/activation_functions/binary_step.py b/neural_network/activation_functions/binary_step.py
new file mode 100644
index 000000000..8f8f4d405
--- /dev/null
+++ b/neural_network/activation_functions/binary_step.py
@@ -0,0 +1,36 @@
+"""
+This script demonstrates the implementation of the Binary Step function.
+
+It's an activation function in which the neuron is activated if the input is positive
+or 0, else it is deactivated
+
+It's a simple activation function which is mentioned in this wikipedia article:
+https://en.wikipedia.org/wiki/Activation_function
+"""
+
+
+import numpy as np
+
+
+def binary_step(vector: np.ndarray) -> np.ndarray:
+    """
+    Implements the binary step function
+
+    Parameters:
+        vector (ndarray): A vector that consists of numeric values
+
+    Returns:
+        vector (ndarray): Input vector after applying binary step function
+
+    >>> vector = np.array([-1.2, 0, 2, 1.45, -3.7, 0.3])
+    >>> binary_step(vector)
+    array([0, 1, 1, 1, 0, 1])
+    """
+
+    return np.where(vector >= 0, 1, 0)
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/neural_network/activation_functions/mish.py b/neural_network/activation_functions/mish.py
index e4f98307f..e51655df8 100644
--- a/neural_network/activation_functions/mish.py
+++ b/neural_network/activation_functions/mish.py
@@ -7,6 +7,7 @@ https://en.wikipedia.org/wiki/Rectifier_(neural_networks)#Mish
 """
 
 import numpy as np
+from softplus import softplus
 
 
 def mish(vector: np.ndarray) -> np.ndarray:
@@ -30,7 +31,7 @@ def mish(vector: np.ndarray) -> np.ndarray:
     array([-0.00092952, -0.15113318,  0.33152014, -0.04745745])
 
     """
-    return vector * np.tanh(np.log(1 + np.exp(vector)))
+    return vector * np.tanh(softplus(vector))
 
 
 if __name__ == "__main__":
diff --git a/neural_network/activation_functions/soboleva_modified_hyperbolic_tangent.py b/neural_network/activation_functions/soboleva_modified_hyperbolic_tangent.py
new file mode 100644
index 000000000..603ac0b7e
--- /dev/null
+++ b/neural_network/activation_functions/soboleva_modified_hyperbolic_tangent.py
@@ -0,0 +1,49 @@
+"""
+This script implements the Soboleva Modified Hyperbolic Tangent function.
+
+The function applies the Soboleva Modified Hyperbolic Tangent function
+to each element of the vector.
+
+More details about the activation function can be found on:
+https://en.wikipedia.org/wiki/Soboleva_modified_hyperbolic_tangent
+"""
+
+
+import numpy as np
+
+
+def soboleva_modified_hyperbolic_tangent(
+    vector: np.ndarray, a_value: float, b_value: float, c_value: float, d_value: float
+) -> np.ndarray:
+    """
+    Implements the Soboleva Modified Hyperbolic Tangent function
+
+    Parameters:
+        vector (ndarray): A vector that consists of numeric values
+        a_value (float): parameter a of the equation
+        b_value (float): parameter b of the equation
+        c_value (float): parameter c of the equation
+        d_value (float): parameter d of the equation
+
+    Returns:
+        vector (ndarray): Input array after applying SMHT function
+
+    >>> vector = np.array([5.4, -2.4, 6.3, -5.23, 3.27, 0.56])
+    >>> soboleva_modified_hyperbolic_tangent(vector, 0.2, 0.4, 0.6, 0.8)
+    array([ 0.11075085, -0.28236685,  0.07861169, -0.1180085 ,  0.22999056,
+            0.1566043 ])
+    """
+
+    # Separate the numerator and denominator for simplicity
+    # Calculate the numerator and denominator element-wise
+    numerator = np.exp(a_value * vector) - np.exp(-b_value * vector)
+    denominator = np.exp(c_value * vector) + np.exp(-d_value * vector)
+
+    # Calculate and return the final result element-wise
+    return numerator / denominator
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/neural_network/activation_functions/squareplus.py b/neural_network/activation_functions/squareplus.py
new file mode 100644
index 000000000..40fa800d6
--- /dev/null
+++ b/neural_network/activation_functions/squareplus.py
@@ -0,0 +1,38 @@
+"""
+Squareplus Activation Function
+
+Use Case: Squareplus designed to enhance positive values and suppress negative values.
+For more detailed information, you can refer to the following link:
+https://en.wikipedia.org/wiki/Rectifier_(neural_networks)#Squareplus
+"""
+
+import numpy as np
+
+
+def squareplus(vector: np.ndarray, beta: float) -> np.ndarray:
+    """
+    Implements the SquarePlus activation function.
+
+    Parameters:
+        vector (np.ndarray): The input array for the SquarePlus activation.
+        beta (float): size of the curved region
+
+    Returns:
+        np.ndarray: The input array after applying the SquarePlus activation.
+
+    Formula: f(x) = ( x + sqrt(x^2 + b) ) / 2
+
+    Examples:
+    >>> squareplus(np.array([2.3, 0.6, -2, -3.8]), beta=2)
+    array([2.5       , 1.06811457, 0.22474487, 0.12731349])
+
+    >>> squareplus(np.array([-9.2, -0.3, 0.45, -4.56]), beta=3)
+    array([0.0808119 , 0.72891979, 1.11977651, 0.15893419])
+    """
+    return (vector + np.sqrt(vector**2 + beta)) / 2
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/other/magicdiamondpattern.py b/other/magicdiamondpattern.py
index 89b973bb4..58889280a 100644
--- a/other/magicdiamondpattern.py
+++ b/other/magicdiamondpattern.py
@@ -4,52 +4,76 @@
 # Function to print upper half of diamond (pyramid)
 def floyd(n):
     """
-        Parameters:
-    n : size of pattern
+    Print the upper half of a diamond pattern with '*' characters.
+
+    Args:
+        n (int): Size of the pattern.
+
+    Examples:
+        >>> floyd(3)
+        '  * \\n * * \\n* * * \\n'
+
+        >>> floyd(5)
+        '    * \\n   * * \\n  * * * \\n * * * * \\n* * * * * \\n'
     """
+    result = ""
     for i in range(n):
         for _ in range(n - i - 1):  # printing spaces
-            print(" ", end="")
+            result += " "
         for _ in range(i + 1):  # printing stars
-            print("* ", end="")
-        print()
+            result += "* "
+        result += "\n"
+    return result
 
 
 # Function to print lower half of diamond (pyramid)
 def reverse_floyd(n):
     """
-        Parameters:
-    n : size of pattern
+    Print the lower half of a diamond pattern with '*' characters.
+
+    Args:
+        n (int): Size of the pattern.
+
+    Examples:
+        >>> reverse_floyd(3)
+        '* * * \\n * * \\n  * \\n   '
+
+        >>> reverse_floyd(5)
+        '* * * * * \\n * * * * \\n  * * * \\n   * * \\n    * \\n     '
     """
+    result = ""
     for i in range(n, 0, -1):
         for _ in range(i, 0, -1):  # printing stars
-            print("* ", end="")
-        print()
+            result += "* "
+        result += "\n"
         for _ in range(n - i + 1, 0, -1):  # printing spaces
-            print(" ", end="")
+            result += " "
+    return result
 
 
 # Function to print complete diamond pattern of "*"
 def pretty_print(n):
     """
-        Parameters:
-    n : size of pattern
+    Print a complete diamond pattern with '*' characters.
+
+    Args:
+        n (int): Size of the pattern.
+
+    Examples:
+        >>> pretty_print(0)
+        '       ...       ....        nothing printing :('
+
+        >>> pretty_print(3)
+        '  * \\n * * \\n* * * \\n* * * \\n * * \\n  * \\n   '
     """
     if n <= 0:
-        print("       ...       ....        nothing printing :(")
-        return
-    floyd(n)  # upper half
-    reverse_floyd(n)  # lower half
+        return "       ...       ....        nothing printing :("
+    upper_half = floyd(n)  # upper half
+    lower_half = reverse_floyd(n)  # lower half
+    return upper_half + lower_half
 
 
 if __name__ == "__main__":
-    print(r"| /\ | |- |  |-  |--| |\  /| |-")
-    print(r"|/  \| |- |_ |_  |__| | \/ | |_")
-    K = 1
-    while K:
-        user_number = int(input("enter the number and , and see the magic : "))
-        print()
-        pretty_print(user_number)
-        K = int(input("press 0 to exit... and 1 to continue..."))
+    import doctest
 
-    print("Good Bye...")
+    doctest.testmod()
diff --git a/physics/reynolds_number.py b/physics/reynolds_number.py
new file mode 100644
index 000000000..dffe690f8
--- /dev/null
+++ b/physics/reynolds_number.py
@@ -0,0 +1,63 @@
+"""
+Title : computing the Reynolds number to find
+        out the type of flow (laminar or turbulent)
+
+Reynolds number is a dimensionless quantity that is used to determine
+the type of flow pattern as laminar or turbulent while flowing through a
+pipe. Reynolds number is defined by the ratio of inertial forces to that of
+viscous forces.
+
+R = Inertial Forces / Viscous Forces
+R = (ρ * V * D)/μ
+
+where :
+ρ = Density of fluid (in Kg/m^3)
+D = Diameter of pipe through which fluid flows (in m)
+V = Velocity of flow of the fluid (in m/s)
+μ = Viscosity of the fluid (in Ns/m^2)
+
+If the Reynolds number calculated is high (greater than 2000), then the
+flow through the pipe is said to be turbulent. If Reynolds number is low
+(less than 2000), the flow is said to be laminar. Numerically, these are
+acceptable values, although in general the laminar and turbulent flows
+are classified according to a range. Laminar flow falls below Reynolds
+number of 1100 and turbulent falls in a range greater than 2200.
+Laminar flow is the type of flow in which the fluid travels smoothly in
+regular paths. Conversely, turbulent flow isn't smooth and follows an
+irregular path with lots of mixing.
+
+Reference : https://byjus.com/physics/reynolds-number/
+"""
+
+
+def reynolds_number(
+    density: float, velocity: float, diameter: float, viscosity: float
+) -> float:
+    """
+    >>> reynolds_number(900, 2.5, 0.05, 0.4)
+    281.25
+    >>> reynolds_number(450, 3.86, 0.078, 0.23)
+    589.0695652173912
+    >>> reynolds_number(234, -4.5, 0.3, 0.44)
+    717.9545454545454
+    >>> reynolds_number(-90, 2, 0.045, 1)
+    Traceback (most recent call last):
+        ...
+    ValueError: please ensure that density, diameter and viscosity are positive
+    >>> reynolds_number(0, 2, -0.4, -2)
+    Traceback (most recent call last):
+        ...
+    ValueError: please ensure that density, diameter and viscosity are positive
+    """
+
+    if density <= 0 or diameter <= 0 or viscosity <= 0:
+        raise ValueError(
+            "please ensure that density, diameter and viscosity are positive"
+        )
+    return (density * abs(velocity) * diameter) / viscosity
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()
diff --git a/project_euler/problem_005/sol2.py b/project_euler/problem_005/sol2.py
index 1b3e5e130..4558e21fd 100644
--- a/project_euler/problem_005/sol2.py
+++ b/project_euler/problem_005/sol2.py
@@ -1,3 +1,5 @@
+from maths.greatest_common_divisor import greatest_common_divisor
+
 """
 Project Euler Problem 5: https://projecteuler.net/problem=5
 
@@ -16,23 +18,6 @@ References:
 """
 
 
-def greatest_common_divisor(x: int, y: int) -> int:
-    """
-    Euclidean Greatest Common Divisor algorithm
-
-    >>> greatest_common_divisor(0, 0)
-    0
-    >>> greatest_common_divisor(23, 42)
-    1
-    >>> greatest_common_divisor(15, 33)
-    3
-    >>> greatest_common_divisor(12345, 67890)
-    15
-    """
-
-    return x if y == 0 else greatest_common_divisor(y, x % y)
-
-
 def lcm(x: int, y: int) -> int:
     """
     Least Common Multiple.
diff --git a/searches/median_of_medians.py b/searches/median_of_medians.py
new file mode 100644
index 000000000..a8011a34a
--- /dev/null
+++ b/searches/median_of_medians.py
@@ -0,0 +1,107 @@
+"""
+A Python implementation of the Median of Medians algorithm
+to select pivots for quick_select, which is efficient for
+calculating the value that would appear in the index of a
+list if it would be sorted, even if it is not already
+sorted. Search in time complexity O(n) at any rank
+deterministically
+https://en.wikipedia.org/wiki/Median_of_medians
+"""
+
+
+def median_of_five(arr: list) -> int:
+    """
+    Return the median of the input list
+    :param arr: Array to find median of
+    :return: median of arr
+
+    >>> median_of_five([2, 4, 5, 7, 899])
+    5
+    >>> median_of_five([5, 7, 899, 54, 32])
+    32
+    >>> median_of_five([5, 4, 3, 2])
+    4
+    >>> median_of_five([3, 5, 7, 10, 2])
+    5
+    """
+    arr = sorted(arr)
+    return arr[len(arr) // 2]
+
+
+def median_of_medians(arr: list) -> int:
+    """
+    Return a pivot to partition data on by calculating
+    Median of medians of input data
+    :param arr: The data to be checked (a list)
+    :return: median of medians of input array
+
+    >>> median_of_medians([2, 4, 5, 7, 899, 54, 32])
+    54
+    >>> median_of_medians([5, 7, 899, 54, 32])
+    32
+    >>> median_of_medians([5, 4, 3, 2])
+    4
+    >>> median_of_medians([3, 5, 7, 10, 2, 12])
+    12
+    """
+
+    if len(arr) <= 5:
+        return median_of_five(arr)
+    medians = []
+    i = 0
+    while i < len(arr):
+        if (i + 4) <= len(arr):
+            medians.append(median_of_five(arr[i:].copy()))
+        else:
+            medians.append(median_of_five(arr[i : i + 5].copy()))
+        i += 5
+    return median_of_medians(medians)
+
+
+def quick_select(arr: list, target: int) -> int:
+    """
+    Two way partition the data into smaller and greater lists,
+    in relationship to the pivot
+    :param arr: The data to be searched (a list)
+    :param target: The rank to be searched
+    :return: element at rank target
+
+    >>> quick_select([2, 4, 5, 7, 899, 54, 32], 5)
+    32
+    >>> quick_select([2, 4, 5, 7, 899, 54, 32], 1)
+    2
+    >>> quick_select([5, 4, 3, 2], 2)
+    3
+    >>> quick_select([3, 5, 7, 10, 2, 12], 3)
+    5
+    """
+
+    # Invalid Input
+    if target > len(arr):
+        return -1
+
+    # x is the estimated pivot by median of medians algorithm
+    x = median_of_medians(arr)
+    left = []
+    right = []
+    check = False
+    for i in range(len(arr)):
+        if arr[i] < x:
+            left.append(arr[i])
+        elif arr[i] > x:
+            right.append(arr[i])
+        elif arr[i] == x and not check:
+            check = True
+        else:
+            right.append(arr[i])
+    rank_x = len(left) + 1
+    if rank_x == target:
+        answer = x
+    elif rank_x > target:
+        answer = quick_select(left, target)
+    elif rank_x < target:
+        answer = quick_select(right, target - rank_x)
+    return answer
+
+
+print(median_of_five([5, 4, 3, 2]))
diff --git a/strings/reverse_letters.py b/strings/reverse_letters.py
index 10b8a6d72..4f73f816b 100644
--- a/strings/reverse_letters.py
+++ b/strings/reverse_letters.py
@@ -1,19 +1,24 @@
-def reverse_letters(input_str: str) -> str:
+def reverse_letters(sentence: str, length: int = 0) -> str:
     """
-    Reverses letters in a given string without adjusting the position of the words
-    >>> reverse_letters('The cat in the hat')
-    'ehT tac ni eht tah'
-    >>> reverse_letters('The quick brown fox jumped over the lazy dog.')
-    'ehT kciuq nworb xof depmuj revo eht yzal .god'
-    >>> reverse_letters('Is this true?')
-    'sI siht ?eurt'
-    >>> reverse_letters("I   love       Python")
-    'I evol nohtyP'
+    Reverse all words that are longer than the given length of characters in a sentence.
+    If unspecified, length is taken as 0
+
+    >>> reverse_letters("Hey wollef sroirraw", 3)
+    'Hey fellow warriors'
+    >>> reverse_letters("nohtyP is nohtyP", 2)
+    'Python is Python'
+    >>> reverse_letters("1 12 123 1234 54321 654321", 0)
+    '1 21 321 4321 12345 123456'
+    >>> reverse_letters("racecar")
+    'racecar'
     """
-    return " ".join([word[::-1] for word in input_str.split()])
+    return " ".join(
+        "".join(word[::-1]) if len(word) > length else word for word in sentence.split()
+    )
 
 
 if __name__ == "__main__":
     import doctest
 
     doctest.testmod()
+    print(reverse_letters("Hey wollef sroirraw"))
diff --git a/strings/reverse_long_words.py b/strings/reverse_long_words.py
deleted file mode 100644
index 39ef11513..000000000
--- a/strings/reverse_long_words.py
+++ /dev/null
@@ -1,21 +0,0 @@
-def reverse_long_words(sentence: str) -> str:
-    """
-    Reverse all words that are longer than 4 characters in a sentence.
-
-    >>> reverse_long_words("Hey wollef sroirraw")
-    'Hey fellow warriors'
-    >>> reverse_long_words("nohtyP is nohtyP")
-    'Python is Python'
-    >>> reverse_long_words("1 12 123 1234 54321 654321")
-    '1 12 123 1234 12345 123456'
-    """
-    return " ".join(
-        "".join(word[::-1]) if len(word) > 4 else word for word in sentence.split()
-    )
-
-
-if __name__ == "__main__":
-    import doctest
-
-    doctest.testmod()
-    print(reverse_long_words("Hey wollef sroirraw"))
diff --git a/strings/strip.py b/strings/strip.py
new file mode 100644
index 000000000..d4f901f0c
--- /dev/null
+++ b/strings/strip.py
@@ -0,0 +1,33 @@
+def strip(user_string: str, characters: str = " \t\n\r") -> str:
+    """
+    Remove leading and trailing characters (whitespace by default) from a string.
+
+    Args:
+        user_string (str): The input string to be stripped.
+        characters (str, optional): Optional characters to be removed
+                (default is whitespace).
+
+    Returns:
+        str: The stripped string.
+
+    Examples:
+        >>> strip("   hello   ")
+        'hello'
+        >>> strip("...world...", ".")
+        'world'
+        >>> strip("123hello123", "123")
+        'hello'
+        >>> strip("")
+        ''
+    """
+
+    start = 0
+    end = len(user_string)
+
+    while start < end and user_string[start] in characters:
+        start += 1
+
+    while end > start and user_string[end - 1] in characters:
+        end -= 1
+
+    return user_string[start:end]