Included is a recipe for performing the Rabin-Miller probabilistic test for a composite witness. As provided by Paul Miller in the comments (not the Miller in Rabin-Miller), Rabin-Miller can only tell us if a value is definitely composite. In the case where a test value is not a witness for the compositeness of a potential prime, it can only lie with a probability of at most 1/4.
With this, we can attempt to catch a liar over some number of trials, and the probability of us not catching at least one liar after k trials (if the number is not actually prime) is at most 4**-k.
Included is an algorithm for generating a number of b bits for which no composite witness was found after k trials. Removing mathematical rigor will suggest that the probability of the value being prime after k trials is at least 1-1/4**k.
Download
Copy to clipboard1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 | '''
Rabin-Miller algorithm for testing the primality of a given number, and
an associated algorithm for generating a b-bit integer that is probably
prime.
Algorithm described in various texts, among them Algorithm Design by
Goodrich and Tamassia.
Copyleft 2005, Josiah Carlson
Use this recipe as you see fit.
'''
import random
import sys
DEBUG = False
def ipow(a,b,n):
#calculates (a**b)%n via binary exponentiation, yielding itermediate
#results as Rabin-Miller requires
A = a = long(a%n)
yield A
t = 1L
while t <= b:
t <<= 1
#t = 2**k, and t > b
t >>= 2
while t:
A = (A * A)%n
if t & b:
A = (A * a) % n
yield A
t >>= 1
def RabinMillerWitness(test, possible):
#Using Rabin-Miller witness test, will return True if possible is
#definitely not prime (composite), False if it may be prime.
return 1 not in ipow(test, possible-1, possible)
smallprimes = (3,5,7,11,13,17,19,23,29,31,37,41,43,
47,53,59,61,67,71,73,79,83,89,97)
def generate_prime(b, k=None):
#Will generate an integer of b bits that is probably prime.
#Reasonably fast on current hardware for values of up to around 512 bits.
bits = int(b)
assert bits > 1
if k is None:
k = 2*bits
k = int(k)
if k < 64:
k = 64
if DEBUG: print "(b=%i, k=%i)"%(bits,k),
good = 0
while not good:
possible = random.randrange(2**(bits-1)+1, 2**bits)|1
good = 1
if DEBUG: sys.stdout.write(';');sys.stdout.flush()
for i in smallprimes:
if possible%i == 0:
good = 0
break
else:
for i in xrange(k):
test = random.randrange(2, possible)|1
if RabinMillerWitness(test, possible):
good = 0
break
if DEBUG: sys.stdout.write('.');sys.stdout.flush()
if DEBUG: print
return possible
if __name__ == '__main__':
DEBUG = True
print "; - trying new possible prime"
print ". - failed to find a composite witness for this trial"
print
bits = 32
while bits < 547:
print hex(generate_prime(bits))
bits = int(bits*1.5)
|
If you are thinking "hey, it would be nice if I had a big prime number, or at least a number that is probably prime", you probably already have a use for them. Most uses for large primes involve cryptography.
I personally have no use for them, but I thought it would be nice to at least have a tester and generator, because you just never know when one would be useful to have.
From what I understand (very little), Rabin Miller is one of the better methods for testing for whether a number is composite. If you are curious about other methods, search for "Jacobi", "Solovay-Strassen test", and/or "Carmichael numbers".
1-1/2k ?** Note that Miller-Rabin is actually a compositeness test more than a primality test. That is, it can only produce a certificate of compositeness for a given integer rather than a certificate of primality. The guarantee on Miller-Rabin is actually that if n is the input to the algorithm, at most 1/4 of the possible witnesses to the compositeness of n will be "liars."
Of course, it's a mistake to say that n is "probably prime," since it's either prime or composite, but letting down our guard a bit on rigor, we can say that if n passes k trials of Miller-Rabin, it is prime with "probability" 1 - 1/4*k, rather than 1-1/2 * k.
I've updated the various portions to be more accurate. Thank you.
Probability of success. These statements are actually incorrect; in practise RM is actually much better, especially for larger bit sizes. See 4.4 on page 146 in the Handbook of Applied Cryptography, online at http://www.cacr.math.uwaterloo.ca/hac/
Implementation. I am a newby in Python but an expert on RM, and from the source I couldn't understand how the algorithm actually works. Or maybe it is a variation that I don't know. The algo I know writes N-1 as 2^u*v with v odd, calculates w^v (modN) and squares this value u times to see if a -1 appears; if not w is a witness of compositeness. According to my Pythonbook pow(x,y,z) = pow(x,y)%z; I think that in combination with this algorithm it would be faster.
My own implementation of RM. Successfully tested on this problem: http://spoj.sphere.pl/problems/PON/
</pre>
Sorry... My prev. comment got doubled and cut. Due to... hmm... imo, due to a free nature of the site. That's it... the free (incl. linux) just sucks. And those people derise at M$ software. Silly.
Plz delete both these my comments.
... Comments cannot be deleted, and crapping on free software randomly, on a site that gives away free software, is pretty foolish.
re: Implementation. I implemented a variant of the Rabin Miller algorithm as provided in Algorithm Design by Goodrich and Tamassia, which is stated in the comment preceeding the content of the recipe.
According to tests I ran initially, repeated calls to pow(x,,z) is slower than ipow(x,,z) for certain large values of x, *, and z, because while ipow generates only those values necessary for the primality test, one ends up recalculating them with the pow solution, unless one uses pow instead of the multiplication/modulo that I use, but then you end up with function call overhead, which is slower than the standard mathematical operations.
I went and read the page (and the few pages before and after), and from what I understand (which isn't a whole lot, given the minimal time at 2AM I spent reading), the probability of failure can be caluclated, and that the actual probability of failure tends to be far lower than (1/4)t, (1/2)(4.4*t) or lower when testing 100-binary digit possible primes with 10 tests.
One thing to note is that I end up doing a HUGE number of tests (I had used 'k' as the number of tests, but that is generally used as # of bits) in comparison to what is generally suggested as being 'good enough'. At least 64 in fact. Take, for example, a 512 bit possible prime, and the minimum number of tests that I do, 64. Plugging that into 4.48 (iii) gets us approximately (1/2)136, which is less than (1/4)64. Using the default number of tests as bits2, puts us generally into rule (iii), which gives failure probability that tends far lower than (1/4)*t, for most possible primes which would be interesting.
Certainly there are diminishing returns for each additional test one does (meaning that we could probably be happy with doing I went and read the page (and the few pages before and after), and from what I understand (which isn't a whole lot, given the minimal time at 2AM I spent reading), the probability of failure can be caluclated, and that the actual probability of failure tends to be far lower than (1/4)t, (1/2)(4.4*t) or lower when testing 100-binary digit possible primes with 10 tests.
One thing to note is that I end up doing a HUGE number of tests (I had used 'k' as the number of tests, but that is generally used as # of bits) in comparison to what is generally suggested as being 'good enough'. At least 64 in fact. Take, for example, a 512 bit possible prime, and the minimum number of tests that I do, 64. Plugging that into 4.48 (iii) gets us approximately (1/2)136, which is less than (1/4)64. Using the default number of tests as bits2, puts us generally into rule (iii), which gives failure probability that tends far lower than (1/4)*t, for most possible primes which would be interesting.
Certainly there are diminishing returns for each additional test one does (meaning that we could probably be happy with doing
Darn it, I should have previewed that entry. Use of < without proper HTML escapes killed it.
In any case, to finish the post, I could probably get away with doing fewer than 20 tests, regardless of the size of the integer desired.
PLEASE any one have rabin miller test for prime number with MPI source code c? please help and send it?! or email it