Tiawichi Research

I stumbled across a post today that demonstrated some straightforward python code for determining whether a number is prime:

def is_prime(n):
    n = abs(n)
    i = 2
    while i < n:
        if n % i == 0:
            return False
        i += 1
    return True

Clean, but slow, reaaaal slow. For the 7th (n=19) Mersenne prime, 524287, it determines primeness in a jiffy, 0.17 s. For the 8th (n=31), 2147483647, determining primeness requires 0:12:09.719914.

I think there are few areas in computer science that will bring out the ‘this is sooo sub-optimal’ comments as this subject as it has been studied quite thoroughly. So to impale myself on my own sword, here goes.

First cut is from code jerked from activestate here and here.

#!/usr/bin/python
import math

def isprime(aNumber):
	'''return True if the number is prime, false otherwise
           jacked from: http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/410662'''
	if aNumber < 2: return False
	if aNumber == 2: return True
	if (( aNumber / 2 ) * 2 == aNumber) :
		return False
	else:
		klist = primes(int(math.sqrt(aNumber+1)))
		for k in klist[1:]:
			if (( aNumber / k ) * k == aNumber ): return False
		return True

def primes(n):
        '''jacked from http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/366178
          Basically the well-known sieve of Eratosthenes but avoiding the sqrt function'''
	if n==2: return [2]
	elif n<2: return []
	s=range(3,n+1,2)
	mroot = n ** 0.5
	half=(n+1)/2-1
	i=0
	m=3
	while m <= mroot:
		if s[i]:
			j=(m*m-3)/2
			s[j]=0
			while j<half:
				s[j]=0
				j+=m
		i=i+1
		m=2*i+3
	return [2]+[x for x in s if x]

if __name__ == '__main__':
    import sys
    if isprime(int(sys.argv[1])):
        print 'yes'
    else:
        print 'no'

For the 8th (n=31), 2147483647, determining primeness requires 0:00:00.016553. Pretty darn quick. The 9th Mersenne prime, 2305843009213693951, throws an exception – out of memory.

Google some more and I come up with the following:

"""Prime number generation and number factorization."""

import bisect, itertools, random, sys

_primes = [2]

_miller_rabin_limit = 48611    # 5000th prime
_miller_rabin_security = 7

def modpow (a, b, c):
    """Efficiently compute (a^b)%c where a, b and c are positive integers."""
    if b == 1: return a % c
    d = b//2
    x = modpow (a, d, c)
    x = x*x%c
    if b % 2 == 1: x = x*a%c
    return x

def miller_rabin (n, t = _miller_rabin_security):
    """Apply Miller-Rabin primality test to detect whether n is prime or not.
    The test is run t times."""
    if n % 2 == 0: return False
    r = (n-1)//2
    for s in itertools.count (1):
        if r % 2: break
        r //= 2
    for tt in xrange (t):
        a = int (random.uniform (2, n-1))
        y = modpow (a, r, n)
        if y != 1 and y != n-1:
            for j in xrange (1, s):
                y = y**2 % n
                if y == 1: return False
                if y == n-1: break
            if y != n-1: return False
    return True

def _is_prime (i):
    if i > _miller_rabin_limit: return miller_rabin (i)
    s = int (i**0.5)
    for j in _primes:
        if i % j == 0: return False
        if j > s: return True
    return True

def _gen_primes ():
    for i in xrange (3, sys.maxint, 2):
        if _is_prime (i):
            _primes.append (i)
            yield i

_primary = _gen_primes ()

def _gen_primes (minval):
    i = 0
    if minval:
        l = _primes[-1]
        if minval > l:
            while minval > l: l =_primary.next ()
            i = len (_primes) - 1
        else:
            i = bisect.bisect_left (_primes, minval)
    for i in itertools.count (i):
        if i == len (_primes): yield _primary.next ()
        else: yield _primes[i]

def gen_primes (maxval = None, minval = None):
    """Prime numbers generator. If minval is given, only primes greater or
    equal to minval are returned. If maxval is given, only primes smaller
    or equal to maxval are returned.

    >>> for p in gen_primes(5): print p
    ...
    2
    3
    5
    7
    11"""
    if maxval:
        return itertools.takewhile (lambda x: x<=maxval, _gen_primes (minval))
    return _gen_primes (minval)

def gen_factors (n, duplicates = True):
    """Generator for factors of n (n > 1). If duplicates is False, do not
    send the same factor more than once."""
    assert n > 1
    if n > _miller_rabin_limit and miller_rabin (n):
        yield n
        return
    s = int (n**0.5)
    for i in gen_primes ():
        if i > s:
            yield n
            return
        if n % i == 0:
            yield i
            n //= i
            while n % i == 0:
                if duplicates: yield i
                n //= i
            if n == 1: return
            if n > _miller_rabin_limit and miller_rabin (n):
                yield n
                return
            s = int (n**0.5)

def factors (n):
    """Factorize n (n > 1) into its prime factors. Return a dictionary where
    keys are prime factors and values are powers.

    >>> factors(18)
    {2: 1, 3: 2}"""
    l = {}
    for i in gen_factors (n):
        try: l[i] += 1
        except KeyError: l[i] = 1
    return l

def factorslist (n):
    """Return a list of prime factors of n (n > 1).

    >>> factorslist(18)
    [2, 3, 3]"""
    return list (gen_factors (n))

def is_prime (n):
    """Check whether n is prime or not."""
    if n < 2: return False
    return gen_factors(n).next() == n

def ufactors (n):
    """Return the list of unique prime factors of n (n > 1).

    >>> ufactors(100)
    [2, 5]"""
    return list (gen_factors(n, duplicates = False))

def nfactors (n):
    """Return the number of unique prime factors of n."""
    return len (ufactors (n))

def totient (n):
    """Euler's totient function. Returns the number of integers between
    1 and n-1 relatively prime to n (n > 0).

    >>> totient(6)
    2
    (6 is relatively prime to 1 and 5)"""
    if n == 1: return 0
    num = den = 1
    for p in gen_factors (n, duplicates = False):
        num *= (p-1)
        den *= p
    return n * num // den

_s = {1: 1}
def s (n):
    """Sum of divisors of n (n > 0).

    >>> s(6)
    12
    (divisors of 6 are 1, 2, 3 and 6, summing to 12)"""
    if not _s.has_key (n):
        t = 1
        for (p, k) in factors(n).items():
            t *= (p**(k+1)-1)//(p-1)
        _s[n] = t
    return _s[n]

if __name__ == '__main__':
    import sys
    start = datetime.datetime.now()
    print start
    if is_prime(int(sys.argv[1])):
        print 'yes'
    else:
        print 'no'
    end = datetime.datetime.now()
    print end
    print(end - start)

For the 8th (n=31), 2147483647, determining primeness requires 0:00:00.001291. Runs like a scalded dog. Let’s try the 9th Mersenne prime, 2305843009213693951 – 0:00:00.002065. Fast, eh?

Let’s try the 6th Fermat number (18446744073709551617 = 18446744073709551617 = 274177 × 67280421310721). Fermat conjectured they were all prime but Euler factored #5 in 1792. Time to determine not prime: 0:00:10.347881.