Geek asks and answers

Implementing the Rabin-Miller Strong Pseudoprime Test will not work

Today we tried to implement the powerful Podudoprime Rabin-Miller test.

Used Wolfram Mathworld as a reference, lines 3-5 summarize my code pretty much.

However, when I run the program, it says (sometimes) that prime numbers (even as low as 5, 7, 11) are not prime. I looked at the code for a very long time and I can not understand what is wrong.

For reference, I looked at this site in the same way as many other sites, but most of them use a different definition (maybe the same thing, but since I'm new to this math, I don't see the same obvious connection).

My code is:

import random

def RabinMiller(n, k):

    # obviously not prime
    if n < 2 or n % 2 == 0:
        return False

    # special case        
    if n == 2:
        return True

    s = 0
    r = n - 1

    # factor n - 1 as 2^(r)*s
    while r % 2 == 0:
        s = s + 1
        r = r // 2  # floor

    # k = accuracy
    for i in range(k):
        a = random.randrange(1, n)

        # a^(s) mod n = 1?
        if pow(a, s, n) == 1:
            return True

        # a^(2^(j) * s) mod n = -1 mod n?
        for j in range(r):
            if pow(a, 2**j*s, n) == -1 % n:
                return True

    return False

print(RabinMiller(7, 5))

How does this differ from the definition given in Mathworld?

+5
5

1.

, , , , , .

  • s = 0
r = n - 1

# factor n - 1 as 2^(r)*s
while r % 2 == 0:
    s = s + 1
    r = r // 2  # floor

    r s swapped: n-1 2 s r. MathWorld, r s :

    # factor n - 1 as 2^(r)*s, where s is odd.
r, s = 0, n - 1
while s % 2 == 0:
    r += 1
    s //= 2


  • for i in range(k):

    i : _.

  • 1 n - 1 :

    a = random.randrange(1, n)

    , MathWorld, . 1, 1 s= 1 (mod n), . n - 1, s (n - 1) s= -1 (mod n). , , :

    a = random.randrange(2, n - 1)

    ( 4, , True , n = 3, , n = 2.)

  • , MathWorld. , "n ", , "n a". , . , , a s= 1 (mod n), .

    # a^(s) = 1 (mod n)?
x = pow(a, s, n)
if x == 1:
    continue

  • . x, , 2 0 s (mod n). , :

    # a^(s) = ±1 (mod n)?
x = pow(a, s, n)
if x == 1 or x == n - 1:
    continue

  • , a 2 j s (mod n), ( n). , . :

    # a^(2^(j) * s) = -1 (mod n)?
for _ in range(r - 1):
    x = pow(x, 2, n)
    if x == n - 1:
        break
else:
    return False

  • , -. , 1977 :

    . p < N, , , N = 1000.

2.

:

from random import randrange

small_primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31] # etc.

def probably_prime(n, k):
    """Return True if n passes k rounds of the Miller-Rabin primality
    test (and is probably prime). Return False if n is proved to be
    composite.

    """
    if n < 2: return False
    for p in small_primes:
        if n < p * p: return True
        if n % p == 0: return False
    r, s = 0, n - 1
    while s % 2 == 0:
        r += 1
        s //= 2
    for _ in range(k):
        a = randrange(2, n - 1)
        x = pow(a, s, n)
        if x == 1 or x == n - 1:
            continue
        for _ in range(r - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break
        else:
            return False
    return True
+6

, Omri Barel, for. true, a, . , a n .

+4

:

# factor n - 1 as 2^(r)*s
while r % 2 == 0:
    s = s + 1
    r = r // 2  # floor

n = 7. , n - 1 = 6. n - 1 2^1 * 3. r = 1 s = 3.

- . r = 6, r % 2 == 0. s = 0, s = 1 r = 3. r % 2 != 0 .

s = 1 r = 3, : 2^r * s = 8.

s . , 2 ( r), s. n = 7, n - 1 = 6 ( r = 1), 3 (so s = 3).

+2

:

# miller-rabin pseudoprimality checker

from random import randrange

def isStrongPseudoprime(n, a):
    d, s = n-1, 0
    while d % 2 == 0:
        d, s = d/2, s+1
    t = pow(a, d, n)
    if t == 1:
        return True
    while s > 0:
        if t == n - 1:
            return True
        t, s = pow(t, 2, n), s - 1
    return False

def isPrime(n, k):
    if n % 2 == 0:
        return n == 2
    for i in range(1, k):
        a = randrange(2, n)
        if not isStrongPseudoprime(n, a):
            return False
    return True

, .

+2

Wikipedia, "" .

  • n < 1,373,653, a = 2 3;
  • if n <9,080,191, it is enough to check a = 31 and 73;
  • if n <4,759,123,141, it is enough to check a = 2, 7 and 61;
  • if n <2,152,302,898,747, it is enough to check a = 2, 3, 5, 7 and 11;
  • if n <3,474,749,660,383, it is enough to check a = 2, 3, 5, 7, 11 and 13;
  • if n <341,550,071,728,321, it is enough to check a = 2, 3, 5, 7, 11, 13 and 17;
+1
source

More articles:


All Articles