diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e915029 --- /dev/null +++ b/.gitignore @@ -0,0 +1,5 @@ +################################################################################ +# This .gitignore file was automatically created by Microsoft(R) Visual Studio. +################################################################################ + +/.vs diff --git a/Readme.md b/Readme.md index 2f7d2d1..4d095d0 100644 --- a/Readme.md +++ b/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) diff --git a/primalityTest/millerRabinMethod/MillerRabinMethod.cs b/primalityTest/millerRabinMethod/MillerRabinMethod.cs new file mode 100644 index 0000000..aef5e30 --- /dev/null +++ b/primalityTest/millerRabinMethod/MillerRabinMethod.cs @@ -0,0 +1,163 @@ +using System.Numerics; +using System.Security.Cryptography; + +namespace Algorithms.primalityTest +{ + + public static partial class PrimalityTest + { + + /// + /// Generate a random BigInteger between and , inclusive. + /// + /// + /// Random BigInteger implementation from this StackOverflow answer + /// Uses an offset on the function to return a random BigInteger in a range. + /// + /// BigInter lower bound. + /// BigInteger upper bound. + /// + /// A random BigInteger in [, ]. + /// + 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; + } + + /// + /// Generate a random BigInteger smaller or equal to . + /// + /// + /// Random BigInteger implementation from this StackOverflow answer + /// Returns a random BigInteger lower than derived from a random byte array" /> + /// + /// BigInteger upper bound. + /// + /// A random BigInteger in [0, ]. + /// + 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; + } + + /// + /// Miller-Rabin primality test + /// + /// + /// Test if could be prime using the Miller-Rabin primality test with rounds. + /// A return value of false means is definitely composite, while true means it is probably prime. + /// The higher is, the more accurate the test is. + /// + /// The number to be tested for primality. + /// How many rounds to use in the test. + /// A bool indicating if the number could be prime or not. + 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}"); + } + } + } + +} diff --git a/primalityTest/millerRabinMethod/MillerRabinMethod.java b/primalityTest/millerRabinMethod/MillerRabinMethod.java new file mode 100644 index 0000000..ee47f3f --- /dev/null +++ b/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 this StackOverflow answer + * 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 this StackOverflow* answer + * 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); + } +} \ No newline at end of file diff --git a/primalityTest/millerRabinMethod/miller_rabin_method.py b/primalityTest/millerRabinMethod/miller_rabin_method.py new file mode 100644 index 0000000..2ea81a0 --- /dev/null +++ b/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}") \ No newline at end of file