Prime number testing

[Mel Wilson]

In article <9d234d68.0407102202.47b5cb43 at posting.google.com>,
codelykos at yahoo.com (Rahul) wrote:
>HI.
>Python , with its support of arbit precision integers, can be great
>for number theory. So i tried writing a program for testing whether a
>number is prime or not. But then found my function painfully slow.Here
>is the function :
>
>from math import sqrt
>def isPrime(x):
>    if x%2==0 or x%3==0:
>       return 0
>    i=5
>    sqrt=sqrt(x)
>    sqrt=sqrt(x)+1
>    while i<sqrt:
>       if x%i==0:
>          return 0
>       i=i+2
>       if x%i==0:
>          return 0
>       i=i+4
>    return 1
>
>Try testing some number like 2**61-1.Its very slow.In comparision
>Mathematica not only tests but finds factors in an instant by
>Divisors[x].
>
>Does anybody have a faster way to test for primality.(and not just
>pseudoprime).

   How long does it take to count to
100,000,000,000,000,000,000 by 6?  If Mathematica is that
fast, I'd think it's using probabilistic techniques, for
example, from Knuth:




import whrandom


def bigrandom (n=1000000L):
    return long (whrandom.uniform (2, n))


# Primality tests from Knuth Vol.II 4.5.4
# Algorithm P, etc.
#
# Return false if n is composite,  true if n may be prime
#
def maybeprime (n):
    if (n & 1) == 0:
        return 0        # n is even.. definitely not prime

    # find q and k, such that n = q*(2**k)+1
    q = n - 1
    k = 0
    while (q & 1) == 0:
        k += 1
        q >>= 1

    x = bigrandom (n)   # get a random number for this test

    y = pow (x, q, n)
    if y == 1  or  y == n-1:
        return 1        # n may be prime

    for j in xrange (1, k+1):
        y = pow (y, 2, n)
        if y == 1:      return 0        # n is definitely not prime
        if y == n-1:    return 1        # n may be prime

    return 0        # n is definitely not prime
# end maybeprime


# Return false if n is composite, true if n is very probably prime.
# The reliablity of the test is increased by increasing t
#
def isprime (n, t=25):
    for j in range (t):
        if not maybeprime (n):
            return 0    # n is known to be composite
    return 1    # n has survived t primality tests, is probably prime
# end isprime



        Regards.        Mel.