Add Solovay-Strassen Primality test

maths/solovay_strassen_primality_test.py

@@ -0,0 +1,107 @@
+"""
+This script implements the Solovay-Strassen Primality test.
+
+This probabilistic primality test is based on Euler's criterion. It is similar
+to the Fermat test but uses quadratic residues. It can quickly identify
+composite numbers but may occasionally classify composite numbers as prime.
+
+More details and concepts about this can be found on:
+https://en.wikipedia.org/wiki/Solovay%E2%80%93Strassen_primality_test
+"""
+
+
+import random
+
+
+def jacobi_symbol(random_a: int, number: int) -> int:
+    """
+    Calculate the Jacobi symbol. The Jacobi symbol is a generalization
+    of the Legendre symbol, which can be used to simplify computations involving
+    quadratic residues. The Jacobi symbol is used in primality tests, like the
+    Solovay-Strassen test, because it helps determine if an integer is a
+    quadratic residue modulo a given modulus, providing valuable information
+    about the number's potential primality or compositeness.
+
+    Parameters:
+        random_a: A randomly chosen integer from 2 to n-2 (inclusive)
+        number: The number that is tested for primality
+
+    Returns:
+        jacobi_symbol: The Jacobi symbol is a mathematical function
+        used to determine whether an integer is a quadratic residue modulo
+        another integer (usually prime) or not.
+
+    >>> jacobi_symbol(2, 13)
+    -1
+    >>> jacobi_symbol(5, 19)
+    1
+    >>> jacobi_symbol(7, 14)
+    0
+    """
+
+    if random_a in (0, 1):
+        return random_a
+
+    random_a %= number
+    t = 1
+
+    while random_a != 0:
+        while random_a % 2 == 0:
+            random_a //= 2
+            r = number % 8
+            if r in (3, 5):
+                t = -t
+
+        random_a, number = number, random_a
+
+        if random_a % 4 == number % 4 == 3:
+            t = -t
+
+        random_a %= number
+
+    return t if number == 1 else 0
+
+
+def solovay_strassen(number: int, iterations: int) -> bool:
+    """
+    Check whether the input number is prime or not using
+    the Solovay-Strassen Primality test
+
+    Parameters:
+        number: The number that is tested for primality
+        iterations: The number of times that the test is run
+        which effects the accuracy
+
+    Returns:
+        result: True if number is probably prime and false
+        if not
+
+    >>> random.seed(10)
+    >>> solovay_strassen(13, 5)
+    True
+    >>> solovay_strassen(9, 10)
+    False
+    >>> solovay_strassen(17, 15)
+    True
+    """
+
+    if number <= 1:
+        return False
+    if number <= 3:
+        return True
+
+    for _ in range(iterations):
+        a = random.randint(2, number - 2)
+        x = jacobi_symbol(a, number)
+        y = pow(a, (number - 1) // 2, number)
+
+        if x == 0 or y != x % number:
+            return False
+
+    return True
+
+
+if __name__ == "__main__":
+    import doctest
+
+    doctest.testmod()