Merge pull request #113 from 96Octavian/master

This PR implements the code for Miller-Rabin primality test in C#, Java, Python referred in issue #79

.gitignore

@@ -0,0 +1,5 @@
+################################################################################
+# This .gitignore file was automatically created by Microsoft(R) Visual Studio.
+################################################################################
+
+/.vs

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/MillerRabinMethod.cs

@@ -0,0 +1,163 @@
+using System.Numerics;
+using System.Security.Cryptography;
+
+namespace Algorithms.primalityTest
+{
+
+    public static partial class PrimalityTest
+    {
+
+        /// <summary>
+        /// Generate a random BigInteger between <paramref name="min"/> and <paramref name="max"/>, inclusive.
+        /// </summary>
+        /// <remarks>
+        /// Random BigInteger implementation from <see href="https://stackoverflow.com/a/48855115/12771343">this StackOverflow answer</see>
+        /// Uses an offset on the <see cref="RandomInRangeFromZeroToPositive(BigInteger)"/> function to return a random BigInteger in a range.
+        /// </remarks>
+        /// <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>
+        /// <remarks>
+        /// Random BigInteger implementation from <see href="https://stackoverflow.com/a/48855115/12771343">this StackOverflow answer</see>
+        /// Returns a random BigInteger lower than <paramref name="max"/> derived from a random byte array" />
+        /// </remarks>
+        /// <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;
+        }
+
+        /// <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 = 40)
+        {
+            // 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.IsEven)
+            {
+                d >>= 1;
+                ++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;
+        }
+
+
+        public static class MillerRabinMethod
+        {
+            static void Main(string[] args)
+            {
+                int count = 0;
+                int upper_bound = 1000;
+                Console.WriteLine($"Prime numbers lower than {upper_bound}:");
+                for (int i = 1; i < upper_bound; ++i)
+                {
+                    if (PrimalityTest.MillerRabin(i))
+                    {
+                        Console.WriteLine($"\t{i}");
+                        ++count;
+                    }
+                }
+                Console.WriteLine($"Total: {count}");
+            }
+        }
+    }
+
+}

primalityTest/millerRabinMethod/MillerRabinMethod.java

@@ -0,0 +1,148 @@
+import java.math.BigInteger;
+import java.security.SecureRandom;
+
+class MillerRabinMethod {
+
+    /**
+     *
+     * Generate a random BigInteger between {@code min} and {@code max}, inclusive.
+     * Random BigInteger implementation from <a href="https://stackoverflow.com/a/48855115/12771343">this StackOverflow answer</a>
+     * Uses an offset on the {@link #RandomInRangeFromZeroToPositive(BigInteger)}
+     * function to return a random BigInteger in a range.
+     *
+     * @param min lower bound
+     * @param max upper bound
+     * @return A random BigInteger in [{@code min}, {@code max}].
+     */
+    public static BigInteger RandomInRange(BigInteger min, BigInteger max) {
+
+        // If min > max we swap the values
+        if (min.compareTo(max) == 1) {
+            max = max.add(min);
+            min = max.subtract(min);
+            max = max.subtract(min);
+        }
+
+        // Offset to set min = 0
+        max = max.subtract(min);
+
+        // Retrieve the random number and offset to get the desired interval
+        BigInteger value = RandomInRangeFromZeroToPositive(max).add(min);
+        return value;
+    }
+
+    /**
+     *
+     * Generate a random BigInteger smaller or equal to {@code max}. Random
+     * BigInteger implementation from <a href="https://stackoverflow.com/a/48855115/12771343">this StackOverflow* answer</a>
+     * Returns a random BigInteger lower than {@code max} derived from a
+     * random byte array" />
+     *
+     * @param max upper bound
+     * @return a random BigInteger in [0, {@code max}].
+     */
+    private static BigInteger RandomInRangeFromZeroToPositive(BigInteger max) {
+
+        SecureRandom rng = new SecureRandom();
+        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[0];
+
+        // 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.nextBytes(bytes);
+
+            // Set most significant bits to 0 (because if any of these bits is 1 `value >
+            // max`)
+            bytes[0] &= 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.compareTo(max) == 1);
+
+        return value;
+    }
+
+    /**
+     *
+     * Miller-Rabin primality test
+     * Test if {@code number} could be prime using the Miller-Rabin primality test
+     * with {@code round} rounds. A return value of false means {@code number} is
+     * definitely composite, while true means it is probably prime.
+     * The higher {@code round} is, the more accurate the test is.
+     *
+     * @param number the number to be tested for primality
+     * @param rounds how many rounds to use in the test
+     * @return a Boolean indicating if the number could be prime or not
+     */
+    public static Boolean MillerRabin(BigInteger number, int rounds) {
+
+        // Handle corner cases
+        if (number.equals(BigInteger.valueOf(1)))
+            return false;
+        if (number.equals(BigInteger.valueOf(2)))
+            return true;
+
+        // Factor out the powers of 2 from {number - 1} and save the result
+        BigInteger d = number.subtract(BigInteger.valueOf(1));
+        int r = 0;
+        while ((d.getLowestSetBit() & 1) == 0) {
+            d = d.shiftLeft(1);
+            ++r;
+        }
+
+        BigInteger x, a;
+        // Cycle at most {round} times
+        for (; rounds > 0; --rounds) {
+            a = RandomInRange(BigInteger.valueOf(2), number.subtract(BigInteger.valueOf(1)));
+            x = a.modPow(d, number);
+            if (x.equals(BigInteger.valueOf(1)) || x.equals(number.subtract(BigInteger.valueOf(1))))
+                continue;
+            // Cycle at most {r - 1} times
+            for (int tmp = 0; tmp < r - 1; ++tmp) {
+                x = x.modPow(BigInteger.valueOf(2), number);
+                if (x.equals(number.subtract(BigInteger.valueOf(1))))
+                    break;
+            }
+            if (x.equals(number.subtract(BigInteger.valueOf(1))))
+                continue;
+            return false;
+        }
+        return true;
+    }
+
+    public static Boolean MillerRabin(BigInteger number) {
+        return MillerRabin(number, 40);
+    }
+
+    public static void main(String[] args) {
+        int count = 0;
+        int upper_bound = 1000;
+        System.out.println("Prime numbers lower than " + upper_bound + ":");
+        for (int i = 1; i < upper_bound; ++i) {
+            if (MillerRabin(BigInteger.valueOf(i))) {
+                System.out.println("\t" + i);
+                ++count;
+            }
+        }
+        System.out.println("Total: " + count);
+    }
+}

primalityTest/millerRabinMethod/miller_rabin_method.py

@@ -0,0 +1,56 @@
+import random
+
+
+def miller_rabin_method(number: int, rounds: int = 40) -> 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
+    if number == 3:
+        return True
+
+    # Factor out the powers of 2 from {number - 1} and save the result
+    d = number - 1
+    r = 0
+    while not d & 1:
+        d = d >> 1
+        r += 1
+
+    # 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
+
+
+if __name__ == "__main__":
+    count = 0
+    upper_bound = 1000
+    print(f"Prime numbers lower than {upper_bound}:")
+    for i in range(1, 1000):
+        if miller_rabin_method(i):
+            print(f"\t{i}")
+            count += 1
+    print(f"Total: {count}")