Add the Miller-Rabin primality test

Author: appgurueu

.spec/math/prime/miller_rabin_test_spec.lua

@@ -0,0 +1,25 @@
+describe("Miller-Rabin probabilistic primality test", function()
+	local test_prime = require("math.prime.miller_rabin_test")
+	it("works for small primes & composite numbers", function()
+		local is_prime = require("math.prime.sieve_of_eratosthenes")(1e6)
+		for number = 1, 100 do -- covers all edge cases
+			assert.equal(is_prime[number], test_prime(number, 20))
+		end
+		for _ = 1, 1e5 do
+			local number = math.random(101, 1e6)
+			assert.equal(is_prime[number], test_prime(number, 20))
+		end
+	end)
+	it("works for selected large primes", function()
+		for _, prime in pairs({ 6199, 7867, 2946901, 39916801 }) do
+			assert.truthy(test_prime(prime))
+		end
+	end)
+	it("works for random composite numbers", function()
+		for _ = 1, 1e3 do
+			-- Care is taken to stay well within double (and int) bounds
+			local composite = math.random(2 ^ 7, 2 ^ 13) * math.random(2 ^ 7, 2 ^ 13)
+			assert.falsy(test_prime(composite))
+		end
+	end)
+end)

src/math/prime/miller_rabin_test.lua

@@ -0,0 +1,60 @@
+-- Miller-Rabin primality test; nondeterministic variant
+-- Be careful when using this with doubles as precision issues may lead to incorrect results
+
+-- Simple integer exponentiation by squaring mod some number
+-- Apply mod after (almost) every operation to not run into issues with double precision
+local function modpow(base, exp, mod)
+	if exp == 1 then
+		return base
+	end
+	if exp % 2 == 1 then
+		return (modpow(base, exp - 1, mod) * base) % mod
+	end
+	return modpow(base, exp / 2, mod) ^ 2 % mod
+end
+
+return function(
+	n, -- number to test for primality; may not exceed (2^52)^.5 = 2^26 = 67108864 due to double limitations
+	rounds -- rounds determine accuracy: probability that a composite is considered probably prime does not exceed 4^-k
+)
+	-- Handle edge cases
+	if n == 1 then
+		return false
+	end
+	if n % 2 == 0 then
+		return n == 2
+	end
+	if n == 3 then
+		return true
+	end
+
+	rounds = rounds or 100 -- decent default for a false positive probability < 1e-60
+
+	-- Write n as d*2^r + 1
+	local d = n - 1
+	local r = 0
+	while d % 2 == 0 do
+		r = r + 1
+		d = d / 2
+	end
+
+	for _ = 1, rounds do
+		local a = math.random(2, n - 2)
+		local x = modpow(a, d, n) -- a^d % n
+		if x ~= 1 and x ~= n - 1 then
+			local composite = true
+			for _ = 1, rounds - 1 do
+				x = x ^ 2 % n
+				if x == n - 1 then
+					composite = false
+					break
+				end
+			end
+			if composite then
+				return false -- certainly not prime
+			end
+		end
+	end
+
+	return true -- likely prime (confidence based on rounds)
+end