60 lines
1.3 KiB
Lua
60 lines
1.3 KiB
Lua
function MRIsPrime(n, k)
|
|
-- If n is prime, returns true (without fail).
|
|
-- If n is not prime, then returns false with probability ≥ 4^(-k), true otherwise.
|
|
-- First, deal with 1 and multiples of 2.
|
|
if n == 1 then
|
|
return false
|
|
elseif n == 2 then
|
|
return true
|
|
elseif n%2 == 0 then
|
|
return false
|
|
end
|
|
-- Now n is odd and greater than 1.
|
|
-- Find the unique expression n = 1 + t*2^h for n, where t is odd and h≥1.
|
|
t = (n-1)/2
|
|
h = 1
|
|
while t%2 == 0 do
|
|
t = t/2
|
|
h = h + 1
|
|
end
|
|
for i = 1, k do
|
|
-- Generate a random integer between 1 and n-1 inclusive.
|
|
a = math.random(n-1)
|
|
-- Test whether a is an element of the set L, and return false if not.
|
|
if not IsInL(n, a, t, h) then
|
|
return false
|
|
end
|
|
end
|
|
-- All generated a were in the set L; thus with high probability, n is prime.
|
|
return true
|
|
end
|
|
|
|
function IsInL(n, a, t, h)
|
|
local b = PowerMod(a, t, n)
|
|
if b == 1 then
|
|
return true
|
|
end
|
|
for j = 0, h-1 do
|
|
if b == n-1 then
|
|
return true
|
|
elseif b == 1 then
|
|
return false
|
|
end
|
|
b = (b^2)%n
|
|
end
|
|
return false
|
|
end
|
|
|
|
function PowerMod(x, y, m)
|
|
-- Computes x^y mod m.
|
|
local z = 1
|
|
while y > 0 do
|
|
if y%2 == 0 then
|
|
x, y, z = (x^2)%m, y//2, z
|
|
else
|
|
x, y, z = (x^2)%m, y//2, (x*z)%m
|
|
end
|
|
end
|
|
return z
|
|
end
|