Prime number testing

[Rahul]

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).

rahul

 

[Paul Rubin]

Does anybody have a faster way to test for primality.(and not just
pseudoprime).

Yes, but the math gets a little fancy. It's a bit much for a news
post. You might look in Knuth vol. 2.

 

[Mel Wilson]

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.

 

[Paul Miller]

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.

Actually I find this function runs quite fast, although it terminates
with an UnboundLocalError.

What you probably meant was something like this:

import math

def isPrime (x):
if x % 2 == 0:
return 0
limit = int(math.sqrt (x))
for i in xrange (3,limit+1,2):
if x % i == 0:
return 0
return 1

I know your function attempted to sieve out multiples of 3 as well,
but that only improves performance in case x is composite.

Looking at this and assuming that 2**61 -1 is prime, you find that
this will do (2**61-1) or about 2**60 divisions. If you assume you
can perform 10**9 divisions per sec, then this will take about 36
years. That's why your function is slow.

You could improve this a bit by dividing by only prime numbers.
However, there are approximately 55997641048889848 +- 183279217216552
prime numbers <= 2**61-1, which means you'll now "only" be performing
about 2**55 divisions, a speedup of about 32 times, or approximately 1
year 2 months. The catch is you need a lot of memory to store all
those primes. [417214937 gigabytes, or so, provided you don't
introduce a compression scheme, which would slow you down a bit]

So, basically, trial division is hopeless for finding any but small
prime numbers. (See http://www.utm.edu/research/primes/howmany.shtml
if you'd like to check my numbers.)

In comparision
Mathematica not only tests but finds factors in an instant by
Divisors[x].

According to http://mathworld.wolfram.com/Rabin-MillerStrongPseudoprimeTest.html,
Mathematica implements primality testing thusly:

"Mathematica versions 2.2 and later have implemented the multiple
Rabin-Miller test combined with a Lucas pseudoprime test as the
primality test used by the function PrimeQ[n]. As of 1997, this
procedure is known to be correct only for all , but no counterexamples
are known and if any exist, they are expected to occur with extremely
small probability (i.e., much less than the probability of a hardware
error on a computer performing the test)."

Does anybody have a faster way to test for primality.(and not just
pseudoprime).

Yes. See above.

 

[Rahul]

Hi.
Thanks for all the info from all people.
Sorry for posting incorrect code.
Change some lines.
import math
(instead of importing sqrt from math)

and these two also

sqrt=math.sqrt(x)
sqrt=sqrt+1

Some clarifications :

Looking at this and assuming that 2**61 -1 is prime, you find that
this will do (2**61-1) or about 2**60 divisions. If you assume you
can perform 10**9 divisions per sec, then this will take about 36
years. That's why your function is slow.

I am diving upto square root.That is upto about 2**30 or around one
billion times.Assuming a modest speed of 0.1 million divisons per
second this should take 100000 secs or three hours.Though still slow.


I will look up knuth and also the links suggested.
Thanks for the illumination

Rahul

 

[Paul Miller]

import math

def isPrime (x):
if x % 2 == 0:
return 0
limit = int(math.sqrt (x))
for i in xrange (3,limit+1,2):
if x % i == 0:
return 0
return 1

I know your function attempted to sieve out multiples of 3 as well,
but that only improves performance in case x is composite.

Looking at this and assuming that 2**61 -1 is prime, you find that
this will do (2**61-1) or about 2**60 divisions. If you assume you
can perform 10**9 divisions per sec, then this will take about 36
years. That's why your function is slow.

This should of course be (2**61-1)/2 divisions.

The catch is you need a lot of memory to store all
those primes. [417214937 gigabytes, or so, provided you don't
introduce a compression scheme, which would slow you down a bit]

Coincidentally, BTW, this would be about the time Longhorn comes out,
and it's not far off the minimum spec for the OS, either.