Add Solovay-Strassen Primality test

maths/solovay_strassen_primality_test.py

@@ -0,0 +1,105 @@
+"""
+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 mathematical function
+    used to determine whether an integer is a quadratic residue modulo
+    another integer (usually prime) or not.
+
+    Parameters:
+        random_a (int): A randomly chosen integer from 2 to n-2 (inclusive)
+        number (int): The number that is tested for primality
+
+    Returns:
+        jocobi_symbol (int): The jacobi symbol
+
+    >>> jacobi_symbol(2, 13)
+    -1
+    >>> jacobi_symbol(5, 19)
+    1
+    >>> jacobi_symbol(7, 14)
+    0
+    """
+
+    if random_a == 0:
+        return 0
+    if random_a == 1:
+        return 1
+
+    random_a = random_a % number
+    t = 1
+
+    while random_a != 0:
+        while random_a % 2 == 0:
+            random_a = 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 = random_a % number
+
+    if number == 1:
+        return t
+    else:
+        return 0
+
+
+def solovay_strassen(number: int, accuracy: int) -> bool:
+    """
+    Check whether the input number is prime or not using
+    the Solovay-Strassen Primality test
+
+    Parameters:
+        number (int): The number that is tested for primality
+        accuracy (int): The accuracy for the test
+
+    Returns:
+        result (bool): True if number is probably prime and false
+        if composite
+
+    >>> 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(accuracy):
+        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()