Add Miller-Rabin primality test C# and Python implementation

Author: 96Octavian

Readme.md

@@ -61,7 +61,7 @@ Basically these are supposed to be my notes, but feel free to use them as you wi
 ### [Primality Test](primalityTest)
 - [Optimised School Method](primalityTest/optimisedSchoolMethod) (Check factors till √n)
 - [Fermat Method](primalityTest/fermatMethod)
-- Miller-Rabin Method
+- [Miller-Rabin Method](primalityTest/millerRabinMethod)
 - [Solovay-Strassen Method](primalityTest/solovayStrassenMethod)
 
 ### [Classic Alogorithms](ClassicalAlgos)

primalityTest/millerRabinMethod/AlgorithmsTools.cs

@@ -0,0 +1,90 @@
+using System.Numerics;
+using System.Security.Cryptography;
+
+namespace Algorithms.primalityTest
+{
+
+    /// <summary>
+    /// Collection of support tools for primality test algorithms
+    /// Random BigInteger implementation from https://stackoverflow.com/a/48855115/12771343
+    /// </summary>
+    public static partial class Algorithms
+    {
+
+        /// <summary>
+        /// Generate a random BigInteger between <paramref name="min"/> and <paramref name="max"/>, inclusive.
+        /// </summary>
+        /// <param name="min">BigInter lower bound.</param>
+        /// <param name="max">BigInteger upper bound.</param>
+        /// <returns>
+        /// A random BigInteger in [<paramref name="min"/>, <paramref name="max"/>].
+        /// </returns>
+        public static BigInteger RandomInRange(BigInteger min, BigInteger max)
+        {
+
+            // If min > max we swap the values
+            if (min > max)
+            {
+                max += min;
+                min = max - min;
+                max -= min;
+            }
+
+            // Offset to set min = 0
+            max -= min;
+
+            // Retrieve the random number and offset to get the desired interval
+            BigInteger value = RandomInRangeFromZeroToPositive(max) + min;
+            return value;
+        }
+
+        /// <summary>
+        /// Generate a random BigInteger smaller or equal to <paramref name="max"/>.
+        /// </summary>
+        /// <param name="max">BigInteger upper bound.</param>
+        /// <returns>
+        /// A random BigInteger in [0, <paramref name="max"/>].
+        /// </returns>
+        private static BigInteger RandomInRangeFromZeroToPositive(BigInteger max)
+        {
+            using RandomNumberGenerator rng = RandomNumberGenerator.Create();
+            BigInteger value;
+            byte[] bytes = max.ToByteArray();
+
+            // NOTE: sign bit is always 0 because `max` must always be positive
+            byte zeroBitsMask = 0b00000000;
+
+            // Count how many bits of the most significant byte are 0
+            var mostSignificantByte = bytes[bytes.Length - 1];
+
+            // Try to set to 0 as many bits as there are in the most significant byte, starting from the left (most significant bits first)
+            // NOTE: `i` starts from 7 because the sign bit is always 0
+            for (var i = 7; i >= 0; i--)
+            {
+                // Keep iterating until the most significant non-zero bit
+                if ((mostSignificantByte & (0b1 << i)) != 0)
+                {
+                    var zeroBits = 7 - i;
+                    zeroBitsMask = (byte)(0b11111111 >> zeroBits);
+                    break;
+                }
+            }
+
+            do
+            {
+                rng.GetBytes(bytes);
+
+                // Set most significant bits to 0 (because if any of these bits is 1 `value > max`)
+                bytes[bytes.Length - 1] &= zeroBitsMask;
+
+                value = new BigInteger(bytes);
+
+                // `value > max` 50% of the times, in which case the fastest way to keep the distribution uniform is to try again
+            } while (value > max);
+
+            return value;
+        }
+
+
+    }
+}

primalityTest/millerRabinMethod/PrimalityTest.cs

@@ -0,0 +1,59 @@
+using System.Numerics;
+
+namespace Algorithms.primalityTest
+{
+
+    public static partial class Algorithms
+    {
+
+        /// <summary>
+        /// Miller-Rabin primality test
+        /// </summary>
+        /// <remarks>
+        /// Test if <paramref name="number"/> could be prime using the Miller-Rabin Primality Test with <paramref name="rounds"/> rounds.
+        /// A return value of false means <paramref name="number"/> is definitely composite, while true means it is probably prime.
+        /// The higher <paramref name="rounds"/> is, the more accurate the test is.
+        /// </remarks>
+        /// <param name="number">The number to be tested for primality.</param>
+        /// <param name="rounds">How many rounds to use in the test.</param>
+        /// <returns>A bool indicating if the number could be prime or not.</returns>
+        public static bool MillerRabin(BigInteger number, int rounds)
+        {
+            // Handle corner cases
+            if (number == 1)
+                return false;
+            if (number == 2)
+                return true;
+
+            // Factor out the powers of 2 from {number - 1} and save the result
+            BigInteger d = number - 1;
+            int r = 0;
+            while (d % 2 == 0)
+            {
+                d /= 2;
+                ++r;
+            }
+
+            BigInteger x, a;
+            // Cycle at most {round} times
+            for (; rounds > 0; --rounds)
+            {
+                a = RandomInRange(2, (number - 1));
+                x = BigInteger.ModPow(a, d, number);
+                if (x == 1 || x == number - 1)
+                    continue;
+                // Cycle at most {r - 1} times
+                for (int tmp = 0; tmp < r - 1; ++tmp)
+                {
+                    x = BigInteger.ModPow(x, 2, number);
+                    if (x == number - 1)
+                        break;
+                }
+                if (x == number - 1)
+                    continue;
+                return false;
+            }
+            return true;
+        }
+    }
+}

primalityTest/millerRabinMethod/miller_rabin.py

@@ -0,0 +1,44 @@
+import random
+
+
+def miller_rabin(number: int, rounds: int) -> bool:
+    """Miller-Rabin primality test
+
+    Test if number could be prime using the Miller-Rabin Primality Test with rounds rounds.
+    A return value of false means number is definitely composite, while true means it is probably prime.
+    The higher rounds is, the more accurate the test is.
+    :param int number: The number to be tested for primality.
+    :param int rounds: How many rounds to use in the test.
+    :return: A bool indicating if the number could be prime or not.
+    :rtype: bool
+
+    """
+    # Handle corner cases
+    if number == 1:
+        return False
+    if number == 2:
+        return True
+
+    # Factor out the powers of 2 from {number - 1} and save the result
+    d = number - 1
+    r = 0
+    while d % 2 == 0:
+        d /= 2
+        r += 1
+    d = int(d)
+
+    # Cycle at most {round} times
+    for i in range(rounds + 1):
+        a = random.randint(2, number - 2)
+        x = pow(a, d, number)
+        if x == 1 or x == number - 1:
+            continue
+        # Cycle at most {r - 1} times
+        for e in range(r):
+            x = x * x % number
+            if x == number - 1:
+                break
+        if x == number - 1:
+            continue
+        return False
+    return True