Нахождение премьер с Миллер Рабин

Question

У меня есть то, что я считаю правильной реализацией алгоритма Миллера-Рабина с использованием Lua, и я пытаюсь получить последовательный возврат для простых чисел.Кажется, моя реализация работает только половину времени.Хотя, если я попытаюсь реализовать подобный код в Python, этот код работает 100% времени.Может ли кто-нибудь указать мне правильное направление?

--decompose n-1 as (2^s)*d
local function decompose(negOne)
  exponent, remainder = 0, negOne
  while (remainder%2) == 0 do 
    exponent = exponent+1
    remainder = remainder/2
  end
  assert((2^exponent)*remainder == negOne and ((remainder%2) == 1), "Error setting up s and d value")
  return exponent, remainder
end

local function isNotWitness(n, possibleWitness, exponent, remainder)
  witness = (possibleWitness^remainder)%n

  if (witness == 1) or (witness == n-1) then
    return false
  end

  for _=0, exponent do
    witness = (witness^2)%n
    if witness == (n-1) then
      return false
    end
  end

  return true
end

--using miller-rabin primality testing
--n the integer to be tested, k the accuracy of the test
local function isProbablyPrime(n, accuracy)
  if n <= 3 then
    return n == 2 or n == 3
  end
  if (n%2) == 0 then
    return false
  end

  exponent, remainder = decompose(n-1)

  --checks if it is composite
  for i=0, accuracy do
    math.randomseed(os.time())
    witness = math.random(2, n - 2)
    if isNotWitness(n, witness, exponent, remainder) then
      return false
    end
  end

  --probably prime
  return true
end

if isProbablyPrime(31, 30) then
  print("prime")
else
  print("nope")
end

Answer 1

Python имеет произвольные длины целых чисел.Lua этого не делает.
Проблема в witness = (possibleWitness^remainder)%n.
Lua не может напрямую рассчитать точный результат 29^15 % 31.
Для чисел n < sqrt(2^53):

* существует обходной путь1008 *

где

local function mulmod(a, e, m)
   local result = 1
   while e > 0 do
      if e % 2 == 1 then 
         result = result * a % m
         e = e - 1
      end
      e = e / 2
      a = a * a % m
   end
   return result
end